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