|
boblsturm@0
|
1 function [F,D] = ifgram(X, N, W, H, SR)
|
|
boblsturm@0
|
2 % [F,D] = ifgram(X, N, W, H, SR) Instantaneous frequency by phase deriv.
|
|
boblsturm@0
|
3 % X is a 1-D signal. Process with N-point FFTs applying a W-point
|
|
boblsturm@0
|
4 % window, stepping by H points; return (N/2)+1 channels with the
|
|
boblsturm@0
|
5 % instantaneous frequency (as a proportion of the sampling rate)
|
|
boblsturm@0
|
6 % obtained as the time-derivative of the phase of the complex spectrum
|
|
boblsturm@0
|
7 % as described by Toshihiro Abe et al in ICASSP'95, Eurospeech'97
|
|
boblsturm@0
|
8 % Same arguments and some common code as dpwebox/stft.m.
|
|
boblsturm@0
|
9 % Calculates regular STFT as side effect - returned in D.
|
|
boblsturm@0
|
10 % after 1998may02 dpwe@icsi.berkeley.edu
|
|
boblsturm@0
|
11 % 2001-03-05 dpwe@ee.columbia.edu revised version
|
|
boblsturm@0
|
12 % 2001-12-13 dpwe@ee.columbia.edu Fixed to work when N != W
|
|
boblsturm@0
|
13 % $Header: $
|
|
boblsturm@0
|
14
|
|
boblsturm@0
|
15 % Copyright (c) 2006 Columbia University.
|
|
boblsturm@0
|
16 %
|
|
boblsturm@0
|
17 % This file is part of LabROSA-coversongID
|
|
boblsturm@0
|
18 %
|
|
boblsturm@0
|
19 % LabROSA-coversongID is free software; you can redistribute it and/or modify
|
|
boblsturm@0
|
20 % it under the terms of the GNU General Public License version 2 as
|
|
boblsturm@0
|
21 % published by the Free Software Foundation.
|
|
boblsturm@0
|
22 %
|
|
boblsturm@0
|
23 % LabROSA-coversongID is distributed in the hope that it will be useful, but
|
|
boblsturm@0
|
24 % WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
boblsturm@0
|
25 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
boblsturm@0
|
26 % General Public License for more details.
|
|
boblsturm@0
|
27 %
|
|
boblsturm@0
|
28 % You should have received a copy of the GNU General Public License
|
|
boblsturm@0
|
29 % along with LabROSA-coversongID; if not, write to the Free Software
|
|
boblsturm@0
|
30 % Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
|
|
boblsturm@0
|
31 % 02110-1301 USA
|
|
boblsturm@0
|
32 %
|
|
boblsturm@0
|
33 % See the file "COPYING" for the text of the license.
|
|
boblsturm@0
|
34
|
|
boblsturm@0
|
35 if nargin < 2; N = 256; end
|
|
boblsturm@0
|
36 if nargin < 3; W = N; end
|
|
boblsturm@0
|
37 if nargin < 4; H = W/2; end
|
|
boblsturm@0
|
38 if nargin < 5; SR = 1; end
|
|
boblsturm@0
|
39
|
|
boblsturm@0
|
40 s = length(X);
|
|
boblsturm@0
|
41 % Make sure it's a single row
|
|
boblsturm@0
|
42 if size(X,1) > 1
|
|
boblsturm@0
|
43 X = X';
|
|
boblsturm@0
|
44 end
|
|
boblsturm@0
|
45
|
|
boblsturm@0
|
46 %win = [0,hanning(W-1)'];
|
|
boblsturm@0
|
47 win = 0.5*(1-cos([0:(W-1)]/W*2*pi));
|
|
boblsturm@0
|
48
|
|
boblsturm@0
|
49 % Window for discrete differentiation
|
|
boblsturm@0
|
50 T = W/SR;
|
|
boblsturm@0
|
51 dwin = -pi / T * sin([0:(W-1)]/W*2*pi);
|
|
boblsturm@0
|
52
|
|
boblsturm@0
|
53 % sum(win) takes out integration due to window, 2 compensates for neg frq
|
|
boblsturm@0
|
54 norm = 2/sum(win);
|
|
boblsturm@0
|
55
|
|
boblsturm@0
|
56 % How many complete windows?
|
|
boblsturm@0
|
57 nhops = 1 + floor((s - W)/H);
|
|
boblsturm@0
|
58
|
|
boblsturm@0
|
59 F = zeros(1 + N/2, nhops);
|
|
boblsturm@0
|
60 D = zeros(1 + N/2, nhops);
|
|
boblsturm@0
|
61
|
|
boblsturm@0
|
62 nmw1 = floor( (N-W)/2 );
|
|
boblsturm@0
|
63 nmw2 = N-W - nmw1;
|
|
boblsturm@0
|
64
|
|
boblsturm@0
|
65 ww = 2*pi*[0:(N-1)]*SR/N;
|
|
boblsturm@0
|
66
|
|
boblsturm@0
|
67 for h = 1:nhops
|
|
boblsturm@0
|
68 u = X((h-1)*H + [1:W]);
|
|
boblsturm@0
|
69 % if(h==0)
|
|
boblsturm@0
|
70 % plot(u)
|
|
boblsturm@0
|
71 % end
|
|
boblsturm@0
|
72 % Apply windows now, while the length is right
|
|
boblsturm@0
|
73 wu = win.*u;
|
|
boblsturm@0
|
74 du = dwin.*u;
|
|
boblsturm@0
|
75
|
|
boblsturm@0
|
76 % Pad or truncate samples if N != W
|
|
boblsturm@0
|
77 if N > W
|
|
boblsturm@0
|
78 wu = [zeros(1,nmw1),wu,zeros(1,nmw2)];
|
|
boblsturm@0
|
79 du = [zeros(1,nmw1),du,zeros(1,nmw2)];
|
|
boblsturm@0
|
80 end
|
|
boblsturm@0
|
81 if N < W
|
|
boblsturm@0
|
82 wu = wu(-nmw1+[1:N]);
|
|
boblsturm@0
|
83 du = du(-nmw1+[1:N]);
|
|
boblsturm@0
|
84 end
|
|
boblsturm@0
|
85 % FFTs of straight samples plus differential-weighted ones
|
|
boblsturm@0
|
86 t1 = fft(fftshift(du));
|
|
boblsturm@0
|
87 t2 = fft(fftshift(wu));
|
|
boblsturm@0
|
88 % Scale down to factor out length & window effects
|
|
boblsturm@0
|
89 D(:,h) = t2(1:(1 + N/2))'*norm;
|
|
boblsturm@0
|
90
|
|
boblsturm@0
|
91 % Calculate instantaneous frequency from phase of differential spectrum
|
|
boblsturm@0
|
92 t = t1 + j*(ww.*t2);
|
|
boblsturm@0
|
93 a = real(t2);
|
|
boblsturm@0
|
94 b = imag(t2);
|
|
boblsturm@0
|
95 da = real(t);
|
|
boblsturm@0
|
96 db = imag(t);
|
|
boblsturm@0
|
97 instf = (1/(2*pi))*(a.*db - b.*da)./((a.*a + b.*b)+(abs(t2)==0));
|
|
boblsturm@0
|
98 % 1/2pi converts rad/s into cycles/s
|
|
boblsturm@0
|
99 % sampling rate already factored in as constant in dwin & ww
|
|
boblsturm@0
|
100 % so result is in Hz
|
|
boblsturm@0
|
101
|
|
boblsturm@0
|
102 F(:,h) = instf(1:(1 + N/2))';
|
|
boblsturm@0
|
103
|
|
boblsturm@0
|
104 end;
|
|
boblsturm@0
|
105
|
