boblsturm@0: function wts = fft2chromamx(nfft,nbins,sr,A440,ctroct,octwidth)
boblsturm@0: % wts = fft2chromamx(nfft,nbins,sr,A440,ctroct,octwidth)
boblsturm@0: %     Create a wts matrix to convert FFT to Chroma
boblsturm@0: %     A440 is optional ref frq for A
boblsturm@0: %     ctroct, octwidth specify a dominance window - Gaussian
boblsturm@0: %     weighting centered on ctroct (in octs, re A0 = 27.5Hz) and 
boblsturm@0: %     with a gaussian half-width of octwidth.  Defaults to
boblsturm@0: %     halfwidth = inf i.e. flat.
boblsturm@0: %     2006-06-29 dpwe@ee.columbia.edu
boblsturm@0: 
boblsturm@0: if nargin < 2;   nbins = 12;  end
boblsturm@0: if nargin < 3;   sr = 22050; end
boblsturm@0: if nargin < 4;   A440 = 440; end
boblsturm@0: if nargin < 5;   ctroct = 5; end
boblsturm@0: if nargin < 6;   octwidth = 0; end
boblsturm@0: 
boblsturm@0: wts = zeros(nbins, nfft);
boblsturm@0: 
boblsturm@0: fftfrqbins = nbins*hz2octs([1:(nfft-1)]/nfft*sr,A440);
boblsturm@0: 
boblsturm@0: % make up a value for the 0 Hz bin = 1.5 octaves below bin 1
boblsturm@0: % (so chroma is 50% rotated from bin 1, and bin width is broad)
boblsturm@0: fftfrqbins = [fftfrqbins(1)-1.5*nbins,fftfrqbins];
boblsturm@0: 
boblsturm@0: binwidthbins = [max(1, fftfrqbins(2:nfft) - fftfrqbins(1:(nfft-1))), 1];
boblsturm@0: 
boblsturm@0: D = repmat(fftfrqbins,nbins,1) - repmat([0:(nbins-1)]',1,nfft);
boblsturm@0: 
boblsturm@0: nbins2 = round(nbins/2);
boblsturm@0: 
boblsturm@0: % Project into range -nbins/2 .. nbins/2
boblsturm@0: % add on fixed offset of 10*nbins to ensure all values passed to rem are +ve
boblsturm@0: D = rem(D + nbins2 + 10*nbins, nbins) - nbins2;
boblsturm@0: 
boblsturm@0: % Gaussian bumps - 2*D to make them narrower
boblsturm@0: wts = exp(-0.5*(2*D./repmat(binwidthbins,nbins,1)).^2);
boblsturm@0: 
boblsturm@0: % normalize each column
boblsturm@0: wts = wts./repmat(sqrt(sum(wts.^2)),nbins,1);
boblsturm@0: 
boblsturm@0: % remove aliasing columns
boblsturm@0: wts(:,[(nfft/2+2):nfft]) = 0;
boblsturm@0: 
boblsturm@0: % Maybe apply scaling for fft bins
boblsturm@0: if octwidth > 0
boblsturm@0:   wts = wts.*repmat(exp(-0.5*(((fftfrqbins/nbins - ctroct)/octwidth).^2)), nbins, 1);
boblsturm@0: end
boblsturm@0: 
boblsturm@0: %wts = binwidthbins;
boblsturm@0: %wts = fftfrqbins;
boblsturm@0: 
boblsturm@0: function octs = hz2octs(freq, A440)
boblsturm@0: % octs = hz2octs(freq, A440)
boblsturm@0: % Convert a frequency in Hz into a real number counting 
boblsturm@0: % the octaves above A0. So hz2octs(440) = 4.0
boblsturm@0: % Optional A440 specifies the Hz to be treated as middle A (default 440).
boblsturm@0: % 2006-06-29 dpwe@ee.columbia.edu for fft2chromamx
boblsturm@0: 
boblsturm@0: %if nargin < 2;   A440 = 440; end
boblsturm@0: 
boblsturm@0: % A4 = A440 = 440 Hz, so A0 = 440/16 Hz
boblsturm@0: octs = log(freq./(A440/16))./log(2);
boblsturm@0: 
boblsturm@0: