source: Sophya/trunk/Poubelle/archTOI.old/fft_mayer.cc@ 2353

Last change on this file since 2353 was 394, checked in by ansari, 26 years ago

Integration detecteur d'etoiles DY
File size: 6.5 KB
Line 
1/*
2** FFT and FHT routines
3** Copyright 1988, 1993; Ron Mayer
4**
5** fht(fz,n);
6** Does a hartley transform of "n" points in the array "fz".
7** fft(n,real,imag)
8** Does a fourier transform of "n" points of the "real" and
9** "imag" arrays.
10** ifft(n,real,imag)
11** Does an inverse fourier transform of "n" points of the "real"
12** and "imag" arrays.
13** 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** 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#include "mayer_fft.h"
47
48#define GOOD_TRIG
49#include "trigtbl.h"
50
51char fht_version[] = "Brcwl-Hrtly-Ron-dbld";
52
53#define SQRT2_2 0.70710678118654752440084436210484
54#define SQRT2 2*0.70710678118654752440084436210484
55
56void fht(REAL *fz,int n)
57{
58 REAL a,b;
59 REAL c1,s1,s2,c2,s3,c3,s4,c4;
60 REAL f0,g0,f1,g1,f2,g2,f3,g3;
61 int i,k,k1,k2,k3,k4,kx;
62 REAL *fi,*fn,*gi;
63 TRIG_VARS;
64
65 for (k1=1,k2=0;k1<n;k1++)
66 {
67 REAL a;
68 for (k=n>>1; (!((k2^=k)&k)); k>>=1);
69 if (k1>k2)
70 {
71 a=fz[k1];fz[k1]=fz[k2];fz[k2]=a;
72 }
73 }
74 for ( k=0 ; (1<<k)<n ; k++ );
75 k &= 1;
76 if (k==0)
77 {
78 for (fi=fz,fn=fz+n;fi<fn;fi+=4)
79 {
80 REAL f0,f1,f2,f3;
81 f1 = fi[0 ]-fi[1 ];
82 f0 = fi[0 ]+fi[1 ];
83 f3 = fi[2 ]-fi[3 ];
84 f2 = fi[2 ]+fi[3 ];
85 fi[2 ] = (f0-f2);
86 fi[0 ] = (f0+f2);
87 fi[3 ] = (f1-f3);
88 fi[1 ] = (f1+f3);
89 }
90 }
91 else
92 {
93 for (fi=fz,fn=fz+n,gi=fi+1;fi<fn;fi+=8,gi+=8)
94 {
95 REAL s1,c1,s2,c2,s3,c3,s4,c4,g0,f0,f1,g1,f2,g2,f3,g3;
96 c1 = fi[0 ] - gi[0 ];
97 s1 = fi[0 ] + gi[0 ];
98 c2 = fi[2 ] - gi[2 ];
99 s2 = fi[2 ] + gi[2 ];
100 c3 = fi[4 ] - gi[4 ];
101 s3 = fi[4 ] + gi[4 ];
102 c4 = fi[6 ] - gi[6 ];
103 s4 = fi[6 ] + gi[6 ];
104 f1 = (s1 - s2);
105 f0 = (s1 + s2);
106 g1 = (c1 - c2);
107 g0 = (c1 + c2);
108 f3 = (s3 - s4);
109 f2 = (s3 + s4);
110 g3 = SQRT2*c4;
111 g2 = SQRT2*c3;
112 fi[4 ] = f0 - f2;
113 fi[0 ] = f0 + f2;
114 fi[6 ] = f1 - f3;
115 fi[2 ] = f1 + f3;
116 gi[4 ] = g0 - g2;
117 gi[0 ] = g0 + g2;
118 gi[6 ] = g1 - g3;
119 gi[2 ] = g1 + g3;
120 }
121 }
122 if (n<16) return;
123
124 do
125 {
126 REAL s1,c1;
127 k += 2;
128 k1 = 1 << k;
129 k2 = k1 << 1;
130 k4 = k2 << 1;
131 k3 = k2 + k1;
132 kx = k1 >> 1;
133 fi = fz;
134 gi = fi + kx;
135 fn = fz + n;
136 do
137 {
138 REAL g0,f0,f1,g1,f2,g2,f3,g3;
139 f1 = fi[0 ] - fi[k1];
140 f0 = fi[0 ] + fi[k1];
141 f3 = fi[k2] - fi[k3];
142 f2 = fi[k2] + fi[k3];
143 fi[k2] = f0 - f2;
144 fi[0 ] = f0 + f2;
145 fi[k3] = f1 - f3;
146 fi[k1] = f1 + f3;
147 g1 = gi[0 ] - gi[k1];
148 g0 = gi[0 ] + gi[k1];
149 g3 = SQRT2 * gi[k3];
150 g2 = SQRT2 * gi[k2];
151 gi[k2] = g0 - g2;
152 gi[0 ] = g0 + g2;
153 gi[k3] = g1 - g3;
154 gi[k1] = g1 + g3;
155 gi += k4;
156 fi += k4;
157 } while (fi<fn);
158 TRIG_INIT(k,c1,s1);
159 for (i=1;i<kx;i++)
160 {
161 REAL c2,s2;
162 TRIG_NEXT(k,c1,s1);
163 c2 = c1*c1 - s1*s1;
164 s2 = 2*(c1*s1);
165 fn = fz + n;
166 fi = fz +i;
167 gi = fz +k1-i;
168 do
169 {
170 REAL a,b,g0,f0,f1,g1,f2,g2,f3,g3;
171 b = s2*fi[k1] - c2*gi[k1];
172 a = c2*fi[k1] + s2*gi[k1];
173 f1 = fi[0 ] - a;
174 f0 = fi[0 ] + a;
175 g1 = gi[0 ] - b;
176 g0 = gi[0 ] + b;
177 b = s2*fi[k3] - c2*gi[k3];
178 a = c2*fi[k3] + s2*gi[k3];
179 f3 = fi[k2] - a;
180 f2 = fi[k2] + a;
181 g3 = gi[k2] - b;
182 g2 = gi[k2] + b;
183 b = s1*f2 - c1*g3;
184 a = c1*f2 + s1*g3;
185 fi[k2] = f0 - a;
186 fi[0 ] = f0 + a;
187 gi[k3] = g1 - b;
188 gi[k1] = g1 + b;
189 b = c1*g2 - s1*f3;
190 a = s1*g2 + c1*f3;
191 gi[k2] = g0 - a;
192 gi[0 ] = g0 + a;
193 fi[k3] = f1 - b;
194 fi[k1] = f1 + b;
195 gi += k4;
196 fi += k4;
197 } while (fi<fn);
198 }
199 TRIG_RESET(k,c1,s1);
200 } while (k4<n);
201}
202
203
204void ifft(int n, double *real, double *imag)
205{
206 double a,b,c,d;
207 double q,r,s,t;
208 int i,j,k;
209 fht(real,n);
210 fht(imag,n);
211 for (i=1,j=n-1,k=n/2;i<k;i++,j--) {
212 a = real[i]; b = real[j]; q=a+b; r=a-b;
213 c = imag[i]; d = imag[j]; s=c+d; t=c-d;
214 imag[i] = (s+r)*0.5; imag[j] = (s-r)*0.5;
215 real[i] = (q-t)*0.5; real[j] = (q+t)*0.5;
216 }
217}
218
219void realfft(int n, double *real)
220{
221 double a,b,c,d;
222 int i,j,k;
223 fht(real,n);
224 for (i=1,j=n-1,k=n/2;i<k;i++,j--) {
225 a = real[i];
226 b = real[j];
227 real[j] = (a-b)*0.5;
228 real[i] = (a+b)*0.5;
229 }
230}
231
232void fft(int n, double *real,double *imag)
233{
234 double a,b,c,d;
235 double q,r,s,t;
236 int i,j,k;
237 for (i=1,j=n-1,k=n/2;i<k;i++,j--) {
238 a = real[i]; b = real[j]; q=a+b; r=a-b;
239 c = imag[i]; d = imag[j]; s=c+d; t=c-d;
240 real[i] = (q+t)*.5; real[j] = (q-t)*.5;
241 imag[i] = (s-r)*.5; imag[j] = (s+r)*.5;
242 }
243 fht(real,n);
244 fht(imag,n);
245}
246
247void realifft(int n,double *real)
248{
249 double a,b,c,d;
250 int i,j,k;
251 for (i=1,j=n-1,k=n/2;i<k;i++,j--) {
252 a = real[i];
253 b = real[j];
254 real[j] = (a-b);
255 real[i] = (a+b);
256 }
257 fht(real,n);
258}
Note: See TracBrowser for help on using the repository browser.