Matlab for Finding Interpolated Spectral Peaks

Matlab listing: findpeaks.m

function [peakamps,peaklocs,peakwidths,resid] = ...
  findpeaks(data,npeaks,minwidth,maxwidth,minpeak,debug);
%FINDPEAKS Find up to npeaks interpolated peaks in data.
%
%  [peakamps,peaklocs,peakwidths,resid] =
%     findpeaks(data,npeaks,minwidth,maxwidth,minpeak,debug);
%
%       finds up to npeaks interpolated peaks in real data.
%       A peak is rejected if its width is
%         less than minwidth samples wide(default=1), or
%         less than minpeak in magnitude (default=min(data)).
%       Quadratic interpolation is used for peak interpolation.
%       Left-over data with peaks removed is returned in resid.
%       Peaks are returned in order of decreasing amplitude.
%         They can be sorted in order of location as follows:
%           [peaklocs,sorter] = sort(peaklocs);
%           amps = zeros(size(peakamps));
%           npeaks = length(peaklocs);
%           for col=1:npeaks,
%             amps(:,col) = peakamps(:,sorter(col));
%           end;
%         They can be sorted in order of width as follows:
%           [peakwidths,sorter] = sort(peakwidths);
%           peakamps = peakamps(sorter);
%           peaklocs = peaklocs(sorter);

len = length(data);

if nargin<6,
  debug = 0;
end;

if nargin<5,
 minpeak = min(data);
end;

if nargin<4,
 maxwidth = 0;
end;

if nargin<3,
 minwidth = 1;
end;

if nargin<2,
 npeaks = len;
end;

peakamps = zeros(1,npeaks);
peaklocs = zeros(1,npeaks);
peakwidths = zeros(1,npeaks);
if debug, peaksarr = 1.1*max(data)*ones(size(data)); end;
if debug, orgdata = data; end;
if debug, npeaks, end

nrej = 0;
pdebug=debug;
ipeak=0;
while ipeak<npeaks
 [ploc, pamp, pcurv] = maxr(data);
 if (pamp==minpeak), warning('findpeaks:min amp reached');
        break;
 end
 plocq = round(ploc);
 ulim = min(len,plocq+1);
 camp = pamp;
 %
 % Follow peak down to determine its width
 %
 drange = max(data) - minpeak; % data dynamic range
 tol = drange * 0.01;
 dmin = camp;
 while ((ulim<len) & (data(ulim)<=dmin+tol)),
   camp = data(ulim);
   ulim = ulim + 1;
   if (camp<dmin), dmin=camp; end
 end;
 ulim = ulim - 1;
 lamp = camp;

 llim = max(1,plocq-1);
 camp = pamp;
 dmin = camp;
 while ((llim>1) & (data(llim)<=dmin+tol)),
   camp = data(llim);
   llim = llim - 1;
   if (camp<dmin), dmin=camp; end
 end;
 llim = llim + 1;
 uamp = camp;
 %
 % Remove the peak
 %
 data(llim:ulim) = min(lamp,uamp) * ones(1,ulim-llim+1);
 %
 % Reject peaks which are too narrow (indicated by zero loc and amp)
 %
 pwid = ulim - llim + 1;
 if ~(pwid < minwidth),
   ipeak = ipeak + 1;
   peaklocs(ipeak) = ploc;
   peakamps(ipeak) = pamp;
   peakwidths(ipeak) = - 1/pcurv;
   nrej = 0;
   if pdebug>1
      peaksarr(plocq) = pamp;
      maxloc = min(len,2*round(max(peaklocs)));
      ttl = sprintf(...
        'Peak %d = %0.2f at %0.2f, width %d',ipeak,pamp,ploc,pwid);
      if x == -1
        pdebug = 0;
      end
   end
 else
   nrej = nrej + 1;
   if (nrej >= 10),
     warning('*** findpeaks: giving up (10 rejected peaks in a row)');
     break;
   end;
 end;
end;
if (ipeak<npeaks),
 warning(sprintf(...
   '*** peaks.m: only %d peaks found instead of %d',ipeak,npeaks));
 peakamps = peakamps(1:ipeak);
 peaklocs = peaklocs(1:ipeak);
 peakwidths = peakwidths(1:ipeak);
end;
resid = data;

Matlab listing: maxr.m

function [xi,yi,hc] = maxr(a)
%MAXR   Find interpolated maximizer(s) and max value(s)
%       for (each column of) a.
%
%               [xi,yi,hc] = maxr(a)
%
%  Calls max() followed by qint() for quadratic interpolation.
%
   [m,n] = size(a);
   if m==1, a=a'; t=m; m=n; n=t; end;
   [y,x] = max(a);
   xi=x;    % vector of maximizer locations, one per col of a
   yi=y;    % vector of maximum values, one per column of a
   if nargout>2, hc = zeros(1,n); end
   for j=1:n,   % loop over columns
     if x(j)>1, % only support interior maxima
       if x(j)<m,
         [xdelta,yij,cj] = qint(a(x(j)-1,j),y(j),a(x(j)+1,j));
         xi(j) = x(j) + xdelta;
         if nargout>2, hc(j) = cj; end
         if (nargout>1), yi(j) = yij; end
       end;
     end;
   end;

Matlab listing: qint.m

function [p,y,a] = qint(ym1,y0,yp1)
%QINT   Quadratic interpolation of 3 uniformly spaced samples
%
%               [p,y,a] = qint(ym1,y0,yp1)
%
%       returns extremum-location p, height y, and half-curvature a
%       of a parabolic fit through three points.
%       The parabola is given by y(x) = a*(x-p)^2+b,
%       where y(-1)=ym1, y(0)=y0, y(1)=yp1.

   p = (yp1 - ym1)/(2*(2*y0 - yp1 - ym1));
   if nargout>1
     y = y0 - 0.25*(ym1-yp1)*p;
   end;
   if nargout>2
     a = 0.5*(ym1 - 2*y0 + yp1);
   end;

Matlab listing: zpfmin.m

function [zpmin] = zpfmin(window,MT,freqBiasMax);
% [zpmin] = zpfmin(window,MT,freqBiasMax)
% returns the minimum zero-padding factor to achieve prescribed
% tolerances in sinusoidal peak parameters measured using the QIFFT
% method with the  window.
%
% INPUTS:
%   window      = 'rect', 'hamming', 'hann', or 'blackman'
%   MT          = window length (seconds)
%   freqBiasMax = max frequency bias allowed in peak frequency (Hz)
%   ampBiasMax  = max peak-magnitude bias allowed (relative fraction)
%
%   If ampBiasMaxPct is not specified, amplitude bias is not considered.
%
% RETURNED:
%   zpmin = minimum zero-padding factor (1 means no zero padding)
%
% EXAMPLE
%   Assuming the window contains one cycle of a sinusoid whose
%   frequency f must be estimated to within 0.1% accuracy, we have
%   zpfmin('blackman',1/f,0.001*f) = 1.9 for min zero-padding factor.
%   Similarly,
%     2 cycles under the Blackman window yields zpf > 1.5
%     4 cycles under the Blackman window yields zpf > 1.2
%     6 cycles under the Blackman window yields zpf > 1.02
%
% Reference:
%    "Design Criteria for Simple Sinusoidal Parameter Estimation
%     based on Quadratic Interpolation of FFT Magnitude Peaks", by
%    Mototsugu Abe and Julius O. Smith,
%    Audio Engineering Society Convention preprint 6256", NY, 2004.

%!test
%! assert(zpfmin('blackman',1,0.001),2.3609,0.001);

% disp(sprintf('%s window',window));
switch window
 case 'rect'
  c0 = 0.4467;
  r0 = -0.3218;
  c1 = 0.8560;
  r1 = -0.2366;
 case 'hamming'
  c0 = 0.2456;
  r0 = -0.3282;
  c1 = 0.4381;
  r1 = -0.2451;
 case 'hann'
  c0 = 0.2436;
  r0 = -0.3288;
  c1 = 0.4149;
  r1 = -0.2456;
 case 'blackman'
  c0 = 0.1868;
  r0 = -0.3307;
  c1 = 0.3156;
  r1 = -0.2475;
 otherwise
  error(sprintf('Window type %s unknown - say help zpfmin',window));
end

zpmin = c0*(MT*freqBiasMax)^r0;

if nargin>3
  zpminamp = c1*(MT*freqBiasMax)^r1;
  zpmin = max(zpmin,zpminamp);
end

if strcmp(window,'rect')
  zpmin=max(1.7,zpmin); % See Abe & Smith 2004, Fig. 2 re. why this
else
  zpmin=max(1,zpmin);
end

Matlab listing: oboeanal.m

The following Matlab script illustrates use of the findpeaks function above to determine the pitch of an oboe tone (given the general location of the correct spectral peak) and configure a spectrum analysis using the rectangular, Hamming, and Blackman windows. This script was used to create Figures 3.16 through 3.18.

diary './dia/oboeanal.dia';
diary on;

zpf = 5; % desired min zero-padding factor in all FFTs
fn = '../wav/oboe.ff.C4B4.wav';
[x,fs,nbits] = wavread(fn);
nx = length(x);
dur = nx/fs;
disp(sprintf(['Read %s, fs = %f, nbits = %d, ',...
              'length = %d samples = %0.1f sec'],...
             fn,fs,nbits,nx,dur));
soundsc(x,fs);
M = fs/10; % 100 ms window for pitch estimation
N = 2^nextpow2(M*zpf) % zero pad
Xzp = fft(x,N);

figure(K-1);
fmax = 500; % maximum freq in Hz
kmax = fmax*N/fs; % maximum freq in bins
prange = 1:kmax+1;
Xdb = dbn(Xzp(prange)); % dB normalized to 0 max
f = fs*[0:N-1]/N/1000; % frequency axis in kHz
plot(f(prange),Xdb);
ylabel('Magnitude (dB)');
xlabel('Normalized Frequency (cycles/sample)');
[amps,bins] = findpeaks(Xdb,1);
k0 = bins(1)-1;
f0 = k0*fs/N;
nP = fs/f0; % samples per period (used below)
disp(sprintf('Estimated pitch = %f Hz',f0));

% Test pitch estimate:
pause(ceil(dur)); % let previous sound finish
T = 1/fs;
t = 0:T:dur;
y = sawtooth(2*pi*f0*t);
sound(y,fs)

winnames={'none','boxcar','none','hamming','none','blackman'};

% Configure window analysis

for K=[2 4 6]
  wintype = winnames{K};
  M = ceil(K*nP); % min Blackman length
  cmd = sprintf('w = %s(M);',wintype);
  eval(cmd);
  Nw = 2^nextpow2(M*zpf) % zero pad
  xw = x(1:M) .* w;
  Xwzp = fft(xw,Nw);
  %fmax = 5*f0; % maximum freq in Hz
  fmax = fs/4;
  kmax = fmax*Nw/fs; % maximum freq in bins
  prange = 1:kmax+1;
  Xwdb = dbn(Xwzp(prange)); % dB normalized to 0 max

  figure(K);
  subplot(2,1,1);
  plot(xw,'-k');
  axis tight;
  title(sprintf('Oboe data - %s window, K=%d',wintype,K));
  ylabel('Amplitude');
  xlabel('Time (samples)');
  subplot(2,1,2);
  f = fs*[0:Nw-1]/Nw/1000; % frequency axis in kHz
  plot(f(prange),Xwdb,'-k'); grid on;
  axis([f(prange(1)) f(prange(end)) -90 0]);
  ylabel('Magnitude (dB)');
  xlabel('Frequency (kHz)');

  plotfile = sprintf('../eps/oboe%s.eps',wintype);
  saveplot(plotfile);

end

diary off;

Sign in