changeset 93:908be1d06bd2

More work on inverse

--- a/cpp-qm-dsp/CQInverse.cpp	Fri May 09 11:15:56 2014 +0100
+++ b/cpp-qm-dsp/CQInverse.cpp	Fri May 09 11:41:44 2014 +0100
@@ -173,9 +173,9 @@
     // and stack these so as to achieve a list, for each octave, of
     // taller columns of height binsPerOctave * atomsPerFrame
     //
-    // 3. For each column, take the product with the inverse CQ kernel
-    // (which is the conjugate of the forward kernel) and perform an
-    // inverse FFT
+    // 3. For each taller column, take the product with the inverse CQ
+    // kernel (which is the conjugate of the forward kernel) and
+    // perform an inverse FFT
     //
     // 4. Overlap-add each octave's resynthesised blocks (unwindowed)
     //
@@ -225,19 +225,43 @@
             ("Columns in octave must be a multiple of atoms per frame");
     }
 
-    ComplexBlock reshaped;
     for (int i = 0; i < ncols; i += m_p.atomsPerFrame) {
+
         ComplexColumn tallcol;
         for (int b = 0; b < m_binsPerOctave; ++b) {
             for (int a = 0; a < m_p.atomsPerFrame; ++a) {
                 tallcol.push_back(columns[i + a][b]);
             }
         }
-        reshaped.push_back(tallcol);
+        
+        processOctaveColumn(octave, tallcol);
+    }
+}
+
+void
+CQInverse::processOctaveColumn(int octave, const ComplexColumn &column)
+{
+    // 3. For each taller column, take the product with the inverse CQ
+    // kernel (which is the conjugate of the forward kernel) and
+    // perform an inverse FFT
+
+    ComplexSequence transformed = m_kernel->processInverse(column);
+
+    int halfLen = m_p.fftSize/2 + 1;
+
+    RealSequence ri(halfLen, 0);
+    RealSequence ii(halfLen, 0);
+
+    for (int i = 0; i < halfLen; ++i) {
+        ri[i] = transformed[i].real();
+        ii[i] = transformed[i].imag();
     }
 
-    
+    RealSequence timeDomain(m_p.fftSize, 0);
 
+    m_fft->inverse(ri.data(), ii.data(), timeDomain.data());
+
+
+    //...
 }
 
-

--- a/cpp-qm-dsp/CQInverse.h	Fri May 09 11:15:56 2014 +0100
+++ b/cpp-qm-dsp/CQInverse.h	Fri May 09 11:41:44 2014 +0100
@@ -83,6 +83,7 @@
     
     void initialise();
     void processOctave(int octave, const ComplexBlock &block);
+    void processOctaveColumn(int octave, const ComplexColumn &column);
 };
 
 #endif

--- a/cpp-qm-dsp/CQKernel.cpp	Fri May 09 11:15:56 2014 +0100
+++ b/cpp-qm-dsp/CQKernel.cpp	Fri May 09 11:41:44 2014 +0100
@@ -241,9 +241,10 @@
     
     cerr << "weight = " << weight << endl;
 
-    // apply normalisation weight, make sparse, and store conjugates
-    // (our multiplication order means we will effectively be using
-    // the adjoint or conjugate transpose of the kernel matrix)
+    // apply normalisation weight, make sparse, and store conjugate
+    // (we use the adjoint or conjugate transpose of the kernel matrix
+    // for the forward transform, the plain kernel for the inverse
+    // which we expect to be less common)
 
     KernelMatrix sk;
 
@@ -274,14 +275,14 @@
 }
 
 vector<C>
-CQKernel::process(const vector<C> &cv)
+CQKernel::processForward(const vector<C> &cv)
 {
-    // matrix multiply m_kernel.data by in, converting in to complex
-    // as we go
+    // straightforward matrix multiply (taking into account m_kernel's
+    // slightly-sparse representation)
 
     int nrows = m_p.binsPerOctave * m_p.atomsPerFrame;
 
-    vector<C> rv(nrows, C(0, 0));
+    vector<C> rv(nrows, C());
 
     for (int i = 0; i < nrows; ++i) {
         for (int j = 0; j < (int)m_kernel.data[i].size(); ++j) {
@@ -292,4 +293,28 @@
     return rv;
 }
 
+vector<C>
+CQKernel::processInverse(const vector<C> &cv)
+{
+    // matrix multiply by conjugate transpose of m_kernel. This is
+    // actually the original kernel as calculated, we just stored the
+    // conjugate-transpose of the kernel because we expect to be doing
+    // more forward transforms than inverse ones.
 
+    int ncols = m_p.binsPerOctave * m_p.atomsPerFrame;
+    int nrows = m_p.fftSize;
+
+    vector<C> rv(nrows, C());
+
+    for (int j = 0; j < ncols; ++j) {
+        int i0 = m_kernel.origin[j];
+        int i1 = i0 + m_kernel.data[j].size();
+        for (int i = i0; i < i1; ++i) {
+            rv[i] += cv[j] * conj(m_kernel.data[j][i - i0]);
+        }
+    }
+
+    return rv;
+}
+
+

--- a/cpp-qm-dsp/CQKernel.h	Fri May 09 11:15:56 2014 +0100
+++ b/cpp-qm-dsp/CQKernel.h	Fri May 09 11:41:44 2014 +0100
@@ -59,7 +59,10 @@
 
     Properties getProperties() const { return m_p; }
 
-    std::vector<std::complex<double> > process
+    std::vector<std::complex<double> > processForward
+        (const std::vector<std::complex<double> > &);
+
+    std::vector<std::complex<double> > processInverse
         (const std::vector<std::complex<double> > &);
 
 private:

--- a/cpp-qm-dsp/ConstantQ.cpp	Fri May 09 11:15:56 2014 +0100
+++ b/cpp-qm-dsp/ConstantQ.cpp	Fri May 09 11:41:44 2014 +0100
@@ -305,7 +305,7 @@
         cv.push_back(Complex(ro[i], io[i]));
     }
 
-    ComplexSequence cqrowvec = m_kernel->process(cv);
+    ComplexSequence cqrowvec = m_kernel->processForward(cv);
 
     // Reform into a column matrix
     ComplexBlock cqblock;