yading@11: function [iploc, ipmag, ipphase] = peakinterp(mX, pX, ploc)
yading@11: % Parabolic interpolation of spectral peaks
yading@11: % mX: magnitude spectrum, pX: phase spectrum, ploc: locations of peaks
yading@11: % iploc, ipmag, ipphase: interpolated values
yading@11: % note that ploc values are assumed to be between 2 and length(mX)-1
yading@11: val = mX(ploc); % magnitude of peak bin
yading@11: lval = mX(ploc-1); % magnitude of bin at left
yading@11: rval= mX(ploc+1); % magnitude of bin at right
yading@11: iploc = ploc + .5*(lval-rval)./(lval-2*val+rval); % center of parabola
yading@11: ipmag = val-.25*(lval-rval).*(iploc-ploc); % magnitude of peaks
yading@11: ipphase = interp1(1:length(pX),pX,iploc,'linear'); % phase of peaks