## 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.4819    0.8352

r_ud1 =

0.5308   -0.4509

```

## 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.8372
-0.7863

w_WHm =

0.8362
-0.7853

```

## 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.8342    0.8372    0.8362
-0.7834   -0.7863   -0.7853

```

## 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

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