diff options
Diffstat (limited to 'lib/qm-dsp/dsp/chromagram/ConstantQ.cpp')
-rw-r--r-- | lib/qm-dsp/dsp/chromagram/ConstantQ.cpp | 387 |
1 files changed, 131 insertions, 256 deletions
diff --git a/lib/qm-dsp/dsp/chromagram/ConstantQ.cpp b/lib/qm-dsp/dsp/chromagram/ConstantQ.cpp index a5b1c12105..f2fad31f69 100644 --- a/lib/qm-dsp/dsp/chromagram/ConstantQ.cpp +++ b/lib/qm-dsp/dsp/chromagram/ConstantQ.cpp @@ -14,57 +14,16 @@ #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); -} - -static double squaredModule(const double & xx, const double & yy) { - return xx*xx + yy*yy; -} - //---------------------------------------------------------------------------- -ConstantQ::ConstantQ( CQConfig Config ) : +ConstantQ::ConstantQ( CQConfig config ) : m_sparseKernel(0) { - initialise( Config ); + initialise(config); } ConstantQ::~ConstantQ() @@ -72,179 +31,130 @@ ConstantQ::~ConstantQ() deInitialise(); } -//---------------------------------------------------------------------------- +static double squaredModule(const double & xx, const double & yy) { + return xx*xx + yy*yy; +} + 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 + 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; - //generates spectral kernel matrix (upside down?) - // initialise temporal kernel with zeros, twice length to deal w. complex numbers + 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; + } - double* hammingWindowRe = new double [ m_FFTLength ]; - double* hammingWindowIm = new double [ m_FFTLength ]; - double* transfHammingWindowRe = new double [ m_FFTLength ]; - double* transfHammingWindowIm = new double [ m_FFTLength ]; + // Compute a complex sinusoid windowed with a hamming window + // of the right length - for (unsigned u=0; u < m_FFTLength; u++) - { - hammingWindowRe[u] = 0; - hammingWindowIm[u] = 0; - } + const double samplesPerCycle = + m_FS / (m_FMin * pow(2, (double)j / (double)m_BPO)); + int windowLength = (int)ceil(m_dQ * samplesPerCycle); - // 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 squareThreshold = m_CQThresh * m_CQThresh; + int origin = m_FFTLength/2 - windowLength/2; - 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; + for (int i = 0; i < windowLength; ++i) { + double angle = (2.0 * M_PI * i) / samplesPerCycle; + windowRe[origin + i] = cos(angle); + windowIm[origin + i] = sin(angle); } - - const double samplesPerCycle = - m_FS / (m_FMin * pow(2, (double)k / (double)m_BPO)); - // Computing a hamming window - const unsigned hammingLength = (int) ceil( - m_dQ * samplesPerCycle); + // Shape with hamming window + Window<double> hamming(HammingWindow, windowLength); + hamming.cut(windowRe + origin); + hamming.cut(windowIm + origin); - unsigned origin = m_FFTLength/2 - hammingLength/2; - - for (unsigned i=0; i<hammingLength; i++) { - const double angle = 2*PI*i/samplesPerCycle; - 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; + // Scale + for (int i = 0; i < windowLength; ++i) { + windowRe[origin + i] /= windowLength; + } + for (int i = 0; i < windowLength; ++i) { + windowIm[origin + i] /= windowLength; } - 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; + // 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 ); + fft.process(false, windowRe, windowIm, transfWindowRe, transfWindowIm); + // convert to sparse form + for (int i = 0; i < m_FFTLength; i++) { - for (unsigned j=0; j<( m_FFTLength ); j++) { // perform thresholding - const double squaredBin = squaredModule( transfHammingWindowRe[ j ], transfHammingWindowIm[ j ]); - if (squaredBin <= squareThreshold) { - continue; - } + double mag = squaredModule(transfWindowRe[i], transfWindowIm[i]); + if (mag <= squareThreshold) continue; + // Insert non-zero position indexes - sk->is.push_back(j); - sk->js.push_back(k); + sk->is.push_back(i); + sk->js.push_back(j); - // take conjugate, normalise and add to array sparkernel - sk->real.push_back( transfHammingWindowRe[ j ]/m_FFTLength); - sk->imag.push_back(-transfHammingWindowIm[ j ]/m_FFTLength); + // 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 [] hammingWindowRe; - delete [] hammingWindowIm; - delete [] transfHammingWindowRe; - delete [] transfHammingWindowIm; + delete [] windowRe; + delete [] windowIm; + delete [] transfWindowRe; + delete [] transfWindowIm; -/* - 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; +} + +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; } //----------------------------------------------------------------------------- @@ -257,26 +167,24 @@ double* ConstantQ::process( const double* fftdata ) SparseKernel *sk = m_sparseKernel; - for (unsigned row=0; row<2*m_uK; row++) { + for (int 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]; - 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 ]; + 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); @@ -285,37 +193,6 @@ double* ConstantQ::process( const double* fftdata ) 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) { @@ -326,27 +203,25 @@ void ConstantQ::process(const double *FFTRe, const double* FFTIm, SparseKernel *sk = m_sparseKernel; - for (unsigned row=0; row<m_uK; row++) { + for (int 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]; - 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 ]; + 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); |