summaryrefslogtreecommitdiffstats
path: root/lib/qm-dsp/dsp/chromagram/ConstantQ.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'lib/qm-dsp/dsp/chromagram/ConstantQ.cpp')
-rw-r--r--lib/qm-dsp/dsp/chromagram/ConstantQ.cpp387
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);
generated by cgit v1.2.3 (git 2.25.1) at 2024-09-10 15:37:32 +0000