comparison data/fft/FFTapi.cpp @ 227:e919a2b97c2a

* Fix in-house FFT provision
author Chris Cannam
date Mon, 12 Feb 2007 11:46:31 +0000
parents a867be73b638
children db946591a391
comparison
equal deleted inserted replaced
226:a867be73b638 227:e919a2b97c2a
17 #include "FFTapi.h" 17 #include "FFTapi.h"
18 18
19 #ifndef HAVE_FFTW3F 19 #ifndef HAVE_FFTW3F
20 20
21 #include <cmath> 21 #include <cmath>
22 #include <iostream>
22 23
23 void 24 void
24 fft(unsigned int n, bool inverse, double *ri, double *ii, double *ro, double *io) 25 fft(unsigned int n, bool inverse, double *ri, double *ii, double *ro, double *io)
25 { 26 {
26 if (!ri || !ro || !io) return; 27 if (!ri || !ro || !io) return;
42 bits = i; 43 bits = i;
43 break; 44 break;
44 } 45 }
45 } 46 }
46 47
47 static unsigned int tableSize = 0; 48 int *table = new int[n];
48 static int *table = 0; 49
49 50 for (i = 0; i < n; ++i) {
50 if (tableSize != n) {
51
52 delete[] table;
53
54 table = new int[n];
55
56 for (i = 0; i < n; ++i) {
57 51
58 m = i; 52 m = i;
59 53
60 for (j = k = 0; j < bits; ++j) { 54 for (j = k = 0; j < bits; ++j) {
61 k = (k << 1) | (m & 1); 55 k = (k << 1) | (m & 1);
62 m >>= 1; 56 m >>= 1;
63 } 57 }
64 58
65 table[i] = k; 59 table[i] = k;
66 }
67
68 tableSize = n;
69 } 60 }
70 61
71 if (ii) { 62 if (ii) {
72 for (i = 0; i < n; ++i) { 63 for (i = 0; i < n; ++i) {
73 ro[table[i]] = ri[i]; 64 ro[table[i]] = ri[i];
123 } 114 }
124 115
125 blockEnd = blockSize; 116 blockEnd = blockSize;
126 } 117 }
127 118
119 /* fftw doesn't normalise, so nor will we
120
128 if (inverse) { 121 if (inverse) {
129 122
130 double denom = (double)n; 123 double denom = (double)n;
131 124
132 for (i = 0; i < n; i++) { 125 for (i = 0; i < n; i++) {
133 ro[i] /= denom; 126 ro[i] /= denom;
134 io[i] /= denom; 127 io[i] /= denom;
135 } 128 }
136 } 129 }
130 */
131 delete[] table;
137 } 132 }
138 133
139 struct fftf_plan_ { 134 struct fftf_plan_ {
140 int size; 135 int size;
141 int inverse; 136 int inverse;
178 } 173 }
179 174
180 void 175 void
181 fftf_execute(const fftf_plan p) 176 fftf_execute(const fftf_plan p)
182 { 177 {
183 //!!! doesn't work yet for inverse transforms
184
185 float *real = p->real; 178 float *real = p->real;
186 fftf_complex *cplx = p->cplx; 179 fftf_complex *cplx = p->cplx;
187 int n = p->size; 180 int n = p->size;
188 int inv = p->inverse; 181 int forward = !p->inverse;
189 182
190 double *ri = new double[n]; 183 double *ri = new double[n];
191 double *ro = new double[n]; 184 double *ro = new double[n];
192 double *io = new double[n]; 185 double *io = new double[n];
193 186
194 double *ii = 0; 187 double *ii = 0;
195 if (inv) ii = new double[n]; 188 if (!forward) ii = new double[n];
196 189
197 if (!inv) { 190 if (forward) {
198 for (int i = 0; i < n; ++i) { 191 for (int i = 0; i < n; ++i) {
199 ri[i] = real[i]; 192 ri[i] = real[i];
200 } 193 }
201 } else { 194 } else {
202 for (int i = 0; i < n/2+1; ++i) { 195 for (int i = 0; i < n/2+1; ++i) {
207 ii[n-i] = -ii[i]; 200 ii[n-i] = -ii[i];
208 } 201 }
209 } 202 }
210 } 203 }
211 204
212 fft(n, inv, ri, ii, ro, io); 205 fft(n, !forward, ri, ii, ro, io);
213 206
214 if (!inv) { 207 if (forward) {
215 for (int i = 0; i < n/2+1; ++i) { 208 for (int i = 0; i < n/2+1; ++i) {
216 cplx[i][0] = ro[i]; 209 cplx[i][0] = ro[i];
217 cplx[i][1] = io[i]; 210 cplx[i][1] = io[i];
218 } 211 }
219 } else { 212 } else {
223 } 216 }
224 217
225 delete[] ri; 218 delete[] ri;
226 delete[] ro; 219 delete[] ro;
227 delete[] io; 220 delete[] io;
228 delete[] ii; 221 if (ii) delete[] ii;
229 } 222 }
230 223
231 #endif 224 #endif