www.pudn.com > LPCToolbox.rar > GETPEAKS.M

%--------------------------------- 
%  [estpk,rawpk] = getpeaks(a, fftlen, Fs) 
%    a is the list of linear prediction coefficients. GETPEAKS 
%    estimates the locations of the peaks and their bandwidths in the 
%    spectrum of 'a', using a peak-picking + parabolic interpolation 
%    method. See Markel & Gray, 1976, Sec. 7.2.2 for details. 
%     
%    a : lpc coeffs 
%    fftlen: zero-padded length for fft 
%    Fs: sampling frequency 
% 
%   rawpk, estpk   : Nx2 matrices, N=# of peaks 
%           pk(:,1) = locations of peaks (w.r.t. Fs) 
%           pk(:,2) = corr. bandwidths  
% 
%--------------------------------- 
function [estpk,rawpk] = getpeaks(a, fftlen, Fs) 
 
F = fft(a,fftlen);          % zero-padded to fftlen 
 
fx = powerspectrum(F);      % all-zero power spectrum, with only 
                            % positive frequencies. 
 
ts = Fs/fftlen;             % frequency resolution 
 
% find the extrema, and do the 3pt parabolic interpolation 
% to find the "true" location of each valley. 
 
l=minima(fx); 
   
rawpk = []; estpk = []; 
for i=1:length(l), 
  j = l(i);        % j = location of the i'th peak 
  [fn,fi,bw] = threepointinterp(fx, j, ts);  
  rawpk(i,:) = floor([fn  0]); 
  estpk(i,:) = floor([fi bw]); 
end 
   
rawpk = [rawpk; 0 0; 0 0; 0 0]; % in case there are too few formants 
estpk = [estpk; 0 0; 0 0; 0 0]; 
   
%================================================================ 
 
%--------------------------------- 
% [fn, fi] = threepointinterp(f, j, ts) 
% 
%   Estimate the "true peak" using three-point parabolic  
%   interpolation. See Markel & Gray, 1976, Sec 7.2.2.,  
%   eqns. 7.4 & 7.5. 
% 
%   f = spectrum, j = location of the i'th peak 
%   ts = frequency quantum 
% 
%   fn : quantal (or raw) peak 
%   fi : estimated peak 
%   bw : est. bandwidth of estimated peak 
% 
%--------------------------------- 
 
function [fn, fi, bw] = threepointinterp(f, j, ts) 
 
 % We (along with ILS) use (j-1)*ts since the array indices 
 % start with 1. So j=1 => 0 Hz, j=2 => 1*ts Hz, etc. 
 
 fn = (j-1) * ts; 
 
 fj = 1 / f(j);       % these inversions are for the 
 fjm1 = 1 / f(j-1);   % all-zero -> all-pole adjustment 
 fjp1 = 1 / f(j+1);   % (so that we get the right values 
                      % for the amplitudes [pk]). 
 
 %------------ 
 % The Markel & Gray algorithm is implemented in a 
 % somewhat convoluted way in ILS. This has been 
 % duplicated for authenticity. For comparision, 
 % here is a more understandable implementation: 
 % 
 % aa = (fjm1+fjp1)/2 - fj; 
 % bb = (fjp1-fjm1)/2; 
 % cc = fj; 
 % delta = (-bb/(2*aa)) * ts; 
 % fi = delta + fn; % fi=location of the peak 
 %------------ 
 
 % The following code is the convoluted implementation,  
 % which is pretty much copied directly from ILS. 
 
 tss = ts*ts; 
 
 %------ location of peak  
 
 aa = (fjm1 - 2*fj) + fjp1; 
 aa = aa * (0.5/tss); 
 
 bb = (fjp1-fjm1) * (0.5/ts); 
 fi = - bb*(0.5/aa) + fn;       % fi=location of the peak 
 
 %------ bandwidth 
 
 cc = fj; 
 ee = cc - (bb*bb)*(.25/aa); 
 dd = cc - 0.5*ee; 
 gg = (bb*bb) - 4*(aa*dd);  
 
 bw = - sqrt(gg) / aa;      % bw=bandwidth 
 pk = ee;                   % estimated amplitude of the peak  
 
%=============================================================== 
% ps = powerspectrum(F) 
%  F is the a complex spectrum (the output of FFT). Throw out  
%  the negative frequencies, and convert to a power spectrum (ps). 
 
function  ps = powerspectrum(F) 
 
n = length(F);    
i = [1:ceil(n/2)]; 
 
ps = abs(F(i)) .^ 2; 
 
%    f  : digital frequencies for ps. Multiply by Fs to 
%         get the analog frequencies. 
% 
% f = (i-1) * (1/n); 
 
%=============================================================== 
% l = minima(x) 
% 
% find when der(x) goes from neg to pos. 
% output l (list of indices) is sorted. 
 
function l = minima(x) 
 
% x(3:end)    a b|c d e f g ... x y z| 
% x(2:end-1)    a|b c d e f g ... x y|z 
% x(1:end-2)     |a b c d e f g ... x|y z 
 
d1 = x(2:end-1) - x(1:end-2); 
d2 = x(3:end)   - x(2:end-1); 
i = find((d1 < 0) & (d2 > 0)); 
 
% d1 and d2 are indexed w.r.t. x(2:end-1). Update 
% the locations of the minima so that they are indexed 
% w.r.t. x(1:end). 
 
l = i+1;