Mercurial > hg > svcore
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 |