Mercurial > hg > constant-q-cpp

changeset 23:4f2f28d5df58

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fix frequency calculation in cpp, and some diagnostics
author Chris Cannam <c.cannam@qmul.ac.uk>
date Fri, 01 Nov 2013 17:58:39 +0000
parents 701900c371b0
children 3973be829352
files cpp-qm-dsp/CQKernel.cpp cpp-qm-dsp/test.cpp yeti/cqtkernel.yeti
diffstat 3 files changed, 67 insertions(+), 19 deletions(-) [+]
line wrap: on
line diff
--- a/cpp-qm-dsp/CQKernel.cpp	Fri Nov 01 16:13:22 2013 +0000
+++ b/cpp-qm-dsp/CQKernel.cpp	Fri Nov 01 17:58:39 2013 +0000
@@ -8,8 +8,11 @@
 #include <cmath>
 #include <cassert>
 #include <vector>
+#include <iostream>
 
 using std::vector;
+using std::cerr;
+using std::endl;
 
 CQKernel::CQKernel(double sampleRate, double maxFreq, int binsPerOctave)
 {
@@ -70,18 +73,18 @@
 	for (int i = 0; i < (int)win.size(); ++i) {
 	    win[i] = sqrt(win[i]) / nk;
 	}
-	
-	double fk = m_p.minFrequency * pow(2, ((k-1.0) * bpo));
+
+	double fk = m_p.minFrequency * pow(2, ((k-1.0) / bpo));
 
 	vector<double> reals, imags;
 	
 	for (int i = 0; i < nk; ++i) {
-	    double arg = 2.0 * M_PI * fk * i / m_p.sampleRate;
-	    reals[i] = win[i] * cos(arg);
-	    imags[i] = win[i] * sin(arg);
+	    double arg = (2.0 * M_PI * fk * i) / m_p.sampleRate;
+	    reals.push_back(win[i] * cos(arg));
+	    imags.push_back(win[i] * sin(arg));
 	}
 
-	int atomOffset = m_p.firstCentre - ceil(nk/2);
+	int atomOffset = m_p.firstCentre - int(ceil(nk/2.0));
 
 	for (int i = 0; i < m_p.atomsPerFrame; ++i) {
 
@@ -92,12 +95,12 @@
 
 	    for (int j = 0; j < nk; ++j) {
 		rin[j + shift] = reals[j];
-		iin[i + shift] = imags[j];
+		iin[j + shift] = imags[j];
 	    }
 
 	    vector<double> rout(m_p.fftSize, 0.0);
 	    vector<double> iout(m_p.fftSize, 0.0);
-		
+
 	    m_fft->process(false,
 			   rin.data(), iin.data(),
 			   rout.data(), iout.data());
@@ -105,9 +108,11 @@
 	    int firstNZ = -1, lastNZ = -1;
 
 	    for (int j = 0; j < m_p.fftSize; ++j) {
-		if (sqrt(rout[i] * rout[i] + iout[i] * iout[i]) >= thresh) {
+		if (sqrt(rout[j] * rout[j] + iout[j] * iout[j]) >= thresh) {
 		    lastNZ = j;
 		    if (firstNZ < 0) firstNZ = j;
+		} else {
+		    rout[j] = iout[j] = 0;
 		}
 	    }
 
@@ -128,11 +133,52 @@
 	}
     }
 
-    assert((int)m_kernel.offsets.size() == m_p.binsPerOctave);
-    assert((int)m_kernel.real.size() == m_p.binsPerOctave);
-    assert((int)m_kernel.imag.size() == m_p.binsPerOctave);
+    assert((int)m_kernel.offsets.size() == m_p.binsPerOctave * m_p.atomsPerFrame);
+    assert((int)m_kernel.real.size() == m_p.binsPerOctave * m_p.atomsPerFrame);
+    assert((int)m_kernel.imag.size() == m_p.binsPerOctave * m_p.atomsPerFrame);
+
+    // print density as diagnostic
+
+    int nnz = 0;
+    for (int i = 0; i < m_kernel.offsets.size(); ++i) {
+	assert(m_kernel.real[i].size() == m_kernel.imag[i].size());
+	for (int j = 0; j < m_kernel.real[i].size(); ++j) {
+	    if (m_kernel.real[i][j] != 0.0 ||
+		m_kernel.imag[i][j] != 0.0) {
+		++nnz;
+	    }
+	}
+    }
+
+    cerr << "density = " << double(nnz) / double(m_p.binsPerOctave * m_p.atomsPerFrame * m_p.fftSize) << " (" << nnz << " of " << m_p.binsPerOctave * m_p.atomsPerFrame * m_p.fftSize << ")" << endl;
 
     //!!! and normalise
+
+/* far too laborious
+
+    int wx1 = -1, wx2 = -1;
+    double wm1 = 0, wm2 = 0;
+
+    for (int i = 0; i < m_kernel.real[0].size(); ++i) {
+	m = m_kernel.real[0][i] * m_kernel.real[0][i] +
+	    m_kernel.imag[0][i] * m_kernel.imag[0][i];
+	if (wx1 == -1 || m > wm1) {
+	    wx1 = i + m_kernel.offsets[0];
+	    wm1 = m;
+	}
+    }
+    
+    int n = m_kernel.offsets.size() - 1;
+    for (int i = 0; i < m_kernel.real[n].size(); ++i) {
+	m = m_kernel.real[n][i] * m_kernel.real[n][i] +
+	    m_kernel.imag[n][i] * m_kernel.imag[n][i];
+	if (wx2 == -1 || m > wm1) {
+	    wx2 = i + m_kernel.offsets[n];
+	    wm2 = m;
+	}
+    }
+
+*/    
 }
 
 
--- a/cpp-qm-dsp/test.cpp	Fri Nov 01 16:13:22 2013 +0000
+++ b/cpp-qm-dsp/test.cpp	Fri Nov 01 17:58:39 2013 +0000
@@ -5,7 +5,7 @@
 
 int main(int argc, char **argv)
 {
-    CQKernel k(96000, 48000, 24);
+    CQKernel k(48000, 24000, 24);
 
     std::cerr << "Q = " << k.getProperties().Q << std::endl;
 
--- a/yeti/cqtkernel.yeti	Fri Nov 01 16:13:22 2013 +0000
+++ b/yeti/cqtkernel.yeti	Fri Nov 01 17:58:39 2013 +0000
@@ -47,16 +47,16 @@
     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]));
@@ -74,11 +74,11 @@
                (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];
@@ -87,7 +87,7 @@
        (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
     
@@ -102,6 +102,8 @@
     weight = (fftHop / fftSize) / (bf.mean (bf.abs wK));
     weight = sqrt(weight);
 
+    eprintln "weight = \(weight)";
+    
     {
         kernel = cm.scaled weight kmat,
         fftSize,