Chris@2: function cqtKernel = genCQTkernel(fmax, bins, fs, varargin) Chris@2: %Calculating the CQT Kernel for one octave. All atoms are center-stacked. Chris@2: %Atoms are placed so that the stacks of lower octaves are centered at the Chris@2: %same positions in time, however, their amount is reduced by factor two for Chris@2: %each octave down. Chris@2: % Chris@2: %INPUT: Chris@2: % fmax ... highest frequency of interest Chris@2: % bins ... number of bins per octave Chris@2: % fs ... sampling frequency Chris@2: % Chris@2: %optional input parameters (parameter name/value pairs): Chris@2: % Chris@2: % 'q' ... Q scaling factor. Default: 1. Chris@2: % 'atomHopFactor' ... relative hop size corresponding to the shortest Chris@2: % temporal atom. Default: 0.25. Chris@2: % 'thresh' ... values smaller than 'tresh' in the spectral kernel are rounded to Chris@2: % zero. Default: 0.0005. 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: % 'perfRast' ... if set to 1 the kernel is designed in order to Chris@2: % enable perfect rasterization using the function Chris@2: % cqtPerfectRast() (Default: perRast=0). See documentation of Chris@2: % 'cqtPerfectRast' for further information. Chris@2: % Chris@2: %OUTPUT: Chris@2: % cqtKernel ... Structure that contains the spectral kernel 'fKernel' Chris@2: % additional design parameters used in cqt(), cqtPerfectRast() and icqt(). Chris@2: % Chris@2: %Christian Schörkhuber, Anssi Klapuri 2010-06 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'; %default value Chris@2: perfRast = 0; %default value Chris@2: Chris@2: for ain = 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},'win'), winFlag = varargin{ain+1}; end; Chris@2: if strcmp(varargin{ain},'perfRast'), perfRast = varargin{ain+1}; end; Chris@2: end Chris@2: Chris@2: %% define Chris@2: fmin = (fmax/2)*2^(1/bins); Chris@2: Q = 1/(2^(1/bins)-1); Chris@2: Q = Q*q; Chris@2: Nk_max = Q * fs / fmin; Chris@2: Nk_max = round(Nk_max); %length of the largest atom [samples] Chris@2: Chris@2: Chris@2: %% Compute FFT size, FFT hop, atom hop, ... Chris@2: Nk_min = round( Q * fs / (fmin*2^((bins-1)/bins)) ); %length of the shortest atom [samples] Chris@2: atomHOP = round(Nk_min*atomHopFactor); %atom hop size Chris@2: first_center = ceil(Nk_max/2); %first possible center position within the frame Chris@2: first_center = atomHOP * ceil(first_center/atomHOP); %lock the first center to an integer multiple of the atom hop size Chris@2: FFTLen = 2^nextpow2(first_center+ceil(Nk_max/2)); %use smallest possible FFT size (increase sparsity) Chris@2: Chris@2: if perfRast Chris@2: winNr = floor((FFTLen-ceil(Nk_max/2)-first_center)/atomHOP); %number of temporal atoms per FFT Frame Chris@2: if winNr == 0 Chris@2: FFTLen = FFTLen * 2; Chris@2: winNr = floor((FFTLen-ceil(Nk_max/2)-first_center)/atomHOP); Chris@2: end Chris@2: else Chris@2: winNr = floor((FFTLen-ceil(Nk_max/2)-first_center)/atomHOP)+1; %number of temporal atoms per FFT Frame Chris@2: end Chris@2: Chris@2: last_center = first_center + (winNr-1)*atomHOP; Chris@2: fftHOP = (last_center + atomHOP) - first_center; % hop size of FFT frames Chris@2: fftOLP = (FFTLen-fftHOP/FFTLen)*100; %overlap of FFT frames in percent ***AK:needed? Chris@2: Chris@2: %% init variables Chris@2: tempKernel= zeros(1,FFTLen); Chris@2: sparKernel= []; Chris@2: Chris@2: %% Compute kernel Chris@2: atomInd = 0; Chris@2: for k = 1:bins Chris@2: Chris@2: Nk = round( Q * fs / (fmin*2^((k-1)/bins)) ); %N[k] = (fs/fk)*Q. Rounding will be omitted in future versions Chris@2: Chris@2: switch winFlag Chris@2: case 'sqrt_blackmanharris' Chris@2: winFct = sqrt(blackmanharris(Nk)); Chris@2: case 'blackmanharris' Chris@2: winFct = blackmanharris(Nk); Chris@2: case 'sqrt_hann' Chris@2: winFct = sqrt(hann(Nk,'periodic')); Chris@2: case 'hann' Chris@2: winFct = hann(Nk,'periodic'); Chris@2: case 'sqrt_blackman' Chris@2: winFct = sqrt(hann(blackman,'periodic')); Chris@2: case 'blackman' Chris@2: winFct = blackman(Nk,'periodic'); Chris@2: otherwise Chris@2: winFct = sqrt(blackmanharris(Nk)); Chris@2: if k==1, warning('CQT:INPUT','Non-existing window function. Default window is used!'); end; Chris@2: end Chris@2: Chris@2: fk = fmin*2^((k-1)/bins); Chris@2: tempKernelBin = (winFct/Nk) .* exp(2*pi*1i*fk*(0:Nk-1)'/fs); Chris@2: atomOffset = first_center - ceil(Nk/2); Chris@2: Chris@2: for i = 1:winNr Chris@2: shift = atomOffset + ((i-1) * atomHOP); Chris@2: tempKernel(1+shift:Nk+shift) = tempKernelBin; Chris@2: atomInd = atomInd+1; Chris@2: specKernel= fft(tempKernel); Chris@2: specKernel(abs(specKernel)<=thresh)= 0; Chris@2: sparKernel= sparse([sparKernel; specKernel]); Chris@2: tempKernel = zeros(1,FFTLen); %reset window Chris@2: end Chris@2: end Chris@2: sparKernel = (sparKernel.')/FFTLen; Chris@2: Chris@2: %% Normalize the magnitudes of the atoms Chris@2: [ignore,wx1]=max(sparKernel(:,1)); Chris@2: [ignore,wx2]=max(sparKernel(:,end)); Chris@2: wK=sparKernel(wx1:wx2,:); Chris@2: wK = diag(wK * wK'); Chris@2: wK = wK(round(1/q)+1:(end-round(1/q)-2)); Chris@2: weight = 1./mean(abs(wK)); Chris@2: weight = weight.*(fftHOP/FFTLen); Chris@2: weight = sqrt(weight); %sqrt because the same weight is applied in icqt again Chris@2: sparKernel = weight.*sparKernel; Chris@2: Chris@2: %% return Chris@2: cqtKernel = struct('fKernel',sparKernel,'fftLEN',FFTLen,'fftHOP',fftHOP,'fftOverlap',fftOLP,'perfRast',perfRast,... Chris@2: 'bins',bins,'firstcenter',first_center,'atomHOP',atomHOP,'atomNr',winNr,'Nk_max',Nk_max,'Q',Q,'fmin',fmin);