1 #include "mcwfft.h"
2
3 /*--------------------------------------------------------------------
4 Complex-to-complex 2D FFT in place:
5 modex = -1 or +1
6 modey = -1 or +1
7 idx = dimension in x (power of 2)
8 idy = dimension in y (power of 2)
9 zc = input/output array
10 ----------------------------------------------------------------------*/
11
FFT_2dcx(int modex,int modey,int idx,int idy,complex * zc)12 void FFT_2dcx( int modex , int modey , int idx , int idy , complex * zc )
13 {
14 int ii , jj ;
15 static complex * tc = NULL ;
16 static int ntc_old = -1 ;
17
18 /*** make workspace, if needed ***/
19
20 if( idy > ntc_old ){
21 if( tc != NULL ) free(tc) ;
22 tc = (complex *) malloc( sizeof(complex) * idy ) ;
23 if( tc == NULL ){
24 fprintf(stderr,"\n***FFT_2dcx cannot malloc workspace!\n") ;
25 exit(-1) ;
26 }
27 ntc_old = idy ;
28 }
29
30 /*** x FFTs ***/
31
32 for( jj=0 ; jj < idy ; jj++ )
33 FFT_1dcx( modex , idx , zc + (jj*idx) ) ;
34
35 /*** y FFTs ***/
36
37 for( ii=0 ; ii < idx ; ii++ ){
38 for( jj=0 ; jj < idy ; jj++ ) tc[jj] = zc[ii+jj*idx] ;
39 FFT_1dcx( modey , idy , tc ) ;
40 for( jj=0 ; jj < idy ; jj++ ) zc[ii+jj*idx] = tc[jj] ;
41 }
42
43 return ;
44 }
45
46 /*--------------------------------------------------------------------
47 Complex-to-complex 1D FFT in place:
48 mode = -1 or +1 [phase factor in exponent: NO SCALING ON INVERSE!]
49 idim = dimension (power of 2)
50 xc = input/output array
51 Automagically re-initializes itself when idim changes from
52 previous call. By AJJesmanowicz, modified by RWCox.
53 ----------------------------------------------------------------------*/
54
FFT_1dcx(int mode,int idim,complex * xc)55 void FFT_1dcx( int mode , int idim , complex * xc )
56 {
57 register unsigned int m, n, i0, i1, i2, i3, k;
58 register complex *r0, *r1, *csp;
59 register float co, si, f0, f1, f2, f3, f4;
60 double al;
61
62 static complex * csplus = NULL , * csminus = NULL ; /* trig consts */
63 static int nold = -666 ;
64
65 /** perhaps initialize **/
66
67 if( nold != idim ){
68
69 if( idim > nold ){
70 if( csplus != NULL ){ free(csplus) ; free(csminus) ; }
71 csplus = (complex *) malloc( sizeof(complex) * idim ) ;
72 csminus = (complex *) malloc( sizeof(complex) * idim ) ;
73 if( csplus == NULL || csminus == NULL ){
74 fprintf(stderr,"\n*** FFT_1dcx cannot malloc workspace! ***\n");
75 exit(-1) ;
76 }
77 }
78 nold = n = idim ;
79
80 f1 = 1.0 ; /* csplus init */
81 m = 1;
82 k = 0;
83 while (n > m) {
84 i3 = m << 1;
85 f2 = m;
86 al = f1*PI/f2;
87 co = cos(al);
88 si = sin(al);
89 (csplus + k)->r = 1.;
90 (csplus + k)->i = 0.;
91 for (i0=0; i0 < m; i0++) {
92 k++;
93 csp = csplus + k;
94 r0 = csp - 1;
95 csp->r = r0->r * co - r0->i * si;
96 csp->i = r0->i * co + r0->r * si;
97 }
98 m = i3;
99 }
100
101 f1 = -1.0 ; /* csminus init */
102 m = 1;
103 k = 0;
104 while (n > m) {
105 i3 = m << 1;
106 f2 = m;
107 al = f1*PI/f2;
108 co = cos(al);
109 si = sin(al);
110 (csminus + k)->r = 1.;
111 (csminus + k)->i = 0.;
112 for (i0=0; i0 < m; i0++) {
113 k++;
114 csp = csminus + k;
115 r0 = csp - 1;
116 csp->r = r0->r * co - r0->i * si;
117 csp->i = r0->i * co + r0->r * si;
118 }
119 m = i3;
120 }
121 }
122
123 /** Main loop starts here **/
124
125 n = idim;
126 i2 = idim >> 1;
127 i1 = 0;
128 csp = (mode > 0) ? csplus : csminus ; /* choose const array */
129
130 for (i0=0; i0 < n; i0 ++) {
131 if ( i1 > i0 ) {
132 r0 = xc + i0;
133 r1 = xc + i1;
134 f1 = r0->r;
135 f2 = r0->i;
136 r0->r = r1->r;
137 r0->i = r1->i;
138 r1->r = f1;
139 r1->i = f2;
140 }
141 m = i2;
142 while ( m && !(i1 < m) ) {
143 i1 -= m;
144 m >>= 1;
145 }
146 i1 += m;
147 }
148
149 m = 1;
150 k = 0;
151 while (n > m) {
152 i3 = m << 1;
153 for (i0=0; i0 < m; i0 ++) {
154 co = (csp + k)->r;
155 si = (csp + k)->i;
156 for (i1=i0; i1 < n; i1 += i3) {
157 r0 = xc + i1;
158 r1 = r0 + m;
159 f1 = r1->r * co;
160 f2 = r1->i * si;
161 f3 = r1->r * si;
162 f4 = r1->i * co;
163 f1 -= f2;
164 f3 += f4;
165 f2 = r0->r;
166 f4 = r0->i;
167 r1->r = f2 - f1;
168 r1->i = f4 - f3;
169 r0->r = f2 + f1;
170 r0->i = f4 + f3;
171 }
172 k++;
173 }
174 m = i3;
175 }
176
177 #ifdef SCALE_INVERSE
178 if (mode > 0) {
179 f0 = 1.0 / idim ;
180 i0 = 0;
181 i1 = 1;
182 while (i0 < n) {
183 r0 = xc + i0;
184 r1 = xc + i1;
185 f1 = r0->r;
186 f2 = r0->i;
187 f3 = r1->r;
188 f4 = r1->i;
189 f1 *= f0;
190 f2 *= f0;
191 f3 *= f0;
192 f4 *= f0;
193 r0->r = f1;
194 r0->i = f2;
195 r1->r = f3;
196 r1->i = f4;
197 i0 += 2;
198 i1 += 2;
199 }
200 }
201 #endif
202
203 return ;
204 }
205
206 /*----------------------------------------------------------------------
207 SUM SUM x exp( -2*Pi*I * (lp+mq)*alpha - 2*Pi*I*(mp-lq)*beta )
208 l=0..N-1 lm
209 m=0..N-1
210
211 where N = idim and x = xc ; result overwrites array xc.
212 alpha and beta are arbitrary floats; N is arbitrary integer > 3.
213 ------------------------------------------------------------------------*/
214
FFT_2dchirpz(float alpha,float beta,int idim,complex * xc)215 void FFT_2dchirpz( float alpha , float beta , int idim , complex * xc )
216 {
217 int nn , ll , ll2 , pp,qq , ll2q , use_old , qloff , qnoff ;
218
219 static complex * yc = NULL , * zc = NULL , * zhatc = NULL ;
220 static int nn_old = -666 , ll_old = -666 ;
221 static float alpha_old = -987.654 , beta_old = 314.1592 ;
222
223 /*** if needed, create new workspace arrays ***/
224
225 nn = idim ;
226 if( nn == nn_old ){ /* same dimensions as before */
227 ll = ll_old ;
228 ll2 = 2*ll ;
229 ll2q = ll2 * ll2 ;
230
231 use_old = ( alpha == alpha_old && beta == beta_old ) ;
232
233 } else { /* new dimensions */
234
235 nn_old = nn ;
236 for( ll=4 ; ll < nn ; ll *= 2 ) ; /* until power of 2 >= nn */
237 ll2 = 2*ll ;
238 ll2q = ll2 * ll2 ;
239
240 use_old = 0 ;
241
242 if( ll != ll_old ){
243 ll_old = ll ;
244 if( yc != NULL ){ free(yc) ; free(zc) ; free(zhatc) ; }
245
246 yc = (complex *) malloc( sizeof(complex) * ll2q ) ;
247 zc = (complex *) malloc( sizeof(complex) * ll2q ) ;
248 zhatc = (complex *) malloc( sizeof(complex) * ll2q ) ;
249
250 if( yc == NULL || zc == NULL || zhatc == NULL ){
251 fprintf(stderr,"\n*** FFT_2dchirpz cannot malloc workspace!\n") ;
252 exit(-1) ;
253 }
254 }
255 }
256
257 /*** if needed, fill "z" workspace arrays ***/
258
259 if( ! use_old ){
260 float aq , bq , t , api ;
261
262 for( pp=0 ; pp < ll2q ; pp++ ) zc[pp].r = zc[pp].i = 0.0 ;
263
264 api = alpha * PI ;
265
266 for( qq=-(nn-1) ; qq <= (nn-1) ; qq++ ){
267
268 qloff = (qq < 0) ? ((qq+ll2)*ll2) : (qq*ll2) ;
269 aq = alpha * PI * qq * qq ;
270 bq = 2.0 * beta * PI * qq ;
271
272 for( pp=0 ; pp < 2*nn ; pp++ ){
273 t = aq - (api * pp + bq ) * pp ;
274 zc[pp+qloff].r = cos(t) ;
275 zc[pp+qloff].i = sin(t) ;
276 }
277 }
278
279 t = 1.0 / ll2q ;
280 for( pp=0 ; pp < ll2q ; pp++ ){
281 zhatc[pp].r = t * zc[pp].r ;
282 zhatc[pp].i = t * zc[pp].i ;
283 }
284 FFT_2dcx( 1 , 1 , ll2 , ll2 , zhatc ) ;
285 }
286
287 /*** fill "yc" workspace ***/
288
289 for( pp=0 ; pp < ll2q ; pp++ ) yc[pp].r = yc[pp].i = 0.0 ;
290
291 for( qq=0 ; qq < nn ; qq++ ){
292 qloff = qq*ll2 ;
293 qnoff = qq*nn ;
294 for( pp=0 ; pp < nn ; pp++ )
295 yc[pp] = CJMULT( xc[pp+qnoff] , zc[pp+qloff] ) ;
296 }
297
298 /*** 2D FFT, multiply, inverse 2D FFT ***/
299
300 FFT_2dcx( -1 , -1 , ll2 , ll2 , yc ) ;
301 for( pp=0 ; pp < ll2q ; pp++ ) yc[pp] = CMULT( yc[pp] , zhatc[pp] ) ;
302 FFT_2dcx( -1 , 1 , ll2 , ll2 , yc ) ;
303
304 /*** produce output ***/
305
306 for( qq=0 ; qq < nn ; qq++ ){
307 qloff = qq*ll2 ;
308 qnoff = qq*nn ;
309 for( pp=0 ; pp < nn ; pp++ )
310 xc[pp+qnoff] = CMULT( zc[pp+qloff] , yc[pp+qloff] ) ;
311 }
312
313 return ;
314 }
315