CLAM-Development  1.4.0
ConstantQTransform.cxx
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2001-2006 MUSIC TECHNOLOGY GROUP (MTG)
3  * UNIVERSITAT POMPEU FABRA
4  *
5  *
6  * This program is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * This program is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with this program; if not, write to the Free Software
18  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19  *
20  */
21 
22 #include <iostream>
23 #include "ConstantQTransform.hxx"
24 #include "FourierTransform.hxx"
25 #include <cmath>
26 
27 namespace Simac
28 {
29 //---------------------------------------------------------------------------
30 // nextpow2 returns the smallest integer n such that 2^n >= x.
31 static double nextpow2(double x)
32 {
33  return ceil(std::log(x)/std::log(2.0));
34 }
35 //----------------------------------------------------------------------------
36 static double squaredModule(const double & xx, const double & yy)
37 {
38  return xx*xx + yy*yy;
39 }
40 //----------------------------------------------------------------------------
41 
42 ConstantQTransform::ConstantQTransform(unsigned iFS, double ifmin, double ifmax, unsigned binsPerOctave)
43  : _binsPerOctave(binsPerOctave)
44 {
45  FS = iFS;
46  fmin = ifmin; // min freq
47  fmax = ifmax; // max freq
48 
49  Q = 1/(std::pow(2.,(1/(double)_binsPerOctave))-1); // Work out Q value for filter bank
50  K = (int) ceil(_binsPerOctave * std::log(fmax/fmin)/std::log(2.0)); // No. of constant Q bins
51 
52  // work out length of fft required for this constant Q filter bank
53  mSpectrumSize = (int) std::pow(2., nextpow2(ceil(Q*FS/fmin)));
54 
55  // allocate memory for cqdata
56  cqdata.resize(2*K);
57 }
58 
60 {
61 }
62 
64 {
65  //generates spectral kernel matrix (upside down?)
66  // initialise temporal kernel with zeros, twice length to deal w. complex numbers
67  double* hammingWindow = new double [2*mSpectrumSize];
68  for (unsigned mm=0; mm<2*mSpectrumSize; mm++) {
69  hammingWindow[mm] = 0;
70  }
71  // Here, fftleng*2 is a guess of the number of sparse cells in the matrix
72  // The matrix K x mSpectrumSize but the non-zero cells are an antialiased
73  // square root function. So mostly is a line, with some grey point.
74  mSparseKernelIs.reserve(mSpectrumSize*2);
75  mSparseKernelJs.reserve(mSpectrumSize*2);
76  mSparseKernelRealValues.reserve(mSpectrumSize*2);
77  mSparseKernelImagValues.reserve(mSpectrumSize*2);
78 
79  // for each bin value K, calculate temporal kernel, take its fft to
80  //calculate the spectral kernel then threshold it to make it sparse and
81  //add it to the sparse kernels matrix
82  double squareThreshold = thresh * thresh;
83  FourierTransform fft(mSpectrumSize, 1, true);
84  for (unsigned k=K; k--; ) {
85  // Computing a hamming window
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;
94  }
95 
96  //do fft of hammingWindow
97  fft.doIt(hammingWindow);
98  const double * transformedHamming = &fft.spectrum()[0];
99  // 2 steps because they are complex
100  for (unsigned j=0; j<(2*mSpectrumSize); j+=2) {
101  // perform thresholding
102  const double squaredBin = squaredModule(transformedHamming[j], transformedHamming[j+1]);
103  if (squaredBin <= squareThreshold) continue;
104  // Insert non-zero position indexes, doubled because they are floats
105  mSparseKernelIs.push_back(j);
106  mSparseKernelJs.push_back(k*2);
107  // take conjugate, normalise and add to array sparkernel
108  mSparseKernelRealValues.push_back( transformedHamming[j ]/mSpectrumSize);
109  mSparseKernelImagValues.push_back(-transformedHamming[j+1]/mSpectrumSize);
110  }
111  }
112  delete [] hammingWindow;
113 }
114 
115 //-----------------------------------------------------------------------------
116 void ConstantQTransform::doIt(const std::vector<double> & fftdata)
117 {
118  // Multiply by sparse kernels matrix and store in cqdata
119  // N.B. complex data
120  // rows = mSpectrumSize
121  // columns = K
122  double * constantQSpectrum = &cqdata[0];
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;
127  }
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++)
134  {
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];
141  // add the multiplication
142  constantQSpectrum[row ] += r1*r2 - i1*i2;
143  constantQSpectrum[row+1] += r1*i2 + i1*r2;
144  }
145 }
146 
147 }
148 
149