Contents

clear all

Repeat the 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 the needed correlations

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

       1.0865      0.48187      0.83517


r_ud1 =

      0.53078     -0.45089

Wiener_Hopf design using estimated correlations from data

Estimated correlation matrix

R1 = [r_u1(1) r_u1(2)
    r_u1(2) r_u1(1)];
% Estimated cross-correlation vector
pd = r_ud1(:);
w_WHd = R1\pd

R_u=[1.1 0.5 ; % Define autocovariance matrix
0.5 1.1 ];
p= [0.5271 ; -0.4458]; % Define cross-correlation vector
w_WHm = R_u\p
w_WHd =

      0.83724
      -0.7863


w_WHm =

      0.83616
     -0.78534

Steepest descent algorithm

R_u=[1.1 0.5 ; % Define autocovariance matrix
0.5 1.1 ];
p= [0.5271 ; -0.4458]; % Define cross-correlation vector
W=[-1;-1];
mu= 0.01; % Initial values for the adaptation algorithm
%mu= 1;
Wt=W; % Wt will record the evolution of vector W
for k=1:1000 % Iterate 1000 times the adaptation step
W=W +mu*(p-R_u*W); % Adaptation Equation ! Quite simple!
Wt=[Wt W]; % Wt records the evolution of vector W
end % Is W the Wiener Filter?

[Wt(:,end) w_WHd  w_WHm]
ans =

      0.83418      0.83724      0.83616
     -0.78337      -0.7863     -0.78534

Plot the trajectory of the algorithm

figure(1),clf
plot(Wt(1,:),Wt(2,:),'-b.',w_WHm(1),w_WHm(2),'ro')
hold on
grid

% Plot the curves where the quadratic form J(w) = sigma_d^2 -2*p'*w + w'*R*w = const
% or equivalently, where (w-w_o)'*R*(w-w_o) = const

L = chol(R_u,'lower');
theta = 0:0.0001:2*pi;
for rho = 0.1:0.1:2
wconst = w_WHm*ones(1,length(theta)) + rho*inv(L)*[cos(theta); sin(theta)];
plot(wconst(1,:),wconst(2,:),'b-')
end
axis equal

figure(2)