Home Original page

Python, MATLAB, Julia, R code: Chapter 10

Intro to Probability for Data Science

Chapters

Lecture Videos + Slides

Selected Exercise Videos

Code and Data

Problem Solutions (Instructor Only)

Purdue Courses

Acknowledgement: The Julia code is written by the contributors listed here.

Acknowledgement: The R code is written by contributors listed here.

Chapter 10.2 Mean and correlation functions

Example 10.5

% MATLAB code for Example 10.5
x = zeros(1000,20);
t = linspace(-2,2,1000);
for i=1:20
    X(:,i) = rand(1)*cos(2*pi*t);
end
plot(t, X, 'LineWidth', 2, 'Color', [0.8 0.8 0.8]); hold on;
plot(t, 0.5*cos(2*pi*t), 'LineWidth', 4, 'Color', [0.6 0 0]);

x = np.zeros((1000,20))
t = np.linspace(-2,2,1000)
for i in range(20):
  x[:,i] = np.random.rand(1)*np.cos(2*np.pi*t)
plt.plot(t,x,color='gray')
plt.plot(t,0.5*np.cos(2*np.pi*t),color='red')
plt.show()

using Plots

x = zeros(1000, 20)
t = range(-2, 2, length=1000)

for col in eachcol(x)
    col .= rand() * cos.(2*pi * t)
end

plot(t, x, linewidth=2, color=:gray80, legend=false)
plot!(t, 0.5 * cos.(2*pi * t), linewidth=4, color=:darkred)

x = matrix(0, 1000, 20)
t = seq(-2, 2, (2+2)/999)

for (i in 1:20) {
    x[,i] = runif(1) * cos(2 * pi * t)
}

matplot(t, x, lwd=2, col="gray")
grid()
lines(t, 0.5*cos(2*pi*t), lwd=4, col="darkred")

Example 10.6

% MATLAB code for Example 10.6
x = zeros(1000,20);
t = linspace(-2,2,1000);
for i=1:20
X(:,i) = cos(2*pi*t+2*pi*rand(1));
end
plot(t, X, 'LineWidth', 2, 'Color', [0.8 0.8 0.8]); hold on;
plot(t, 0*cos(2*pi*t), 'LineWidth', 4, 'Color', [0.6 0 0]);

x = np.zeros((1000,20))
t = np.linspace(-2,2,1000)
for i in range(20):
  Theta = 2*np.pi*(np.random.rand(1))
  x[:,i] = np.cos(2*np.pi*t+Theta)
plt.plot(t,x,color='gray')
plt.plot(t,np.zeros((1000,1)),color='red')
plt.show()

using Plots

x = zeros(1000, 20)
t = range(-2, 2, length=1000)

for col in eachcol(x)
    col .= cos.(2*pi*t .+ 2*pi*rand())
end

plot(t, x, linewidth=2, color=:gray80, legend=false)
hline!([0], linewidth=4, color=:darkred)

x = matrix(0, 1000, 20)
t = seq(-2, 2, (2+2)/999)

for (i in 1:20) {
    x[,i] = cos(2*pi*t+2*pi*runif(1))
}

matplot(t, x, lwd=2, col="gray")
grid()
lines(t, 0*cos(2*pi*t), lwd=4, col="darkred")

Example 10.7

% MATLAB code for Example 10.7
t = 0:20;
for i=1:20
    X(:,i) = rand(1).^t;
end
stem(t, X, 'LineWidth', 2, 'Color', [0.8 0.8 0.8]); hold on;
stem(t, 1./(t+1), 'LineWidth', 2, 'MarkerSize', 8);

using Plots

x = zeros(21, 20)
t = 0:20

for col in eachcol(x)
    col .= rand() .^ t
end

plot(t, x, line=:stem, linewidth=2, shape=:circle, color=:gray80, legend=false)
plot!(t, 1 ./ (t .+ 1), line=:stem, linewidth=2, shape=:circle, markersize=6, color=:darkred)

x = matrix(0, 21, 20)
t = 0:20

for (i in 1:20) {
    x[,i] = runif(1) ^ t
}

stem <- function(x,y,pch=1,...){
    if (missing(y)){
        y = x
        x = 1:length(x)
    }
    for (i in 1:length(x)){
        lines(c(x[i],x[i]), c(0,y[i]), ...)
    }
    lines(c(x[1]-2,x[length(x)]+2), c(0,0), ...)
}

matplot(t, x, pch=1, col="gray", lwd=2)
grid()

for (i in 1:ncol(x)) {
    stem(t, x[,i], col="gray", lwd=2)
}
stem(t, 1/(t+1), col="darkred", lwd=2, pch=1, type="o")

Example 10.11

% MATLAB code for Example 10.11: Plot the time function
t = linspace(-2,2,1000);
for i=1:20
X(:,i) = rand(1)*cos(2*pi*t);
end
plot(t, X, 'LineWidth', 2, 'Color', [0.8 0.8 0.8]); hold on;
plot(t, 0.5*cos(2*pi*t), 'LineWidth', 4, 'Color', [0.6 0 0]);
scatter(zeros(20,1), X(501,:)','LineWidth',2,'MarkerEdgeColor',[1 0.5 0]);
scatter(0.5*ones(20,1), X(626,:)','x','LineWidth',2,'SizeData',70,'MarkerEdgeColor',[0 0.5 1]);
grid on;
xticks(-2:0.5:2);
axis([-2 2 -1.2 1.2]);
set(gcf, 'Position', [100, 100, 600, 400]);
set(gca,'FontWeight','bold','fontsize',14);

% MATLAB code for Example 10.11: Plot the autocorrelation function
t = linspace(-1,1,1000);
R = (1/3)*cos(2*pi*t(:)).*cos(2*pi*t);
imagesc(t,t,R);

import numpy as np
import matplotlib.pyplot as plt
t = np.linspace(-1,1,1000)
R = (1/3)*np.outer(np.cos(2*np.pi*t), np.cos(2*np.pi*t))
plt.imshow(R, extent=[-1, 1, -1, 1])
plt.show()

using Plots

x = zeros(1000, 20)
t = range(-2, 2, length=1000)

for col in eachcol(x)
    col .= rand() * cos.(2*pi * t)
end

plot(t, x, linewidth=2, color=:gray80, legend=false)
plot!(t, 0.5 * cos.(2*pi * t), linewidth=4, color=:darkred)
scatter!(zeros(20), x[501, :], markersize=4, color=:orange)
scatter!(0.5 * ones(20), x[626, :], markersize=4, color=:blue)


t = range(-1, 1, length=1000)
R = 1/3 * cos.(2*pi * t) * cos.(2*pi * t')

heatmap(t, t, R, yticks= -1:0.25:1, xticks= -1:0.25:1, legend=false, yflip=true)

t = seq(-2, 2, len = 1000)
x = matrix(0, 1000, 20)

for (i in 1:20) {
    x[,i] = runif(1) * cos(2*pi*t)
}

matplot(t, x, lwd=2, col="gray", type="l", lty="solid", xaxp=c(-2,2, 8))
grid()
lines(t, 0.5*cos(2*pi*t), lwd=5, col="darkred")
points(numeric(20), x[501,], lwd=2, col="darkorange")
points(numeric(20) + 0.5, x[626,], lwd=2, col="blue", pch=4)


t = seq(-1, 1, len=1000)
R = (1/3)*outer(cos(2*pi*t), cos(2*pi*t))
image(t, t, R, col=topo.colors(255), xlab="t_1", ylab="t_2")

Example 10.12

% MATLAB code for Example 10.12
t = linspace(-1,1,1000);
R = toeplitz(0.5*cos(2*pi*t(:)));
imagesc(t,t,R);
grid on;
xticks(-1:0.25:1);
yticks(-1:0.25:1);

using ToeplitzMatrices
using Plots

x = zeros(1000, 20)
t = range(-2, 2, length=1000)

for col in eachcol(x)
    col .= cos.(2*pi*t .+ 2*pi*rand())
end

plot(t, x, linewidth=2, color=:gray80, legend=false)
plot!(t, 0 * cos.(2*pi * t), linewidth=4, color=:darkred)
scatter!(zeros(20), x[501, :], markersize=4, color=:orange)
scatter!(0.5 * ones(20), x[626, :], markersize=4, color=:blue)


t = range(-1, 1, length=1000)
R = SymmetricToeplitz(0.5 * cos.(2*pi * t))

heatmap(t, t, R, yticks= -1:0.25:1, xticks= -1:0.25:1, legend=false, yflip=true)

t = seq(-2, 2, len = 1000)
x = matrix(0, 1000, 20)

for (i in 1:20) {
    x[,i] = cos((2*pi*t) + (2*pi*runif(1)))
}

matplot(t, x, lwd=2, col="gray", type="l", lty="solid", xaxp=c(-2,2, 8))
grid()
lines(t, 0*cos(2*pi*t), lwd=5, col="darkred")
points(numeric(20), x[501,], lwd=2, col="darkorange")
points(numeric(20) + 0.5, x[626,], lwd=2, col="blue", pch=4)


t = seq(-1, 1, len=1000)
R = toeplitz(0.5*(cos(2*pi*t)))
image(t, t, R, col=topo.colors(255), ylim=c(1,-1), xlab="t_1", ylab="t_2")

Chapter 10.3 Wide sense stationary processes

Example 10.14

% MATLAB code to demonstrate autocorrelation
N     = 1000; % number of sample paths
T     = 1000; % number of time stamps
X     = 1*randn(N,T);
xc    = zeros(N,2*T-1);
for i=1:N
    xc(i,:) = xcorr(X(i,:))/T;
end
plot(xc(1,:),'b:', 'LineWidth', 2); hold on;
plot(xc(2,:),'k:', 'LineWidth', 2);

N = 1000
T = 1000
X = np.random.randn(N,T)
xc= np.zeros((N,2*T-1))
for i in range(N):
  xc[i,:] = np.correlate(X[i,:],X[i,:],mode='full')/T
plt.plot(xc[0,:],'b:')
plt.plot(xc[1,:],'k:')
plt.show()

using Plots
using DSP


X1 = randn(1000)
X2 = randn(1000)

plot(X1, linewidth=2, linestyle=:dot, color=:blue, legend=false)
plot!(X2, linewidth=2, linestyle=:dot, color=:black)


N  = 1000 
T  = 1000 
X  = randn(N, T)
xc = zeros(N, 2T - 1)

for i in 1:N
    xc[i, :] .= xcorr(X[i, :], X[i, :]) / T
end

plot(xc[1, :], linewidth=2, linestyle=:dot, color=:blue, label="correlation of sample 1")
plot!(xc[2, :], linewidth=2, linestyle=:dot, color=:black, label="correlation of sample 2")



Xa = rnorm(1000)
Xa2 = rnorm(1000)

plot(Xa, type="l", lwd=2, col="blue")
grid()
lines(Xa2, lwd=2)


N = 1000
T = 1000
X = matrix(rnorm(N*T), N, T)
xc = matrix(0, N, 2*T-1)

for (i in 1:N) {
    xc[i,] = ccf(X[i,], X[i,], lag.max=2*T-1, pl=FALSE)$acf/T
}

plot(xc[1,], type="l", lwd=2, col="darkblue")
lines(xc[2,], lwd=2)

Chapter 10.5 Wide sense stationary processes

Example 10.15

% MATLAB code for Example 10.15
t = -10:0.001:10;
L = length(t);
X = randn(1,L);
h = 10*max(1-abs(t),0)/1000;
Y = imfilter(X,h,'circular');

figure(1);
plot(t, X, 'LineWidth', 0.5, 'Color', [0.6 0.6 0.6]); hold on;
plot(t, zeros(1,L), 'LineWidth', 4, 'Color', [0.2 0.2 0.2]);
plot(t, Y, 'LineWidth', 2, 'Color', [0.85 0.3 0]);
plot(t, zeros(1,L), ':', 'LineWidth', 4, 'Color', [1 0.8 0]);
legend('X(t)','\mu_X(t)','Y(t)','\mu_Y(t)');
grid on;
axis([-10 10 -4 4]);
set(gcf, 'Position', [100, 100, 600, 400]);
set(gca,'FontWeight','bold','fontsize',14);

figure(2);
h2 = conv(h,h);
Rx = zeros(1,40001); Rx(20001) = 0.2;
plot(-20:0.001:20, Rx, 'LineWidth', 2, 'Color', [0.6 0.6 0.6]); hold on;
plot(-20:0.001:20, h2, 'LineWidth', 2, 'Color', [0.85 0.3 0]);
% axis([-2 2 -0.05 0.2]);
grid on;
legend('R_X(t)','R_Y(t)');
set(gcf, 'Position', [100, 100, 600, 400]);
set(gca,'FontWeight','bold','fontsize',14);

using Plots, LaTeXStrings
using Images
using DSP


t = -10:0.001:10
L = length(t)
X = randn(L)
h = 10 * max.(1 .- abs.(t), 0) / 1000
Y = imfilter(X, h, "circular")

plot(t, X, linewidth=0.5, color=:gray60, label=L"X(t)")
hline!([0], linewidth=4, color=:gray20, label=L"\mu_X(t)")
plot!(t, Y, linewidth=2, color=:red, label=L"Y(t)")
hline!([0], linewidth=4, linestyle=:dot, color=:yellow, label=L"\mu_Y(t)")
xlims!((-10, 10))


h2 = conv(h, h)
Rx = zeros(40001)
Rx[20001] = 0.2

plot(-20:0.001:20, Rx, linewidth=2, color=:gray60, label=L"R_X(t)")
plot!(-20:0.001:20, h2, linewidth=2, color=:red, label=L"R_Y(t)")
xlims!((-2, 2))
ylims!((-0.05, 0.2))


t = seq(-10, 10, by=0.001)
L = length(t)
X = rnorm(L)
h = 10 * sapply((1 - abs(t)), max, 0) / 1000
Y = convolve(X, h, type=c("circular"))


plot(t, X, lwd=1, col="gray", type="l")
grid()
legend(6, 3.8, legend=c("X(t)", "μ_x(t)", "Y(t)", "μ_y(t)"), col=c("gray", "black", "darkorange", "yellow"), lty=c(1, 1, 1, 3), lwd=c(1, 4, 3, 4), bg="white")
abline(h=0, lwd=4, lty=1, col="yellow")
lines(t, Y, lwd=3, col="darkorange")
abline(h=0, lwd=4, lty=3)


h2 = convolve(h, h, type="open")
Rx = numeric(40001)
Rx[20001] = 0.2
plot(seq(-20, 20, by=0.001), Rx, lwd=2, col="gray", type="l", xlim=c(-2,2), ylim=c(-0.05, 0.2))
grid()
legend(-2, 0.19, legend=c("R_x(t)", "R_y(t)"), col=c("gray", "darkorange"), lty=1:1, lwd=2)
lines(seq(-20, 20, by=0.001), h2, lwd=2, col="darkorange")

Chapter 10.6 Optimal linear filter

Solve the Yule Walker equation

Data download: ch10_LPC_data.txt (4KB)

% MATLAB code to solve the Yule Walker Equation
y = load('ch10_LPC_data.txt');
K = 10;
N = 320;
y_corr = xcorr(y);
R      = toeplitz(y_corr(N+[0:K-1]));
lhs    = y_corr(N+[1:K]);
h      = R\lhs;

y = np.loadtxt('./ch10_LPC_data.txt')
K = 10
N = 320
y_corr = np.correlate(y,y,mode='full')
R      = lin.toeplitz(y_corr[N-1:N+K-1]) 
lhs    = y_corr[N:N+K]
h      = np.linalg.lstsq(R,lhs,rcond = None)[0]

using Plots, LaTeXStrings
using ToeplitzMatrices
using DelimitedFiles
using DSP

y = vec(readdlm("ch10_LPC_data.txt"))
K = 10
N = length(y)
y_corr = xcorr(y, y)
R = SymmetricToeplitz(y_corr[N:N+K-1])
lhs = y_corr[N+1:N+K]
h = R \ lhs


plot(y, linewidth=4, color=:blue, label=L"Y[n]")


plot(y_corr, linewidth=4, color=:black, label=L"R_Y[k]")


y = scan("./ch10_LPC_data.txt")
K = 10
N = 320
y_corr = ccf(y, y, lag.max=2*length(y)-1, pl=FALSE)$acf
R = toeplitz(y_corr[N + 1:K-1])
lhs = y_corr[N + 1:K]
h = solve(R, lhs)


plot(y, lwd = 4, col="blue", type="l")
grid()
legend(10, 0.1, legend=c("Y[n]"), col=c("blue"), lty=1:1, lwd=4)


plot(y_corr, lwd = 4, type="l")
grid()
legend(10, 0.9, legend=c("R_y[k]"), col=c("black"), lty=1:1, lwd=4)

Predict sample

Data download: ch10_LPC_data_02.txt (4KB)

% MATLAB code to predict the samples
z    = y(311:320);
yhat = zeros(340,1);
yhat(1:320) = y;
for t = 1:20
    predict     = z'*h;
    z           = [z(2:10); predict];
    yhat(320+t) = predict;
end
plot(yhat, 'r', 'LineWidth', 3); hold on;
plot(y,    'k', 'LineWidth', 4);

z = y[310:320]
yhat = np.zeros((340,1))
yhat[0:320,0] = y
for t in range(20):
  predict = np.inner(np.reshape(z,(1,10)),h)
  z = np.concatenate((z[1:10], predict))
  yhat[320+t,0] = predict
plt.plot(yhat,'r')
plt.plot(y,'k')
plt.show()

using ToeplitzMatrices
using DelimitedFiles
using Plots
using DSP

y = vec(readdlm("ch10_LPC_data_02.txt"))
K = 10
N = length(y)

y_corr = xcorr(y, y)
R = SymmetricToeplitz(y_corr[N:N+K-1])
lhs = y_corr[N+1:N+K]
h = R \ lhs

z = y[311:320]
yhat = zeros(340)
yhat[1:320] .= y

for t in 1:20
    pred = z' * h
    z = [z[2:10]; pred]
    yhat[320+t] = pred
end

plot(yhat, linewidth=3, color=:red, label="Prediction", legend=:bottomleft)
plot!(y, linewidth=4, linestyle=:dash, color=:gray60, label="Input")


y = scan("./ch10_LPC_data_02.txt")
K = 10
N = length(y)

y_corr = ccf(y, y, lag.max=2*N-1)$acf[,,1]
R = toeplitz(y_corr[N + 1:K-1])
lhs = y_corr[N + 1:K]
h = solve(R, lhs)

z = y[311:320]
yhat = numeric(340)
yhat[1:320] = y

for (t in 1:20) {
    predict = t(z) * h
    z = c(z[2:10], predict)
    yhat[320+t] = predict
}

plot(yhat, lwd=3, col="red", type="l")
grid()
legend(10, 0.9, legend=c("Prediction", "Input"), col=c("red", "gray"), lty=c(1, 2), lwd=2)
lines(y, lwd=4, col="gray", lty=3)

Wiener filter

% MATLAB code for Wiener filtering
w = 0.05*randn(320,1);
x = y + w;

Ry = xcorr(y);
Rw = xcorr(w);
Sy = fft(Ry);
Sw = fft(Rw);
H = Sy./(Sy + Sw);
Yhat = H.*fft(x, 639);
yhat = real(ifft(Yhat));

plot(x, 'LineWidth', 4, 'Color', [0.7, 0.7, 0.7]); hold on;
plot(yhat(1:320), 'r', 'LineWidth', 2);
plot(y, 'k:', 'LineWidth', 2);

from scipy.fft import fft, ifft
w = 0.05*np.random.randn(320)
x = y + w

Ry = np.correlate(y,y,mode='full')
Rw = np.correlate(w,w,mode='full')
Sy = fft(Ry)
Sw = fft(Rw)
H  = Sy / (Sy+Sw)
Yhat = H * fft(x, 639)
yhat = np.real(ifft(Yhat))

plt.plot(x,color='gray')
plt.plot(yhat[0:320],'r')
plt.plot(y,'k:')

using DSP, FFTW
using Plots

w = 0.05 * randn(320)
x = y + w

Ry = xcorr(y, y)
Rw = xcorr(w, w)
Sy = fft(Ry)
Sw = fft(Rw)
H = Sy ./ (Sy + Sw)

Yhat = H .* fft_len(x, 639)
yhat = real.(ifft(Yhat))


plot(x, linewidth=4, color=:gray, label="Noise Input X[n]", legend=:bottomleft)
plot!(yhat[1:320], linewidth=2, color=:red, label="Wiener Filtered Yhat[n]")
plot!(y, linewidth=2, linestyle=:dot, color=:black, label="Ground Truth Y[n]")


plot(Rw, linewidth=4, color=:blue, label="h[n]")

library(pracma)
y = scan("./ch10_LPC_data_02.txt")
w = 0.05 * rnorm(320)
x = y + w

Ry = ccf(y, y, lag.max=2*length(y)-1, pl=FALSE)$acf
Rw = ccf(w, w, lag.max=2*length(w)-1, pl=FALSE)$acf
Sy = fft(Ry)
Sw = fft(Rw)
H = Sy / (Sy + Sw)

a = x[1:639]
a[is.na(a)] = 0
Yhat = H * fft(a[1:639])
yhat = Re(ifft(as.vector(Yhat)))


plot(x, lwd=5, col="gray", type="l")
grid()
legend(10, -0.7, legend=c("Noisy Input X[n]", "Wiener Filtered Yhat[n]", "Ground Truth Y[n]"), col=c("gray", "red", "black"), lty=c(1, 1, 2), lwd=2)
lines(yhat[1:320], lwd=2, col="red")
lines(y, lwd=2, lty=3)


plot(Rw, lwd=4, col="blue", label="h[n]", type="l")
legend(500, 0.85, legend=c("h[n]"), col=c("blue"), lty=c(1), lwd=c(4))
grid()

Wiener deblurring

Data download: ch10_wiener_deblur_data.txt (4KB)

% MATLAB code to solve the Wiener deconvolution problem
load('ch10_wiener_deblur_data');
g = ones(32,1)/32;
w = 0.02*randn(320,1);
x = conv(y,g,'same') + w;

Ry = xcorr(y);
Rw = xcorr(w);
Sy = fft(Ry);
Sw = fft(Rw);
G  = fft(g,639);

H = (conj(G).*Sy)./(abs(G).^2.*Sy + Sw);
Yhat = H.*fft(x, 639);
yhat = real(ifft(Yhat));

figure;
plot(x, 'LineWidth', 4, 'Color', [0.5, 0.5, 0.5]); hold on;
plot(16:320+15, yhat(1:320), 'r', 'LineWidth', 2);
plot(1:320, y, 'k:', 'LineWidth', 2);

y = np.loadtxt('./ch10_wiener_deblur_data.txt')
g = np.ones(64)/64
w = 0.02*np.random.randn(320)
x = np.convolve(y,g,mode='same') + w

Ry = np.correlate(y,y,mode='full')
Rw = np.correlate(w,w,mode='full')
Sy = fft(Ry)
Sw = fft(Rw)
G  = fft(g,639)
H = (np.conj(G)*Sy)/( np.power(np.abs(G),2)*Sy + Sw )

Yhat = H * fft(x, 639)
yhat = np.real(ifft(Yhat))

plt.plot(x,color='gray')
plt.plot(np.arange(32,320+32),yhat[0:320],'r')
plt.plot(y,'k:')

using DelimitedFiles
using DSP, FFTW
using Plots

y = vec(readdlm("ch10_wiener_deblur_data.txt"))
g = ones(32) / 32
w = 0.02 * randn(320)
x = conv_same(y, g) + w

Ry = xcorr(y, y)
Rw = xcorr(w, w)
Sy = fft(Ry)
Sw = fft(Rw)
G  = fft_len(g, 639)

H = @. (conj(G) * Sy) / (abs(G)^2 * Sy + Sw)
Yhat = H .* fft_len(x, 639)
yhat = real.(ifft(Yhat))

plot(x, linewidth=4, color=:gray, label="Noise Input X[n]")
plot!(16:320+15, yhat[1:320], linewidth=2, color=:red, label="Wiener Filtered Yhat[n]")
plot!(y, linewidth=2, linestyle=:dot, color=:black, label="Ground Truth Y[n]")

library(pracma)

conv_same = function(x, y) {
    s = length(y) / 2
    e = length(x) + s - 1
    return (convolve(x, y, type="open")[s:e])
}

y = scan("./ch10_wiener_deblur_data.txt")
g = ones(32, 1)/32
w = 0.02 * rnorm(320)
s = length(y)/2
e = length(g) + s - 1
x = conv_same(y, g) + w

Ry = ccf(y, y, lag.max=2*length(y)-1, pl=FALSE)$acf
Rw = ccf(w, w, lag.max=2*length(w)-1, pl=FALSE)$acf
Sy= fft(Ry)
Sw = fft(Rw)
a = g[1:639]
a[is.na(a)] = 0
G = fft(a[1:639])

H = (Conj(G) * Sy) / (abs(G) ^ 2 * Sy + Sw)
b = x[1:639]
b[is.na(b)] = 0
Yhat = H * fft(b[1:639])
yhat = Re(ifft(as.vector(Yhat)))

plot(x, lwd=4, col="gray", type="l", ylim=c(-0.7, 0.7))
grid()
legend(150, -0.45, legend=c("Noisy Input X[n]", "Wiener Filtered Yhat[n]", "Ground Truth Y[n]"), col=c("gray", "red", "black"), lty=c(1, 1, 2), lwd=2)
lines(16:(320+15), yhat[1:320], col="red", lwd=2)
lines(1:320, y, lwd=2, lty=3)

Auxiliary functions for Julia

function fft_len(x::AbstractVector, len::Int)
    if len > length(x)
        return fft([x; zeros(len - length(x))])
    else
        return fft(x[1:len])
    end
end

function conv_same(x::AbstractVector, y::AbstractVector)
    s = div(length(y), 2)
    e = length(x) + s - 1
    return conv(x, y)[s:e]
end