1 /* 2 * Copyright(c)1995,97 Mark Olesen <olesen (at) 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 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 nTotal = 1; 176 if (ndim && dims [0]) 177 { 178 for (i = 0; i < ndim; i++) 179 { 180 if (dims [i] <= 0) 181 { 182 return -1; 183 } 184 nTotal *= dims [i]; 185 } 186 } 187 else 188 { 189 ndim = 0; 190 for (i = 0; dims [i]; i++) 191 { 192 if (dims [i] <= 0) 193 { 194 return -1; 195 } 196 nTotal *= dims [i]; 197 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 max_factors = max_perm = 1; 208 for (i = 0; i < ndim; i++) 209 { 210 nSpan = dims [i]; 211 if ((int)nSpan > max_factors) 212 { 213 max_factors = (int)nSpan; 214 } 215 if ((int)nSpan > max_perm) 216 { 217 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 nPass = 1; 227 for (i = 0; i < ndim; i++) 228 { 229 nSpan = dims [i]; 230 nPass *= nSpan; 231 ret = FFTRADIX (Re, Im, nTotal, nSpan, nPass, iSign, 232 max_factors, max_perm, fftstate); 233 /* exit, clean-up already done */ 234 if (ret) 235 return ret; 236 } 237 238 /* Divide through by the normalizing constant: */ 239 if (scaling && scaling != 1.0) 240 { 241 if (iSign < 0) iSign = -iSign; 242 if (scaling < 0.0) 243 { 244 scaling = (double)nTotal; 245 if (scaling < -1.0) 246 scaling = sqrt (scaling); 247 } 248 scaling = 1.0 / scaling; /* multiply is often faster */ 249 for (i = 0; i < nTotal; i += iSign) 250 { 251 Re [i] *= scaling; 252 Im [i] *= scaling; 253 } 254 } 255 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 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 REAL *Rtmp = NULL; /* temp space for real part*/ 284 REAL *Itmp = NULL; /* temp space for imaginary part */ 285 REAL *Cos = NULL; /* Cosine values */ 286 REAL *Sin = NULL; /* Sine values */ 287 288 REAL s60 = SIN60; /* sin(60 deg) */ 289 REAL c72 = COS72; /* cos(72 deg) */ 290 REAL s72 = SIN72; /* sin(72 deg) */ 291 REAL pi2 = M_PI; /* use PI first, 2 PI later */ 292 293 294 fftstate->SpaceAlloced = 0; 295 fftstate->MaxPermAlloced = 0; 296 297 298 // initialize to avoid warnings 299 k3 = c2 = c3 = s2 = s3 = 0.0; 300 301 if (nPass < 2) 302 return 0; 303 304 /* allocate storage */ 305 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 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 max_factors = fftstate->SpaceAlloced / sizeof (REAL); 324 } 325 if (fftstate->MaxPermAlloced < max_perm) 326 { 327 #ifdef SUN_BROKEN_REALLOC 328 if (!fftstate->MaxPermAlloced) /* first time */ 329 else 330 #endif 331 fftstate->MaxPermAlloced = max_perm; 332 } 333 else 334 { 335 /* allow full use of alloc'd space */ 336 max_perm = fftstate->MaxPermAlloced; 337 } 338 if (fftstate->Tmp0 == NULL || fftstate->Tmp1 == NULL || fftstate->Tmp2 == NULL || fftstate->Tmp3 == NULL 339 || fftstate->Perm == NULL) { 340 return -1; 341 } 342 343 /* assign pointers */ 344 Rtmp = (REAL *) fftstate->Tmp0; 345 Itmp = (REAL *) fftstate->Tmp1; 346 Cos = (REAL *) fftstate->Tmp2; 347 Sin = (REAL *) fftstate->Tmp3; 348 349 /* 350 * Function Body 351 */ 352 inc = iSign; 353 if (iSign < 0) { 354 s72 = -s72; 355 s60 = -s60; 356 pi2 = -pi2; 357 inc = -inc; /* absolute value */ 358 } 359 360 /* adjust for strange increments */ 361 nt = inc * (int)nTotal; 362 ns = inc * (int)nSpan; 363 kspan = ns; 364 365 nn = nt - inc; 366 jc = ns / (int)nPass; 367 radf = pi2 * (double) jc; 368 pi2 *= 2.0; /* use 2 PI from here on */ 369 370 ii = 0; 371 jf = 0; 372 /* determine the factors of n */ 373 mfactor = 0; 374 k = (int)nPass; 375 while (k % 16 == 0) { 376 mfactor++; 377 fftstate->factor [mfactor - 1] = 4; 378 k /= 16; 379 } 380 j = 3; 381 jj = 9; 382 do { 383 while (k % jj == 0) { 384 mfactor++; 385 fftstate->factor [mfactor - 1] = j; 386 k /= jj; 387 } 388 j += 2; 389 jj = j * j; 390 } while (jj <= k); 391 if (k <= 4) { 392 kt = mfactor; 393 fftstate->factor [mfactor] = k; 394 if (k != 1) 395 mfactor++; 396 } else { 397 if (k - (k / 4 << 2) == 0) { 398 mfactor++; 399 fftstate->factor [mfactor - 1] = 2; 400 k /= 4; 401 } 402 kt = mfactor; 403 j = 2; 404 do { 405 if (k % j == 0) { 406 mfactor++; 407 fftstate->factor [mfactor - 1] = j; 408 k /= j; 409 } 410 j = ((j + 1) / 2 << 1) + 1; 411 } while (j <= k); 412 } 413 if (kt) { 414 j = kt; 415 do { 416 mfactor++; 417 fftstate->factor [mfactor - 1] = fftstate->factor [j - 1]; 418 j--; 419 } while (j); 420 } 421 422 /* test that mfactors is in range */ 423 if (mfactor > NFACTOR) 424 { 425 return -1; 426 } 427 428 /* compute fourier transform */ 429 for (;;) { 430 sd = radf / (double) kspan; 431 cd = sin(sd); 432 cd = 2.0 * cd * cd; 433 sd = sin(sd + sd); 434 kk = 0; 435 ii++; 436 437 switch (fftstate->factor [ii - 1]) { 438 case 2: 439 /* transform for factor of 2 (including rotation factor) */ 440 kspan /= 2; 441 k1 = kspan + 2; 442 do { 443 do { 444 k2 = kk + kspan; 445 ak = Re [k2]; 446 bk = Im [k2]; 447 Re [k2] = Re [kk] - ak; 448 Im [k2] = Im [kk] - bk; 449 Re [kk] += ak; 450 Im [kk] += bk; 451 kk = k2 + kspan; 452 } while (kk < nn); 453 kk -= nn; 454 } while (kk < jc); 455 if (kk >= kspan) 456 goto Permute_Results_Label; /* exit infinite loop */ 457 do { 458 c1 = 1.0 - cd; 459 s1 = sd; 460 do { 461 do { 462 do { 463 k2 = kk + kspan; 464 ak = Re [kk] - Re [k2]; 465 bk = Im [kk] - Im [k2]; 466 Re [kk] += Re [k2]; 467 Im [kk] += Im [k2]; 468 Re [k2] = c1 * ak - s1 * bk; 469 Im [k2] = s1 * ak + c1 * bk; 470 kk = k2 + kspan; 471 } while (kk < (nt-1)); 472 k2 = kk - nt; 473 c1 = -c1; 474 kk = k1 - k2; 475 } while (kk > k2); 476 ak = c1 - (cd * c1 + sd * s1); 477 s1 = sd * c1 - cd * s1 + s1; 478 c1 = 2.0 - (ak * ak + s1 * s1); 479 s1 *= c1; 480 c1 *= ak; 481 kk += jc; 482 } while (kk < k2); 483 k1 += inc + inc; 484 kk = (k1 - kspan + 1) / 2 + jc - 1; 485 } while (kk < (jc + jc)); 486 break; 487 488 case 4: /* transform for factor of 4 */ 489 ispan = kspan; 490 kspan /= 4; 491 492 do { 493 c1 = 1.0; 494 s1 = 0.0; 495 do { 496 do { 497 k1 = kk + kspan; 498 k2 = k1 + kspan; 499 k3 = k2 + kspan; 500 akp = Re [kk] + Re [k2]; 501 akm = Re [kk] - Re [k2]; 502 ajp = Re [k1] + Re [k3]; 503 ajm = Re [k1] - Re [k3]; 504 bkp = Im [kk] + Im [k2]; 505 bkm = Im [kk] - Im [k2]; 506 bjp = Im [k1] + Im [k3]; 507 bjm = Im [k1] - Im [k3]; 508 Re [kk] = akp + ajp; 509 Im [kk] = bkp + bjp; 510 ajp = akp - ajp; 511 bjp = bkp - bjp; 512 if (iSign < 0) { 513 akp = akm + bjm; 514 bkp = bkm - ajm; 515 akm -= bjm; 516 bkm += ajm; 517 } else { 518 akp = akm - bjm; 519 bkp = bkm + ajm; 520 akm += bjm; 521 bkm -= ajm; 522 } 523 /* avoid useless multiplies */ 524 if (s1 == 0.0) { 525 Re [k1] = akp; 526 Re [k2] = ajp; 527 Re [k3] = akm; 528 Im [k1] = bkp; 529 Im [k2] = bjp; 530 Im [k3] = bkm; 531 } else { 532 Re [k1] = akp * c1 - bkp * s1; 533 Re [k2] = ajp * c2 - bjp * s2; 534 Re [k3] = akm * c3 - bkm * s3; 535 Im [k1] = akp * s1 + bkp * c1; 536 Im [k2] = ajp * s2 + bjp * c2; 537 Im [k3] = akm * s3 + bkm * c3; 538 } 539 kk = k3 + kspan; 540 } while (kk < nt); 541 542 c2 = c1 - (cd * c1 + sd * s1); 543 s1 = sd * c1 - cd * s1 + s1; 544 c1 = 2.0 - (c2 * c2 + s1 * s1); 545 s1 *= c1; 546 c1 *= c2; 547 /* values of c2, c3, s2, s3 that will get used next time */ 548 c2 = c1 * c1 - s1 * s1; 549 s2 = 2.0 * c1 * s1; 550 c3 = c2 * c1 - s2 * s1; 551 s3 = c2 * s1 + s2 * c1; 552 kk = kk - nt + jc; 553 } while (kk < kspan); 554 kk = kk - kspan + inc; 555 } while (kk < jc); 556 if (kspan == jc) 557 goto Permute_Results_Label; /* exit infinite loop */ 558 break; 559 560 default: 561 /* transform for odd factors */ 562 #ifdef FFT_RADIX4 563 return -1; 564 break; 565 #else /* FFT_RADIX4 */ 566 k = fftstate->factor [ii - 1]; 567 ispan = kspan; 568 kspan /= k; 569 570 switch (k) { 571 case 3: /* transform for factor of 3 (optional code) */ 572 do { 573 do { 574 k1 = kk + kspan; 575 k2 = k1 + kspan; 576 ak = Re [kk]; 577 bk = Im [kk]; 578 aj = Re [k1] + Re [k2]; 579 bj = Im [k1] + Im [k2]; 580 Re [kk] = ak + aj; 581 Im [kk] = bk + bj; 582 ak -= 0.5 * aj; 583 bk -= 0.5 * bj; 584 aj = (Re [k1] - Re [k2]) * s60; 585 bj = (Im [k1] - Im [k2]) * s60; 586 Re [k1] = ak - bj; 587 Re [k2] = ak + bj; 588 Im [k1] = bk + aj; 589 Im [k2] = bk - aj; 590 kk = k2 + kspan; 591 } while (kk < (nn - 1)); 592 kk -= nn; 593 } while (kk < kspan); 594 break; 595 596 case 5: /* transform for factor of 5 (optional code) */ 597 c2 = c72 * c72 - s72 * s72; 598 s2 = 2.0 * c72 * s72; 599 do { 600 do { 601 k1 = kk + kspan; 602 k2 = k1 + kspan; 603 k3 = k2 + kspan; 604 k4 = k3 + kspan; 605 akp = Re [k1] + Re [k4]; 606 akm = Re [k1] - Re [k4]; 607 bkp = Im [k1] + Im [k4]; 608 bkm = Im [k1] - Im [k4]; 609 ajp = Re [k2] + Re [k3]; 610 ajm = Re [k2] - Re [k3]; 611 bjp = Im [k2] + Im [k3]; 612 bjm = Im [k2] - Im [k3]; 613 aa = Re [kk]; 614 bb = Im [kk]; 615 Re [kk] = aa + akp + ajp; 616 Im [kk] = bb + bkp + bjp; 617 ak = akp * c72 + ajp * c2 + aa; 618 bk = bkp * c72 + bjp * c2 + bb; 619 aj = akm * s72 + ajm * s2; 620 bj = bkm * s72 + bjm * s2; 621 Re [k1] = ak - bj; 622 Re [k4] = ak + bj; 623 Im [k1] = bk + aj; 624 Im [k4] = bk - aj; 625 ak = akp * c2 + ajp * c72 + aa; 626 bk = bkp * c2 + bjp * c72 + bb; 627 aj = akm * s2 - ajm * s72; 628 bj = bkm * s2 - bjm * s72; 629 Re [k2] = ak - bj; 630 Re [k3] = ak + bj; 631 Im [k2] = bk + aj; 632 Im [k3] = bk - aj; 633 kk = k4 + kspan; 634 } while (kk < (nn-1)); 635 kk -= nn; 636 } while (kk < kspan); 637 break; 638 639 default: 640 if (k != jf) { 641 jf = k; 642 s1 = pi2 / (double) k; 643 c1 = cos(s1); 644 s1 = sin(s1); 645 if (jf > max_factors){ 646 return -1; 647 } 648 Cos [jf - 1] = 1.0; 649 Sin [jf - 1] = 0.0; 650 j = 1; 651 do { 652 Cos [j - 1] = Cos [k - 1] * c1 + Sin [k - 1] * s1; 653 Sin [j - 1] = Cos [k - 1] * s1 - Sin [k - 1] * c1; 654 k--; 655 Cos [k - 1] = Cos [j - 1]; 656 Sin [k - 1] = -Sin [j - 1]; 657 j++; 658 } while (j < k); 659 } 660 do { 661 do { 662 k1 = kk; 663 k2 = kk + ispan; 664 ak = aa = Re [kk]; 665 bk = bb = Im [kk]; 666 j = 1; 667 k1 += kspan; 668 do { 669 k2 -= kspan; 670 j++; 671 Rtmp [j - 1] = Re [k1] + Re [k2]; 672 ak += Rtmp [j - 1]; 673 Itmp [j - 1] = Im [k1] + Im [k2]; 674 bk += Itmp [j - 1]; 675 j++; 676 Rtmp [j - 1] = Re [k1] - Re [k2]; 677 Itmp [j - 1] = Im [k1] - Im [k2]; 678 k1 += kspan; 679 } while (k1 < k2); 680 Re [kk] = ak; 681 Im [kk] = bk; 682 k1 = kk; 683 k2 = kk + ispan; 684 j = 1; 685 do { 686 k1 += kspan; 687 k2 -= kspan; 688 jj = j; 689 ak = aa; 690 bk = bb; 691 aj = 0.0; 692 bj = 0.0; 693 k = 1; 694 do { 695 k++; 696 ak += Rtmp [k - 1] * Cos [jj - 1]; 697 bk += Itmp [k - 1] * Cos [jj - 1]; 698 k++; 699 aj += Rtmp [k - 1] * Sin [jj - 1]; 700 bj += Itmp [k - 1] * Sin [jj - 1]; 701 jj += j; 702 if (jj > jf) { 703 jj -= jf; 704 } 705 } while (k < jf); 706 k = jf - j; 707 Re [k1] = ak - bj; 708 Im [k1] = bk + aj; 709 Re [k2] = ak + bj; 710 Im [k2] = bk - aj; 711 j++; 712 } while (j < k); 713 kk += ispan; 714 } while (kk < nn); 715 kk -= nn; 716 } while (kk < kspan); 717 break; 718 } 719 720 /* multiply by rotation factor (except for factors of 2 and 4) */ 721 if (ii == mfactor) 722 goto Permute_Results_Label; /* exit infinite loop */ 723 kk = jc; 724 do { 725 c2 = 1.0 - cd; 726 s1 = sd; 727 do { 728 c1 = c2; 729 s2 = s1; 730 kk += kspan; 731 do { 732 do { 733 ak = Re [kk]; 734 Re [kk] = c2 * ak - s2 * Im [kk]; 735 Im [kk] = s2 * ak + c2 * Im [kk]; 736 kk += ispan; 737 } while (kk < nt); 738 ak = s1 * s2; 739 s2 = s1 * c2 + c1 * s2; 740 c2 = c1 * c2 - ak; 741 kk = kk - nt + kspan; 742 } while (kk < ispan); 743 c2 = c1 - (cd * c1 + sd * s1); 744 s1 += sd * c1 - cd * s1; 745 c1 = 2.0 - (c2 * c2 + s1 * s1); 746 s1 *= c1; 747 c2 *= c1; 748 kk = kk - ispan + jc; 749 } while (kk < kspan); 750 kk = kk - kspan + jc + inc; 751 } while (kk < (jc + jc)); 752 break; 753 #endif /* FFT_RADIX4 */ 754 } 755 } 756 757 /* permute the results to normal order---done in two stages */ 758 /* permutation for square factors of n */ 759 Permute_Results_Label: 760 fftstate->Perm [0] = ns; 761 if (kt) { 762 k = kt + kt + 1; 763 if (mfactor < k) 764 k--; 765 j = 1; 766 fftstate->Perm [k] = jc; 767 do { 768 fftstate->Perm [j] = fftstate->Perm [j - 1] / fftstate->factor [j - 1]; 769 fftstate->Perm [k - 1] = fftstate->Perm [k] * fftstate->factor [j - 1]; 770 j++; 771 k--; 772 } while (j < k); 773 k3 = fftstate->Perm [k]; 774 kspan = fftstate->Perm [1]; 775 kk = jc; 776 k2 = kspan; 777 j = 1; 778 if (nPass != nTotal) { 779 /* permutation for multivariate transform */ 780 Permute_Multi_Label: 781 do { 782 do { 783 k = kk + jc; 784 do { 785 /* swap Re [kk] <> Re [k2], Im [kk] <> Im [k2] */ 786 ak = Re [kk]; Re [kk] = Re [k2]; Re [k2] = ak; 787 bk = Im [kk]; Im [kk] = Im [k2]; Im [k2] = bk; 788 kk += inc; 789 k2 += inc; 790 } while (kk < (k-1)); 791 kk += ns - jc; 792 k2 += ns - jc; 793 } while (kk < (nt-1)); 794 k2 = k2 - nt + kspan; 795 kk = kk - nt + jc; 796 } while (k2 < (ns-1)); 797 do { 798 do { 799 k2 -= fftstate->Perm [j - 1]; 800 j++; 801 k2 = fftstate->Perm [j] + k2; 802 } while (k2 > fftstate->Perm [j - 1]); 803 j = 1; 804 do { 805 if (kk < (k2-1)) 806 goto Permute_Multi_Label; 807 kk += jc; 808 k2 += kspan; 809 } while (k2 < (ns-1)); 810 } while (kk < (ns-1)); 811 } else { 812 /* permutation for single-variate transform (optional code) */ 813 Permute_Single_Label: 814 do { 815 /* swap Re [kk] <> Re [k2], Im [kk] <> Im [k2] */ 816 ak = Re [kk]; Re [kk] = Re [k2]; Re [k2] = ak; 817 bk = Im [kk]; Im [kk] = Im [k2]; Im [k2] = bk; 818 kk += inc; 819 k2 += kspan; 820 } while (k2 < (ns-1)); 821 do { 822 do { 823 k2 -= fftstate->Perm [j - 1]; 824 j++; 825 k2 = fftstate->Perm [j] + k2; 826 } while (k2 >= fftstate->Perm [j - 1]); 827 j = 1; 828 do { 829 if (kk < k2) 830 goto Permute_Single_Label; 831 kk += inc; 832 k2 += kspan; 833 } while (k2 < (ns-1)); 834 } while (kk < (ns-1)); 835 } 836 jc = k3; 837 } 838 839 if ((kt << 1) + 1 >= mfactor) 840 return 0; 841 ispan = fftstate->Perm [kt]; 842 /* permutation for square-free factors of n */ 843 j = mfactor - kt; 844 fftstate->factor [j] = 1; 845 do { 846 fftstate->factor [j - 1] *= fftstate->factor [j]; 847 j--; 848 } while (j != kt); 849 kt++; 850 nn = fftstate->factor [kt - 1] - 1; 851 if (nn > (int) max_perm) { 852 return -1; 853 } 854 j = jj = 0; 855 for (;;) { 856 k = kt + 1; 857 k2 = fftstate->factor [kt - 1]; 858 kk = fftstate->factor [k - 1]; 859 j++; 860 if (j > nn) 861 break; /* exit infinite loop */ 862 jj += kk; 863 while (jj >= k2) { 864 jj -= k2; 865 k2 = kk; 866 k++; 867 kk = fftstate->factor [k - 1]; 868 jj += kk; 869 } 870 fftstate->Perm [j - 1] = jj; 871 } 872 /* determine the permutation cycles of length greater than 1 */ 873 j = 0; 874 for (;;) { 875 do { 876 j++; 877 kk = fftstate->Perm [j - 1]; 878 } while (kk < 0); 879 if (kk != j) { 880 do { 881 k = kk; 882 kk = fftstate->Perm [k - 1]; 883 fftstate->Perm [k - 1] = -kk; 884 } while (kk != j); 885 k3 = kk; 886 } else { 887 fftstate->Perm [j - 1] = -j; 888 if (j == nn) 889 break; /* exit infinite loop */ 890 } 891 } 892 max_factors *= inc; 893 /* reorder a and b, following the permutation cycles */ 894 for (;;) { 895 j = k3 + 1; 896 nt -= ispan; 897 ii = nt - inc + 1; 898 if (nt < 0) 899 break; /* exit infinite loop */ 900 do { 901 do { 902 j--; 903 } while (fftstate->Perm [j - 1] < 0); 904 jj = jc; 905 do { 906 kspan = jj; 907 if (jj > max_factors) { 908 kspan = max_factors; 909 } 910 jj -= kspan; 911 k = fftstate->Perm [j - 1]; 912 kk = jc * k + ii + jj; 913 k1 = kk + kspan - 1; 914 k2 = 0; 915 do { 916 k2++; 917 Rtmp [k2 - 1] = Re [k1]; 918 Itmp [k2 - 1] = Im [k1]; 919 k1 -= inc; 920 } while (k1 != (kk-1)); 921 do { 922 k1 = kk + kspan - 1; 923 k2 = k1 - jc * (k + fftstate->Perm [k - 1]); 924 k = -fftstate->Perm [k - 1]; 925 do { 926 Re [k1] = Re [k2]; 927 Im [k1] = Im [k2]; 928 k1 -= inc; 929 k2 -= inc; 930 } while (k1 != (kk-1)); 931 kk = k2 + 1; 932 } while (k != j); 933 k1 = kk + kspan - 1; 934 k2 = 0; 935 do { 936 k2++; 937 Re [k1] = Rtmp [k2 - 1]; 938 Im [k1] = Itmp [k2 - 1]; 939 k1 -= inc; 940 } while (k1 != (kk-1)); 941 } while (jj); 942 } while (j != 1); 943 } 944 return 0; /* exit point here */ 945 } 946 /* ---------------------- end-of-file (c source) ---------------------- */ 947 948
