# HG changeset patch
# User Chris Cannam <c.cannam@qmul.ac.uk>
# Date 1383725433 0
# Node ID ba648b8672bde390b79db96463aa6c284ac5000e
# Parent  01a3e110bf8d2fa386390c1ede54821d23aa2d94
Much of the actual CQ processing logic

diff -r 01a3e110bf8d -r ba648b8672bd cpp-qm-dsp/CQKernel.cpp
--- a/cpp-qm-dsp/CQKernel.cpp	Tue Nov 05 16:58:18 2013 +0000
+++ b/cpp-qm-dsp/CQKernel.cpp	Wed Nov 06 08:10:33 2013 +0000
@@ -143,6 +143,11 @@
         }
     }
 
+    cerr << "size = " << m_kernel.data.size() << "*" << m_kernel.data[0].size() << " (fft size = " << m_p.fftSize << ")" << endl;
+
+    assert(m_kernel.data.size() == m_p.binsPerOctave * m_p.atomsPerFrame);
+    assert(m_kernel.data[0].size() == m_p.fftSize);
+
     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;
 
     finaliseKernel();
@@ -235,4 +240,24 @@
     m_kernel = sk;
 }
 
+vector<C>
+CQKernel::process(const vector<C> &cv)
+{
+    // matrix multiply m_kernel.data by in, converting in to complex
+    // as we go
 
+    int ncols = m_p.fftSize;
+    int nrows = m_p.binsPerOctave * m_p.atomsPerFrame;
+
+    vector<C> rv(nrows, C(0, 0));
+
+    for (int i = 0; i < nrows; ++i) {
+        for (int j = 0; j < m_kernel.data[i].size(); ++j) {
+            rv[i] += cv[j + m_kernel.origin[i]] * m_kernel.data[i][j];
+        }
+    }
+
+    return rv;
+}
+
+
diff -r 01a3e110bf8d -r ba648b8672bd cpp-qm-dsp/CQKernel.h
--- a/cpp-qm-dsp/CQKernel.h	Tue Nov 05 16:58:18 2013 +0000
+++ b/cpp-qm-dsp/CQKernel.h	Wed Nov 06 08:10:33 2013 +0000
@@ -29,6 +29,9 @@
 
     Properties getProperties() const { return m_p; }
 
+    std::vector<std::complex<double> > process
+        (const std::vector<std::complex<double> > &);
+
 private:
     Properties m_p;
     FFT *m_fft;
diff -r 01a3e110bf8d -r ba648b8672bd cpp-qm-dsp/ConstantQ.cpp
--- a/cpp-qm-dsp/ConstantQ.cpp	Tue Nov 05 16:58:18 2013 +0000
+++ b/cpp-qm-dsp/ConstantQ.cpp	Wed Nov 06 08:10:33 2013 +0000
@@ -17,6 +17,8 @@
 using std::cerr;
 using std::endl;
 
+typedef std::complex<double> C;
+
 ConstantQ::ConstantQ(double sampleRate,
 		     double minFreq,
 		     double maxFreq,
@@ -81,7 +83,7 @@
 	m_buffers.push_back(vector<double>(m_extraLatencies[i], 0.0));
     }
 
-    m_fft = new FFT(m_p.fftSize);
+    m_fft = new FFTReal(m_p.fftSize);
     m_bigBlockSize = m_p.fftSize * pow(2, m_octaves) / 2;
 
     cerr << "m_bigBlockSize = " << m_bigBlockSize << " for " << m_octaves << " octaves" << endl;
@@ -97,9 +99,72 @@
 	m_buffers[i].insert(m_buffers[i].end(), dec.begin(), dec.end());
     }
 
-    //!!! do the work!
+    vector<vector<double> > out;
 
-    vector<vector<double> > out;
+    //!!!! need some mechanism for handling remaining samples at the end
+
+    while (m_buffers[0].size() >= m_bigBlockSize) {
+	vector<vector<double> > chunk = processBigBlock();
+	for (int i = 0; i < m_octaves; ++i) {
+	    vector<double> v;
+	    v.insert(v.end(),
+		     m_buffers[i].begin() + m_bigBlockSize/pow(2,i), m_buffers[i].end());
+	    m_buffers[i] = v;
+	}
+	out.insert(out.end(), chunk.begin(), chunk.end());
+    }
+    
     return out;
 }
 
+vector<vector<double> > 
+ConstantQ::processBigBlock()
+{
+    vector<vector<double> > out;
+
+    for (int i = 0; i < pow(2, m_octaves - 1) * m_p.atomsPerFrame; ++i) {
+	out.push_back(vector<double>());
+    }
+    cerr << "returning " << out.size() << " cols per big block" << endl;
+
+    for (int i = 0; i < m_octaves; ++i) {
+	int n = pow(2, (m_octaves - i - 1));
+	cerr << "octave " << i+1 << " of " << m_octaves << ", n = " << n << endl;
+	for (int j = 0; j < n; ++j) {
+	    int origin = j * m_p.fftSize; //!!! no -- it's the hop but how does that add up?
+	    vector<vector<double> > block = 
+		processOctaveBlock(m_buffers[i].data() + origin);
+
+	    //...
+	}
+    }
+
+    return out;
+}
+	
+vector<vector<double> >
+ConstantQ::processOctaveBlock(const double *data)
+{
+    vector<double> ro(m_p.fftSize, 0.0);
+    vector<double> io(m_p.fftSize, 0.0);
+    m_fft->forward(data, ro.data(), io.data());
+
+    vector<C> cv;
+    for (int i = 0; i < m_p.fftSize; ++i) {
+	cv.push_back(C(ro[i], io[i]));
+    }
+
+    vector<C> cqrowvec = m_kernel->process(cv);
+
+    vector<vector<double> > cqblock;
+    for (int i = 0; i < m_p.binsPerOctave; ++i) {
+	cqblock.push_back(vector<double>());
+	for (int j = 0; j < m_p.atomsPerFrame; ++j) {
+	    cqblock[i].push_back(abs(cqrowvec[i * m_p.atomsPerFrame + j]));
+	}
+    }
+
+    return cqblock;
+}
+
+
diff -r 01a3e110bf8d -r ba648b8672bd cpp-qm-dsp/ConstantQ.h
--- a/cpp-qm-dsp/ConstantQ.h	Tue Nov 05 16:58:18 2013 +0000
+++ b/cpp-qm-dsp/ConstantQ.h	Wed Nov 06 08:10:33 2013 +0000
@@ -8,7 +8,7 @@
 #include <vector>
 
 class Resampler;
-class FFT;
+class FFTReal;
 
 class ConstantQ
 {
@@ -37,9 +37,11 @@
     int m_totalLatency;
     std::vector<int> m_extraLatencies; // per resampler, to make up to total
 
-    FFT *m_fft;
+    FFTReal *m_fft;
 
     void initialise();
+    std::vector<std::vector<double> > processBigBlock();
+    std::vector<std::vector<double> > processOctaveBlock(const double *data);
 };
 
 #endif
diff -r 01a3e110bf8d -r ba648b8672bd cpp-qm-dsp/test.cpp
--- a/cpp-qm-dsp/test.cpp	Tue Nov 05 16:58:18 2013 +0000
+++ b/cpp-qm-dsp/test.cpp	Wed Nov 06 08:10:33 2013 +0000
@@ -12,7 +12,7 @@
 {
     ConstantQ k(48000, 50, 24000, 24);
 
-    vector<double> in(1024, 0.0);
+    vector<double> in(65536*4, 0.0);
     vector<vector<double> > out = k.process(in);