annotate mirex2012-matlab/genCQTkernel.m @ 116:91bb029a847a timing

Reorder the calculations to match the series of vector operations in the most recent bqvec code, just in case it's the order of vector calculations that is saving the time rather than the avoidance of std::vector

function cqtKernel = genCQTkernel(fmax, bins, fs, varargin)
%Calculating the CQT Kernel for one octave. All atoms are center-stacked.
%Atoms are placed so that the stacks of lower octaves are centered at the
%same positions in time, however, their amount is reduced by factor two for
%each octave down.
%
%INPUT:
%   fmax          ... highest frequency of interest
%   bins          ... number of bins per octave
%   fs            ... sampling frequency
%
%optional input parameters (parameter name/value pairs):
%
%   'q'             ... Q scaling factor. Default: 1.
%   'atomHopFactor' ... relative hop size corresponding to the shortest
%                       temporal atom. Default: 0.25.
%   'thresh'        ... values smaller than 'tresh' in the spectral kernel are rounded to
%                       zero. Default: 0.0005.
%   'win'           ... defines which window will be used for the CQT. Valid
%                       values are: 'blackman','hann' and 'blackmanharris'. To
%                       use the square root of each window use the prefix 'sqrt_'
%                      (i.e. 'sqrt_blackman'). Default: 'sqrt_blackmanharris'
%   'perfRast'      ... if set to 1 the kernel is designed in order to
%                       enable perfect rasterization using the function
%                       cqtPerfectRast() (Default: perRast=0). See documentation of
%                       'cqtPerfectRast' for further information.
%
%OUTPUT:
%   cqtKernel   ... Structure that contains the spectral kernel 'fKernel'
%                   additional design parameters used in cqt(), cqtPerfectRast() and icqt().
%
%Christian Schörkhuber, Anssi Klapuri 2010-06

%% input parameters
q = 1; %default value
atomHopFactor = 0.25; %default value
thresh = 0.0005; %default value
winFlag = 'sqrt_blackmanharris'; %default value
perfRast = 0; %default value

for ain = 1:length(varargin)
    if strcmp(varargin{ain},'q'), q = varargin{ain+1}; end;
    if strcmp(varargin{ain},'atomHopFactor'), atomHopFactor = varargin{ain+1}; end;
    if strcmp(varargin{ain},'thresh'), thresh = varargin{ain+1}; end;
    if strcmp(varargin{ain},'win'), winFlag = varargin{ain+1}; end;
    if strcmp(varargin{ain},'perfRast'), perfRast = varargin{ain+1}; end;
end

%% define
fmin = (fmax/2)*2^(1/bins);
Q = 1/(2^(1/bins)-1);
Q = Q*q;
Nk_max = Q * fs / fmin; 
Nk_max = round(Nk_max); %length of the largest atom [samples]


%% Compute FFT size, FFT hop, atom hop, ...
Nk_min = round( Q * fs / (fmin*2^((bins-1)/bins)) ); %length of the shortest atom [samples]
atomHOP = round(Nk_min*atomHopFactor); %atom hop size
first_center = ceil(Nk_max/2); %first possible center position within the frame
first_center = atomHOP * ceil(first_center/atomHOP); %lock the first center to an integer multiple of the atom hop size
FFTLen = 2^nextpow2(first_center+ceil(Nk_max/2)); %use smallest possible FFT size (increase sparsity)

if perfRast
    winNr = floor((FFTLen-ceil(Nk_max/2)-first_center)/atomHOP); %number of temporal atoms per FFT Frame
    if winNr == 0
        FFTLen = FFTLen * 2;
        winNr = floor((FFTLen-ceil(Nk_max/2)-first_center)/atomHOP);
    end
else
    winNr = floor((FFTLen-ceil(Nk_max/2)-first_center)/atomHOP)+1; %number of temporal atoms per FFT Frame
end

last_center = first_center + (winNr-1)*atomHOP;
fftHOP = (last_center + atomHOP) - first_center; % hop size of FFT frames
fftOLP = (FFTLen-fftHOP/FFTLen)*100; %overlap of FFT frames in percent ***AK:needed?

%% init variables
tempKernel= zeros(1,FFTLen); 
sparKernel= []; 

%% Compute kernel
atomInd = 0;
for k = 1:bins
    
   Nk = round( Q * fs / (fmin*2^((k-1)/bins)) ); %N[k] = (fs/fk)*Q. Rounding will be omitted in future versions
   
   switch winFlag
       case 'sqrt_blackmanharris'
            winFct = sqrt(blackmanharris(Nk));
       case 'blackmanharris'
            winFct = blackmanharris(Nk);
       case 'sqrt_hann'
            winFct = sqrt(hann(Nk,'periodic'));
       case 'hann'
            winFct = hann(Nk,'periodic');
       case 'sqrt_blackman'
            winFct = sqrt(hann(blackman,'periodic'));
       case 'blackman'
            winFct = blackman(Nk,'periodic');
       otherwise
            winFct = sqrt(blackmanharris(Nk));
            if k==1, warning('CQT:INPUT','Non-existing window function. Default window is used!'); end;
   end
   
   fk = fmin*2^((k-1)/bins);
   tempKernelBin = (winFct/Nk) .* exp(2*pi*1i*fk*(0:Nk-1)'/fs);
   atomOffset = first_center - ceil(Nk/2);

   for i = 1:winNr 
       shift = atomOffset + ((i-1) * atomHOP);
       tempKernel(1+shift:Nk+shift) = tempKernelBin; 
       atomInd = atomInd+1;
       specKernel= fft(tempKernel);
       specKernel(abs(specKernel)<=thresh)= 0;   
       sparKernel= sparse([sparKernel; specKernel]);
       tempKernel = zeros(1,FFTLen); %reset window     
   end
end
sparKernel = (sparKernel.')/FFTLen;

%% Normalize the magnitudes of the atoms
[ignore,wx1]=max(sparKernel(:,1));
[ignore,wx2]=max(sparKernel(:,end));
wK=sparKernel(wx1:wx2,:);
wK = diag(wK * wK');
wK = wK(round(1/q)+1:(end-round(1/q)-2));
weight = 1./mean(abs(wK));
weight = weight.*(fftHOP/FFTLen); 
weight = sqrt(weight); %sqrt because the same weight is applied in icqt again
sparKernel = weight.*sparKernel;

%% return
cqtKernel = struct('fKernel',sparKernel,'fftLEN',FFTLen,'fftHOP',fftHOP,'fftOverlap',fftOLP,'perfRast',perfRast,...
    'bins',bins,'firstcenter',first_center,'atomHOP',atomHOP,'atomNr',winNr,'Nk_max',Nk_max,'Q',Q,'fmin',fmin);