qm-dsp
1.8
|
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 }