1 /* This is the FFT routine taken from PureData, a great piece of
2 software by Miller S. Puckette.
3 http://crca.ucsd.edu/~msp/software.html */
4
5 /*
6 ** FFT and FHT routines
7 ** Copyright 1988, 1993; Ron Mayer
8 **
9 ** mayer_fht(fz,n);
10 ** Does a hartley transform of "n" points in the array "fz".
11 ** mayer_fft(n,real,imag)
12 ** Does a fourier transform of "n" points of the "real" and
13 ** "imag" arrays.
14 ** mayer_ifft(n,real,imag)
15 ** Does an inverse fourier transform of "n" points of the "real"
16 ** and "imag" arrays.
17 ** mayer_realfft(n,real)
18 ** Does a real-valued fourier transform of "n" points of the
19 ** "real" array. The real part of the transform ends
20 ** up in the first half of the array and the imaginary part of the
21 ** transform ends up in the second half of the array.
22 ** mayer_realifft(n,real)
23 ** The inverse of the realfft() routine above.
24 **
25 **
26 ** NOTE: This routine uses at least 2 patented algorithms, and may be
27 ** under the restrictions of a bunch of different organizations.
28 ** Although I wrote it completely myself, it is kind of a derivative
29 ** of a routine I once authored and released under the GPL, so it
30 ** may fall under the free software foundation's restrictions;
31 ** it was worked on as a Stanford Univ project, so they claim
32 ** some rights to it; it was further optimized at work here, so
33 ** I think this company claims parts of it. The patents are
34 ** held by R. Bracewell (the FHT algorithm) and O. Buneman (the
35 ** trig generator), both at Stanford Univ.
36 ** If it were up to me, I'd say go do whatever you want with it;
37 ** but it would be polite to give credit to the following people
38 ** if you use this anywhere:
39 ** Euler - probable inventor of the fourier transform.
40 ** Gauss - probable inventor of the FFT.
41 ** Hartley - probable inventor of the hartley transform.
42 ** Buneman - for a really cool trig generator
43 ** Mayer(me) - for authoring this particular version and
44 ** including all the optimizations in one package.
45 ** Thanks,
46 ** Ron Mayer; mayer@acuson.com
47 **
48 */
49
50 /* This is a slightly modified version of Mayer's contribution; write
51 * msp@ucsd.edu for the original code. Kudos to Mayer for a fine piece
52 * of work. -msp
53 */
54
55 #define REAL float
56 #define GOOD_TRIG
57
58 #ifdef GOOD_TRIG
59 #else
60 #define FAST_TRIG
61 #endif
62
63 #if defined(GOOD_TRIG)
64 #define FHT_SWAP(a,b,t) {(t)=(a);(a)=(b);(b)=(t);}
65 #define TRIG_VARS \
66 int t_lam=0;
67 #define TRIG_INIT(k,c,s) \
68 { \
69 int i; \
70 for (i=2 ; i<=k ; i++) \
71 {coswrk[i]=costab[i];sinwrk[i]=sintab[i];} \
72 t_lam = 0; \
73 c = 1; \
74 s = 0; \
75 }
76 #define TRIG_NEXT(k,c,s) \
77 { \
78 int i,j; \
79 (t_lam)++; \
80 for (i=0 ; !((1<<i)&t_lam) ; i++); \
81 i = k-i; \
82 s = sinwrk[i]; \
83 c = coswrk[i]; \
84 if (i>1) \
85 { \
86 for (j=k-i+2 ; (1<<j)&t_lam ; j++); \
87 j = k - j; \
88 sinwrk[i] = halsec[i] * (sinwrk[i-1] + sinwrk[j]); \
89 coswrk[i] = halsec[i] * (coswrk[i-1] + coswrk[j]); \
90 } \
91 }
92 #define TRIG_RESET(k,c,s)
93 #endif
94
95 #if defined(FAST_TRIG)
96 #define TRIG_VARS \
97 REAL t_c,t_s;
98 #define TRIG_INIT(k,c,s) \
99 { \
100 t_c = costab[k]; \
101 t_s = sintab[k]; \
102 c = 1; \
103 s = 0; \
104 }
105 #define TRIG_NEXT(k,c,s) \
106 { \
107 REAL t = c; \
108 c = t*t_c - s*t_s; \
109 s = t*t_s + s*t_c; \
110 }
111 #define TRIG_RESET(k,c,s)
112 #endif
113
114 static REAL halsec[20]=
115 {
116 0,
117 0,
118 .54119610014619698439972320536638942006107206337801,
119 .50979557910415916894193980398784391368261849190893,
120 .50241928618815570551167011928012092247859337193963,
121 .50060299823519630134550410676638239611758632599591,
122 .50015063602065098821477101271097658495974913010340,
123 .50003765191554772296778139077905492847503165398345,
124 .50000941253588775676512870469186533538523133757983,
125 .50000235310628608051401267171204408939326297376426,
126 .50000058827484117879868526730916804925780637276181,
127 .50000014706860214875463798283871198206179118093251,
128 .50000003676714377807315864400643020315103490883972,
129 .50000000919178552207366560348853455333939112569380,
130 .50000000229794635411562887767906868558991922348920,
131 .50000000057448658687873302235147272458812263401372
132 };
133 static REAL costab[20]=
134 {
135 .00000000000000000000000000000000000000000000000000,
136 .70710678118654752440084436210484903928483593768847,
137 .92387953251128675612818318939678828682241662586364,
138 .98078528040323044912618223613423903697393373089333,
139 .99518472667219688624483695310947992157547486872985,
140 .99879545620517239271477160475910069444320361470461,
141 .99969881869620422011576564966617219685006108125772,
142 .99992470183914454092164649119638322435060646880221,
143 .99998117528260114265699043772856771617391725094433,
144 .99999529380957617151158012570011989955298763362218,
145 .99999882345170190992902571017152601904826792288976,
146 .99999970586288221916022821773876567711626389934930,
147 .99999992646571785114473148070738785694820115568892,
148 .99999998161642929380834691540290971450507605124278,
149 .99999999540410731289097193313960614895889430318945,
150 .99999999885102682756267330779455410840053741619428
151 };
152 static REAL sintab[20]=
153 {
154 1.0000000000000000000000000000000000000000000000000,
155 .70710678118654752440084436210484903928483593768846,
156 .38268343236508977172845998403039886676134456248561,
157 .19509032201612826784828486847702224092769161775195,
158 .09801714032956060199419556388864184586113667316749,
159 .04906767432741801425495497694268265831474536302574,
160 .02454122852291228803173452945928292506546611923944,
161 .01227153828571992607940826195100321214037231959176,
162 .00613588464915447535964023459037258091705788631738,
163 .00306795676296597627014536549091984251894461021344,
164 .00153398018628476561230369715026407907995486457522,
165 .00076699031874270452693856835794857664314091945205,
166 .00038349518757139558907246168118138126339502603495,
167 .00019174759731070330743990956198900093346887403385,
168 .00009587379909597734587051721097647635118706561284,
169 .00004793689960306688454900399049465887274686668768
170 };
171 static REAL coswrk[20]=
172 {
173 .00000000000000000000000000000000000000000000000000,
174 .70710678118654752440084436210484903928483593768847,
175 .92387953251128675612818318939678828682241662586364,
176 .98078528040323044912618223613423903697393373089333,
177 .99518472667219688624483695310947992157547486872985,
178 .99879545620517239271477160475910069444320361470461,
179 .99969881869620422011576564966617219685006108125772,
180 .99992470183914454092164649119638322435060646880221,
181 .99998117528260114265699043772856771617391725094433,
182 .99999529380957617151158012570011989955298763362218,
183 .99999882345170190992902571017152601904826792288976,
184 .99999970586288221916022821773876567711626389934930,
185 .99999992646571785114473148070738785694820115568892,
186 .99999998161642929380834691540290971450507605124278,
187 .99999999540410731289097193313960614895889430318945,
188 .99999999885102682756267330779455410840053741619428
189 };
190 static REAL sinwrk[20]=
191 {
192 1.0000000000000000000000000000000000000000000000000,
193 .70710678118654752440084436210484903928483593768846,
194 .38268343236508977172845998403039886676134456248561,
195 .19509032201612826784828486847702224092769161775195,
196 .09801714032956060199419556388864184586113667316749,
197 .04906767432741801425495497694268265831474536302574,
198 .02454122852291228803173452945928292506546611923944,
199 .01227153828571992607940826195100321214037231959176,
200 .00613588464915447535964023459037258091705788631738,
201 .00306795676296597627014536549091984251894461021344,
202 .00153398018628476561230369715026407907995486457522,
203 .00076699031874270452693856835794857664314091945205,
204 .00038349518757139558907246168118138126339502603495,
205 .00019174759731070330743990956198900093346887403385,
206 .00009587379909597734587051721097647635118706561284,
207 .00004793689960306688454900399049465887274686668768
208 };
209
210
211 #define SQRT2_2 0.70710678118654752440084436210484
212 #define SQRT2 2*0.70710678118654752440084436210484
213
mayer_fht(REAL * fz,int n)214 void mayer_fht(REAL *fz, int n)
215 {
216 /* REAL a,b;
217 REAL c1,s1,s2,c2,s3,c3,s4,c4;
218 REAL f0,g0,f1,g1,f2,g2,f3,g3; */
219 int k,k1,k2,k3,k4,kx;
220 REAL *fi,*fn,*gi;
221 TRIG_VARS;
222
223 for (k1=1,k2=0;k1<n;k1++)
224 {
225 REAL aa;
226 for (k=n>>1; (!((k2^=k)&k)); k>>=1);
227 if (k1>k2)
228 {
229 aa=fz[k1];fz[k1]=fz[k2];fz[k2]=aa;
230 }
231 }
232 for ( k=0 ; (1<<k)<n ; k++ );
233 k &= 1;
234 if (k==0)
235 {
236 for (fi=fz,fn=fz+n;fi<fn;fi+=4)
237 {
238 REAL f0,f1,f2,f3;
239 f1 = fi[0 ]-fi[1 ];
240 f0 = fi[0 ]+fi[1 ];
241 f3 = fi[2 ]-fi[3 ];
242 f2 = fi[2 ]+fi[3 ];
243 fi[2 ] = (f0-f2);
244 fi[0 ] = (f0+f2);
245 fi[3 ] = (f1-f3);
246 fi[1 ] = (f1+f3);
247 }
248 }
249 else
250 {
251 for (fi=fz,fn=fz+n,gi=fi+1;fi<fn;fi+=8,gi+=8)
252 {
253 REAL bs1,bc1,bs2,bc2,bs3,bc3,bs4,bc4,
254 bg0,bf0,bf1,bg1,bf2,bg2,bf3,bg3;
255 bc1 = fi[0 ] - gi[0 ];
256 bs1 = fi[0 ] + gi[0 ];
257 bc2 = fi[2 ] - gi[2 ];
258 bs2 = fi[2 ] + gi[2 ];
259 bc3 = fi[4 ] - gi[4 ];
260 bs3 = fi[4 ] + gi[4 ];
261 bc4 = fi[6 ] - gi[6 ];
262 bs4 = fi[6 ] + gi[6 ];
263 bf1 = (bs1 - bs2);
264 bf0 = (bs1 + bs2);
265 bg1 = (bc1 - bc2);
266 bg0 = (bc1 + bc2);
267 bf3 = (bs3 - bs4);
268 bf2 = (bs3 + bs4);
269 bg3 = SQRT2*bc4;
270 bg2 = SQRT2*bc3;
271 fi[4 ] = bf0 - bf2;
272 fi[0 ] = bf0 + bf2;
273 fi[6 ] = bf1 - bf3;
274 fi[2 ] = bf1 + bf3;
275 gi[4 ] = bg0 - bg2;
276 gi[0 ] = bg0 + bg2;
277 gi[6 ] = bg1 - bg3;
278 gi[2 ] = bg1 + bg3;
279 }
280 }
281 if (n<16) return;
282
283 do
284 {
285 REAL s1,c1;
286 int ii;
287 k += 2;
288 k1 = 1 << k;
289 k2 = k1 << 1;
290 k4 = k2 << 1;
291 k3 = k2 + k1;
292 kx = k1 >> 1;
293 fi = fz;
294 gi = fi + kx;
295 fn = fz + n;
296 do
297 {
298 REAL g0,f0,f1,g1,f2,g2,f3,g3;
299 f1 = fi[0 ] - fi[k1];
300 f0 = fi[0 ] + fi[k1];
301 f3 = fi[k2] - fi[k3];
302 f2 = fi[k2] + fi[k3];
303 fi[k2] = f0 - f2;
304 fi[0 ] = f0 + f2;
305 fi[k3] = f1 - f3;
306 fi[k1] = f1 + f3;
307 g1 = gi[0 ] - gi[k1];
308 g0 = gi[0 ] + gi[k1];
309 g3 = SQRT2 * gi[k3];
310 g2 = SQRT2 * gi[k2];
311 gi[k2] = g0 - g2;
312 gi[0 ] = g0 + g2;
313 gi[k3] = g1 - g3;
314 gi[k1] = g1 + g3;
315 gi += k4;
316 fi += k4;
317 } while (fi<fn);
318 TRIG_INIT(k,c1,s1);
319 for (ii=1;ii<kx;ii++)
320 {
321 REAL c2,s2;
322 TRIG_NEXT(k,c1,s1);
323 c2 = c1*c1 - s1*s1;
324 s2 = 2*(c1*s1);
325 fn = fz + n;
326 fi = fz +ii;
327 gi = fz +k1-ii;
328 do
329 {
330 REAL a,b,g0,f0,f1,g1,f2,g2,f3,g3;
331 b = s2*fi[k1] - c2*gi[k1];
332 a = c2*fi[k1] + s2*gi[k1];
333 f1 = fi[0 ] - a;
334 f0 = fi[0 ] + a;
335 g1 = gi[0 ] - b;
336 g0 = gi[0 ] + b;
337 b = s2*fi[k3] - c2*gi[k3];
338 a = c2*fi[k3] + s2*gi[k3];
339 f3 = fi[k2] - a;
340 f2 = fi[k2] + a;
341 g3 = gi[k2] - b;
342 g2 = gi[k2] + b;
343 b = s1*f2 - c1*g3;
344 a = c1*f2 + s1*g3;
345 fi[k2] = f0 - a;
346 fi[0 ] = f0 + a;
347 gi[k3] = g1 - b;
348 gi[k1] = g1 + b;
349 b = c1*g2 - s1*f3;
350 a = s1*g2 + c1*f3;
351 gi[k2] = g0 - a;
352 gi[0 ] = g0 + a;
353 fi[k3] = f1 - b;
354 fi[k1] = f1 + b;
355 gi += k4;
356 fi += k4;
357 } while (fi<fn);
358 }
359 TRIG_RESET(k,c1,s1);
360 } while (k4<n);
361 }
362
mayer_fft(int n,REAL * real,REAL * imag)363 void mayer_fft(int n, REAL *real, REAL *imag)
364 {
365 REAL a,b,c,d;
366 REAL q,r,s,t;
367 int i,j,k;
368 for (i=1,j=n-1,k=n/2;i<k;i++,j--) {
369 a = real[i]; b = real[j]; q=a+b; r=a-b;
370 c = imag[i]; d = imag[j]; s=c+d; t=c-d;
371 real[i] = (q+t)*.5; real[j] = (q-t)*.5;
372 imag[i] = (s-r)*.5; imag[j] = (s+r)*.5;
373 }
374 mayer_fht(real,n);
375 mayer_fht(imag,n);
376 }
377
mayer_ifft(int n,REAL * real,REAL * imag)378 void mayer_ifft(int n, REAL *real, REAL *imag)
379 {
380 REAL a,b,c,d;
381 REAL q,r,s,t;
382 int i,j,k;
383 mayer_fht(real,n);
384 mayer_fht(imag,n);
385 for (i=1,j=n-1,k=n/2;i<k;i++,j--) {
386 a = real[i]; b = real[j]; q=a+b; r=a-b;
387 c = imag[i]; d = imag[j]; s=c+d; t=c-d;
388 imag[i] = (s+r)*0.5; imag[j] = (s-r)*0.5;
389 real[i] = (q-t)*0.5; real[j] = (q+t)*0.5;
390 }
391 }
392
mayer_realfft(int n,REAL * real)393 void mayer_realfft(int n, REAL *real)
394 {
395 REAL a,b;
396 int i,j,k;
397
398 mayer_fht(real,n);
399 for (i=1,j=n-1,k=n/2;i<k;i++,j--) {
400 a = real[i];
401 b = real[j];
402 real[j] = (a-b)*0.5;
403 real[i] = (a+b)*0.5;
404 }
405 }
406
mayer_realifft(int n,REAL * real)407 void mayer_realifft(int n, REAL *real)
408 {
409 REAL a,b;
410 int i,j,k;
411
412 for (i=1,j=n-1,k=n/2;i<k;i++,j--) {
413 a = real[i];
414 b = real[j];
415 real[j] = (a-b);
416 real[i] = (a+b);
417 }
418 mayer_fht(real,n);
419 }
420