boblsturm@0: function [F,D] = ifgram(X, N, W, H, SR)
boblsturm@0: % [F,D] = ifgram(X, N, W, H, SR)       Instantaneous frequency by phase deriv.
boblsturm@0: %    X is a 1-D signal.  Process with N-point FFTs applying a W-point 
boblsturm@0: %    window, stepping by H points; return (N/2)+1 channels with the 
boblsturm@0: %    instantaneous frequency (as a proportion of the sampling rate) 
boblsturm@0: %    obtained as the time-derivative of the phase of the complex spectrum
boblsturm@0: %    as described by Toshihiro Abe et al in ICASSP'95, Eurospeech'97
boblsturm@0: %    Same arguments and some common code as dpwebox/stft.m.
boblsturm@0: %    Calculates regular STFT as side effect - returned in D.
boblsturm@0: % after 1998may02 dpwe@icsi.berkeley.edu
boblsturm@0: % 2001-03-05 dpwe@ee.columbia.edu  revised version
boblsturm@0: % 2001-12-13 dpwe@ee.columbia.edu  Fixed to work when N != W
boblsturm@0: % $Header: $
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: if nargin < 2;  N = 256; end
boblsturm@0: if nargin < 3;  W = N;   end
boblsturm@0: if nargin < 4;  H = W/2; end
boblsturm@0: if nargin < 5;  SR = 1;  end
boblsturm@0: 
boblsturm@0: s = length(X);
boblsturm@0: % Make sure it's a single row
boblsturm@0: if size(X,1) > 1
boblsturm@0:   X = X';
boblsturm@0: end
boblsturm@0: 
boblsturm@0: %win = [0,hanning(W-1)'];
boblsturm@0: win = 0.5*(1-cos([0:(W-1)]/W*2*pi));
boblsturm@0: 
boblsturm@0: % Window for discrete differentiation
boblsturm@0: T = W/SR;
boblsturm@0: dwin = -pi / T * sin([0:(W-1)]/W*2*pi);
boblsturm@0: 
boblsturm@0: % sum(win) takes out integration due to window, 2 compensates for neg frq
boblsturm@0: norm = 2/sum(win);
boblsturm@0: 
boblsturm@0: % How many complete windows?
boblsturm@0: nhops = 1 + floor((s - W)/H);
boblsturm@0: 
boblsturm@0: F = zeros(1 + N/2, nhops);
boblsturm@0: D = zeros(1 + N/2, nhops);
boblsturm@0: 
boblsturm@0: nmw1 = floor( (N-W)/2 );
boblsturm@0: nmw2 = N-W - nmw1;
boblsturm@0: 
boblsturm@0: ww = 2*pi*[0:(N-1)]*SR/N;
boblsturm@0: 
boblsturm@0: for h = 1:nhops
boblsturm@0:   u = X((h-1)*H + [1:W]);
boblsturm@0: %  if(h==0)
boblsturm@0: %	plot(u)
boblsturm@0: %  end
boblsturm@0:   % Apply windows now, while the length is right
boblsturm@0:   wu = win.*u;
boblsturm@0:   du = dwin.*u;
boblsturm@0:   
boblsturm@0:   % Pad or truncate samples if N != W
boblsturm@0:   if N > W
boblsturm@0:     wu = [zeros(1,nmw1),wu,zeros(1,nmw2)];
boblsturm@0:     du = [zeros(1,nmw1),du,zeros(1,nmw2)];
boblsturm@0:   end
boblsturm@0:   if N < W
boblsturm@0:     wu = wu(-nmw1+[1:N]);
boblsturm@0:     du = du(-nmw1+[1:N]);
boblsturm@0:   end
boblsturm@0:   % FFTs of straight samples plus differential-weighted ones
boblsturm@0:   t1 = fft(fftshift(du));
boblsturm@0:   t2 = fft(fftshift(wu));
boblsturm@0:   % Scale down to factor out length & window effects
boblsturm@0:   D(:,h) = t2(1:(1 + N/2))'*norm;
boblsturm@0: 
boblsturm@0:   % Calculate instantaneous frequency from phase of differential spectrum
boblsturm@0:   t = t1 + j*(ww.*t2);
boblsturm@0:   a = real(t2);
boblsturm@0:   b = imag(t2);
boblsturm@0:   da = real(t);
boblsturm@0:   db = imag(t);
boblsturm@0:   instf = (1/(2*pi))*(a.*db - b.*da)./((a.*a + b.*b)+(abs(t2)==0));
boblsturm@0:   % 1/2pi converts rad/s into cycles/s
boblsturm@0:   % sampling rate already factored in as constant in dwin & ww
boblsturm@0:   % so result is in Hz
boblsturm@0:   
boblsturm@0:   F(:,h) = instf(1:(1 + N/2))';
boblsturm@0:     
boblsturm@0: end;
boblsturm@0: