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