function stsa_wlr(filename,outfile)

%
%  Implements the Bayesian estimator based on the weighted likelihood ratio
%  distortion measure [1, Eq. 37].
% 
%  Usage:  stsa_wlr(noisyFile, outputFile)
%           
%         infile - noisy speech file in .wav format
%         outputFile - enhanced output file in .wav format
%  
%
%  Example call:  stsa_wlr('sp04_babble_sn10.wav','out_wlr.wav');
%
%  References:
%   [1] Loizou, P. (2005). Speech enhancement based on perceptually motivated 
%       Bayesian estimators of the speech magnitude spectrum. IEEE Trans. on Speech 
%       and Audio Processing, 13(5), 857-869.
%   
% Author: Philipos C. Loizou
%
% Copyright (c) 2006 by Philipos C. Loizou
% $Revision: 0.0 $  $Date: 10/09/2006 $
%-------------------------------------------------------------------------

if nargin<2
   fprintf('Usage: stsa_wlr inFile outFile.wav \n\n');
   return;
end

[x, Srate, bits]= wavread( filename);	


% =============== Initialize variables ===============
%
len=floor(20*Srate/1000); % Frame size in samples
if rem(len,2)==1, len=len+1; end;
PERC=50; % window overlap in percent of frame size
len1=floor(len*PERC/100);
len2=len-len1; 


win=hanning(len); %tukey(len,PERC);  % define window


% Noise magnitude calculations - assuming that the first 6 frames is noise/silence
%
nFFT=len;
nFFT2=len/2;
noise_mean=zeros(nFFT,1);
j=1;
for k=1:5
   noise_mean=noise_mean+abs(fft(win.*x(j:j+len-1),nFFT));
   j=j+len;
end
noise_mu=noise_mean/5;
noise_mu2=noise_mu.^2;

%--- allocate memory and initialize various variables
   
img=sqrt(-1);
x_old=zeros(len1,1);
Nframes=floor(length(x)/len2)-1;
xfinal=zeros(Nframes*len2,1);
xinterv=0.001:0.01:10;
k=1;
aa=0.98;

%===============================  Start Processing =======================================================
%
fprintf('This might take some time ...\n')
for n=1:Nframes 
   
  
   insign=win.*x(k:k+len-1);
    
   %--- Take fourier transform of  frame
   
   spec=fft(insign,nFFT);   
   sig=abs(spec); % compute the magnitude
   sig2=sig.^2;
   
       gammak=min(sig2./noise_mu2,40);  % post SNR. Limit it to avoid overflows
              
       if n==1
           ksi=aa+(1-aa)*max(gammak-1,0);
       else
           ksi=aa*Xk_prev./noise_mu2 + (1-aa)*max(gammak-1,0);     % a priori SNR   
       end
     
       vk=ksi.*gammak./(1+ksi);
      
       
       xx=solve_wlr(vk,gammak,sig,xinterv); % solves Eq. 37 in [1]
       
       sig_hat=xx;
       Xk_prev=sig_hat.^2;
       
       xi_w= ifft( sig_hat.* exp(img*angle(spec))); 
	   xi_w= real( xi_w);
	  
      
	% --- Overlap and add ---------------
    %
    xfinal(k:k+ len2-1)= x_old+ xi_w(1:len1);
	x_old= xi_w(len1+ 1: len);

    if rem(n,20)==0, fprintf('Frame: %d Percent completed:%4.2f \n',n,n*100/Nframes); end;
    
 k=k+len2;
end
%========================================================================================




wavwrite(xfinal,Srate,16,outfile);


%==========================================================================
function x=solve_wlr(vk,gammak,Yk,xx);

% solves non-linear Eq. 37 in [1]
%

Len=length(vk);
L2=Len/2+1;

lk05=sqrt(vk).*Yk./gammak;
Ex=gamma(1.5)*lk05.*confhyperg(-0.5,1,-vk,100);
Elogx=1-0.5*(2*log(lk05)+log(vk)+expint(vk));

x=zeros(Len,1);
  
for n=1:L2

    a=Elogx(n);
    b=Ex(n);
    ff=sprintf('log(x)+%f - %f/x',a,b);
    y=log(xx)+a-b./xx;
    bet=xx(1); tox=200;
    if y(1)<0
            ind=find(y>0);
            bet=xx(1)/2;
            tox=xx(ind(1));
             
            [x(n),fval,flag]=fzero(inline(ff),[bet tox]);
            if flag<0
                x(n)=x(n-1);
           end
    else
        ind=find(y<0);
        if ~isempty(ind)
            bet=xx(1);
            tox=xx(ind(1));
            [x(n),fval]=fzero(inline(ff),[bet tox]);
            
        else
   
            x(n)=0.001;  % spectral floor
   
        end
    end    
       

end

x(L2+1:Len)=flipud(x(2:L2-1));