CLAM-Development  1.4.0
FourierTransform.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 <cmath>
23 #include "FourierTransform.hxx"
24 #include <iostream>
25 
26 #define SWAP(a,b) do {double swaptemp=(a); (a)=(b); (b)=swaptemp;} while(false)
27 
28 // constructor
29 // input points to the input data. n is twice the frame size
30 // complex = 0 -> augmentation with 0s is necessary
31 FourierTransform::FourierTransform(unsigned long int framesize, double datanorm, bool isComplex)
32  : mFrameSize(framesize)
33  , mIsComplex(isComplex)
34 {
35  _spectrum.resize(2*mFrameSize);
36  #if USE_FFTW3
37  _realInput = (double*) fftw_malloc(sizeof(double) * mFrameSize);
38  _complexOutput = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * mFrameSize);
39 
40  fftw_import_system_wisdom();
41  if (isComplex)
42  _plan = fftw_plan_dft_1d(framesize, _complexOutput, _complexOutput, 1, FFTW_ESTIMATE);
43  else
44  _plan = fftw_plan_dft_r2c_1d(framesize, _realInput, _complexOutput, FFTW_ESTIMATE);
45  #endif
46 }
47 
49 {
50  #if USE_FFTW3
51  fftw_destroy_plan(_plan);
52  fftw_free(_realInput);
53  fftw_free(_complexOutput);
54  #endif
55 }
56 
57 //----------------------------------------------------------------------------------------
58 // FourierTransform transform code adapted from: Chapter 12, Fast Fourier Transform from Numerical
59 // Recipes in C: The Art of Scientific Computing. Numerical Recipes Software (1988-1992)
60 // Cambridge University Press, Cambridge, accessed 25.08.04,
61 // available at http://www.library.cornell.edu/nr/cbookcpdf.html
62 
63 void doFourierTransform(double * data, unsigned frameSize)
64 {
65  // Matlab indeces, data[0] is out of bounds TOREMOVE?
66  const unsigned long n=frameSize<<1; //n is the data array size
67  for (unsigned long i=1, j=1; i<n; i+=2) {
68  if (j>i) {
69  SWAP(data[j], data[i]);
70  SWAP(data[j+1], data[i+1]);
71  }
72  unsigned long m=frameSize;
73  while (m >= 2 && j>m) {
74  j-= m;
75  m >>= 1;
76  }
77  j += (m);
78  }
79  //now Danielson-Lanczos
80  unsigned long istep;
81  for (unsigned long mmax = 2; mmax<n; mmax = istep)
82  {
83  istep = mmax << 1;
84  const double theta =(6.28318530717959/mmax); // multiply by -1 here for ifft
85  double wtemp = sin(0.5*theta);
86  double wpr = -2.0*wtemp*wtemp;
87  double wpi = sin(theta);
88  double wr=1.0;
89  double wi=0.0;
90  for(unsigned long m=1; m<mmax; m+=2) {
91  for (unsigned long i=m; i<=n; i+=istep) {
92  unsigned long j = i+mmax;
93  double tempr = float (wr*data[j] - wi*data[j+1]);
94  double tempi = float (wr*data[j+1] + wi*data[j]);
95  data[j] = data[i] - tempr;
96  data[j+1] = data[i+1] - tempi;
97  data[i] += tempr;
98  data[i+1] += tempi;
99  }
100  wtemp=wr;
101  wr = wtemp*wpr - wi*wpi + wr; //trigonometric recurrence
102  //wr = (wtemp=wr)*wpr - wi*wpi + wr; //trigonometric recurrence
103  wi = wi*wpr + wtemp*wpi + wi;
104  }
105  }
106 }
107 
108 template <typename T>
109 void FourierTransform::doItGeneric(const T * input)
110 {
111  double * spectrum = &_spectrum[0];
112 #if USE_FFTW3
113  fftw_complex * complexOutput = &_complexOutput[0];
114  if (mIsComplex)
115  {
116  for (unsigned long i=0; i<mFrameSize; i++)
117  {
118  complexOutput[i][0] = input[i<<1];
119  complexOutput[i][1] = input[(i<<1)+1];
120  }
121  fftw_execute(_plan);
122  for (unsigned long i=0; i<mFrameSize*2; i+=2)
123  {
124  spectrum[i] = complexOutput[i/2][0];
125  spectrum[i+1] = complexOutput[i/2][1];
126  }
127  }
128  else
129  {
130  for (unsigned long i=0; i<mFrameSize; i++)
131  {
132  _realInput[i] = input[i];
133  }
134  fftw_execute(_plan);
135  for (int i=0; i<mFrameSize; i+=2)
136  {
137  spectrum[i] = complexOutput[i/2][0];
138  spectrum[i+1] = - complexOutput[i/2][1];
139  }
140  for (int i=1; i<mFrameSize; i+=2)
141  {
142  unsigned j = mFrameSize*2-i;
143  spectrum[j] = complexOutput[i/2][0];
144  spectrum[j+1] = complexOutput[i/2][1];
145  }
146  }
147 #else
148  if (mIsComplex)
149  {
150  for (unsigned long i=0; i<mFrameSize*2; i++)
151  spectrum[i] = input[i];
152  }
153  else
154  {
155  for (unsigned long i=0; i<mFrameSize*2; i+=2)
156  {
157  spectrum[i] = input[i/2];
158  spectrum[i+1] = 0.0;
159  }
160  }
161  double * data = spectrum-1; // Hack to use Matlab indeces TOREMOVE?
162  doFourierTransform(data, mFrameSize);
163 #endif
164 }
165 
166 void FourierTransform::doIt(const float * input)
167 {
168  doItGeneric(input);
169 }
170 void FourierTransform::doIt(const double * input)
171 {
172  doItGeneric(input);
173 }
174 
175 //-------------------------------------------------------------------------
176