qm-dsp  1.8
ConstantQ.cpp
Go to the documentation of this file.
00001 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
00002 /*
00003     QM DSP Library
00004 
00005     Centre for Digital Music, Queen Mary, University of London.
00006     This file 2005-2006 Christian Landone.
00007 
00008     This program is free software; you can redistribute it and/or
00009     modify it under the terms of the GNU General Public License as
00010     published by the Free Software Foundation; either version 2 of the
00011     License, or (at your option) any later version.  See the file
00012     COPYING included with this distribution for more information.
00013 */
00014 
00015 #include "ConstantQ.h"
00016 #include "dsp/transforms/FFT.h"
00017 
00018 #include <iostream>
00019 
00020 #ifdef NOT_DEFINED
00021 // see note in CQprecalc
00022 
00023 #include "CQprecalc.cpp"
00024 
00025 static bool push_precalculated(int uk, int fftlength,
00026                                std::vector<unsigned> &is,
00027                                std::vector<unsigned> &js,
00028                                std::vector<double> &real,
00029                                std::vector<double> &imag)
00030 {
00031     if (uk == 76 && fftlength == 16384) {
00032         push_76_16384(is, js, real, imag);
00033         return true;
00034     }
00035     if (uk == 144 && fftlength == 4096) {
00036         push_144_4096(is, js, real, imag);
00037         return true;
00038     }
00039     if (uk == 65 && fftlength == 2048) {
00040         push_65_2048(is, js, real, imag);
00041         return true;
00042     }
00043     if (uk == 84 && fftlength == 65536) {
00044         push_84_65536(is, js, real, imag);
00045         return true;
00046     }
00047     return false;
00048 }
00049 #endif
00050 
00051 //---------------------------------------------------------------------------
00052 // nextpow2 returns the smallest integer n such that 2^n >= x.
00053 static double nextpow2(double x) {
00054     double y = ceil(log(x)/log(2.0));
00055     return(y);
00056 }
00057 
00058 static double squaredModule(const double & xx, const double & yy) {
00059     return xx*xx + yy*yy;
00060 }
00061 
00062 //----------------------------------------------------------------------------
00063 
00064 ConstantQ::ConstantQ( CQConfig Config ) :
00065     m_sparseKernel(0)
00066 {
00067     initialise( Config );
00068 }
00069 
00070 ConstantQ::~ConstantQ()
00071 {
00072     deInitialise();
00073 }
00074 
00075 //----------------------------------------------------------------------------
00076 void ConstantQ::sparsekernel()
00077 {
00078 //    std::cerr << "ConstantQ: initialising sparse kernel, uK = " << m_uK << ", FFTLength = " << m_FFTLength << "...";
00079 
00080     SparseKernel *sk = new SparseKernel();
00081 
00082 #ifdef NOT_DEFINED
00083     if (push_precalculated(m_uK, m_FFTLength,
00084                            sk->is, sk->js, sk->real, sk->imag)) {
00085 //        std::cerr << "using precalculated kernel" << std::endl;
00086         m_sparseKernel = sk;
00087         return;
00088     }
00089 #endif
00090 
00091     //generates spectral kernel matrix (upside down?)
00092     // initialise temporal kernel with zeros, twice length to deal w. complex numbers
00093 
00094     double* hammingWindowRe = new double [ m_FFTLength ];
00095     double* hammingWindowIm = new double [ m_FFTLength ];
00096     double* transfHammingWindowRe = new double [ m_FFTLength ];
00097     double* transfHammingWindowIm = new double [ m_FFTLength ];
00098 
00099     for (unsigned u=0; u < m_FFTLength; u++) 
00100     {
00101         hammingWindowRe[u] = 0;
00102         hammingWindowIm[u] = 0;
00103     }
00104 
00105     // Here, fftleng*2 is a guess of the number of sparse cells in the matrix
00106     // The matrix K x fftlength but the non-zero cells are an antialiased
00107     // square root function. So mostly is a line, with some grey point.
00108     sk->is.reserve( m_FFTLength*2 );
00109     sk->js.reserve( m_FFTLength*2 );
00110     sk->real.reserve( m_FFTLength*2 );
00111     sk->imag.reserve( m_FFTLength*2 );
00112         
00113     // for each bin value K, calculate temporal kernel, take its fft to
00114     //calculate the spectral kernel then threshold it to make it sparse and 
00115     //add it to the sparse kernels matrix
00116     double squareThreshold = m_CQThresh * m_CQThresh;
00117 
00118     FFT m_FFT(m_FFTLength);
00119         
00120     for (unsigned k = m_uK; k--; ) 
00121     {
00122         for (unsigned u=0; u < m_FFTLength; u++) 
00123         {
00124             hammingWindowRe[u] = 0;
00125             hammingWindowIm[u] = 0;
00126         }
00127         
00128         // Computing a hamming window
00129         const unsigned hammingLength = (int) ceil( m_dQ * m_FS / ( m_FMin * pow(2,((double)(k))/(double)m_BPO)));
00130 
00131         unsigned origin = m_FFTLength/2 - hammingLength/2;
00132 
00133         for (unsigned i=0; i<hammingLength; i++) 
00134         {
00135             const double angle = 2*PI*m_dQ*i/hammingLength;
00136             const double real = cos(angle);
00137             const double imag = sin(angle);
00138             const double absol = hamming(hammingLength, i)/hammingLength;
00139             hammingWindowRe[ origin + i ] = absol*real;
00140             hammingWindowIm[ origin + i ] = absol*imag;
00141         }
00142 
00143         for (unsigned i = 0; i < m_FFTLength/2; ++i) {
00144             double temp = hammingWindowRe[i];
00145             hammingWindowRe[i] = hammingWindowRe[i + m_FFTLength/2];
00146             hammingWindowRe[i + m_FFTLength/2] = temp;
00147             temp = hammingWindowIm[i];
00148             hammingWindowIm[i] = hammingWindowIm[i + m_FFTLength/2];
00149             hammingWindowIm[i + m_FFTLength/2] = temp;
00150         }
00151     
00152         //do fft of hammingWindow
00153         m_FFT.process( 0, hammingWindowRe, hammingWindowIm, transfHammingWindowRe, transfHammingWindowIm );
00154 
00155                 
00156         for (unsigned j=0; j<( m_FFTLength ); j++) 
00157         {
00158             // perform thresholding
00159             const double squaredBin = squaredModule( transfHammingWindowRe[ j ], transfHammingWindowIm[ j ]);
00160             if (squaredBin <= squareThreshold) continue;
00161                 
00162             // Insert non-zero position indexes, doubled because they are floats
00163             sk->is.push_back(j);
00164             sk->js.push_back(k);
00165 
00166             // take conjugate, normalise and add to array sparkernel
00167             sk->real.push_back( transfHammingWindowRe[ j ]/m_FFTLength);
00168             sk->imag.push_back(-transfHammingWindowIm[ j ]/m_FFTLength);
00169         }
00170 
00171     }
00172 
00173     delete [] hammingWindowRe;
00174     delete [] hammingWindowIm;
00175     delete [] transfHammingWindowRe;
00176     delete [] transfHammingWindowIm;
00177 
00178 /*
00179     using std::cout;
00180     using std::endl;
00181 
00182     cout.precision(28);
00183 
00184     int n = sk->is.size();
00185     int w = 8;
00186     cout << "static unsigned int sk_i_" << m_uK << "_" << m_FFTLength << "[" << n << "] = {" << endl;
00187     for (int i = 0; i < n; ++i) {
00188         if (i % w == 0) cout << "    ";
00189         cout << sk->is[i];
00190         if (i + 1 < n) cout << ", ";
00191         if (i % w == w-1) cout << endl;
00192     };
00193     if (n % w != 0) cout << endl;
00194     cout << "};" << endl;
00195 
00196     n = sk->js.size();
00197     cout << "static unsigned int sk_j_" << m_uK << "_" << m_FFTLength << "[" << n << "] = {" << endl;
00198     for (int i = 0; i < n; ++i) {
00199         if (i % w == 0) cout << "    ";
00200         cout << sk->js[i];
00201         if (i + 1 < n) cout << ", ";
00202         if (i % w == w-1) cout << endl;
00203     };
00204     if (n % w != 0) cout << endl;
00205     cout << "};" << endl;
00206 
00207     w = 2;
00208     n = sk->real.size();
00209     cout << "static double sk_real_" << m_uK << "_" << m_FFTLength << "[" << n << "] = {" << endl;
00210     for (int i = 0; i < n; ++i) {
00211         if (i % w == 0) cout << "    ";
00212         cout << sk->real[i];
00213         if (i + 1 < n) cout << ", ";
00214         if (i % w == w-1) cout << endl;
00215     };
00216     if (n % w != 0) cout << endl;
00217     cout << "};" << endl;
00218 
00219     n = sk->imag.size();
00220     cout << "static double sk_imag_" << m_uK << "_" << m_FFTLength << "[" << n << "] = {" << endl;
00221     for (int i = 0; i < n; ++i) {
00222         if (i % w == 0) cout << "    ";
00223         cout << sk->imag[i];
00224         if (i + 1 < n) cout << ", ";
00225         if (i % w == w-1) cout << endl;
00226     };
00227     if (n % w != 0) cout << endl;
00228     cout << "};" << endl;
00229 
00230     cout << "static void push_" << m_uK << "_" << m_FFTLength << "(vector<unsigned int> &is, vector<unsigned int> &js, vector<double> &real, vector<double> &imag)" << endl;
00231     cout << "{\n    is.reserve(" << n << ");\n";
00232     cout << "    js.reserve(" << n << ");\n";
00233     cout << "    real.reserve(" << n << ");\n";
00234     cout << "    imag.reserve(" << n << ");\n";
00235     cout << "    for (int i = 0; i < " << n << "; ++i) {" << endl;
00236     cout << "        is.push_back(sk_i_" << m_uK << "_" << m_FFTLength << "[i]);" << endl;
00237     cout << "        js.push_back(sk_j_" << m_uK << "_" << m_FFTLength << "[i]);" << endl;
00238     cout << "        real.push_back(sk_real_" << m_uK << "_" << m_FFTLength << "[i]);" << endl;
00239     cout << "        imag.push_back(sk_imag_" << m_uK << "_" << m_FFTLength << "[i]);" << endl;
00240     cout << "    }" << endl;
00241     cout << "}" << endl;
00242 */
00243 //    std::cerr << "done\n -> is: " << sk->is.size() << ", js: " << sk->js.size() << ", reals: " << sk->real.size() << ", imags: " << sk->imag.size() << std::endl;
00244     
00245     m_sparseKernel = sk;
00246     return;
00247 }
00248 
00249 //-----------------------------------------------------------------------------
00250 double* ConstantQ::process( const double* fftdata )
00251 {
00252     if (!m_sparseKernel) {
00253         std::cerr << "ERROR: ConstantQ::process: Sparse kernel has not been initialised" << std::endl;
00254         return m_CQdata;
00255     }
00256 
00257     SparseKernel *sk = m_sparseKernel;
00258 
00259     for (unsigned row=0; row<2*m_uK; row++) 
00260     {
00261         m_CQdata[ row ] = 0;
00262         m_CQdata[ row+1 ] = 0;
00263     }
00264     const unsigned *fftbin = &(sk->is[0]);
00265     const unsigned *cqbin  = &(sk->js[0]);
00266     const double   *real   = &(sk->real[0]);
00267     const double   *imag   = &(sk->imag[0]);
00268     const unsigned int sparseCells = sk->real.size();
00269         
00270     for (unsigned i = 0; i<sparseCells; i++)
00271     {
00272         const unsigned row = cqbin[i];
00273         const unsigned col = fftbin[i];
00274         const double & r1  = real[i];
00275         const double & i1  = imag[i];
00276         const double & r2  = fftdata[ (2*m_FFTLength) - 2*col - 2 ];
00277         const double & i2  = fftdata[ (2*m_FFTLength) - 2*col - 2 + 1 ];
00278         // add the multiplication
00279         m_CQdata[ 2*row  ] += (r1*r2 - i1*i2);
00280         m_CQdata[ 2*row+1] += (r1*i2 + i1*r2);
00281     }
00282 
00283     return m_CQdata;
00284 }
00285 
00286 
00287 void ConstantQ::initialise( CQConfig Config )
00288 {
00289     m_FS = Config.FS;
00290     m_FMin = Config.min;                // min freq
00291     m_FMax = Config.max;                // max freq
00292     m_BPO = Config.BPO;         // bins per octave
00293     m_CQThresh = Config.CQThresh;// ConstantQ threshold for kernel generation
00294 
00295     m_dQ = 1/(pow(2,(1/(double)m_BPO))-1);      // Work out Q value for Filter bank
00296     m_uK = (unsigned int) ceil(m_BPO * log(m_FMax/m_FMin)/log(2.0));    // No. of constant Q bins
00297 
00298 //    std::cerr << "ConstantQ::initialise: rate = " << m_FS << ", fmin = " << m_FMin << ", fmax = " << m_FMax << ", bpo = " << m_BPO << ", K = " << m_uK << ", Q = " << m_dQ << std::endl;
00299 
00300     // work out length of fft required for this constant Q Filter bank
00301     m_FFTLength = (int) pow(2, nextpow2(ceil( m_dQ*m_FS/m_FMin )));
00302 
00303     m_hop = m_FFTLength/8; // <------ hop size is window length divided by 32
00304 
00305 //    std::cerr << "ConstantQ::initialise: -> fft length = " << m_FFTLength << ", hop = " << m_hop << std::endl;
00306 
00307     // allocate memory for cqdata
00308     m_CQdata = new double [2*m_uK];
00309 }
00310 
00311 void ConstantQ::deInitialise()
00312 {
00313     delete [] m_CQdata;
00314     delete m_sparseKernel;
00315 }
00316 
00317 void ConstantQ::process(const double *FFTRe, const double* FFTIm,
00318                         double *CQRe, double *CQIm)
00319 {
00320     if (!m_sparseKernel) {
00321         std::cerr << "ERROR: ConstantQ::process: Sparse kernel has not been initialised" << std::endl;
00322         return;
00323     }
00324 
00325     SparseKernel *sk = m_sparseKernel;
00326 
00327     for (unsigned row=0; row<m_uK; row++) 
00328     {
00329         CQRe[ row ] = 0;
00330         CQIm[ row ] = 0;
00331     }
00332 
00333     const unsigned *fftbin = &(sk->is[0]);
00334     const unsigned *cqbin  = &(sk->js[0]);
00335     const double   *real   = &(sk->real[0]);
00336     const double   *imag   = &(sk->imag[0]);
00337     const unsigned int sparseCells = sk->real.size();
00338         
00339     for (unsigned i = 0; i<sparseCells; i++)
00340     {
00341         const unsigned row = cqbin[i];
00342         const unsigned col = fftbin[i];
00343         const double & r1  = real[i];
00344         const double & i1  = imag[i];
00345         const double & r2  = FFTRe[ m_FFTLength - col - 1 ];
00346         const double & i2  = FFTIm[ m_FFTLength - col - 1 ];
00347         // add the multiplication
00348         CQRe[ row ] += (r1*r2 - i1*i2);
00349         CQIm[ row ] += (r1*i2 + i1*r2);
00350     }
00351 }