diff src/matlab/ifgram.m @ 0:c52bc3e8d3ad tip

user: boblsturm branch 'default' added README.md added assets/.DS_Store added assets/playButton.jpg added assets/stopButton.png added assets/swapButton.jpg added data/.DS_Store added data/fiveoctaves.mp3 added data/glock2.wav added data/sinScale.mp3 added data/speech_female.mp3 added data/sweep.wav added nimfks.m.lnk added src/.DS_Store added src/matlab/.DS_Store added src/matlab/AnalysisCache.m added src/matlab/CSS.m added src/matlab/DataHash.m added src/matlab/ExistsInCache.m added src/matlab/KLDivCost.m added src/matlab/LoadFromCache.m added src/matlab/SA_B_NMF.m added src/matlab/SaveInCache.m added src/matlab/Sound.m added src/matlab/SynthesisCache.m added src/matlab/chromagram_E.m added src/matlab/chromagram_IF.m added src/matlab/chromagram_P.m added src/matlab/chromsynth.m added src/matlab/computeSTFTFeat.m added src/matlab/controller.m added src/matlab/decibelSliderReleaseCallback.m added src/matlab/drawClickCallBack.m added src/matlab/fft2chromamx.m added src/matlab/hz2octs.m added src/matlab/ifgram.m added src/matlab/ifptrack.m added src/matlab/istft.m added src/matlab/nimfks.fig added src/matlab/nimfks.m added src/matlab/nmfFn.m added src/matlab/nmf_beta.m added src/matlab/nmf_divergence.m added src/matlab/nmf_euclidean.m added src/matlab/prune_corpus.m added src/matlab/rot_kernel.m added src/matlab/templateAdditionResynth.m added src/matlab/templateDelCb.m added src/matlab/templateScrollCb.m
author boblsturm
date Sun, 18 Jun 2017 06:26:13 -0400
parents
children
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/matlab/ifgram.m	Sun Jun 18 06:26:13 2017 -0400
@@ -0,0 +1,105 @@
+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;
+