function[yw,w,en,xn]=aawienernoisecancelor(dn,a1,a2,v,M,N)
%  
% Input: 
%   dn=desired signal;                                              
%   a1=first order IIR coefficient for v_1(n)                                 
%   a2=first order IIR coefficient for v_2(n) 
%   v=noise;                                                        
%   M=number of Wiener filter coefficients;                         
%   N=number of sequence elemets of dn(desired signal) and v(noise);
% 
% Output: 
%   yn=
%   w=Wiener filter coefficients;      
%   en=our estimated signal \hat{d}(n) time series after removing an estimate of v_1(n) from x(n)
%   xn=corrupted signal = d(n) + v_1(n);               
% 
v1(1)=0;v2(1)=0;
for n=2:N
    v1(n)=a1*v1(n-1)+v(n-1);
    v2(n)=a2*v2(n-1)+v(n-1);
end;
xn=dn+v1;                                 % <- the noised signal ... with v_1(n)

v2autoc=aasamplebiasedautoc(v2,M);        % <- we observe v_2(n) however ...
Rv2 = toeplitz(v2autoc);
R   = Rv2;

p1=xcorr(xn,v2,'biased');
%p1=xcorr(v1,v2,'biased');
if M>N
    disp(['error:M must be less than N']);
end;
p=p1((length(p1)+1)/2:(length(p1)+1)/2+M-1);
%w=inv(R)*p';
w=R\(p');                                 % <- the optimal coefficients to estimate v_1(n) 
                                          %    from the observed v_2(n) 
yw=filter(w,1,v2);                        % <- estimate v_1(n) using our MA(M) filter coeffients 

en=xn-yw;                                 % <- remove \hat{v}_1(n) from x(n) to get \hat{d}(n)