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@1: fs = 48000;
c@1: 
c@1: fmax = fs/2;
c@1: 
c@1: bins = 24;
c@1: 
c@1: q = 1;
c@1: 
c@1: atomHopFactor = 0.25;
c@1: 
c@1: thresh = 0.0005;
c@1: 
c@1: fmin = (fmax/2) * (pow 2 (1/bins));
c@1: 
c@1: bigQ = q / ((pow 2 (1/bins)) - 1);
c@1: 
c@1: nk_max = round(bigQ * fs / fmin);
c@1: 
c@1: nk_min = round(bigQ * fs / (fmin * (pow 2 ((bins-1)/bins))));
c@1: 
c@1: atomHop = round(nk_min * atomHopFactor);
c@1: 
c@1: first_center = atomHop * Math#ceil(Math#ceil(nk_max/2) / atomHop);
c@1: 
c@1: fftLen = nextPowerOfTwo (first_center + Math#ceil(nk_max/2));
c@1: 
c@1: println "fs = \(fs), fmax = \(fmax), bins = \(bins), q = \(q), atomHopFactor = \(atomHopFactor), thresh = \(thresh)";
c@1: 
c@1: println "fmin = \(fmin), bigQ = \(bigQ), nk_max = \(nk_max), nk_min = \(nk_min), atomHop = \(atomHop), first_center = \(first_center), fftLen = \(fftLen)";
c@1: 
c@2: winNr = floor((fftLen - ceil(nk_max/2) - first_center) / atomHop) + 1;
c@2: 
c@2: last_center = first_center + (winNr - 1) * atomHop;
c@2: 
c@2: fftHop = (last_center + atomHop) - first_center;
c@2: 
c@2: fftOverlap = ((fftLen - fftHop) / fftLen) * 100; // as % -- why? just for diagnostics?
c@2: 
c@2: println "winNr = \(winNr), last_center = \(last_center), fftHop = \(fftHop), fftOverlap = \(fftOverlap)%";
c@2: 
c@3: fftFunc = fft.forward fftLen;
c@2: 
c@3: // Note the MATLAB uses exp(2*pi*1i*x) for a complex generating
c@3: // function. We can't do that here; we need to generate real and imag
c@3: // parts separately as real = cos(2*pi*x), imag = sin(2*pi*x).
c@2: 
c@3: kernels = map do k:
c@3: 
c@3:     nk = round(bigQ * fs / (fmin * (pow 2 ((k-1)/bins))));
c@3:     println "k = \(k) -> nk = \(nk)";
c@3:     
c@3:     win = bf.divideBy nk (bf.sqrt (window.blackmanHarris nk));
c@3:     
c@3:     fk = fmin * (pow 2 ((k-1)/bins));
c@3: 
c@3:     println "fk = \(fk)";
c@3: 
c@3:     genKernel f = bf.multiply win
c@3:        (vec.fromList (map do i: f (2 * pi * fk * i / fs) done [0..nk-1]));
c@3: 
c@3:     reals = genKernel cos;
c@3:     imags = genKernel sin;
c@3: 
c@3:     atomOffset = first_center - ceil(nk/2);
c@3: 
c@3:     map do i:
c@3: 
c@3:         shift = atomOffset + ((i-1) * atomHop);
c@3: 
c@4:         println "shift = \(shift)";
c@4: 
c@6:         shiftedReals = vec.concat [vec.zeros shift, reals];
c@6:         shiftedImags = vec.concat [vec.zeros shift, imags];
c@6: 
c@6: //        println "shiftedReals = \(vec.list shiftedReals)";
c@6: //        println "shiftedImags = \(vec.list shiftedImags)";
c@6: 
c@3:         specKernel = fftFunc
c@3:            (complex.complexArray
c@3:                (vec.concat [vec.zeros shift, reals])
c@3:                (vec.concat [vec.zeros shift, imags]));
c@3: 
c@3:         map do c:
c@6:             if complex.magnitude c < thresh then complex.zero else c fi
c@3:             done specKernel;
c@4: 
c@3:     done [1..winNr];
c@3: 
c@3: done [1..bins];
c@3: 
c@6: kmat = cm.toSparse
c@6:    (cm.scaled (1/fftLen)
c@6:        (cm.newComplexMatrix (ColumnMajor()) (concat kernels))); // or row major?
c@2: 
c@6: println "density = \(cm.density kmat)";
c@1: 
c@6: kmat;
c@5: 
c@6: 
c@6: