Mercurial > hg > qm-dsp
changeset 495:1bea13b8f951
Style fixes in constant-Q: avoid unsigned, reuse our Window class, fix
comments
| author | Chris Cannam <cannam@all-day-breakfast.com> |
|---|---|
| date | Fri, 31 May 2019 18:25:31 +0100 |
| parents | 3f649fbb1172 |
| children | 97008bc9420f |
| files | dsp/chromagram/Chromagram.cpp dsp/chromagram/ConstantQ.cpp dsp/chromagram/ConstantQ.h dsp/segmentation/ClusterMeltSegmenter.cpp |
| diffstat | 4 files changed, 150 insertions(+), 169 deletions(-) [+] |
line wrap: on
line diff
--- a/dsp/chromagram/Chromagram.cpp Fri May 31 18:15:59 2019 +0100 +++ b/dsp/chromagram/Chromagram.cpp Fri May 31 18:25:31 2019 +0100 @@ -57,8 +57,8 @@ m_uK = m_ConstantQ->getK(); // Initialise working arrays - m_frameSize = m_ConstantQ->getfftlength(); - m_hopSize = m_ConstantQ->gethop(); + m_frameSize = m_ConstantQ->getFFTLength(); + m_hopSize = m_ConstantQ->getHop(); // Initialise FFT object m_FFT = new FFTReal(m_frameSize);
--- a/dsp/chromagram/ConstantQ.cpp Fri May 31 18:15:59 2019 +0100 +++ b/dsp/chromagram/ConstantQ.cpp Fri May 31 18:25:31 2019 +0100 @@ -14,26 +14,16 @@ #include "ConstantQ.h" #include "dsp/transforms/FFT.h" +#include "base/Window.h" #include <iostream> -//--------------------------------------------------------------------------- -// 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() @@ -41,100 +31,129 @@ 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(); - //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 ); + 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 + // 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); + FFT 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 j = m_uK; j >= 0; --j) { - // 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*M_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 (int i = 0; i < m_FFTLength; ++i) { + windowRe[i] = 0; + windowIm[i] = 0; } - 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; + // 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++) { + fft.process(false, windowRe, windowIm, transfWindowRe, transfWindowIm); + + // convert to sparse form + for (int i = 0; i < m_FFTLength; i++) { + // 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; -// 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; } //----------------------------------------------------------------------------- @@ -147,24 +166,24 @@ 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(); + 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 (unsigned i = 0; i<sparseCells; i++) { - const unsigned row = cqbin[i]; - const unsigned col = fftbin[i]; + 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 ]; + 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); @@ -173,37 +192,6 @@ 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) { @@ -214,25 +202,25 @@ 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(); + 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 (unsigned i = 0; i<sparseCells; i++) { - const unsigned row = cqbin[i]; - const unsigned col = fftbin[i]; + 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 ]; + 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);
--- a/dsp/chromagram/ConstantQ.h Fri May 31 18:15:59 2019 +0100 +++ b/dsp/chromagram/ConstantQ.h Fri May 31 18:25:31 2019 +0100 @@ -31,29 +31,23 @@ class ConstantQ { public: - void process( const double* FFTRe, const double* FFTIm, - double* CQRe, double* CQIm ); - - ConstantQ( CQConfig Config ); + ConstantQ(CQConfig config); ~ConstantQ(); - double* process( const double* FFTData ); + void process(const double* FFTRe, const double* FFTIm, + double* CQRe, double* CQIm); + + double* process(const double* FFTData); void sparsekernel(); - double hamming(int len, int n) { - double out = 0.54 - 0.46*cos(2*M_PI*n/len); - return(out); - } - - int getnumwin() { return m_numWin;} - double getQ() { return m_dQ;} - int getK() {return m_uK ;} - int getfftlength() { return m_FFTLength;} - int gethop() { return m_hop;} + double getQ() { return m_dQ; } + int getK() { return m_uK; } + int getFFTLength() { return m_FFTLength; } + int getHop() { return m_hop; } private: - void initialise( CQConfig Config ); + void initialise(CQConfig config); void deInitialise(); double* m_CQdata; @@ -62,15 +56,14 @@ double m_FMax; double m_dQ; double m_CQThresh; - unsigned int m_numWin; - unsigned int m_hop; - unsigned int m_BPO; - unsigned int m_FFTLength; - unsigned int m_uK; + int m_hop; + int m_BPO; + int m_FFTLength; + int m_uK; struct SparseKernel { - std::vector<unsigned> is; - std::vector<unsigned> js; + std::vector<int> is; + std::vector<int> js; std::vector<double> imag; std::vector<double> real; };
--- a/dsp/segmentation/ClusterMeltSegmenter.cpp Fri May 31 18:15:59 2019 +0100 +++ b/dsp/segmentation/ClusterMeltSegmenter.cpp Fri May 31 18:25:31 2019 +0100 @@ -82,7 +82,7 @@ ncoeff = constq->getK(); - fft = new FFTReal(constq->getfftlength()); + fft = new FFTReal(constq->getFFTLength()); } else if (featureType == FEATURE_TYPE_MFCC) { @@ -156,7 +156,7 @@ return; } - int fftsize = constq->getfftlength(); + int fftsize = constq->getFFTLength(); if (!window || window->getSize() != fftsize) { delete window;