boblsturm@0: function [p,m,S] = ifptrack(d,w,sr,fminl,fminu,fmaxl,fmaxu)
boblsturm@0: % [p,m,S] = ifptrack(d,w,sr,fminl,fminu,fmaxl,fmaxu)
boblsturm@0: %     Pitch track based on inst freq.
boblsturm@0: %     Look for adjacent bins with same inst freq.
boblsturm@0: %     d is the input waveform.  sr is its sample rate
boblsturm@0: %     w is the basic STFT DFT length (window is half, hop is 1/4)
boblsturm@0: %     S returns the underlying complex STFT.
boblsturm@0: %     fmin,fmax define ramps at edge of sensitivity
boblsturm@0: % 2006-05-03 dpwe@ee.columbia.edu
boblsturm@0: 
boblsturm@0: %   Copyright (c) 2006 Columbia University.
boblsturm@0: % 
boblsturm@0: %   This file is part of LabROSA-coversongID
boblsturm@0: % 
boblsturm@0: %   LabROSA-coversongID is free software; you can redistribute it and/or modify
boblsturm@0: %   it under the terms of the GNU General Public License version 2 as
boblsturm@0: %   published by the Free Software Foundation.
boblsturm@0: % 
boblsturm@0: %   LabROSA-coversongID is distributed in the hope that it will be useful, but
boblsturm@0: %   WITHOUT ANY WARRANTY; without even the implied warranty of
boblsturm@0: %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
boblsturm@0: %   General Public License for more details.
boblsturm@0: % 
boblsturm@0: %   You should have received a copy of the GNU General Public License
boblsturm@0: %   along with LabROSA-coversongID; if not, write to the Free Software
boblsturm@0: %   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
boblsturm@0: %   02110-1301 USA
boblsturm@0: % 
boblsturm@0: %   See the file "COPYING" for the text of the license.
boblsturm@0: 
boblsturm@0: % downweight fundamentals below here
boblsturm@0: if nargin < 4; fminl = 150; end
boblsturm@0: if nargin < 5; fminu = 300; end
boblsturm@0: % highest frequency we look to
boblsturm@0: if nargin < 6; fmaxl = 2000; end
boblsturm@0: if nargin < 7; fmaxu = 4000; end
boblsturm@0: 
boblsturm@0: 
boblsturm@0: % Calculate the inst freq gram
boblsturm@0: [I,S] = ifgram(d,w,w/2,w/4,sr);
boblsturm@0: 
boblsturm@0: % Only look at bins up to 2 kHz
boblsturm@0: maxbin = round(fmaxu * (w/sr) );
boblsturm@0: %maxbin = size(I,1)
boblsturm@0: minbin = round(fminl * (w/sr) );
boblsturm@0: 
boblsturm@0: % Find plateaus in ifgram - stretches where delta IF is < thr
boblsturm@0: ddif = [I(2:maxbin, :);I(maxbin,:)] - [I(1,:);I(1:(maxbin-1),:)];
boblsturm@0: 
boblsturm@0: % expected increment per bin = sr/w, threshold at 3/4 that
boblsturm@0: dgood = abs(ddif) < .75*sr/w;
boblsturm@0: 
boblsturm@0: % delete any single bins (both above and below are zero);
boblsturm@0: dgood = dgood .* ([dgood(2:maxbin,:);dgood(maxbin,:)] >  0 | [dgood(1,:);dgood(1:(maxbin-1),:)] > 0);
boblsturm@0: 
boblsturm@0: % check it out
boblsturm@0: %p = dgood;
boblsturm@0: 
boblsturm@0: % reconstruct just pitchy cells?
boblsturm@0: %r = istft(p.*S,w,w/2,w/4);
boblsturm@0: 
boblsturm@0: p = 0*dgood;
boblsturm@0: m = 0*dgood;
boblsturm@0: 
boblsturm@0: % For each frame, extract all harmonic freqs & magnitudes
boblsturm@0: for t = 1:size(I,2)
boblsturm@0:   ds = dgood(:,t)';
boblsturm@0:   lds = length(ds);
boblsturm@0:   % find nonzero regions in this vector
boblsturm@0:   st = find(([0,ds(1:(lds-1))]==0) & (ds > 0));
boblsturm@0:   en = find((ds > 0) & ([ds(2:lds),0] == 0));
boblsturm@0:   npks = length(st);
boblsturm@0:   frqs = zeros(1,npks);
boblsturm@0:   mags = zeros(1,npks);
boblsturm@0:   for i = 1:length(st)
boblsturm@0:     bump = abs(S(st(i):en(i),t));
boblsturm@0:     frqs(i) = (bump'*I(st(i):en(i),t))/(sum(bump)+(sum(bump)==0));
boblsturm@0:     mags(i) = sum(bump);
boblsturm@0:     if frqs(i) > fmaxu
boblsturm@0:       mags(i) = 0;
boblsturm@0:       frqs(i) = 0;
boblsturm@0:     elseif frqs(i) > fmaxl
boblsturm@0:       mags(i) = mags(i) * max(0, (fmaxu - frqs(i))/(fmaxu-fmaxl));
boblsturm@0:     end
boblsturm@0:     % downweight magnitudes below? 200 Hz
boblsturm@0:     if frqs(i) < fminl
boblsturm@0:       mags(i) = 0;
boblsturm@0:       frqs(i) = 0;
boblsturm@0:     elseif frqs(i) < fminu
boblsturm@0:       % 1 octave fade-out
boblsturm@0:       mags(i) = mags(i) * (frqs(i) - fminl)/(fminu-fminl);
boblsturm@0:     end
boblsturm@0:     if frqs(i) < 0 
boblsturm@0:       mags(i) = 0;
boblsturm@0:       frqs(i) = 0;
boblsturm@0:     end
boblsturm@0:     
boblsturm@0:   end
boblsturm@0: 
boblsturm@0: % then just keep the largest at each frame (for now)
boblsturm@0: %  [v,ix] = max(mags);
boblsturm@0: %  p(t) = frqs(ix);
boblsturm@0: %  m(t) = mags(ix);
boblsturm@0:   % No, keep them all
boblsturm@0:   %bin = st;
boblsturm@0:   bin = round((st+en)/2);
boblsturm@0:   p(bin,t) = frqs;
boblsturm@0:   m(bin,t) = mags;
boblsturm@0: end
boblsturm@0: 
boblsturm@0: %% Pull out the max in each column
boblsturm@0: %[mm,ix] = max(m);
boblsturm@0: %% idiom to retrieve different element from each column
boblsturm@0: %[nr,nc] = size(p);
boblsturm@0: %pp = p((nr*[0:(nc-1)])+ix);
boblsturm@0: %mm = m((nr*[0:(nc-1)])+ix);
boblsturm@0: % r = synthtrax(pp,mm,sr,w/4);
boblsturm@0: 
boblsturm@0: %p = pp;
boblsturm@0: %m = mm;
boblsturm@0: