dwiener4wlt1(sigwlt,nu,dn)

function w = dwiener4wlt1(sigwlt,nu,dn)

%----------------------------------------------------
% Compute the Wiener filter coefficients
%
% INPUT:
% sigwlt - sampled discrete signal wavelet (length n >= M filter
% length);
% nu - noise power ratio to signal power;
% dn - desired signal at n points;
%
% OUTPUT:
% w - Wiener filter in time domain;
%
% USAGE:
% w = dwiener4wlt1(sigwlt,nu,dn);
% Dependency: external "xcorr" and "toeplitz" functions
%


M=length(sigwlt);
%sigt=sigwlt+nswlt;

xc=xcorr(double(sigwlt), M-1);
%xcn=xcorr(nswlt, M-1);

%%[ys,is]=max(sigwlt);
%%[ydn,idn]=max(dn);
%%xc=xcorr([sigwlt;sigwlt], M-1); % symmetric reflection
%%xcn=xcorr([nswlt;nswlt], M-1); % symmetric reflection

R=toeplitz(xc(M:size(xc,1)));
Rn=sum(double(xc(M:size(xc,1))))*diag(ones(length(R),1),0);

R=R+nu*Rn;

p=xcorr(double(dn), double(sigwlt), M-1);

%%p=xcorr([dn;dn], [sigwlt;sigwlt], M-1); % symmetric reflection

if size(p,1)<size(p,2)
p=double(p(M:2*M-1)');
else p=double(p(M:2*M-1));
end;

%%%% Start Levinson Solution for general p-vector

%R=toep([1,.5,.2]) % example R
%p=-[-4,1,-3]' % example right hand side

%% Normalize equation by r0
p=double(p)/double(R(1,1)); R=double(R)/double(R(1,1));
p;

 rvec=double(R(2:end,1)); w=double(p)*0; 
%% Intial conditions
w(1)=double(p(1)); y=double(R(2:end,1));y(1)=-rvec(1); beta=1; alfa=-rvec(1);

for k=1:M-1
  beta=(1-double(alfa)^2)*double(beta);
  mu=(double(p(k+1))-double((rvec(1:k))'*w(k:-1:1)))/double(beta);
  
  v(1:k)=double(w(1:k))+double(mu)*double(y(k:-1:1)); 
  dimv=size(v); if dimv(2)>dimv(1) v=v'; end;
 
  w(1:k+1)=[double(v(1:k));double(mu)];
 
  if k<M-1
   alfa=-(double(rvec(k+1))+double((rvec(1:k))'*y(k:-1:1)))/double(beta);
   z(1:k)=double(y(1:k))+double(alfa*y(k:-1:1)); 
   dimz=size(z); if dimz(2)>dimz(1) z=z'; end;
 
   y(1:k+1)=[double(z(1:k));double(alfa)];
 
  end;  
end;

%iR=inv(R); %% conventional inversion

%w=iR*p; %% filter coefficients after inversion

%w=[w;w];
%sn=xcorr(w(size(w,1):-1:1), sigt, M-1);
%sn=xcorr(w,sigt,M-1);

%sn=sn(size(sn,1):-1:size(sn,1)-M+1);


return;