comparison dsp/chromagram/ConstantQ.cpp @ 276:4c901426b9f3

* Do not calculate CQ sparse kernel when chromagram is constructed: only
  when it's actually used
* Pre-calculate CQ sparse kernels in the sizes required for the default
  configurations of some of our transforms

 
 #include "ConstantQ.h"
 #include "dsp/transforms/FFT.h"
 
 #include <iostream>
+
+#include "CQprecalc.cpp"
+
+static bool push_precalculated(int uk, int fftlength,
+                               std::vector<unsigned> &is,
+                               std::vector<unsigned> &js,
+                               std::vector<double> &real,
+                               std::vector<double> &imag)
+{
+    if (uk == 76 && fftlength == 16384) {
+        push_76_16384(is, js, real, imag);
+        return true;
+    }
+    if (uk == 144 && fftlength == 4096) {
+        push_144_4096(is, js, real, imag);
+        return true;
+    }
+    if (uk == 65 && fftlength == 2048) {
+        push_65_2048(is, js, real, imag);
+        return true;
+    }
+    if (uk == 84 && fftlength == 65536) {
+        push_84_65536(is, js, real, imag);
+        return true;
+    }
+    return false;
+}
 
 //---------------------------------------------------------------------------
 // nextpow2 returns the smallest integer n such that 2^n >= x.
 static double nextpow2(double x) {
     double y = ceil(log(x)/log(2.0));

     return xx*xx + yy*yy;
 }
 
 //----------------------------------------------------------------------------
 
-ConstantQ::ConstantQ( CQConfig Config ) 
+ConstantQ::ConstantQ( CQConfig Config ) :
+    m_sparseKernel(0)
 {
     initialise( Config );
 }
 
 ConstantQ::~ConstantQ()

 }
 
 //----------------------------------------------------------------------------
 void ConstantQ::sparsekernel()
 {
+//    std::cerr << "ConstantQ: initialising sparse kernel, uK = " << m_uK << ", FFTLength = " << m_FFTLength << "...";
+
+    SparseKernel *sk = new SparseKernel();
+
+    if (push_precalculated(m_uK, m_FFTLength,
+                           sk->is, sk->js, sk->real, sk->imag)) {
+        m_sparseKernel = sk;
+        return;
+    }
+    
     //generates spectral kernel matrix (upside down?)
     // initialise temporal kernel with zeros, twice length to deal w. complex numbers
 
     double* hammingWindowRe = new double [ m_FFTLength ];
     double* hammingWindowIm = new double [ m_FFTLength ];

     {
 	hammingWindowRe[u] = 0;
 	hammingWindowIm[u] = 0;
     }
 
-
     // Here, fftleng*2 is a guess of the number of sparse cells in the matrix
     // The matrix K x fftlength but the non-zero cells are an antialiased
     // square root function. So mostly is a line, with some grey point.
-    m_sparseKernelIs.reserve( m_FFTLength*2 );
-    m_sparseKernelJs.reserve( m_FFTLength*2 );
-    m_sparseKernelRealValues.reserve( m_FFTLength*2 );
-    m_sparseKernelImagValues.reserve( m_FFTLength*2 );
+    sk->is.reserve( m_FFTLength*2 );
+    sk->js.reserve( m_FFTLength*2 );
+    sk->real.reserve( m_FFTLength*2 );
+    sk->imag.reserve( m_FFTLength*2 );
 	
     // for each bin value K, calculate temporal kernel, take its fft to
     //calculate the spectral kernel then threshold it to make it sparse and 
     //add it to the sparse kernels matrix
     double squareThreshold = m_CQThresh * m_CQThresh;

 	    // perform thresholding
 	    const double squaredBin = squaredModule( transfHammingWindowRe[ j ], transfHammingWindowIm[ j ]);
 	    if (squaredBin <= squareThreshold) continue;
 		
 	    // Insert non-zero position indexes, doubled because they are floats
-	    m_sparseKernelIs.push_back(j);
-	    m_sparseKernelJs.push_back(k);
+	    sk->is.push_back(j);
+	    sk->js.push_back(k);
 
 	    // take conjugate, normalise and add to array sparkernel
-	    m_sparseKernelRealValues.push_back( transfHammingWindowRe[ j ]/m_FFTLength);
-	    m_sparseKernelImagValues.push_back(-transfHammingWindowIm[ j ]/m_FFTLength);
+	    sk->real.push_back( transfHammingWindowRe[ j ]/m_FFTLength);
+	    sk->imag.push_back(-transfHammingWindowIm[ j ]/m_FFTLength);
 	}
 
     }
 
     delete [] hammingWindowRe;
     delete [] hammingWindowIm;
     delete [] transfHammingWindowRe;
     delete [] transfHammingWindowIm;
 
+/*
+    using std::cout;
+    using std::endl;
+
+    cout.precision(28);
+
+    int n = sk->is.size();
+    int w = 8;
+    cout << "static unsigned int sk_i_" << m_uK << "_" << m_FFTLength << "[" << n << "] = {" << endl;
+    for (int i = 0; i < n; ++i) {
+        if (i % w == 0) cout << "    ";
+        cout << sk->is[i];
+        if (i + 1 < n) cout << ", ";
+        if (i % w == w-1) cout << endl;
+    };
+    if (n % w != 0) cout << endl;
+    cout << "};" << endl;
+
+    n = sk->js.size();
+    cout << "static unsigned int sk_j_" << m_uK << "_" << m_FFTLength << "[" << n << "] = {" << endl;
+    for (int i = 0; i < n; ++i) {
+        if (i % w == 0) cout << "    ";
+        cout << sk->js[i];
+        if (i + 1 < n) cout << ", ";
+        if (i % w == w-1) cout << endl;
+    };
+    if (n % w != 0) cout << endl;
+    cout << "};" << endl;
+
+    w = 2;
+    n = sk->real.size();
+    cout << "static double sk_real_" << m_uK << "_" << m_FFTLength << "[" << n << "] = {" << endl;
+    for (int i = 0; i < n; ++i) {
+        if (i % w == 0) cout << "    ";
+        cout << sk->real[i];
+        if (i + 1 < n) cout << ", ";
+        if (i % w == w-1) cout << endl;
+    };
+    if (n % w != 0) cout << endl;
+    cout << "};" << endl;
+
+    n = sk->imag.size();
+    cout << "static double sk_imag_" << m_uK << "_" << m_FFTLength << "[" << n << "] = {" << endl;
+    for (int i = 0; i < n; ++i) {
+        if (i % w == 0) cout << "    ";
+        cout << sk->imag[i];
+        if (i + 1 < n) cout << ", ";
+        if (i % w == w-1) cout << endl;
+    };
+    if (n % w != 0) cout << endl;
+    cout << "};" << endl;
+
+    cout << "static void push_" << m_uK << "_" << m_FFTLength << "(vector<unsigned int> &is, vector<unsigned int> &js, vector<double> &real, vector<double> &imag)" << endl;
+    cout << "{\n    is.reserve(" << n << ");\n";
+    cout << "    js.reserve(" << n << ");\n";
+    cout << "    real.reserve(" << n << ");\n";
+    cout << "    imag.reserve(" << n << ");\n";
+    cout << "    for (int i = 0; i < " << n << "; ++i) {" << endl;
+    cout << "        is.push_back(sk_i_" << m_uK << "_" << m_FFTLength << "[i]);" << endl;
+    cout << "        js.push_back(sk_j_" << m_uK << "_" << m_FFTLength << "[i]);" << endl;
+    cout << "        real.push_back(sk_real_" << m_uK << "_" << m_FFTLength << "[i]);" << endl;
+    cout << "        imag.push_back(sk_imag_" << m_uK << "_" << m_FFTLength << "[i]);" << endl;
+    cout << "    }" << endl;
+    cout << "}" << endl;
+*/
+//    std::cerr << "done\n -> is: " << sk->is.size() << ", js: " << sk->js.size() << ", reals: " << sk->real.size() << ", imags: " << sk->imag.size() << std::endl;
+    
+    m_sparseKernel = sk;
+    return;
 }
 
 //-----------------------------------------------------------------------------
 double* ConstantQ::process( const double* fftdata )
 {
+    if (!m_sparseKernel) {
+        std::cerr << "ERROR: ConstantQ::process: Sparse kernel has not been initialised" << std::endl;
+        return m_CQdata;
+    }
+
+    SparseKernel *sk = m_sparseKernel;
+
     for (unsigned row=0; row<2*m_uK; row++) 
     {
 	m_CQdata[ row ] = 0;
 	m_CQdata[ row+1 ] = 0;
     }
-    const unsigned *fftbin = &(m_sparseKernelIs[0]);
-    const unsigned *cqbin  = &(m_sparseKernelJs[0]);
-    const double   *real   = &(m_sparseKernelRealValues[0]);
-    const double   *imag   = &(m_sparseKernelImagValues[0]);
-    const unsigned int sparseCells = m_sparseKernelRealValues.size();
+    const unsigned *fftbin = &(sk->is[0]);
+    const unsigned *cqbin  = &(sk->js[0]);
+    const double   *real   = &(sk->real[0]);
+    const double   *imag   = &(sk->imag[0]);
+    const unsigned int sparseCells = sk->real.size();
 	
     for (unsigned i = 0; i<sparseCells; i++)
     {
 	const unsigned row = cqbin[i];
 	const unsigned col = fftbin[i];

 }
 
 void ConstantQ::deInitialise()
 {
     delete [] m_CQdata;
+    delete m_sparseKernel;
 }
 
 void ConstantQ::process(const double *FFTRe, const double* FFTIm,
                         double *CQRe, double *CQIm)
 {
+    if (!m_sparseKernel) {
+        std::cerr << "ERROR: ConstantQ::process: Sparse kernel has not been initialised" << std::endl;
+        return;
+    }
+
+    SparseKernel *sk = m_sparseKernel;
+
     for (unsigned row=0; row<m_uK; row++) 
     {
 	CQRe[ row ] = 0;
 	CQIm[ row ] = 0;
     }
 
-    const unsigned *fftbin = &(m_sparseKernelIs[0]);
-    const unsigned *cqbin  = &(m_sparseKernelJs[0]);
-    const double   *real   = &(m_sparseKernelRealValues[0]);
-    const double   *imag   = &(m_sparseKernelImagValues[0]);
-    const unsigned int sparseCells = m_sparseKernelRealValues.size();
+    const unsigned *fftbin = &(sk->is[0]);
+    const unsigned *cqbin  = &(sk->js[0]);
+    const double   *real   = &(sk->real[0]);
+    const double   *imag   = &(sk->imag[0]);
+    const unsigned int sparseCells = sk->real.size();
 	
     for (unsigned i = 0; i<sparseCells; i++)
     {
 	const unsigned row = cqbin[i];
 	const unsigned col = fftbin[i];