% Demonstration of Wiener filter for noise cancellation.
% Rich Kozick, Spring 1997

% Desired signal
s = auread('/home/pleiades/ACM/snd/clint_eastwood.au');
Ls = length(s);

% Interference + random noise
n = auread('/home/pleiades/ACM/snd/spacemusic.au');
Sign = 0.1;
n = n(1:Ls) + Sign*randn(Ls,1);

y = s + n;

Sigx = 0.01;
bx = [1];	% bx and ax are filtering on n to produce x
ax = [1 0.9];
Dx = 10;	% Delay of n that appears in x
x = filter(bx, ax, n);
x = [x(Dx:(Ls)); zeros(Dx-1,1)] + Sigx*randn(Ls,1);

N = 20;		% Length of FIR Wiener filter

% Compute Wiener filter coefficients

P = zeros(N,1);
R = zeros(N,N);

for k=N:Ls
  xk = x(k:-1:(k-N+1));
  P = P + y(k)*xk;
  R = R + xk*xk';
end

w = inv(R)*P;

nhat = filter(w, 1, x);

shat = y - nhat;
