00001
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 #ifdef NATIVE_FFT
00037
00038 #include <cstring>
00039 #include <cmath>
00040 #include "native_fft.h"
00041 #include "log.h"
00042
00043 #include <stdio.h>
00044 #include <stdlib.h>
00045 #include <math.h>
00046
00047
00048 #include "util.h"
00049
00050 using namespace EMAN;
00051
00052 int Nativefft::fftmcf_(float *a, float *b, integer *ntot,
00053 integer *n, integer *nspan, integer *isn)
00054 {
00055
00056 integer i__1;
00057 float r__1, r__2;
00058 int status = -1;
00059
00060
00061
00062
00063
00064 static integer nfac[11];
00065 static float radf;
00066 static integer maxf, maxp, i__, j, k, m, kspan;
00067 static float c1, c2, c3;
00068 static integer kspnn, k1, k2, k3, k4;
00069 static float s1, s2, s3, aa, bb, cd, aj, c72;
00070 static integer jc;
00071 static float ck[23], ak;
00072 static integer jf;
00073 static float bk, bj;
00074 static integer jj;
00075 static float at[23], bt[23], sd;
00076 static integer kk;
00077 static float s72;
00078 static integer nn, np[209];
00079 static float sk[23];
00080 static integer ks, kt, nt;
00081 static float s120, rad, ajm, akm;
00082 static integer inc;
00083 static float ajp, akp, bkp, bkm, bjp, bjm;
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094 --b;
00095 --a;
00096
00097
00098 maxf = 23;
00099 maxp = 209;
00100 if (*n < 2) {
00101 return 0;
00102 }
00103 inc = *isn;
00104 rad = atan((float)1.) * (float)8.;
00105 s72 = rad / (float)5.;
00106 c72 = cos(s72);
00107 s72 = sin(s72);
00108 s120 = sqrt((float).75);
00109 if (*isn >= 0) {
00110 goto L10;
00111 }
00112 s72 = -s72;
00113 s120 = -s120;
00114 rad = -rad;
00115 inc = -inc;
00116 L10:
00117 nt = inc * *ntot;
00118 ks = inc * *nspan;
00119 kspan = ks;
00120 nn = nt - inc;
00121 jc = ks / *n;
00122 radf = rad * (float) jc * (float).5;
00123 i__ = 0;
00124 jf = 0;
00125
00126
00127
00128 m = 0;
00129 k = *n;
00130 goto L20;
00131 L15:
00132 ++m;
00133 nfac[m - 1] = 4;
00134 k /= 16;
00135 L20:
00136 if (k - (k / 16 << 4) == 0) {
00137 goto L15;
00138 }
00139 j = 3;
00140 jj = 9;
00141 goto L30;
00142 L25:
00143 ++m;
00144 nfac[m - 1] = j;
00145 k /= jj;
00146 L30:
00147 if (k % jj == 0) {
00148 goto L25;
00149 }
00150 j += 2;
00151
00152 i__1 = j;
00153 jj = i__1 * i__1;
00154 if (jj <= k) {
00155 goto L30;
00156 }
00157 if (k > 4) {
00158 goto L40;
00159 }
00160 kt = m;
00161 nfac[m] = k;
00162 if (k != 1) {
00163 ++m;
00164 }
00165 goto L80;
00166 L40:
00167 if (k - (k / 4 << 2) != 0) {
00168 goto L50;
00169 }
00170 ++m;
00171 nfac[m - 1] = 2;
00172 k /= 4;
00173 L50:
00174 kt = m;
00175 j = 2;
00176 L60:
00177 if (k % j != 0) {
00178 goto L70;
00179 }
00180 ++m;
00181 nfac[m - 1] = j;
00182 k /= j;
00183 L70:
00184 j = ((j + 1) / 2 << 1) + 1;
00185 if (j <= k) {
00186 goto L60;
00187 }
00188 L80:
00189 if (kt == 0) {
00190 goto L100;
00191 }
00192 j = kt;
00193 L90:
00194 ++m;
00195 nfac[m - 1] = nfac[j - 1];
00196 --j;
00197 if (j != 0) {
00198 goto L90;
00199 }
00200
00201
00202
00203 L100:
00204 sd = radf / (float) kspan;
00205
00206 r__1 = sin(sd);
00207 cd = r__1 * r__1 * (float)2.;
00208 sd = sin(sd + sd);
00209 kk = 1;
00210 ++i__;
00211 if (nfac[i__ - 1] != 2) {
00212 goto L400;
00213 }
00214
00215
00216
00217 kspan /= 2;
00218 k1 = kspan + 2;
00219 L210:
00220 k2 = kk + kspan;
00221 ak = a[k2];
00222 bk = b[k2];
00223 a[k2] = a[kk] - ak;
00224 b[k2] = b[kk] - bk;
00225 a[kk] += ak;
00226 b[kk] += bk;
00227 kk = k2 + kspan;
00228 if (kk <= nn) {
00229 goto L210;
00230 }
00231 kk -= nn;
00232 if (kk <= jc) {
00233 goto L210;
00234 }
00235 if (kk > kspan) {
00236 goto L800;
00237 }
00238 L220:
00239 c1 = (float)1. - cd;
00240 s1 = sd;
00241 L230:
00242 k2 = kk + kspan;
00243 ak = a[kk] - a[k2];
00244 bk = b[kk] - b[k2];
00245 a[kk] += a[k2];
00246 b[kk] += b[k2];
00247 a[k2] = c1 * ak - s1 * bk;
00248 b[k2] = s1 * ak + c1 * bk;
00249 kk = k2 + kspan;
00250 if (kk < nt) {
00251 goto L230;
00252 }
00253 k2 = kk - nt;
00254 c1 = -c1;
00255 kk = k1 - k2;
00256 if (kk > k2) {
00257 goto L230;
00258 }
00259 ak = c1 - (cd * c1 + sd * s1);
00260 s1 = sd * c1 - cd * s1 + s1;
00261
00262
00263
00264
00265
00266
00267 r__1 = ak;
00268
00269 r__2 = s1;
00270 c1 = (float).5 / (r__1 * r__1 + r__2 * r__2) + (float).5;
00271 s1 = c1 * s1;
00272 c1 *= ak;
00273 kk += jc;
00274 if (kk < k2) {
00275 goto L230;
00276 }
00277 k1 = k1 + inc + inc;
00278 kk = (k1 - kspan) / 2 + jc;
00279 if (kk <= jc + jc) {
00280 goto L220;
00281 }
00282 goto L100;
00283
00284
00285
00286 L320:
00287 k1 = kk + kspan;
00288 k2 = k1 + kspan;
00289 ak = a[kk];
00290 bk = b[kk];
00291 aj = a[k1] + a[k2];
00292 bj = b[k1] + b[k2];
00293 a[kk] = ak + aj;
00294 b[kk] = bk + bj;
00295 ak = aj * (float)-.5 + ak;
00296 bk = bj * (float)-.5 + bk;
00297 aj = (a[k1] - a[k2]) * s120;
00298 bj = (b[k1] - b[k2]) * s120;
00299 a[k1] = ak - bj;
00300 b[k1] = bk + aj;
00301 a[k2] = ak + bj;
00302 b[k2] = bk - aj;
00303 kk = k2 + kspan;
00304 if (kk < nn) {
00305 goto L320;
00306 }
00307 kk -= nn;
00308 if (kk <= kspan) {
00309 goto L320;
00310 }
00311 goto L700;
00312
00313
00314
00315 L400:
00316 if (nfac[i__ - 1] != 4) {
00317 goto L600;
00318 }
00319 kspnn = kspan;
00320 kspan /= 4;
00321 L410:
00322 c1 = (float)1.;
00323 s1 = (float)0.;
00324 L420:
00325 k1 = kk + kspan;
00326 k2 = k1 + kspan;
00327 k3 = k2 + kspan;
00328 akp = a[kk] + a[k2];
00329 akm = a[kk] - a[k2];
00330 ajp = a[k1] + a[k3];
00331 ajm = a[k1] - a[k3];
00332 a[kk] = akp + ajp;
00333 ajp = akp - ajp;
00334 bkp = b[kk] + b[k2];
00335 bkm = b[kk] - b[k2];
00336 bjp = b[k1] + b[k3];
00337 bjm = b[k1] - b[k3];
00338 b[kk] = bkp + bjp;
00339 bjp = bkp - bjp;
00340 if (*isn < 0) {
00341 goto L450;
00342 }
00343 akp = akm - bjm;
00344 akm += bjm;
00345 bkp = bkm + ajm;
00346 bkm -= ajm;
00347 if (s1 == (float)0.) {
00348 goto L460;
00349 }
00350 L430:
00351 a[k1] = akp * c1 - bkp * s1;
00352 b[k1] = akp * s1 + bkp * c1;
00353 a[k2] = ajp * c2 - bjp * s2;
00354 b[k2] = ajp * s2 + bjp * c2;
00355 a[k3] = akm * c3 - bkm * s3;
00356 b[k3] = akm * s3 + bkm * c3;
00357 kk = k3 + kspan;
00358 if (kk <= nt) {
00359 goto L420;
00360 }
00361 L440:
00362 c2 = c1 - (cd * c1 + sd * s1);
00363 s1 = sd * c1 - cd * s1 + s1;
00364
00365
00366
00367
00368
00369
00370 r__1 = c2;
00371
00372 r__2 = s1;
00373 c1 = (float).5 / (r__1 * r__1 + r__2 * r__2) + (float).5;
00374 s1 = c1 * s1;
00375 c1 *= c2;
00376
00377 r__1 = c1;
00378
00379 r__2 = s1;
00380 c2 = r__1 * r__1 - r__2 * r__2;
00381 s2 = c1 * (float)2. * s1;
00382 c3 = c2 * c1 - s2 * s1;
00383 s3 = c2 * s1 + s2 * c1;
00384 kk = kk - nt + jc;
00385 if (kk <= kspan) {
00386 goto L420;
00387 }
00388 kk = kk - kspan + inc;
00389 if (kk <= jc) {
00390 goto L410;
00391 }
00392 if (kspan == jc) {
00393 goto L800;
00394 }
00395 goto L100;
00396 L450:
00397 akp = akm + bjm;
00398 akm -= bjm;
00399 bkp = bkm - ajm;
00400 bkm += ajm;
00401 if (s1 != (float)0.) {
00402 goto L430;
00403 }
00404 L460:
00405 a[k1] = akp;
00406 b[k1] = bkp;
00407 a[k2] = ajp;
00408 b[k2] = bjp;
00409 a[k3] = akm;
00410 b[k3] = bkm;
00411 kk = k3 + kspan;
00412 if (kk <= nt) {
00413 goto L420;
00414 }
00415 goto L440;
00416
00417
00418
00419 L510:
00420
00421 r__1 = c72;
00422
00423 r__2 = s72;
00424 c2 = r__1 * r__1 - r__2 * r__2;
00425 s2 = c72 * (float)2. * s72;
00426 L520:
00427 k1 = kk + kspan;
00428 k2 = k1 + kspan;
00429 k3 = k2 + kspan;
00430 k4 = k3 + kspan;
00431 akp = a[k1] + a[k4];
00432 akm = a[k1] - a[k4];
00433 bkp = b[k1] + b[k4];
00434 bkm = b[k1] - b[k4];
00435 ajp = a[k2] + a[k3];
00436 ajm = a[k2] - a[k3];
00437 bjp = b[k2] + b[k3];
00438 bjm = b[k2] - b[k3];
00439 aa = a[kk];
00440 bb = b[kk];
00441 a[kk] = aa + akp + ajp;
00442 b[kk] = bb + bkp + bjp;
00443 ak = akp * c72 + ajp * c2 + aa;
00444 bk = bkp * c72 + bjp * c2 + bb;
00445 aj = akm * s72 + ajm * s2;
00446 bj = bkm * s72 + bjm * s2;
00447 a[k1] = ak - bj;
00448 a[k4] = ak + bj;
00449 b[k1] = bk + aj;
00450 b[k4] = bk - aj;
00451 ak = akp * c2 + ajp * c72 + aa;
00452 bk = bkp * c2 + bjp * c72 + bb;
00453 aj = akm * s2 - ajm * s72;
00454 bj = bkm * s2 - bjm * s72;
00455 a[k2] = ak - bj;
00456 a[k3] = ak + bj;
00457 b[k2] = bk + aj;
00458 b[k3] = bk - aj;
00459 kk = k4 + kspan;
00460 if (kk < nn) {
00461 goto L520;
00462 }
00463 kk -= nn;
00464 if (kk <= kspan) {
00465 goto L520;
00466 }
00467 goto L700;
00468
00469
00470
00471 L600:
00472 k = nfac[i__ - 1];
00473 kspnn = kspan;
00474 kspan /= k;
00475 if (k == 3) {
00476 goto L320;
00477 }
00478 if (k == 5) {
00479 goto L510;
00480 }
00481 if (k == jf) {
00482 goto L640;
00483 }
00484 jf = k;
00485 s1 = rad / (float) k;
00486 c1 = cos(s1);
00487 s1 = sin(s1);
00488 if (jf > maxf) {
00489 goto L998;
00490 }
00491 ck[jf - 1] = (float)1.;
00492 sk[jf - 1] = (float)0.;
00493 j = 1;
00494 L630:
00495 ck[j - 1] = ck[k - 1] * c1 + sk[k - 1] * s1;
00496 sk[j - 1] = ck[k - 1] * s1 - sk[k - 1] * c1;
00497 --k;
00498 ck[k - 1] = ck[j - 1];
00499 sk[k - 1] = -sk[j - 1];
00500 ++j;
00501 if (j < k) {
00502 goto L630;
00503 }
00504 L640:
00505 k1 = kk;
00506 k2 = kk + kspnn;
00507 aa = a[kk];
00508 bb = b[kk];
00509 ak = aa;
00510 bk = bb;
00511 j = 1;
00512 k1 += kspan;
00513 L650:
00514 k2 -= kspan;
00515 ++j;
00516 at[j - 1] = a[k1] + a[k2];
00517 ak = at[j - 1] + ak;
00518 bt[j - 1] = b[k1] + b[k2];
00519 bk = bt[j - 1] + bk;
00520 ++j;
00521 at[j - 1] = a[k1] - a[k2];
00522 bt[j - 1] = b[k1] - b[k2];
00523 k1 += kspan;
00524 if (k1 < k2) {
00525 goto L650;
00526 }
00527 a[kk] = ak;
00528 b[kk] = bk;
00529 k1 = kk;
00530 k2 = kk + kspnn;
00531 j = 1;
00532 L660:
00533 k1 += kspan;
00534 k2 -= kspan;
00535 jj = j;
00536 ak = aa;
00537 bk = bb;
00538 aj = (float)0.;
00539 bj = (float)0.;
00540 k = 1;
00541 L670:
00542 ++k;
00543 ak = at[k - 1] * ck[jj - 1] + ak;
00544 bk = bt[k - 1] * ck[jj - 1] + bk;
00545 ++k;
00546 aj = at[k - 1] * sk[jj - 1] + aj;
00547 bj = bt[k - 1] * sk[jj - 1] + bj;
00548 jj += j;
00549 if (jj > jf) {
00550 jj -= jf;
00551 }
00552 if (k < jf) {
00553 goto L670;
00554 }
00555 k = jf - j;
00556 a[k1] = ak - bj;
00557 b[k1] = bk + aj;
00558 a[k2] = ak + bj;
00559 b[k2] = bk - aj;
00560 ++j;
00561 if (j < k) {
00562 goto L660;
00563 }
00564 kk += kspnn;
00565 if (kk <= nn) {
00566 goto L640;
00567 }
00568 kk -= nn;
00569 if (kk <= kspan) {
00570 goto L640;
00571 }
00572
00573
00574
00575 L700:
00576 if (i__ == m) {
00577 goto L800;
00578 }
00579 kk = jc + 1;
00580 L710:
00581 c2 = (float)1. - cd;
00582 s1 = sd;
00583 L720:
00584 c1 = c2;
00585 s2 = s1;
00586 kk += kspan;
00587 L730:
00588 ak = a[kk];
00589 a[kk] = c2 * ak - s2 * b[kk];
00590 b[kk] = s2 * ak + c2 * b[kk];
00591 kk += kspnn;
00592 if (kk <= nt) {
00593 goto L730;
00594 }
00595 ak = s1 * s2;
00596 s2 = s1 * c2 + c1 * s2;
00597 c2 = c1 * c2 - ak;
00598 kk = kk - nt + kspan;
00599 if (kk <= kspnn) {
00600 goto L730;
00601 }
00602 c2 = c1 - (cd * c1 + sd * s1);
00603 s1 += sd * c1 - cd * s1;
00604
00605
00606
00607
00608
00609
00610 r__1 = c2;
00611
00612 r__2 = s1;
00613 c1 = (float).5 / (r__1 * r__1 + r__2 * r__2) + (float).5;
00614 s1 = c1 * s1;
00615 c2 = c1 * c2;
00616 kk = kk - kspnn + jc;
00617 if (kk <= kspan) {
00618 goto L720;
00619 }
00620 kk = kk - kspan + jc + inc;
00621 if (kk <= jc + jc) {
00622 goto L710;
00623 }
00624 goto L100;
00625
00626
00627
00628
00629 L800:
00630 np[0] = ks;
00631 if (kt == 0) {
00632 goto L890;
00633 }
00634 k = kt + kt + 1;
00635 if (m < k) {
00636 --k;
00637 }
00638 j = 1;
00639 np[k] = jc;
00640 L810:
00641 np[j] = np[j - 1] / nfac[j - 1];
00642 np[k - 1] = np[k] * nfac[j - 1];
00643 ++j;
00644 --k;
00645 if (j < k) {
00646 goto L810;
00647 }
00648 k3 = np[k];
00649 kspan = np[1];
00650 kk = jc + 1;
00651 j = 1;
00652 k2 = kspan + 1;
00653 if (*n != *ntot) {
00654 goto L850;
00655 }
00656
00657
00658
00659 L820:
00660 ak = a[kk];
00661 a[kk] = a[k2];
00662 a[k2] = ak;
00663 bk = b[kk];
00664 b[kk] = b[k2];
00665 b[k2] = bk;
00666 kk += inc;
00667 k2 = kspan + k2;
00668 if (k2 < ks) {
00669 goto L820;
00670 }
00671 L830:
00672 k2 -= np[j - 1];
00673 ++j;
00674 k2 = np[j] + k2;
00675 if (k2 > np[j - 1]) {
00676 goto L830;
00677 }
00678 j = 1;
00679 L840:
00680 if (kk < k2) {
00681 goto L820;
00682 }
00683 kk += inc;
00684 k2 = kspan + k2;
00685 if (k2 < ks) {
00686 goto L840;
00687 }
00688 if (kk < ks) {
00689 goto L830;
00690 }
00691 jc = k3;
00692 goto L890;
00693
00694
00695
00696 L850:
00697 k = kk + jc;
00698 L860:
00699 ak = a[kk];
00700 a[kk] = a[k2];
00701 a[k2] = ak;
00702 bk = b[kk];
00703 b[kk] = b[k2];
00704 b[k2] = bk;
00705 kk += inc;
00706 k2 += inc;
00707 if (kk < k) {
00708 goto L860;
00709 }
00710 kk = kk + ks - jc;
00711 k2 = k2 + ks - jc;
00712 if (kk < nt) {
00713 goto L850;
00714 }
00715 k2 = k2 - nt + kspan;
00716 kk = kk - nt + jc;
00717 if (k2 < ks) {
00718 goto L850;
00719 }
00720 L870:
00721 k2 -= np[j - 1];
00722 ++j;
00723 k2 = np[j] + k2;
00724 if (k2 > np[j - 1]) {
00725 goto L870;
00726 }
00727 j = 1;
00728 L880:
00729 if (kk < k2) {
00730 goto L850;
00731 }
00732 kk += jc;
00733 k2 = kspan + k2;
00734 if (k2 < ks) {
00735 goto L880;
00736 }
00737 if (kk < ks) {
00738 goto L870;
00739 }
00740 jc = k3;
00741 L890:
00742 if ((kt << 1) + 1 >= m) {
00743 return 0;
00744 }
00745 kspnn = np[kt];
00746
00747
00748
00749 j = m - kt;
00750 nfac[j] = 1;
00751 L900:
00752 nfac[j - 1] *= nfac[j];
00753 --j;
00754 if (j != kt) {
00755 goto L900;
00756 }
00757 ++kt;
00758 nn = nfac[kt - 1] - 1;
00759 if (nn > maxp) {
00760 goto L998;
00761 }
00762 jj = 0;
00763 j = 0;
00764 goto L906;
00765 L902:
00766 jj -= k2;
00767 k2 = kk;
00768 ++k;
00769 kk = nfac[k - 1];
00770 L904:
00771 jj = kk + jj;
00772 if (jj >= k2) {
00773 goto L902;
00774 }
00775 np[j - 1] = jj;
00776 L906:
00777 k2 = nfac[kt - 1];
00778 k = kt + 1;
00779 kk = nfac[k - 1];
00780 ++j;
00781 if (j <= nn) {
00782 goto L904;
00783 }
00784
00785
00786
00787 j = 0;
00788 goto L914;
00789 L910:
00790 k = kk;
00791 kk = np[k - 1];
00792 np[k - 1] = -kk;
00793 if (kk != j) {
00794 goto L910;
00795 }
00796 k3 = kk;
00797 L914:
00798 ++j;
00799 kk = np[j - 1];
00800 if (kk < 0) {
00801 goto L914;
00802 }
00803 if (kk != j) {
00804 goto L910;
00805 }
00806 np[j - 1] = -j;
00807 if (j != nn) {
00808 goto L914;
00809 }
00810 maxf = inc * maxf;
00811
00812
00813
00814 goto L950;
00815 L924:
00816 --j;
00817 if (np[j - 1] < 0) {
00818 goto L924;
00819 }
00820 jj = jc;
00821 L926:
00822 kspan = jj;
00823 if (jj > maxf) {
00824 kspan = maxf;
00825 }
00826 jj -= kspan;
00827 k = np[j - 1];
00828 kk = jc * k + i__ + jj;
00829 k1 = kk + kspan;
00830 k2 = 0;
00831 L928:
00832 ++k2;
00833 at[k2 - 1] = a[k1];
00834 bt[k2 - 1] = b[k1];
00835 k1 -= inc;
00836 if (k1 != kk) {
00837 goto L928;
00838 }
00839 L932:
00840 k1 = kk + kspan;
00841 k2 = k1 - jc * (k + np[k - 1]);
00842 k = -np[k - 1];
00843 L936:
00844 a[k1] = a[k2];
00845 b[k1] = b[k2];
00846 k1 -= inc;
00847 k2 -= inc;
00848 if (k1 != kk) {
00849 goto L936;
00850 }
00851 kk = k2;
00852 if (k != j) {
00853 goto L932;
00854 }
00855 k1 = kk + kspan;
00856 k2 = 0;
00857 L940:
00858 ++k2;
00859 a[k1] = at[k2 - 1];
00860 b[k1] = bt[k2 - 1];
00861 k1 -= inc;
00862 if (k1 != kk) {
00863 goto L940;
00864 }
00865 if (jj != 0) {
00866 goto L926;
00867 }
00868 if (j != 1) {
00869 goto L924;
00870 }
00871 L950:
00872 j = k3 + 1;
00873 nt -= kspnn;
00874 i__ = nt - inc + 1;
00875 if (nt >= 0) {
00876 goto L924;
00877 }
00878 return 0;
00879
00880
00881
00882 L998:
00883 *isn = 0;
00884 return status;
00885 }
00886
00887 #define work(i) work[(i)-1]
00888 #define x(i) x[(i)-1]
00889 #define y(i) y[(i)-1]
00890
00891
00892 int Nativefft::fmrs_1rf(float *x, float *work, int nsam)
00893 {
00894
00895
00896
00897 int i, n, status;
00898 int inv=1;
00899
00900 n = nsam;
00901
00902 for (i=1; i<=n; i++) work(i) = 0.0;
00903 status = fftmcf_(x,work,&n,&n,&n,&inv);
00904
00905 if (status == -1) {
00906 fprintf(stderr, "Native FFT cannot be performed on images with one of ");
00907 fprintf(stderr, "the dimensions set to n = %d\n", nsam);
00908 exit(1);
00909 }
00910
00911
00912 if (n%2 == 0) {
00913 for (i=n+1;i>=3;i=i-2) {
00914 x(i) = x((i+1)/2);
00915 x(i+1) = work((i+1)/2);
00916 }
00917 x(2) = 0.0;
00918 x(n+2) = 0.0;
00919 }
00920 else {
00921 for (i=n;i>=3;i=i-2) {
00922 x(i) = x(i/2+1);
00923 x(i+1) = work(i/2+1);
00924 }
00925 x(2)=0.0;
00926 }
00927 return status;
00928 }
00929
00930
00931 int Nativefft::fmrs_1rb(float *x, float *work, int nsam)
00932 {
00933
00934
00935
00936 int i, n, status;
00937 int inv=-1;
00938
00939 n = nsam;
00940
00941 for (i=2;i<=n/2+1;i++) {
00942 work(i) = x(2*i)/n;
00943 work(n-i+2) = -work(i);
00944 }
00945 work(1) = 0.0;
00946
00947 for (i=1;i<=n/2+1;i++) x(i) = x(2*i-1)/n;
00948 for (i=n;i>=n/2+2;i--) x(i) = x(n-i+2);
00949
00950 status = fftmcf_(x,work,&n,&n,&n,&inv);
00951
00952 if (status == -1) {
00953 fprintf(stderr, "Native IFT cannot be performed on images with one of ");
00954 fprintf(stderr, "the dimensions set to n = %d\n", nsam);
00955 exit(1);
00956 }
00957
00958 return status;
00959 }
00960
00961 #undef x
00962 #undef y
00963 #undef work
00964
00965
00966
00967
00968 #define x(i,j) x[((j)-1)*lda + (i)-1]
00969 #define y(i,j) y[((j)-1)*lda + (i)-1]
00970
00971
00972 int Nativefft::fmrs_2rf(float *x, float *work, int lda, int nsam, int nrow)
00973 {
00974
00975
00976
00977 int ins, status=0, i, j;
00978
00979 ins = lda;
00980
00981
00982
00983
00984
00985
00986
00987
00988
00989 for (i=1;i<=lda;i=i+2){
00990 status = fftmcf_(&x(i,1),&x(i+1,1),&nrow,&nrow,&nrow,&ins);
00991 if (status == -1) {
00992 fprintf(stderr, "Native FFT cannot be performed on images with one of ");
00993 fprintf(stderr, "the dimensions set to n = %d\n", nrow);
00994 exit(1);
00995 }
00996 }
00997
00998 return status;
00999 }
01000
01001
01002
01003 int Nativefft::fmrs_2rb(float *y, float *work, int lda, int nsam, int nrow)
01004 {
01005
01006
01007
01008 int ins, status=0, i, j;
01009 float q;
01010
01011 ins=-lda;
01012
01013 for (i=1; i<=lda; i=i+2) {
01014 fftmcf_(&y(i,1),&y(i+1,1),&nrow,&nrow,&nrow,&ins);
01015 if (status == -1) {
01016 fprintf(stderr, "Native IFT cannot be performed on images with one of ");
01017 fprintf(stderr, "the dimensions set to n = %d\n", nrow);
01018 exit(1);
01019 }
01020 }
01021
01022
01023 q = 1.0/(float)(nrow);
01024 for (j=1; j<=nrow; j++) for (i=1; i<=lda; i++) y(i,j)*=q;
01025
01026
01027 for (j=1; j<=nrow; j++) status = fmrs_1rb(&y(1,j), work, nsam);
01028
01029 return status;
01030 }
01031 #undef x
01032 #undef y
01033
01034
01035
01036 #define complexd(i,j) complexd[(j-1)*lda + i-1]
01037 #define reald(i,j) reald[(j-1)*nsam + i-1]
01038
01039
01040 int Nativefft::ftp_2rf(float *reald, float *complexd, int lda, int nsam, int nrow)
01041 {
01042
01043
01044
01045
01046 int ins, status=0, i, j;
01047
01048 float * work = (float*) malloc((2*nsam+15)*sizeof(float));
01049 if (!work) {
01050 fprintf(stderr,"real_to_complex_nd(2df): failed to allocate work\n");
01051 LOGERR("real_to_complex_nd(2df): failed to allocate work\n");
01052 }
01053
01054 rffti(nsam, work);
01055
01056 for (j=1; j<=nrow; j++) {
01057 memcpy(&complexd(2,j), &reald(1,j), nsam * sizeof(float));
01058 rfftf(nsam, &complexd(2,j), work);
01059 complexd(1,j) = complexd(2,j) ;
01060 complexd(2,j) = 0.0f ;
01061 if (nsam%2 == 0) complexd(nsam+2,j) = 0.0f ;
01062 }
01063
01064 free(work);
01065
01066 ins = lda;
01067
01068 for (i=1; i<=lda; i=i+2) {
01069 status = fftmcf_(&complexd(i,1),&complexd(i+1,1),&nrow,&nrow,&nrow,&ins);
01070 if (status == -1) {
01071 fprintf(stderr, "Native FFT cannot be performed on images with one of ");
01072 fprintf(stderr, "the dimensions set to n = %d\n", nrow);
01073 exit(1);
01074 }
01075 }
01076
01077 return status;
01078 }
01079
01080
01081 int Nativefft::ftp_2rb(float *complexd, int lda, int nsam, int nrow)
01082 {
01083
01084
01085
01086
01087 int ins, status=0, i, j;
01088
01089 ins = -lda;
01090
01091 for (i=1; i<=lda; i=i+2) {
01092 status = fftmcf_(&complexd(i,1),&complexd(i+1,1),&nrow,&nrow,&nrow,&ins);
01093 if (status == -1) {
01094 fprintf(stderr, "Native FFT cannot be performed on images with one of ");
01095 fprintf(stderr, "the dimensions set to n = %d\n", nrow);
01096 exit(1);
01097 }
01098 }
01099
01100
01101 float * work = (float*) malloc((2*nsam+15)*sizeof(float));
01102 if (!work) {
01103 fprintf(stderr,"real_to_complex_nd(2df): failed to allocate work\n");
01104 LOGERR("real_to_complex_nd(2df): failed to allocate work\n");
01105 }
01106
01107 rffti(nsam, work);
01108
01109 for (j=1; j<=nrow; j++) {
01110 memmove(&complexd(2,j), &complexd(3,j), (nsam-1) * sizeof(float));
01111 rfftb(nsam, &complexd(1,j), work);
01112 }
01113 free(work);
01114
01115 float nrm = 1.0f/float(nsam*nrow);
01116 for (int j = 1; j<=nrow; j++) for (int i = 1; i<=nsam; i++) complexd(i,j) *= nrm;
01117
01118 return status;
01119 }
01120 #undef complexd
01121 #undef reald
01122
01123 #define reald(i,j,k) reald[i-1 + (j-1+(k-1)*nrow)*nsam]
01124 #define complexd(i,j,k) complexd[i-1 + (j-1+(k-1)*nrow)*lda]
01125
01126 int Nativefft::ftp_3rf(float *reald, float *complexd, int lda, int nsam, int nrow, int nslice)
01127 {
01128
01129
01130
01131 int ndr, i, j, k, status=0;
01132
01133 for (k=1; k<=nslice;k ++) status = ftp_2rf(&reald(1,1,k), &complexd(1,1,k), lda, nsam, nrow);
01134
01135 ndr=lda*nrow;
01136
01137 for (j=1; j<=nrow; j++) {
01138 for (i=1; i<=lda-1; i=i+2) {
01139 status = fftmcf_(&complexd(i,j,1),&complexd(i+1,j,1),&nslice,&nslice,&nslice,&ndr);
01140 if (status == -1) {
01141 fprintf(stderr, "Native FFT cannot be performed on images with one of ");
01142 fprintf(stderr, "the dimensions set to n = %d\n", nslice);
01143 exit(1);
01144 }
01145 }
01146 }
01147
01148 return status;
01149 }
01150 #undef reald
01151
01152
01153 int Nativefft::ftp_3rb(float *complexd, int lda, int nsam, int nrow, int nslice)
01154 {
01155
01156
01157
01158 int ndr, i, j, k, status=0;
01159 float q;
01160
01161 ndr=-lda*nrow;
01162
01163 for (j=1; j<=nrow; j++) {
01164 for (i=1; i<=lda-1; i=i+2) {
01165 status = fftmcf_(&complexd(i,j,1),&complexd(i+1,j,1),&nslice,&nslice,&nslice,&ndr);
01166 if (status == -1) {
01167 fprintf(stderr, "Native IFT cannot be performed on images with one of ");
01168 fprintf(stderr, "the dimensions set to n = %d\n", nslice);
01169 exit(1);
01170 }
01171 }
01172 }
01173
01174
01175 q=1.0/(float)(nslice);
01176 for (k=1; k<=nslice; k++) for (j=1;j<=nrow;j++) for (i=1;i<=lda;i++) complexd(i,j,k) *= q;
01177
01178 for (k=1; k<=nslice; k++) status = ftp_2rb(&complexd(1,1,k), lda, nsam, nrow);
01179
01180 return status;
01181 }
01182 #undef complexd
01183
01184
01185 #define b(i,j,k) b[((k)-1)*lda*nrow + ((j)-1)*lda + (i) - 1]
01186
01187
01188 int Nativefft::fmrs_3rf(float *b, float *work, int lda, int nsam, int nrow, int nslice)
01189 {
01190
01191
01192
01193 int ndr, i, j, k, status=0;
01194
01195 ndr=lda*nrow;
01196
01197 for (k=1; k<=nslice; k++) status = fmrs_2rf(&b(1,1,k), work, lda, nsam, nrow);
01198
01199 for (j=1; j<=nrow; j++) {
01200 for (i=1; i<=lda-1; i=i+2) {
01201 status = fftmcf_(&b(i,j,1),&b(i+1,j,1),&nslice,&nslice,&nslice,&ndr);
01202 if (status == -1) {
01203 fprintf(stderr, "Native FFT cannot be performed on images with one of ");
01204 fprintf(stderr, "the dimensions set to n = %d\n", nslice);
01205 exit(1);
01206 }
01207 }
01208 }
01209
01210
01211 return status;
01212 }
01213
01214
01215 int Nativefft::fmrs_3rb(float *b, float *work, int lda, int nsam, int nrow, int nslice)
01216 {
01217
01218
01219
01220 int ndr, i, j, k, status=0;
01221 float q;
01222
01223 ndr=-lda*nrow;
01224
01225 for (j=1;j<=nrow;j++) {
01226 for (i=1;i<=lda-1;i=i+2) {
01227 status = fftmcf_(&b(i,j,1),&b(i+1,j,1),&nslice,&nslice,&nslice,&ndr);
01228 if (status == -1) {
01229 fprintf(stderr, "Native IFT cannot be performed on images with one of ");
01230 fprintf(stderr, "the dimensions set to n = %d\n", nslice);
01231 exit(1);
01232 }
01233 }
01234 }
01235
01236
01237
01238
01239 q=1.0/(float)(nslice);
01240 for (k=1;k<=nslice;k++) for (j=1;j<=nrow;j++) for (i=1;i<=lda;i++) b(i,j,k)*=q;
01241
01242 for (k=1; k<=nslice; k++) status = fmrs_2rb(&b(1,1,k), work, lda, nsam, nrow);
01243 return status;
01244 }
01245 #undef b
01246
01247
01248
01249
01250
01251
01252
01253
01254
01255
01256
01257 #include <math.h>
01258 #include <stdio.h>
01259
01260
01261 #ifdef DOUBLE
01262 #define Treal double
01263 #else
01264 #define Treal float
01265 #endif
01266
01267
01268 #define ref(u,a) u[a]
01269
01270 #define MAXFAC 13
01271 #define NSPECIAL 4
01272
01273
01274
01275
01276
01277
01278
01279
01280
01281
01282
01283 static void passf2(int ido, int l1, const Treal cc[], Treal ch[], const Treal wa1[], int isign)
01284
01285 {
01286 int i, k, ah, ac;
01287 Treal ti2, tr2;
01288 if (ido <= 2) {
01289 for (k=0; k<l1; k++) {
01290 ah = k*ido;
01291 ac = 2*k*ido;
01292 ch[ah] = ref(cc,ac) + ref(cc,ac + ido);
01293 ch[ah + ido*l1] = ref(cc,ac) - ref(cc,ac + ido);
01294 ch[ah+1] = ref(cc,ac+1) + ref(cc,ac + ido + 1);
01295 ch[ah + ido*l1 + 1] = ref(cc,ac+1) - ref(cc,ac + ido + 1);
01296 }
01297 } else {
01298 for (k=0; k<l1; k++) {
01299 for (i=0; i<ido-1; i+=2) {
01300 ah = i + k*ido;
01301 ac = i + 2*k*ido;
01302 ch[ah] = ref(cc,ac) + ref(cc,ac + ido);
01303 tr2 = ref(cc,ac) - ref(cc,ac + ido);
01304 ch[ah+1] = ref(cc,ac+1) + ref(cc,ac + 1 + ido);
01305 ti2 = ref(cc,ac+1) - ref(cc,ac + 1 + ido);
01306 ch[ah+l1*ido+1] = wa1[i]*ti2 + isign*wa1[i+1]*tr2;
01307 ch[ah+l1*ido] = wa1[i]*tr2 - isign*wa1[i+1]*ti2;
01308 }
01309 }
01310 }
01311 }
01312
01313
01314 static void passf3(int ido, int l1, const Treal cc[], Treal ch[],
01315 const Treal wa1[], const Treal wa2[], int isign)
01316
01317 {
01318 static const Treal taur = -0.5;
01319 static const Treal taui = 0.866025403784439;
01320 int i, k, ac, ah;
01321 Treal ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
01322 if (ido == 2) {
01323 for (k=1; k<=l1; k++) {
01324 ac = (3*k - 2)*ido;
01325 tr2 = ref(cc,ac) + ref(cc,ac + ido);
01326 cr2 = ref(cc,ac - ido) + taur*tr2;
01327 ah = (k - 1)*ido;
01328 ch[ah] = ref(cc,ac - ido) + tr2;
01329
01330 ti2 = ref(cc,ac + 1) + ref(cc,ac + ido + 1);
01331 ci2 = ref(cc,ac - ido + 1) + taur*ti2;
01332 ch[ah + 1] = ref(cc,ac - ido + 1) + ti2;
01333
01334 cr3 = isign*taui*(ref(cc,ac) - ref(cc,ac + ido));
01335 ci3 = isign*taui*(ref(cc,ac + 1) - ref(cc,ac + ido + 1));
01336 ch[ah + l1*ido] = cr2 - ci3;
01337 ch[ah + 2*l1*ido] = cr2 + ci3;
01338 ch[ah + l1*ido + 1] = ci2 + cr3;
01339 ch[ah + 2*l1*ido + 1] = ci2 - cr3;
01340 }
01341 } else {
01342 for (k=1; k<=l1; k++) {
01343 for (i=0; i<ido-1; i+=2) {
01344 ac = i + (3*k - 2)*ido;
01345 tr2 = ref(cc,ac) + ref(cc,ac + ido);
01346 cr2 = ref(cc,ac - ido) + taur*tr2;
01347 ah = i + (k-1)*ido;
01348 ch[ah] = ref(cc,ac - ido) + tr2;
01349 ti2 = ref(cc,ac + 1) + ref(cc,ac + ido + 1);
01350 ci2 = ref(cc,ac - ido + 1) + taur*ti2;
01351 ch[ah + 1] = ref(cc,ac - ido + 1) + ti2;
01352 cr3 = isign*taui*(ref(cc,ac) - ref(cc,ac + ido));
01353 ci3 = isign*taui*(ref(cc,ac + 1) - ref(cc,ac + ido + 1));
01354 dr2 = cr2 - ci3;
01355 dr3 = cr2 + ci3;
01356 di2 = ci2 + cr3;
01357 di3 = ci2 - cr3;
01358 ch[ah + l1*ido + 1] = wa1[i]*di2 + isign*wa1[i+1]*dr2;
01359 ch[ah + l1*ido] = wa1[i]*dr2 - isign*wa1[i+1]*di2;
01360 ch[ah + 2*l1*ido + 1] = wa2[i]*di3 + isign*wa2[i+1]*dr3;
01361 ch[ah + 2*l1*ido] = wa2[i]*dr3 - isign*wa2[i+1]*di3;
01362 }
01363 }
01364 }
01365 }
01366
01367
01368 static void passf4(int ido, int l1, const Treal cc[], Treal ch[],
01369 const Treal wa1[], const Treal wa2[], const Treal wa3[], int isign)
01370
01371 {
01372 int i, k, ac, ah;
01373 Treal ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
01374 if (ido == 2) {
01375 for (k=0; k<l1; k++) {
01376 ac = 4*k*ido + 1;
01377 ti1 = ref(cc,ac) - ref(cc,ac + 2*ido);
01378 ti2 = ref(cc,ac) + ref(cc,ac + 2*ido);
01379 tr4 = ref(cc,ac + 3*ido) - ref(cc,ac + ido);
01380 ti3 = ref(cc,ac + ido) + ref(cc,ac + 3*ido);
01381 tr1 = ref(cc,ac - 1) - ref(cc,ac + 2*ido - 1);
01382 tr2 = ref(cc,ac - 1) + ref(cc,ac + 2*ido - 1);
01383 ti4 = ref(cc,ac + ido - 1) - ref(cc,ac + 3*ido - 1);
01384 tr3 = ref(cc,ac + ido - 1) + ref(cc,ac + 3*ido - 1);
01385 ah = k*ido;
01386 ch[ah] = tr2 + tr3;
01387 ch[ah + 2*l1*ido] = tr2 - tr3;
01388 ch[ah + 1] = ti2 + ti3;
01389 ch[ah + 2*l1*ido + 1] = ti2 - ti3;
01390 ch[ah + l1*ido] = tr1 + isign*tr4;
01391 ch[ah + 3*l1*ido] = tr1 - isign*tr4;
01392 ch[ah + l1*ido + 1] = ti1 + isign*ti4;
01393 ch[ah + 3*l1*ido + 1] = ti1 - isign*ti4;
01394 }
01395 } else {
01396 for (k=0; k<l1; k++) {
01397 for (i=0; i<ido-1; i+=2) {
01398 ac = i + 1 + 4*k*ido;
01399 ti1 = ref(cc,ac) - ref(cc,ac + 2*ido);
01400 ti2 = ref(cc,ac) + ref(cc,ac + 2*ido);
01401 ti3 = ref(cc,ac + ido) + ref(cc,ac + 3*ido);
01402 tr4 = ref(cc,ac + 3*ido) - ref(cc,ac + ido);
01403 tr1 = ref(cc,ac - 1) - ref(cc,ac + 2*ido - 1);
01404 tr2 = ref(cc,ac - 1) + ref(cc,ac + 2*ido - 1);
01405 ti4 = ref(cc,ac + ido - 1) - ref(cc,ac + 3*ido - 1);
01406 tr3 = ref(cc,ac + ido - 1) + ref(cc,ac + 3*ido - 1);
01407 ah = i + k*ido;
01408 ch[ah] = tr2 + tr3;
01409 cr3 = tr2 - tr3;
01410 ch[ah + 1] = ti2 + ti3;
01411 ci3 = ti2 - ti3;
01412 cr2 = tr1 + isign*tr4;
01413 cr4 = tr1 - isign*tr4;
01414 ci2 = ti1 + isign*ti4;
01415 ci4 = ti1 - isign*ti4;
01416 ch[ah + l1*ido] = wa1[i]*cr2 - isign*wa1[i + 1]*ci2;
01417 ch[ah + l1*ido + 1] = wa1[i]*ci2 + isign*wa1[i + 1]*cr2;
01418 ch[ah + 2*l1*ido] = wa2[i]*cr3 - isign*wa2[i + 1]*ci3;
01419 ch[ah + 2*l1*ido + 1] = wa2[i]*ci3 + isign*wa2[i + 1]*cr3;
01420 ch[ah + 3*l1*ido] = wa3[i]*cr4 -isign*wa3[i + 1]*ci4;
01421 ch[ah + 3*l1*ido + 1] = wa3[i]*ci4 + isign*wa3[i + 1]*cr4;
01422 }
01423 }
01424 }
01425 }
01426
01427
01428 static void passf5(int ido, int l1, const Treal cc[], Treal ch[],
01429 const Treal wa1[], const Treal wa2[], const Treal wa3[], const Treal wa4[], int isign)
01430
01431 {
01432 static const Treal tr11 = 0.309016994374947;
01433 static const Treal ti11 = 0.951056516295154;
01434 static const Treal tr12 = -0.809016994374947;
01435 static const Treal ti12 = 0.587785252292473;
01436 int i, k, ac, ah;
01437 Treal ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4, ti2, ti3,
01438 ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;
01439 if (ido == 2) {
01440 for (k = 1; k <= l1; ++k) {
01441 ac = (5*k - 4)*ido + 1;
01442 ti5 = ref(cc,ac) - ref(cc,ac + 3*ido);
01443 ti2 = ref(cc,ac) + ref(cc,ac + 3*ido);
01444 ti4 = ref(cc,ac + ido) - ref(cc,ac + 2*ido);
01445 ti3 = ref(cc,ac + ido) + ref(cc,ac + 2*ido);
01446 tr5 = ref(cc,ac - 1) - ref(cc,ac + 3*ido - 1);
01447 tr2 = ref(cc,ac - 1) + ref(cc,ac + 3*ido - 1);
01448 tr4 = ref(cc,ac + ido - 1) - ref(cc,ac + 2*ido - 1);
01449 tr3 = ref(cc,ac + ido - 1) + ref(cc,ac + 2*ido - 1);
01450 ah = (k - 1)*ido;
01451 ch[ah] = ref(cc,ac - ido - 1) + tr2 + tr3;
01452 ch[ah + 1] = ref(cc,ac - ido) + ti2 + ti3;
01453 cr2 = ref(cc,ac - ido - 1) + tr11*tr2 + tr12*tr3;
01454 ci2 = ref(cc,ac - ido) + tr11*ti2 + tr12*ti3;
01455 cr3 = ref(cc,ac - ido - 1) + tr12*tr2 + tr11*tr3;
01456 ci3 = ref(cc,ac - ido) + tr12*ti2 + tr11*ti3;
01457 cr5 = isign*(ti11*tr5 + ti12*tr4);
01458 ci5 = isign*(ti11*ti5 + ti12*ti4);
01459 cr4 = isign*(ti12*tr5 - ti11*tr4);
01460 ci4 = isign*(ti12*ti5 - ti11*ti4);
01461 ch[ah + l1*ido] = cr2 - ci5;
01462 ch[ah + 4*l1*ido] = cr2 + ci5;
01463 ch[ah + l1*ido + 1] = ci2 + cr5;
01464 ch[ah + 2*l1*ido + 1] = ci3 + cr4;
01465 ch[ah + 2*l1*ido] = cr3 - ci4;
01466 ch[ah + 3*l1*ido] = cr3 + ci4;
01467 ch[ah + 3*l1*ido + 1] = ci3 - cr4;
01468 ch[ah + 4*l1*ido + 1] = ci2 - cr5;
01469 }
01470 } else {
01471 for (k=1; k<=l1; k++) {
01472 for (i=0; i<ido-1; i+=2) {
01473 ac = i + 1 + (k*5 - 4)*ido;
01474 ti5 = ref(cc,ac) - ref(cc,ac + 3*ido);
01475 ti2 = ref(cc,ac) + ref(cc,ac + 3*ido);
01476 ti4 = ref(cc,ac + ido) - ref(cc,ac + 2*ido);
01477 ti3 = ref(cc,ac + ido) + ref(cc,ac + 2*ido);
01478 tr5 = ref(cc,ac - 1) - ref(cc,ac + 3*ido - 1);
01479 tr2 = ref(cc,ac - 1) + ref(cc,ac + 3*ido - 1);
01480 tr4 = ref(cc,ac + ido - 1) - ref(cc,ac + 2*ido - 1);
01481 tr3 = ref(cc,ac + ido - 1) + ref(cc,ac + 2*ido - 1);
01482 ah = i + (k - 1)*ido;
01483 ch[ah] = ref(cc,ac - ido - 1) + tr2 + tr3;
01484 ch[ah + 1] = ref(cc,ac - ido) + ti2 + ti3;
01485 cr2 = ref(cc,ac - ido - 1) + tr11*tr2 + tr12*tr3;
01486
01487 ci2 = ref(cc,ac - ido) + tr11*ti2 + tr12*ti3;
01488 cr3 = ref(cc,ac - ido - 1) + tr12*tr2 + tr11*tr3;
01489
01490 ci3 = ref(cc,ac - ido) + tr12*ti2 + tr11*ti3;
01491 cr5 = isign*(ti11*tr5 + ti12*tr4);
01492 ci5 = isign*(ti11*ti5 + ti12*ti4);
01493 cr4 = isign*(ti12*tr5 - ti11*tr4);
01494 ci4 = isign*(ti12*ti5 - ti11*ti4);
01495 dr3 = cr3 - ci4;
01496 dr4 = cr3 + ci4;
01497 di3 = ci3 + cr4;
01498 di4 = ci3 - cr4;
01499 dr5 = cr2 + ci5;
01500 dr2 = cr2 - ci5;
01501 di5 = ci2 - cr5;
01502 di2 = ci2 + cr5;
01503 ch[ah + l1*ido] = wa1[i]*dr2 - isign*wa1[i+1]*di2;
01504 ch[ah + l1*ido + 1] = wa1[i]*di2 + isign*wa1[i+1]*dr2;
01505 ch[ah + 2*l1*ido] = wa2[i]*dr3 - isign*wa2[i+1]*di3;
01506 ch[ah + 2*l1*ido + 1] = wa2[i]*di3 + isign*wa2[i+1]*dr3;
01507 ch[ah + 3*l1*ido] = wa3[i]*dr4 - isign*wa3[i+1]*di4;
01508 ch[ah + 3*l1*ido + 1] = wa3[i]*di4 + isign*wa3[i+1]*dr4;
01509 ch[ah + 4*l1*ido] = wa4[i]*dr5 - isign*wa4[i+1]*di5;
01510 ch[ah + 4*l1*ido + 1] = wa4[i]*di5 + isign*wa4[i+1]*dr5;
01511 }
01512 }
01513 }
01514 }
01515
01516
01517 static void passf(int *nac, int ido, int ip, int l1, int idl1,
01518 Treal cc[], Treal ch[],
01519 const Treal wa[], int isign)
01520
01521 {
01522 int idij, idlj, idot, ipph, i, j, k, l, jc, lc, ik, nt, idj, idl, inc,idp;
01523 Treal wai, war;
01524
01525 idot = ido / 2;
01526 nt = ip*idl1;
01527 ipph = (ip + 1) / 2;
01528 idp = ip*ido;
01529 if (ido >= l1) {
01530 for (j=1; j<ipph; j++) {
01531 jc = ip - j;
01532 for (k=0; k<l1; k++) {
01533 for (i=0; i<ido; i++) {
01534 ch[i + (k + j*l1)*ido] =
01535 ref(cc,i + (j + k*ip)*ido) + ref(cc,i + (jc + k*ip)*ido);
01536 ch[i + (k + jc*l1)*ido] =
01537 ref(cc,i + (j + k*ip)*ido) - ref(cc,i + (jc + k*ip)*ido);
01538 }
01539 }
01540 }
01541 for (k=0; k<l1; k++)
01542 for (i=0; i<ido; i++)
01543 ch[i + k*ido] = ref(cc,i + k*ip*ido);
01544 } else {
01545 for (j=1; j<ipph; j++) {
01546 jc = ip - j;
01547 for (i=0; i<ido; i++) {
01548 for (k=0; k<l1; k++) {
01549 ch[i + (k + j*l1)*ido] = ref(cc,i + (j + k*ip)*ido) + ref(cc,i + (jc + k*
01550 ip)*ido);
01551 ch[i + (k + jc*l1)*ido] = ref(cc,i + (j + k*ip)*ido) - ref(cc,i + (jc + k*
01552 ip)*ido);
01553 }
01554 }
01555 }
01556 for (i=0; i<ido; i++)
01557 for (k=0; k<l1; k++)
01558 ch[i + k*ido] = ref(cc,i + k*ip*ido);
01559 }
01560
01561 idl = 2 - ido;
01562 inc = 0;
01563 for (l=1; l<ipph; l++) {
01564 lc = ip - l;
01565 idl += ido;
01566 for (ik=0; ik<idl1; ik++) {
01567 cc[ik + l*idl1] = ch[ik] + wa[idl - 2]*ch[ik + idl1];
01568 cc[ik + lc*idl1] = isign*wa[idl-1]*ch[ik + (ip-1)*idl1];
01569 }
01570 idlj = idl;
01571 inc += ido;
01572 for (j=2; j<ipph; j++) {
01573 jc = ip - j;
01574 idlj += inc;
01575 if (idlj > idp) idlj -= idp;
01576 war = wa[idlj - 2];
01577 wai = wa[idlj-1];
01578 for (ik=0; ik<idl1; ik++) {
01579 cc[ik + l*idl1] += war*ch[ik + j*idl1];
01580 cc[ik + lc*idl1] += isign*wai*ch[ik + jc*idl1];
01581 }
01582 }
01583 }
01584 for (j=1; j<ipph; j++)
01585 for (ik=0; ik<idl1; ik++)
01586 ch[ik] += ch[ik + j*idl1];
01587 for (j=1; j<ipph; j++) {
01588 jc = ip - j;
01589 for (ik=1; ik<idl1; ik+=2) {
01590 ch[ik - 1 + j*idl1] = cc[ik - 1 + j*idl1] - cc[ik + jc*idl1];
01591 ch[ik - 1 + jc*idl1] = cc[ik - 1 + j*idl1] + cc[ik + jc*idl1];
01592 ch[ik + j*idl1] = cc[ik + j*idl1] + cc[ik - 1 + jc*idl1];
01593 ch[ik + jc*idl1] = cc[ik + j*idl1] - cc[ik - 1 + jc*idl1];
01594 }
01595 }
01596 *nac = 1;
01597 if (ido == 2) return;
01598 *nac = 0;
01599 for (ik=0; ik<idl1; ik++)
01600 cc[ik] = ch[ik];
01601 for (j=1; j<ip; j++) {
01602 for (k=0; k<l1; k++) {
01603 cc[(k + j*l1)*ido + 0] = ch[(k + j*l1)*ido + 0];
01604 cc[(k + j*l1)*ido + 1] = ch[(k + j*l1)*ido + 1];
01605 }
01606 }
01607 if (idot <= l1) {
01608 idij = 0;
01609 for (j=1; j<ip; j++) {
01610 idij += 2;
01611 for (i=3; i<ido; i+=2) {
01612 idij += 2;
01613 for (k=0; k<l1; k++) {
01614 cc[i - 1 + (k + j*l1)*ido] =
01615 wa[idij - 2]*ch[i - 1 + (k + j*l1)*ido] -
01616 isign*wa[idij-1]*ch[i + (k + j*l1)*ido];
01617 cc[i + (k + j*l1)*ido] =
01618 wa[idij - 2]*ch[i + (k + j*l1)*ido] +
01619 isign*wa[idij-1]*ch[i - 1 + (k + j*l1)*ido];
01620 }
01621 }
01622 }
01623 } else {
01624 idj = 2 - ido;
01625 for (j=1; j<ip; j++) {
01626 idj += ido;
01627 for (k = 0; k < l1; k++) {
01628 idij = idj;
01629 for (i=3; i<ido; i+=2) {
01630 idij += 2;
01631 cc[i - 1 + (k + j*l1)*ido] =
01632 wa[idij - 2]*ch[i - 1 + (k + j*l1)*ido] -
01633 isign*wa[idij-1]*ch[i + (k + j*l1)*ido];
01634 cc[i + (k + j*l1)*ido] =
01635 wa[idij - 2]*ch[i + (k + j*l1)*ido] +
01636 isign*wa[idij-1]*ch[i - 1 + (k + j*l1)*ido];
01637 }
01638 }
01639 }
01640 }
01641 }
01642
01643
01644
01645
01646
01647
01648
01649 static void radf2(int ido, int l1, const Treal cc[], Treal ch[], const Treal wa1[])
01650 {
01651 int i, k, ic;
01652 Treal ti2, tr2;
01653 for (k=0; k<l1; k++) {
01654 ch[2*k*ido] =
01655 ref(cc,k*ido) + ref(cc,(k + l1)*ido);
01656 ch[(2*k+1)*ido + ido-1] =
01657 ref(cc,k*ido) - ref(cc,(k + l1)*ido);
01658 }
01659 if (ido < 2) return;
01660 if (ido != 2) {
01661 for (k=0; k<l1; k++) {
01662 for (i=2; i<ido; i+=2) {
01663 ic = ido - i;
01664 tr2 = wa1[i - 2]*ref(cc, i-1 + (k + l1)*ido) + wa1[i - 1]*ref(cc, i + (k + l1)*ido);
01665 ti2 = wa1[i - 2]*ref(cc, i + (k + l1)*ido) - wa1[i - 1]*ref(cc, i-1 + (k + l1)*ido);
01666 ch[i + 2*k*ido] = ref(cc,i + k*ido) + ti2;
01667 ch[ic + (2*k+1)*ido] = ti2 - ref(cc,i + k*ido);
01668 ch[i - 1 + 2*k*ido] = ref(cc,i - 1 + k*ido) + tr2;
01669 ch[ic - 1 + (2*k+1)*ido] = ref(cc,i - 1 + k*ido) - tr2;
01670 }
01671 }
01672 if (ido % 2 == 1) return;
01673 }
01674 for (k=0; k<l1; k++) {
01675 ch[(2*k+1)*ido] = -ref(cc,ido-1 + (k + l1)*ido);
01676 ch[ido-1 + 2*k*ido] = ref(cc,ido-1 + k*ido);
01677 }
01678 }
01679
01680
01681 static void radb2(int ido, int l1, const Treal cc[], Treal ch[], const Treal wa1[])
01682 {
01683 int i, k, ic;
01684 Treal ti2, tr2;
01685 for (k=0; k<l1; k++) {
01686 ch[k*ido] =
01687 ref(cc,2*k*ido) + ref(cc,ido-1 + (2*k+1)*ido);
01688 ch[(k + l1)*ido] =
01689 ref(cc,2*k*ido) - ref(cc,ido-1 + (2*k+1)*ido);
01690 }
01691 if (ido < 2) return;
01692 if (ido != 2) {
01693 for (k = 0; k < l1; ++k) {
01694 for (i = 2; i < ido; i += 2) {
01695 ic = ido - i;
01696 ch[i-1 + k*ido] =
01697 ref(cc,i-1 + 2*k*ido) + ref(cc,ic-1 + (2*k+1)*ido);
01698 tr2 = ref(cc,i-1 + 2*k*ido) - ref(cc,ic-1 + (2*k+1)*ido);
01699 ch[i + k*ido] =
01700 ref(cc,i + 2*k*ido) - ref(cc,ic + (2*k+1)*ido);
01701 ti2 = ref(cc,i + (2*k)*ido) + ref(cc,ic + (2*k+1)*ido);
01702 ch[i-1 + (k + l1)*ido] =
01703 wa1[i - 2]*tr2 - wa1[i - 1]*ti2;
01704 ch[i + (k + l1)*ido] =
01705 wa1[i - 2]*ti2 + wa1[i - 1]*tr2;
01706 }
01707 }
01708 if (ido % 2 == 1) return;
01709 }
01710 for (k = 0; k < l1; k++) {
01711 ch[ido-1 + k*ido] = 2*ref(cc,ido-1 + 2*k*ido);
01712 ch[ido-1 + (k + l1)*ido] = -2*ref(cc,(2*k+1)*ido);
01713 }
01714 }
01715
01716
01717 static void radf3(int ido, int l1, const Treal cc[], Treal ch[],
01718 const Treal wa1[], const Treal wa2[])
01719 {
01720 static const Treal taur = -0.5;
01721 static const Treal taui = 0.866025403784439;
01722 int i, k, ic;
01723 Treal ci2, di2, di3, cr2, dr2, dr3, ti2, ti3, tr2, tr3;
01724 for (k=0; k<l1; k++) {
01725 cr2 = ref(cc,(k + l1)*ido) + ref(cc,(k + 2*l1)*ido);
01726 ch[3*k*ido] = ref(cc,k*ido) + cr2;
01727 ch[(3*k+2)*ido] = taui*(ref(cc,(k + l1*2)*ido) - ref(cc,(k + l1)*ido));
01728 ch[ido-1 + (3*k + 1)*ido] = ref(cc,k*ido) + taur*cr2;
01729 }
01730 if (ido == 1) return;
01731 for (k=0; k<l1; k++) {
01732 for (i=2; i<ido; i+=2) {
01733 ic = ido - i;
01734 dr2 = wa1[i - 2]*ref(cc,i - 1 + (k + l1)*ido) +
01735 wa1[i - 1]*ref(cc,i + (k + l1)*ido);
01736 di2 = wa1[i - 2]*ref(cc,i + (k + l1)*ido) - wa1[i - 1]*ref(cc,i - 1 + (k + l1)*ido);
01737 dr3 = wa2[i - 2]*ref(cc,i - 1 + (k + l1*2)*ido) + wa2[i - 1]*ref(cc,i + (k + l1*2)*ido);
01738 di3 = wa2[i - 2]*ref(cc,i + (k + l1*2)*ido) - wa2[i - 1]*ref(cc,i - 1 + (k + l1*2)*ido);
01739 cr2 = dr2 + dr3;
01740 ci2 = di2 + di3;
01741 ch[i - 1 + 3*k*ido] = ref(cc,i - 1 + k*ido) + cr2;
01742 ch[i + 3*k*ido] = ref(cc,i + k*ido) + ci2;
01743 tr2 = ref(cc,i - 1 + k*ido) + taur*cr2;
01744 ti2 = ref(cc,i + k*ido) + taur*ci2;
01745 tr3 = taui*(di2 - di3);
01746 ti3 = taui*(dr3 - dr2);
01747 ch[i - 1 + (3*k + 2)*ido] = tr2 + tr3;
01748 ch[ic - 1 + (3*k + 1)*ido] = tr2 - tr3;
01749 ch[i + (3*k + 2)*ido] = ti2 + ti3;
01750 ch[ic + (3*k + 1)*ido] = ti3 - ti2;
01751 }
01752 }
01753 }
01754
01755
01756 static void radb3(int ido, int l1, const Treal cc[], Treal ch[],
01757 const Treal wa1[], const Treal wa2[])
01758 {
01759 static const Treal taur = -0.5;
01760 static const Treal taui = 0.866025403784439;
01761 int i, k, ic;
01762 Treal ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
01763 for (k=0; k<l1; k++) {
01764 tr2 = 2*ref(cc,ido-1 + (3*k + 1)*ido);
01765 cr2 = ref(cc,3*k*ido) + taur*tr2;
01766 ch[k*ido] = ref(cc,3*k*ido) + tr2;
01767 ci3 = 2*taui*ref(cc,(3*k + 2)*ido);
01768 ch[(k + l1)*ido] = cr2 - ci3;
01769 ch[(k + 2*l1)*ido] = cr2 + ci3;
01770 }
01771 if (ido == 1) return;
01772 for (k=0; k<l1; k++) {
01773 for (i=2; i<ido; i+=2) {
01774 ic = ido - i;
01775 tr2 = ref(cc,i - 1 + (3*k + 2)*ido) + ref(cc,ic - 1 + (3*k + 1)*ido);
01776 cr2 = ref(cc,i - 1 + 3*k*ido) + taur*tr2;
01777 ch[i - 1 + k*ido] = ref(cc,i - 1 + 3*k*ido) + tr2;
01778 ti2 = ref(cc,i + (3*k + 2)*ido) - ref(cc,ic + (3*k + 1)*ido);
01779 ci2 = ref(cc,i + 3*k*ido) + taur*ti2;
01780 ch[i + k*ido] = ref(cc,i + 3*k*ido) + ti2;
01781 cr3 = taui*(ref(cc,i - 1 + (3*k + 2)*ido) - ref(cc,ic - 1 + (3*k + 1)*ido));
01782 ci3 = taui*(ref(cc,i + (3*k + 2)*ido) + ref(cc,ic + (3*k + 1)*ido));
01783 dr2 = cr2 - ci3;
01784 dr3 = cr2 + ci3;
01785 di2 = ci2 + cr3;
01786 di3 = ci2 - cr3;
01787 ch[i - 1 + (k + l1)*ido] = wa1[i - 2]*dr2 - wa1[i - 1]*di2;
01788 ch[i + (k + l1)*ido] = wa1[i - 2]*di2 + wa1[i - 1]*dr2;
01789 ch[i - 1 + (k + 2*l1)*ido] = wa2[i - 2]*dr3 - wa2[i - 1]*di3;
01790 ch[i + (k + 2*l1)*ido] = wa2[i - 2]*di3 + wa2[i - 1]*dr3;
01791 }
01792 }
01793 }
01794
01795
01796 static void radf4(int ido, int l1, const Treal cc[], Treal ch[],
01797 const Treal wa1[], const Treal wa2[], const Treal wa3[])
01798 {
01799 static const Treal hsqt2 = 0.7071067811865475;
01800 int i, k, ic;
01801 Treal ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
01802 for (k=0; k<l1; k++) {
01803 tr1 = ref(cc,(k + l1)*ido) + ref(cc,(k + 3*l1)*ido);
01804 tr2 = ref(cc,k*ido) + ref(cc,(k + 2*l1)*ido);
01805 ch[4*k*ido] = tr1 + tr2;
01806 ch[ido-1 + (4*k + 3)*ido] = tr2 - tr1;
01807 ch[ido-1 + (4*k + 1)*ido] = ref(cc,k*ido) - ref(cc,(k + 2*l1)*ido);
01808 ch[(4*k + 2)*ido] = ref(cc,(k + 3*l1)*ido) - ref(cc,(k + l1)*ido);
01809 }
01810 if (ido < 2) return;
01811 if (ido != 2) {
01812 for (k=0; k<l1; k++) {
01813 for (i=2; i<ido; i += 2) {
01814 ic = ido - i;
01815 cr2 = wa1[i - 2]*ref(cc,i - 1 + (k + l1)*ido) + wa1[i - 1]*ref(cc,i + (k + l1)*ido);
01816 ci2 = wa1[i - 2]*ref(cc,i + (k + l1)*ido) - wa1[i - 1]*ref(cc,i - 1 + (k + l1)*ido);
01817 cr3 = wa2[i - 2]*ref(cc,i - 1 + (k + 2*l1)*ido) + wa2[i - 1]*ref(cc,i + (k + 2*l1)*
01818 ido);
01819 ci3 = wa2[i - 2]*ref(cc,i + (k + 2*l1)*ido) - wa2[i - 1]*ref(cc,i - 1 + (k + 2*l1)*
01820 ido);
01821 cr4 = wa3[i - 2]*ref(cc,i - 1 + (k + 3*l1)*ido) + wa3[i - 1]*ref(cc,i + (k + 3*l1)*
01822 ido);
01823 ci4 = wa3[i - 2]*ref(cc,i + (k + 3*l1)*ido) - wa3[i - 1]*ref(cc,i - 1 + (k + 3*l1)*
01824 ido);
01825 tr1 = cr2 + cr4;
01826 tr4 = cr4 - cr2;
01827 ti1 = ci2 + ci4;
01828 ti4 = ci2 - ci4;
01829 ti2 = ref(cc,i + k*ido) + ci3;
01830 ti3 = ref(cc,i + k*ido) - ci3;
01831 tr2 = ref(cc,i - 1 + k*ido) + cr3;
01832 tr3 = ref(cc,i - 1 + k*ido) - cr3;
01833 ch[i - 1 + 4*k*ido] = tr1 + tr2;
01834 ch[ic - 1 + (4*k + 3)*ido] = tr2 - tr1;
01835 ch[i + 4*k*ido] = ti1 + ti2;
01836 ch[ic + (4*k + 3)*ido] = ti1 - ti2;
01837 ch[i - 1 + (4*k + 2)*ido] = ti4 + tr3;
01838 ch[ic - 1 + (4*k + 1)*ido] = tr3 - ti4;
01839 ch[i + (4*k + 2)*ido] = tr4 + ti3;
01840 ch[ic + (4*k + 1)*ido] = tr4 - ti3;
01841 }
01842 }
01843 if (ido % 2 == 1) return;
01844 }
01845 for (k=0; k<l1; k++) {
01846 ti1 = -hsqt2*(ref(cc,ido-1 + (k + l1)*ido) + ref(cc,ido-1 + (k + 3*l1)*ido));
01847 tr1 = hsqt2*(ref(cc,ido-1 + (k + l1)*ido) - ref(cc,ido-1 + (k + 3*l1)*ido));
01848 ch[ido-1 + 4*k*ido] = tr1 + ref(cc,ido-1 + k*ido);
01849 ch[ido-1 + (4*k + 2)*ido] = ref(cc,ido-1 + k*ido) - tr1;
01850 ch[(4*k + 1)*ido] = ti1 - ref(cc,ido-1 + (k + 2*l1)*ido);
01851 ch[(4*k + 3)*ido] = ti1 + ref(cc,ido-1 + (k + 2*l1)*ido);
01852 }
01853 }
01854
01855
01856 static void radb4(int ido, int l1, const Treal cc[], Treal ch[],
01857 const Treal wa1[], const Treal wa2[], const Treal wa3[])
01858 {
01859 static const Treal sqrt2 = 1.414213562373095;
01860 int i, k, ic;
01861 Treal ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
01862 for (k = 0; k < l1; k++) {
01863 tr1 = ref(cc,4*k*ido) - ref(cc,ido-1 + (4*k + 3)*ido);
01864 tr2 = ref(cc,4*k*ido) + ref(cc,ido-1 + (4*k + 3)*ido);
01865 tr3 = ref(cc,ido-1 + (4*k + 1)*ido) + ref(cc,ido-1 + (4*k + 1)*ido);
01866 tr4 = ref(cc,(4*k + 2)*ido) + ref(cc,(4*k + 2)*ido);
01867 ch[k*ido] = tr2 + tr3;
01868 ch[(k + l1)*ido] = tr1 - tr4;
01869 ch[(k + 2*l1)*ido] = tr2 - tr3;
01870 ch[(k + 3*l1)*ido] = tr1 + tr4;
01871 }
01872 if (ido < 2) return;
01873 if (ido != 2) {
01874 for (k = 0; k < l1; ++k) {
01875 for (i = 2; i < ido; i += 2) {
01876 ic = ido - i;
01877 ti1 = ref(cc,i + 4*k*ido) + ref(cc,ic + (4*k + 3)*ido);
01878 ti2 = ref(cc,i + 4*k*ido) - ref(cc,ic + (4*k + 3)*ido);
01879 ti3 = ref(cc,i + (4*k + 2)*ido) - ref(cc,ic + (4*k + 1)*ido);
01880 tr4 = ref(cc,i + (4*k + 2)*ido) + ref(cc,ic + (4*k + 1)*ido);
01881 tr1 = ref(cc,i - 1 + 4*k*ido) - ref(cc,ic - 1 + (4*k + 3)*ido);
01882 tr2 = ref(cc,i - 1 + 4*k*ido) + ref(cc,ic - 1 + (4*k + 3)*ido);
01883 ti4 = ref(cc,i - 1 + (4*k + 2)*ido) - ref(cc,ic - 1 + (4*k + 1)*ido);
01884 tr3 = ref(cc,i - 1 + (4*k + 2)*ido) + ref(cc,ic - 1 + (4*k + 1)*ido);
01885 ch[i - 1 + k*ido] = tr2 + tr3;
01886 cr3 = tr2 - tr3;
01887 ch[i + k*ido] = ti2 + ti3;
01888 ci3 = ti2 - ti3;
01889 cr2 = tr1 - tr4;
01890 cr4 = tr1 + tr4;
01891 ci2 = ti1 + ti4;
01892 ci4 = ti1 - ti4;
01893 ch[i - 1 + (k + l1)*ido] = wa1[i - 2]*cr2 - wa1[i - 1]*ci2;
01894 ch[i + (k + l1)*ido] = wa1[i - 2]*ci2 + wa1[i - 1]*cr2;
01895 ch[i - 1 + (k + 2*l1)*ido] = wa2[i - 2]*cr3 - wa2[i - 1]*ci3;
01896 ch[i + (k + 2*l1)*ido] = wa2[i - 2]*ci3 + wa2[i - 1]*cr3;
01897 ch[i - 1 + (k + 3*l1)*ido] = wa3[i - 2]*cr4 - wa3[i - 1]*ci4;
01898 ch[i + (k + 3*l1)*ido] = wa3[i - 2]*ci4 + wa3[i - 1]*cr4;
01899 }
01900 }
01901 if (ido % 2 == 1) return;
01902 }
01903 for (k = 0; k < l1; k++) {
01904 ti1 = ref(cc,(4*k + 1)*ido) + ref(cc,(4*k + 3)*ido);
01905 ti2 = ref(cc,(4*k + 3)*ido) - ref(cc,(4*k + 1)*ido);
01906 tr1 = ref(cc,ido-1 + 4*k*ido) - ref(cc,ido-1 + (4*k + 2)*ido);
01907 tr2 = ref(cc,ido-1 + 4*k*ido) + ref(cc,ido-1 + (4*k + 2)*ido);
01908 ch[ido-1 + k*ido] = tr2 + tr2;
01909 ch[ido-1 + (k + l1)*ido] = sqrt2*(tr1 - ti1);
01910 ch[ido-1 + (k + 2*l1)*ido] = ti2 + ti2;
01911 ch[ido-1 + (k + 3*l1)*ido] = -sqrt2*(tr1 + ti1);
01912 }
01913 }
01914
01915
01916 static void radf5(int ido, int l1, const Treal cc[], Treal ch[],
01917 const Treal wa1[], const Treal wa2[], const Treal wa3[], const Treal wa4[])
01918 {
01919 static const Treal tr11 = 0.309016994374947;
01920 static const Treal ti11 = 0.951056516295154;
01921 static const Treal tr12 = -0.809016994374947;
01922 static const Treal ti12 = 0.587785252292473;
01923 int i, k, ic;
01924 Treal ci2, di2, ci4, ci5, di3, di4, di5, ci3, cr2, cr3, dr2, dr3, dr4, dr5,
01925 cr5, cr4, ti2, ti3, ti5, ti4, tr2, tr3, tr4, tr5;
01926 for (k = 0; k < l1; k++) {
01927 cr2 = ref(cc,(k + 4*l1)*ido) + ref(cc,(k + l1)*ido);
01928 ci5 = ref(cc,(k + 4*l1)*ido) - ref(cc,(k + l1)*ido);
01929 cr3 = ref(cc,(k + 3*l1)*ido) + ref(cc,(k + 2*l1)*ido);
01930 ci4 = ref(cc,(k + 3*l1)*ido) - ref(cc,(k + 2*l1)*ido);
01931 ch[5*k*ido] = ref(cc,k*ido) + cr2 + cr3;
01932 ch[ido-1 + (5*k + 1)*ido] = ref(cc,k*ido) + tr11*cr2 + tr12*cr3;
01933 ch[(5*k + 2)*ido] = ti11*ci5 + ti12*ci4;
01934 ch[ido-1 + (5*k + 3)*ido] = ref(cc,k*ido) + tr12*cr2 + tr11*cr3;
01935 ch[(5*k + 4)*ido] = ti12*ci5 - ti11*ci4;
01936 }
01937 if (ido == 1) return;
01938 for (k = 0; k < l1; ++k) {
01939 for (i = 2; i < ido; i += 2) {
01940 ic = ido - i;
01941 dr2 = wa1[i - 2]*ref(cc,i - 1 + (k + l1)*ido) + wa1[i - 1]*ref(cc,i + (k + l1)*ido);
01942 di2 = wa1[i - 2]*ref(cc,i + (k + l1)*ido) - wa1[i - 1]*ref(cc,i - 1 + (k + l1)*ido);
01943 dr3 = wa2[i - 2]*ref(cc,i - 1 + (k + 2*l1)*ido) + wa2[i - 1]*ref(cc,i + (k + 2*l1)*ido);
01944 di3 = wa2[i - 2]*ref(cc,i + (k + 2*l1)*ido) - wa2[i - 1]*ref(cc,i - 1 + (k + 2*l1)*ido);
01945 dr4 = wa3[i - 2]*ref(cc,i - 1 + (k + 3*l1)*ido) + wa3[i - 1]*ref(cc,i + (k + 3*l1)*ido);
01946 di4 = wa3[i - 2]*ref(cc,i + (k + 3*l1)*ido) - wa3[i - 1]*ref(cc,i - 1 + (k + 3*l1)*ido);
01947 dr5 = wa4[i - 2]*ref(cc,i - 1 + (k + 4*l1)*ido) + wa4[i - 1]*ref(cc,i + (k + 4*l1)*ido);
01948 di5 = wa4[i - 2]*ref(cc,i + (k + 4*l1)*ido) - wa4[i - 1]*ref(cc,i - 1 + (k + 4*l1)*ido);
01949 cr2 = dr2 + dr5;
01950 ci5 = dr5 - dr2;
01951 cr5 = di2 - di5;
01952 ci2 = di2 + di5;
01953 cr3 = dr3 + dr4;
01954 ci4 = dr4 - dr3;
01955 cr4 = di3 - di4;
01956 ci3 = di3 + di4;
01957 ch[i - 1 + 5*k*ido] = ref(cc,i - 1 + k*ido) + cr2 + cr3;
01958 ch[i + 5*k*ido] = ref(cc,i + k*ido) + ci2 + ci3;
01959 tr2 = ref(cc,i - 1 + k*ido) + tr11*cr2 + tr12*cr3;
01960 ti2 = ref(cc,i + k*ido) + tr11*ci2 + tr12*ci3;
01961 tr3 = ref(cc,i - 1 + k*ido) + tr12*cr2 + tr11*cr3;
01962 ti3 = ref(cc,i + k*ido) + tr12*ci2 + tr11*ci3;
01963 tr5 = ti11*cr5 + ti12*cr4;
01964 ti5 = ti11*ci5 + ti12*ci4;
01965 tr4 = ti12*cr5 - ti11*cr4;
01966 ti4 = ti12*ci5 - ti11*ci4;
01967 ch[i - 1 + (5*k + 2)*ido] = tr2 + tr5;
01968 ch[ic - 1 + (5*k + 1)*ido] = tr2 - tr5;
01969 ch[i + (5*k + 2)*ido] = ti2 + ti5;
01970 ch[ic + (5*k + 1)*ido] = ti5 - ti2;
01971 ch[i - 1 + (5*k + 4)*ido] = tr3 + tr4;
01972 ch[ic - 1 + (5*k + 3)*ido] = tr3 - tr4;
01973 ch[i + (5*k + 4)*ido] = ti3 + ti4;
01974 ch[ic + (5*k + 3)*ido] = ti4 - ti3;
01975 }
01976 }
01977 }
01978
01979
01980 static void radb5(int ido, int l1, const Treal cc[], Treal ch[],
01981 const Treal wa1[], const Treal wa2[], const Treal wa3[], const Treal wa4[])
01982 {
01983 static const Treal tr11 = 0.309016994374947;
01984 static const Treal ti11 = 0.951056516295154;
01985 static const Treal tr12 = -0.809016994374947;
01986 static const Treal ti12 = 0.587785252292473;
01987 int i, k, ic;
01988 Treal ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4, ti2, ti3,
01989 ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;
01990 for (k = 0; k < l1; k++) {
01991 ti5 = 2*ref(cc,(5*k + 2)*ido);
01992 ti4 = 2*ref(cc,(5*k + 4)*ido);
01993 tr2 = 2*ref(cc,ido-1 + (5*k + 1)*ido);
01994 tr3 = 2*ref(cc,ido-1 + (5*k + 3)*ido);
01995 ch[k*ido] = ref(cc,5*k*ido) + tr2 + tr3;
01996 cr2 = ref(cc,5*k*ido) + tr11*tr2 + tr12*tr3;
01997 cr3 = ref(cc,5*k*ido) + tr12*tr2 + tr11*tr3;
01998 ci5 = ti11*ti5 + ti12*ti4;
01999 ci4 = ti12*ti5 - ti11*ti4;
02000 ch[(k + l1)*ido] = cr2 - ci5;
02001 ch[(k + 2*l1)*ido] = cr3 - ci4;
02002 ch[(k + 3*l1)*ido] = cr3 + ci4;
02003 ch[(k + 4*l1)*ido] = cr2 + ci5;
02004 }
02005 if (ido == 1) return;
02006 for (k = 0; k < l1; ++k) {
02007 for (i = 2; i < ido; i += 2) {
02008 ic = ido - i;
02009 ti5 = ref(cc,i + (5*k + 2)*ido) + ref(cc,ic + (5*k + 1)*ido);
02010 ti2 = ref(cc,i + (5*k + 2)*ido) - ref(cc,ic + (5*k + 1)*ido);
02011 ti4 = ref(cc,i + (5*k + 4)*ido) + ref(cc,ic + (5*k + 3)*ido);
02012 ti3 = ref(cc,i + (5*k + 4)*ido) - ref(cc,ic + (5*k + 3)*ido);
02013 tr5 = ref(cc,i - 1 + (5*k + 2)*ido) - ref(cc,ic - 1 + (5*k + 1)*ido);
02014 tr2 = ref(cc,i - 1 + (5*k + 2)*ido) + ref(cc,ic - 1 + (5*k + 1)*ido);
02015 tr4 = ref(cc,i - 1 + (5*k + 4)*ido) - ref(cc,ic - 1 + (5*k + 3)*ido);
02016 tr3 = ref(cc,i - 1 + (5*k + 4)*ido) + ref(cc,ic - 1 + (5*k + 3)*ido);
02017 ch[i - 1 + k*ido] = ref(cc,i - 1 + 5*k*ido) + tr2 + tr3;
02018 ch[i + k*ido] = ref(cc,i + 5*k*ido) + ti2 + ti3;
02019 cr2 = ref(cc,i - 1 + 5*k*ido) + tr11*tr2 + tr12*tr3;
02020
02021 ci2 = ref(cc,i + 5*k*ido) + tr11*ti2 + tr12*ti3;
02022 cr3 = ref(cc,i - 1 + 5*k*ido) + tr12*tr2 + tr11*tr3;
02023
02024 ci3 = ref(cc,i + 5*k*ido) + tr12*ti2 + tr11*ti3;
02025 cr5 = ti11*tr5 + ti12*tr4;
02026 ci5 = ti11*ti5 + ti12*ti4;
02027 cr4 = ti12*tr5 - ti11*tr4;
02028 ci4 = ti12*ti5 - ti11*ti4;
02029 dr3 = cr3 - ci4;
02030 dr4 = cr3 + ci4;
02031 di3 = ci3 + cr4;
02032 di4 = ci3 - cr4;
02033 dr5 = cr2 + ci5;
02034 dr2 = cr2 - ci5;
02035 di5 = ci2 - cr5;
02036 di2 = ci2 + cr5;
02037 ch[i - 1 + (k + l1)*ido] = wa1[i - 2]*dr2 - wa1[i - 1]*di2;
02038 ch[i + (k + l1)*ido] = wa1[i - 2]*di2 + wa1[i - 1]*dr2;
02039 ch[i - 1 + (k + 2*l1)*ido] = wa2[i - 2]*dr3 - wa2[i - 1]*di3;
02040 ch[i + (k + 2*l1)*ido] = wa2[i - 2]*di3 + wa2[i - 1]*dr3;
02041 ch[i - 1 + (k + 3*l1)*ido] = wa3[i - 2]*dr4 - wa3[i - 1]*di4;
02042 ch[i + (k + 3*l1)*ido] = wa3[i - 2]*di4 + wa3[i - 1]*dr4;
02043 ch[i - 1 + (k + 4*l1)*ido] = wa4[i - 2]*dr5 - wa4[i - 1]*di5;
02044 ch[i + (k + 4*l1)*ido] = wa4[i - 2]*di5 + wa4[i - 1]*dr5;
02045 }
02046 }
02047 }
02048
02049
02050 static void radfg(int ido, int ip, int l1, int idl1,
02051 Treal cc[], Treal ch[], const Treal wa[])
02052 {
02053 static const Treal twopi = 6.28318530717959;
02054 int idij, ipph, i, j, k, l, j2, ic, jc, lc, ik, is, nbd;
02055 Treal dc2, ai1, ai2, ar1, ar2, ds2, dcp, arg, dsp, ar1h, ar2h;
02056 arg = twopi / ip;
02057 dcp = cos(arg);
02058 dsp = sin(arg);
02059 ipph = (ip + 1) / 2;
02060 nbd = (ido - 1) / 2;
02061 if (ido != 1) {
02062 for (ik=0; ik<idl1; ik++) ch[ik] = cc[ik];
02063 for (j=1; j<ip; j++)
02064 for (k=0; k<l1; k++)
02065 ch[(k + j*l1)*ido] = cc[(k + j*l1)*ido];
02066 if (nbd <= l1) {
02067 is = -ido;
02068 for (j=1; j<ip; j++) {
02069 is += ido;
02070 idij = is-1;
02071 for (i=2; i<ido; i+=2) {
02072 idij += 2;
02073 for (k=0; k<l1; k++) {
02074 ch[i - 1 + (k + j*l1)*ido] =
02075 wa[idij - 1]*cc[i - 1 + (k + j*l1)*ido] + wa[idij]*cc[i + (k + j*l1)*ido];
02076 ch[i + (k + j*l1)*ido] =
02077 wa[idij - 1]*cc[i + (k + j*l1)*ido] - wa[idij]*cc[i - 1 + (k + j*l1)*ido];
02078 }
02079 }
02080 }
02081 } else {
02082 is = -ido;
02083 for (j=1; j<ip; j++) {
02084 is += ido;
02085 for (k=0; k<l1; k++) {
02086 idij = is-1;
02087 for (i=2; i<ido; i+=2) {
02088 idij += 2;
02089 ch[i - 1 + (k + j*l1)*ido] =
02090 wa[idij - 1]*cc[i - 1 + (k + j*l1)*ido] + wa[idij]*cc[i + (k + j*l1)*ido];
02091 ch[i + (k + j*l1)*ido] =
02092 wa[idij - 1]*cc[i + (k + j*l1)*ido] - wa[idij]*cc[i - 1 + (k + j*l1)*ido];
02093 }
02094 }
02095 }
02096 }
02097 if (nbd >= l1) {
02098 for (j=1; j<ipph; j++) {
02099 jc = ip - j;
02100 for (k=0; k<l1; k++) {
02101 for (i=2; i<ido; i+=2) {
02102 cc[i - 1 + (k + j*l1)*ido] = ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido];
02103 cc[i - 1 + (k + jc*l1)*ido] = ch[i + (k + j*l1)*ido] - ch[i + (k + jc*l1)*ido];
02104 cc[i + (k + j*l1)*ido] = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido];
02105 cc[i + (k + jc*l1)*ido] = ch[i - 1 + (k + jc*l1)*ido] - ch[i - 1 + (k + j*l1)*ido];
02106 }
02107 }
02108 }
02109 } else {
02110 for (j=1; j<ipph; j++) {
02111 jc = ip - j;
02112 for (i=2; i<ido; i+=2) {
02113 for (k=0; k<l1; k++) {
02114 cc[i - 1 + (k + j*l1)*ido] =
02115 ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido];
02116 cc[i - 1 + (k + jc*l1)*ido] = ch[i + (k + j*l1)*ido] - ch[i + (k + jc*l1)*ido];
02117 cc[i + (k + j*l1)*ido] = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido];
02118 cc[i + (k + jc*l1)*ido] = ch[i - 1 + (k + jc*l1)*ido] - ch[i - 1 + (k + j*l1)*ido];
02119 }
02120 }
02121 }
02122 }
02123 } else {
02124 for (ik=0; ik<idl1; ik++) cc[ik] = ch[ik];
02125 }
02126 for (j=1; j<ipph; j++) {
02127 jc = ip - j;
02128 for (k=0; k<l1; k++) {
02129 cc[(k + j*l1)*ido] = ch[(k + j*l1)*ido] + ch[(k + jc*l1)*ido];
02130 cc[(k + jc*l1)*ido] = ch[(k + jc*l1)*ido] - ch[(k + j*l1)*ido];
02131 }
02132 }
02133
02134 ar1 = 1;
02135 ai1 = 0;
02136 for (l=1; l<ipph; l++) {
02137 lc = ip - l;
02138 ar1h = dcp*ar1 - dsp*ai1;
02139 ai1 = dcp*ai1 + dsp*ar1;
02140 ar1 = ar1h;
02141 for (ik=0; ik<idl1; ik++) {
02142 ch[ik + l*idl1] = cc[ik] + ar1*cc[ik + idl1];
02143 ch[ik + lc*idl1] = ai1*cc[ik + (ip-1)*idl1];
02144 }
02145 dc2 = ar1;
02146 ds2 = ai1;
02147 ar2 = ar1;
02148 ai2 = ai1;
02149 for (j=2; j<ipph; j++) {
02150 jc = ip - j;
02151 ar2h = dc2*ar2 - ds2*ai2;
02152 ai2 = dc2*ai2 + ds2*ar2;
02153 ar2 = ar2h;
02154 for (ik=0; ik<idl1; ik++) {
02155 ch[ik + l*idl1] += ar2*cc[ik + j*idl1];
02156 ch[ik + lc*idl1] += ai2*cc[ik + jc*idl1];
02157 }
02158 }
02159 }
02160 for (j=1; j<ipph; j++)
02161 for (ik=0; ik<idl1; ik++)
02162 ch[ik] += cc[ik + j*idl1];
02163
02164 if (ido >= l1) {
02165 for (k=0; k<l1; k++) {
02166 for (i=0; i<ido; i++) {
02167 ref(cc,i + k*ip*ido) = ch[i + k*ido];
02168 }
02169 }
02170 } else {
02171 for (i=0; i<ido; i++) {
02172 for (k=0; k<l1; k++) {
02173 ref(cc,i + k*ip*ido) = ch[i + k*ido];
02174 }
02175 }
02176 }
02177 for (j=1; j<ipph; j++) {
02178 jc = ip - j;
02179 j2 = 2*j;
02180 for (k=0; k<l1; k++) {
02181 ref(cc,ido-1 + (j2 - 1 + k*ip)*ido) =
02182 ch[(k + j*l1)*ido];
02183 ref(cc,(j2 + k*ip)*ido) =
02184 ch[(k + jc*l1)*ido];
02185 }
02186 }
02187 if (ido == 1) return;
02188 if (nbd >= l1) {
02189 for (j=1; j<ipph; j++) {
02190 jc = ip - j;
02191 j2 = 2*j;
02192 for (k=0; k<l1; k++) {
02193 for (i=2; i<ido; i+=2) {
02194 ic = ido - i;
02195 ref(cc,i - 1 + (j2 + k*ip)*ido) = ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido];
02196 ref(cc,ic - 1 + (j2 - 1 + k*ip)*ido) = ch[i - 1 + (k + j*l1)*ido] - ch[i - 1 + (k + jc*l1)*ido];
02197 ref(cc,i + (j2 + k*ip)*ido) = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido];
02198 ref(cc,ic + (j2 - 1 + k*ip)*ido) = ch[i + (k + jc*l1)*ido] - ch[i + (k + j*l1)*ido];
02199 }
02200 }
02201 }
02202 } else {
02203 for (j=1; j<ipph; j++) {
02204 jc = ip - j;
02205 j2 = 2*j;
02206 for (i=2; i<ido; i+=2) {
02207 ic = ido - i;
02208 for (k=0; k<l1; k++) {
02209 ref(cc,i - 1 + (j2 + k*ip)*ido) = ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido];
02210 ref(cc,ic - 1 + (j2 - 1 + k*ip)*ido) = ch[i - 1 + (k + j*l1)*ido] - ch[i - 1 + (k + jc*l1)*ido];
02211 ref(cc,i + (j2 + k*ip)*ido) = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido];
02212 ref(cc,ic + (j2 - 1 + k*ip)*ido) = ch[i + (k + jc*l1)*ido] - ch[i + (k + j*l1)*ido];
02213 }
02214 }
02215 }
02216 }
02217 }
02218
02219
02220 static void radbg(int ido, int ip, int l1, int idl1,
02221 Treal cc[], Treal ch[], const Treal wa[])
02222 {
02223 static const Treal twopi = 6.28318530717959;
02224 int idij, ipph, i, j, k, l, j2, ic, jc, lc, ik, is;
02225 Treal dc2, ai1, ai2, ar1, ar2, ds2;
02226 int nbd;
02227 Treal dcp, arg, dsp, ar1h, ar2h;
02228 arg = twopi / ip;
02229 dcp = cos(arg);
02230 dsp = sin(arg);
02231 nbd = (ido - 1) / 2;
02232 ipph = (ip + 1) / 2;
02233 if (ido >= l1) {
02234 for (k=0; k<l1; k++) {
02235 for (i=0; i<ido; i++) {
02236 ch[i + k*ido] = ref(cc,i + k*ip*ido);
02237 }
02238 }
02239 } else {
02240 for (i=0; i<ido; i++) {
02241 for (k=0; k<l1; k++) {
02242 ch[i + k*ido] = ref(cc,i + k*ip*ido);
02243 }
02244 }
02245 }
02246 for (j=1; j<ipph; j++) {
02247 jc = ip - j;
02248 j2 = 2*j;
02249 for (k=0; k<l1; k++) {
02250 ch[(k + j*l1)*ido] = ref(cc,ido-1 + (j2 - 1 + k*ip)*ido) + ref(cc,ido-1 + (j2 - 1 + k*ip)*
02251 ido);
02252 ch[(k + jc*l1)*ido] = ref(cc,(j2 + k*ip)*ido) + ref(cc,(j2 + k*ip)*ido);
02253 }
02254 }
02255
02256 if (ido != 1) {
02257 if (nbd >= l1) {
02258 for (j=1; j<ipph; j++) {
02259 jc = ip - j;
02260 for (k=0; k<l1; k++) {
02261 for (i=2; i<ido; i+=2) {
02262 ic = ido - i;
02263 ch[i - 1 + (k + j*l1)*ido] = ref(cc,i - 1 + (2*j + k*ip)*ido) + ref(cc,
02264 ic - 1 + (2*j - 1 + k*ip)*ido);
02265 ch[i - 1 + (k + jc*l1)*ido] = ref(cc,i - 1 + (2*j + k*ip)*ido) -
02266 ref(cc,ic - 1 + (2*j - 1 + k*ip)*ido);
02267 ch[i + (k + j*l1)*ido] = ref(cc,i + (2*j + k*ip)*ido) - ref(cc,ic
02268 + (2*j - 1 + k*ip)*ido);
02269 ch[i + (k + jc*l1)*ido] = ref(cc,i + (2*j + k*ip)*ido) + ref(cc,ic
02270 + (2*j - 1 + k*ip)*ido);
02271 }
02272 }
02273 }
02274 } else {
02275 for (j=1; j<ipph; j++) {
02276 jc = ip - j;
02277 for (i=2; i<ido; i+=2) {
02278 ic = ido - i;
02279 for (k=0; k<l1; k++) {
02280 ch[i - 1 + (k + j*l1)*ido] = ref(cc,i - 1 + (2*j + k*ip)*ido) + ref(cc,
02281 ic - 1 + (2*j - 1 + k*ip)*ido);
02282 ch[i - 1 + (k + jc*l1)*ido] = ref(cc,i - 1 + (2*j + k*ip)*ido) -
02283 ref(cc,ic - 1 + (2*j - 1 + k*ip)*ido);
02284 ch[i + (k + j*l1)*ido] = ref(cc,i + (2*j + k*ip)*ido) - ref(cc,ic
02285 + (2*j - 1 + k*ip)*ido);
02286 ch[i + (k + jc*l1)*ido] = ref(cc,i + (2*j + k*ip)*ido) + ref(cc,ic
02287 + (2*j - 1 + k*ip)*ido);
02288 }
02289 }
02290 }
02291 }
02292 }
02293
02294 ar1 = 1;
02295 ai1 = 0;
02296 for (l=1; l<ipph; l++) {
02297 lc = ip - l;
02298 ar1h = dcp*ar1 - dsp*ai1;
02299 ai1 = dcp*ai1 + dsp*ar1;
02300 ar1 = ar1h;
02301 for (ik=0; ik<idl1; ik++) {
02302 cc[ik + l*idl1] = ch[ik] + ar1*ch[ik + idl1];
02303 cc[ik + lc*idl1] = ai1*ch[ik + (ip-1)*idl1];
02304 }
02305 dc2 = ar1;
02306 ds2 = ai1;
02307 ar2 = ar1;
02308 ai2 = ai1;
02309 for (j=2; j<ipph; j++) {
02310 jc = ip - j;
02311 ar2h = dc2*ar2 - ds2*ai2;
02312 ai2 = dc2*ai2 + ds2*ar2;
02313 ar2 = ar2h;
02314 for (ik=0; ik<idl1; ik++) {
02315 cc[ik + l*idl1] += ar2*ch[ik + j*idl1];
02316 cc[ik + lc*idl1] += ai2*ch[ik + jc*idl1];
02317 }
02318 }
02319 }
02320 for (j=1; j<ipph; j++) {
02321 for (ik=0; ik<idl1; ik++) {
02322 ch[ik] += ch[ik + j*idl1];
02323 }
02324 }
02325 for (j=1; j<ipph; j++) {
02326 jc = ip - j;
02327 for (k=0; k<l1; k++) {
02328 ch[(k + j*l1)*ido] = cc[(k + j*l1)*ido] - cc[(k + jc*l1)*ido];
02329 ch[(k + jc*l1)*ido] = cc[(k + j*l1)*ido] + cc[(k + jc*l1)*ido];
02330 }
02331 }
02332
02333 if (ido == 1) return;
02334 if (nbd >= l1) {
02335 for (j=1; j<ipph; j++) {
02336 jc = ip - j;
02337 for (k=0; k<l1; k++) {
02338 for (i=2; i<ido; i+=2) {
02339 ch[i - 1 + (k + j*l1)*ido] = cc[i - 1 + (k + j*l1)*ido] - cc[i + (k + jc*l1)*ido];
02340 ch[i - 1 + (k + jc*l1)*ido] = cc[i - 1 + (k + j*l1)*ido] + cc[i + (k + jc*l1)*ido];
02341 ch[i + (k + j*l1)*ido] = cc[i + (k + j*l1)*ido] + cc[i - 1 + (k + jc*l1)*ido];
02342 ch[i + (k + jc*l1)*ido] = cc[i + (k + j*l1)*ido] - cc[i - 1 + (k + jc*l1)*ido];
02343 }
02344 }
02345 }
02346 } else {
02347 for (j=1; j<ipph; j++) {
02348 jc = ip - j;
02349 for (i=2; i<ido; i+=2) {
02350 for (k=0; k<l1; k++) {
02351 ch[i - 1 + (k + j*l1)*ido] = cc[i - 1 + (k + j*l1)*ido] - cc[i + (k + jc*l1)*ido];
02352 ch[i - 1 + (k + jc*l1)*ido] = cc[i - 1 + (k + j *l1)*ido] + cc[i + (k + jc*l1)*ido];
02353 ch[i + (k + j*l1)*ido] = cc[i + (k + j*l1)*ido] + cc[i - 1 + (k + jc*l1)*ido];
02354 ch[i + (k + jc*l1)*ido] = cc[i + (k + j*l1)*ido] - cc[i - 1 + (k + jc*l1)*ido];
02355 }
02356 }
02357 }
02358 }
02359 for (ik=0; ik<idl1; ik++) cc[ik] = ch[ik];
02360 for (j=1; j<ip; j++)
02361 for (k=0; k<l1; k++)
02362 cc[(k + j*l1)*ido] = ch[(k + j*l1)*ido];
02363 if (nbd <= l1) {
02364 is = -ido;
02365 for (j=1; j<ip; j++) {
02366 is += ido;
02367 idij = is-1;
02368 for (i=2; i<ido; i+=2) {
02369 idij += 2;
02370 for (k=0; k<l1; k++) {
02371 cc[i - 1 + (k + j*l1)*ido] = wa[idij - 1]*ch[i - 1 + (k + j*l1)*ido] - wa[idij]*
02372 ch[i + (k + j*l1)*ido];
02373 cc[i + (k + j*l1)*ido] = wa[idij - 1]*ch[i + (k + j*l1)*ido] + wa[idij]*ch[i - 1 + (k + j*l1)*ido];
02374 }
02375 }
02376 }
02377 } else {
02378 is = -ido;
02379 for (j=1; j<ip; j++) {
02380 is += ido;
02381 for (k=0; k<l1; k++) {
02382 idij = is - 1;
02383 for (i=2; i<ido; i+=2) {
02384 idij += 2;
02385 cc[i - 1 + (k + j*l1)*ido] = wa[idij-1]*ch[i - 1 + (k + j*l1)*ido] - wa[idij]*
02386 ch[i + (k + j*l1)*ido];
02387 cc[i + (k + j*l1)*ido] = wa[idij-1]*ch[i + (k + j*l1)*ido] + wa[idij]*ch[i - 1 + (k + j*l1)*ido];
02388 }
02389 }
02390 }
02391 }
02392 }
02393
02394
02395
02396
02397
02398 static void cfftf1(int n, Treal c[], Treal ch[], const Treal wa[], const int ifac[MAXFAC+2], int isign)
02399 {
02400 int idot, i;
02401 int k1, l1, l2;
02402 int na, nf, ip, iw, ix2, ix3, ix4, nac, ido, idl1;
02403 Treal *cinput, *coutput;
02404 nf = ifac[1];
02405 na = 0;
02406 l1 = 1;
02407 iw = 0;
02408 for (k1=2; k1<=nf+1; k1++) {
02409 ip = ifac[k1];
02410 l2 = ip*l1;
02411 ido = n / l2;
02412 idot = ido + ido;
02413 idl1 = idot*l1;
02414 if (na) {
02415 cinput = ch;
02416 coutput = c;
02417 } else {
02418 cinput = c;
02419 coutput = ch;
02420 }
02421 switch (ip) {
02422 case 4:
02423 ix2 = iw + idot;
02424 ix3 = ix2 + idot;
02425 passf4(idot, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], isign);
02426 na = !na;
02427 break;
02428 case 2:
02429 passf2(idot, l1, cinput, coutput, &wa[iw], isign);
02430 na = !na;
02431 break;
02432 case 3:
02433 ix2 = iw + idot;
02434 passf3(idot, l1, cinput, coutput, &wa[iw], &wa[ix2], isign);
02435 na = !na;
02436 break;
02437 case 5:
02438 ix2 = iw + idot;
02439 ix3 = ix2 + idot;
02440 ix4 = ix3 + idot;
02441 passf5(idot, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign);
02442 na = !na;
02443 break;
02444 default:
02445 passf(&nac, idot, ip, l1, idl1, cinput, coutput, &wa[iw], isign);
02446 if (nac != 0) na = !na;
02447 }
02448 l1 = l2;
02449 iw += (ip - 1)*idot;
02450 }
02451 if (na == 0) return;
02452 for (i=0; i<2*n; i++) c[i] = ch[i];
02453 }
02454
02455
02456 void cfftf(int n, Treal c[], Treal wsave[])
02457 {
02458 int iw1, iw2;
02459 if (n == 1) return;
02460 iw1 = 2*n;
02461 iw2 = iw1 + 2*n;
02462 cfftf1(n, c, wsave, wsave+iw1, (int*)(wsave+iw2), -1);
02463 }
02464
02465
02466 void cfftb(int n, Treal c[], Treal wsave[])
02467 {
02468 int iw1, iw2;
02469 if (n == 1) return;
02470 iw1 = 2*n;
02471 iw2 = iw1 + 2*n;
02472 cfftf1(n, c, wsave, wsave+iw1, (int*)(wsave+iw2), +1);
02473 }
02474
02475
02476 static void factorize(int n, int ifac[MAXFAC+2], const int ntryh[NSPECIAL])
02477
02478
02479
02480 {
02481 int ntry=3, i, j=0, ib, nf=0, nl=n, nq, nr;
02482 startloop:
02483 if (j < NSPECIAL)
02484 ntry = ntryh[j];
02485 else
02486 ntry+= 2;
02487 j++;
02488 do {
02489 nq = nl / ntry;
02490 nr = nl - ntry*nq;
02491 if (nr != 0) goto startloop;
02492 nf++;
02493 ifac[nf + 1] = ntry;
02494 nl = nq;
02495 if (ntry == 2 && nf != 1) {
02496 for (i=2; i<=nf; i++) {
02497 ib = nf - i + 2;
02498 ifac[ib + 1] = ifac[ib];
02499 }
02500 ifac[2] = 2;
02501 }
02502 } while (nl != 1);
02503 ifac[0] = n;
02504 ifac[1] = nf;
02505 }
02506
02507
02508 static void cffti1(int n, Treal wa[], int ifac[MAXFAC+2])
02509 {
02510 static const Treal twopi = 6.28318530717959;
02511 Treal arg, argh, argld, fi;
02512 int idot, i, j;
02513 int i1, k1, l1, l2;
02514 int ld, ii, nf, ip;
02515 int ido, ipm;
02516
02517 static const int ntryh[NSPECIAL] = {
02518 3,4,2,5 };
02519
02520 factorize(n,ifac,ntryh);
02521 nf = ifac[1];
02522 argh = twopi/(Treal)n;
02523 i = 1;
02524 l1 = 1;
02525 for (k1=1; k1<=nf; k1++) {
02526 ip = ifac[k1+1];
02527 ld = 0;
02528 l2 = l1*ip;
02529 ido = n / l2;
02530 idot = ido + ido + 2;
02531 ipm = ip - 1;
02532 for (j=1; j<=ipm; j++) {
02533 i1 = i;
02534 wa[i-1] = 1;
02535 wa[i] = 0;
02536 ld += l1;
02537 fi = 0;
02538 argld = ld*argh;
02539 for (ii=4; ii<=idot; ii+=2) {
02540 i+= 2;
02541 fi+= 1;
02542 arg = fi*argld;
02543 wa[i-1] = cos(arg);
02544 wa[i] = sin(arg);
02545 }
02546 if (ip > 5) {
02547 wa[i1-1] = wa[i-1];
02548 wa[i1] = wa[i];
02549 }
02550 }
02551 l1 = l2;
02552 }
02553 }
02554
02555
02556 void cffti(int n, Treal wsave[])
02557 {
02558 int iw1, iw2;
02559 if (n == 1) return;
02560 iw1 = 2*n;
02561 iw2 = iw1 + 2*n;
02562 cffti1(n, wsave+iw1, (int*)(wsave+iw2));
02563 }
02564
02565
02566
02567
02568
02569 static void rfftf1(int n, Treal c[], Treal ch[], const Treal wa[], const int ifac[MAXFAC+2])
02570 {
02571 int i;
02572 int k1, l1, l2, na, kh, nf, ip, iw, ix2, ix3, ix4, ido, idl1;
02573 Treal *cinput, *coutput;
02574 nf = ifac[1];
02575 na = 1;
02576 l2 = n;
02577 iw = n-1;
02578 for (k1 = 1; k1 <= nf; ++k1) {
02579 kh = nf - k1;
02580 ip = ifac[kh + 2];
02581 l1 = l2 / ip;
02582 ido = n / l2;
02583 idl1 = ido*l1;
02584 iw -= (ip - 1)*ido;
02585 na = !na;
02586 if (na) {
02587 cinput = ch;
02588 coutput = c;
02589 } else {
02590 cinput = c;
02591 coutput = ch;
02592 }
02593 switch (ip) {
02594 case 4:
02595 ix2 = iw + ido;
02596 ix3 = ix2 + ido;
02597 radf4(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3]);
02598 break;
02599 case 2:
02600 radf2(ido, l1, cinput, coutput, &wa[iw]);
02601 break;
02602 case 3:
02603 ix2 = iw + ido;
02604 radf3(ido, l1, cinput, coutput, &wa[iw], &wa[ix2]);
02605 break;
02606 case 5:
02607 ix2 = iw + ido;
02608 ix3 = ix2 + ido;
02609 ix4 = ix3 + ido;
02610 radf5(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4]);
02611 break;
02612 default:
02613 if (ido == 1)
02614 na = !na;
02615 if (na == 0) {
02616 radfg(ido, ip, l1, idl1, c, ch, &wa[iw]);
02617 na = 1;
02618 } else {
02619 radfg(ido, ip, l1, idl1, ch, c, &wa[iw]);
02620 na = 0;
02621 }
02622 }
02623 l2 = l1;
02624 }
02625 if (na == 1) return;
02626 for (i = 0; i < n; i++) c[i] = ch[i];
02627 }
02628
02629
02630 void rfftb1(int n, Treal c[], Treal ch[], const Treal wa[], const int ifac[MAXFAC+2])
02631 {
02632 int i;
02633 int k1, l1, l2, na, nf, ip, iw, ix2, ix3, ix4, ido, idl1;
02634 Treal *cinput, *coutput;
02635 nf = ifac[1];
02636 na = 0;
02637 l1 = 1;
02638 iw = 0;
02639 for (k1=1; k1<=nf; k1++) {
02640 ip = ifac[k1 + 1];
02641 l2 = ip*l1;
02642 ido = n / l2;
02643 idl1 = ido*l1;
02644 if (na) {
02645 cinput = ch;
02646 coutput = c;
02647 } else {
02648 cinput = c;
02649 coutput = ch;
02650 }
02651 switch (ip) {
02652 case 4:
02653 ix2 = iw + ido;
02654 ix3 = ix2 + ido;
02655 radb4(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3]);
02656 na = !na;
02657 break;
02658 case 2:
02659 radb2(ido, l1, cinput, coutput, &wa[iw]);
02660 na = !na;
02661 break;
02662 case 3:
02663 ix2 = iw + ido;
02664 radb3(ido, l1, cinput, coutput, &wa[iw], &wa[ix2]);
02665 na = !na;
02666 break;
02667 case 5:
02668 ix2 = iw + ido;
02669 ix3 = ix2 + ido;
02670 ix4 = ix3 + ido;
02671 radb5(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4]);
02672 na = !na;
02673 break;
02674 default:
02675 radbg(ido, ip, l1, idl1, cinput, coutput, &wa[iw]);
02676 if (ido == 1) na = !na;
02677 }
02678 l1 = l2;
02679 iw += (ip - 1)*ido;
02680 }
02681 if (na == 0) return;
02682 for (i=0; i<n; i++) c[i] = ch[i];
02683 }
02684
02685
02686 void EMAN::rfftf(int n, Treal r[], Treal wsave[])
02687 {
02688 if (n == 1) return;
02689 rfftf1(n, r, wsave, wsave+n, (int*)(wsave+2*n));
02690 }
02691
02692
02693 void EMAN::rfftb(int n, Treal r[], Treal wsave[])
02694 {
02695 if (n == 1) return;
02696 rfftb1(n, r, wsave, wsave+n, (int*)(wsave+2*n));
02697 }
02698
02699
02700 static void rffti1(int n, Treal wa[], int ifac[MAXFAC+2])
02701 {
02702 static const Treal twopi = 6.28318530717959;
02703 Treal arg, argh, argld, fi;
02704 int i, j;
02705 int k1, l1, l2;
02706 int ld, ii, nf, ip, is;
02707 int ido, ipm, nfm1;
02708 static const int ntryh[NSPECIAL] = {
02709 4,2,3,5 };
02710 factorize(n,ifac,ntryh);
02711 nf = ifac[1];
02712 argh = twopi / n;
02713 is = 0;
02714 nfm1 = nf - 1;
02715 l1 = 1;
02716 if (nfm1 == 0) return;
02717 for (k1 = 1; k1 <= nfm1; k1++) {
02718 ip = ifac[k1 + 1];
02719 ld = 0;
02720 l2 = l1*ip;
02721 ido = n / l2;
02722 ipm = ip - 1;
02723 for (j = 1; j <= ipm; ++j) {
02724 ld += l1;
02725 i = is;
02726 argld = (Treal) ld*argh;
02727 fi = 0;
02728 for (ii = 3; ii <= ido; ii += 2) {
02729 i += 2;
02730 fi += 1;
02731 arg = fi*argld;
02732 wa[i - 2] = cos(arg);
02733 wa[i - 1] = sin(arg);
02734 }
02735 is += ido;
02736 }
02737 l1 = l2;
02738 }
02739 }
02740
02741
02742 void EMAN::rffti(int n, Treal wsave[])
02743 {
02744 if (n == 1) return;
02745 rffti1(n, wsave+n, (int*)(wsave+2*n));
02746 }
02747
02748
02749
02750
02751
02752 #endif //NATIVE_FFT