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: $
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);
% 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;