function [F,D] = ifgram(X, N, W, H, SR)
% [F,D] = ifgram(X, N, W, H, SR) Instantaneous frequency by phase deriv.
% X is a 1-D signal. Process with N-point FFTs applying a W-point
% window, stepping by H points; return (N/2)+1 channels with the
% instantaneous frequency (as a proportion of the sampling rate)
% obtained as the time-derivative of the phase of the complex spectrum
% as described by Toshihiro Abe et al in ICASSP'95, Eurospeech'97
% Same arguments and some common code as dpwebox/stft.m.
% Calculates regular STFT as side effect - returned in D.
% after 1998may02 dpwe@icsi.berkeley.edu
% 2001-03-05 dpwe@ee.columbia.edu revised version
% 2001-12-13 dpwe@ee.columbia.edu Fixed to work when N != W
% $Header: $
% Copyright (c) 2006 Columbia University.
%
% This file is part of LabROSA-coversongID
%
% LabROSA-coversongID is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License version 2 as
% published by the Free Software Foundation.
%
% LabROSA-coversongID is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with LabROSA-coversongID; if not, write to the Free Software
% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
% 02110-1301 USA
%
% See the file "COPYING" for the text of the license.
if nargin < 2; N = 256; end
if nargin < 3; W = N; end
if nargin < 4; H = W/2; end
if nargin < 5; SR = 1; end
s = length(X);
% Make sure it's a single row
if size(X,1) > 1
X = X';
end
%win = [0,hanning(W-1)'];
win = 0.5*(1-cos([0:(W-1)]/W*2*pi));
% Window for discrete differentiation
T = W/SR;
dwin = -pi / T * sin([0:(W-1)]/W*2*pi);
% sum(win) takes out integration due to window, 2 compensates for neg frq
norm = 2/sum(win);
% How many complete windows?
nhops = 1 + floor((s - W)/H);
F = zeros(1 + N/2, nhops);
D = zeros(1 + N/2, nhops);
nmw1 = floor( (N-W)/2 );
nmw2 = N-W - nmw1;
ww = 2*pi*[0:(N-1)]*SR/N;
for h = 1:nhops
u = X((h-1)*H + [1:W]);
% if(h==0)
% plot(u)
% end
% Apply windows now, while the length is right
wu = win.*u;
du = dwin.*u;
% Pad or truncate samples if N != W
if N > W
wu = [zeros(1,nmw1),wu,zeros(1,nmw2)];
du = [zeros(1,nmw1),du,zeros(1,nmw2)];
end
if N < W
wu = wu(-nmw1+[1:N]);
du = du(-nmw1+[1:N]);
end
% FFTs of straight samples plus differential-weighted ones
t1 = fft(fftshift(du));
t2 = fft(fftshift(wu));
% Scale down to factor out length & window effects
D(:,h) = t2(1:(1 + N/2))'*norm;
% Calculate instantaneous frequency from phase of differential spectrum
t = t1 + j*(ww.*t2);
a = real(t2);
b = imag(t2);
da = real(t);
db = imag(t);
instf = (1/(2*pi))*(a.*db - b.*da)./((a.*a + b.*b)+(abs(t2)==0));
% 1/2pi converts rad/s into cycles/s
% sampling rate already factored in as constant in dwin & ww
% so result is in Hz
F(:,h) = instf(1:(1 + N/2))';
end;