view yeti/cqtkernel.yeti @ 5:6b59109469c4

Return sparse complex matrix

module cqtkernel;

vec = load may.vector;
bf = load may.vector.blockfuncs;
complex = load may.complex;
window = load may.signal.window;
fft = load may.transform.fft;
pl = load may.plot;
complexmatrix = load may.matrix.complex;

{ pow, round, floor, ceil, nextPowerOfTwo } = load may.mathmisc;

fs = 48000;

fmax = fs/2;

bins = 24;

q = 1;

atomHopFactor = 0.25;

thresh = 0.0005;

fmin = (fmax/2) * (pow 2 (1/bins));

bigQ = q / ((pow 2 (1/bins)) - 1);

nk_max = round(bigQ * fs / fmin);

nk_min = round(bigQ * fs / (fmin * (pow 2 ((bins-1)/bins))));

atomHop = round(nk_min * atomHopFactor);

first_center = atomHop * Math#ceil(Math#ceil(nk_max/2) / atomHop);

fftLen = nextPowerOfTwo (first_center + Math#ceil(nk_max/2));

println "fs = \(fs), fmax = \(fmax), bins = \(bins), q = \(q), atomHopFactor = \(atomHopFactor), thresh = \(thresh)";

println "fmin = \(fmin), bigQ = \(bigQ), nk_max = \(nk_max), nk_min = \(nk_min), atomHop = \(atomHop), first_center = \(first_center), fftLen = \(fftLen)";

winNr = floor((fftLen - ceil(nk_max/2) - first_center) / atomHop) + 1;

last_center = first_center + (winNr - 1) * atomHop;

fftHop = (last_center + atomHop) - first_center;

fftOverlap = ((fftLen - fftHop) / fftLen) * 100; // as % -- why? just for diagnostics?

println "winNr = \(winNr), last_center = \(last_center), fftHop = \(fftHop), fftOverlap = \(fftOverlap)%";

fftFunc = fft.forward fftLen;

// Note the MATLAB uses exp(2*pi*1i*x) for a complex generating
// function. We can't do that here; we need to generate real and imag
// parts separately as real = cos(2*pi*x), imag = sin(2*pi*x).

kernels = map do k:

    nk = round(bigQ * fs / (fmin * (pow 2 ((k-1)/bins))));
    println "k = \(k) -> nk = \(nk)";
    
    win = bf.divideBy nk (bf.sqrt (window.blackmanHarris nk));
    
    fk = fmin * (pow 2 ((k-1)/bins));

    println "fk = \(fk)";

    genKernel f = bf.multiply win
       (vec.fromList (map do i: f (2 * pi * fk * i / fs) done [0..nk-1]));

    reals = genKernel cos;
    imags = genKernel sin;

    atomOffset = first_center - ceil(nk/2);

    map do i:

        shift = atomOffset + ((i-1) * atomHop);

        println "shift = \(shift)";

        specKernel = fftFunc
           (complex.complexArray
               (vec.concat [vec.zeros shift, reals])
               (vec.concat [vec.zeros shift, imags]));

        map do c:
            if complex.magnitude c < thresh then complex.zero
            else complex.scale (1/fftLen) c fi
            done specKernel;

    done [1..winNr];

done [1..bins];

kmat = complexmatrix.newComplexMatrix (ColumnMajor()) (concat kernels); // or row major?
println "density = \(complexmatrix.density kmat)";
complexmatrix.toSparse kmat;