function[yw,w,en,xn]=aaWNC_with_leaking_signal(dn,alpha,a1,a2,v,M,N)
%
% Input: 
%   dn=desired signal;                              
%   alpha=fraction of signal d(n) that is observed mixed into the signal v_2(n) 
%   a1=first order IIR coefficient                                  
%   a2=first order IIR coefficient;                                 
%   v=noise;                                                        
%   M=number of Wiener filter coefficients;                         
%   N=number of sequence elemets of dn(desired signal) and v(noise);
% 
% Output: 
%   yw=estimated noise v_1(n)
%   w=optimal Wiener filter coefficients for estimating the noise;
%   en=the observed signal x(n) minus our estimate of v_1(n) i.e. our estimate of d(n) 
%   xn=corrupted and observed signal = d(n) + v_1(n)
%   

if M>N
    disp(['error:M must be less than N']);
end;

% compute the two versions of the noise
%
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 ... noised with with v_1(n)
yn=v2+alpha*dn;                           % <- this is the observed signal with the additional d(n) 

% Get R_y the autocovariance matrix of y(n):
%
yautoc=aasamplebiasedautoc(yn,M);         % <- we can get y(n) covariance by sampling 
Ry    =toeplitz(yautoc); 

% Get p_{v1,y} the cross-correlation of v1(n) and y(n)
%
pxy    = xcorr(xn,yn,'biased');             % <- this we can get from sampling ... 
pxy2   = pxy((length(pxy)+1)/2:(length(pxy)+1)/2+M-1);
dautoc = aasamplebiasedautoc(dn,M);         % <- this we would have to be obtained from earlier training ...  
p1     = pxy2 - alpha * dautoc;             

% Invert to get the optimal filter coefficients: 
%
w=(Ry)\(p1');                             % <- the optimal coefficients to estimate v_1(n) 
                                          %    from the observed y(n) 
                                          
yw=filter(w,1,yn);                        % <- estimate v_1(n) using our MA(M) filter coeffients 

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