summaryrefslogtreecommitdiff
path: root/apps/plugins/pdbox/PDa/src/d_mayer_fft.c
diff options
context:
space:
mode:
Diffstat (limited to 'apps/plugins/pdbox/PDa/src/d_mayer_fft.c')
-rw-r--r--apps/plugins/pdbox/PDa/src/d_mayer_fft.c838
1 files changed, 838 insertions, 0 deletions
diff --git a/apps/plugins/pdbox/PDa/src/d_mayer_fft.c b/apps/plugins/pdbox/PDa/src/d_mayer_fft.c
new file mode 100644
index 0000000000..7c29c6e7c8
--- /dev/null
+++ b/apps/plugins/pdbox/PDa/src/d_mayer_fft.c
@@ -0,0 +1,838 @@
1/*
2** FFT and FHT routines
3** Copyright 1988, 1993; Ron Mayer
4**
5** mayer_fht(fz,n);
6** Does a hartley transform of "n" points in the array "fz".
7** mayer_fft(n,real,imag)
8** Does a fourier transform of "n" points of the "real" and
9** "imag" arrays.
10** mayer_ifft(n,real,imag)
11** Does an inverse fourier transform of "n" points of the "real"
12** and "imag" arrays.
13** mayer_realfft(n,real)
14** Does a real-valued fourier transform of "n" points of the
15** "real" array. The real part of the transform ends
16** up in the first half of the array and the imaginary part of the
17** transform ends up in the second half of the array.
18** mayer_realifft(n,real)
19** The inverse of the realfft() routine above.
20**
21**
22** NOTE: This routine uses at least 2 patented algorithms, and may be
23** under the restrictions of a bunch of different organizations.
24** Although I wrote it completely myself, it is kind of a derivative
25** of a routine I once authored and released under the GPL, so it
26** may fall under the free software foundation's restrictions;
27** it was worked on as a Stanford Univ project, so they claim
28** some rights to it; it was further optimized at work here, so
29** I think this company claims parts of it. The patents are
30** held by R. Bracewell (the FHT algorithm) and O. Buneman (the
31** trig generator), both at Stanford Univ.
32** If it were up to me, I'd say go do whatever you want with it;
33** but it would be polite to give credit to the following people
34** if you use this anywhere:
35** Euler - probable inventor of the fourier transform.
36** Gauss - probable inventor of the FFT.
37** Hartley - probable inventor of the hartley transform.
38** Buneman - for a really cool trig generator
39** Mayer(me) - for authoring this particular version and
40** including all the optimizations in one package.
41** Thanks,
42** Ron Mayer; mayer@acuson.com
43**
44*/
45
46/* This is a slightly modified version of Mayer's contribution; write
47* msp@ucsd.edu for the original code. Kudos to Mayer for a fine piece
48* of work. -msp
49*/
50
51#ifdef MSW
52#pragma warning( disable : 4305 ) /* uncast const double to float */
53#pragma warning( disable : 4244 ) /* uncast double to float */
54#pragma warning( disable : 4101 ) /* unused local variables */
55#endif
56
57#define REAL float
58#define GOOD_TRIG
59
60#ifdef GOOD_TRIG
61#else
62#define FAST_TRIG
63#endif
64
65#if defined(GOOD_TRIG)
66#define FHT_SWAP(a,b,t) {(t)=(a);(a)=(b);(b)=(t);}
67#define TRIG_VARS \
68 int t_lam=0;
69#define TRIG_INIT(k,c,s) \
70 { \
71 int i; \
72 for (i=2 ; i<=k ; i++) \
73 {coswrk[i]=costab[i];sinwrk[i]=sintab[i];} \
74 t_lam = 0; \
75 c = 1; \
76 s = 0; \
77 }
78#define TRIG_NEXT(k,c,s) \
79 { \
80 int i,j; \
81 (t_lam)++; \
82 for (i=0 ; !((1<<i)&t_lam) ; i++); \
83 i = k-i; \
84 s = sinwrk[i]; \
85 c = coswrk[i]; \
86 if (i>1) \
87 { \
88 for (j=k-i+2 ; (1<<j)&t_lam ; j++); \
89 j = k - j; \
90 sinwrk[i] = halsec[i] * (sinwrk[i-1] + sinwrk[j]); \
91 coswrk[i] = halsec[i] * (coswrk[i-1] + coswrk[j]); \
92 } \
93 }
94#define TRIG_RESET(k,c,s)
95#endif
96
97#if defined(FAST_TRIG)
98#define TRIG_VARS \
99 REAL t_c,t_s;
100#define TRIG_INIT(k,c,s) \
101 { \
102 t_c = costab[k]; \
103 t_s = sintab[k]; \
104 c = 1; \
105 s = 0; \
106 }
107#define TRIG_NEXT(k,c,s) \
108 { \
109 REAL t = c; \
110 c = t*t_c - s*t_s; \
111 s = t*t_s + s*t_c; \
112 }
113#define TRIG_RESET(k,c,s)
114#endif
115
116static REAL halsec[20]=
117 {
118 0,
119 0,
120 .54119610014619698439972320536638942006107206337801,
121 .50979557910415916894193980398784391368261849190893,
122 .50241928618815570551167011928012092247859337193963,
123 .50060299823519630134550410676638239611758632599591,
124 .50015063602065098821477101271097658495974913010340,
125 .50003765191554772296778139077905492847503165398345,
126 .50000941253588775676512870469186533538523133757983,
127 .50000235310628608051401267171204408939326297376426,
128 .50000058827484117879868526730916804925780637276181,
129 .50000014706860214875463798283871198206179118093251,
130 .50000003676714377807315864400643020315103490883972,
131 .50000000919178552207366560348853455333939112569380,
132 .50000000229794635411562887767906868558991922348920,
133 .50000000057448658687873302235147272458812263401372
134 };
135static REAL costab[20]=
136 {
137 .00000000000000000000000000000000000000000000000000,
138 .70710678118654752440084436210484903928483593768847,
139 .92387953251128675612818318939678828682241662586364,
140 .98078528040323044912618223613423903697393373089333,
141 .99518472667219688624483695310947992157547486872985,
142 .99879545620517239271477160475910069444320361470461,
143 .99969881869620422011576564966617219685006108125772,
144 .99992470183914454092164649119638322435060646880221,
145 .99998117528260114265699043772856771617391725094433,
146 .99999529380957617151158012570011989955298763362218,
147 .99999882345170190992902571017152601904826792288976,
148 .99999970586288221916022821773876567711626389934930,
149 .99999992646571785114473148070738785694820115568892,
150 .99999998161642929380834691540290971450507605124278,
151 .99999999540410731289097193313960614895889430318945,
152 .99999999885102682756267330779455410840053741619428
153 };
154static REAL sintab[20]=
155 {
156 1.0000000000000000000000000000000000000000000000000,
157 .70710678118654752440084436210484903928483593768846,
158 .38268343236508977172845998403039886676134456248561,
159 .19509032201612826784828486847702224092769161775195,
160 .09801714032956060199419556388864184586113667316749,
161 .04906767432741801425495497694268265831474536302574,
162 .02454122852291228803173452945928292506546611923944,
163 .01227153828571992607940826195100321214037231959176,
164 .00613588464915447535964023459037258091705788631738,
165 .00306795676296597627014536549091984251894461021344,
166 .00153398018628476561230369715026407907995486457522,
167 .00076699031874270452693856835794857664314091945205,
168 .00038349518757139558907246168118138126339502603495,
169 .00019174759731070330743990956198900093346887403385,
170 .00009587379909597734587051721097647635118706561284,
171 .00004793689960306688454900399049465887274686668768
172 };
173static REAL coswrk[20]=
174 {
175 .00000000000000000000000000000000000000000000000000,
176 .70710678118654752440084436210484903928483593768847,
177 .92387953251128675612818318939678828682241662586364,
178 .98078528040323044912618223613423903697393373089333,
179 .99518472667219688624483695310947992157547486872985,
180 .99879545620517239271477160475910069444320361470461,
181 .99969881869620422011576564966617219685006108125772,
182 .99992470183914454092164649119638322435060646880221,
183 .99998117528260114265699043772856771617391725094433,
184 .99999529380957617151158012570011989955298763362218,
185 .99999882345170190992902571017152601904826792288976,
186 .99999970586288221916022821773876567711626389934930,
187 .99999992646571785114473148070738785694820115568892,
188 .99999998161642929380834691540290971450507605124278,
189 .99999999540410731289097193313960614895889430318945,
190 .99999999885102682756267330779455410840053741619428
191 };
192static REAL sinwrk[20]=
193 {
194 1.0000000000000000000000000000000000000000000000000,
195 .70710678118654752440084436210484903928483593768846,
196 .38268343236508977172845998403039886676134456248561,
197 .19509032201612826784828486847702224092769161775195,
198 .09801714032956060199419556388864184586113667316749,
199 .04906767432741801425495497694268265831474536302574,
200 .02454122852291228803173452945928292506546611923944,
201 .01227153828571992607940826195100321214037231959176,
202 .00613588464915447535964023459037258091705788631738,
203 .00306795676296597627014536549091984251894461021344,
204 .00153398018628476561230369715026407907995486457522,
205 .00076699031874270452693856835794857664314091945205,
206 .00038349518757139558907246168118138126339502603495,
207 .00019174759731070330743990956198900093346887403385,
208 .00009587379909597734587051721097647635118706561284,
209 .00004793689960306688454900399049465887274686668768
210 };
211
212
213#define SQRT2_2 0.70710678118654752440084436210484
214#define SQRT2 2*0.70710678118654752440084436210484
215
216void mayer_fht(REAL *fz, int n)
217{
218/* REAL a,b;
219REAL c1,s1,s2,c2,s3,c3,s4,c4;
220 REAL f0,g0,f1,g1,f2,g2,f3,g3; */
221 int k,k1,k2,k3,k4,kx;
222 REAL *fi,*fn,*gi;
223 TRIG_VARS;
224
225 for (k1=1,k2=0;k1<n;k1++)
226 {
227 REAL aa;
228 for (k=n>>1; (!((k2^=k)&k)); k>>=1);
229 if (k1>k2)
230 {
231 aa=fz[k1];fz[k1]=fz[k2];fz[k2]=aa;
232 }
233 }
234 for ( k=0 ; (1<<k)<n ; k++ );
235 k &= 1;
236 if (k==0)
237 {
238 for (fi=fz,fn=fz+n;fi<fn;fi+=4)
239 {
240 REAL f0,f1,f2,f3;
241 f1 = fi[0 ]-fi[1 ];
242 f0 = fi[0 ]+fi[1 ];
243 f3 = fi[2 ]-fi[3 ];
244 f2 = fi[2 ]+fi[3 ];
245 fi[2 ] = (f0-f2);
246 fi[0 ] = (f0+f2);
247 fi[3 ] = (f1-f3);
248 fi[1 ] = (f1+f3);
249 }
250 }
251 else
252 {
253 for (fi=fz,fn=fz+n,gi=fi+1;fi<fn;fi+=8,gi+=8)
254 {
255 REAL bs1,bc1,bs2,bc2,bs3,bc3,bs4,bc4,
256 bg0,bf0,bf1,bg1,bf2,bg2,bf3,bg3;
257 bc1 = fi[0 ] - gi[0 ];
258 bs1 = fi[0 ] + gi[0 ];
259 bc2 = fi[2 ] - gi[2 ];
260 bs2 = fi[2 ] + gi[2 ];
261 bc3 = fi[4 ] - gi[4 ];
262 bs3 = fi[4 ] + gi[4 ];
263 bc4 = fi[6 ] - gi[6 ];
264 bs4 = fi[6 ] + gi[6 ];
265 bf1 = (bs1 - bs2);
266 bf0 = (bs1 + bs2);
267 bg1 = (bc1 - bc2);
268 bg0 = (bc1 + bc2);
269 bf3 = (bs3 - bs4);
270 bf2 = (bs3 + bs4);
271 bg3 = SQRT2*bc4;
272 bg2 = SQRT2*bc3;
273 fi[4 ] = bf0 - bf2;
274 fi[0 ] = bf0 + bf2;
275 fi[6 ] = bf1 - bf3;
276 fi[2 ] = bf1 + bf3;
277 gi[4 ] = bg0 - bg2;
278 gi[0 ] = bg0 + bg2;
279 gi[6 ] = bg1 - bg3;
280 gi[2 ] = bg1 + bg3;
281 }
282 }
283 if (n<16) return;
284
285 do
286 {
287 REAL s1,c1;
288 int ii;
289 k += 2;
290 k1 = 1 << k;
291 k2 = k1 << 1;
292 k4 = k2 << 1;
293 k3 = k2 + k1;
294 kx = k1 >> 1;
295 fi = fz;
296 gi = fi + kx;
297 fn = fz + n;
298 do
299 {
300 REAL g0,f0,f1,g1,f2,g2,f3,g3;
301 f1 = fi[0 ] - fi[k1];
302 f0 = fi[0 ] + fi[k1];
303 f3 = fi[k2] - fi[k3];
304 f2 = fi[k2] + fi[k3];
305 fi[k2] = f0 - f2;
306 fi[0 ] = f0 + f2;
307 fi[k3] = f1 - f3;
308 fi[k1] = f1 + f3;
309 g1 = gi[0 ] - gi[k1];
310 g0 = gi[0 ] + gi[k1];
311 g3 = SQRT2 * gi[k3];
312 g2 = SQRT2 * gi[k2];
313 gi[k2] = g0 - g2;
314 gi[0 ] = g0 + g2;
315 gi[k3] = g1 - g3;
316 gi[k1] = g1 + g3;
317 gi += k4;
318 fi += k4;
319 } while (fi<fn);
320 TRIG_INIT(k,c1,s1);
321 for (ii=1;ii<kx;ii++)
322 {
323 REAL c2,s2;
324 TRIG_NEXT(k,c1,s1);
325 c2 = c1*c1 - s1*s1;
326 s2 = 2*(c1*s1);
327 fn = fz + n;
328 fi = fz +ii;
329 gi = fz +k1-ii;
330 do
331 {
332 REAL a,b,g0,f0,f1,g1,f2,g2,f3,g3;
333 b = s2*fi[k1] - c2*gi[k1];
334 a = c2*fi[k1] + s2*gi[k1];
335 f1 = fi[0 ] - a;
336 f0 = fi[0 ] + a;
337 g1 = gi[0 ] - b;
338 g0 = gi[0 ] + b;
339 b = s2*fi[k3] - c2*gi[k3];
340 a = c2*fi[k3] + s2*gi[k3];
341 f3 = fi[k2] - a;
342 f2 = fi[k2] + a;
343 g3 = gi[k2] - b;
344 g2 = gi[k2] + b;
345 b = s1*f2 - c1*g3;
346 a = c1*f2 + s1*g3;
347 fi[k2] = f0 - a;
348 fi[0 ] = f0 + a;
349 gi[k3] = g1 - b;
350 gi[k1] = g1 + b;
351 b = c1*g2 - s1*f3;
352 a = s1*g2 + c1*f3;
353 gi[k2] = g0 - a;
354 gi[0 ] = g0 + a;
355 fi[k3] = f1 - b;
356 fi[k1] = f1 + b;
357 gi += k4;
358 fi += k4;
359 } while (fi<fn);
360 }
361 TRIG_RESET(k,c1,s1);
362 } while (k4<n);
363}
364
365void mayer_fft(int n, REAL *real, REAL *imag)
366{
367 REAL a,b,c,d;
368 REAL q,r,s,t;
369 int i,j,k;
370 for (i=1,j=n-1,k=n/2;i<k;i++,j--) {
371 a = real[i]; b = real[j]; q=a+b; r=a-b;
372 c = imag[i]; d = imag[j]; s=c+d; t=c-d;
373 real[i] = (q+t)*.5; real[j] = (q-t)*.5;
374 imag[i] = (s-r)*.5; imag[j] = (s+r)*.5;
375 }
376 mayer_fht(real,n);
377 mayer_fht(imag,n);
378}
379
380void mayer_ifft(int n, REAL *real, REAL *imag)
381{
382 REAL a,b,c,d;
383 REAL q,r,s,t;
384 int i,j,k;
385 mayer_fht(real,n);
386 mayer_fht(imag,n);
387 for (i=1,j=n-1,k=n/2;i<k;i++,j--) {
388 a = real[i]; b = real[j]; q=a+b; r=a-b;
389 c = imag[i]; d = imag[j]; s=c+d; t=c-d;
390 imag[i] = (s+r)*0.5; imag[j] = (s-r)*0.5;
391 real[i] = (q-t)*0.5; real[j] = (q+t)*0.5;
392 }
393}
394
395void mayer_realfft(int n, REAL *real)
396{
397 REAL a,b,c,d;
398 int i,j,k;
399 mayer_fht(real,n);
400 for (i=1,j=n-1,k=n/2;i<k;i++,j--) {
401 a = real[i];
402 b = real[j];
403 real[j] = (a-b)*0.5;
404 real[i] = (a+b)*0.5;
405 }
406}
407
408void mayer_realifft(int n, REAL *real)
409{
410 REAL a,b,c,d;
411 int i,j,k;
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/*
421** FFT and FHT routines
422** Copyright 1988, 1993; Ron Mayer
423**
424** mayer_fht(fz,n);
425** Does a hartley transform of "n" points in the array "fz".
426** mayer_fft(n,real,imag)
427** Does a fourier transform of "n" points of the "real" and
428** "imag" arrays.
429** mayer_ifft(n,real,imag)
430** Does an inverse fourier transform of "n" points of the "real"
431** and "imag" arrays.
432** mayer_realfft(n,real)
433** Does a real-valued fourier transform of "n" points of the
434** "real" array. The real part of the transform ends
435** up in the first half of the array and the imaginary part of the
436** transform ends up in the second half of the array.
437** mayer_realifft(n,real)
438** The inverse of the realfft() routine above.
439**
440**
441** NOTE: This routine uses at least 2 patented algorithms, and may be
442** under the restrictions of a bunch of different organizations.
443** Although I wrote it completely myself, it is kind of a derivative
444** of a routine I once authored and released under the GPL, so it
445** may fall under the free software foundation's restrictions;
446** it was worked on as a Stanford Univ project, so they claim
447** some rights to it; it was further optimized at work here, so
448** I think this company claims parts of it. The patents are
449** held by R. Bracewell (the FHT algorithm) and O. Buneman (the
450** trig generator), both at Stanford Univ.
451** If it were up to me, I'd say go do whatever you want with it;
452** but it would be polite to give credit to the following people
453** if you use this anywhere:
454** Euler - probable inventor of the fourier transform.
455** Gauss - probable inventor of the FFT.
456** Hartley - probable inventor of the hartley transform.
457** Buneman - for a really cool trig generator
458** Mayer(me) - for authoring this particular version and
459** including all the optimizations in one package.
460** Thanks,
461** Ron Mayer; mayer@acuson.com
462**
463*/
464
465/* This is a slightly modified version of Mayer's contribution; write
466* msp@ucsd.edu for the original code. Kudos to Mayer for a fine piece
467* of work. -msp
468*/
469
470#ifdef MSW
471#pragma warning( disable : 4305 ) /* uncast const double to float */
472#pragma warning( disable : 4244 ) /* uncast double to float */
473#pragma warning( disable : 4101 ) /* unused local variables */
474#endif
475
476#define REAL float
477#define GOOD_TRIG
478
479#ifdef GOOD_TRIG
480#else
481#define FAST_TRIG
482#endif
483
484#if defined(GOOD_TRIG)
485#define FHT_SWAP(a,b,t) {(t)=(a);(a)=(b);(b)=(t);}
486#define TRIG_VARS \
487 int t_lam=0;
488#define TRIG_INIT(k,c,s) \
489 { \
490 int i; \
491 for (i=2 ; i<=k ; i++) \
492 {coswrk[i]=costab[i];sinwrk[i]=sintab[i];} \
493 t_lam = 0; \
494 c = 1; \
495 s = 0; \
496 }
497#define TRIG_NEXT(k,c,s) \
498 { \
499 int i,j; \
500 (t_lam)++; \
501 for (i=0 ; !((1<<i)&t_lam) ; i++); \
502 i = k-i; \
503 s = sinwrk[i]; \
504 c = coswrk[i]; \
505 if (i>1) \
506 { \
507 for (j=k-i+2 ; (1<<j)&t_lam ; j++); \
508 j = k - j; \
509 sinwrk[i] = halsec[i] * (sinwrk[i-1] + sinwrk[j]); \
510 coswrk[i] = halsec[i] * (coswrk[i-1] + coswrk[j]); \
511 } \
512 }
513#define TRIG_RESET(k,c,s)
514#endif
515
516#if defined(FAST_TRIG)
517#define TRIG_VARS \
518 REAL t_c,t_s;
519#define TRIG_INIT(k,c,s) \
520 { \
521 t_c = costab[k]; \
522 t_s = sintab[k]; \
523 c = 1; \
524 s = 0; \
525 }
526#define TRIG_NEXT(k,c,s) \
527 { \
528 REAL t = c; \
529 c = t*t_c - s*t_s; \
530 s = t*t_s + s*t_c; \
531 }
532#define TRIG_RESET(k,c,s)
533#endif
534
535static REAL halsec[20]=
536 {
537 0,
538 0,
539 .54119610014619698439972320536638942006107206337801,
540 .50979557910415916894193980398784391368261849190893,
541 .50241928618815570551167011928012092247859337193963,
542 .50060299823519630134550410676638239611758632599591,
543 .50015063602065098821477101271097658495974913010340,
544 .50003765191554772296778139077905492847503165398345,
545 .50000941253588775676512870469186533538523133757983,
546 .50000235310628608051401267171204408939326297376426,
547 .50000058827484117879868526730916804925780637276181,
548 .50000014706860214875463798283871198206179118093251,
549 .50000003676714377807315864400643020315103490883972,
550 .50000000919178552207366560348853455333939112569380,
551 .50000000229794635411562887767906868558991922348920,
552 .50000000057448658687873302235147272458812263401372
553 };
554static REAL costab[20]=
555 {
556 .00000000000000000000000000000000000000000000000000,
557 .70710678118654752440084436210484903928483593768847,
558 .92387953251128675612818318939678828682241662586364,
559 .98078528040323044912618223613423903697393373089333,
560 .99518472667219688624483695310947992157547486872985,
561 .99879545620517239271477160475910069444320361470461,
562 .99969881869620422011576564966617219685006108125772,
563 .99992470183914454092164649119638322435060646880221,
564 .99998117528260114265699043772856771617391725094433,
565 .99999529380957617151158012570011989955298763362218,
566 .99999882345170190992902571017152601904826792288976,
567 .99999970586288221916022821773876567711626389934930,
568 .99999992646571785114473148070738785694820115568892,
569 .99999998161642929380834691540290971450507605124278,
570 .99999999540410731289097193313960614895889430318945,
571 .99999999885102682756267330779455410840053741619428
572 };
573static REAL sintab[20]=
574 {
575 1.0000000000000000000000000000000000000000000000000,
576 .70710678118654752440084436210484903928483593768846,
577 .38268343236508977172845998403039886676134456248561,
578 .19509032201612826784828486847702224092769161775195,
579 .09801714032956060199419556388864184586113667316749,
580 .04906767432741801425495497694268265831474536302574,
581 .02454122852291228803173452945928292506546611923944,
582 .01227153828571992607940826195100321214037231959176,
583 .00613588464915447535964023459037258091705788631738,
584 .00306795676296597627014536549091984251894461021344,
585 .00153398018628476561230369715026407907995486457522,
586 .00076699031874270452693856835794857664314091945205,
587 .00038349518757139558907246168118138126339502603495,
588 .00019174759731070330743990956198900093346887403385,
589 .00009587379909597734587051721097647635118706561284,
590 .00004793689960306688454900399049465887274686668768
591 };
592static REAL coswrk[20]=
593 {
594 .00000000000000000000000000000000000000000000000000,
595 .70710678118654752440084436210484903928483593768847,
596 .92387953251128675612818318939678828682241662586364,
597 .98078528040323044912618223613423903697393373089333,
598 .99518472667219688624483695310947992157547486872985,
599 .99879545620517239271477160475910069444320361470461,
600 .99969881869620422011576564966617219685006108125772,
601 .99992470183914454092164649119638322435060646880221,
602 .99998117528260114265699043772856771617391725094433,
603 .99999529380957617151158012570011989955298763362218,
604 .99999882345170190992902571017152601904826792288976,
605 .99999970586288221916022821773876567711626389934930,
606 .99999992646571785114473148070738785694820115568892,
607 .99999998161642929380834691540290971450507605124278,
608 .99999999540410731289097193313960614895889430318945,
609 .99999999885102682756267330779455410840053741619428
610 };
611static REAL sinwrk[20]=
612 {
613 1.0000000000000000000000000000000000000000000000000,
614 .70710678118654752440084436210484903928483593768846,
615 .38268343236508977172845998403039886676134456248561,
616 .19509032201612826784828486847702224092769161775195,
617 .09801714032956060199419556388864184586113667316749,
618 .04906767432741801425495497694268265831474536302574,
619 .02454122852291228803173452945928292506546611923944,
620 .01227153828571992607940826195100321214037231959176,
621 .00613588464915447535964023459037258091705788631738,
622 .00306795676296597627014536549091984251894461021344,
623 .00153398018628476561230369715026407907995486457522,
624 .00076699031874270452693856835794857664314091945205,
625 .00038349518757139558907246168118138126339502603495,
626 .00019174759731070330743990956198900093346887403385,
627 .00009587379909597734587051721097647635118706561284,
628 .00004793689960306688454900399049465887274686668768
629 };
630
631
632#define SQRT2_2 0.70710678118654752440084436210484
633#define SQRT2 2*0.70710678118654752440084436210484
634
635void mayer_fht(REAL *fz, int n)
636{
637/* REAL a,b;
638REAL c1,s1,s2,c2,s3,c3,s4,c4;
639 REAL f0,g0,f1,g1,f2,g2,f3,g3; */
640 int k,k1,k2,k3,k4,kx;
641 REAL *fi,*fn,*gi;
642 TRIG_VARS;
643
644 for (k1=1,k2=0;k1<n;k1++)
645 {
646 REAL aa;
647 for (k=n>>1; (!((k2^=k)&k)); k>>=1);
648 if (k1>k2)
649 {
650 aa=fz[k1];fz[k1]=fz[k2];fz[k2]=aa;
651 }
652 }
653 for ( k=0 ; (1<<k)<n ; k++ );
654 k &= 1;
655 if (k==0)
656 {
657 for (fi=fz,fn=fz+n;fi<fn;fi+=4)
658 {
659 REAL f0,f1,f2,f3;
660 f1 = fi[0 ]-fi[1 ];
661 f0 = fi[0 ]+fi[1 ];
662 f3 = fi[2 ]-fi[3 ];
663 f2 = fi[2 ]+fi[3 ];
664 fi[2 ] = (f0-f2);
665 fi[0 ] = (f0+f2);
666 fi[3 ] = (f1-f3);
667 fi[1 ] = (f1+f3);
668 }
669 }
670 else
671 {
672 for (fi=fz,fn=fz+n,gi=fi+1;fi<fn;fi+=8,gi+=8)
673 {
674 REAL bs1,bc1,bs2,bc2,bs3,bc3,bs4,bc4,
675 bg0,bf0,bf1,bg1,bf2,bg2,bf3,bg3;
676 bc1 = fi[0 ] - gi[0 ];
677 bs1 = fi[0 ] + gi[0 ];
678 bc2 = fi[2 ] - gi[2 ];
679 bs2 = fi[2 ] + gi[2 ];
680 bc3 = fi[4 ] - gi[4 ];
681 bs3 = fi[4 ] + gi[4 ];
682 bc4 = fi[6 ] - gi[6 ];
683 bs4 = fi[6 ] + gi[6 ];
684 bf1 = (bs1 - bs2);
685 bf0 = (bs1 + bs2);
686 bg1 = (bc1 - bc2);
687 bg0 = (bc1 + bc2);
688 bf3 = (bs3 - bs4);
689 bf2 = (bs3 + bs4);
690 bg3 = SQRT2*bc4;
691 bg2 = SQRT2*bc3;
692 fi[4 ] = bf0 - bf2;
693 fi[0 ] = bf0 + bf2;
694 fi[6 ] = bf1 - bf3;
695 fi[2 ] = bf1 + bf3;
696 gi[4 ] = bg0 - bg2;
697 gi[0 ] = bg0 + bg2;
698 gi[6 ] = bg1 - bg3;
699 gi[2 ] = bg1 + bg3;
700 }
701 }
702 if (n<16) return;
703
704 do
705 {
706 REAL s1,c1;
707 int ii;
708 k += 2;
709 k1 = 1 << k;
710 k2 = k1 << 1;
711 k4 = k2 << 1;
712 k3 = k2 + k1;
713 kx = k1 >> 1;
714 fi = fz;
715 gi = fi + kx;
716 fn = fz + n;
717 do
718 {
719 REAL g0,f0,f1,g1,f2,g2,f3,g3;
720 f1 = fi[0 ] - fi[k1];
721 f0 = fi[0 ] + fi[k1];
722 f3 = fi[k2] - fi[k3];
723 f2 = fi[k2] + fi[k3];
724 fi[k2] = f0 - f2;
725 fi[0 ] = f0 + f2;
726 fi[k3] = f1 - f3;
727 fi[k1] = f1 + f3;
728 g1 = gi[0 ] - gi[k1];
729 g0 = gi[0 ] + gi[k1];
730 g3 = SQRT2 * gi[k3];
731 g2 = SQRT2 * gi[k2];
732 gi[k2] = g0 - g2;
733 gi[0 ] = g0 + g2;
734 gi[k3] = g1 - g3;
735 gi[k1] = g1 + g3;
736 gi += k4;
737 fi += k4;
738 } while (fi<fn);
739 TRIG_INIT(k,c1,s1);
740 for (ii=1;ii<kx;ii++)
741 {
742 REAL c2,s2;
743 TRIG_NEXT(k,c1,s1);
744 c2 = c1*c1 - s1*s1;
745 s2 = 2*(c1*s1);
746 fn = fz + n;
747 fi = fz +ii;
748 gi = fz +k1-ii;
749 do
750 {
751 REAL a,b,g0,f0,f1,g1,f2,g2,f3,g3;
752 b = s2*fi[k1] - c2*gi[k1];
753 a = c2*fi[k1] + s2*gi[k1];
754 f1 = fi[0 ] - a;
755 f0 = fi[0 ] + a;
756 g1 = gi[0 ] - b;
757 g0 = gi[0 ] + b;
758 b = s2*fi[k3] - c2*gi[k3];
759 a = c2*fi[k3] + s2*gi[k3];
760 f3 = fi[k2] - a;
761 f2 = fi[k2] + a;
762 g3 = gi[k2] - b;
763 g2 = gi[k2] + b;
764 b = s1*f2 - c1*g3;
765 a = c1*f2 + s1*g3;
766 fi[k2] = f0 - a;
767 fi[0 ] = f0 + a;
768 gi[k3] = g1 - b;
769 gi[k1] = g1 + b;
770 b = c1*g2 - s1*f3;
771 a = s1*g2 + c1*f3;
772 gi[k2] = g0 - a;
773 gi[0 ] = g0 + a;
774 fi[k3] = f1 - b;
775 fi[k1] = f1 + b;
776 gi += k4;
777 fi += k4;
778 } while (fi<fn);
779 }
780 TRIG_RESET(k,c1,s1);
781 } while (k4<n);
782}
783
784void mayer_fft(int n, REAL *real, REAL *imag)
785{
786 REAL a,b,c,d;
787 REAL q,r,s,t;
788 int i,j,k;
789 for (i=1,j=n-1,k=n/2;i<k;i++,j--) {
790 a = real[i]; b = real[j]; q=a+b; r=a-b;
791 c = imag[i]; d = imag[j]; s=c+d; t=c-d;
792 real[i] = (q+t)*.5; real[j] = (q-t)*.5;
793 imag[i] = (s-r)*.5; imag[j] = (s+r)*.5;
794 }
795 mayer_fht(real,n);
796 mayer_fht(imag,n);
797}
798
799void mayer_ifft(int n, REAL *real, REAL *imag)
800{
801 REAL a,b,c,d;
802 REAL q,r,s,t;
803 int i,j,k;
804 mayer_fht(real,n);
805 mayer_fht(imag,n);
806 for (i=1,j=n-1,k=n/2;i<k;i++,j--) {
807 a = real[i]; b = real[j]; q=a+b; r=a-b;
808 c = imag[i]; d = imag[j]; s=c+d; t=c-d;
809 imag[i] = (s+r)*0.5; imag[j] = (s-r)*0.5;
810 real[i] = (q-t)*0.5; real[j] = (q+t)*0.5;
811 }
812}
813
814void mayer_realfft(int n, REAL *real)
815{
816 REAL a,b,c,d;
817 int i,j,k;
818 mayer_fht(real,n);
819 for (i=1,j=n-1,k=n/2;i<k;i++,j--) {
820 a = real[i];
821 b = real[j];
822 real[j] = (a-b)*0.5;
823 real[i] = (a+b)*0.5;
824 }
825}
826
827void mayer_realifft(int n, REAL *real)
828{
829 REAL a,b,c,d;
830 int i,j,k;
831 for (i=1,j=n-1,k=n/2;i<k;i++,j--) {
832 a = real[i];
833 b = real[j];
834 real[j] = (a-b);
835 real[i] = (a+b);
836 }
837 mayer_fht(real,n);
838}
generated by cgit v1.2.3 (git 2.39.1) at 2024-09-12 17:50:10 +0000