comparison yeti/cqtkernel.yeti @ 23:4f2f28d5df58

Fix frequency calculation in cpp, and some diagnostics

     // parts separately as real = cos(2*pi*x), imag = sin(2*pi*x).
     
     kernels = map do k:
     
         nk = round(bigQ * sampleRate / (minFreq * (pow 2 ((k-1)/binsPerOctave))));
-    
+
         // the cq MATLAB toolbox uses a symmetric window for
         // blackmanharris -- which is odd because it uses a periodic one
         // for other types. Oh well
         win = bf.divideBy nk
            (bf.sqrt
                (window.windowFunction (BlackmanHarris ()) [Symmetric true] nk));
-        
+
         fk = minFreq * (pow 2 ((k-1)/binsPerOctave));
-    
+
         genKernel f = bf.multiply win
            (vec.fromList
                (map do i: f (2 * pi * fk * i / sampleRate) done [0..nk-1]));
     
         reals = genKernel cos;

     
             specKernel = fftFunc
                (complex.complexArray
                    (vec.concat [shift, reals])
                    (vec.concat [shift, imags]));
-    
+
             map do c:
                 if complex.magnitude c <= thresh then complex.zero else c fi
                 done specKernel;
-    
+
         done [1..winNr];
     
     done [1..binsPerOctave];
     
     kmat = cm.toSparse
        (cm.scaled (1/fftSize)
            (cm.newComplexMatrix (RowMajor()) (concat kernels)));
     
-    eprintln "density = \(cm.density kmat)";
+    eprintln "density = \(cm.density kmat) (\(cm.nonZeroValues kmat) of \(cm.width kmat * cm.height kmat))";
     
     // Normalisation
     
     wx1 = bf.maxindex (complex.magnitudes (cm.getRow 0 kmat));
     wx2 = bf.maxindex (complex.magnitudes (cm.getRow (cm.height kmat - 1) kmat));

     wK = vec.slice diag (round(1/q)) (vec.length diag - round(1/q) - 2);
     
     weight = (fftHop / fftSize) / (bf.mean (bf.abs wK));
     weight = sqrt(weight);
 
+    eprintln "weight = \(weight)";
+    
     {
         kernel = cm.scaled weight kmat,
         fftSize,
         fftHop,
         binsPerOctave,