Line data Source code
1 : /*
2 : * Copyright(c)1995,97 Mark Olesen <olesen@me.QueensU.CA>
3 : * Queen's Univ at Kingston (Canada)
4 : *
5 : * Permission to use, copy, modify, and distribute this software for
6 : * any purpose without fee is hereby granted, provided that this
7 : * entire notice is included in all copies of any software which is
8 : * or includes a copy or modification of this software and in all
9 : * copies of the supporting documentation for such software.
10 : *
11 : * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR
12 : * IMPLIED WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR QUEEN'S
13 : * UNIVERSITY AT KINGSTON MAKES ANY REPRESENTATION OR WARRANTY OF ANY
14 : * KIND CONCERNING THE MERCHANTABILITY OF THIS SOFTWARE OR ITS
15 : * FITNESS FOR ANY PARTICULAR PURPOSE.
16 : *
17 : * All of which is to say that you can do what you like with this
18 : * source code provided you don't try to sell it as your own and you
19 : * include an unaltered copy of this message (including the
20 : * copyright).
21 : *
22 : * It is also implicitly understood that bug fixes and improvements
23 : * should make their way back to the general Internet community so
24 : * that everyone benefits.
25 : *
26 : * Changes:
27 : * Trivial type modifications by the WebRTC authors.
28 : */
29 :
30 :
31 : /*
32 : * File:
33 : * WebRtcIsac_Fftn.c
34 : *
35 : * Public:
36 : * WebRtcIsac_Fftn / fftnf ();
37 : *
38 : * Private:
39 : * WebRtcIsac_Fftradix / fftradixf ();
40 : *
41 : * Descript:
42 : * multivariate complex Fourier transform, computed in place
43 : * using mixed-radix Fast Fourier Transform algorithm.
44 : *
45 : * Fortran code by:
46 : * RC Singleton, Stanford Research Institute, Sept. 1968
47 : *
48 : * translated by f2c (version 19950721).
49 : *
50 : * int WebRtcIsac_Fftn (int ndim, const int dims[], REAL Re[], REAL Im[],
51 : * int iSign, double scaling);
52 : *
53 : * NDIM = the total number dimensions
54 : * DIMS = a vector of array sizes
55 : * if NDIM is zero then DIMS must be zero-terminated
56 : *
57 : * RE and IM hold the real and imaginary components of the data, and return
58 : * the resulting real and imaginary Fourier coefficients. Multidimensional
59 : * data *must* be allocated contiguously. There is no limit on the number
60 : * of dimensions.
61 : *
62 : * ISIGN = the sign of the complex exponential (ie, forward or inverse FFT)
63 : * the magnitude of ISIGN (normally 1) is used to determine the
64 : * correct indexing increment (see below).
65 : *
66 : * SCALING = normalizing constant by which the final result is *divided*
67 : * if SCALING == -1, normalize by total dimension of the transform
68 : * if SCALING < -1, normalize by the square-root of the total dimension
69 : *
70 : * example:
71 : * tri-variate transform with Re[n1][n2][n3], Im[n1][n2][n3]
72 : *
73 : * int dims[3] = {n1,n2,n3}
74 : * WebRtcIsac_Fftn (3, dims, Re, Im, 1, scaling);
75 : *
76 : *-----------------------------------------------------------------------*
77 : * int WebRtcIsac_Fftradix (REAL Re[], REAL Im[], size_t nTotal, size_t nPass,
78 : * size_t nSpan, int iSign, size_t max_factors,
79 : * size_t max_perm);
80 : *
81 : * RE, IM - see above documentation
82 : *
83 : * Although there is no limit on the number of dimensions, WebRtcIsac_Fftradix() must
84 : * be called once for each dimension, but the calls may be in any order.
85 : *
86 : * NTOTAL = the total number of complex data values
87 : * NPASS = the dimension of the current variable
88 : * NSPAN/NPASS = the spacing of consecutive data values while indexing the
89 : * current variable
90 : * ISIGN - see above documentation
91 : *
92 : * example:
93 : * tri-variate transform with Re[n1][n2][n3], Im[n1][n2][n3]
94 : *
95 : * WebRtcIsac_Fftradix (Re, Im, n1*n2*n3, n1, n1, 1, maxf, maxp);
96 : * WebRtcIsac_Fftradix (Re, Im, n1*n2*n3, n2, n1*n2, 1, maxf, maxp);
97 : * WebRtcIsac_Fftradix (Re, Im, n1*n2*n3, n3, n1*n2*n3, 1, maxf, maxp);
98 : *
99 : * single-variate transform,
100 : * NTOTAL = N = NSPAN = (number of complex data values),
101 : *
102 : * WebRtcIsac_Fftradix (Re, Im, n, n, n, 1, maxf, maxp);
103 : *
104 : * The data can also be stored in a single array with alternating real and
105 : * imaginary parts, the magnitude of ISIGN is changed to 2 to give correct
106 : * indexing increment, and data [0] and data [1] used to pass the initial
107 : * addresses for the sequences of real and imaginary values,
108 : *
109 : * example:
110 : * REAL data [2*NTOTAL];
111 : * WebRtcIsac_Fftradix ( &data[0], &data[1], NTOTAL, nPass, nSpan, 2, maxf, maxp);
112 : *
113 : * for temporary allocation:
114 : *
115 : * MAX_FACTORS >= the maximum prime factor of NPASS
116 : * MAX_PERM >= the number of prime factors of NPASS. In addition,
117 : * if the square-free portion K of NPASS has two or more prime
118 : * factors, then MAX_PERM >= (K-1)
119 : *
120 : * storage in FACTOR for a maximum of 15 prime factors of NPASS. if NPASS
121 : * has more than one square-free factor, the product of the square-free
122 : * factors must be <= 210 array storage for maximum prime factor of 23 the
123 : * following two constants should agree with the array dimensions.
124 : *
125 : *----------------------------------------------------------------------*/
126 : #include "fft.h"
127 :
128 : #include <stdlib.h>
129 : #include <math.h>
130 :
131 :
132 :
133 : /* double precision routine */
134 : static int
135 : WebRtcIsac_Fftradix (double Re[], double Im[],
136 : size_t nTotal, size_t nPass, size_t nSpan, int isign,
137 : int max_factors, unsigned int max_perm,
138 : FFTstr *fftstate);
139 :
140 :
141 :
142 : #ifndef M_PI
143 : # define M_PI 3.14159265358979323846264338327950288
144 : #endif
145 :
146 : #ifndef SIN60
147 : # define SIN60 0.86602540378443865 /* sin(60 deg) */
148 : # define COS72 0.30901699437494742 /* cos(72 deg) */
149 : # define SIN72 0.95105651629515357 /* sin(72 deg) */
150 : #endif
151 :
152 : # define REAL double
153 : # define FFTN WebRtcIsac_Fftn
154 : # define FFTNS "fftn"
155 : # define FFTRADIX WebRtcIsac_Fftradix
156 : # define FFTRADIXS "fftradix"
157 :
158 :
159 0 : int WebRtcIsac_Fftns(unsigned int ndim, const int dims[],
160 : double Re[],
161 : double Im[],
162 : int iSign,
163 : double scaling,
164 : FFTstr *fftstate)
165 : {
166 :
167 : size_t nSpan, nPass, nTotal;
168 : unsigned int i;
169 : int ret, max_factors, max_perm;
170 :
171 : /*
172 : * tally the number of elements in the data array
173 : * and determine the number of dimensions
174 : */
175 0 : nTotal = 1;
176 0 : if (ndim && dims [0])
177 : {
178 0 : for (i = 0; i < ndim; i++)
179 : {
180 0 : if (dims [i] <= 0)
181 : {
182 0 : return -1;
183 : }
184 0 : nTotal *= dims [i];
185 : }
186 : }
187 : else
188 : {
189 0 : ndim = 0;
190 0 : for (i = 0; dims [i]; i++)
191 : {
192 0 : if (dims [i] <= 0)
193 : {
194 0 : return -1;
195 : }
196 0 : nTotal *= dims [i];
197 0 : ndim++;
198 : }
199 : }
200 :
201 : /* determine maximum number of factors and permuations */
202 : #if 1
203 : /*
204 : * follow John Beale's example, just use the largest dimension and don't
205 : * worry about excess allocation. May be someone else will do it?
206 : */
207 0 : max_factors = max_perm = 1;
208 0 : for (i = 0; i < ndim; i++)
209 : {
210 0 : nSpan = dims [i];
211 0 : if ((int)nSpan > max_factors)
212 : {
213 0 : max_factors = (int)nSpan;
214 : }
215 0 : if ((int)nSpan > max_perm)
216 : {
217 0 : max_perm = (int)nSpan;
218 : }
219 : }
220 : #else
221 : /* use the constants used in the original Fortran code */
222 : max_factors = 23;
223 : max_perm = 209;
224 : #endif
225 : /* loop over the dimensions: */
226 0 : nPass = 1;
227 0 : for (i = 0; i < ndim; i++)
228 : {
229 0 : nSpan = dims [i];
230 0 : nPass *= nSpan;
231 0 : ret = FFTRADIX (Re, Im, nTotal, nSpan, nPass, iSign,
232 : max_factors, max_perm, fftstate);
233 : /* exit, clean-up already done */
234 0 : if (ret)
235 0 : return ret;
236 : }
237 :
238 : /* Divide through by the normalizing constant: */
239 0 : if (scaling && scaling != 1.0)
240 : {
241 0 : if (iSign < 0) iSign = -iSign;
242 0 : if (scaling < 0.0)
243 : {
244 0 : scaling = (double)nTotal;
245 0 : if (scaling < -1.0)
246 0 : scaling = sqrt (scaling);
247 : }
248 0 : scaling = 1.0 / scaling; /* multiply is often faster */
249 0 : for (i = 0; i < nTotal; i += iSign)
250 : {
251 0 : Re [i] *= scaling;
252 0 : Im [i] *= scaling;
253 : }
254 : }
255 0 : return 0;
256 : }
257 :
258 : /*
259 : * singleton's mixed radix routine
260 : *
261 : * could move allocation out to WebRtcIsac_Fftn(), but leave it here so that it's
262 : * possible to make this a standalone function
263 : */
264 :
265 0 : static int FFTRADIX (REAL Re[],
266 : REAL Im[],
267 : size_t nTotal,
268 : size_t nPass,
269 : size_t nSpan,
270 : int iSign,
271 : int max_factors,
272 : unsigned int max_perm,
273 : FFTstr *fftstate)
274 : {
275 : int ii, mfactor, kspan, ispan, inc;
276 : int j, jc, jf, jj, k, k1, k2, k3, k4, kk, kt, nn, ns, nt;
277 :
278 :
279 : REAL radf;
280 : REAL c1, c2, c3, cd, aa, aj, ak, ajm, ajp, akm, akp;
281 : REAL s1, s2, s3, sd, bb, bj, bk, bjm, bjp, bkm, bkp;
282 :
283 0 : REAL *Rtmp = NULL; /* temp space for real part*/
284 0 : REAL *Itmp = NULL; /* temp space for imaginary part */
285 0 : REAL *Cos = NULL; /* Cosine values */
286 0 : REAL *Sin = NULL; /* Sine values */
287 :
288 0 : REAL s60 = SIN60; /* sin(60 deg) */
289 0 : REAL c72 = COS72; /* cos(72 deg) */
290 0 : REAL s72 = SIN72; /* sin(72 deg) */
291 0 : REAL pi2 = M_PI; /* use PI first, 2 PI later */
292 :
293 :
294 0 : fftstate->SpaceAlloced = 0;
295 0 : fftstate->MaxPermAlloced = 0;
296 :
297 :
298 : // initialize to avoid warnings
299 0 : k3 = c2 = c3 = s2 = s3 = 0.0;
300 :
301 0 : if (nPass < 2)
302 0 : return 0;
303 :
304 : /* allocate storage */
305 0 : if (fftstate->SpaceAlloced < max_factors * sizeof (REAL))
306 : {
307 : #ifdef SUN_BROKEN_REALLOC
308 : if (!fftstate->SpaceAlloced) /* first time */
309 : {
310 : fftstate->SpaceAlloced = max_factors * sizeof (REAL);
311 : }
312 : else
313 : {
314 : #endif
315 0 : fftstate->SpaceAlloced = max_factors * sizeof (REAL);
316 : #ifdef SUN_BROKEN_REALLOC
317 : }
318 : #endif
319 : }
320 : else
321 : {
322 : /* allow full use of alloc'd space */
323 0 : max_factors = fftstate->SpaceAlloced / sizeof (REAL);
324 : }
325 0 : if (fftstate->MaxPermAlloced < max_perm)
326 : {
327 : #ifdef SUN_BROKEN_REALLOC
328 : if (!fftstate->MaxPermAlloced) /* first time */
329 : else
330 : #endif
331 0 : fftstate->MaxPermAlloced = max_perm;
332 : }
333 : else
334 : {
335 : /* allow full use of alloc'd space */
336 0 : max_perm = fftstate->MaxPermAlloced;
337 : }
338 :
339 : /* assign pointers */
340 0 : Rtmp = (REAL *) fftstate->Tmp0;
341 0 : Itmp = (REAL *) fftstate->Tmp1;
342 0 : Cos = (REAL *) fftstate->Tmp2;
343 0 : Sin = (REAL *) fftstate->Tmp3;
344 :
345 : /*
346 : * Function Body
347 : */
348 0 : inc = iSign;
349 0 : if (iSign < 0) {
350 0 : s72 = -s72;
351 0 : s60 = -s60;
352 0 : pi2 = -pi2;
353 0 : inc = -inc; /* absolute value */
354 : }
355 :
356 : /* adjust for strange increments */
357 0 : nt = inc * (int)nTotal;
358 0 : ns = inc * (int)nSpan;
359 0 : kspan = ns;
360 :
361 0 : nn = nt - inc;
362 0 : jc = ns / (int)nPass;
363 0 : radf = pi2 * (double) jc;
364 0 : pi2 *= 2.0; /* use 2 PI from here on */
365 :
366 0 : ii = 0;
367 0 : jf = 0;
368 : /* determine the factors of n */
369 0 : mfactor = 0;
370 0 : k = (int)nPass;
371 0 : while (k % 16 == 0) {
372 0 : mfactor++;
373 0 : fftstate->factor [mfactor - 1] = 4;
374 0 : k /= 16;
375 : }
376 0 : j = 3;
377 0 : jj = 9;
378 : do {
379 0 : while (k % jj == 0) {
380 0 : mfactor++;
381 0 : fftstate->factor [mfactor - 1] = j;
382 0 : k /= jj;
383 : }
384 0 : j += 2;
385 0 : jj = j * j;
386 0 : } while (jj <= k);
387 0 : if (k <= 4) {
388 0 : kt = mfactor;
389 0 : fftstate->factor [mfactor] = k;
390 0 : if (k != 1)
391 0 : mfactor++;
392 : } else {
393 0 : if (k - (k / 4 << 2) == 0) {
394 0 : mfactor++;
395 0 : fftstate->factor [mfactor - 1] = 2;
396 0 : k /= 4;
397 : }
398 0 : kt = mfactor;
399 0 : j = 2;
400 : do {
401 0 : if (k % j == 0) {
402 0 : mfactor++;
403 0 : fftstate->factor [mfactor - 1] = j;
404 0 : k /= j;
405 : }
406 0 : j = ((j + 1) / 2 << 1) + 1;
407 0 : } while (j <= k);
408 : }
409 0 : if (kt) {
410 0 : j = kt;
411 : do {
412 0 : mfactor++;
413 0 : fftstate->factor [mfactor - 1] = fftstate->factor [j - 1];
414 0 : j--;
415 0 : } while (j);
416 : }
417 :
418 : /* test that mfactors is in range */
419 0 : if (mfactor > NFACTOR)
420 : {
421 0 : return -1;
422 : }
423 :
424 : /* compute fourier transform */
425 : for (;;) {
426 0 : sd = radf / (double) kspan;
427 0 : cd = sin(sd);
428 0 : cd = 2.0 * cd * cd;
429 0 : sd = sin(sd + sd);
430 0 : kk = 0;
431 0 : ii++;
432 :
433 0 : switch (fftstate->factor [ii - 1]) {
434 : case 2:
435 : /* transform for factor of 2 (including rotation factor) */
436 0 : kspan /= 2;
437 0 : k1 = kspan + 2;
438 : do {
439 : do {
440 0 : k2 = kk + kspan;
441 0 : ak = Re [k2];
442 0 : bk = Im [k2];
443 0 : Re [k2] = Re [kk] - ak;
444 0 : Im [k2] = Im [kk] - bk;
445 0 : Re [kk] += ak;
446 0 : Im [kk] += bk;
447 0 : kk = k2 + kspan;
448 0 : } while (kk < nn);
449 0 : kk -= nn;
450 0 : } while (kk < jc);
451 0 : if (kk >= kspan)
452 0 : goto Permute_Results_Label; /* exit infinite loop */
453 : do {
454 0 : c1 = 1.0 - cd;
455 0 : s1 = sd;
456 : do {
457 : do {
458 : do {
459 0 : k2 = kk + kspan;
460 0 : ak = Re [kk] - Re [k2];
461 0 : bk = Im [kk] - Im [k2];
462 0 : Re [kk] += Re [k2];
463 0 : Im [kk] += Im [k2];
464 0 : Re [k2] = c1 * ak - s1 * bk;
465 0 : Im [k2] = s1 * ak + c1 * bk;
466 0 : kk = k2 + kspan;
467 0 : } while (kk < (nt-1));
468 0 : k2 = kk - nt;
469 0 : c1 = -c1;
470 0 : kk = k1 - k2;
471 0 : } while (kk > k2);
472 0 : ak = c1 - (cd * c1 + sd * s1);
473 0 : s1 = sd * c1 - cd * s1 + s1;
474 0 : c1 = 2.0 - (ak * ak + s1 * s1);
475 0 : s1 *= c1;
476 0 : c1 *= ak;
477 0 : kk += jc;
478 0 : } while (kk < k2);
479 0 : k1 += inc + inc;
480 0 : kk = (k1 - kspan + 1) / 2 + jc - 1;
481 0 : } while (kk < (jc + jc));
482 0 : break;
483 :
484 : case 4: /* transform for factor of 4 */
485 0 : ispan = kspan;
486 0 : kspan /= 4;
487 :
488 : do {
489 0 : c1 = 1.0;
490 0 : s1 = 0.0;
491 : do {
492 : do {
493 0 : k1 = kk + kspan;
494 0 : k2 = k1 + kspan;
495 0 : k3 = k2 + kspan;
496 0 : akp = Re [kk] + Re [k2];
497 0 : akm = Re [kk] - Re [k2];
498 0 : ajp = Re [k1] + Re [k3];
499 0 : ajm = Re [k1] - Re [k3];
500 0 : bkp = Im [kk] + Im [k2];
501 0 : bkm = Im [kk] - Im [k2];
502 0 : bjp = Im [k1] + Im [k3];
503 0 : bjm = Im [k1] - Im [k3];
504 0 : Re [kk] = akp + ajp;
505 0 : Im [kk] = bkp + bjp;
506 0 : ajp = akp - ajp;
507 0 : bjp = bkp - bjp;
508 0 : if (iSign < 0) {
509 0 : akp = akm + bjm;
510 0 : bkp = bkm - ajm;
511 0 : akm -= bjm;
512 0 : bkm += ajm;
513 : } else {
514 0 : akp = akm - bjm;
515 0 : bkp = bkm + ajm;
516 0 : akm += bjm;
517 0 : bkm -= ajm;
518 : }
519 : /* avoid useless multiplies */
520 0 : if (s1 == 0.0) {
521 0 : Re [k1] = akp;
522 0 : Re [k2] = ajp;
523 0 : Re [k3] = akm;
524 0 : Im [k1] = bkp;
525 0 : Im [k2] = bjp;
526 0 : Im [k3] = bkm;
527 : } else {
528 0 : Re [k1] = akp * c1 - bkp * s1;
529 0 : Re [k2] = ajp * c2 - bjp * s2;
530 0 : Re [k3] = akm * c3 - bkm * s3;
531 0 : Im [k1] = akp * s1 + bkp * c1;
532 0 : Im [k2] = ajp * s2 + bjp * c2;
533 0 : Im [k3] = akm * s3 + bkm * c3;
534 : }
535 0 : kk = k3 + kspan;
536 0 : } while (kk < nt);
537 :
538 0 : c2 = c1 - (cd * c1 + sd * s1);
539 0 : s1 = sd * c1 - cd * s1 + s1;
540 0 : c1 = 2.0 - (c2 * c2 + s1 * s1);
541 0 : s1 *= c1;
542 0 : c1 *= c2;
543 : /* values of c2, c3, s2, s3 that will get used next time */
544 0 : c2 = c1 * c1 - s1 * s1;
545 0 : s2 = 2.0 * c1 * s1;
546 0 : c3 = c2 * c1 - s2 * s1;
547 0 : s3 = c2 * s1 + s2 * c1;
548 0 : kk = kk - nt + jc;
549 0 : } while (kk < kspan);
550 0 : kk = kk - kspan + inc;
551 0 : } while (kk < jc);
552 0 : if (kspan == jc)
553 0 : goto Permute_Results_Label; /* exit infinite loop */
554 0 : break;
555 :
556 : default:
557 : /* transform for odd factors */
558 : #ifdef FFT_RADIX4
559 : return -1;
560 : break;
561 : #else /* FFT_RADIX4 */
562 0 : k = fftstate->factor [ii - 1];
563 0 : ispan = kspan;
564 0 : kspan /= k;
565 :
566 0 : switch (k) {
567 : case 3: /* transform for factor of 3 (optional code) */
568 : do {
569 : do {
570 0 : k1 = kk + kspan;
571 0 : k2 = k1 + kspan;
572 0 : ak = Re [kk];
573 0 : bk = Im [kk];
574 0 : aj = Re [k1] + Re [k2];
575 0 : bj = Im [k1] + Im [k2];
576 0 : Re [kk] = ak + aj;
577 0 : Im [kk] = bk + bj;
578 0 : ak -= 0.5 * aj;
579 0 : bk -= 0.5 * bj;
580 0 : aj = (Re [k1] - Re [k2]) * s60;
581 0 : bj = (Im [k1] - Im [k2]) * s60;
582 0 : Re [k1] = ak - bj;
583 0 : Re [k2] = ak + bj;
584 0 : Im [k1] = bk + aj;
585 0 : Im [k2] = bk - aj;
586 0 : kk = k2 + kspan;
587 0 : } while (kk < (nn - 1));
588 0 : kk -= nn;
589 0 : } while (kk < kspan);
590 0 : break;
591 :
592 : case 5: /* transform for factor of 5 (optional code) */
593 0 : c2 = c72 * c72 - s72 * s72;
594 0 : s2 = 2.0 * c72 * s72;
595 : do {
596 : do {
597 0 : k1 = kk + kspan;
598 0 : k2 = k1 + kspan;
599 0 : k3 = k2 + kspan;
600 0 : k4 = k3 + kspan;
601 0 : akp = Re [k1] + Re [k4];
602 0 : akm = Re [k1] - Re [k4];
603 0 : bkp = Im [k1] + Im [k4];
604 0 : bkm = Im [k1] - Im [k4];
605 0 : ajp = Re [k2] + Re [k3];
606 0 : ajm = Re [k2] - Re [k3];
607 0 : bjp = Im [k2] + Im [k3];
608 0 : bjm = Im [k2] - Im [k3];
609 0 : aa = Re [kk];
610 0 : bb = Im [kk];
611 0 : Re [kk] = aa + akp + ajp;
612 0 : Im [kk] = bb + bkp + bjp;
613 0 : ak = akp * c72 + ajp * c2 + aa;
614 0 : bk = bkp * c72 + bjp * c2 + bb;
615 0 : aj = akm * s72 + ajm * s2;
616 0 : bj = bkm * s72 + bjm * s2;
617 0 : Re [k1] = ak - bj;
618 0 : Re [k4] = ak + bj;
619 0 : Im [k1] = bk + aj;
620 0 : Im [k4] = bk - aj;
621 0 : ak = akp * c2 + ajp * c72 + aa;
622 0 : bk = bkp * c2 + bjp * c72 + bb;
623 0 : aj = akm * s2 - ajm * s72;
624 0 : bj = bkm * s2 - bjm * s72;
625 0 : Re [k2] = ak - bj;
626 0 : Re [k3] = ak + bj;
627 0 : Im [k2] = bk + aj;
628 0 : Im [k3] = bk - aj;
629 0 : kk = k4 + kspan;
630 0 : } while (kk < (nn-1));
631 0 : kk -= nn;
632 0 : } while (kk < kspan);
633 0 : break;
634 :
635 : default:
636 0 : if (k != jf) {
637 0 : jf = k;
638 0 : s1 = pi2 / (double) k;
639 0 : c1 = cos(s1);
640 0 : s1 = sin(s1);
641 0 : if (jf > max_factors){
642 0 : return -1;
643 : }
644 0 : Cos [jf - 1] = 1.0;
645 0 : Sin [jf - 1] = 0.0;
646 0 : j = 1;
647 : do {
648 0 : Cos [j - 1] = Cos [k - 1] * c1 + Sin [k - 1] * s1;
649 0 : Sin [j - 1] = Cos [k - 1] * s1 - Sin [k - 1] * c1;
650 0 : k--;
651 0 : Cos [k - 1] = Cos [j - 1];
652 0 : Sin [k - 1] = -Sin [j - 1];
653 0 : j++;
654 0 : } while (j < k);
655 : }
656 : do {
657 : do {
658 0 : k1 = kk;
659 0 : k2 = kk + ispan;
660 0 : ak = aa = Re [kk];
661 0 : bk = bb = Im [kk];
662 0 : j = 1;
663 0 : k1 += kspan;
664 : do {
665 0 : k2 -= kspan;
666 0 : j++;
667 0 : Rtmp [j - 1] = Re [k1] + Re [k2];
668 0 : ak += Rtmp [j - 1];
669 0 : Itmp [j - 1] = Im [k1] + Im [k2];
670 0 : bk += Itmp [j - 1];
671 0 : j++;
672 0 : Rtmp [j - 1] = Re [k1] - Re [k2];
673 0 : Itmp [j - 1] = Im [k1] - Im [k2];
674 0 : k1 += kspan;
675 0 : } while (k1 < k2);
676 0 : Re [kk] = ak;
677 0 : Im [kk] = bk;
678 0 : k1 = kk;
679 0 : k2 = kk + ispan;
680 0 : j = 1;
681 : do {
682 0 : k1 += kspan;
683 0 : k2 -= kspan;
684 0 : jj = j;
685 0 : ak = aa;
686 0 : bk = bb;
687 0 : aj = 0.0;
688 0 : bj = 0.0;
689 0 : k = 1;
690 : do {
691 0 : k++;
692 0 : ak += Rtmp [k - 1] * Cos [jj - 1];
693 0 : bk += Itmp [k - 1] * Cos [jj - 1];
694 0 : k++;
695 0 : aj += Rtmp [k - 1] * Sin [jj - 1];
696 0 : bj += Itmp [k - 1] * Sin [jj - 1];
697 0 : jj += j;
698 0 : if (jj > jf) {
699 0 : jj -= jf;
700 : }
701 0 : } while (k < jf);
702 0 : k = jf - j;
703 0 : Re [k1] = ak - bj;
704 0 : Im [k1] = bk + aj;
705 0 : Re [k2] = ak + bj;
706 0 : Im [k2] = bk - aj;
707 0 : j++;
708 0 : } while (j < k);
709 0 : kk += ispan;
710 0 : } while (kk < nn);
711 0 : kk -= nn;
712 0 : } while (kk < kspan);
713 0 : break;
714 : }
715 :
716 : /* multiply by rotation factor (except for factors of 2 and 4) */
717 0 : if (ii == mfactor)
718 0 : goto Permute_Results_Label; /* exit infinite loop */
719 0 : kk = jc;
720 : do {
721 0 : c2 = 1.0 - cd;
722 0 : s1 = sd;
723 : do {
724 0 : c1 = c2;
725 0 : s2 = s1;
726 0 : kk += kspan;
727 : do {
728 : do {
729 0 : ak = Re [kk];
730 0 : Re [kk] = c2 * ak - s2 * Im [kk];
731 0 : Im [kk] = s2 * ak + c2 * Im [kk];
732 0 : kk += ispan;
733 0 : } while (kk < nt);
734 0 : ak = s1 * s2;
735 0 : s2 = s1 * c2 + c1 * s2;
736 0 : c2 = c1 * c2 - ak;
737 0 : kk = kk - nt + kspan;
738 0 : } while (kk < ispan);
739 0 : c2 = c1 - (cd * c1 + sd * s1);
740 0 : s1 += sd * c1 - cd * s1;
741 0 : c1 = 2.0 - (c2 * c2 + s1 * s1);
742 0 : s1 *= c1;
743 0 : c2 *= c1;
744 0 : kk = kk - ispan + jc;
745 0 : } while (kk < kspan);
746 0 : kk = kk - kspan + jc + inc;
747 0 : } while (kk < (jc + jc));
748 0 : break;
749 : #endif /* FFT_RADIX4 */
750 : }
751 : }
752 :
753 : /* permute the results to normal order---done in two stages */
754 : /* permutation for square factors of n */
755 : Permute_Results_Label:
756 0 : fftstate->Perm [0] = ns;
757 0 : if (kt) {
758 0 : k = kt + kt + 1;
759 0 : if (mfactor < k)
760 0 : k--;
761 0 : j = 1;
762 0 : fftstate->Perm [k] = jc;
763 : do {
764 0 : fftstate->Perm [j] = fftstate->Perm [j - 1] / fftstate->factor [j - 1];
765 0 : fftstate->Perm [k - 1] = fftstate->Perm [k] * fftstate->factor [j - 1];
766 0 : j++;
767 0 : k--;
768 0 : } while (j < k);
769 0 : k3 = fftstate->Perm [k];
770 0 : kspan = fftstate->Perm [1];
771 0 : kk = jc;
772 0 : k2 = kspan;
773 0 : j = 1;
774 0 : if (nPass != nTotal) {
775 : /* permutation for multivariate transform */
776 : Permute_Multi_Label:
777 : do {
778 : do {
779 0 : k = kk + jc;
780 : do {
781 : /* swap Re [kk] <> Re [k2], Im [kk] <> Im [k2] */
782 0 : ak = Re [kk]; Re [kk] = Re [k2]; Re [k2] = ak;
783 0 : bk = Im [kk]; Im [kk] = Im [k2]; Im [k2] = bk;
784 0 : kk += inc;
785 0 : k2 += inc;
786 0 : } while (kk < (k-1));
787 0 : kk += ns - jc;
788 0 : k2 += ns - jc;
789 0 : } while (kk < (nt-1));
790 0 : k2 = k2 - nt + kspan;
791 0 : kk = kk - nt + jc;
792 0 : } while (k2 < (ns-1));
793 : do {
794 : do {
795 0 : k2 -= fftstate->Perm [j - 1];
796 0 : j++;
797 0 : k2 = fftstate->Perm [j] + k2;
798 0 : } while (k2 > fftstate->Perm [j - 1]);
799 0 : j = 1;
800 : do {
801 0 : if (kk < (k2-1))
802 0 : goto Permute_Multi_Label;
803 0 : kk += jc;
804 0 : k2 += kspan;
805 0 : } while (k2 < (ns-1));
806 0 : } while (kk < (ns-1));
807 : } else {
808 : /* permutation for single-variate transform (optional code) */
809 : Permute_Single_Label:
810 : do {
811 : /* swap Re [kk] <> Re [k2], Im [kk] <> Im [k2] */
812 0 : ak = Re [kk]; Re [kk] = Re [k2]; Re [k2] = ak;
813 0 : bk = Im [kk]; Im [kk] = Im [k2]; Im [k2] = bk;
814 0 : kk += inc;
815 0 : k2 += kspan;
816 0 : } while (k2 < (ns-1));
817 : do {
818 : do {
819 0 : k2 -= fftstate->Perm [j - 1];
820 0 : j++;
821 0 : k2 = fftstate->Perm [j] + k2;
822 0 : } while (k2 >= fftstate->Perm [j - 1]);
823 0 : j = 1;
824 : do {
825 0 : if (kk < k2)
826 0 : goto Permute_Single_Label;
827 0 : kk += inc;
828 0 : k2 += kspan;
829 0 : } while (k2 < (ns-1));
830 0 : } while (kk < (ns-1));
831 : }
832 0 : jc = k3;
833 : }
834 :
835 0 : if ((kt << 1) + 1 >= mfactor)
836 0 : return 0;
837 0 : ispan = fftstate->Perm [kt];
838 : /* permutation for square-free factors of n */
839 0 : j = mfactor - kt;
840 0 : fftstate->factor [j] = 1;
841 : do {
842 0 : fftstate->factor [j - 1] *= fftstate->factor [j];
843 0 : j--;
844 0 : } while (j != kt);
845 0 : kt++;
846 0 : nn = fftstate->factor [kt - 1] - 1;
847 0 : if (nn > (int) max_perm) {
848 0 : return -1;
849 : }
850 0 : j = jj = 0;
851 : for (;;) {
852 0 : k = kt + 1;
853 0 : k2 = fftstate->factor [kt - 1];
854 0 : kk = fftstate->factor [k - 1];
855 0 : j++;
856 0 : if (j > nn)
857 0 : break; /* exit infinite loop */
858 0 : jj += kk;
859 0 : while (jj >= k2) {
860 0 : jj -= k2;
861 0 : k2 = kk;
862 0 : k++;
863 0 : kk = fftstate->factor [k - 1];
864 0 : jj += kk;
865 : }
866 0 : fftstate->Perm [j - 1] = jj;
867 : }
868 : /* determine the permutation cycles of length greater than 1 */
869 0 : j = 0;
870 0 : for (;;) {
871 : do {
872 0 : j++;
873 0 : kk = fftstate->Perm [j - 1];
874 0 : } while (kk < 0);
875 0 : if (kk != j) {
876 : do {
877 0 : k = kk;
878 0 : kk = fftstate->Perm [k - 1];
879 0 : fftstate->Perm [k - 1] = -kk;
880 0 : } while (kk != j);
881 0 : k3 = kk;
882 : } else {
883 0 : fftstate->Perm [j - 1] = -j;
884 0 : if (j == nn)
885 0 : break; /* exit infinite loop */
886 : }
887 : }
888 0 : max_factors *= inc;
889 : /* reorder a and b, following the permutation cycles */
890 : for (;;) {
891 0 : j = k3 + 1;
892 0 : nt -= ispan;
893 0 : ii = nt - inc + 1;
894 0 : if (nt < 0)
895 0 : break; /* exit infinite loop */
896 : do {
897 : do {
898 0 : j--;
899 0 : } while (fftstate->Perm [j - 1] < 0);
900 0 : jj = jc;
901 : do {
902 0 : kspan = jj;
903 0 : if (jj > max_factors) {
904 0 : kspan = max_factors;
905 : }
906 0 : jj -= kspan;
907 0 : k = fftstate->Perm [j - 1];
908 0 : kk = jc * k + ii + jj;
909 0 : k1 = kk + kspan - 1;
910 0 : k2 = 0;
911 : do {
912 0 : k2++;
913 0 : Rtmp [k2 - 1] = Re [k1];
914 0 : Itmp [k2 - 1] = Im [k1];
915 0 : k1 -= inc;
916 0 : } while (k1 != (kk-1));
917 : do {
918 0 : k1 = kk + kspan - 1;
919 0 : k2 = k1 - jc * (k + fftstate->Perm [k - 1]);
920 0 : k = -fftstate->Perm [k - 1];
921 : do {
922 0 : Re [k1] = Re [k2];
923 0 : Im [k1] = Im [k2];
924 0 : k1 -= inc;
925 0 : k2 -= inc;
926 0 : } while (k1 != (kk-1));
927 0 : kk = k2 + 1;
928 0 : } while (k != j);
929 0 : k1 = kk + kspan - 1;
930 0 : k2 = 0;
931 : do {
932 0 : k2++;
933 0 : Re [k1] = Rtmp [k2 - 1];
934 0 : Im [k1] = Itmp [k2 - 1];
935 0 : k1 -= inc;
936 0 : } while (k1 != (kk-1));
937 0 : } while (jj);
938 0 : } while (j != 1);
939 : }
940 0 : return 0; /* exit point here */
941 : }
942 : /* ---------------------- end-of-file (c source) ---------------------- */
943 :
|
