boblsturm@0: function Y = chromagram_IF(d,sr,fftlen,nbin,f_ctr,f_sd)
boblsturm@0: % Y = chromagram_IF(d,sr,fftlen,nbin,f_ctr,f_sd)
boblsturm@0: %  Calculate a "chromagram" of the sound in d (at sampling rate sr)
boblsturm@0: %  Use windows of fftlen points, hopped by ffthop points
boblsturm@0: %  Divide the octave into nbin steps
boblsturm@0: %  Weight with center frequency f_ctr (in Hz) and gaussian SD f_sd
boblsturm@0: %  (in octaves)
boblsturm@0: %  Use instantaneous frequency to keep only real harmonics.
boblsturm@0: % 2006-09-26 dpwe@ee.columbia.edu
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 < 3;   fftlen = 2048; end
boblsturm@0: if nargin < 4;   nbin = 12; end
boblsturm@0: if nargin < 5;   f_ctr = 1000; end
boblsturm@0: if nargin < 6;   f_sd = 1; end
boblsturm@0: 
boblsturm@0: A0 = 27.5; % Hz
boblsturm@0: A440 = 440; % Hz
boblsturm@0: f_ctr_log = log(f_ctr/A0) / log(2);
boblsturm@0: 
boblsturm@0: fminl = octs2hz(hz2octs(f_ctr)-2*f_sd);
boblsturm@0: fminu = octs2hz(hz2octs(f_ctr)-f_sd);
boblsturm@0: fmaxl = octs2hz(hz2octs(f_ctr)+f_sd);
boblsturm@0: fmaxu = octs2hz(hz2octs(f_ctr)+2*f_sd);
boblsturm@0: 
boblsturm@0: ffthop = fftlen/4;
boblsturm@0: nchr = 12;
boblsturm@0: 
boblsturm@0: % Calculate spectrogram and IF gram pitch tracks...
boblsturm@0: [p,m]=ifptrack(d,fftlen,sr,fminl,fminu,fmaxl,fmaxu); 
boblsturm@0: 
boblsturm@0: [nbins,ncols] = size(p);
boblsturm@0: 
boblsturm@0: %disp(['ncols = ',num2str(ncols)]);
boblsturm@0: 
boblsturm@0: % chroma-quantized IF sinusoids
boblsturm@0: Pocts = hz2octs(p+(p==0));
boblsturm@0: Pocts(p(:)==0) = 0;
boblsturm@0: % Figure best tuning alignment
boblsturm@0: nzp = find(p(:)>0);
boblsturm@0: %hist(nchr*Pmapo(nzp)-round(nchr*Pmapo(nzp)),100)
boblsturm@0: [hn,hx] = hist(nchr*Pocts(nzp)-round(nchr*Pocts(nzp)),100);
boblsturm@0: centsoff = hx(find(hn == max(hn)));
boblsturm@0: % Adjust tunings to align better with chroma
boblsturm@0: Pocts(nzp) = Pocts(nzp) - centsoff(1)/nchr;
boblsturm@0: 
boblsturm@0: % Quantize to chroma bins
boblsturm@0: PoctsQ = Pocts;
boblsturm@0: PoctsQ(nzp) = round(nchr*Pocts(nzp))/nchr;
boblsturm@0: 
boblsturm@0: % map IF pitches to chroma bins
boblsturm@0: Pmapc = round(nchr*(PoctsQ - floor(PoctsQ)));
boblsturm@0: Pmapc(p(:) == 0) = -1; 
boblsturm@0: Pmapc(Pmapc(:) == nchr) = 0;
boblsturm@0: 
boblsturm@0: Y = zeros(nchr,ncols);
boblsturm@0: for t = 1:ncols;
boblsturm@0:   Y(:,t)=(repmat([0:(nchr-1)]',1,size(Pmapc,1))==repmat(Pmapc(:,t)',nchr,1))*m(:,t);
boblsturm@0: end