plp_audw.m

这是一个用于语音信号处理的工具箱

% Computes auditory weighting functions           
% function [nfilt, cb, ibegen]=  plp_audw(npoint, sf)
%
% npoint (in): number of points in the fft spectrum                 
% nfilt  (in): number of samples of auditory spectrum               
%              equally spaced on the bark scale;                     
%              1st filter at dc and last at the Nyquist frequency   
%                                                                    
% cb    (out): array of weighting coefficients to simulate               
%              critical band spectral resolution                        
%              and equal loudness sensitivity of hearing                
%              on npoint speech power spectrum                          
% ibegen(out): three-dimensional array which indicates where         
%              to begin and where to end integration                
%              of speech spectrum and where the given            
%              weighting function starts in array cb                
%              get Nyquist frequency in Bark                
%
function [nfilt, cb, ibegen]=  plp_audw(npoint, sf)

fnqbar = 6.0*log((sf/1200)+sqrt((sf/1200)^2+1.0));

% compute number of filters for less than 1 Bark spacing 
nfilt = floor(fnqbar)+2;

% frequency -> fft spectral sample conversion            
f2samp = (npoint-1)/(sf/2.0);

% compute filter step in Bark 
zdel = fnqbar/(nfilt-1);    
    
icount = 1;                       
% loop over all weighting functions  
for j = 2:(nfilt-1)
	ibegen(j,3) = icount;

	% get center frequency of the j-th filter in Bark        
	z0 = zdel*(j-1);
                  
	% get center frequency of the j-th filter in Hz          
	f0 = 600.0* (exp(z0/6.0) - exp(-z0/6.0))/2.0;
	
	% get low-cut frequency of j-th filter in Hz             
	fl = 600.0* (exp((z0 - 2.5)/6.0) - exp(-(z0 - 2.5)/6.0))/2.0;
                  
	ibegen(j,1) = round(fl*f2samp)+1;
	if (ibegen(j,1) < 1)
	    	ibegen(j,1) = 1;
	end
	
	% get high-cut frequency of j-th filter in Hz           
	fh = 600.0*(exp((z0 + 1.3)/6.0) - exp(-( z0 + 1.3)/6.0)) /2.0;
                   
	ibegen(j,2) = round(fh*f2samp) + 1;  
	if (ibegen(j,2) > npoint)
	    	ibegen(j,2) = npoint;
	end
	
	for i = ibegen(j,1):ibegen(j,2),
		% get frequency of j-th spectral point in Hz   
	    	freq = (i-1)/f2samp;

		% get frequency of j-th spectral point in Bark 
	    	x = freq / 600.0;
                            
	    	z = 6.0*log(x+sqrt(x^2+1.0));
	    
		% normalize by center frequency in barks:      
	    	z = z-z0;	    

		% compute weighting 
	    	if z <= -0.5                            
			cb(icount) = 10^(z+0.5);
	    	elseif z >= 0.5
			cb(icount) = 10^((-2.5)*(z-0.5));
	    	else
			cb(icount) = 1.0;
		end

		%     calculate the 40 db equal-loudness curve        
		%     at center frequency and add it to the weighting 

		fsq = f0^2;
		rsss = fsq^2*(fsq+1200^2)/(((fsq+400^2)^2)*(fsq+3100^2));
	   	 
		cb(icount) = rsss*cb(icount);
                       
		% add the qual-loundness curve to the weighting    
		icount = icount+1;
	end 
end