26 #define SWAP(a,b) do {double swaptemp=(a); (a)=(b); (b)=swaptemp;} while(false)
32 : mFrameSize(framesize)
33 , mIsComplex(isComplex)
35 _spectrum.resize(2*mFrameSize);
37 _realInput = (
double*) fftw_malloc(
sizeof(
double) * mFrameSize);
38 _complexOutput = (fftw_complex*) fftw_malloc(
sizeof(fftw_complex) * mFrameSize);
40 fftw_import_system_wisdom();
42 _plan = fftw_plan_dft_1d(framesize, _complexOutput, _complexOutput, 1, FFTW_ESTIMATE);
44 _plan = fftw_plan_dft_r2c_1d(framesize, _realInput, _complexOutput, FFTW_ESTIMATE);
51 fftw_destroy_plan(_plan);
52 fftw_free(_realInput);
53 fftw_free(_complexOutput);
66 const unsigned long n=frameSize<<1;
67 for (
unsigned long i=1, j=1; i<n; i+=2) {
69 SWAP(data[j], data[i]);
70 SWAP(data[j+1], data[i+1]);
72 unsigned long m=frameSize;
73 while (m >= 2 && j>m) {
81 for (
unsigned long mmax = 2; mmax<n; mmax = istep)
84 const double theta =(6.28318530717959/mmax);
85 double wtemp = sin(0.5*theta);
86 double wpr = -2.0*wtemp*wtemp;
87 double wpi = sin(theta);
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;
101 wr = wtemp*wpr - wi*wpi + wr;
103 wi = wi*wpr + wtemp*wpi + wi;
108 template <
typename T>
109 void FourierTransform::doItGeneric(
const T * input)
113 fftw_complex * complexOutput = &_complexOutput[0];
116 for (
unsigned long i=0; i<mFrameSize; i++)
118 complexOutput[i][0] = input[i<<1];
119 complexOutput[i][1] = input[(i<<1)+1];
122 for (
unsigned long i=0; i<mFrameSize*2; i+=2)
124 spectrum[i] = complexOutput[i/2][0];
125 spectrum[i+1] = complexOutput[i/2][1];
130 for (
unsigned long i=0; i<mFrameSize; i++)
132 _realInput[i] = input[i];
135 for (
int i=0; i<mFrameSize; i+=2)
137 spectrum[i] = complexOutput[i/2][0];
138 spectrum[i+1] = - complexOutput[i/2][1];
140 for (
int i=1; i<mFrameSize; i+=2)
142 unsigned j = mFrameSize*2-i;
143 spectrum[j] = complexOutput[i/2][0];
144 spectrum[j+1] = complexOutput[i/2][1];
150 for (
unsigned long i=0; i<mFrameSize*2; i++)
151 spectrum[i] = input[i];
155 for (
unsigned long i=0; i<mFrameSize*2; i+=2)
157 spectrum[i] = input[i/2];
161 double * data = spectrum-1;