function outfile= mt_mask( noisy_file, outfile) 

%
%  Implements a psychoacoustically motivated algorithm [1].
% 
%  Usage:  mt_mask(noisyFile, outputFile)
%           
%         infile - noisy speech file in .wav format
%         outputFile - enhanced output file in .wav format
%
%         
%  Example call:  mt_mask('sp04_babble_sn10.wav','out_mask.wav');
%
%  References:
%   [1] Hu, Y. and Loizou, P. (2004). Incorporating a psychoacoustical model in 
%   	frequency domain speech enhancement. IEEE Signal Processing Letters, 11(2), 
%       270-273.
%   
% Authors: Yi Hu and Philipos C. Loizou
%
% Copyright (c) 2006 by Philipos C. Loizou
% $Revision: 0.0 $  $Date: 10/09/2006 $
%-------------------------------------------------------------------------

if nargin<2
   fprintf('Usage: mt_mask(noisyfile.wav,outFile.wav) \n\n');
   return;
end



% Initialize wavelet parameters (see also wiener_wt.m)
wavname='db4';
thre_type='ds';thre_func_type='s';q_0=5;
taper_num=16;

%------------------get the noisy speech data
[noisy_speech, Srate, NBITS]= wavread( noisy_file);

%===========initiate the parameters=======================
frame_dur= 20; %unit is milli-second
len= floor( Srate* frame_dur/ 1000);    
if rem( len, 2)~= 0
	len= len+ 1;
end
NFFT= len;  %number of FFT points
tapers= sine_taper( taper_num, NFFT);
diga= digamma( taper_num)- log( taper_num);

win= hamming( len);
% win= win/ norm( win);
PERC= 50; % window overlap in percent of frame size
len1=floor(len* PERC/ 100);
len2= len- len1; 
L120= floor( 120* Srate/ 1000);
bfl=0.002; % spectral floor


k= 1;   %k is starting point of each frame

%================================================

q= ceil( log2( len)); 
M= 2^ q;

sigma_eta_square= trigamma( taper_num);
N_autoc= sigma_eta_square* ( 1- ( 0: taper_num+ 1)/ ( taper_num+ 1));
N_autoc( M/ 2+ 1)= 0;
Sigma_N_firstrow= [N_autoc( 1: M/ 2+ 1), fliplr( N_autoc( 2: M/ 2))];
noise_stat= real( fft( Sigma_N_firstrow));

[wfilter( 1, :), wfilter( 2, :), wfilter( 3, :), wfilter( 4, :)]= ...
	wfilters( wavname);
%------get the wavelet/scaling filter for decomposition/reconstruction

noise= noisy_speech( 1: L120);
noise_ps= psd_mt_sine( noise, tapers);
log_noise_ps= log( noise_ps)- diga;
den_log_noise_ps= thre_wavelet( log_noise_ps, noise_stat, thre_type, ...
    thre_func_type, wfilter, q_0);
den_log_noise_ps= [den_log_noise_ps( 1: len/ 2+ 1); ...
    flipud( den_log_noise_ps( 2: len/ 2))];
noise_ps= exp( den_log_noise_ps);
%=================

mu_vad= 0.98; % smoothing factor in noise spectrum update
aa= 0.98; % smoothing factor in priori update
eta= 0.15; % VAD threshold

%=================

Nframes= floor( length( noisy_speech)/ len2)- 1;
x_old= zeros( len1, 1);
xfinal= zeros( Nframes* len2, 1);

%===============================  Start Processing ==========

for n= 1: Nframes 
	
	insign= noisy_speech( k: k+ len- 1);    
	insign_spec= fft( insign.* win, NFFT);
	
	%========estimate the noisy speech power spectrum
	ns_ps= psd_mt_sine( insign, tapers);    
	
    log_ns_ps= log( ns_ps)- diga;
	den_log_ns_ps= thre_wavelet( log_ns_ps, noise_stat, thre_type, ...
		thre_func_type, wfilter, q_0);
	den_log_ns_ps= [den_log_ns_ps( 1: NFFT/ 2+ 1); ...
			flipud( den_log_ns_ps( 2: NFFT/ 2))];
	ns_ps= exp( den_log_ns_ps); 
	%=================================================	   
    
    gammak= abs( insign_spec).^ 2/ (norm( win)^2)./ noise_ps;
    if n==1
        ksi=aa+(1-aa)*max(gammak-1,0);
    else
        ksi=aa*Xk_prev./noise_ps + (1-aa)*max(gammak-1,0);  
        % decision-direct estimate of a priori SNR
    end

    log_sigma_k= gammak.* ksi./ (1+ ksi)- log(1+ ksi);    
    vad_decision(n)= sum( log_sigma_k)/ len;    
    if (vad_decision(n)< eta) 
        % noise only frame found
        noise_ps= mu_vad* noise_ps+ (1- mu_vad)* ns_ps;
        vad( k: k+ len- 1)= 0;
    else
        vad( k: k+ len- 1)= 1;
    end
    % ===end of vad===
	
	%========estimate the clean speech power spectrum
	cl_ps= ns_ps- noise_ps; 
	cl_ps= max( cl_ps, bfl* ns_ps);
	%--providing a spectral floor
	%========
	
	%compute the masking threshold 
	mask_thre= mask( cl_ps( 1: NFFT/ 2+ 1), NFFT, Srate, 16);
	mask_thre= [mask_thre; flipud( mask_thre( 2: NFFT/ 2))];
	%expand it to NFFT length
	
	noise_mask_ratio= noise_ps./ mask_thre;
	%=======two methods to compute g_wi
	%   get the mu_k by u= max( sqrt( Sn/ alpha- 1), 0) * Sx/ Sn
	%aprioSNR= cl_ps./ noise_ps;
	%mu( :, n)= max( sqrt( noise_mask_ratio)-1, 0).* aprioSNR;
	%g_wi= aprioSNR./ ( aprioSNR+ mu_n);
    tmp= max( sqrt( noise_mask_ratio)-1, 0);        
	g_wi= 1./ (1+ tmp);
    
    xi_freq= g_wi.* insign_spec;
    Xk_prev= abs( xi_freq).^ 2; 
	
	xi_w= ifft( xi_freq); 
	xi_w= real( xi_w);
	
	xfinal( k: k+ len2- 1)= x_old+ xi_w( 1: len1);
	x_old= xi_w( len1+ 1: len);
	k= k+ len2;
	
end
%========================================================================================

wavwrite( xfinal, Srate, 16, outfile);
 
    
%========================================================================================

function after_thre= thre_wavelet( before_thre, noise_stat, ...
	thre_type, thre_func_type, wfilter, q_0)

%this function implements the wavelet thresholding technique
%   refer to the paper by Walden/1998, Donoho/1995, Johnstone/1997

%note on the parameters
%   before_thre: data before thresholding
%   noise_stat: the power spectrum of the noise (i.e., noise statistics), 
%       DFT of the first row of Sigma_N, refer to Eq. (8) in Walden's paper
%   thre_type: threshold type, scale-dependent Universal ('d'), 
%		scale-independent Universal ('i'), scale-dependent SURE ('ds'), 
%		scale-independent SURE ('is'), or scale-dependent Generalized 
%		Corss-Validation ('dg')
%   thre_func_type: threshold function type: soft ('s') or hard ('h');
%   wfilter: wavelet low pass and high pass decomposition/reconstruction filters [lo_d, hi_d, lo_r, hi_r]
%       the 1st row is lo_d, the 2nd row is hi_d, the 3rd row is lo_r, and the 4th row is hi_r
%   q_0 is the decomposition level

%   after_thre: data after thresholding

s= size( before_thre);  
before_thre= before_thre( :)';   %make it a row vector
noise_stat= noise_stat( :)';

N= length( before_thre);    %length of before-thresholded data
q= ceil( log2( N)); 
M= 2^ q;

%==get the low pass and high pass decomposition/reconstruction filters from wfilter
lo_d= wfilter( 1, :);   %low pass decomposition filter/ scaling filter
hi_d= wfilter( 2, :);   %high pass decomposition filter/ wavelet filter
lo_r= wfilter( 3, :);   %low pass reconstruction filter/ scaling filter
hi_r= wfilter( 4, :);   %high pass reconstruction filter/ wavelet filter

%==refer to pp. 3155 in Walden's paper
H= zeros( q_0, M);
H( 1, :)= fft( hi_d, M);   %frequency response of wavelet filter
G( 1, :)= fft( lo_d, M);  %frequency response of scaling filter
for i= 2: q_0- 1
    G( i, :)= G( 1, rem( (2^ (i- 1) )* (0: M- 1), M)+ 1);
end

for j= 2: q_0
    H( j, :)= prod( [G( 1: j- 1, :); H( 1, rem( (2^ (j- 1) )* (0: M- 1), M)+ 1)], 1);
end

[y_coeff, len_info]= wavedec( before_thre, q_0, lo_d, hi_d);

% --decompose before_thre into q_0 levels using wavelet filter hi_d and scaling filter lo_d
% --where y_coeff contains the coefficients and len_info contains the length information
% --different segments of y_coeff correspond approximation and detail coefficients;
% -- length of len_info should be q_0+ 2

%===============processing according to 'thre_type'
%-------with 'd'--scale-dependent thresholding, threshold has to be computed for each level
%-------with 'i'--scale-independent thresholding, threshold is set to a fixed level

if thre_type== 'i' %scale-independent universal thresholding
    sigma_square= mean( noise_stat);
    thre= sqrt( sigma_square* 2* log( M)) ;    %mean( noise_stat) is sigma_eta_square in Eq. (6)
    y_coeff( len_info( 1)+ 1: end)= ...
        wthresh( y_coeff( len_info( 1)+ 1: end), thre_func_type, thre);
    
elseif thre_type== 'd' %scale-dependent universal thresholding
    %------first we need to compute the energy level of each scale from j= 1: q_0
    for i= 1: q_0   %refer to Eq. (9) in Walden's paper
        sigma_j_square( i)= mean( noise_stat.* (abs( H( i, :)).^ 2), 2);   %average along the row          
    end

    for i= 2: q_0+ 1    %thresholding for each scale
        
        sp= sum( len_info( 1: i- 1), 2)+ 1; %starting point
        ep= sp+ len_info( i)- 1;
        thre= sqrt( sigma_j_square( q_0- i+ 2)* 2* log( len_info( i)));
        y_coeff( sp: ep)= wthresh( y_coeff( sp: ep), thre_func_type, thre);
        
    end
    
elseif thre_type== 'ds' %scale-dependent SURE thresholding
    
%=======use Eq. (9) in Walden's paper to get sigma_j, MDA estimate seems to be better
%     for i= 1: q_0   
%         sigma_j_square( i)= mean( noise_stat.* (abs( H( i, :)).^ 2), 2);   %average along the row 
%         sigma_j( i)= sqrt( sigma_j_square( i));
%     end

%======MDA estimate of sigma_j
    sigma_j= wnoisest( y_coeff, len_info, 1: q_0);  
    
    for i= 2: q_0+ 1    %thresholding for each scale
        
        sp= sum( len_info( 1: i- 1), 2)+ 1; %starting point
        ep= sp+ len_info( i)- 1;    %ending point
        if sigma_j( q_0- i+ 2)< sqrt( eps)* max(  y_coeff( sp: ep));
            thre= 0;
        else
            thre= sigma_j( q_0- i+ 2)* thselect( y_coeff( sp: ep)/ ...
                sigma_j( q_0- i+ 2), 'heursure');
        end
        
        %fprintf( 1, 'sigma_j is %6.2f, thre is %6.2f\n', sigma_j, thre);
        y_coeff( sp: ep)= wthresh( y_coeff( sp: ep), thre_func_type, thre);
        
    end
    
elseif thre_type== 'dn' %new risk function defined in Xiao-ping Zhang's paper
    
    sigma_j= wnoisest( y_coeff, len_info, 1: q_0);  
    sigma_j_square= sigma_j.^ 2;
    
    for i= 2: q_0+ 1    %thresholding for each scale
        
        sp= sum( len_info( 1: i- 1), 2)+ 1; %starting point
        ep= sp+ len_info( i)- 1;    %ending point       
        if sigma_j( q_0- i+ 2)< sqrt( eps)* max(  y_coeff( sp: ep));
            thre= 0;
        else    
            
            %based on some evidece, the following theme let thre vary with SNR
            %   with ultra low SNR indicating low probability of signal presence, 
            %       hence using universal threshold
            %   and very high SNR indicates high probability of signal presence,
            %       hence using SURE threshold
            
            thre_max= sigma_j( q_0- i+ 2)* sqrt( 2* log( len_info( i))); %thre with SNRlog< -5dB
            thre_min= sigma_j( q_0- i+ 2)* fminbnd( @riskfunc, 0, sqrt(2* log( ep- sp+ 1)), ...
                optimset( 'MaxFunEvals',1000,'MaxIter',1000), ...
                y_coeff( sp: ep)/ sigma_j( q_0- i+ 2), 3);   %thre with SNRlog> 20dB
            slope= (thre_max- thre_min)/ 25;
            thre_0= thre_min+ 20* slope;        

            SNRlog= 10* log10( mean( max( y_coeff( sp: ep).^ 2/ sigma_j_square( q_0- i+ 2)- 1, 0)));            
            if SNRlog>= 20
                thre= thre_min;  %actually this corresponds to SURE threshold
            elseif ( SNRlog< 20) & ( SNRlog>= -5)
                thre= thre_0- SNRlog* slope;
            else
                thre= thre_max;   %this corresponds to oversmooth threshold
            end
            
            %the theme below is similar to the option 'heursure' in the function 'thselect'
%             univ_thr = sqrt(2* log( len_info( i)));  %universal thresholding
%             eta = (norm( y_coeff( sp: ep)/ sigma_j( q_0- i+ 2)).^2)/ ( len_info( i))- 1;
%             crit = (log2( len_info( i)))^(1.5)/ sqrt( len_info( i));
%             if 1%eta > crit   %high probility that speech exists
%                 thre= sigma_j( q_0- i+ 2)* fminbnd( @riskfunc, 0, sqrt(2* log( ep- sp+ 1)), ...
%                     optimset( 'MaxFunEvals',1000,'MaxIter',1000), ...
%                     y_coeff( sp: ep)/ sigma_j( q_0- i+ 2), 3); 
%             else
%                 thre = sigma_j( q_0- i+ 2)* univ_thr;                    
%             end

        end
        
        y_coeff( sp: ep)= wthresh( y_coeff( sp: ep), thre_func_type, thre);
        
    end

elseif thre_type== 'dg' %scale-dependent Generalized Cross Validation thresholding
    
    for i= 2: q_0+ 1    %thresholding for each scale
        
        sp= sum( len_info( 1: i- 1), 2)+ 1; %starting point
        ep= sp+ len_info( i)- 1;    %ending point       
        [y_coeff( sp: ep), thre]= mingcv( y_coeff( sp: ep), thre_func_type);
        
    end   
   
else 
    error( 'wrong thresholding type');
end

%--reconstruct the thresholded coefficients
after_thre= waverec( y_coeff, len_info, lo_r, hi_r);

if s(1)>1 
    after_thre= after_thre'; 
end
%fprintf( 1, 'thre is %f\n', thre);



function mt_psd= psd_mt_sine( data, sine_tapers)

% this function uses sine tapers to get multitaper power spectrum estimation
% 'x' is the incoming data, 'sine_tapers' is a matrix with each column being
% sine taper, sine_tapers can be obtained using the function sine_taper

[frame_len, taper_num]= size( sine_tapers);

eigen_spectra= zeros( frame_len, taper_num);

data= data( :);
data_len= length( data);
data_hankel= hankel( data( 1: frame_len), data( frame_len: data_len));

x_mt_psd= zeros( frame_len, data_len- frame_len+ 1);

for pp= 1: data_len- frame_len+ 1
    for index= 1: taper_num
        x_taperd= sine_tapers( :, index).* data_hankel( :, pp);
        x_taperd_spec= fft( x_taperd);
        eigen_spectra( :, index)= abs( x_taperd_spec).^ 2;
    end
    x_mt_psd(:, pp)= mean( eigen_spectra, 2);
end

mt_psd= mean( x_mt_psd, 2);



function tapers= sine_taper( L, N)

% this function is used to generate the sine tapers proposed by Riedel et
% al in IEEE Transactions on Signal Processing, pp. 188- 195, Jan. 1995

% there are two parameters, 'L' is the number of the sine tapers generated,
% and 'N' is the length of each sine taper; the returned value 'tapers' is
% a N-by-L matrix with each column being sine taper

tapers= zeros( N, L);

for index= 1: L
    tapers( :, index)= sqrt( 2/ (N+ 1))* sin (pi* index* (1: N)'/ (N+ 1));
end





function y = trigamma(z,method,debug)

%  y = trigamma(z)   ... Trigamma-Function for real positive z
%
%  trigamma(z) = (d/dz)^2 log(gamma(z)) = d/dz digamma(z)
%
%  if 'z' is a matrix, then the digamma-function is evaluated for
%  each element. Results are inaccurate for real arguments < 10 which are
%  neither integers nor half-integers.
%
%  y = trigamma(z,method)
%
%  possible values for optional argument 'method':
%    method = 1           : quick asymptotic series expansion (approximate)
%    method = 2           : finite recursion for integer values (exact)
%    method = 3           : finite recursion for half-integer values (exact)
%    method = 4 (default) : automatic selection of 1,2 or 3 for individual
%                           elements in z whichever is appropriate.
%
%  see also: digamma, gamma, gammaln, gammainc, specfun


%  reference: Abramowitz & Stegun, "Handbook of Mathematical Functions"
%             Chapter "Gamma Function and Related Functions" :
%  implemented by: Christoph Mecklenbraeuker
%             (email: cfm@sth.ruhr-uni-bochum.de), July 4, 1995.


dim = size(z); 				% save original matrix dimension
z = reshape(z,dim(1)*dim(2),1); 	% make a column vector
I1 = ones(length(z),1); 		% auxiliary vector of ones

if(nargin==1)
    method=4; debug=0;
elseif(nargin==2)
    debug=0;
end;


if(debug == 1) 				% if debug==1: track recursion
    [m,n] =size(z);
    fprintf(1,'trigamma: method =	%d, size(z)=[%d %d],\t min(z)=%f, max(z)=%f\n',...
        method,m,n,min(min(z)),max(max(z)));
end;

if(method==1) 				% use 9th order asymptotic expansion
    if(any(z<1))
        fprintf(1,'Warning: some elements in argument of "trigamma(z,1)" are < 1\n');
        fprintf(1,'minimal argument = %g: trigamma-result is inaccurate!\n',min(min(z)));
    end

    % calculate powers of 1/z :
    w1 = 1./z; w2 = w1.*w1; w3 = w1.*w2; w5 = w2.*w3; w7 = w2.*w5; w9 = w2.*w7;
    % generate coefficients of expansion: matrix with constant columns
    a = [ I1   I1/2   I1/6  -I1/30  I1/42 -I1/30];
    % make vector of powers of 1/z:
    w = [ w1   w2     w3     w5      w7    w9];
    % calculate expansion by summing the ROWS of (a .* w) :
    y = sum((a.*w).').';
elseif(method==2)
    zmax = max(max(floor(z)));
    ytab = zeros(zmax,1);
    ytab(1) = pi^2/6; 			% = psi'(1)
    for n=1:zmax-1;
        ytab(n+1) = ytab(n) - 1/n^2; 	% generate lookup table
    end;
    y = ytab(z);
elseif(method==3)
    zmax = max(max(floor(z)));
    ytab = zeros(zmax+1,1);
    ytab(1) = pi^2/2; 			% = psi'(1/2)
    for n=1:zmax;
        ytab(n+1) = ytab(n) - 4/(2*n-1)^2; % generate lookup table
    end;
    y = ytab(z+0.5);
elseif(method==4) 			% decide here which method to use
    Less0 = find(z<0); 			% negative arguments evaluated by reflexion formula
    Less1 = find(z>0 & z<1); 		% values between 0 and 1.
    fraction = rem(z,1); 			% fractional part of arguments
    f2 = rem(2*fraction,1);
    Integers = find(fraction==0 & z>0); 	% Index set of positive integer arguments
    NegInts  = find(fraction==0 & z<=0); 	% Index set of positive integer arguments
    HalfInts = find(abs(fraction-0.5)<1e-7 & z>0); % Index set of positive half-integers
    Reals    = find(f2>1e-7 & z>1); 	% Index set of all other arguments > 1
    if(~isempty(Reals)) y(Reals)    = trigamma(z(Reals),1,debug);    end;
    if(~isempty(Less1)) y(Less1)    = trigamma(z(Less1)+2,1,debug) + ...
            1./z(Less1).^2+1./(z(Less1)+1).^2;end;
    % reflexion formula:
    if(~isempty(Less0)) y(Less0)= -trigamma(1-z(Less0),1,debug)+(pi./sin(pi*z(Less0))).^2; end;
    % integers:
    if(~isempty(Integers)) y(Integers) = trigamma(z(Integers),2,debug); end;
    % half-integers:
    if(~isempty(HalfInts)) y(HalfInts) = trigamma(z(HalfInts),3,debug); end;
    % negative integers:
    if(~isempty(NegInts))  y(NegInts)  = Inf * NegInts; end;
end

y = reshape(y,dim(1),dim(2));
return;




function psi = digamma(z,method,debug)
%
%  psi = digamma(z)   ... Digamma-Function for real argument z.
%
%  digamma(z) = d/dz log(gamma(z)) = gamma'(z)/gamma(z)
%
%  if 'z' is a matrix, then the digamma-function is evaluated for
%  each element. Results may be inaccurate for real arguments < 10
%  which are neither integers nor half-integers.
%
%  psi = digamma(z,method)
%
%  possible values for optional argument 'method':
%    method = 1           : quick asymptotic series expansion (approximate)
%    method = 2           : finite recursion for integer values (exact)
%    method = 3           : finite recursion for half-integer values (exact)
%    method = 4 (default) : automatic selection of 1,2 or 3 for individual
%                           elements in z whichever is appropriate.
%
%  see also:  trigamma, gamma, gammaln, gammainc, specfun

%  reference: Abramowitz & Stegun, "Handbook of Mathematical Functions"
%             Chapter "Gamma Function and Related Functions" :
%  implemented by: Christoph Mecklenbraeuker
%             (email: cfm@sth.ruhr-uni-bochum.de), July 1, 1995.


dim = size(z); 				% save original matrix dimension
z = reshape(z,dim(1)*dim(2),1); 	% make a column vector
I1 = ones(length(z),1); 		% auxiliary vector of ones

if(nargin==1)
    method=4; debug=0;
elseif(nargin==2)
    debug=0;
end;

if(debug == 1) 				% if debug==1: track recursion
    [m,n] = size(z);
    fprintf(1,'digamma: method = %d, size(z)=[%d %d],\t min(z)=%f, max(z)=%f\n',...
        method,m,n,min(min(z)),max(max(z)));
end;


if(method==1) 				% use 8th order asymptotic expansion
    if(any(z<1))
        fprintf(1,'Warning: some elements in argument of "digamma(z,1)" are < 1\n');
        fprintf(1,'minimal argument = %g: digamma-result is inaccurate!\n',min(min(z)));
    end
    % calculate powers of 1/z :
    w1 = 1./z; w2 = w1.*w1; w4 = w2.*w2; w6 = w2.*w4; w8 = w4.*w4;
    % generate coefficients of expansion: matrix with constant columns
    a = [ -I1/2   -I1/12   I1/120  -I1/252  I1/240 ];
    % make vector of powers of 1/z:
    w = [  w1      w2       w4       w6       w8   ];
    % calculate expansion by summing the ROWS of (a .* w) :
    psi = log(z) + sum((a.*w).').';
elseif(method==2)
    zmax = max(max(floor(z)));
    psitab = zeros(zmax,1);
    psitab(1) = -0.5772156649015328606;
    for n=1:zmax-1;
        psitab(n+1) = psitab(n) + 1/n; 	% generate lookup table
    end;
    psi = psitab(z);
elseif(method==3)
    zmax = max(max(floor(z)));
    psitab = zeros(zmax+1,1);
    psitab(1) = -0.5772156649015328606 - 2*log(2);  % = psi(1/2)
    for n=1:zmax;
        psitab(n+1) = psitab(n) + 2/(2*n-1); % generate lookup table
    end;
    psi = psitab(z+0.5);
elseif(method==4) 			% decide here which method to use
    Less0 = find(z<0); 			% negative arguments evaluated by reflexion formula
    Less1 = find(z>0 & z<1); 		% values between 0 and 1.
    fraction = rem(z,1); 			% fractional part of arguments
    f2 = rem(2*fraction,1);
    Integers = find(fraction==0 & z>0); 	% Index set of positive integer arguments
    NegInts  = find(fraction==0 & z<=0); 	% Index set of positive integer arguments
    HalfInts = find(abs(fraction-0.5)<1e-7 & z>0); % Index set of positive half-integers
    Reals    = find(f2>1e-7 & z>1); 	% Index set of all other arguments > 1
    if(~isempty(Reals)) psi(Reals)    = digamma(z(Reals),1,debug);    end;
    if(~isempty(Less1)) psi(Less1)    = digamma(z(Less1)+2,1,debug) - ...
            1./z(Less1)-1./(z(Less1)+1);end;
    % reflexion formula:
    if(~isempty(Less0)) psi(Less0) = digamma(1-z(Less0),1,debug) - pi./tan(pi*z(Less0)); end;
    if(~isempty(Integers)) psi(Integers) = digamma(z(Integers),2,debug); end;
    if(~isempty(HalfInts)) psi(HalfInts) = digamma(z(HalfInts),3,debug); end;
    if(~isempty(NegInts))  psi(NegInts)  = Inf * NegInts; end;
end

psi = reshape(psi,dim(1),dim(2));

return;


% Author: Patrick J. Wolfe
%         Signal Processing Group
%         Cambridge University Engineering Department
%         p.wolfe@ieee.org
% Johnston perceptual model initialisation
function M= mask( Sx, dft_length, Fs, nbits)

frame_overlap= dft_length/ 2;    
freq_val = (0:Fs/dft_length:Fs/2)';
half_lsb = (1/(2^nbits-1))^2/dft_length;

freq= freq_val;
thresh= half_lsb;
crit_band_ends = [0;100;200;300;400;510;630;770;920;1080;1270;...
        1480;1720;2000;2320;2700;3150;3700;4400;5300;6400;7700;...
        9500;12000;15500;Inf];

% Maximum Bark frequency
%
imax = max(find(crit_band_ends < freq(end)));

% Normalised (to 0 dB) threshold of hearing values (Fletcher, 1929) 
% as used  by Johnston.  First and last thresholds are corresponding 
% critical band endpoint values, elsewhere means of interpolated 
% critical band endpoint threshold values are used.
%
abs_thr = 10.^([38;31;22;18.5;15.5;13;11;9.5;8.75;7.25;4.75;2.75;...
        1.5;0.5;0;0;0;0;2;7;12;15.5;18;24;29]./10);
ABSOLUTE_THRESH = thresh.*abs_thr(1:imax);

% Calculation of tone-masking-noise offset ratio in dB
%
OFFSET_RATIO_DB = 9+ (1:imax)';

% Initialisation of matrices for bark/linear frequency conversion
% (loop increments i to the proper critical band)
%
num_bins = length(freq);
LIN_TO_BARK = zeros(imax,num_bins);
i = 1;
for j = 1:num_bins
    while ~((freq(j) >= crit_band_ends(i)) & ...
            (freq(j) < crit_band_ends(i+1))),
        i = i+1;
    end
    LIN_TO_BARK(i,j) = 1;
end

% Calculation of spreading function (Schroeder et al., 82)

spreading_fcn = zeros(imax);
summ = 0.474:imax;
spread = 10.^((15.81+7.5.*summ-17.5.*sqrt(1+summ.^2))./10);
for i = 1:imax
    for j = 1:imax
        spreading_fcn(i,j) = spread(abs(j-i)+1);
    end
end

% Calculation of excitation pattern function

EX_PAT = spreading_fcn* LIN_TO_BARK;

% Calculation of DC gain due to spreading function

DC_GAIN = spreading_fcn* ones(imax,1);


%Sx = X.* conj(X);

C = EX_PAT* Sx;

% Calculation of spectral flatness measure SFM_dB
%
[num_bins num_frames] = size(Sx);
k = 1/num_bins;
SFM_dB = 10.*log10((prod(Sx).^k)./(k.*sum(Sx)+eps)+ eps);

% Calculation of tonality coefficient and masked threshold offset
%
alpha = min(1,SFM_dB./-60);
O_dB = OFFSET_RATIO_DB(:,ones(1,num_frames)).*...
    alpha(ones(length(OFFSET_RATIO_DB),1),:) + 5.5;

% Threshold calculation and renormalisation, accounting for absolute 
% thresholds

T = C./10.^(O_dB./10);
T = T./DC_GAIN(:,ones(1,num_frames));
T = max( T, ABSOLUTE_THRESH(:, ones(1, num_frames)));

% Reconversion to linear frequency scale 

%M = 1.* sqrt((LIN_TO_BARK')*T);
M= LIN_TO_BARK'* T;