Chris@2: function Xcqt = cqt(x,fmin,fmax,bins,fs,varargin) 
Chris@2: %Xcqt = cqt(x,fmin,fmax,bins,fs,varargin) calculates the constant-Q transform of the input signal x.
Chris@2: %
Chris@2: %INPUT:
Chris@2: %   fmin      ... lowest frequency of interest
Chris@2: %   fmax      ... highest frequency of interest
Chris@2: %   bins      ... frequency bins per octave
Chris@2: %   fs        ... sampling rate
Chris@2: %
Chris@2: %   optional input parameters (parameter name/value pairs):
Chris@2: %
Chris@2: %   'atomHopFactor'   ... overlap of temporal atoms in percent. Default: 0.25.
Chris@2: %    
Chris@2: %   'q'       ... the maximum value for optimal reconstruction is q=1.
Chris@2: %                 For values smaller than 1 the bandwidths of the spectral
Chris@2: %                 atoms (filter) are increased retaining their center
Chris@2: %                 frequencies (frequency 'smearing', frequency domain redundancy 
Chris@2: %                 increases, time resolutin improves). Default: 1.
Chris@2: %   'thresh'  ... all values in the cqt kernel smaller than tresh are
Chris@2: %                 rounded to zero. A high value for thresh yields a
Chris@2: %                 very sparse kernel (fast) but introduces a bigger error. 
Chris@2: %                 The default value is chosen so that the error due to rounding is negligible.
Chris@2: %   'kernel'  ... if the cqt kernel structure has been precomputed
Chris@2: %                 (using function 'genCQTkernel'), the computation of the kernel
Chris@2: %                 will be by-passed below).
Chris@2: %   'win'     ... defines which window will be used for the CQT. Valid
Chris@2: %                 values are: 'blackman','hann' and 'blackmanharris'. To
Chris@2: %                 use the square root of each window use the prefix 'sqrt_'
Chris@2: %                 (i.e. 'sqrt_blackman'). Default: 'sqrt_blackmanharris'
Chris@2: %   'coeffB',
Chris@2: %   'coeffA'  ... Filter coefficients for the anti-aliasing filter, where
Chris@2: %                 'coeffB' is the numerator and 'coeffA' is the
Chris@2: %                 denominator (listed in descending powers of z). 
Chris@2: %                                                  
Chris@2: %OUTPUT:
Chris@2: %   Xcqt      ... struct that comprises various fields: 
Chris@2: %              spCQT: CQT coefficients in the form of a sparse matrix 
Chris@2: %                    (rasterized, not interpolated)
Chris@2: %              fKernel: spectral Kernel 
Chris@2: %              fmin: frequency of the lowest bin
Chris@2: %              fmax: frequency of the hiqhest bin
Chris@2: %              octaveNr: number of octaves processed
Chris@2: %              bins: number of bins per octave
Chris@2: %              intParams: structure containing additional parameters for the inverse transform   
Chris@2: %
Chris@2: %Christian Schörkhuber, Anssi Klapuri 2010-06
Chris@2: 
Chris@2: %% input checking
Chris@2: if size(x,2) > 1 && size(x,1) > 1, error('cqt requires one-dimensional input!'); end;
Chris@2: if size(x,2) > 1, x = x(:); end; %column vector
Chris@2: 
Chris@2: %% input parameters
Chris@2: q = 1; %default value
Chris@2: atomHopFactor = 0.25; %default value
Chris@2: thresh = 0.0005; %default value
Chris@2: winFlag = 'sqrt_blackmanharris';
Chris@2: 
Chris@2: for ain = 1:1:length(varargin)
Chris@2:     if strcmp(varargin{ain},'q'), q = varargin{ain+1}; end;
Chris@2:     if strcmp(varargin{ain},'atomHopFactor'), atomHopFactor = varargin{ain+1}; end;
Chris@2:     if strcmp(varargin{ain},'thresh'), thresh = varargin{ain+1}; end;
Chris@2:     if strcmp(varargin{ain},'kernel'), cqtKernel = varargin{ain+1}; end;
Chris@2:     if strcmp(varargin{ain},'win'), winFlag = varargin{ain+1}; end;
Chris@2:     if strcmp(varargin{ain},'coeffB'), B = varargin{ain+1}; end;
Chris@2:     if strcmp(varargin{ain},'coeffA'), A = varargin{ain+1}; end;    
Chris@2: end
Chris@2: 
Chris@2: %% define
Chris@2: octaveNr = ceil(log2(fmax/fmin));
Chris@2: xlen_init = length(x);
Chris@2: 
Chris@2: %% design lowpass filter
Chris@2: if ~exist('B','var') || ~exist('A','var')
Chris@2:     LPorder = 6; %order of the anti-aliasing filter
Chris@2:     cutoff = 0.5;
Chris@2:     [B A] = butter(LPorder,cutoff,'low'); %design f_nyquist/2-lowpass filter
Chris@2: end
Chris@2: 
Chris@2: %% design kernel for one octave 
Chris@2: if ~exist('cqtKernel','var')
Chris@2:     cqtKernel = genCQTkernel(fmax, bins,fs,'q',q,'atomHopFactor',atomHopFactor,'thresh',thresh,'win',winFlag);
Chris@2: end
Chris@2: 
Chris@2: %% calculate CQT
Chris@2: cellCQT = cell(1,octaveNr);
Chris@2: maxBlock = cqtKernel.fftLEN * 2^(octaveNr-1); %largest FFT Block (virtual)
Chris@2: suffixZeros = maxBlock; 
Chris@2: prefixZeros = maxBlock;
Chris@2: x = [zeros(prefixZeros,1); x; zeros(suffixZeros,1)]; %zeropadding
Chris@2: OVRLP = cqtKernel.fftLEN - cqtKernel.fftHOP;
Chris@2: K = cqtKernel.fKernel'; %conjugate spectral kernel for cqt transformation
Chris@2: 
Chris@2: for i = 1:octaveNr
Chris@2:     xx = buffer(x,cqtKernel.fftLEN, OVRLP,'nodelay'); %generating FFT blocks
Chris@2:     XX = fft(xx); %applying fft to each column (each FFT frame)
Chris@2:     cellCQT{i} = K*XX; %calculating cqt coefficients for all FFT frames for this octave
Chris@2: 
Chris@2:     if i~=octaveNr
Chris@2:         x = filtfilt(B,A,x); %anti aliasing filter
Chris@2:         x = x(1:2:end); %drop samplerate by 2
Chris@2:     end
Chris@2: end
Chris@2: 
Chris@2: %% map to sparse matrix representation
Chris@2: spCQT = cell2sparse(cellCQT,octaveNr,bins,cqtKernel.firstcenter,cqtKernel.atomHOP,cqtKernel.atomNr);
Chris@2: 
Chris@2: %% return
Chris@2: intParam = struct('sufZeros',suffixZeros,'preZeros',prefixZeros,'xlen_init',xlen_init,'fftLEN',cqtKernel.fftLEN,'fftHOP',cqtKernel.fftHOP,...
Chris@2:     'q',q,'filtCoeffA',A,'filtCoeffB',B,'firstcenter',cqtKernel.firstcenter,'atomHOP',cqtKernel.atomHOP,...
Chris@2:     'atomNr',cqtKernel.atomNr,'Nk_max',cqtKernel.Nk_max,'Q',cqtKernel.Q,'rast',0);
Chris@2: 
Chris@2: Xcqt = struct('spCQT',spCQT,'fKernel',cqtKernel.fKernel,'fmax',fmax,'fmin',fmin*2^(1/bins),'octaveNr',octaveNr,'bins',cqtKernel.bins,'intParams',intParam);
Chris@2: 
Chris@2: