## Contents

clear all


## Experiment setup for Wiener filter example

N = 100000;
d = zeros(N,1);
x = zeros(N,1);


## Generate independent random variable v_1

Generate N i.i.d (independent and identically distributed) random variables with distribution N(0,0.27)

v1= sqrt(0.27)*randn(N,1);


## Generate the variable d

a = 0.8458;
d(1) = 0;
for n= 2:N
d(n) = -a*d(n-1) +v1(n);
end


## Generate the variable x

b = -0.9458;
x(1) = 0;
for n= 2:N
x(n) = -b*x(n-1) +d(n);
end


## Generate the variable u

v2 = sqrt(0.1)*randn(N,1);
u = x(:) + v2(:);


## Compute correlations

k = (N/4):N;  % interval over which to compute the correlations
for tau = 0:20
r_u1(tau+1) = sum(u(k).*u(k-tau))/length(k);
end
r_u1
for tau = 0:20
r_ud1(tau+1) = sum(d(k).*u(k-tau))/length(k);
end
r_ud1

r_u1 =

Columns 1 through 7

1.0941    0.4903    0.8442    0.4759    0.7228    0.4522    0.6234

Columns 8 through 14

0.4246    0.5408    0.3966    0.4707    0.3656    0.4133    0.3364

Columns 15 through 21

0.3660    0.3087    0.3235    0.2823    0.2894    0.2551    0.2587

r_ud1 =

Columns 1 through 7

0.5300   -0.4499    0.3798   -0.3220    0.2725   -0.2303    0.1947

Columns 8 through 14

-0.1645    0.1385   -0.1152    0.0962   -0.0812    0.0684   -0.0552

Columns 15 through 21

0.0483   -0.0384    0.0331   -0.0258    0.0233   -0.0194    0.0168



## Compute the correlation values from the Yule_walker equation

See the matrix form at the bottom of page 19

a1 = -0.1;
a2 = -0.8;
sigma2_v = 0.27;
sigma2_v2 = 0.1;
B=  [1 a1 a2
a1 1+a2 0
a2 a1 1];
C = [sigma2_v; 0; 0];
r_x = B\C

r_u(1:2,1) = r_x(1:2,1) + [sigma2_v2; 0];

p(1,1) = r_x(1) + b*r_x(2);
p(2,1) = r_x(2) + b*r_x(1)

r_x =

1.0000
0.5000
0.8500

p =

0.5271
-0.4458



## Wiener_Hopf design

w = [r_u(1) r_u(2); r_u(2) r_u(1)]\[p(1); p(2)]

w =

0.8362
-0.7853



## Estimate by the Yule_Walker equation 20 values of r_x

This is just one extension of YW to compute more values of r_x

B=  [1 a1 a2 zeros(1,17)
a1 1+a2 0 zeros(1,17)
a2 a1 1 zeros(1,17)
0 a2 a1 1 zeros(1,16)
zeros(1,2) a2 a1 1 zeros(1,15)];
for ii = 6:20
B = [B; zeros(1,ii-3) a2 a1 1 zeros(1,20-ii)];
end
C = [sigma2_v; 0; 0; zeros(17,1)];
r_x_long = B\C;


## Plot the waveforms

figure(1),plot(v1(1:200))
title('v1 are normal distributed i.i.d samples')
xlabel('index k')
ylabel('the value u(k)')
grid
figure(2),plot(d(1:200))
title('d is a first order AR(1) process')
xlabel('index k')
ylabel('the value d(k)')
grid
figure(3),plot(x(1:200))
title('x is a second order AR(2) process')
xlabel('index k')
ylabel('the value x(k)')
grid
figure(4),plot(u(1:200))
title('u is a second order AR(2) process plus noise')
xlabel('index k')
ylabel('the value u(k)')
grid
figure(5),plot(0:20, r_u1,'or-',0:19, r_x_long,'ob-')
title('r_u is the correlation of a second order AR(2) process plus noise')
xlabel('index \tau')
ylabel('the value r_u(\tau)')
legend('r_u(\tau)','r_x(\tau)')
grid