Mercurial > hg > constant-q-cpp

changeset 26:a5f71b5c9f85

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Minor updates
author Chris Cannam <c.cannam@qmul.ac.uk>
date Mon, 04 Nov 2013 17:51:27 +0000
parents 5ca24ff67566
children 5435f00a4103
files cpp-qm-dsp/CQKernel.cpp cpp-qm-dsp/CQKernel.h
diffstat 2 files changed, 106 insertions(+), 89 deletions(-) [+]
line wrap: on
line diff
--- a/cpp-qm-dsp/CQKernel.cpp	Mon Nov 04 14:20:29 2013 +0000
+++ b/cpp-qm-dsp/CQKernel.cpp	Mon Nov 04 17:51:27 2013 +0000
@@ -1,3 +1,4 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
 
 #include "CQKernel.h"
 
@@ -9,12 +10,15 @@
 #include <cassert>
 #include <vector>
 #include <iostream>
+#include <algorithm>
 
 using std::vector;
 using std::complex;
 using std::cerr;
 using std::endl;
 
+typedef std::complex<double> C;
+
 CQKernel::CQKernel(double sampleRate, double maxFreq, int binsPerOctave)
 {
     m_p.sampleRate = sampleRate;
@@ -42,16 +46,16 @@
 
     double maxNK = round(m_p.Q * m_p.sampleRate / m_p.minFrequency);
     double minNK = round
-	(m_p.Q * m_p.sampleRate /
-	 (m_p.minFrequency * pow(2, (bpo - 1.0) / bpo)));
+        (m_p.Q * m_p.sampleRate /
+         (m_p.minFrequency * pow(2, (bpo - 1.0) / bpo)));
 
     m_p.atomSpacing = round(minNK * atomHopFactor);
     m_p.firstCentre = m_p.atomSpacing * ceil(ceil(maxNK / 2.0) / m_p.atomSpacing);
     m_p.fftSize = MathUtilities::nextPowerOfTwo
-	(m_p.firstCentre + ceil(maxNK / 2.0));
+        (m_p.firstCentre + ceil(maxNK / 2.0));
 
     m_p.atomsPerFrame = floor
-	(1.0 + (m_p.fftSize - ceil(maxNK / 2.0) - m_p.firstCentre) / m_p.atomSpacing);
+        (1.0 + (m_p.fftSize - ceil(maxNK / 2.0) - m_p.firstCentre) / m_p.atomSpacing);
 
     int lastCentre = m_p.firstCentre + (m_p.atomsPerFrame - 1) * m_p.atomSpacing;
 
@@ -60,110 +64,120 @@
     m_fft = new FFT(m_p.fftSize);
 
     for (int k = 1; k <= m_p.binsPerOctave; ++k) {
-	
-	int nk = round(m_p.Q * m_p.sampleRate /
-		       (m_p.minFrequency * pow(2, ((k-1.0) / bpo))));
+        
+        int nk = round(m_p.Q * m_p.sampleRate /
+                       (m_p.minFrequency * pow(2, ((k-1.0) / bpo))));
 
-	// The MATLAB version uses a symmetric window, but our windows
-	// are periodic. A symmetric window of size N is a periodic
-	// one of size N-1 with the first element stuck on the end
-	Window<double> w(BlackmanHarrisWindow, nk-1);
-	vector<double> win = w.getWindowData();
-	win.push_back(win[0]);
+        // The MATLAB version uses a symmetric window, but our windows
+        // are periodic. A symmetric window of size N is a periodic
+        // one of size N-1 with the first element stuck on the end
+        Window<double> w(BlackmanHarrisWindow, nk-1);
+        vector<double> win = w.getWindowData();
+        win.push_back(win[0]);
 
-	for (int i = 0; i < (int)win.size(); ++i) {
-	    win[i] = sqrt(win[i]) / nk;
-	}
+        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.push_back(win[i] * cos(arg));
-	    imags.push_back(win[i] * sin(arg));
-	}
+        vector<double> reals, imags;
+        
+        for (int i = 0; i < nk; ++i) {
+            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 - int(ceil(nk/2.0));
+        int atomOffset = m_p.firstCentre - int(ceil(nk/2.0));
 
-	for (int i = 0; i < m_p.atomsPerFrame; ++i) {
+        for (int i = 0; i < m_p.atomsPerFrame; ++i) {
 
-	    int shift = atomOffset + (i * m_p.atomSpacing);
+            int shift = atomOffset + (i * m_p.atomSpacing);
 
-	    vector<double> rin(m_p.fftSize, 0.0);
-	    vector<double> iin(m_p.fftSize, 0.0);
+            vector<double> rin(m_p.fftSize, 0.0);
+            vector<double> iin(m_p.fftSize, 0.0);
 
-	    for (int j = 0; j < nk; ++j) {
-		rin[j + shift] = reals[j];
-		iin[j + shift] = imags[j];
-	    }
+            for (int j = 0; j < nk; ++j) {
+                rin[j + shift] = reals[j];
+                iin[j + shift] = imags[j];
+            }
 
-	    vector<double> rout(m_p.fftSize, 0.0);
-	    vector<double> iout(m_p.fftSize, 0.0);
+            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());
+            m_fft->process(false,
+                           rin.data(), iin.data(),
+                           rout.data(), iout.data());
 
-	    int firstNZ = -1, lastNZ = -1;
+            // Keep this dense for the moment (until after
+            // normalisation calculations)
 
-	    vector<complex<double> > row;
+            vector<C> row;
 
-	    for (int j = 0; j < m_p.fftSize; ++j) {
-		if (sqrt(rout[j] * rout[j] + iout[j] * iout[j]) < thresh) {
-		    row.push_back(complex<double>(0, 0));
-		} else {
-		    row.push_back(complex<double>(rout[j], iout[j]));
-		}
-	    }
+            for (int j = 0; j < m_p.fftSize; ++j) {
+                if (sqrt(rout[j] * rout[j] + iout[j] * iout[j]) < thresh) {
+                    row.push_back(C(0, 0));
+                } else {
+                    row.push_back(C(rout[j], iout[j]));
+                }
+            }
 
-	    m_kernel.row.push_back(row);
-	}
+            m_kernel.origin.push_back(0);
+            m_kernel.data.push_back(row);
+        }
     }
 
-    assert((int)m_kernel.row.size() == m_p.binsPerOctave * m_p.atomsPerFrame);
+    assert((int)m_kernel.data.size() == m_p.binsPerOctave * m_p.atomsPerFrame);
 
     // print density as diagnostic
 
     int nnz = 0;
-    for (int i = 0; i < m_kernel.row.size(); ++i) {
-	for (int j = 0; j < m_kernel.row[i].size(); ++j) {
-	    if (m_kernel.row[i][j] != complex<double>(0, 0)) {
-		++nnz;
-	    }
-	}
+    for (int i = 0; i < m_kernel.data.size(); ++i) {
+        for (int j = 0; j < m_kernel.data[i].size(); ++j) {
+            if (m_kernel.data[i][j] != C(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
+    normaliseKernel();
+}
 
-/* far too laborious
+static bool ccomparator(C &c1, C &c2)
+{
+    return abs(c1) < abs(c2);
+}
 
-    int wx1 = -1, wx2 = -1;
-    double wm1 = 0, wm2 = 0;
+static int maxidx(vector<C> &v)
+{
+    return std::max_element(v.begin(), v.end(), ccomparator) - v.begin();
+}
 
-    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;
-	}
+void
+CQKernel::normaliseKernel()
+{
+    // and normalise
+
+    int wx1 = maxidx(m_kernel.data[0]);
+    int wx2 = maxidx(m_kernel.data[m_kernel.data.size()-1]);
+
+    vector<vector<C> > subset(m_kernel.data.size());
+    for (int j = wx1; j <= wx2; ++j) {
+        for (int i = 0; i < m_kernel.data.size(); ++i) {
+            subset[i].push_back(m_kernel.data[i][j]);
+        }
     }
 
-*/    
+    vector<vector<C> > square;
+    // conjugate transpose of subset * subset
+    for (int j = 0; j < subset[0].size(); ++j) {
+
+
+    }
+        
 }
 
 
--- a/cpp-qm-dsp/CQKernel.h	Mon Nov 04 14:20:29 2013 +0000
+++ b/cpp-qm-dsp/CQKernel.h	Mon Nov 04 17:51:27 2013 +0000
@@ -1,3 +1,4 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
 
 #ifndef CQ_KERNEL_H
 #define CQ_KERNEL_H
@@ -14,16 +15,16 @@
     ~CQKernel();
     
     struct Properties {
-	double sampleRate;
-	double maxFrequency;
-	double minFrequency;
-	int binsPerOctave;
-	int fftSize;
-	int fftHop;
-	int atomsPerFrame;
-	int atomSpacing;
-	int firstCentre;
-	double Q;
+        double sampleRate;
+        double maxFrequency;
+        double minFrequency;
+        int binsPerOctave;
+        int fftSize;
+        int fftHop;
+        int atomsPerFrame;
+        int atomSpacing;
+        int firstCentre;
+        double Q;
     };
 
     Properties getProperties() const { return m_p; }
@@ -33,11 +34,13 @@
     FFT *m_fft;
 
     struct KernelMatrix {
-	std::vector<std::vector<std::complex<double> > > row;
+        std::vector<int> origin;
+        std::vector<std::vector<std::complex<double> > > data;
     };
     KernelMatrix m_kernel;
 
     void generateKernel();
+    void normaliseKernel();
 };
 
 #endif