?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
bump version to 3.5beta1
[supercollider.git] / common / fftlib.c
blobcd2220288b2e4f22127199b783be0932bd56420d
1 /* FFT library */
2 /* Library of in-place fast fourier transforms */
3 /* Forward and inverse complex transforms */
4 /* and real forward transform */
5 /* Optimized for RISC processors with many registers */
6 /* Version 1.1 John Green NUWC New London CT January 96 */
7 /* Version 1.1 renamed as fftlib from fftbig */
8 /* Version 1.1 removed (float *)((char *) ptr) optimization */
9 /* Version 1.0 John Green NUWC New London CT December 95 */
10 /* (John Green) green_jt@vsdec.nl.nuwc.navy.mil */
11 /* green_jt@vsdec.nl.nuwc.navy.mil */
12
13 #include <math.h>
14 #include "fftlib.h"
15
16 #define MAXMROOT 9 /* 2^(MAXMROOT-1) = # of elements in BRcnt */
17
18 /* Bit reversed counter */
19 static const unsigned char BRcnt[256] = {
20 0, 128, 64, 192, 32, 160, 96, 224, 16, 144, 80, 208,
21 48, 176, 112, 240, 8, 136, 72, 200, 40, 168, 104, 232,
22 24, 152, 88, 216, 56, 184, 120, 248, 4, 132, 68, 196,
23 36, 164, 100, 228, 20, 148, 84, 212, 52, 180, 116, 244,
24 12, 140, 76, 204, 44, 172, 108, 236, 28, 156, 92, 220,
25 60, 188, 124, 252, 2, 130, 66, 194, 34, 162, 98, 226,
26 18, 146, 82, 210, 50, 178, 114, 242, 10, 138, 74, 202,
27 42, 170, 106, 234, 26, 154, 90, 218, 58, 186, 122, 250,
28 6, 134, 70, 198, 38, 166, 102, 230, 22, 150, 86, 214,
29 54, 182, 118, 246, 14, 142, 78, 206, 46, 174, 110, 238,
30 30, 158, 94, 222, 62, 190, 126, 254, 1, 129, 65, 193,
31 33, 161, 97, 225, 17, 145, 81, 209, 49, 177, 113, 241,
32 9, 137, 73, 201, 41, 169, 105, 233, 25, 153, 89, 217,
33 57, 185, 121, 249, 5, 133, 69, 197, 37, 165, 101, 229,
34 21, 149, 85, 213, 53, 181, 117, 245, 13, 141, 77, 205,
35 45, 173, 109, 237, 29, 157, 93, 221, 61, 189, 125, 253,
36 3, 131, 67, 195, 35, 163, 99, 227, 19, 147, 83, 211,
37 51, 179, 115, 243, 11, 139, 75, 203, 43, 171, 107, 235,
38 27, 155, 91, 219, 59, 187, 123, 251, 7, 135, 71, 199,
39 39, 167, 103, 231, 23, 151, 87, 215, 55, 183, 119, 247,
40 15, 143, 79, 207, 47, 175, 111, 239, 31, 159, 95, 223,
41 63, 191, 127, 255 };
42
43 /* Table of powers of 2 */
44 static const long Ntbl[21] = {1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048,
45 4096, 8192, 16384, 32768, 65536, 131072, 262144, 524288, 1048576};
46
47 long FFTInit(long *fftMptr, long fftN, float *Utbl){
48 /* Compute cosine table and check size for complex ffts */
49 /* INPUTS */
50 /* fftN = size of fft */
51 /* OUTPUTS */
52 /* *fftMptr = log2 of fft size */
53 /* *Utbl = cosine table with fftN/4 + 1 entries (angles = 0 to pi/2 inclusive) */
54 /* RETURNS */
55 /* 1 if fftN is invalid, 0 otherwise */
56 long i1, ErrVal;
57 ErrVal = 0;
58 *fftMptr = (long)(log(fftN)/log(2.0) + 0.5);
59 if ((*fftMptr >= 3) & (*fftMptr <= 19) & (int)(pow(2,*fftMptr)+.5) == fftN)
60 for (i1 = 0; i1 <= fftN/4; i1++)
61 Utbl[i1] = cos( ( 3.1415926535897932384626433832795 * 2.0 * i1 ) / fftN );
62 else
63 ErrVal = 1;
64 return ErrVal;
65 }
66
67 long rFFTInit(long *fftMptr, long fftN, float *Utbl){
68 /* Compute cosine table and check size for a real input fft */
69 /* INPUTS */
70 /* fftN = size of fft */
71 /* OUTPUTS */
72 /* *fftMptr = log2 of fft size */
73 /* *Utbl = cosine table with fftN/4 + 1 entries (angles = 0 to pi/2 inclusive) */
74 /* RETURNS */
75 /* 1 if fftN is invalid, 0 otherwise */
76 long i1, ErrVal;
77 ErrVal = 0;
78 *fftMptr = (long)(log(fftN)/log(2.0) + 0.5);
79 if ((*fftMptr >= 4) & (*fftMptr <= 20) & (int)(pow(2,*fftMptr)+.5) == fftN)
80
81 for (i1 = 0; i1 <= fftN/4; i1++)
82 Utbl[i1] = cos( ( 3.1415926535897932384626433832795 * 2.0 * i1 ) / fftN );
83 else
84 ErrVal = 1;
85 return ErrVal;
86 }
87
88 void ffts(float *ioptr, long M, long Rows, float *Utbl){
89 /* Compute in-place complex fft on the rows of the input array */
90 /* INPUTS */
91 /* M = log2 of fft size */
92 /* *ioptr = input data array */
93 /* *Utbl = cosine table */
94 /* Rows = number of rows in ioptr array (use Rows of 1 if ioptr is a 1 dimensional array) */
95 /* OUTPUTS */
96 /* *ioptr = output data array */
97
98 long Flyinc;
99 long FlyOffsetA;
100 long FlyOffsetAIm;
101 long FlyOffsetB;
102 long FlyOffsetBIm;
103 long NsameU1;
104 long NsameU2;
105 long NsameU4;
106 long diffUcnt;
107 long LoopCnt;
108
109 float scale;
110 float fly0r;
111 float fly0i;
112 float fly1r;
113 float fly1i;
114 float fly2r;
115 float fly2i;
116 float fly3r;
117 float fly3i;
118 float fly4r;
119 float fly4i;
120 float fly5r;
121 float fly5i;
122 float fly6r;
123 float fly6i;
124 float fly7r;
125 float fly7i;
126 float U0r;
127 float U0i;
128 float U1r;
129 float U1i;
130 float U2r;
131 float U2i;
132 float U3r;
133 float U3i;
134 float t0r;
135 float t0i;
136 float t1r;
137 float t1i;
138
139 float *fly0P;
140 float *fly1P;
141 float *fly2P;
142 float *fly3P;
143
144 float *U0rP;
145 float *U0iP;
146 float *U1rP;
147 float *U1iP;
148 float *U2rP;
149 float *U2iP;
150 float *IOP;
151 long U3offset;
152
153 long stage;
154 long NdiffU;
155 long LoopN;
156
157 const long BRshift = MAXMROOT - (M>>1);
158 const long Nrems2 = Ntbl[M-(M>>1)+1];
159 const long Nroot_1_ColInc = (Ntbl[M-1]-Ntbl[M-(M>>1)])*2;
160
161 scale = 2.0;
162 for (;Rows>0;Rows--){
163
164 FlyOffsetA = Ntbl[M] * (2/2);
165 FlyOffsetAIm = FlyOffsetA + 1;
166 FlyOffsetB = FlyOffsetA + 2;
167 FlyOffsetBIm = FlyOffsetB + 1;
168
169 /* BitrevR2 ** bit reverse and first radix 2 stage ******/
170
171 for (stage = 0; stage < Ntbl[M-(M>>1)]*2; stage += Ntbl[M>>1]*2){
172 for (LoopN = (Ntbl[(M>>1)-1]-1); LoopN >= 0; LoopN--){
173 LoopCnt = (Ntbl[(M>>1)-1]-1);
174 fly0P = ioptr+ Nroot_1_ColInc + ((int)BRcnt[LoopN] >> BRshift)*(2*2) + stage;
175 IOP = ioptr + (LoopN<<(M+1)/2) * 2 + stage;
176 fly1P = IOP + ((int)BRcnt[LoopCnt] >> BRshift)*(2*2);
177 fly0r = *(fly0P);
178 fly0i = *(fly0P+1);
179 fly1r = *(fly0P+FlyOffsetA);
180 fly1i = *(fly0P+FlyOffsetAIm);
181 for (; LoopCnt > LoopN;){
182 fly2r = *(fly0P+2);
183 fly2i = *(fly0P+(2+1));
184 fly3r = *(fly0P+FlyOffsetB);
185 fly3i = *(fly0P+FlyOffsetBIm);
186 fly4r = *(fly1P);
187 fly4i = *(fly1P+1);
188 fly5r = *(fly1P+FlyOffsetA);
189 fly5i = *(fly1P+FlyOffsetAIm);
190 fly6r = *(fly1P+2);
191 fly6i = *(fly1P+(2+1));
192 fly7r = *(fly1P+FlyOffsetB);
193 fly7i = *(fly1P+FlyOffsetBIm);
194
195 t0r = fly0r + fly1r;
196 t0i = fly0i + fly1i;
197 fly1r = fly0r - fly1r;
198 fly1i = fly0i - fly1i;
199 t1r = fly2r + fly3r;
200 t1i = fly2i + fly3i;
201 fly3r = fly2r - fly3r;
202 fly3i = fly2i - fly3i;
203 fly0r = fly4r + fly5r;
204 fly0i = fly4i + fly5i;
205 fly5r = fly4r - fly5r;
206 fly5i = fly4i - fly5i;
207 fly2r = fly6r + fly7r;
208 fly2i = fly6i + fly7i;
209 fly7r = fly6r - fly7r;
210 fly7i = fly6i - fly7i;
211
212 *(fly1P) = t0r;
213 *(fly1P+1) = t0i;
214 *(fly1P+2) = fly1r;
215 *(fly1P+(2+1)) = fly1i;
216 *(fly1P+FlyOffsetA) = t1r;
217 *(fly1P+FlyOffsetAIm) = t1i;
218 *(fly1P+FlyOffsetB) = fly3r;
219 *(fly1P+FlyOffsetBIm) = fly3i;
220 *(fly0P) = fly0r;
221 *(fly0P+1) = fly0i;
222 *(fly0P+2) = fly5r;
223 *(fly0P+(2+1)) = fly5i;
224 *(fly0P+FlyOffsetA) = fly2r;
225 *(fly0P+FlyOffsetAIm) = fly2i;
226 *(fly0P+FlyOffsetB) = fly7r;
227 *(fly0P+FlyOffsetBIm) = fly7i;
228
229 fly0P -= Nrems2;
230 fly0r = *(fly0P);
231 fly0i = *(fly0P+1);
232 fly1r = *(fly0P+FlyOffsetA);
233 fly1i = *(fly0P+FlyOffsetAIm);
234 LoopCnt -= 1;
235 fly1P = (IOP + ((int)BRcnt[LoopCnt] >> BRshift)*(2*2));
236 };
237 fly2r = *(fly0P+2);
238 fly2i = *(fly0P+(2+1));
239 fly3r = *(fly0P+FlyOffsetB);
240 fly3i = *(fly0P+FlyOffsetBIm);
241
242 t0r = fly0r + fly1r;
243 t0i = fly0i + fly1i;
244 fly1r = fly0r - fly1r;
245 fly1i = fly0i - fly1i;
246 t1r = fly2r + fly3r;
247 t1i = fly2i + fly3i;
248 fly3r = fly2r - fly3r;
249 fly3i = fly2i - fly3i;
250
251 *(fly0P) = t0r;
252 *(fly0P+1) = t0i;
253 *(fly0P+2) = fly1r;
254 *(fly0P+(2+1)) = fly1i;
255 *(fly0P+FlyOffsetA) = t1r;
256 *(fly0P+FlyOffsetAIm) = t1i;
257 *(fly0P+FlyOffsetB) = fly3r;
258 *(fly0P+FlyOffsetBIm) = fly3i;
259
260 };
261 };
262
263
264
265 /**** FFTC **************/
266
267 NdiffU = 2;
268 Flyinc = (NdiffU);
269
270 NsameU4 = Ntbl[M]/4;
271 LoopN = Ntbl[M-3];
272
273 stage = ((M-1)/3);
274 if ( (M-1-(stage * 3)) != 0 ){
275 FlyOffsetA = Flyinc << 2;
276 FlyOffsetB = FlyOffsetA << 1;
277 FlyOffsetAIm = FlyOffsetA + 1;
278 FlyOffsetBIm = FlyOffsetB + 1;
279 if ( (M-1-(stage * 3)) == 1 ){
280 /* 1 radix 2 stage */
281
282 IOP = ioptr;
283 fly0P = IOP;
284 fly1P = (IOP+Flyinc);
285 fly2P = (fly1P+Flyinc);
286 fly3P = (fly2P+Flyinc);
287
288 /* Butterflys */
289 /*
290 t0 - - t0
291 t1 - - t1
292 f2 - 1- f5
293 f1 - -i- f7
294 f4 - - f4
295 f0 - - f0
296 f6 - 1- f2
297 f3 - -i- f1
298 */
299
300 for (LoopCnt = LoopN; LoopCnt > 0 ; LoopCnt--){
301 t0r = *(fly0P);
302 t0i = *(fly0P + 1);
303 t1r = *(fly1P);
304 t1i = *(fly1P + 1);
305 fly2r = *(fly2P);
306 fly2i = *(fly2P + 1);
307 fly1r = *(fly3P);
308 fly1i = *(fly3P + 1);
309 fly4r = *(fly0P + FlyOffsetA);
310 fly4i = *(fly0P + FlyOffsetAIm);
311 fly0r = *(fly1P + FlyOffsetA);
312 fly0i = *(fly1P + FlyOffsetAIm);
313 fly6r = *(fly2P + FlyOffsetA);
314 fly6i = *(fly2P + FlyOffsetAIm);
315 fly3r = *(fly3P + FlyOffsetA);
316 fly3i = *(fly3P + FlyOffsetAIm);
317
318 fly5r = t0r - fly2r;
319 fly5i = t0i - fly2i;
320 t0r = t0r + fly2r;
321 t0i = t0i + fly2i;
322
323 fly7r = t1r - fly1i;
324 fly7i = t1i + fly1r;
325 t1r = t1r + fly1i;
326 t1i = t1i - fly1r;
327
328 fly2r = fly4r - fly6r;
329 fly2i = fly4i - fly6i;
330 fly4r = fly4r + fly6r;
331 fly4i = fly4i + fly6i;
332
333 fly1r = fly0r - fly3i;
334 fly1i = fly0i + fly3r;
335 fly0r = fly0r + fly3i;
336 fly0i = fly0i - fly3r;
337
338 *(fly2P) = fly5r;
339 *(fly2P + 1) = fly5i;
340 *(fly0P) = t0r;
341 *(fly0P + 1) = t0i;
342 *(fly3P) = fly7r;
343 *(fly3P + 1) = fly7i;
344 *(fly1P) = t1r;
345 *(fly1P + 1) = t1i;
346 *(fly2P + FlyOffsetA) = fly2r;
347 *(fly2P + FlyOffsetAIm) = fly2i;
348 *(fly0P + FlyOffsetA) = fly4r;
349 *(fly0P + FlyOffsetAIm) = fly4i;
350 *(fly3P + FlyOffsetA) = fly1r;
351 *(fly3P + FlyOffsetAIm) = fly1i;
352 *(fly1P + FlyOffsetA) = fly0r;
353 *(fly1P + FlyOffsetAIm) = fly0i;
354
355 fly0P = (fly0P + FlyOffsetB);
356 fly1P = (fly1P + FlyOffsetB);
357 fly2P = (fly2P + FlyOffsetB);
358 fly3P = (fly3P + FlyOffsetB);
359 };
360
361 NsameU4 >>= 1;
362 LoopN >>= 1;
363 NdiffU <<= 1;
364 Flyinc = Flyinc << 1;
365 }
366 else{
367 /* 1 radix 4 stage */
368 IOP = ioptr;
369
370 U3r = 0.7071067811865475244008443621; /* sqrt(0.5); */
371 U3i = U3r;
372 fly0P = IOP;
373 fly1P = (IOP+Flyinc);
374 fly2P = (fly1P+Flyinc);
375 fly3P = (fly2P+Flyinc);
376
377 /* Butterflys */
378 /*
379 t0 - - t0 - - t0
380 t1 - - t1 - - t1
381 f2 - 1- f5 - - f5
382 f1 - -i- f7 - - f7
383 f4 - - f4 - 1- f6
384 f0 - - f0 -U3 - f3
385 f6 - 1- f2 - -i- f4
386 f3 - -i- f1 -U3a- f2
387 */
388
389 for (LoopCnt = LoopN; LoopCnt > 0 ; LoopCnt--){
390 t0r = *(fly0P);
391 t0i = *(fly0P + 1);
392 t1r = *(fly1P);
393 t1i = *(fly1P + 1);
394 fly2r = *(fly2P);
395 fly2i = *(fly2P + 1);
396 fly1r = *(fly3P);
397 fly1i = *(fly3P + 1);
398 fly4r = *(fly0P + FlyOffsetA);
399 fly4i = *(fly0P + FlyOffsetAIm);
400 fly0r = *(fly1P + FlyOffsetA);
401 fly0i = *(fly1P + FlyOffsetAIm);
402 fly6r = *(fly2P + FlyOffsetA);
403 fly6i = *(fly2P + FlyOffsetAIm);
404 fly3r = *(fly3P + FlyOffsetA);
405 fly3i = *(fly3P + FlyOffsetAIm);
406
407 fly5r = t0r - fly2r;
408 fly5i = t0i - fly2i;
409 t0r = t0r + fly2r;
410 t0i = t0i + fly2i;
411
412 fly7r = t1r - fly1i;
413 fly7i = t1i + fly1r;
414 t1r = t1r + fly1i;
415 t1i = t1i - fly1r;
416
417 fly2r = fly4r - fly6r;
418 fly2i = fly4i - fly6i;
419 fly4r = fly4r + fly6r;
420 fly4i = fly4i + fly6i;
421
422 fly1r = fly0r - fly3i;
423 fly1i = fly0i + fly3r;
424 fly0r = fly0r + fly3i;
425 fly0i = fly0i - fly3r;
426
427
428 fly6r = t0r - fly4r;
429 fly6i = t0i - fly4i;
430 t0r = t0r + fly4r;
431 t0i = t0i + fly4i;
432
433 fly3r = fly5r - fly2i;
434 fly3i = fly5i + fly2r;
435 fly5r = fly5r + fly2i;
436 fly5i = fly5i - fly2r;
437
438 fly4r = t1r - U3r * fly0r;
439 fly4r = fly4r - U3i * fly0i;
440 fly4i = t1i + U3i * fly0r;
441 fly4i = fly4i - U3r * fly0i;
442 t1r = scale * t1r - fly4r;
443 t1i = scale * t1i - fly4i;
444
445 fly2r = fly7r + U3i * fly1r;
446 fly2r = fly2r - U3r * fly1i;
447 fly2i = fly7i + U3r * fly1r;
448 fly2i = fly2i + U3i * fly1i;
449 fly7r = scale * fly7r - fly2r;
450 fly7i = scale * fly7i - fly2i;
451
452 *(fly0P + FlyOffsetA) = fly6r;
453 *(fly0P + FlyOffsetAIm) = fly6i;
454 *(fly0P) = t0r;
455 *(fly0P + 1) = t0i;
456 *(fly2P + FlyOffsetA) = fly3r;
457 *(fly2P + FlyOffsetAIm) = fly3i;
458 *(fly2P) = fly5r;
459 *(fly2P + 1) = fly5i;
460 *(fly1P + FlyOffsetA) = fly4r;
461 *(fly1P + FlyOffsetAIm) = fly4i;
462 *(fly1P) = t1r;
463 *(fly1P + 1) = t1i;
464 *(fly3P + FlyOffsetA) = fly2r;
465 *(fly3P + FlyOffsetAIm) = fly2i;
466 *(fly3P) = fly7r;
467 *(fly3P + 1) = fly7i;
468
469 fly0P = (fly0P + FlyOffsetB);
470 fly1P = (fly1P + FlyOffsetB);
471 fly2P = (fly2P + FlyOffsetB);
472 fly3P = (fly3P + FlyOffsetB);
473
474 };
475
476 NsameU4 >>= 2;
477 LoopN >>= 2;
478 NdiffU <<= 2;
479 Flyinc = Flyinc << 2;
480 };
481 };
482
483 NsameU2 = NsameU4 >> 1;
484 NsameU1 = NsameU2 >> 1;
485 Flyinc <<= 1;
486 FlyOffsetA = Flyinc << 2;
487 FlyOffsetB = FlyOffsetA << 1;
488 FlyOffsetAIm = FlyOffsetA + 1;
489 FlyOffsetBIm = FlyOffsetB + 1;
490 LoopN >>= 1;
491 /* ****** RADIX 8 Stages */
492 for (stage = stage<<1; stage > 0 ; stage--){
493
494 /* an fft stage is done in two parts to ease use of the single quadrant cos table */
495
496 /* fftcalc1(iobuf, Utbl, N, NdiffU, LoopN); */
497 if(!(stage&1)){
498 U0rP = &Utbl[0];
499 U0iP = &Utbl[Ntbl[M-2]];
500 U1rP = U0rP;
501 U1iP = U0iP;
502 U2rP = U0rP;
503 U2iP = U0iP;
504 U3offset = (Ntbl[M]) / 8;
505
506 IOP = ioptr;
507
508 U0r = *U0rP,
509 U0i = *U0iP;
510 U1r = *U1rP,
511 U1i = *U1iP;
512 U2r = *U2rP,
513 U2i = *U2iP;
514 U3r = *( U2rP + U3offset);
515 U3i = *( U2iP - U3offset);
516 }
517
518 fly0P = IOP;
519 fly1P = (IOP+Flyinc);
520 fly2P = (fly1P+Flyinc);
521 fly3P = (fly2P+Flyinc);
522
523 for (diffUcnt = (NdiffU)>>1; diffUcnt != 0; diffUcnt--){
524
525 /* Butterflys */
526 /*
527 f0 - - t0 - - t0 - - t0
528 f1 -U0 - t1 - - t1 - - t1
529 f2 - - f2 -U1 - f5 - - f3
530 f3 -U0 - f1 -U1a- f7 - - f7
531 f4 - - f4 - - f4 -U2 - f6
532 f5 -U0 - f0 - - f0 -U3 - f4
533 f6 - - f6 -U1 - f2 -U2a- f2
534 f7 -U0 - f3 -U1a- f1 -U3a- f5
535 */
536
537 fly0r = *(IOP);
538 fly0i = *(IOP+1);
539 fly1r = *(fly1P);
540 fly1i = *(fly1P+1);
541
542 for (LoopCnt = LoopN-1; LoopCnt > 0 ; LoopCnt--){
543
544 fly2r = *(fly2P);
545 fly2i = *(fly2P + 1);
546 fly3r = *(fly3P);
547 fly3i = *(fly3P + 1);
548 fly4r = *(fly0P + FlyOffsetA);
549 fly4i = *(fly0P + FlyOffsetAIm);
550 fly5r = *(fly1P + FlyOffsetA);
551 fly5i = *(fly1P + FlyOffsetAIm);
552 fly6r = *(fly2P + FlyOffsetA);
553 fly6i = *(fly2P + FlyOffsetAIm);
554 fly7r = *(fly3P + FlyOffsetA);
555 fly7i = *(fly3P + FlyOffsetAIm);
556
557 t1r = fly0r - U0r * fly1r;
558 t1r = t1r - U0i * fly1i;
559 t1i = fly0i + U0i * fly1r;
560 t1i = t1i - U0r * fly1i;
561 t0r = scale * fly0r - t1r;
562 t0i = scale * fly0i - t1i;
563
564 fly1r = fly2r - U0r * fly3r;
565 fly1r = fly1r - U0i * fly3i;
566 fly1i = fly2i + U0i * fly3r;
567 fly1i = fly1i - U0r * fly3i;
568 fly2r = scale * fly2r - fly1r;
569 fly2i = scale * fly2i - fly1i;
570
571 fly0r = fly4r - U0r * fly5r;
572 fly0r = fly0r - U0i * fly5i;
573 fly0i = fly4i + U0i * fly5r;
574 fly0i = fly0i - U0r * fly5i;
575 fly4r = scale * fly4r - fly0r;
576 fly4i = scale * fly4i - fly0i;
577
578 fly3r = fly6r - U0r * fly7r;
579 fly3r = fly3r - U0i * fly7i;
580 fly3i = fly6i + U0i * fly7r;
581 fly3i = fly3i - U0r * fly7i;
582 fly6r = scale * fly6r - fly3r;
583 fly6i = scale * fly6i - fly3i;
584
585
586 fly5r = t0r - U1r * fly2r;
587 fly5r = fly5r - U1i * fly2i;
588 fly5i = t0i + U1i * fly2r;
589 fly5i = fly5i - U1r * fly2i;
590 t0r = scale * t0r - fly5r;
591 t0i = scale * t0i - fly5i;
592
593 fly7r = t1r + U1i * fly1r;
594 fly7r = fly7r - U1r * fly1i;
595 fly7i = t1i + U1r * fly1r;
596 fly7i = fly7i + U1i * fly1i;
597 t1r = scale * t1r - fly7r;
598 t1i = scale * t1i - fly7i;
599
600 fly2r = fly4r - U1r * fly6r;
601 fly2r = fly2r - U1i * fly6i;
602 fly2i = fly4i + U1i * fly6r;
603 fly2i = fly2i - U1r * fly6i;
604 fly4r = scale * fly4r - fly2r;
605 fly4i = scale * fly4i - fly2i;
606
607 fly1r = fly0r + U1i * fly3r;
608 fly1r = fly1r - U1r * fly3i;
609 fly1i = fly0i + U1r * fly3r;
610 fly1i = fly1i + U1i * fly3i;
611 fly0r = scale * fly0r - fly1r;
612 fly0i = scale * fly0i - fly1i;
613
614 fly6r = t0r - U2r * fly4r;
615 fly6r = fly6r - U2i * fly4i;
616 fly6i = t0i + U2i * fly4r;
617 fly6i = fly6i - U2r * fly4i;
618 t0r = scale * t0r - fly6r;
619 t0i = scale * t0i - fly6i;
620
621 fly3r = fly5r - U2i * fly2r;
622 fly3r = fly3r + U2r * fly2i;
623 fly3i = fly5i - U2r * fly2r;
624 fly3i = fly3i - U2i * fly2i;
625 fly2r = scale * fly5r - fly3r;
626 fly2i = scale * fly5i - fly3i;
627
628 fly4r = t1r - U3r * fly0r;
629 fly4r = fly4r - U3i * fly0i;
630 fly4i = t1i + U3i * fly0r;
631 fly4i = fly4i - U3r * fly0i;
632 t1r = scale * t1r - fly4r;
633 t1i = scale * t1i - fly4i;
634
635 fly5r = fly7r + U3i * fly1r;
636 fly5r = fly5r - U3r * fly1i;
637 fly5i = fly7i + U3r * fly1r;
638 fly5i = fly5i + U3i * fly1i;
639 fly7r = scale * fly7r - fly5r;
640 fly7i = scale * fly7i - fly5i;
641
642 *(fly0P + FlyOffsetA) = fly6r;
643 *(fly0P + FlyOffsetAIm) = fly6i;
644 *(fly0P) = t0r;
645 *(fly0P + 1) = t0i;
646 *(fly2P) = fly3r;
647 *(fly2P + 1) = fly3i;
648 *(fly2P + FlyOffsetA) = fly2r;
649 *(fly2P + FlyOffsetAIm) = fly2i;
650
651 fly0r = *(fly0P + FlyOffsetB);
652 fly0i = *(fly0P + FlyOffsetBIm);
653
654 *(fly1P + FlyOffsetA) = fly4r;
655 *(fly1P + FlyOffsetAIm) = fly4i;
656 *(fly1P) = t1r;
657 *(fly1P + 1) = t1i;
658
659 fly1r = *(fly1P + FlyOffsetB);
660 fly1i = *(fly1P + FlyOffsetBIm);
661
662 *(fly3P + FlyOffsetA) = fly5r;
663 *(fly3P + FlyOffsetAIm) = fly5i;
664 *(fly3P) = fly7r;
665 *(fly3P + 1) = fly7i;
666
667 fly0P = (fly0P + FlyOffsetB);
668 fly1P = (fly1P + FlyOffsetB);
669 fly2P = (fly2P + FlyOffsetB);
670 fly3P = (fly3P + FlyOffsetB);
671 };
672 fly2r = *(fly2P);
673 fly2i = *(fly2P + 1);
674 fly3r = *(fly3P);
675 fly3i = *(fly3P + 1);
676 fly4r = *(fly0P + FlyOffsetA);
677 fly4i = *(fly0P + FlyOffsetAIm);
678 fly5r = *(fly1P + FlyOffsetA);
679 fly5i = *(fly1P + FlyOffsetAIm);
680 fly6r = *(fly2P + FlyOffsetA);
681 fly6i = *(fly2P + FlyOffsetAIm);
682 fly7r = *(fly3P + FlyOffsetA);
683 fly7i = *(fly3P + FlyOffsetAIm);
684
685 t1r = fly0r - U0r * fly1r;
686 t1r = t1r - U0i * fly1i;
687 t1i = fly0i + U0i * fly1r;
688 t1i = t1i - U0r * fly1i;
689 t0r = scale * fly0r - t1r;
690 t0i = scale * fly0i - t1i;
691
692 fly1r = fly2r - U0r * fly3r;
693 fly1r = fly1r - U0i * fly3i;
694 fly1i = fly2i + U0i * fly3r;
695 fly1i = fly1i - U0r * fly3i;
696 fly2r = scale * fly2r - fly1r;
697 fly2i = scale * fly2i - fly1i;
698
699 fly0r = fly4r - U0r * fly5r;
700 fly0r = fly0r - U0i * fly5i;
701 fly0i = fly4i + U0i * fly5r;
702 fly0i = fly0i - U0r * fly5i;
703 fly4r = scale * fly4r - fly0r;
704 fly4i = scale * fly4i - fly0i;
705
706 fly3r = fly6r - U0r * fly7r;
707 fly3r = fly3r - U0i * fly7i;
708 fly3i = fly6i + U0i * fly7r;
709 fly3i = fly3i - U0r * fly7i;
710 fly6r = scale * fly6r - fly3r;
711 fly6i = scale * fly6i - fly3i;
712
713 fly5r = t0r - U1r * fly2r;
714 fly5r = fly5r - U1i * fly2i;
715 fly5i = t0i + U1i * fly2r;
716 fly5i = fly5i - U1r * fly2i;
717 t0r = scale * t0r - fly5r;
718 t0i = scale * t0i - fly5i;
719
720 fly7r = t1r + U1i * fly1r;
721 fly7r = fly7r - U1r * fly1i;
722 fly7i = t1i + U1r * fly1r;
723 fly7i = fly7i + U1i * fly1i;
724 t1r = scale * t1r - fly7r;
725 t1i = scale * t1i - fly7i;
726
727 fly2r = fly4r - U1r * fly6r;
728 fly2r = fly2r - U1i * fly6i;
729 fly2i = fly4i + U1i * fly6r;
730 fly2i = fly2i - U1r * fly6i;
731 fly4r = scale * fly4r - fly2r;
732 fly4i = scale * fly4i - fly2i;
733
734 fly1r = fly0r + U1i * fly3r;
735 fly1r = fly1r - U1r * fly3i;
736 fly1i = fly0i + U1r * fly3r;
737 fly1i = fly1i + U1i * fly3i;
738 fly0r = scale * fly0r - fly1r;
739 fly0i = scale * fly0i - fly1i;
740
741 U0i = *(U0iP = (U0iP - NsameU4));
742 U0r = *(U0rP = (U0rP + NsameU4));
743 U1r = *(U1rP = (U1rP + NsameU2));
744 U1i = *(U1iP = (U1iP - NsameU2));
745 if(stage&1)
746 U0r = -U0r;
747
748 fly6r = t0r - U2r * fly4r;
749 fly6r = fly6r - U2i * fly4i;
750 fly6i = t0i + U2i * fly4r;
751 fly6i = fly6i - U2r * fly4i;
752 t0r = scale * t0r - fly6r;
753 t0i = scale * t0i - fly6i;
754
755 fly3r = fly5r - U2i * fly2r;
756 fly3r = fly3r + U2r * fly2i;
757 fly3i = fly5i - U2r * fly2r;
758 fly3i = fly3i - U2i * fly2i;
759 fly2r = scale * fly5r - fly3r;
760 fly2i = scale * fly5i - fly3i;
761
762 fly4r = t1r - U3r * fly0r;
763 fly4r = fly4r - U3i * fly0i;
764 fly4i = t1i + U3i * fly0r;
765 fly4i = fly4i - U3r * fly0i;
766 t1r = scale * t1r - fly4r;
767 t1i = scale * t1i - fly4i;
768
769 fly5r = fly7r + U3i * fly1r;
770 fly5r = fly5r - U3r * fly1i;
771 fly5i = fly7i + U3r * fly1r;
772 fly5i = fly5i + U3i * fly1i;
773 fly7r = scale * fly7r - fly5r;
774 fly7i = scale * fly7i - fly5i;
775
776 *(fly0P + FlyOffsetA) = fly6r;
777 *(fly0P + FlyOffsetAIm) = fly6i;
778 *(fly0P) = t0r;
779 *(fly0P + 1) = t0i;
780
781 U2r = *(U2rP = (U2rP + NsameU1));
782 U2i = *(U2iP = (U2iP - NsameU1));
783
784 *(fly2P) = fly3r;
785 *(fly2P + 1) = fly3i;
786 *(fly2P + FlyOffsetA) = fly2r;
787 *(fly2P + FlyOffsetAIm) = fly2i;
788 *(fly1P + FlyOffsetA) = fly4r;
789 *(fly1P + FlyOffsetAIm) = fly4i;
790 *(fly1P) = t1r;
791 *(fly1P + 1) = t1i;
792
793 U3r = *( U2rP + U3offset);
794 U3i = *( U2iP - U3offset);
795
796 *(fly3P + FlyOffsetA) = fly5r;
797 *(fly3P + FlyOffsetAIm) = fly5i;
798 *(fly3P) = fly7r;
799 *(fly3P + 1) = fly7i;
800
801 IOP = IOP + 2;
802 fly0P = IOP;
803 fly1P = (IOP+Flyinc);
804 fly2P = (fly1P+Flyinc);
805 fly3P = (fly2P+Flyinc);
806 };
807 NsameU4 = -NsameU4;
808
809 if(stage&1){
810 LoopN >>= 3;
811 NsameU1 >>= 3;
812 NsameU2 >>= 3;
813 NsameU4 >>= 3;
814 NdiffU <<= 3;
815 Flyinc <<= 3;
816 FlyOffsetA <<= 3;
817 FlyOffsetB <<= 3;
818 FlyOffsetAIm = FlyOffsetA + 1;
819 FlyOffsetBIm = FlyOffsetB + 1;
820 }
821 }
822 ioptr += 2*Ntbl[M];
823 }
824 }
825
826 void iffts(float *ioptr, long M, long Rows, float *Utbl){
827 /* Compute in-place inverse complex fft on the rows of the input array */
828 /* INPUTS */
829 /* M = log2 of fft size */
830 /* *ioptr = input data array */
831 /* *Utbl = cosine table */
832 /* Rows = number of rows in ioptr array (use Rows of 1 if ioptr is a 1 dimensional array) */
833 /* OUTPUTS */
834 /* *ioptr = output data array */
835
836 long Flyinc;
837 long FlyOffsetA;
838 long FlyOffsetAIm;
839 long FlyOffsetB;
840 long FlyOffsetBIm;
841 long NsameU1;
842 long NsameU2;
843 long NsameU4;
844 long diffUcnt;
845 long LoopCnt;
846
847 float scale;
848 float fly0r;
849 float fly0i;
850 float fly1r;
851 float fly1i;
852 float fly2r;
853 float fly2i;
854 float fly3r;
855 float fly3i;
856 float fly4r;
857 float fly4i;
858 float fly5r;
859 float fly5i;
860 float fly6r;
861 float fly6i;
862 float fly7r;
863 float fly7i;
864 float U0r;
865 float U0i;
866 float U1r;
867 float U1i;
868 float U2r;
869 float U2i;
870 float U3r;
871 float U3i;
872 float t0r;
873 float t0i;
874 float t1r;
875 float t1i;
876
877 float *fly0P;
878 float *fly1P;
879 float *fly2P;
880 float *fly3P;
881 float *U0rP;
882 float *U0iP;
883 float *U1rP;
884 float *U1iP;
885 float *U2rP;
886 float *U2iP;
887 float *IOP;
888 long U3offset;
889
890 long stage;
891 long NdiffU;
892 long LoopN;
893
894 const long BRshift = MAXMROOT - (M>>1);
895 const long Nrems2 = Ntbl[M-(M>>1)+1];
896 const long Nroot_1_ColInc = (Ntbl[M-1]-Ntbl[M-(M>>1)])*2;
897
898 for (;Rows>0;Rows--){
899
900 FlyOffsetA = Ntbl[M] * 2/2;
901 FlyOffsetAIm = FlyOffsetA + 1;
902 FlyOffsetB = FlyOffsetA + 2;
903 FlyOffsetBIm = FlyOffsetB + 1;
904
905 /* BitrevR2 ** bit reverse and first radix 2 stage ******/
906
907 scale = 1./Ntbl[M];
908 for (stage = 0; stage < Ntbl[M-(M>>1)]*2; stage += Ntbl[M>>1]*2){
909 for (LoopN = (Ntbl[(M>>1)-1]-1); LoopN >= 0; LoopN--){
910 LoopCnt = (Ntbl[(M>>1)-1]-1);
911 fly0P = ioptr + Nroot_1_ColInc + ((int)BRcnt[LoopN] >> BRshift)*(2*2) + stage;
912 IOP = ioptr + (LoopN<<(M+1)/2) * 2 + stage;
913 fly1P = IOP + ((int)BRcnt[LoopCnt] >> BRshift)*(2*2);
914 fly0r = *(fly0P);
915 fly0i = *(fly0P+1);
916 fly1r = *(fly0P+FlyOffsetA);
917 fly1i = *(fly0P+FlyOffsetAIm);
918 for (; LoopCnt > LoopN;){
919 fly2r = *(fly0P+2);
920 fly2i = *(fly0P+(2+1));
921 fly3r = *(fly0P+FlyOffsetB);
922 fly3i = *(fly0P+FlyOffsetBIm);
923 fly4r = *(fly1P);
924 fly4i = *(fly1P+1);
925 fly5r = *(fly1P+FlyOffsetA);
926 fly5i = *(fly1P+FlyOffsetAIm);
927 fly6r = *(fly1P+2);
928 fly6i = *(fly1P+(2+1));
929 fly7r = *(fly1P+FlyOffsetB);
930 fly7i = *(fly1P+FlyOffsetBIm);
931
932 t0r = fly0r + fly1r;
933 t0i = fly0i + fly1i;
934 fly1r = fly0r - fly1r;
935 fly1i = fly0i - fly1i;
936 t1r = fly2r + fly3r;
937 t1i = fly2i + fly3i;
938 fly3r = fly2r - fly3r;
939 fly3i = fly2i - fly3i;
940 fly0r = fly4r + fly5r;
941 fly0i = fly4i + fly5i;
942 fly5r = fly4r - fly5r;
943 fly5i = fly4i - fly5i;
944 fly2r = fly6r + fly7r;
945 fly2i = fly6i + fly7i;
946 fly7r = fly6r - fly7r;
947 fly7i = fly6i - fly7i;
948
949 *(fly1P) = scale*t0r;
950 *(fly1P+1) = scale*t0i;
951 *(fly1P+2) = scale*fly1r;
952 *(fly1P+(2+1)) = scale*fly1i;
953 *(fly1P+FlyOffsetA) = scale*t1r;
954 *(fly1P+FlyOffsetAIm) = scale*t1i;
955 *(fly1P+FlyOffsetB) = scale*fly3r;
956 *(fly1P+FlyOffsetBIm) = scale*fly3i;
957 *(fly0P) = scale*fly0r;
958 *(fly0P+1) = scale*fly0i;
959 *(fly0P+2) = scale*fly5r;
960 *(fly0P+(2+1)) = scale*fly5i;
961 *(fly0P+FlyOffsetA) = scale*fly2r;
962 *(fly0P+FlyOffsetAIm) = scale*fly2i;
963 *(fly0P+FlyOffsetB) = scale*fly7r;
964 *(fly0P+FlyOffsetBIm) = scale*fly7i;
965
966 fly0P -= Nrems2;
967 fly0r = *(fly0P);
968 fly0i = *(fly0P+1);
969 fly1r = *(fly0P+FlyOffsetA);
970 fly1i = *(fly0P+FlyOffsetAIm);
971 LoopCnt -= 1;
972 fly1P = (IOP + ((int)BRcnt[LoopCnt] >> BRshift)*(2*2));
973 };
974 fly2r = *(fly0P+2);
975 fly2i = *(fly0P+(2+1));
976 fly3r = *(fly0P+FlyOffsetB);
977 fly3i = *(fly0P+FlyOffsetBIm);
978
979 t0r = fly0r + fly1r;
980 t0i = fly0i + fly1i;
981 fly1r = fly0r - fly1r;
982 fly1i = fly0i - fly1i;
983 t1r = fly2r + fly3r;
984 t1i = fly2i + fly3i;
985 fly3r = fly2r - fly3r;
986 fly3i = fly2i - fly3i;
987
988 *(fly0P) = scale*t0r;
989 *(fly0P+1) = scale*t0i;
990 *(fly0P+2) = scale*fly1r;
991 *(fly0P+(2+1)) = scale*fly1i;
992 *(fly0P+FlyOffsetA) = scale*t1r;
993 *(fly0P+FlyOffsetAIm) = scale*t1i;
994 *(fly0P+FlyOffsetB) = scale*fly3r;
995 *(fly0P+FlyOffsetBIm) = scale*fly3i;
996
997 };
998 };
999
1000 /**** FFTC **************/
1001
1002 scale = 2.0;
1003
1004 NdiffU = 2;
1005 Flyinc = (NdiffU);
1006
1007 NsameU4 = Ntbl[M]/4;
1008 LoopN = Ntbl[M-3];
1009
1010 stage = ((M-1)/3);
1011 if ( (M-1-(stage * 3)) != 0 ){
1012 FlyOffsetA = Flyinc << 2;
1013 FlyOffsetB = FlyOffsetA << 1;
1014 FlyOffsetAIm = FlyOffsetA + 1;
1015 FlyOffsetBIm = FlyOffsetB + 1;
1016 if ( (M-1-(stage * 3)) == 1 ){
1017 /* 1 radix 2 stage */
1018
1019 IOP = ioptr;
1020 fly0P = IOP;
1021 fly1P = (IOP+Flyinc);
1022 fly2P = (fly1P+Flyinc);
1023 fly3P = (fly2P+Flyinc);
1024
1025 /* Butterflys */
1026 /*
1027 t0 - - t0
1028 t1 - - t1
1029 f2 - 1- f5
1030 f1 - -i- f7
1031 f4 - - f4
1032 f0 - - f0
1033 f6 - 1- f2
1034 f3 - -i- f1
1035 */
1036
1037 for (LoopCnt = LoopN; LoopCnt > 0 ; LoopCnt--){
1038 t0r = *(fly0P);
1039 t0i = *(fly0P + 1);
1040 t1r = *(fly1P);
1041 t1i = *(fly1P + 1);
1042 fly2r = *(fly2P);
1043 fly2i = *(fly2P + 1);
1044 fly1r = *(fly3P);
1045 fly1i = *(fly3P + 1);
1046 fly4r = *(fly0P + FlyOffsetA);
1047 fly4i = *(fly0P + FlyOffsetAIm);
1048 fly0r = *(fly1P + FlyOffsetA);
1049 fly0i = *(fly1P + FlyOffsetAIm);
1050 fly6r = *(fly2P + FlyOffsetA);
1051 fly6i = *(fly2P + FlyOffsetAIm);
1052 fly3r = *(fly3P + FlyOffsetA);
1053 fly3i = *(fly3P + FlyOffsetAIm);
1054
1055 fly5r = t0r - fly2r;
1056 fly5i = t0i - fly2i;
1057 t0r = t0r + fly2r;
1058 t0i = t0i + fly2i;
1059
1060 fly7r = t1r + fly1i;
1061 fly7i = t1i - fly1r;
1062 t1r = t1r - fly1i;
1063 t1i = t1i + fly1r;
1064
1065 fly2r = fly4r - fly6r;
1066 fly2i = fly4i - fly6i;
1067 fly4r = fly4r + fly6r;
1068 fly4i = fly4i + fly6i;
1069
1070 fly1r = fly0r + fly3i;
1071 fly1i = fly0i - fly3r;
1072 fly0r = fly0r - fly3i;
1073 fly0i = fly0i + fly3r;
1074
1075 *(fly2P) = fly5r;
1076 *(fly2P + 1) = fly5i;
1077 *(fly0P) = t0r;
1078 *(fly0P + 1) = t0i;
1079 *(fly3P) = fly7r;
1080 *(fly3P + 1) = fly7i;
1081 *(fly1P) = t1r;
1082 *(fly1P + 1) = t1i;
1083 *(fly2P + FlyOffsetA) = fly2r;
1084 *(fly2P + FlyOffsetAIm) = fly2i;
1085 *(fly0P + FlyOffsetA) = fly4r;
1086 *(fly0P + FlyOffsetAIm) = fly4i;
1087 *(fly3P + FlyOffsetA) = fly1r;
1088 *(fly3P + FlyOffsetAIm) = fly1i;
1089 *(fly1P + FlyOffsetA) = fly0r;
1090 *(fly1P + FlyOffsetAIm) = fly0i;
1091
1092 fly0P = (fly0P + FlyOffsetB);
1093 fly1P = (fly1P + FlyOffsetB);
1094 fly2P = (fly2P + FlyOffsetB);
1095 fly3P = (fly3P + FlyOffsetB);
1096 };
1097
1098 NsameU4 >>= 1;
1099 LoopN >>= 1;
1100 NdiffU <<= 1;
1101 Flyinc = Flyinc << 1;
1102 }
1103 else{
1104 /* 1 radix 4 stage */
1105 IOP = ioptr;
1106
1107 U3r = 0.7071067811865475244008443621; /* sqrt(0.5); */
1108 U3i = U3r;
1109 fly0P = IOP;
1110 fly1P = (IOP+Flyinc);
1111 fly2P = (fly1P+Flyinc);
1112 fly3P = (fly2P+Flyinc);
1113
1114 /* Butterflys */
1115 /*
1116 t0 - - t0 - - t0
1117 t1 - - t1 - - t1
1118 f2 - 1- f5 - - f5
1119 f1 - -i- f7 - - f7
1120 f4 - - f4 - 1- f6
1121 f0 - - f0 -U3 - f3
1122 f6 - 1- f2 - -i- f4
1123 f3 - -i- f1 -U3a- f2
1124 */
1125
1126 for (LoopCnt = LoopN; LoopCnt > 0 ; LoopCnt--){
1127 t0r = *(fly0P);
1128 t0i = *(fly0P + 1);
1129 t1r = *(fly1P);
1130 t1i = *(fly1P + 1);
1131 fly2r = *(fly2P);
1132 fly2i = *(fly2P + 1);
1133 fly1r = *(fly3P);
1134 fly1i = *(fly3P + 1);
1135 fly4r = *(fly0P + FlyOffsetA);
1136 fly4i = *(fly0P + FlyOffsetAIm);
1137 fly0r = *(fly1P + FlyOffsetA);
1138 fly0i = *(fly1P + FlyOffsetAIm);
1139 fly6r = *(fly2P + FlyOffsetA);
1140 fly6i = *(fly2P + FlyOffsetAIm);
1141 fly3r = *(fly3P + FlyOffsetA);
1142 fly3i = *(fly3P + FlyOffsetAIm);
1143
1144 fly5r = t0r - fly2r;
1145 fly5i = t0i - fly2i;
1146 t0r = t0r + fly2r;
1147 t0i = t0i + fly2i;
1148
1149 fly7r = t1r + fly1i;
1150 fly7i = t1i - fly1r;
1151 t1r = t1r - fly1i;
1152 t1i = t1i + fly1r;
1153
1154 fly2r = fly4r - fly6r;
1155 fly2i = fly4i - fly6i;
1156 fly4r = fly4r + fly6r;
1157 fly4i = fly4i + fly6i;
1158
1159 fly1r = fly0r + fly3i;
1160 fly1i = fly0i - fly3r;
1161 fly0r = fly0r - fly3i;
1162 fly0i = fly0i + fly3r;
1163
1164 fly6r = t0r - fly4r;
1165 fly6i = t0i - fly4i;
1166 t0r = t0r + fly4r;
1167 t0i = t0i + fly4i;
1168
1169 fly3r = fly5r + fly2i;
1170 fly3i = fly5i - fly2r;
1171 fly5r = fly5r - fly2i;
1172 fly5i = fly5i + fly2r;
1173
1174 fly4r = t1r - U3r * fly0r;
1175 fly4r = fly4r + U3i * fly0i;
1176 fly4i = t1i - U3i * fly0r;
1177 fly4i = fly4i - U3r * fly0i;
1178 t1r = scale * t1r - fly4r;
1179 t1i = scale * t1i - fly4i;
1180
1181 fly2r = fly7r + U3i * fly1r;
1182 fly2r = fly2r + U3r * fly1i;
1183 fly2i = fly7i - U3r * fly1r;
1184 fly2i = fly2i + U3i * fly1i;
1185 fly7r = scale * fly7r - fly2r;
1186 fly7i = scale * fly7i - fly2i;
1187
1188 *(fly0P + FlyOffsetA) = fly6r;
1189 *(fly0P + FlyOffsetAIm) = fly6i;
1190 *(fly0P) = t0r;
1191 *(fly0P + 1) = t0i;
1192 *(fly2P + FlyOffsetA) = fly3r;
1193 *(fly2P + FlyOffsetAIm) = fly3i;
1194 *(fly2P) = fly5r;
1195 *(fly2P + 1) = fly5i;
1196 *(fly1P + FlyOffsetA) = fly4r;
1197 *(fly1P + FlyOffsetAIm) = fly4i;
1198 *(fly1P) = t1r;
1199 *(fly1P + 1) = t1i;
1200 *(fly3P + FlyOffsetA) = fly2r;
1201 *(fly3P + FlyOffsetAIm) = fly2i;
1202 *(fly3P) = fly7r;
1203 *(fly3P + 1) = fly7i;
1204
1205 fly0P = (fly0P + FlyOffsetB);
1206 fly1P = (fly1P + FlyOffsetB);
1207 fly2P = (fly2P + FlyOffsetB);
1208 fly3P = (fly3P + FlyOffsetB);
1209
1210 };
1211
1212 NsameU4 >>= 2;
1213 LoopN >>= 2;
1214 NdiffU <<= 2;
1215 Flyinc = Flyinc << 2;
1216 };
1217 };
1218
1219 NsameU2 = NsameU4 >> 1;
1220 NsameU1 = NsameU2 >> 1;
1221 Flyinc <<= 1;
1222 FlyOffsetA = Flyinc << 2;
1223 FlyOffsetB = FlyOffsetA << 1;
1224 FlyOffsetAIm = FlyOffsetA + 1;
1225 FlyOffsetBIm = FlyOffsetB + 1;
1226 LoopN >>= 1;
1227
1228 /* ****** RADIX 8 Stages */
1229 for (stage = stage<<1; stage > 0 ; stage--){
1230
1231 /* an fft stage is done in two parts to ease use of the single quadrant cos table */
1232
1233 /* fftcalc1(iobuf, Utbl, N, NdiffU, LoopN); */
1234 if(!(stage&1)){
1235 U0rP = &Utbl[0];
1236 U0iP = &Utbl[Ntbl[M-2]];
1237 U1rP = U0rP;
1238 U1iP = U0iP;
1239 U2rP = U0rP;
1240 U2iP = U0iP;
1241 U3offset = (Ntbl[M]) / 8;
1242
1243 IOP = ioptr;
1244
1245 U0r = *U0rP,
1246 U0i = *U0iP;
1247 U1r = *U1rP,
1248 U1i = *U1iP;
1249 U2r = *U2rP,
1250 U2i = *U2iP;
1251 U3r = *( U2rP + U3offset);
1252 U3i = *( U2iP - U3offset);
1253 }
1254
1255 fly0P = IOP;
1256 fly1P = (IOP+Flyinc);
1257 fly2P = (fly1P+Flyinc);
1258 fly3P = (fly2P+Flyinc);
1259
1260 for (diffUcnt = (NdiffU)>>1; diffUcnt > 0; diffUcnt--){
1261
1262 /* Butterflys */
1263 /*
1264 f0 - - t0 - - t0 - - t0
1265 f1 -U0 - t1 - - t1 - - t1
1266 f2 - - f2 -U1 - f5 - - f3
1267 f3 -U0 - f1 -U1a- f7 - - f7
1268 f4 - - f4 - - f4 -U2 - f6
1269 f5 -U0 - f0 - - f0 -U3 - f4
1270 f6 - - f6 -U1 - f2 -U2a- f2
1271 f7 -U0 - f3 -U1a- f1 -U3a- f5
1272 */
1273
1274 fly0r = *(IOP);
1275 fly0i = *(IOP+1);
1276 fly1r = *(fly1P);
1277 fly1i = *(fly1P+1);
1278
1279 for (LoopCnt = LoopN-1; LoopCnt > 0 ; LoopCnt--){
1280
1281 fly2r = *(fly2P);
1282 fly2i = *(fly2P + 1);
1283 fly3r = *(fly3P);
1284 fly3i = *(fly3P + 1);
1285 fly4r = *(fly0P + FlyOffsetA);
1286 fly4i = *(fly0P + FlyOffsetAIm);
1287 fly5r = *(fly1P + FlyOffsetA);
1288 fly5i = *(fly1P + FlyOffsetAIm);
1289 fly6r = *(fly2P + FlyOffsetA);
1290 fly6i = *(fly2P + FlyOffsetAIm);
1291 fly7r = *(fly3P + FlyOffsetA);
1292 fly7i = *(fly3P + FlyOffsetAIm);
1293
1294 t1r = fly0r - U0r * fly1r;
1295 t1r = t1r + U0i * fly1i;
1296 t1i = fly0i - U0i * fly1r;
1297 t1i = t1i - U0r * fly1i;
1298 t0r = scale * fly0r - t1r;
1299 t0i = scale * fly0i - t1i;
1300
1301 fly1r = fly2r - U0r * fly3r;
1302 fly1r = fly1r + U0i * fly3i;
1303 fly1i = fly2i - U0i * fly3r;
1304 fly1i = fly1i - U0r * fly3i;
1305 fly2r = scale * fly2r - fly1r;
1306 fly2i = scale * fly2i - fly1i;
1307
1308 fly0r = fly4r - U0r * fly5r;
1309 fly0r = fly0r + U0i * fly5i;
1310 fly0i = fly4i - U0i * fly5r;
1311 fly0i = fly0i - U0r * fly5i;
1312 fly4r = scale * fly4r - fly0r;
1313 fly4i = scale * fly4i - fly0i;
1314
1315 fly3r = fly6r - U0r * fly7r;
1316 fly3r = fly3r + U0i * fly7i;
1317 fly3i = fly6i - U0i * fly7r;
1318 fly3i = fly3i - U0r * fly7i;
1319 fly6r = scale * fly6r - fly3r;
1320 fly6i = scale * fly6i - fly3i;
1321
1322
1323 fly5r = t0r - U1r * fly2r;
1324 fly5r = fly5r + U1i * fly2i;
1325 fly5i = t0i - U1i * fly2r;
1326 fly5i = fly5i - U1r * fly2i;
1327 t0r = scale * t0r - fly5r;
1328 t0i = scale * t0i - fly5i;
1329
1330 fly7r = t1r + U1i * fly1r;
1331 fly7r = fly7r + U1r * fly1i;
1332 fly7i = t1i - U1r * fly1r;
1333 fly7i = fly7i + U1i * fly1i;
1334 t1r = scale * t1r - fly7r;
1335 t1i = scale * t1i - fly7i;
1336
1337 fly2r = fly4r - U1r * fly6r;
1338 fly2r = fly2r + U1i * fly6i;
1339 fly2i = fly4i - U1i * fly6r;
1340 fly2i = fly2i - U1r * fly6i;
1341 fly4r = scale * fly4r - fly2r;
1342 fly4i = scale * fly4i - fly2i;
1343
1344 fly1r = fly0r + U1i * fly3r;
1345 fly1r = fly1r + U1r * fly3i;
1346 fly1i = fly0i - U1r * fly3r;
1347 fly1i = fly1i + U1i * fly3i;
1348 fly0r = scale * fly0r - fly1r;
1349 fly0i = scale * fly0i - fly1i;
1350
1351 fly6r = t0r - U2r * fly4r;
1352 fly6r = fly6r + U2i * fly4i;
1353 fly6i = t0i - U2i * fly4r;
1354 fly6i = fly6i - U2r * fly4i;
1355 t0r = scale * t0r - fly6r;
1356 t0i = scale * t0i - fly6i;
1357
1358 fly3r = fly5r - U2i * fly2r;
1359 fly3r = fly3r - U2r * fly2i;
1360 fly3i = fly5i + U2r * fly2r;
1361 fly3i = fly3i - U2i * fly2i;
1362 fly2r = scale * fly5r - fly3r;
1363 fly2i = scale * fly5i - fly3i;
1364
1365 fly4r = t1r - U3r * fly0r;
1366 fly4r = fly4r + U3i * fly0i;
1367 fly4i = t1i - U3i * fly0r;
1368 fly4i = fly4i - U3r * fly0i;
1369 t1r = scale * t1r - fly4r;
1370 t1i = scale * t1i - fly4i;
1371
1372 fly5r = fly7r + U3i * fly1r;
1373 fly5r = fly5r + U3r * fly1i;
1374 fly5i = fly7i - U3r * fly1r;
1375 fly5i = fly5i + U3i * fly1i;
1376 fly7r = scale * fly7r - fly5r;
1377 fly7i = scale * fly7i - fly5i;
1378
1379 *(fly0P + FlyOffsetA) = fly6r;
1380 *(fly0P + FlyOffsetAIm) = fly6i;
1381 *(fly0P) = t0r;
1382 *(fly0P + 1) = t0i;
1383 *(fly2P) = fly3r;
1384 *(fly2P + 1) = fly3i;
1385 *(fly2P + FlyOffsetA) = fly2r;
1386 *(fly2P + FlyOffsetAIm) = fly2i;
1387
1388 fly0r = *(fly0P + FlyOffsetB);
1389 fly0i = *(fly0P + FlyOffsetBIm);
1390
1391 *(fly1P + FlyOffsetA) = fly4r;
1392 *(fly1P + FlyOffsetAIm) = fly4i;
1393 *(fly1P) = t1r;
1394 *(fly1P + 1) = t1i;
1395
1396 fly1r = *(fly1P + FlyOffsetB);
1397 fly1i = *(fly1P + FlyOffsetBIm);
1398
1399 *(fly3P + FlyOffsetA) = fly5r;
1400 *(fly3P + FlyOffsetAIm) = fly5i;
1401 *(fly3P) = fly7r;
1402 *(fly3P + 1) = fly7i;
1403
1404 fly0P = (fly0P + FlyOffsetB);
1405 fly1P = (fly1P + FlyOffsetB);
1406 fly2P = (fly2P + FlyOffsetB);
1407 fly3P = (fly3P + FlyOffsetB);
1408
1409 };
1410 fly2r = *(fly2P);
1411 fly2i = *(fly2P + 1);
1412 fly3r = *(fly3P);
1413 fly3i = *(fly3P + 1);
1414 fly4r = *(fly0P + FlyOffsetA);
1415 fly4i = *(fly0P + FlyOffsetAIm);
1416 fly5r = *(fly1P + FlyOffsetA);
1417 fly5i = *(fly1P + FlyOffsetAIm);
1418 fly6r = *(fly2P + FlyOffsetA);
1419 fly6i = *(fly2P + FlyOffsetAIm);
1420 fly7r = *(fly3P + FlyOffsetA);
1421 fly7i = *(fly3P + FlyOffsetAIm);
1422
1423 t1r = fly0r - U0r * fly1r;
1424 t1r = t1r + U0i * fly1i;
1425 t1i = fly0i - U0i * fly1r;
1426 t1i = t1i - U0r * fly1i;
1427 t0r = scale * fly0r - t1r;
1428 t0i = scale * fly0i - t1i;
1429
1430 fly1r = fly2r - U0r * fly3r;
1431 fly1r = fly1r + U0i * fly3i;
1432 fly1i = fly2i - U0i * fly3r;
1433 fly1i = fly1i - U0r * fly3i;
1434 fly2r = scale * fly2r - fly1r;
1435 fly2i = scale * fly2i - fly1i;
1436
1437 fly0r = fly4r - U0r * fly5r;
1438 fly0r = fly0r + U0i * fly5i;
1439 fly0i = fly4i - U0i * fly5r;
1440 fly0i = fly0i - U0r * fly5i;
1441 fly4r = scale * fly4r - fly0r;
1442 fly4i = scale * fly4i - fly0i;
1443
1444 fly3r = fly6r - U0r * fly7r;
1445 fly3r = fly3r + U0i * fly7i;
1446 fly3i = fly6i - U0i * fly7r;
1447 fly3i = fly3i - U0r * fly7i;
1448 fly6r = scale * fly6r - fly3r;
1449 fly6i = scale * fly6i - fly3i;
1450
1451 fly5r = t0r - U1r * fly2r;
1452 fly5r = fly5r + U1i * fly2i;
1453 fly5i = t0i - U1i * fly2r;
1454 fly5i = fly5i - U1r * fly2i;
1455 t0r = scale * t0r - fly5r;
1456 t0i = scale * t0i - fly5i;
1457
1458 fly7r = t1r + U1i * fly1r;
1459 fly7r = fly7r + U1r * fly1i;
1460 fly7i = t1i - U1r * fly1r;
1461 fly7i = fly7i + U1i * fly1i;
1462 t1r = scale * t1r - fly7r;
1463 t1i = scale * t1i - fly7i;
1464
1465 fly2r = fly4r - U1r * fly6r;
1466 fly2r = fly2r + U1i * fly6i;
1467 fly2i = fly4i - U1i * fly6r;
1468 fly2i = fly2i - U1r * fly6i;
1469 fly4r = scale * fly4r - fly2r;
1470 fly4i = scale * fly4i - fly2i;
1471
1472 fly1r = fly0r + U1i * fly3r;
1473 fly1r = fly1r + U1r * fly3i;
1474 fly1i = fly0i - U1r * fly3r;
1475 fly1i = fly1i + U1i * fly3i;
1476 fly0r = scale * fly0r - fly1r;
1477 fly0i = scale * fly0i - fly1i;
1478
1479 U0i = *(U0iP = (U0iP - NsameU4));
1480 U0r = *(U0rP = (U0rP + NsameU4));
1481 U1r = *(U1rP = (U1rP + NsameU2));
1482 U1i = *(U1iP = (U1iP - NsameU2));
1483 if(stage&1)
1484 U0r = - U0r;
1485
1486 fly6r = t0r - U2r * fly4r;
1487 fly6r = fly6r + U2i * fly4i;
1488 fly6i = t0i - U2i * fly4r;
1489 fly6i = fly6i - U2r * fly4i;
1490 t0r = scale * t0r - fly6r;
1491 t0i = scale * t0i - fly6i;
1492
1493 fly3r = fly5r - U2i * fly2r;
1494 fly3r = fly3r - U2r * fly2i;
1495 fly3i = fly5i + U2r * fly2r;
1496 fly3i = fly3i - U2i * fly2i;
1497 fly2r = scale * fly5r - fly3r;
1498 fly2i = scale * fly5i - fly3i;
1499
1500 fly4r = t1r - U3r * fly0r;
1501 fly4r = fly4r + U3i * fly0i;
1502 fly4i = t1i - U3i * fly0r;
1503 fly4i = fly4i - U3r * fly0i;
1504 t1r = scale * t1r - fly4r;
1505 t1i = scale * t1i - fly4i;
1506
1507 fly5r = fly7r + U3i * fly1r;
1508 fly5r = fly5r + U3r * fly1i;
1509 fly5i = fly7i - U3r * fly1r;
1510 fly5i = fly5i + U3i * fly1i;
1511 fly7r = scale * fly7r - fly5r;
1512 fly7i = scale * fly7i - fly5i;
1513
1514 *(fly0P + FlyOffsetA) = fly6r;
1515 *(fly0P + FlyOffsetAIm) = fly6i;
1516 *(fly0P) = t0r;
1517 *(fly0P + 1) = t0i;
1518
1519 U2r = *(U2rP = (U2rP + NsameU1));
1520 U2i = *(U2iP = (U2iP - NsameU1));
1521
1522 *(fly2P) = fly3r;
1523 *(fly2P + 1) = fly3i;
1524 *(fly2P + FlyOffsetA) = fly2r;
1525 *(fly2P + FlyOffsetAIm) = fly2i;
1526 *(fly1P + FlyOffsetA) = fly4r;
1527 *(fly1P + FlyOffsetAIm) = fly4i;
1528 *(fly1P) = t1r;
1529 *(fly1P + 1) = t1i;
1530
1531 U3r = *( U2rP + U3offset);
1532 U3i = *( U2iP - U3offset);
1533
1534 *(fly3P + FlyOffsetA) = fly5r;
1535 *(fly3P + FlyOffsetAIm) = fly5i;
1536 *(fly3P) = fly7r;
1537 *(fly3P + 1) = fly7i;
1538
1539 IOP = IOP + 2;
1540 fly0P = IOP;
1541 fly1P = (IOP+Flyinc);
1542 fly2P = (fly1P+Flyinc);
1543 fly3P = (fly2P+Flyinc);
1544 };
1545 NsameU4 = -NsameU4;
1546
1547 if(stage&1){
1548 LoopN >>= 3;
1549 NsameU1 >>= 3;
1550 NsameU2 >>= 3;
1551 NsameU4 >>= 3;
1552 NdiffU <<= 3;
1553 Flyinc = Flyinc << 3;
1554 FlyOffsetA <<= 3;
1555 FlyOffsetB <<= 3;
1556 FlyOffsetAIm = FlyOffsetA + 1;
1557 FlyOffsetBIm = FlyOffsetB + 1;
1558 }
1559 }
1560 ioptr += 2*Ntbl[M];
1561 }
1562 }
1563
1564 /* rffts */
1565 /* compute multiple real input FFTs on 'Rows' consecutively stored arrays */
1566 /* ioptr = pointer to the data in and out */
1567 /* M = log2 of fft size */
1568 /* Rows = number of FFTs to compute */
1569 /* Utbl = Pointer to cosine table */
1570
1571 void rffts(float *ioptr, long M, long Rows, float *Utbl){
1572 /* Compute in-place real fft on the rows of the input array */
1573 /* INPUTS */
1574 /* M = log2 of fft size */
1575 /* *ioptr = real input data array */
1576 /* *Utbl = cosine table */
1577 /* Rows = number of rows in ioptr array (use Rows of 1 if ioptr is a 1 dimensional array) */
1578 /* OUTPUTS */
1579 /* *ioptr = output data array in the following order */
1580 /* Re(x[0]), Re(x[N/2]), Re(x[1]), Im(x[1]), Re(x[2]), Im(x[2]), ... Re(x[N/2-1]), Im(x[N/2-1]). */
1581
1582 long Flyinc;
1583 long FlyOffsetA;
1584 long FlyOffsetAIm;
1585 long FlyOffsetB;
1586 long FlyOffsetBIm;
1587 long NsameU1;
1588 long NsameU2;
1589 long NsameU4;
1590 long diffUcnt;
1591 long LoopCnt;
1592 float scale;
1593 float fly0r;
1594 float fly0i;
1595 float fly1r;
1596 float fly1i;
1597 float fly2r;
1598 float fly2i;
1599 float fly3r;
1600 float fly3i;
1601 float fly4r;
1602 float fly4i;
1603 float fly5r;
1604 float fly5i;
1605 float fly6r;
1606 float fly6i;
1607 float fly7r;
1608 float fly7i;
1609 float U0r;
1610 float U0i;
1611 float U1r;
1612 float U1i;
1613 float U2r;
1614 float U2i;
1615 float U3r;
1616 float U3i;
1617 float t0r;
1618 float t0i;
1619 float t1r;
1620 float t1i;
1621
1622 float *fly0P;
1623 float *fly1P;
1624 float *fly2P;
1625 float *fly3P;
1626
1627 float *U0rP;
1628 float *U0iP;
1629 float *U1rP;
1630 float *U1iP;
1631 float *U2rP;
1632 float *U2iP;
1633 float *IOP;
1634 long U3offset;
1635
1636 long stage;
1637 long NdiffU;
1638 long LoopN;
1639
1640 const long BRshift = MAXMROOT - ((M-1)>>1); /* for RFFT */
1641 const long Nrems2 = Ntbl[(M-1)-((M-1)>>1)+1]; /* for RFFT */
1642 const long Nroot_1_ColInc = (Ntbl[(M-1)-1]-Ntbl[(M-1)-((M-1)>>1)])*2; /* for RFFT */
1643
1644 for (;Rows>0;Rows--){
1645
1646 M=M-1; /* for RFFT */
1647
1648 FlyOffsetA = Ntbl[M] * 2/2;
1649 FlyOffsetAIm = FlyOffsetA + 1;
1650 FlyOffsetB = FlyOffsetA + 2;
1651 FlyOffsetBIm = FlyOffsetB + 1;
1652
1653 /* BitrevR2 ** bit reverse shuffle and first radix 2 stage ******/
1654
1655 scale = 0.5;
1656 for (stage = 0; stage < Ntbl[M-(M>>1)]*2; stage += Ntbl[M>>1]*2){
1657 for (LoopN = (Ntbl[(M>>1)-1]-1); LoopN >= 0; LoopN--){
1658 LoopCnt = (Ntbl[(M>>1)-1]-1);
1659 fly0P = ioptr + Nroot_1_ColInc + ((int)BRcnt[LoopN] >> BRshift)*(2*2) + stage;
1660 IOP = ioptr + (LoopN<<(M+1)/2) * 2 + stage;
1661 fly1P = IOP + ((int)BRcnt[LoopCnt] >> BRshift)*(2*2);
1662 fly0r = *fly0P;
1663 fly0i = *(fly0P+1);
1664 fly1r = *(fly0P+FlyOffsetA);
1665 fly1i = *(fly0P+FlyOffsetAIm);
1666 for (; LoopCnt > LoopN;){
1667 fly2r = *(fly0P+2);
1668 fly2i = *(fly0P+(2+1));
1669 fly3r = *(fly0P+FlyOffsetB);
1670 fly3i = *(fly0P+FlyOffsetBIm);
1671 fly4r = *fly1P;
1672 fly4i = *(fly1P+1);
1673 fly5r = *(fly1P+FlyOffsetA);
1674 fly5i = *(fly1P+FlyOffsetAIm);
1675 fly6r = *(fly1P+2);
1676 fly6i = *(fly1P+(2+1));
1677 fly7r = *(fly1P+FlyOffsetB);
1678 fly7i = *(fly1P+FlyOffsetBIm);
1679
1680 t0r = fly0r + fly1r;
1681 t0i = fly0i + fly1i;
1682 fly1r = fly0r - fly1r;
1683 fly1i = fly0i - fly1i;
1684 t1r = fly2r + fly3r;
1685 t1i = fly2i + fly3i;
1686 fly3r = fly2r - fly3r;
1687 fly3i = fly2i - fly3i;
1688 fly0r = fly4r + fly5r;
1689 fly0i = fly4i + fly5i;
1690 fly5r = fly4r - fly5r;
1691 fly5i = fly4i - fly5i;
1692 fly2r = fly6r + fly7r;
1693 fly2i = fly6i + fly7i;
1694 fly7r = fly6r - fly7r;
1695 fly7i = fly6i - fly7i;
1696
1697 *fly1P = scale*t0r;
1698 *(fly1P+1) = scale*t0i;
1699 *(fly1P+2) = scale*fly1r;
1700 *(fly1P+(2+1)) = scale*fly1i;
1701 *(fly1P+FlyOffsetA) = scale*t1r;
1702 *(fly1P+FlyOffsetAIm) = scale*t1i;
1703 *(fly1P+FlyOffsetB) = scale*fly3r;
1704 *(fly1P+FlyOffsetBIm) = scale*fly3i;
1705 *fly0P = scale*fly0r;
1706 *(fly0P+1) = scale*fly0i;
1707 *(fly0P+2) = scale*fly5r;
1708 *(fly0P+(2+1)) = scale*fly5i;
1709 *(fly0P+FlyOffsetA) = scale*fly2r;
1710 *(fly0P+FlyOffsetAIm) = scale*fly2i;
1711 *(fly0P+FlyOffsetB) = scale*fly7r;
1712 *(fly0P+FlyOffsetBIm) = scale*fly7i;
1713
1714 fly0P -= Nrems2;
1715 fly0r = *fly0P;
1716 fly0i = *(fly0P+1);
1717 fly1r = *(fly0P+FlyOffsetA);
1718 fly1i = *(fly0P+FlyOffsetAIm);
1719 LoopCnt -= 1;
1720 fly1P = (IOP + ((int)BRcnt[LoopCnt] >> BRshift)*(2*2));
1721 };
1722 fly2r = *(fly0P+2);
1723 fly2i = *(fly0P+(2+1));
1724 fly3r = *(fly0P+FlyOffsetB);
1725 fly3i = *(fly0P+FlyOffsetBIm);
1726
1727 t0r = fly0r + fly1r;
1728 t0i = fly0i + fly1i;
1729 fly1r = fly0r - fly1r;
1730 fly1i = fly0i - fly1i;
1731 t1r = fly2r + fly3r;
1732 t1i = fly2i + fly3i;
1733 fly3r = fly2r - fly3r;
1734 fly3i = fly2i - fly3i;
1735
1736 *fly0P = scale*t0r;
1737 *(fly0P+1) = scale*t0i;
1738 *(fly0P+2) = scale*fly1r;
1739 *(fly0P+(2+1)) = scale*fly1i;
1740 *(fly0P+FlyOffsetA) = scale*t1r;
1741 *(fly0P+FlyOffsetAIm) = scale*t1i;
1742 *(fly0P+FlyOffsetB) = scale*fly3r;
1743 *(fly0P+FlyOffsetBIm) = scale*fly3i;
1744
1745 };
1746 };
1747
1748
1749
1750 /**** FFTC **************/
1751
1752 scale = 2.0;
1753
1754 NdiffU = 2;
1755 Flyinc = (NdiffU);
1756
1757 NsameU4 = Ntbl[M+1]/4; /* for RFFT */
1758 LoopN = Ntbl[M-3];
1759
1760 stage = ((M-1)/3);
1761 if ( (M-1-(stage * 3)) != 0 ){
1762 FlyOffsetA = Flyinc << 2;
1763 FlyOffsetB = FlyOffsetA << 1;
1764 FlyOffsetAIm = FlyOffsetA + 1;
1765 FlyOffsetBIm = FlyOffsetB + 1;
1766 if ( (M-1-(stage * 3)) == 1 ){
1767 /* 1 radix 2 stage */
1768
1769 IOP = ioptr;
1770 fly0P = IOP;
1771 fly1P = (IOP+Flyinc);
1772 fly2P = (fly1P+Flyinc);
1773 fly3P = (fly2P+Flyinc);
1774
1775 /* Butterflys */
1776 /*
1777 t0 - - t0
1778 t1 - - t1
1779 f2 - 1- f5
1780 f1 - -i- f7
1781 f4 - - f4
1782 f0 - - f0
1783 f6 - 1- f2
1784 f3 - -i- f1
1785 */
1786
1787 for (LoopCnt = LoopN; LoopCnt > 0 ; LoopCnt--){
1788 t0r = *(fly0P);
1789 t0i = *(fly0P + 1);
1790 t1r = *(fly1P);
1791 t1i = *(fly1P + 1);
1792 fly2r = *(fly2P);
1793 fly2i = *(fly2P + 1);
1794 fly1r = *(fly3P);
1795 fly1i = *(fly3P + 1);
1796 fly4r = *(fly0P + FlyOffsetA);
1797 fly4i = *(fly0P + FlyOffsetAIm);
1798 fly0r = *(fly1P + FlyOffsetA);
1799 fly0i = *(fly1P + FlyOffsetAIm);
1800 fly6r = *(fly2P + FlyOffsetA);
1801 fly6i = *(fly2P + FlyOffsetAIm);
1802 fly3r = *(fly3P + FlyOffsetA);
1803 fly3i = *(fly3P + FlyOffsetAIm);
1804
1805 fly5r = t0r - fly2r;
1806 fly5i = t0i - fly2i;
1807 t0r = t0r + fly2r;
1808 t0i = t0i + fly2i;
1809
1810 fly7r = t1r - fly1i;
1811 fly7i = t1i + fly1r;
1812 t1r = t1r + fly1i;
1813 t1i = t1i - fly1r;
1814
1815 fly2r = fly4r - fly6r;
1816 fly2i = fly4i - fly6i;
1817 fly4r = fly4r + fly6r;
1818 fly4i = fly4i + fly6i;
1819
1820 fly1r = fly0r - fly3i;
1821 fly1i = fly0i + fly3r;
1822 fly0r = fly0r + fly3i;
1823 fly0i = fly0i - fly3r;
1824
1825 *(fly2P) = fly5r;
1826 *(fly2P + 1) = fly5i;
1827 *(fly0P) = t0r;
1828 *(fly0P + 1) = t0i;
1829 *(fly3P) = fly7r;
1830 *(fly3P + 1) = fly7i;
1831 *(fly1P) = t1r;
1832 *(fly1P + 1) = t1i;
1833 *(fly2P + FlyOffsetA) = fly2r;
1834 *(fly2P + FlyOffsetAIm) = fly2i;
1835 *(fly0P + FlyOffsetA) = fly4r;
1836 *(fly0P + FlyOffsetAIm) = fly4i;
1837 *(fly3P + FlyOffsetA) = fly1r;
1838 *(fly3P + FlyOffsetAIm) = fly1i;
1839 *(fly1P + FlyOffsetA) = fly0r;
1840 *(fly1P + FlyOffsetAIm) = fly0i;
1841
1842 fly0P = (fly0P + FlyOffsetB);
1843 fly1P = (fly1P + FlyOffsetB);
1844 fly2P = (fly2P + FlyOffsetB);
1845 fly3P = (fly3P + FlyOffsetB);
1846 };
1847
1848 NsameU4 >>= 1;
1849 LoopN >>= 1;
1850 NdiffU <<= 1;
1851 Flyinc = Flyinc << 1;
1852 }
1853 else{
1854 /* 1 radix 4 stage */
1855 IOP = ioptr;
1856
1857 U3r = 0.7071067811865475244008443621; /* sqrt(0.5); */
1858 U3i = U3r;
1859 fly0P = IOP;
1860 fly1P = (IOP+Flyinc);
1861 fly2P = (fly1P+Flyinc);
1862 fly3P = (fly2P+Flyinc);
1863
1864 /* Butterflys */
1865 /*
1866 t0 - - t0 - - t0
1867 t1 - - t1 - - t1
1868 f2 - 1- f5 - - f5
1869 f1 - -i- f7 - - f7
1870 f4 - - f4 - 1- f6
1871 f0 - - f0 -U3 - f3
1872 f6 - 1- f2 - -i- f4
1873 f3 - -i- f1 -U3a- f2
1874 */
1875
1876 for (LoopCnt = LoopN; LoopCnt > 0 ; LoopCnt--){
1877 t0r = *(fly0P);
1878 t0i = *(fly0P + 1);
1879 t1r = *(fly1P);
1880 t1i = *(fly1P + 1);
1881 fly2r = *(fly2P);
1882 fly2i = *(fly2P + 1);
1883 fly1r = *(fly3P);
1884 fly1i = *(fly3P + 1);
1885 fly4r = *(fly0P + FlyOffsetA);
1886 fly4i = *(fly0P + FlyOffsetAIm);
1887 fly0r = *(fly1P + FlyOffsetA);
1888 fly0i = *(fly1P + FlyOffsetAIm);
1889 fly6r = *(fly2P + FlyOffsetA);
1890 fly6i = *(fly2P + FlyOffsetAIm);
1891 fly3r = *(fly3P + FlyOffsetA);
1892 fly3i = *(fly3P + FlyOffsetAIm);
1893
1894 fly5r = t0r - fly2r;
1895 fly5i = t0i - fly2i;
1896 t0r = t0r + fly2r;
1897 t0i = t0i + fly2i;
1898
1899 fly7r = t1r - fly1i;
1900 fly7i = t1i + fly1r;
1901 t1r = t1r + fly1i;
1902 t1i = t1i - fly1r;
1903
1904 fly2r = fly4r - fly6r;
1905 fly2i = fly4i - fly6i;
1906 fly4r = fly4r + fly6r;
1907 fly4i = fly4i + fly6i;
1908
1909 fly1r = fly0r - fly3i;
1910 fly1i = fly0i + fly3r;
1911 fly0r = fly0r + fly3i;
1912 fly0i = fly0i - fly3r;
1913
1914 fly6r = t0r - fly4r;
1915 fly6i = t0i - fly4i;
1916 t0r = t0r + fly4r;
1917 t0i = t0i + fly4i;
1918
1919 fly3r = fly5r - fly2i;
1920 fly3i = fly5i + fly2r;
1921 fly5r = fly5r + fly2i;
1922 fly5i = fly5i - fly2r;
1923
1924 fly4r = t1r - U3r * fly0r;
1925 fly4r = fly4r - U3i * fly0i;
1926 fly4i = t1i + U3i * fly0r;
1927 fly4i = fly4i - U3r * fly0i;
1928 t1r = scale * t1r - fly4r;
1929 t1i = scale * t1i - fly4i;
1930
1931 fly2r = fly7r + U3i * fly1r;
1932 fly2r = fly2r - U3r * fly1i;
1933 fly2i = fly7i + U3r * fly1r;
1934 fly2i = fly2i + U3i * fly1i;
1935 fly7r = scale * fly7r - fly2r;
1936 fly7i = scale * fly7i - fly2i;
1937
1938 *(fly0P + FlyOffsetA) = fly6r;
1939 *(fly0P + FlyOffsetAIm) = fly6i;
1940 *(fly0P) = t0r;
1941 *(fly0P + 1) = t0i;
1942 *(fly2P + FlyOffsetA) = fly3r;
1943 *(fly2P + FlyOffsetAIm) = fly3i;
1944 *(fly2P) = fly5r;
1945 *(fly2P + 1) = fly5i;
1946 *(fly1P + FlyOffsetA) = fly4r;
1947 *(fly1P + FlyOffsetAIm) = fly4i;
1948 *(fly1P) = t1r;
1949 *(fly1P + 1) = t1i;
1950 *(fly3P + FlyOffsetA) = fly2r;
1951 *(fly3P + FlyOffsetAIm) = fly2i;
1952 *(fly3P) = fly7r;
1953 *(fly3P + 1) = fly7i;
1954
1955 fly0P = (fly0P + FlyOffsetB);
1956 fly1P = (fly1P + FlyOffsetB);
1957 fly2P = (fly2P + FlyOffsetB);
1958 fly3P = (fly3P + FlyOffsetB);
1959
1960 };
1961
1962 NsameU4 >>= 2;
1963 LoopN >>= 2;
1964 NdiffU <<= 2;
1965 Flyinc = Flyinc << 2;
1966 };
1967 };
1968
1969 NsameU2 = NsameU4 >> 1;
1970 NsameU1 = NsameU2 >> 1;
1971 Flyinc <<= 1;
1972 FlyOffsetA = Flyinc << 2;
1973 FlyOffsetB = FlyOffsetA << 1;
1974 FlyOffsetAIm = FlyOffsetA + 1;
1975 FlyOffsetBIm = FlyOffsetB + 1;
1976 LoopN >>= 1;
1977
1978 /* ****** RADIX 8 Stages */
1979 for (stage = stage<<1; stage > 0 ; stage--){
1980
1981 /* an fft stage is done in two parts to ease use of the single quadrant cos table */
1982
1983 /* fftcalc1(iobuf, Utbl, N, NdiffU, LoopN); */
1984 if(!(stage&1)){
1985 U0rP = &Utbl[0];
1986 U0iP = &Utbl[Ntbl[M-1]]; /* for RFFT */
1987 U1rP = U0rP;
1988 U1iP = U0iP;
1989 U2rP = U0rP;
1990 U2iP = U0iP;
1991 U3offset = (Ntbl[M+1]) / 8; /* for RFFT */
1992
1993 IOP = ioptr;
1994
1995 U0r = *U0rP,
1996 U0i = *U0iP;
1997 U1r = *U1rP,
1998 U1i = *U1iP;
1999 U2r = *U2rP,
2000 U2i = *U2iP;
2001 U3r = *( U2rP + U3offset);
2002 U3i = *( U2iP - U3offset);
2003 }
2004
2005 fly0P = IOP;
2006 fly1P = (IOP+Flyinc);
2007 fly2P = (fly1P+Flyinc);
2008 fly3P = (fly2P+Flyinc);
2009
2010 for (diffUcnt = (NdiffU)>>1; diffUcnt > 0; diffUcnt--){
2011
2012 /* Butterflys */
2013 /*
2014 f0 - - t0 - - t0 - - t0
2015 f1 -U0 - t1 - - t1 - - t1
2016 f2 - - f2 -U1 - f5 - - f3
2017 f3 -U0 - f1 -U1a- f7 - - f7
2018 f4 - - f4 - - f4 -U2 - f6
2019 f5 -U0 - f0 - - f0 -U3 - f4
2020 f6 - - f6 -U1 - f2 -U2a- f2
2021 f7 -U0 - f3 -U1a- f1 -U3a- f5
2022 */
2023
2024 fly0r = *(IOP);
2025 fly0i = *(IOP+1);
2026 fly1r = *(fly1P);
2027 fly1i = *(fly1P+1);
2028
2029 for (LoopCnt = LoopN-1; LoopCnt > 0 ; LoopCnt--){
2030
2031 fly2r = *(fly2P);
2032 fly2i = *(fly2P + 1);
2033 fly3r = *(fly3P);
2034 fly3i = *(fly3P + 1);
2035 fly4r = *(fly0P + FlyOffsetA);
2036 fly4i = *(fly0P + FlyOffsetAIm);
2037 fly5r = *(fly1P + FlyOffsetA);
2038 fly5i = *(fly1P + FlyOffsetAIm);
2039 fly6r = *(fly2P + FlyOffsetA);
2040 fly6i = *(fly2P + FlyOffsetAIm);
2041 fly7r = *(fly3P + FlyOffsetA);
2042 fly7i = *(fly3P + FlyOffsetAIm);
2043
2044 t1r = fly0r - U0r * fly1r;
2045 t1r = t1r - U0i * fly1i;
2046 t1i = fly0i + U0i * fly1r;
2047 t1i = t1i - U0r * fly1i;
2048 t0r = scale * fly0r - t1r;
2049 t0i = scale * fly0i - t1i;
2050
2051 fly1r = fly2r - U0r * fly3r;
2052 fly1r = fly1r - U0i * fly3i;
2053 fly1i = fly2i + U0i * fly3r;
2054 fly1i = fly1i - U0r * fly3i;
2055 fly2r = scale * fly2r - fly1r;
2056 fly2i = scale * fly2i - fly1i;
2057
2058 fly0r = fly4r - U0r * fly5r;
2059 fly0r = fly0r - U0i * fly5i;
2060 fly0i = fly4i + U0i * fly5r;
2061 fly0i = fly0i - U0r * fly5i;
2062 fly4r = scale * fly4r - fly0r;
2063 fly4i = scale * fly4i - fly0i;
2064
2065 fly3r = fly6r - U0r * fly7r;
2066 fly3r = fly3r - U0i * fly7i;
2067 fly3i = fly6i + U0i * fly7r;
2068 fly3i = fly3i - U0r * fly7i;
2069 fly6r = scale * fly6r - fly3r;
2070 fly6i = scale * fly6i - fly3i;
2071
2072
2073 fly5r = t0r - U1r * fly2r;
2074 fly5r = fly5r - U1i * fly2i;
2075 fly5i = t0i + U1i * fly2r;
2076 fly5i = fly5i - U1r * fly2i;
2077 t0r = scale * t0r - fly5r;
2078 t0i = scale * t0i - fly5i;
2079
2080 fly7r = t1r + U1i * fly1r;
2081 fly7r = fly7r - U1r * fly1i;
2082 fly7i = t1i + U1r * fly1r;
2083 fly7i = fly7i + U1i * fly1i;
2084 t1r = scale * t1r - fly7r;
2085 t1i = scale * t1i - fly7i;
2086
2087 fly2r = fly4r - U1r * fly6r;
2088 fly2r = fly2r - U1i * fly6i;
2089 fly2i = fly4i + U1i * fly6r;
2090 fly2i = fly2i - U1r * fly6i;
2091 fly4r = scale * fly4r - fly2r;
2092 fly4i = scale * fly4i - fly2i;
2093
2094 fly1r = fly0r + U1i * fly3r;
2095 fly1r = fly1r - U1r * fly3i;
2096 fly1i = fly0i + U1r * fly3r;
2097 fly1i = fly1i + U1i * fly3i;
2098 fly0r = scale * fly0r - fly1r;
2099 fly0i = scale * fly0i - fly1i;
2100
2101 fly6r = t0r - U2r * fly4r;
2102 fly6r = fly6r - U2i * fly4i;
2103 fly6i = t0i + U2i * fly4r;
2104 fly6i = fly6i - U2r * fly4i;
2105 t0r = scale * t0r - fly6r;
2106 t0i = scale * t0i - fly6i;
2107
2108 fly3r = fly5r - U2i * fly2r;
2109 fly3r = fly3r + U2r * fly2i;
2110 fly3i = fly5i - U2r * fly2r;
2111 fly3i = fly3i - U2i * fly2i;
2112 fly2r = scale * fly5r - fly3r;
2113 fly2i = scale * fly5i - fly3i;
2114
2115 fly4r = t1r - U3r * fly0r;
2116 fly4r = fly4r - U3i * fly0i;
2117 fly4i = t1i + U3i * fly0r;
2118 fly4i = fly4i - U3r * fly0i;
2119 t1r = scale * t1r - fly4r;
2120 t1i = scale * t1i - fly4i;
2121
2122 fly5r = fly7r + U3i * fly1r;
2123 fly5r = fly5r - U3r * fly1i;
2124 fly5i = fly7i + U3r * fly1r;
2125 fly5i = fly5i + U3i * fly1i;
2126 fly7r = scale * fly7r - fly5r;
2127 fly7i = scale * fly7i - fly5i;
2128
2129 *(fly0P + FlyOffsetA) = fly6r;
2130 *(fly0P + FlyOffsetAIm) = fly6i;
2131 *(fly0P) = t0r;
2132 *(fly0P + 1) = t0i;
2133 *(fly2P) = fly3r;
2134 *(fly2P + 1) = fly3i;
2135 *(fly2P + FlyOffsetA) = fly2r;
2136 *(fly2P + FlyOffsetAIm) = fly2i;
2137
2138 fly0r = *(fly0P + FlyOffsetB);
2139 fly0i = *(fly0P + FlyOffsetBIm);
2140
2141 *(fly1P + FlyOffsetA) = fly4r;
2142 *(fly1P + FlyOffsetAIm) = fly4i;
2143 *(fly1P) = t1r;
2144 *(fly1P + 1) = t1i;
2145
2146 fly1r = *(fly1P + FlyOffsetB);
2147 fly1i = *(fly1P + FlyOffsetBIm);
2148
2149 *(fly3P + FlyOffsetA) = fly5r;
2150 *(fly3P + FlyOffsetAIm) = fly5i;
2151 *(fly3P) = fly7r;
2152 *(fly3P + 1) = fly7i;
2153
2154 fly0P = (fly0P + FlyOffsetB);
2155 fly1P = (fly1P + FlyOffsetB);
2156 fly2P = (fly2P + FlyOffsetB);
2157 fly3P = (fly3P + FlyOffsetB);
2158 };
2159
2160 fly2r = *(fly2P);
2161 fly2i = *(fly2P + 1);
2162 fly3r = *(fly3P);
2163 fly3i = *(fly3P + 1);
2164 fly4r = *(fly0P + FlyOffsetA);
2165 fly4i = *(fly0P + FlyOffsetAIm);
2166 fly5r = *(fly1P + FlyOffsetA);
2167 fly5i = *(fly1P + FlyOffsetAIm);
2168 fly6r = *(fly2P + FlyOffsetA);
2169 fly6i = *(fly2P + FlyOffsetAIm);
2170 fly7r = *(fly3P + FlyOffsetA);
2171 fly7i = *(fly3P + FlyOffsetAIm);
2172
2173 t1r = fly0r - U0r * fly1r;
2174 t1r = t1r - U0i * fly1i;
2175 t1i = fly0i + U0i * fly1r;
2176 t1i = t1i - U0r * fly1i;
2177 t0r = scale * fly0r - t1r;
2178 t0i = scale * fly0i - t1i;
2179
2180 fly1r = fly2r - U0r * fly3r;
2181 fly1r = fly1r - U0i * fly3i;
2182 fly1i = fly2i + U0i * fly3r;
2183 fly1i = fly1i - U0r * fly3i;
2184 fly2r = scale * fly2r - fly1r;
2185 fly2i = scale * fly2i - fly1i;
2186
2187 fly0r = fly4r - U0r * fly5r;
2188 fly0r = fly0r - U0i * fly5i;
2189 fly0i = fly4i + U0i * fly5r;
2190 fly0i = fly0i - U0r * fly5i;
2191 fly4r = scale * fly4r - fly0r;
2192 fly4i = scale * fly4i - fly0i;
2193
2194 fly3r = fly6r - U0r * fly7r;
2195 fly3r = fly3r - U0i * fly7i;
2196 fly3i = fly6i + U0i * fly7r;
2197 fly3i = fly3i - U0r * fly7i;
2198 fly6r = scale * fly6r - fly3r;
2199 fly6i = scale * fly6i - fly3i;
2200
2201
2202 fly5r = t0r - U1r * fly2r;
2203 fly5r = fly5r - U1i * fly2i;
2204 fly5i = t0i + U1i * fly2r;
2205 fly5i = fly5i - U1r * fly2i;
2206 t0r = scale * t0r - fly5r;
2207 t0i = scale * t0i - fly5i;
2208
2209 fly7r = t1r + U1i * fly1r;
2210 fly7r = fly7r - U1r * fly1i;
2211 fly7i = t1i + U1r * fly1r;
2212 fly7i = fly7i + U1i * fly1i;
2213 t1r = scale * t1r - fly7r;
2214 t1i = scale * t1i - fly7i;
2215
2216 fly2r = fly4r - U1r * fly6r;
2217 fly2r = fly2r - U1i * fly6i;
2218 fly2i = fly4i + U1i * fly6r;
2219 fly2i = fly2i - U1r * fly6i;
2220 fly4r = scale * fly4r - fly2r;
2221 fly4i = scale * fly4i - fly2i;
2222
2223 fly1r = fly0r + U1i * fly3r;
2224 fly1r = fly1r - U1r * fly3i;
2225 fly1i = fly0i + U1r * fly3r;
2226 fly1i = fly1i + U1i * fly3i;
2227 fly0r = scale * fly0r - fly1r;
2228 fly0i = scale * fly0i - fly1i;
2229
2230 U0i = *(U0iP = (U0iP - NsameU4));
2231 U0r = *(U0rP = (U0rP + NsameU4));
2232 U1r = *(U1rP = (U1rP + NsameU2));
2233 U1i = *(U1iP = (U1iP - NsameU2));
2234 if(stage&1)
2235 U0r = -U0r;
2236
2237 fly6r = t0r - U2r * fly4r;
2238 fly6r = fly6r - U2i * fly4i;
2239 fly6i = t0i + U2i * fly4r;
2240 fly6i = fly6i - U2r * fly4i;
2241 t0r = scale * t0r - fly6r;
2242 t0i = scale * t0i - fly6i;
2243
2244 fly3r = fly5r - U2i * fly2r;
2245 fly3r = fly3r + U2r * fly2i;
2246 fly3i = fly5i - U2r * fly2r;
2247 fly3i = fly3i - U2i * fly2i;
2248 fly2r = scale * fly5r - fly3r;
2249 fly2i = scale * fly5i - fly3i;
2250
2251 fly4r = t1r - U3r * fly0r;
2252 fly4r = fly4r - U3i * fly0i;
2253 fly4i = t1i + U3i * fly0r;
2254 fly4i = fly4i - U3r * fly0i;
2255 t1r = scale * t1r - fly4r;
2256 t1i = scale * t1i - fly4i;
2257
2258 fly5r = fly7r + U3i * fly1r;
2259 fly5r = fly5r - U3r * fly1i;
2260 fly5i = fly7i + U3r * fly1r;
2261 fly5i = fly5i + U3i * fly1i;
2262 fly7r = scale * fly7r - fly5r;
2263 fly7i = scale * fly7i - fly5i;
2264
2265 *(fly0P + FlyOffsetA) = fly6r;
2266 *(fly0P + FlyOffsetAIm) = fly6i;
2267 *(fly0P) = t0r;
2268 *(fly0P + 1) = t0i;
2269
2270 U2r = *(U2rP = (U2rP + NsameU1));
2271 U2i = *(U2iP = (U2iP - NsameU1));
2272
2273 *(fly2P) = fly3r;
2274 *(fly2P + 1) = fly3i;
2275 *(fly2P + FlyOffsetA) = fly2r;
2276 *(fly2P + FlyOffsetAIm) = fly2i;
2277 *(fly1P + FlyOffsetA) = fly4r;
2278 *(fly1P + FlyOffsetAIm) = fly4i;
2279 *(fly1P) = t1r;
2280 *(fly1P + 1) = t1i;
2281
2282 U3r = *( U2rP + U3offset);
2283 U3i = *( U2iP - U3offset);
2284
2285 *(fly3P + FlyOffsetA) = fly5r;
2286 *(fly3P + FlyOffsetAIm) = fly5i;
2287 *(fly3P) = fly7r;
2288 *(fly3P + 1) = fly7i;
2289
2290 IOP = IOP + 2;
2291 fly0P = IOP;
2292 fly1P = (IOP+Flyinc);
2293 fly2P = (fly1P+Flyinc);
2294 fly3P = (fly2P+Flyinc);
2295 };
2296 NsameU4 = -NsameU4;
2297
2298 if(stage&1){
2299 LoopN >>= 3;
2300 NsameU1 >>= 3;
2301 NsameU2 >>= 3;
2302 NsameU4 >>= 3;
2303 NdiffU <<= 3;
2304 Flyinc = Flyinc << 3;
2305 FlyOffsetA <<= 3;
2306 FlyOffsetB <<= 3;
2307 FlyOffsetAIm = FlyOffsetA + 1;
2308 FlyOffsetBIm = FlyOffsetB + 1;
2309 }
2310 }
2311
2312 /* Finish RFFT */
2313 M=M+1;
2314
2315 FlyOffsetA = Ntbl[M] * 2/4;
2316 FlyOffsetAIm = Ntbl[M] * 2/4 + 1;
2317
2318 IOP = ioptr;
2319
2320 fly0P = (IOP + Ntbl[M]*2/4);
2321 fly1P = (IOP + Ntbl[M]*2/8);
2322
2323 U0rP = &Utbl[Ntbl[M-3]];
2324
2325 U0r = *U0rP,
2326
2327 fly0r = *(IOP);
2328 fly0i = *(IOP + 1);
2329 fly1r = *(fly0P);
2330 fly1i = *(fly0P + 1);
2331 fly2r = *(fly1P);
2332 fly2i = *(fly1P + 1);
2333 fly3r = *(fly1P + FlyOffsetA);
2334 fly3i = *(fly1P + FlyOffsetAIm);
2335
2336 t0r = scale * fly0r + scale * fly0i; /* compute Re(x[0]) */
2337 t0i = scale * fly0r - scale * fly0i; /* compute Re(x[N/2]) */
2338 t1r = scale * fly1r;
2339 t1i = -scale * fly1i;
2340
2341 fly0r = fly2r + fly3r;
2342 fly0i = fly2i - fly3i;
2343 fly1r = fly2i + fly3i;
2344 fly1i = fly3r - fly2r;
2345
2346 fly2r = fly0r + U0r * fly1r;
2347 fly2r = fly2r + U0r * fly1i;
2348 fly2i = fly0i - U0r * fly1r;
2349 fly2i = fly2i + U0r * fly1i;
2350 fly3r = scale * fly0r - fly2r;
2351 fly3i = fly2i - scale * fly0i;
2352
2353 *(IOP) = t0r;
2354 *(IOP + 1) = t0i;
2355 *(fly0P) = t1r;
2356 *(fly0P + 1) = t1i;
2357 *(fly1P) = fly2r;
2358 *(fly1P + 1) = fly2i;
2359 *(fly1P + FlyOffsetA) = fly3r;
2360 *(fly1P + FlyOffsetAIm) = fly3i;
2361
2362 U0rP = &Utbl[1];
2363 U0iP = &Utbl[Ntbl[M-2]-1];
2364
2365 U0r = *U0rP,
2366 U0i = *U0iP;
2367
2368 fly0P = (IOP + 2);
2369 fly1P = (IOP + (Ntbl[M-2]-1)*2);
2370
2371 /* Butterflys */
2372 /*
2373 f0 - t0 - - f2
2374 f1 - t1 -U0 - f3
2375 f2 - f0 - - t0
2376 f3 - f1 -U0a- t1
2377 */
2378
2379 for (diffUcnt = Ntbl[M-3]-1; diffUcnt > 0 ; diffUcnt--){
2380
2381 fly0r = *(fly0P);
2382 fly0i = *(fly0P + 1);
2383 fly1r = *(fly1P + FlyOffsetA);
2384 fly1i = *(fly1P + FlyOffsetAIm);
2385 fly2r = *(fly1P);
2386 fly2i = *(fly1P + 1);
2387 fly3r = *(fly0P + FlyOffsetA);
2388 fly3i = *(fly0P + FlyOffsetAIm);
2389
2390 t0r = fly0r + fly1r;
2391 t0i = fly0i - fly1i;
2392 t1r = fly0i + fly1i;
2393 t1i = fly1r - fly0r;
2394
2395 fly0r = fly2r + fly3r;
2396 fly0i = fly2i - fly3i;
2397 fly1r = fly2i + fly3i;
2398 fly1i = fly3r - fly2r;
2399
2400 fly2r = t0r + U0r * t1r;
2401 fly2r = fly2r + U0i * t1i;
2402 fly2i = t0i - U0i * t1r;
2403 fly2i = fly2i + U0r * t1i;
2404 fly3r = scale * t0r - fly2r;
2405 fly3i = fly2i - scale * t0i;
2406
2407 t0r = fly0r + U0i * fly1r;
2408 t0r = t0r + U0r * fly1i;
2409 t0i = fly0i - U0r * fly1r;
2410 t0i = t0i + U0i * fly1i;
2411 t1r = scale * fly0r - t0r;
2412 t1i = t0i - scale * fly0i;
2413
2414 *(fly0P) = fly2r;
2415 *(fly0P + 1) = fly2i;
2416 *(fly1P + FlyOffsetA) = fly3r;
2417 *(fly1P + FlyOffsetAIm) = fly3i;
2418
2419 U0r = *(U0rP = (U0rP + 1));
2420 U0i = *(U0iP = (U0iP - 1));
2421
2422 *(fly1P) = t0r;
2423 *(fly1P + 1) = t0i;
2424 *(fly0P + FlyOffsetA) = t1r;
2425 *(fly0P + FlyOffsetAIm) = t1i;
2426
2427 fly0P += 2;
2428 fly1P -= 2;
2429 };
2430
2431 ioptr += Ntbl[M];
2432 }
2433 }
2434
2435 void riffts(float *ioptr, long M, long Rows, float *Utbl){
2436 /* Compute in-place real ifft on the rows of the input array */
2437 /* INPUTS */
2438 /* M = log2 of fft size */
2439 /* *ioptr = input data array in the following order */
2440 /* Re(x[0]), Re(x[N/2]), Re(x[1]), Im(x[1]), Re(x[2]), Im(x[2]), ... Re(x[N/2-1]), Im(x[N/2-1]). */
2441 /* *Utbl = cosine table */
2442 /* Rows = number of rows in ioptr array (use Rows of 1 if ioptr is a 1 dimensional array) */
2443 /* OUTPUTS */
2444 /* *ioptr = real output data array */
2445
2446 long Flyinc;
2447 long FlyOffsetA;
2448 long FlyOffsetAIm;
2449 long FlyOffsetB;
2450 long FlyOffsetBIm;
2451 long NsameU1;
2452 long NsameU2;
2453 long NsameU4;
2454 long diffUcnt;
2455 long LoopCnt;
2456 float scale;
2457 float fly0r;
2458 float fly0i;
2459 float fly1r;
2460 float fly1i;
2461 float fly2r;
2462 float fly2i;
2463 float fly3r;
2464 float fly3i;
2465 float fly4r;
2466 float fly4i;
2467 float fly5r;
2468 float fly5i;
2469 float fly6r;
2470 float fly6i;
2471 float fly7r;
2472 float fly7i;
2473 float U0r;
2474 float U0i;
2475 float U1r;
2476 float U1i;
2477 float U2r;
2478 float U2i;
2479 float U3r;
2480 float U3i;
2481 float t0r;
2482 float t0i;
2483 float t1r;
2484 float t1i;
2485
2486 float *fly0P;
2487 float *fly1P;
2488 float *fly2P;
2489 float *fly3P;
2490
2491 float *U0rP;
2492 float *U0iP;
2493 float *U1rP;
2494 float *U1iP;
2495 float *U2rP;
2496 float *U2iP;
2497 float *IOP;
2498 long U3offset;
2499
2500 long stage;
2501 long NdiffU;
2502 long LoopN;
2503
2504 const long BRshift = MAXMROOT - ((M-1)>>1); /* for RIFFT */
2505 const long Nrems2 = Ntbl[(M-1)-((M-1)>>1)+1]; /* for RIFFT */
2506 const long Nroot_1_ColInc = (Ntbl[(M-1)-1]-Ntbl[(M-1)-((M-1)>>1)])*2; /* for RIFFT */
2507
2508 for (;Rows>0;Rows--){
2509
2510 /* Start RIFFT */
2511
2512 FlyOffsetA = Ntbl[M] * 2/4;
2513 FlyOffsetAIm = Ntbl[M] * 2/4 + 1;
2514
2515 IOP = ioptr;
2516
2517 fly0P = (IOP + Ntbl[M]*2/4);
2518 fly1P = (IOP + Ntbl[M]*2/8);
2519
2520 U0rP = &Utbl[Ntbl[M-3]];
2521
2522 U0r = *U0rP,
2523
2524 fly0r = *(IOP);
2525 fly0i = *(IOP + 1);
2526 fly1r = *(fly0P);
2527 fly1i = *(fly0P + 1);
2528 fly2r = *(fly1P);
2529 fly2i = *(fly1P + 1);
2530 fly3r = *(fly1P + FlyOffsetA);
2531 fly3i = *(fly1P + FlyOffsetAIm);
2532
2533 t0r = fly0r + fly0i;
2534 t0i = fly0r - fly0i;
2535 t1r = fly1r + fly1r;
2536 t1i = -fly1i - fly1i;
2537
2538 fly0r = fly2r + fly3r;
2539 fly0i = fly2i - fly3i;
2540 fly1r = fly2r - fly3r;
2541 fly1i = fly2i + fly3i;
2542
2543 fly3r = fly1r * U0r;
2544 fly3r = fly3r - U0r * fly1i;
2545 fly3i = fly1r * U0r;
2546 fly3i = fly3i + U0r * fly1i;
2547
2548 fly2r = fly0r - fly3i;
2549 fly2i = fly0i + fly3r;
2550 fly1r = fly0r + fly3i;
2551 fly1i = -fly0i + fly3r;
2552
2553 *(IOP) = t0r;
2554 *(IOP + 1) = t0i;
2555 *(fly0P) = t1r;
2556 *(fly0P + 1) = t1i;
2557 *(fly1P) = fly2r;
2558 *(fly1P + 1) = fly2i;
2559 *(fly1P + FlyOffsetA) = fly1r;
2560 *(fly1P + FlyOffsetAIm) = fly1i;
2561
2562 U0rP = &Utbl[1];
2563 U0iP = &Utbl[Ntbl[M-2]-1];
2564
2565 U0r = *U0rP,
2566 U0i = *U0iP;
2567
2568 fly0P = (IOP + 2);
2569 fly1P = (IOP + (Ntbl[M-2]-1)*2);
2570
2571 /* Butterflys */
2572 /*
2573 f0 - t0 - f2
2574 f1 -t1 -U0- f3 - f1
2575 f2 - f0 - t0
2576 f3 -f1 -U0a t1 - f3
2577 */
2578
2579 for (diffUcnt = Ntbl[M-3]-1; diffUcnt > 0 ; diffUcnt--){
2580
2581 fly0r = *(fly0P);
2582 fly0i = *(fly0P + 1);
2583 fly1r = *(fly1P + FlyOffsetA);
2584 fly1i = *(fly1P + FlyOffsetAIm);
2585 fly2r = *(fly1P);
2586 fly2i = *(fly1P + 1);
2587 fly3r = *(fly0P + FlyOffsetA);
2588 fly3i = *(fly0P + FlyOffsetAIm);
2589
2590 t0r = fly0r + fly1r;
2591 t0i = fly0i - fly1i;
2592 t1r = fly0r - fly1r;
2593 t1i = fly0i + fly1i;
2594
2595 fly0r = fly2r + fly3r;
2596 fly0i = fly2i - fly3i;
2597 fly1r = fly2r - fly3r;
2598 fly1i = fly2i + fly3i;
2599
2600
2601 fly3r = t1r * U0r;
2602 fly3r = fly3r - U0i * t1i;
2603 fly3i = t1r * U0i;
2604 fly3i = fly3i + U0r * t1i;
2605
2606 t1r = fly1r * U0i;
2607 t1r = t1r - U0r * fly1i;
2608 t1i = fly1r * U0r;
2609 t1i = t1i + U0i * fly1i;
2610
2611
2612 fly2r = t0r - fly3i;
2613 fly2i = t0i + fly3r;
2614 fly1r = t0r + fly3i;
2615 fly1i = -t0i + fly3r;
2616
2617 t0r = fly0r - t1i;
2618 t0i = fly0i + t1r;
2619 fly3r = fly0r + t1i;
2620 fly3i = -fly0i + t1r;
2621
2622
2623 *(fly0P) = fly2r;
2624 *(fly0P + 1) = fly2i;
2625 *(fly1P + FlyOffsetA) = fly1r;
2626 *(fly1P + FlyOffsetAIm) = fly1i;
2627
2628 U0r = *(U0rP = (U0rP + 1));
2629 U0i = *(U0iP = (U0iP - 1));
2630
2631 *(fly1P) = t0r;
2632 *(fly1P + 1) = t0i;
2633 *(fly0P + FlyOffsetA) = fly3r;
2634 *(fly0P + FlyOffsetAIm) = fly3i;
2635
2636 fly0P += 2;
2637 fly1P -= 2;
2638 };
2639
2640 M=M-1; /* for RIFFT */
2641
2642 FlyOffsetA = Ntbl[M] * 2/2;
2643 FlyOffsetAIm = FlyOffsetA + 1;
2644 FlyOffsetB = FlyOffsetA + 2;
2645 FlyOffsetBIm = FlyOffsetB + 1;
2646
2647 /* BitrevR2 ** bit reverse shuffle and first radix 2 stage ******/
2648
2649 scale = 1./Ntbl[M+1];
2650 for (stage = 0; stage < Ntbl[M-(M>>1)]*2; stage += Ntbl[M>>1]*2){
2651 for (LoopN = (Ntbl[(M>>1)-1]-1); LoopN >= 0; LoopN--){
2652 LoopCnt = (Ntbl[(M>>1)-1]-1);
2653 fly0P = ioptr + Nroot_1_ColInc + ((int)BRcnt[LoopN] >> BRshift)*(2*2) + stage;
2654 IOP = ioptr + (LoopN<<(M+1)/2) * 2 + stage;
2655 fly1P = IOP + ((int)BRcnt[LoopCnt] >> BRshift)*(2*2);
2656 fly0r = *(fly0P);
2657 fly0i = *(fly0P+1);
2658 fly1r = *(fly0P+FlyOffsetA);
2659 fly1i = *(fly0P+FlyOffsetAIm);
2660 for (; LoopCnt > LoopN;){
2661 fly2r = *(fly0P+2);
2662 fly2i = *(fly0P+(2+1));
2663 fly3r = *(fly0P+FlyOffsetB);
2664 fly3i = *(fly0P+FlyOffsetBIm);
2665 fly4r = *(fly1P);
2666 fly4i = *(fly1P+1);
2667 fly5r = *(fly1P+FlyOffsetA);
2668 fly5i = *(fly1P+FlyOffsetAIm);
2669 fly6r = *(fly1P+2);
2670 fly6i = *(fly1P+(2+1));
2671 fly7r = *(fly1P+FlyOffsetB);
2672 fly7i = *(fly1P+FlyOffsetBIm);
2673
2674 t0r = fly0r + fly1r;
2675 t0i = fly0i + fly1i;
2676 fly1r = fly0r - fly1r;
2677 fly1i = fly0i - fly1i;
2678 t1r = fly2r + fly3r;
2679 t1i = fly2i + fly3i;
2680 fly3r = fly2r - fly3r;
2681 fly3i = fly2i - fly3i;
2682 fly0r = fly4r + fly5r;
2683 fly0i = fly4i + fly5i;
2684 fly5r = fly4r - fly5r;
2685 fly5i = fly4i - fly5i;
2686 fly2r = fly6r + fly7r;
2687 fly2i = fly6i + fly7i;
2688 fly7r = fly6r - fly7r;
2689 fly7i = fly6i - fly7i;
2690
2691 *(fly1P) = scale*t0r;
2692 *(fly1P+1) = scale*t0i;
2693 *(fly1P+2) = scale*fly1r;
2694 *(fly1P+(2+1)) = scale*fly1i;
2695 *(fly1P+FlyOffsetA) = scale*t1r;
2696 *(fly1P+FlyOffsetAIm) = scale*t1i;
2697 *(fly1P+FlyOffsetB) = scale*fly3r;
2698 *(fly1P+FlyOffsetBIm) = scale*fly3i;
2699 *(fly0P) = scale*fly0r;
2700 *(fly0P+1) = scale*fly0i;
2701 *(fly0P+2) = scale*fly5r;
2702 *(fly0P+(2+1)) = scale*fly5i;
2703 *(fly0P+FlyOffsetA) = scale*fly2r;
2704 *(fly0P+FlyOffsetAIm) = scale*fly2i;
2705 *(fly0P+FlyOffsetB) = scale*fly7r;
2706 *(fly0P+FlyOffsetBIm) = scale*fly7i;
2707
2708 fly0P -= Nrems2;
2709 fly0r = *(fly0P);
2710 fly0i = *(fly0P+1);
2711 fly1r = *(fly0P+FlyOffsetA);
2712 fly1i = *(fly0P+FlyOffsetAIm);
2713 LoopCnt -= 1;
2714 fly1P = (IOP + ((int)BRcnt[LoopCnt] >> BRshift)*(2*2));
2715 };
2716 fly2r = *(fly0P+2);
2717 fly2i = *(fly0P+(2+1));
2718 fly3r = *(fly0P+FlyOffsetB);
2719 fly3i = *(fly0P+FlyOffsetBIm);
2720
2721 t0r = fly0r + fly1r;
2722 t0i = fly0i + fly1i;
2723 fly1r = fly0r - fly1r;
2724 fly1i = fly0i - fly1i;
2725 t1r = fly2r + fly3r;
2726 t1i = fly2i + fly3i;
2727 fly3r = fly2r - fly3r;
2728 fly3i = fly2i - fly3i;
2729
2730 *(fly0P) = scale*t0r;
2731 *(fly0P+1) = scale*t0i;
2732 *(fly0P+2) = scale*fly1r;
2733 *(fly0P+(2+1)) = scale*fly1i;
2734 *(fly0P+FlyOffsetA) = scale*t1r;
2735 *(fly0P+FlyOffsetAIm) = scale*t1i;
2736 *(fly0P+FlyOffsetB) = scale*fly3r;
2737 *(fly0P+FlyOffsetBIm) = scale*fly3i;
2738
2739 };
2740 };
2741
2742 /**** FFTC **************/
2743
2744 scale = 2.0;
2745
2746 NdiffU = 2;
2747 Flyinc = (NdiffU);
2748
2749 NsameU4 = Ntbl[M+1]/4; /* for RIFFT */
2750 LoopN = Ntbl[M-3];
2751
2752 stage = ((M-1)/3);
2753 if ( (M-1-(stage * 3)) != 0 ){
2754 FlyOffsetA = Flyinc << 2;
2755 FlyOffsetB = FlyOffsetA << 1;
2756 FlyOffsetAIm = FlyOffsetA + 1;
2757 FlyOffsetBIm = FlyOffsetB + 1;
2758 if ( (M-1-(stage * 3)) == 1 ){
2759 /* 1 radix 2 stage */
2760
2761 IOP = ioptr;
2762 fly0P = IOP;
2763 fly1P = (IOP+Flyinc);
2764 fly2P = (fly1P+Flyinc);
2765 fly3P = (fly2P+Flyinc);
2766
2767 /* Butterflys */
2768 /*
2769 t0 - - t0
2770 t1 - - t1
2771 f2 - 1- f5
2772 f1 - -i- f7
2773 f4 - - f4
2774 f0 - - f0
2775 f6 - 1- f2
2776 f3 - -i- f1
2777 */
2778
2779 for (LoopCnt = LoopN; LoopCnt > 0 ; LoopCnt--){
2780 t0r = *(fly0P);
2781 t0i = *(fly0P + 1);
2782 t1r = *(fly1P);
2783 t1i = *(fly1P + 1);
2784 fly2r = *(fly2P);
2785 fly2i = *(fly2P + 1);
2786 fly1r = *(fly3P);
2787 fly1i = *(fly3P + 1);
2788 fly4r = *(fly0P + FlyOffsetA);
2789 fly4i = *(fly0P + FlyOffsetAIm);
2790 fly0r = *(fly1P + FlyOffsetA);
2791 fly0i = *(fly1P + FlyOffsetAIm);
2792 fly6r = *(fly2P + FlyOffsetA);
2793 fly6i = *(fly2P + FlyOffsetAIm);
2794 fly3r = *(fly3P + FlyOffsetA);
2795 fly3i = *(fly3P + FlyOffsetAIm);
2796
2797 fly5r = t0r - fly2r;
2798 fly5i = t0i - fly2i;
2799 t0r = t0r + fly2r;
2800 t0i = t0i + fly2i;
2801
2802 fly7r = t1r + fly1i;
2803 fly7i = t1i - fly1r;
2804 t1r = t1r - fly1i;
2805 t1i = t1i + fly1r;
2806
2807 fly2r = fly4r - fly6r;
2808 fly2i = fly4i - fly6i;
2809 fly4r = fly4r + fly6r;
2810 fly4i = fly4i + fly6i;
2811
2812 fly1r = fly0r + fly3i;
2813 fly1i = fly0i - fly3r;
2814 fly0r = fly0r - fly3i;
2815 fly0i = fly0i + fly3r;
2816
2817 *(fly2P) = fly5r;
2818 *(fly2P + 1) = fly5i;
2819 *(fly0P) = t0r;
2820 *(fly0P + 1) = t0i;
2821 *(fly3P) = fly7r;
2822 *(fly3P + 1) = fly7i;
2823 *(fly1P) = t1r;
2824 *(fly1P + 1) = t1i;
2825 *(fly2P + FlyOffsetA) = fly2r;
2826 *(fly2P + FlyOffsetAIm) = fly2i;
2827 *(fly0P + FlyOffsetA) = fly4r;
2828 *(fly0P + FlyOffsetAIm) = fly4i;
2829 *(fly3P + FlyOffsetA) = fly1r;
2830 *(fly3P + FlyOffsetAIm) = fly1i;
2831 *(fly1P + FlyOffsetA) = fly0r;
2832 *(fly1P + FlyOffsetAIm) = fly0i;
2833
2834 fly0P = (fly0P + FlyOffsetB);
2835 fly1P = (fly1P + FlyOffsetB);
2836 fly2P = (fly2P + FlyOffsetB);
2837 fly3P = (fly3P + FlyOffsetB);
2838 };
2839
2840 NsameU4 >>= 1;
2841 LoopN >>= 1;
2842 NdiffU <<= 1;
2843 Flyinc = Flyinc << 1;
2844 }
2845 else{
2846 /* 1 radix 4 stage */
2847 IOP = ioptr;
2848
2849 U3r = 0.7071067811865475244008443621; /* sqrt(0.5); */
2850 U3i = U3r;
2851 fly0P = IOP;
2852 fly1P = (IOP+Flyinc);
2853 fly2P = (fly1P+Flyinc);
2854 fly3P = (fly2P+Flyinc);
2855
2856 /* Butterflys */
2857 /*
2858 t0 - - t0 - - t0
2859 t1 - - t1 - - t1
2860 f2 - 1- f5 - - f5
2861 f1 - -i- f7 - - f7
2862 f4 - - f4 - 1- f6
2863 f0 - - f0 -U3 - f3
2864 f6 - 1- f2 - -i- f4
2865 f3 - -i- f1 -U3a- f2
2866 */
2867
2868 for (LoopCnt = LoopN; LoopCnt > 0 ; LoopCnt--){
2869 t0r = *(fly0P);
2870 t0i = *(fly0P + 1);
2871 t1r = *(fly1P);
2872 t1i = *(fly1P + 1);
2873 fly2r = *(fly2P);
2874 fly2i = *(fly2P + 1);
2875 fly1r = *(fly3P);
2876 fly1i = *(fly3P + 1);
2877 fly4r = *(fly0P + FlyOffsetA);
2878 fly4i = *(fly0P + FlyOffsetAIm);
2879 fly0r = *(fly1P + FlyOffsetA);
2880 fly0i = *(fly1P + FlyOffsetAIm);
2881 fly6r = *(fly2P + FlyOffsetA);
2882 fly6i = *(fly2P + FlyOffsetAIm);
2883 fly3r = *(fly3P + FlyOffsetA);
2884 fly3i = *(fly3P + FlyOffsetAIm);
2885
2886 fly5r = t0r - fly2r;
2887 fly5i = t0i - fly2i;
2888 t0r = t0r + fly2r;
2889 t0i = t0i + fly2i;
2890
2891 fly7r = t1r + fly1i;
2892 fly7i = t1i - fly1r;
2893 t1r = t1r - fly1i;
2894 t1i = t1i + fly1r;
2895
2896 fly2r = fly4r - fly6r;
2897 fly2i = fly4i - fly6i;
2898 fly4r = fly4r + fly6r;
2899 fly4i = fly4i + fly6i;
2900
2901 fly1r = fly0r + fly3i;
2902 fly1i = fly0i - fly3r;
2903 fly0r = fly0r - fly3i;
2904 fly0i = fly0i + fly3r;
2905
2906 fly6r = t0r - fly4r;
2907 fly6i = t0i - fly4i;
2908 t0r = t0r + fly4r;
2909 t0i = t0i + fly4i;
2910
2911 fly3r = fly5r + fly2i;
2912 fly3i = fly5i - fly2r;
2913 fly5r = fly5r - fly2i;
2914 fly5i = fly5i + fly2r;
2915
2916 fly4r = t1r - U3r * fly0r;
2917 fly4r = fly4r + U3i * fly0i;
2918 fly4i = t1i - U3i * fly0r;
2919 fly4i = fly4i - U3r * fly0i;
2920 t1r = scale * t1r - fly4r;
2921 t1i = scale * t1i - fly4i;
2922
2923 fly2r = fly7r + U3i * fly1r;
2924 fly2r = fly2r + U3r * fly1i;
2925 fly2i = fly7i - U3r * fly1r;
2926 fly2i = fly2i + U3i * fly1i;
2927 fly7r = scale * fly7r - fly2r;
2928 fly7i = scale * fly7i - fly2i;
2929
2930 *(fly0P + FlyOffsetA) = fly6r;
2931 *(fly0P + FlyOffsetAIm) = fly6i;
2932 *(fly0P) = t0r;
2933 *(fly0P + 1) = t0i;
2934 *(fly2P + FlyOffsetA) = fly3r;
2935 *(fly2P + FlyOffsetAIm) = fly3i;
2936 *(fly2P) = fly5r;
2937 *(fly2P + 1) = fly5i;
2938 *(fly1P + FlyOffsetA) = fly4r;
2939 *(fly1P + FlyOffsetAIm) = fly4i;
2940 *(fly1P) = t1r;
2941 *(fly1P + 1) = t1i;
2942 *(fly3P + FlyOffsetA) = fly2r;
2943 *(fly3P + FlyOffsetAIm) = fly2i;
2944 *(fly3P) = fly7r;
2945 *(fly3P + 1) = fly7i;
2946
2947 fly0P = (fly0P + FlyOffsetB);
2948 fly1P = (fly1P + FlyOffsetB);
2949 fly2P = (fly2P + FlyOffsetB);
2950 fly3P = (fly3P + FlyOffsetB);
2951
2952 };
2953
2954 NsameU4 >>= 2;
2955 LoopN >>= 2;
2956 NdiffU <<= 2;
2957 Flyinc = Flyinc << 2;
2958 };
2959 };
2960
2961 NsameU2 = NsameU4 >> 1;
2962 NsameU1 = NsameU2 >> 1;
2963 Flyinc <<= 1;
2964 FlyOffsetA = Flyinc << 2;
2965 FlyOffsetB = FlyOffsetA << 1;
2966 FlyOffsetAIm = FlyOffsetA + 1;
2967 FlyOffsetBIm = FlyOffsetB + 1;
2968 LoopN >>= 1;
2969
2970 /* ****** RADIX 8 Stages */
2971 for (stage = stage<<1; stage > 0 ; stage--){
2972
2973 /* an fft stage is done in two parts to ease use of the single quadrant cos table */
2974
2975 /* fftcalc1(iobuf, Utbl, N, NdiffU, LoopN); */
2976 if(!(stage&1)){
2977 U0rP = &Utbl[0];
2978 U0iP = &Utbl[Ntbl[M-1]]; /* for RIFFT */
2979 U1rP = U0rP;
2980 U1iP = U0iP;
2981 U2rP = U0rP;
2982 U2iP = U0iP;
2983 U3offset = (Ntbl[M+1]) / 8; /* for RIFFT */
2984
2985 IOP = ioptr;
2986
2987 U0r = *U0rP,
2988 U0i = *U0iP;
2989 U1r = *U1rP,
2990 U1i = *U1iP;
2991 U2r = *U2rP,
2992 U2i = *U2iP;
2993 U3r = *( U2rP + U3offset);
2994 U3i = *( U2iP - U3offset);
2995 }
2996
2997 fly0P = IOP;
2998 fly1P = (IOP+Flyinc);
2999 fly2P = (fly1P+Flyinc);
3000 fly3P = (fly2P+Flyinc);
3001
3002 for (diffUcnt = (NdiffU)>>1; diffUcnt > 0; diffUcnt--){
3003
3004 /* Butterflys */
3005 /*
3006 f0 - - t0 - - t0 - - t0
3007 f1 -U0 - t1 - - t1 - - t1
3008 f2 - - f2 -U1 - f5 - - f3
3009 f3 -U0 - f1 -U1a- f7 - - f7
3010 f4 - - f4 - - f4 -U2 - f6
3011 f5 -U0 - f0 - - f0 -U3 - f4
3012 f6 - - f6 -U1 - f2 -U2a- f2
3013 f7 -U0 - f3 -U1a- f1 -U3a- f5
3014 */
3015
3016 fly0r = *(IOP);
3017 fly0i = *(IOP+1);
3018 fly1r = *(fly1P);
3019 fly1i = *(fly1P+1);
3020
3021 for (LoopCnt = LoopN-1; LoopCnt > 0 ; LoopCnt--){
3022
3023 fly2r = *(fly2P);
3024 fly2i = *(fly2P + 1);
3025 fly3r = *(fly3P);
3026 fly3i = *(fly3P + 1);
3027 fly4r = *(fly0P + FlyOffsetA);
3028 fly4i = *(fly0P + FlyOffsetAIm);
3029 fly5r = *(fly1P + FlyOffsetA);
3030 fly5i = *(fly1P + FlyOffsetAIm);
3031 fly6r = *(fly2P + FlyOffsetA);
3032 fly6i = *(fly2P + FlyOffsetAIm);
3033 fly7r = *(fly3P + FlyOffsetA);
3034 fly7i = *(fly3P + FlyOffsetAIm);
3035
3036 t1r = fly0r - U0r * fly1r;
3037 t1r = t1r + U0i * fly1i;
3038 t1i = fly0i - U0i * fly1r;
3039 t1i = t1i - U0r * fly1i;
3040 t0r = scale * fly0r - t1r;
3041 t0i = scale * fly0i - t1i;
3042
3043 fly1r = fly2r - U0r * fly3r;
3044 fly1r = fly1r + U0i * fly3i;
3045 fly1i = fly2i - U0i * fly3r;
3046 fly1i = fly1i - U0r * fly3i;
3047 fly2r = scale * fly2r - fly1r;
3048 fly2i = scale * fly2i - fly1i;
3049
3050 fly0r = fly4r - U0r * fly5r;
3051 fly0r = fly0r + U0i * fly5i;
3052 fly0i = fly4i - U0i * fly5r;
3053 fly0i = fly0i - U0r * fly5i;
3054 fly4r = scale * fly4r - fly0r;
3055 fly4i = scale * fly4i - fly0i;
3056
3057 fly3r = fly6r - U0r * fly7r;
3058 fly3r = fly3r + U0i * fly7i;
3059 fly3i = fly6i - U0i * fly7r;
3060 fly3i = fly3i - U0r * fly7i;
3061 fly6r = scale * fly6r - fly3r;
3062 fly6i = scale * fly6i - fly3i;
3063
3064
3065 fly5r = t0r - U1r * fly2r;
3066 fly5r = fly5r + U1i * fly2i;
3067 fly5i = t0i - U1i * fly2r;
3068 fly5i = fly5i - U1r * fly2i;
3069 t0r = scale * t0r - fly5r;
3070 t0i = scale * t0i - fly5i;
3071
3072 fly7r = t1r + U1i * fly1r;
3073 fly7r = fly7r + U1r * fly1i;
3074 fly7i = t1i - U1r * fly1r;
3075 fly7i = fly7i + U1i * fly1i;
3076 t1r = scale * t1r - fly7r;
3077 t1i = scale * t1i - fly7i;
3078
3079 fly2r = fly4r - U1r * fly6r;
3080 fly2r = fly2r + U1i * fly6i;
3081 fly2i = fly4i - U1i * fly6r;
3082 fly2i = fly2i - U1r * fly6i;
3083 fly4r = scale * fly4r - fly2r;
3084 fly4i = scale * fly4i - fly2i;
3085
3086 fly1r = fly0r + U1i * fly3r;
3087 fly1r = fly1r + U1r * fly3i;
3088 fly1i = fly0i - U1r * fly3r;
3089 fly1i = fly1i + U1i * fly3i;
3090 fly0r = scale * fly0r - fly1r;
3091 fly0i = scale * fly0i - fly1i;
3092
3093 fly6r = t0r - U2r * fly4r;
3094 fly6r = fly6r + U2i * fly4i;
3095 fly6i = t0i - U2i * fly4r;
3096 fly6i = fly6i - U2r * fly4i;
3097 t0r = scale * t0r - fly6r;
3098 t0i = scale * t0i - fly6i;
3099
3100 fly3r = fly5r - U2i * fly2r;
3101 fly3r = fly3r - U2r * fly2i;
3102 fly3i = fly5i + U2r * fly2r;
3103 fly3i = fly3i - U2i * fly2i;
3104 fly2r = scale * fly5r - fly3r;
3105 fly2i = scale * fly5i - fly3i;
3106
3107 fly4r = t1r - U3r * fly0r;
3108 fly4r = fly4r + U3i * fly0i;
3109 fly4i = t1i - U3i * fly0r;
3110 fly4i = fly4i - U3r * fly0i;
3111 t1r = scale * t1r - fly4r;
3112 t1i = scale * t1i - fly4i;
3113
3114 fly5r = fly7r + U3i * fly1r;
3115 fly5r = fly5r + U3r * fly1i;
3116 fly5i = fly7i - U3r * fly1r;
3117 fly5i = fly5i + U3i * fly1i;
3118 fly7r = scale * fly7r - fly5r;
3119 fly7i = scale * fly7i - fly5i;
3120
3121 *(fly0P + FlyOffsetA) = fly6r;
3122 *(fly0P + FlyOffsetAIm) = fly6i;
3123 *(fly0P) = t0r;
3124 *(fly0P + 1) = t0i;
3125 *(fly2P) = fly3r;
3126 *(fly2P + 1) = fly3i;
3127 *(fly2P + FlyOffsetA) = fly2r;
3128 *(fly2P + FlyOffsetAIm) = fly2i;
3129
3130 fly0r = *(fly0P + FlyOffsetB);
3131 fly0i = *(fly0P + FlyOffsetBIm);
3132
3133 *(fly1P + FlyOffsetA) = fly4r;
3134 *(fly1P + FlyOffsetAIm) = fly4i;
3135 *(fly1P) = t1r;
3136 *(fly1P + 1) = t1i;
3137
3138 fly1r = *(fly1P + FlyOffsetB);
3139 fly1i = *(fly1P + FlyOffsetBIm);
3140
3141 *(fly3P + FlyOffsetA) = fly5r;
3142 *(fly3P + FlyOffsetAIm) = fly5i;
3143 *(fly3P) = fly7r;
3144 *(fly3P + 1) = fly7i;
3145
3146 fly0P = (fly0P + FlyOffsetB);
3147 fly1P = (fly1P + FlyOffsetB);
3148 fly2P = (fly2P + FlyOffsetB);
3149 fly3P = (fly3P + FlyOffsetB);
3150
3151 };
3152 fly2r = *(fly2P);
3153 fly2i = *(fly2P + 1);
3154 fly3r = *(fly3P);
3155 fly3i = *(fly3P + 1);
3156 fly4r = *(fly0P + FlyOffsetA);
3157 fly4i = *(fly0P + FlyOffsetAIm);
3158 fly5r = *(fly1P + FlyOffsetA);
3159 fly5i = *(fly1P + FlyOffsetAIm);
3160 fly6r = *(fly2P + FlyOffsetA);
3161 fly6i = *(fly2P + FlyOffsetAIm);
3162 fly7r = *(fly3P + FlyOffsetA);
3163 fly7i = *(fly3P + FlyOffsetAIm);
3164
3165 t1r = fly0r - U0r * fly1r;
3166 t1r = t1r + U0i * fly1i;
3167 t1i = fly0i - U0i * fly1r;
3168 t1i = t1i - U0r * fly1i;
3169 t0r = scale * fly0r - t1r;
3170 t0i = scale * fly0i - t1i;
3171
3172 fly1r = fly2r - U0r * fly3r;
3173 fly1r = fly1r + U0i * fly3i;
3174 fly1i = fly2i - U0i * fly3r;
3175 fly1i = fly1i - U0r * fly3i;
3176 fly2r = scale * fly2r - fly1r;
3177 fly2i = scale * fly2i - fly1i;
3178
3179 fly0r = fly4r - U0r * fly5r;
3180 fly0r = fly0r + U0i * fly5i;
3181 fly0i = fly4i - U0i * fly5r;
3182 fly0i = fly0i - U0r * fly5i;
3183 fly4r = scale * fly4r - fly0r;
3184 fly4i = scale * fly4i - fly0i;
3185
3186 fly3r = fly6r - U0r * fly7r;
3187 fly3r = fly3r + U0i * fly7i;
3188 fly3i = fly6i - U0i * fly7r;
3189 fly3i = fly3i - U0r * fly7i;
3190 fly6r = scale * fly6r - fly3r;
3191 fly6i = scale * fly6i - fly3i;
3192
3193 fly5r = t0r - U1r * fly2r;
3194 fly5r = fly5r + U1i * fly2i;
3195 fly5i = t0i - U1i * fly2r;
3196 fly5i = fly5i - U1r * fly2i;
3197 t0r = scale * t0r - fly5r;
3198 t0i = scale * t0i - fly5i;
3199
3200 fly7r = t1r + U1i * fly1r;
3201 fly7r = fly7r + U1r * fly1i;
3202 fly7i = t1i - U1r * fly1r;
3203 fly7i = fly7i + U1i * fly1i;
3204 t1r = scale * t1r - fly7r;
3205 t1i = scale * t1i - fly7i;
3206
3207 fly2r = fly4r - U1r * fly6r;
3208 fly2r = fly2r + U1i * fly6i;
3209 fly2i = fly4i - U1i * fly6r;
3210 fly2i = fly2i - U1r * fly6i;
3211 fly4r = scale * fly4r - fly2r;
3212 fly4i = scale * fly4i - fly2i;
3213
3214 fly1r = fly0r + U1i * fly3r;
3215 fly1r = fly1r + U1r * fly3i;
3216 fly1i = fly0i - U1r * fly3r;
3217 fly1i = fly1i + U1i * fly3i;
3218 fly0r = scale * fly0r - fly1r;
3219 fly0i = scale * fly0i - fly1i;
3220
3221 U0i = *(U0iP = (U0iP - NsameU4));
3222 U0r = *(U0rP = (U0rP + NsameU4));
3223 U1r = *(U1rP = (U1rP + NsameU2));
3224 U1i = *(U1iP = (U1iP - NsameU2));
3225 if(stage&1)
3226 U0r = - U0r;
3227
3228 fly6r = t0r - U2r * fly4r;
3229 fly6r = fly6r + U2i * fly4i;
3230 fly6i = t0i - U2i * fly4r;
3231 fly6i = fly6i - U2r * fly4i;
3232 t0r = scale * t0r - fly6r;
3233 t0i = scale * t0i - fly6i;
3234
3235 fly3r = fly5r - U2i * fly2r;
3236 fly3r = fly3r - U2r * fly2i;
3237 fly3i = fly5i + U2r * fly2r;
3238 fly3i = fly3i - U2i * fly2i;
3239 fly2r = scale * fly5r - fly3r;
3240 fly2i = scale * fly5i - fly3i;
3241
3242 fly4r = t1r - U3r * fly0r;
3243 fly4r = fly4r + U3i * fly0i;
3244 fly4i = t1i - U3i * fly0r;
3245 fly4i = fly4i - U3r * fly0i;
3246 t1r = scale * t1r - fly4r;
3247 t1i = scale * t1i - fly4i;
3248
3249 fly5r = fly7r + U3i * fly1r;
3250 fly5r = fly5r + U3r * fly1i;
3251 fly5i = fly7i - U3r * fly1r;
3252 fly5i = fly5i + U3i * fly1i;
3253 fly7r = scale * fly7r - fly5r;
3254 fly7i = scale * fly7i - fly5i;
3255
3256 *(fly0P + FlyOffsetA) = fly6r;
3257 *(fly0P + FlyOffsetAIm) = fly6i;
3258 *(fly0P) = t0r;
3259 *(fly0P + 1) = t0i;
3260
3261 U2r = *(U2rP = (U2rP + NsameU1));
3262 U2i = *(U2iP = (U2iP - NsameU1));
3263
3264 *(fly2P) = fly3r;
3265 *(fly2P + 1) = fly3i;
3266 *(fly2P + FlyOffsetA) = fly2r;
3267 *(fly2P + FlyOffsetAIm) = fly2i;
3268 *(fly1P + FlyOffsetA) = fly4r;
3269 *(fly1P + FlyOffsetAIm) = fly4i;
3270 *(fly1P) = t1r;
3271 *(fly1P + 1) = t1i;
3272
3273 U3r = *( U2rP + U3offset);
3274 U3i = *( U2iP - U3offset);
3275
3276 *(fly3P + FlyOffsetA) = fly5r;
3277 *(fly3P + FlyOffsetAIm) = fly5i;
3278 *(fly3P) = fly7r;
3279 *(fly3P + 1) = fly7i;
3280
3281 IOP = IOP + 2;
3282 fly0P = IOP;
3283 fly1P = (IOP+Flyinc);
3284 fly2P = (fly1P+Flyinc);
3285 fly3P = (fly2P+Flyinc);
3286 };
3287 NsameU4 = -NsameU4;
3288
3289 if(stage&1){
3290 LoopN >>= 3;
3291 NsameU1 >>= 3;
3292 NsameU2 >>= 3;
3293 NsameU4 >>= 3;
3294 NdiffU <<= 3;
3295 Flyinc = Flyinc << 3;
3296 FlyOffsetA <<= 3;
3297 FlyOffsetB <<= 3;
3298 FlyOffsetAIm = FlyOffsetA + 1;
3299 FlyOffsetBIm = FlyOffsetB + 1;
3300 }
3301 }
3302 M=M+1; /* for RIFFT */
3303
3304 ioptr += Ntbl[M];
3305 }
3306 }