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