comparison dsp/chromagram/ConstantQ.cpp @ 505:930b5b0f707d

Merge branch 'codestyle-and-tidy'

     COPYING included with this distribution for more information.
 */
 
 #include "ConstantQ.h"
 #include "dsp/transforms/FFT.h"
+#include "base/Window.h"
 
 #include <iostream>
 
-#ifdef NOT_DEFINED
-// see note in CQprecalc
-
-#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;
-}
-#endif
-
-//---------------------------------------------------------------------------
-// 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(y);
+//----------------------------------------------------------------------------
+
+ConstantQ::ConstantQ( CQConfig config ) :
+    m_sparseKernel(0)
+{
+    initialise(config);
+}
+
+ConstantQ::~ConstantQ()
+{
+    deInitialise();
 }
 
 static double squaredModule(const double & xx, const double & yy) {
     return xx*xx + yy*yy;
 }
 
-//----------------------------------------------------------------------------
-
-ConstantQ::ConstantQ( CQConfig Config ) :
-    m_sparseKernel(0)
-{
-    initialise( Config );
-}
-
-ConstantQ::~ConstantQ()
-{
-    deInitialise();
-}
-
-//----------------------------------------------------------------------------
 void ConstantQ::sparsekernel()
 {
-//    std::cerr << "ConstantQ: initialising sparse kernel, uK = " << m_uK << ", FFTLength = " << m_FFTLength << "...";
-
     SparseKernel *sk = new SparseKernel();
 
-#ifdef NOT_DEFINED
-    if (push_precalculated(m_uK, m_FFTLength,
-                           sk->is, sk->js, sk->real, sk->imag)) {
-//        std::cerr << "using precalculated kernel" << std::endl;
-        m_sparseKernel = sk;
-        return;
-    }
-#endif
-
-    //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 ];
-    double* transfHammingWindowRe = new double [ m_FFTLength ];
-    double* transfHammingWindowIm = new double [ m_FFTLength ];
-
-    for (unsigned u=0; u < m_FFTLength; u++) 
-    {
-	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.
-    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* windowRe = new double [ m_FFTLength ];
+    double* windowIm = new double [ m_FFTLength ];
+    double* transfWindowRe = new double [ m_FFTLength ];
+    double* transfWindowIm = new double [ m_FFTLength ];
+        
+    // 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;
 
-    FFT m_FFT(m_FFTLength);
-	
-    for (unsigned k = m_uK; k--; ) 
-    {
-        for (unsigned u=0; u < m_FFTLength; u++) 
-        {
-            hammingWindowRe[u] = 0;
-            hammingWindowIm[u] = 0;
-        }
-        
-	// Computing a hamming window
-	const unsigned hammingLength = (int) ceil( m_dQ * m_FS / ( m_FMin * pow(2,((double)(k))/(double)m_BPO)));
-
-//        cerr << "k = " << k << ", q = " << m_dQ << ", m_FMin = " << m_FMin << ", hammingLength = " << hammingLength << " (rounded up from " << (m_dQ * m_FS / ( m_FMin * pow(2,((double)(k))/(double)m_BPO))) << ")" << endl;
-        
-
-        unsigned origin = m_FFTLength/2 - hammingLength/2;
-
-	for (unsigned i=0; i<hammingLength; i++) 
-	{
-	    const double angle = 2*PI*m_dQ*i/hammingLength;
-	    const double real = cos(angle);
-	    const double imag = sin(angle);
-	    const double absol = hamming(hammingLength, i)/hammingLength;
-	    hammingWindowRe[ origin + i ] = absol*real;
-	    hammingWindowIm[ origin + i ] = absol*imag;
-	}
-
-        for (unsigned i = 0; i < m_FFTLength/2; ++i) {
-            double temp = hammingWindowRe[i];
-            hammingWindowRe[i] = hammingWindowRe[i + m_FFTLength/2];
-            hammingWindowRe[i + m_FFTLength/2] = temp;
-            temp = hammingWindowIm[i];
-            hammingWindowIm[i] = hammingWindowIm[i + m_FFTLength/2];
-            hammingWindowIm[i + m_FFTLength/2] = temp;
+    FFT fft(m_FFTLength);
+        
+    for (int j = m_uK - 1; j >= 0; --j) {
+        
+        for (int i = 0; i < m_FFTLength; ++i) {
+            windowRe[i] = 0;
+            windowIm[i] = 0;
+        }
+
+        // Compute a complex sinusoid windowed with a hamming window
+        // of the right length
+        
+        int windowLength = (int)ceil
+            (m_dQ * m_FS / (m_FMin * pow(2, (double)j / (double)m_BPO)));
+
+        int origin = m_FFTLength/2 - windowLength/2;
+
+        for (int i = 0; i < windowLength; ++i) {
+            double angle = (2.0 * M_PI * m_dQ * i) / windowLength;
+            windowRe[origin + i] = cos(angle);
+            windowIm[origin + i] = sin(angle);
+        }
+
+        // Shape with hamming window
+        Window<double> hamming(HammingWindow, windowLength);
+        hamming.cut(windowRe + origin);
+        hamming.cut(windowIm + origin);
+
+        // Scale
+        for (int i = 0; i < windowLength; ++i) {
+            windowRe[origin + i] /= windowLength;
+        }
+        for (int i = 0; i < windowLength; ++i) {
+            windowIm[origin + i] /= windowLength;
+        }
+
+        // Input is expected to have been fftshifted, so do the
+        // same to the input to the fft that contains the kernel
+        for (int i = 0; i < m_FFTLength/2; ++i) {
+            double temp = windowRe[i];
+            windowRe[i] = windowRe[i + m_FFTLength/2];
+            windowRe[i + m_FFTLength/2] = temp;
+        }
+        for (int i = 0; i < m_FFTLength/2; ++i) {
+            double temp = windowIm[i];
+            windowIm[i] = windowIm[i + m_FFTLength/2];
+            windowIm[i + m_FFTLength/2] = temp;
         }
     
-	//do fft of hammingWindow
-	m_FFT.process( 0, hammingWindowRe, hammingWindowIm, transfHammingWindowRe, transfHammingWindowIm );
-
-		
-	for (unsigned j=0; j<( m_FFTLength ); j++) 
-	{
-	    // perform thresholding
-	    const double squaredBin = squaredModule( transfHammingWindowRe[ j ], transfHammingWindowIm[ j ]);
-	    if (squaredBin <= squareThreshold) continue;
-		
-	    // Insert non-zero position indexes
-	    sk->is.push_back(j);
-	    sk->js.push_back(k);
-
-	    // take conjugate, normalise and add to array sparkernel
-	    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;
-    
+        fft.process(false, windowRe, windowIm, transfWindowRe, transfWindowIm);
+
+        // convert to sparse form
+        for (int i = 0; i < m_FFTLength; i++) {
+
+            // perform thresholding
+            double mag = squaredModule(transfWindowRe[i], transfWindowIm[i]);
+            if (mag <= squareThreshold) continue;
+                
+            // Insert non-zero position indexes
+            sk->is.push_back(i);
+            sk->js.push_back(j);
+
+            // take conjugate, normalise and add to array for sparse kernel
+            sk->real.push_back( transfWindowRe[i] / m_FFTLength);
+            sk->imag.push_back(-transfWindowIm[i] / m_FFTLength);
+        }
+    }
+
+    delete [] windowRe;
+    delete [] windowIm;
+    delete [] transfWindowRe;
+    delete [] transfWindowIm;
+
     m_sparseKernel = sk;
-    return;
+}
+
+void ConstantQ::initialise( CQConfig Config )
+{
+    m_FS = Config.FS;             // Sample rate
+    m_FMin = Config.min;          // Minimum frequency
+    m_FMax = Config.max;          // Maximum frequency
+    m_BPO = Config.BPO;           // Bins per octave
+    m_CQThresh = Config.CQThresh; // Threshold for sparse kernel generation
+
+    // Q value for filter bank
+    m_dQ = 1/(pow(2,(1/(double)m_BPO))-1);
+
+    // No. of constant Q bins
+    m_uK = (int)ceil(m_BPO * log(m_FMax/m_FMin)/log(2.0));
+
+    // Length of fft required for this Constant Q filter bank
+    m_FFTLength = MathUtilities::nextPowerOfTwo(int(ceil(m_dQ * m_FS/m_FMin)));
+
+    // Hop from one frame to next
+    m_hop = m_FFTLength / 8;
+
+    // allocate memory for cqdata
+    m_CQdata = new double [2*m_uK];
+}
+
+void ConstantQ::deInitialise()
+{
+    delete [] m_CQdata;
+    delete m_sparseKernel;
 }
 
 //-----------------------------------------------------------------------------
 double* ConstantQ::process( const double* fftdata )
 {

         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 = &(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];
+    for (int row=0; row < 2 * m_uK; row++) {
+        m_CQdata[ row ] = 0;
+        m_CQdata[ row+1 ] = 0;
+    }
+    const int *fftbin = &(sk->is[0]);
+    const int *cqbin = &(sk->js[0]);
+    const double *real = &(sk->real[0]);
+    const double *imag = &(sk->imag[0]);
+    const int sparseCells = int(sk->real.size());
+        
+    for (int i = 0; i < sparseCells; i++) {
+        const int row = cqbin[i];
+        const int col = fftbin[i];
         if (col == 0) continue;
-	const double & r1  = real[i];
-	const double & i1  = imag[i];
-	const double & r2  = fftdata[ (2*m_FFTLength) - 2*col - 2 ];
-	const double & i2  = fftdata[ (2*m_FFTLength) - 2*col - 2 + 1 ];
-	// add the multiplication
-	m_CQdata[ 2*row  ] += (r1*r2 - i1*i2);
-	m_CQdata[ 2*row+1] += (r1*i2 + i1*r2);
+        const double & r1 = real[i];
+        const double & i1 = imag[i];
+        const double & r2 = fftdata[ (2*m_FFTLength) - 2*col - 2 ];
+        const double & i2 = fftdata[ (2*m_FFTLength) - 2*col - 2 + 1 ];
+        // add the multiplication
+        m_CQdata[ 2*row  ] += (r1*r2 - i1*i2);
+        m_CQdata[ 2*row+1] += (r1*i2 + i1*r2);
     }
 
     return m_CQdata;
-}
-
-
-void ConstantQ::initialise( CQConfig Config )
-{
-    m_FS = Config.FS;
-    m_FMin = Config.min;		// min freq
-    m_FMax = Config.max;		// max freq
-    m_BPO = Config.BPO;		// bins per octave
-    m_CQThresh = Config.CQThresh;// ConstantQ threshold for kernel generation
-
-    m_dQ = 1/(pow(2,(1/(double)m_BPO))-1);	// Work out Q value for Filter bank
-    m_uK = (unsigned int) ceil(m_BPO * log(m_FMax/m_FMin)/log(2.0));	// No. of constant Q bins
-
-//    std::cerr << "ConstantQ::initialise: rate = " << m_FS << ", fmin = " << m_FMin << ", fmax = " << m_FMax << ", bpo = " << m_BPO << ", K = " << m_uK << ", Q = " << m_dQ << std::endl;
-
-    // work out length of fft required for this constant Q Filter bank
-    m_FFTLength = (int) pow(2, nextpow2(ceil( m_dQ*m_FS/m_FMin )));
-
-    m_hop = m_FFTLength/8;
-
-//    std::cerr << "ConstantQ::initialise: -> fft length = " << m_FFTLength << ", hop = " << m_hop << std::endl;
-
-    // allocate memory for cqdata
-    m_CQdata = new double [2*m_uK];
-}
-
-void ConstantQ::deInitialise()
-{
-    delete [] m_CQdata;
-    delete m_sparseKernel;
 }
 
 void ConstantQ::process(const double *FFTRe, const double* FFTIm,
                         double *CQRe, double *CQIm)
 {

         return;
     }
 
     SparseKernel *sk = m_sparseKernel;
 
-    for (unsigned row=0; row<m_uK; row++) 
-    {
-	CQRe[ row ] = 0;
-	CQIm[ row ] = 0;
-    }
-
-    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];
+    for (int row = 0; row < m_uK; row++) {
+        CQRe[ row ] = 0;
+        CQIm[ row ] = 0;
+    }
+
+    const int *fftbin = &(sk->is[0]);
+    const int *cqbin = &(sk->js[0]);
+    const double *real = &(sk->real[0]);
+    const double *imag = &(sk->imag[0]);
+    const int sparseCells = int(sk->real.size());
+        
+    for (int i = 0; i<sparseCells; i++) {
+        const int row = cqbin[i];
+        const int col = fftbin[i];
         if (col == 0) continue;
-	const double & r1  = real[i];
-	const double & i1  = imag[i];
-	const double & r2  = FFTRe[ m_FFTLength - col ];
-	const double & i2  = FFTIm[ m_FFTLength - col ];
-	// add the multiplication
-	CQRe[ row ] += (r1*r2 - i1*i2);
-	CQIm[ row ] += (r1*i2 + i1*r2);
-    }
-}
+        const double & r1 = real[i];
+        const double & i1 = imag[i];
+        const double & r2 = FFTRe[ m_FFTLength - col ];
+        const double & i2 = FFTIm[ m_FFTLength - col ];
+        // add the multiplication
+        CQRe[ row ] += (r1*r2 - i1*i2);
+        CQIm[ row ] += (r1*i2 + i1*r2);
+    }
+}