33 return ceil(std::log(x)/std::log(2.0));
43 : _binsPerOctave(binsPerOctave)
49 Q = 1/(std::pow(2.,(1/(
double)_binsPerOctave))-1);
50 K = (
int) ceil(_binsPerOctave * std::log(fmax/fmin)/std::log(2.0));
53 mSpectrumSize = (
int) std::pow(2.,
nextpow2(ceil(Q*FS/fmin)));
67 double* hammingWindow =
new double [2*mSpectrumSize];
68 for (
unsigned mm=0; mm<2*mSpectrumSize; mm++) {
69 hammingWindow[mm] = 0;
74 mSparseKernelIs.reserve(mSpectrumSize*2);
75 mSparseKernelJs.reserve(mSpectrumSize*2);
76 mSparseKernelRealValues.reserve(mSpectrumSize*2);
77 mSparseKernelImagValues.reserve(mSpectrumSize*2);
82 double squareThreshold = thresh * thresh;
84 for (
unsigned k=K; k--; ) {
86 const unsigned hammingLength = (
int) ceil(Q * FS / (fmin * std::pow(2.,((
double)(k))/(
double)_binsPerOctave)));
87 for (
unsigned i=0; i<hammingLength; i++) {
88 const double angle = 2*M_PI*Q*i/hammingLength;
89 const double real = cos(angle);
90 const double imag = sin(angle);
91 const double absol = Hamming(hammingLength, i)/hammingLength;
92 hammingWindow[2*i ] = absol*real;
93 hammingWindow[2*i+1] = absol*imag;
97 fft.
doIt(hammingWindow);
98 const double * transformedHamming = &fft.
spectrum()[0];
100 for (
unsigned j=0; j<(2*mSpectrumSize); j+=2) {
102 const double squaredBin =
squaredModule(transformedHamming[j], transformedHamming[j+1]);
103 if (squaredBin <= squareThreshold)
continue;
105 mSparseKernelIs.push_back(j);
106 mSparseKernelJs.push_back(k*2);
108 mSparseKernelRealValues.push_back( transformedHamming[j ]/mSpectrumSize);
109 mSparseKernelImagValues.push_back(-transformedHamming[j+1]/mSpectrumSize);
112 delete [] hammingWindow;
123 const double * fftSpectrum = &fftdata[0];
124 for (
unsigned row=0; row<2*K; row+=2) {
125 constantQSpectrum[row ] = 0;
126 constantQSpectrum[row+1] = 0;
128 const unsigned *fftbin = &(mSparseKernelIs[0]);
129 const unsigned *cqbin = &(mSparseKernelJs[0]);
130 const double *real = &(mSparseKernelRealValues[0]);
131 const double *imag = &(mSparseKernelImagValues[0]);
132 const unsigned int nSparseCells = mSparseKernelRealValues.size();
133 for (
unsigned i = 0; i<nSparseCells; i++)
135 const unsigned row = cqbin[i];
136 const unsigned col = fftbin[i];
137 const double & r1 = real[i];
138 const double & i1 = imag[i];
139 const double & r2 = fftSpectrum[col];
140 const double & i2 = fftSpectrum[col+1];
142 constantQSpectrum[row ] += r1*r2 - i1*i2;
143 constantQSpectrum[row+1] += r1*i2 + i1*r2;