c@1: 
c@1: module cqtkernel;
c@1: 
c@3: vec = load may.vector;
c@3: bf = load may.vector.blockfuncs;
c@3: complex = load may.complex;
c@3: window = load may.signal.window;
c@3: fft = load may.transform.fft;
c@4: pl = load may.plot;
c@6: cm = load may.matrix.complex;
c@3: 
c@2: { pow, round, floor, ceil, nextPowerOfTwo } = load may.mathmisc;
c@1: 
c@9: makeKernel { sampleRate, maxFreq, binsPerOctave } =
c@9:    (q = 1;
c@9:     atomHopFactor = 0.25;
c@9:     thresh = 0.0005;
c@9:     minFreq = (maxFreq/2) * (pow 2 (1/binsPerOctave));
c@9:     bigQ = q / ((pow 2 (1/binsPerOctave)) - 1);
c@1: 
c@9:     maxNK = round(bigQ * sampleRate / minFreq);
c@9:     minNK = round(bigQ * sampleRate /
c@9:                   (minFreq * (pow 2 ((binsPerOctave-1) / binsPerOctave))));
c@1: 
c@9:     atomHop = round(minNK * atomHopFactor);
c@9:     
c@9:     firstCentre = atomHop * (ceil ((ceil (maxNK/2)) / atomHop));
c@9:     
c@9:     fftSize = nextPowerOfTwo (firstCentre + ceil (maxNK/2));
c@9:     
c@9:     println "sampleRate = \(sampleRate), maxFreq = \(maxFreq), binsPerOctave = \(binsPerOctave), q = \(q), atomHopFactor = \(atomHopFactor), thresh = \(thresh)";
c@9:     println "minFreq = \(minFreq), bigQ = \(bigQ), maxNK = \(maxNK), minNK = \(minNK), atomHop = \(atomHop), firstCentre = \(firstCentre), fftSize = \(fftSize)";
c@9:     
c@9:     winNr = floor((fftSize - ceil(maxNK/2) - firstCentre) / atomHop) + 1;
c@9:     
c@9:     lastCentre = firstCentre + (winNr - 1) * atomHop;
c@9:     
c@9:     fftHop = (lastCentre + atomHop) - firstCentre;
c@9:     
c@9:     println "winNr = \(winNr), lastCentre = \(lastCentre), fftHop = \(fftHop)";
c@9:     
c@9:     fftFunc = fft.forward fftSize;
c@9:     
c@9:     // Note the MATLAB uses exp(2*pi*1i*x) for a complex generating
c@9:     // function. We can't do that here; we need to generate real and imag
c@9:     // parts separately as real = cos(2*pi*x), imag = sin(2*pi*x).
c@9:     
c@9:     kernels = map do k:
c@9:     
c@9:         nk = round(bigQ * sampleRate / (minFreq * (pow 2 ((k-1)/binsPerOctave))));
c@9:     
c@9:         // the cq MATLAB toolbox uses a symmetric window for
c@9:         // blackmanharris -- which is odd because it uses a periodic one
c@9:         // for other types. Oh well
c@9:         win = bf.divideBy nk
c@9:            (bf.sqrt
c@9:                (window.windowFunction (BlackmanHarris ()) [Symmetric true] nk));
c@9:         
c@9:         fk = minFreq * (pow 2 ((k-1)/binsPerOctave));
c@9:     
c@9:         genKernel f = bf.multiply win
c@9:            (vec.fromList
c@9:                (map do i: f (2 * pi * fk * i / sampleRate) done [0..nk-1]));
c@9:     
c@9:         reals = genKernel cos;
c@9:         imags = genKernel sin;
c@9:     
c@9:         atomOffset = firstCentre - ceil(nk/2);
c@9:     
c@9:         map do i:
c@9:     
c@9:             shift = vec.zeros (atomOffset + ((i-1) * atomHop));
c@9:     
c@9:             specKernel = fftFunc
c@9:                (complex.complexArray
c@9:                    (vec.concat [shift, reals])
c@9:                    (vec.concat [shift, imags]));
c@9:     
c@9:             map do c:
c@9:                 if complex.magnitude c <= thresh then complex.zero else c fi
c@9:                 done specKernel;
c@9:     
c@9:         done [1..winNr];
c@9:     
c@9:     done [1..binsPerOctave];
c@9:     
c@9:     kmat = cm.toSparse
c@9:        (cm.scaled (1/fftSize)
c@9:            (cm.newComplexMatrix (RowMajor()) (concat kernels)));
c@9:     
c@9:     println "density = \(cm.density kmat)";
c@9:     
c@9:     // Normalisation
c@9:     
c@9:     wx1 = bf.maxindex (complex.magnitudes (cm.getRow 0 kmat));
c@9:     wx2 = bf.maxindex (complex.magnitudes (cm.getRow (cm.height kmat - 1) kmat));
c@9:     
c@9:     subset = cm.columnSlice kmat wx1 (wx2+1);
c@9:     square = cm.product (cm.conjugateTransposed subset) subset;
c@9:     diag = complex.magnitudes (cm.getDiagonal 0 square);
c@9:     wK = vec.slice diag (round(1/q)) (vec.length diag - round(1/q) - 2);
c@9:     
c@9:     weight = (fftHop / fftSize) / (bf.mean (bf.abs wK));
c@9:     weight = sqrt(weight);
c@1: 
c@9:     {
c@9:         kernel = cm.scaled weight kmat,
c@9:         fftSize,
c@9:         fftHop,
c@9:         binsPerOctave,
c@12:         atomsPerFrame = winNr,
c@12:         atomSpacing = atomHop,
c@13:         firstCentre,
c@9:         maxFreq,
c@9:         minFreq,
c@9:         bigQ
c@9:     });
c@1: 
c@9: {
c@9:     makeKernel
c@9: }
c@1: