00001
00002
00003
00004
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
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285 void cdft(int n, int isgn, double *a, int *ip, double *w)
00286 {
00287 void makewt(int nw, int *ip, double *w);
00288 void cftfsub(int n, double *a, int *ip, int nw, double *w);
00289 void cftbsub(int n, double *a, int *ip, int nw, double *w);
00290 int nw;
00291
00292 nw = ip[0];
00293 if (n > (nw << 2)) {
00294 nw = n >> 2;
00295 makewt(nw, ip, w);
00296 }
00297 if (isgn >= 0) {
00298 cftfsub(n, a, ip, nw, w);
00299 } else {
00300 cftbsub(n, a, ip, nw, w);
00301 }
00302 }
00303
00304
00305 void rdft(int n, int isgn, double *a, int *ip, double *w)
00306 {
00307 void makewt(int nw, int *ip, double *w);
00308 void makect(int nc, int *ip, double *c);
00309 void cftfsub(int n, double *a, int *ip, int nw, double *w);
00310 void cftbsub(int n, double *a, int *ip, int nw, double *w);
00311 void rftfsub(int n, double *a, int nc, double *c);
00312 void rftbsub(int n, double *a, int nc, double *c);
00313 int nw, nc;
00314 double xi;
00315
00316 nw = ip[0];
00317 if (n > (nw << 2)) {
00318 nw = n >> 2;
00319 makewt(nw, ip, w);
00320 }
00321 nc = ip[1];
00322 if (n > (nc << 2)) {
00323 nc = n >> 2;
00324 makect(nc, ip, w + nw);
00325 }
00326 if (isgn >= 0) {
00327 if (n > 4) {
00328 cftfsub(n, a, ip, nw, w);
00329 rftfsub(n, a, nc, w + nw);
00330 } else if (n == 4) {
00331 cftfsub(n, a, ip, nw, w);
00332 }
00333 xi = a[0] - a[1];
00334 a[0] += a[1];
00335 a[1] = xi;
00336 } else {
00337 a[1] = 0.5 * (a[0] - a[1]);
00338 a[0] -= a[1];
00339 if (n > 4) {
00340 rftbsub(n, a, nc, w + nw);
00341 cftbsub(n, a, ip, nw, w);
00342 } else if (n == 4) {
00343 cftbsub(n, a, ip, nw, w);
00344 }
00345 }
00346 }
00347
00348
00349 void ddct(int n, int isgn, double *a, int *ip, double *w)
00350 {
00351 void makewt(int nw, int *ip, double *w);
00352 void makect(int nc, int *ip, double *c);
00353 void cftfsub(int n, double *a, int *ip, int nw, double *w);
00354 void cftbsub(int n, double *a, int *ip, int nw, double *w);
00355 void rftfsub(int n, double *a, int nc, double *c);
00356 void rftbsub(int n, double *a, int nc, double *c);
00357 void dctsub(int n, double *a, int nc, double *c);
00358 int j, nw, nc;
00359 double xr;
00360
00361 nw = ip[0];
00362 if (n > (nw << 2)) {
00363 nw = n >> 2;
00364 makewt(nw, ip, w);
00365 }
00366 nc = ip[1];
00367 if (n > nc) {
00368 nc = n;
00369 makect(nc, ip, w + nw);
00370 }
00371 if (isgn < 0) {
00372 xr = a[n - 1];
00373 for (j = n - 2; j >= 2; j -= 2) {
00374 a[j + 1] = a[j] - a[j - 1];
00375 a[j] += a[j - 1];
00376 }
00377 a[1] = a[0] - xr;
00378 a[0] += xr;
00379 if (n > 4) {
00380 rftbsub(n, a, nc, w + nw);
00381 cftbsub(n, a, ip, nw, w);
00382 } else if (n == 4) {
00383 cftbsub(n, a, ip, nw, w);
00384 }
00385 }
00386 dctsub(n, a, nc, w + nw);
00387 if (isgn >= 0) {
00388 if (n > 4) {
00389 cftfsub(n, a, ip, nw, w);
00390 rftfsub(n, a, nc, w + nw);
00391 } else if (n == 4) {
00392 cftfsub(n, a, ip, nw, w);
00393 }
00394 xr = a[0] - a[1];
00395 a[0] += a[1];
00396 for (j = 2; j < n; j += 2) {
00397 a[j - 1] = a[j] - a[j + 1];
00398 a[j] += a[j + 1];
00399 }
00400 a[n - 1] = xr;
00401 }
00402 }
00403
00404
00405 void ddst(int n, int isgn, double *a, int *ip, double *w)
00406 {
00407 void makewt(int nw, int *ip, double *w);
00408 void makect(int nc, int *ip, double *c);
00409 void cftfsub(int n, double *a, int *ip, int nw, double *w);
00410 void cftbsub(int n, double *a, int *ip, int nw, double *w);
00411 void rftfsub(int n, double *a, int nc, double *c);
00412 void rftbsub(int n, double *a, int nc, double *c);
00413 void dstsub(int n, double *a, int nc, double *c);
00414 int j, nw, nc;
00415 double xr;
00416
00417 nw = ip[0];
00418 if (n > (nw << 2)) {
00419 nw = n >> 2;
00420 makewt(nw, ip, w);
00421 }
00422 nc = ip[1];
00423 if (n > nc) {
00424 nc = n;
00425 makect(nc, ip, w + nw);
00426 }
00427 if (isgn < 0) {
00428 xr = a[n - 1];
00429 for (j = n - 2; j >= 2; j -= 2) {
00430 a[j + 1] = -a[j] - a[j - 1];
00431 a[j] -= a[j - 1];
00432 }
00433 a[1] = a[0] + xr;
00434 a[0] -= xr;
00435 if (n > 4) {
00436 rftbsub(n, a, nc, w + nw);
00437 cftbsub(n, a, ip, nw, w);
00438 } else if (n == 4) {
00439 cftbsub(n, a, ip, nw, w);
00440 }
00441 }
00442 dstsub(n, a, nc, w + nw);
00443 if (isgn >= 0) {
00444 if (n > 4) {
00445 cftfsub(n, a, ip, nw, w);
00446 rftfsub(n, a, nc, w + nw);
00447 } else if (n == 4) {
00448 cftfsub(n, a, ip, nw, w);
00449 }
00450 xr = a[0] - a[1];
00451 a[0] += a[1];
00452 for (j = 2; j < n; j += 2) {
00453 a[j - 1] = -a[j] - a[j + 1];
00454 a[j] -= a[j + 1];
00455 }
00456 a[n - 1] = -xr;
00457 }
00458 }
00459
00460
00461 void dfct(int n, double *a, double *t, int *ip, double *w)
00462 {
00463 void makewt(int nw, int *ip, double *w);
00464 void makect(int nc, int *ip, double *c);
00465 void cftfsub(int n, double *a, int *ip, int nw, double *w);
00466 void rftfsub(int n, double *a, int nc, double *c);
00467 void dctsub(int n, double *a, int nc, double *c);
00468 int j, k, l, m, mh, nw, nc;
00469 double xr, xi, yr, yi;
00470
00471 nw = ip[0];
00472 if (n > (nw << 3)) {
00473 nw = n >> 3;
00474 makewt(nw, ip, w);
00475 }
00476 nc = ip[1];
00477 if (n > (nc << 1)) {
00478 nc = n >> 1;
00479 makect(nc, ip, w + nw);
00480 }
00481 m = n >> 1;
00482 yi = a[m];
00483 xi = a[0] + a[n];
00484 a[0] -= a[n];
00485 t[0] = xi - yi;
00486 t[m] = xi + yi;
00487 if (n > 2) {
00488 mh = m >> 1;
00489 for (j = 1; j < mh; j++) {
00490 k = m - j;
00491 xr = a[j] - a[n - j];
00492 xi = a[j] + a[n - j];
00493 yr = a[k] - a[n - k];
00494 yi = a[k] + a[n - k];
00495 a[j] = xr;
00496 a[k] = yr;
00497 t[j] = xi - yi;
00498 t[k] = xi + yi;
00499 }
00500 t[mh] = a[mh] + a[n - mh];
00501 a[mh] -= a[n - mh];
00502 dctsub(m, a, nc, w + nw);
00503 if (m > 4) {
00504 cftfsub(m, a, ip, nw, w);
00505 rftfsub(m, a, nc, w + nw);
00506 } else if (m == 4) {
00507 cftfsub(m, a, ip, nw, w);
00508 }
00509 a[n - 1] = a[0] - a[1];
00510 a[1] = a[0] + a[1];
00511 for (j = m - 2; j >= 2; j -= 2) {
00512 a[2 * j + 1] = a[j] + a[j + 1];
00513 a[2 * j - 1] = a[j] - a[j + 1];
00514 }
00515 l = 2;
00516 m = mh;
00517 while (m >= 2) {
00518 dctsub(m, t, nc, w + nw);
00519 if (m > 4) {
00520 cftfsub(m, t, ip, nw, w);
00521 rftfsub(m, t, nc, w + nw);
00522 } else if (m == 4) {
00523 cftfsub(m, t, ip, nw, w);
00524 }
00525 a[n - l] = t[0] - t[1];
00526 a[l] = t[0] + t[1];
00527 k = 0;
00528 for (j = 2; j < m; j += 2) {
00529 k += l << 2;
00530 a[k - l] = t[j] - t[j + 1];
00531 a[k + l] = t[j] + t[j + 1];
00532 }
00533 l <<= 1;
00534 mh = m >> 1;
00535 for (j = 0; j < mh; j++) {
00536 k = m - j;
00537 t[j] = t[m + k] - t[m + j];
00538 t[k] = t[m + k] + t[m + j];
00539 }
00540 t[mh] = t[m + mh];
00541 m = mh;
00542 }
00543 a[l] = t[0];
00544 a[n] = t[2] - t[1];
00545 a[0] = t[2] + t[1];
00546 } else {
00547 a[1] = a[0];
00548 a[2] = t[0];
00549 a[0] = t[1];
00550 }
00551 }
00552
00553
00554 void dfst(int n, double *a, double *t, int *ip, double *w)
00555 {
00556 void makewt(int nw, int *ip, double *w);
00557 void makect(int nc, int *ip, double *c);
00558 void cftfsub(int n, double *a, int *ip, int nw, double *w);
00559 void rftfsub(int n, double *a, int nc, double *c);
00560 void dstsub(int n, double *a, int nc, double *c);
00561 int j, k, l, m, mh, nw, nc;
00562 double xr, xi, yr, yi;
00563
00564 nw = ip[0];
00565 if (n > (nw << 3)) {
00566 nw = n >> 3;
00567 makewt(nw, ip, w);
00568 }
00569 nc = ip[1];
00570 if (n > (nc << 1)) {
00571 nc = n >> 1;
00572 makect(nc, ip, w + nw);
00573 }
00574 if (n > 2) {
00575 m = n >> 1;
00576 mh = m >> 1;
00577 for (j = 1; j < mh; j++) {
00578 k = m - j;
00579 xr = a[j] + a[n - j];
00580 xi = a[j] - a[n - j];
00581 yr = a[k] + a[n - k];
00582 yi = a[k] - a[n - k];
00583 a[j] = xr;
00584 a[k] = yr;
00585 t[j] = xi + yi;
00586 t[k] = xi - yi;
00587 }
00588 t[0] = a[mh] - a[n - mh];
00589 a[mh] += a[n - mh];
00590 a[0] = a[m];
00591 dstsub(m, a, nc, w + nw);
00592 if (m > 4) {
00593 cftfsub(m, a, ip, nw, w);
00594 rftfsub(m, a, nc, w + nw);
00595 } else if (m == 4) {
00596 cftfsub(m, a, ip, nw, w);
00597 }
00598 a[n - 1] = a[1] - a[0];
00599 a[1] = a[0] + a[1];
00600 for (j = m - 2; j >= 2; j -= 2) {
00601 a[2 * j + 1] = a[j] - a[j + 1];
00602 a[2 * j - 1] = -a[j] - a[j + 1];
00603 }
00604 l = 2;
00605 m = mh;
00606 while (m >= 2) {
00607 dstsub(m, t, nc, w + nw);
00608 if (m > 4) {
00609 cftfsub(m, t, ip, nw, w);
00610 rftfsub(m, t, nc, w + nw);
00611 } else if (m == 4) {
00612 cftfsub(m, t, ip, nw, w);
00613 }
00614 a[n - l] = t[1] - t[0];
00615 a[l] = t[0] + t[1];
00616 k = 0;
00617 for (j = 2; j < m; j += 2) {
00618 k += l << 2;
00619 a[k - l] = -t[j] - t[j + 1];
00620 a[k + l] = t[j] - t[j + 1];
00621 }
00622 l <<= 1;
00623 mh = m >> 1;
00624 for (j = 1; j < mh; j++) {
00625 k = m - j;
00626 t[j] = t[m + k] + t[m + j];
00627 t[k] = t[m + k] - t[m + j];
00628 }
00629 t[0] = t[m + mh];
00630 m = mh;
00631 }
00632 a[l] = t[0];
00633 }
00634 a[0] = 0;
00635 }
00636
00637
00638
00639
00640
00641 #include <math.h>
00642
00643 void makewt(int nw, int *ip, double *w)
00644 {
00645 void makeipt(int nw, int *ip);
00646 int j, nwh, nw0, nw1;
00647 double delta, wn4r, wk1r, wk1i, wk3r, wk3i;
00648
00649 ip[0] = nw;
00650 ip[1] = 1;
00651 if (nw > 2) {
00652 nwh = nw >> 1;
00653 delta = atan(1.0) / nwh;
00654 wn4r = cos(delta * nwh);
00655 w[0] = 1;
00656 w[1] = wn4r;
00657 if (nwh == 4) {
00658 w[2] = cos(delta * 2);
00659 w[3] = sin(delta * 2);
00660 } else if (nwh > 4) {
00661 makeipt(nw, ip);
00662 w[2] = 0.5 / cos(delta * 2);
00663 w[3] = 0.5 / cos(delta * 6);
00664 for (j = 4; j < nwh; j += 4) {
00665 w[j] = cos(delta * j);
00666 w[j + 1] = sin(delta * j);
00667 w[j + 2] = cos(3 * delta * j);
00668 w[j + 3] = -sin(3 * delta * j);
00669 }
00670 }
00671 nw0 = 0;
00672 while (nwh > 2) {
00673 nw1 = nw0 + nwh;
00674 nwh >>= 1;
00675 w[nw1] = 1;
00676 w[nw1 + 1] = wn4r;
00677 if (nwh == 4) {
00678 wk1r = w[nw0 + 4];
00679 wk1i = w[nw0 + 5];
00680 w[nw1 + 2] = wk1r;
00681 w[nw1 + 3] = wk1i;
00682 } else if (nwh > 4) {
00683 wk1r = w[nw0 + 4];
00684 wk3r = w[nw0 + 6];
00685 w[nw1 + 2] = 0.5 / wk1r;
00686 w[nw1 + 3] = 0.5 / wk3r;
00687 for (j = 4; j < nwh; j += 4) {
00688 wk1r = w[nw0 + 2 * j];
00689 wk1i = w[nw0 + 2 * j + 1];
00690 wk3r = w[nw0 + 2 * j + 2];
00691 wk3i = w[nw0 + 2 * j + 3];
00692 w[nw1 + j] = wk1r;
00693 w[nw1 + j + 1] = wk1i;
00694 w[nw1 + j + 2] = wk3r;
00695 w[nw1 + j + 3] = wk3i;
00696 }
00697 }
00698 nw0 = nw1;
00699 }
00700 }
00701 }
00702
00703
00704 void makeipt(int nw, int *ip)
00705 {
00706 int j, l, m, m2, p, q;
00707
00708 ip[2] = 0;
00709 ip[3] = 16;
00710 m = 2;
00711 for (l = nw; l > 32; l >>= 2) {
00712 m2 = m << 1;
00713 q = m2 << 3;
00714 for (j = m; j < m2; j++) {
00715 p = ip[j] << 2;
00716 ip[m + j] = p;
00717 ip[m2 + j] = p + q;
00718 }
00719 m = m2;
00720 }
00721 }
00722
00723
00724 void makect(int nc, int *ip, double *c)
00725 {
00726 int j, nch;
00727 double delta;
00728
00729 ip[1] = nc;
00730 if (nc > 1) {
00731 nch = nc >> 1;
00732 delta = atan(1.0) / nch;
00733 c[0] = cos(delta * nch);
00734 c[nch] = 0.5 * c[0];
00735 for (j = 1; j < nch; j++) {
00736 c[j] = 0.5 * cos(delta * j);
00737 c[nc - j] = 0.5 * sin(delta * j);
00738 }
00739 }
00740 }
00741
00742
00743
00744
00745
00746 #ifdef USE_CDFT_PTHREADS
00747 #define USE_CDFT_THREADS
00748 #ifndef CDFT_THREADS_BEGIN_N
00749 #define CDFT_THREADS_BEGIN_N 8192
00750 #endif
00751 #ifndef CDFT_4THREADS_BEGIN_N
00752 #define CDFT_4THREADS_BEGIN_N 65536
00753 #endif
00754 #include <pthread.h>
00755 #include <stdio.h>
00756 #include <stdlib.h>
00757 #define cdft_thread_t pthread_t
00758 #define cdft_thread_create(thp,func,argp) { \
00759 if (pthread_create(thp, NULL, func, (void *) argp) != 0) { \
00760 fprintf(stderr, "cdft thread error\n"); \
00761 exit(1); \
00762 } \
00763 }
00764 #define cdft_thread_wait(th) { \
00765 if (pthread_join(th, NULL) != 0) { \
00766 fprintf(stderr, "cdft thread error\n"); \
00767 exit(1); \
00768 } \
00769 }
00770 #endif
00771
00772
00773 #ifdef USE_CDFT_WINTHREADS
00774 #define USE_CDFT_THREADS
00775 #ifndef CDFT_THREADS_BEGIN_N
00776 #define CDFT_THREADS_BEGIN_N 32768
00777 #endif
00778 #ifndef CDFT_4THREADS_BEGIN_N
00779 #define CDFT_4THREADS_BEGIN_N 524288
00780 #endif
00781 #include <windows.h>
00782 #include <stdio.h>
00783 #include <stdlib.h>
00784 #define cdft_thread_t HANDLE
00785 #define cdft_thread_create(thp,func,argp) { \
00786 DWORD thid; \
00787 *(thp) = CreateThread(NULL, 0, (LPTHREAD_START_ROUTINE) func, (LPVOID) argp, 0, &thid); \
00788 if (*(thp) == 0) { \
00789 fprintf(stderr, "cdft thread error\n"); \
00790 exit(1); \
00791 } \
00792 }
00793 #define cdft_thread_wait(th) { \
00794 WaitForSingleObject(th, INFINITE); \
00795 CloseHandle(th); \
00796 }
00797 #endif
00798
00799
00800 void cftfsub(int n, double *a, int *ip, int nw, double *w)
00801 {
00802 void bitrv2(int n, int *ip, double *a);
00803 void bitrv216(double *a);
00804 void bitrv208(double *a);
00805 void cftf1st(int n, double *a, double *w);
00806 void cftrec4(int n, double *a, int nw, double *w);
00807 void cftleaf(int n, int isplt, double *a, int nw, double *w);
00808 void cftfx41(int n, double *a, int nw, double *w);
00809 void cftf161(double *a, double *w);
00810 void cftf081(double *a, double *w);
00811 void cftf040(double *a);
00812 void cftx020(double *a);
00813 #ifdef USE_CDFT_THREADS
00814 void cftrec4_th(int n, double *a, int nw, double *w);
00815 #endif
00816
00817 if (n > 8) {
00818 if (n > 32) {
00819 cftf1st(n, a, &w[nw - (n >> 2)]);
00820 #ifdef USE_CDFT_THREADS
00821 if (n > CDFT_THREADS_BEGIN_N) {
00822 cftrec4_th(n, a, nw, w);
00823 } else
00824 #endif
00825 if (n > 512) {
00826 cftrec4(n, a, nw, w);
00827 } else if (n > 128) {
00828 cftleaf(n, 1, a, nw, w);
00829 } else {
00830 cftfx41(n, a, nw, w);
00831 }
00832 bitrv2(n, ip, a);
00833 } else if (n == 32) {
00834 cftf161(a, &w[nw - 8]);
00835 bitrv216(a);
00836 } else {
00837 cftf081(a, w);
00838 bitrv208(a);
00839 }
00840 } else if (n == 8) {
00841 cftf040(a);
00842 } else if (n == 4) {
00843 cftx020(a);
00844 }
00845 }
00846
00847
00848 void cftbsub(int n, double *a, int *ip, int nw, double *w)
00849 {
00850 void bitrv2conj(int n, int *ip, double *a);
00851 void bitrv216neg(double *a);
00852 void bitrv208neg(double *a);
00853 void cftb1st(int n, double *a, double *w);
00854 void cftrec4(int n, double *a, int nw, double *w);
00855 void cftleaf(int n, int isplt, double *a, int nw, double *w);
00856 void cftfx41(int n, double *a, int nw, double *w);
00857 void cftf161(double *a, double *w);
00858 void cftf081(double *a, double *w);
00859 void cftb040(double *a);
00860 void cftx020(double *a);
00861 #ifdef USE_CDFT_THREADS
00862 void cftrec4_th(int n, double *a, int nw, double *w);
00863 #endif
00864
00865 if (n > 8) {
00866 if (n > 32) {
00867 cftb1st(n, a, &w[nw - (n >> 2)]);
00868 #ifdef USE_CDFT_THREADS
00869 if (n > CDFT_THREADS_BEGIN_N) {
00870 cftrec4_th(n, a, nw, w);
00871 } else
00872 #endif
00873 if (n > 512) {
00874 cftrec4(n, a, nw, w);
00875 } else if (n > 128) {
00876 cftleaf(n, 1, a, nw, w);
00877 } else {
00878 cftfx41(n, a, nw, w);
00879 }
00880 bitrv2conj(n, ip, a);
00881 } else if (n == 32) {
00882 cftf161(a, &w[nw - 8]);
00883 bitrv216neg(a);
00884 } else {
00885 cftf081(a, w);
00886 bitrv208neg(a);
00887 }
00888 } else if (n == 8) {
00889 cftb040(a);
00890 } else if (n == 4) {
00891 cftx020(a);
00892 }
00893 }
00894
00895
00896 void bitrv2(int n, int *ip, double *a)
00897 {
00898 int j, j1, k, k1, l, m, nh, nm;
00899 double xr, xi, yr, yi;
00900
00901 m = 1;
00902 for (l = n >> 2; l > 8; l >>= 2) {
00903 m <<= 1;
00904 }
00905 nh = n >> 1;
00906 nm = 4 * m;
00907 if (l == 8) {
00908 for (k = 0; k < m; k++) {
00909 for (j = 0; j < k; j++) {
00910 j1 = 4 * j + 2 * ip[m + k];
00911 k1 = 4 * k + 2 * ip[m + j];
00912 xr = a[j1];
00913 xi = a[j1 + 1];
00914 yr = a[k1];
00915 yi = a[k1 + 1];
00916 a[j1] = yr;
00917 a[j1 + 1] = yi;
00918 a[k1] = xr;
00919 a[k1 + 1] = xi;
00920 j1 += nm;
00921 k1 += 2 * nm;
00922 xr = a[j1];
00923 xi = a[j1 + 1];
00924 yr = a[k1];
00925 yi = a[k1 + 1];
00926 a[j1] = yr;
00927 a[j1 + 1] = yi;
00928 a[k1] = xr;
00929 a[k1 + 1] = xi;
00930 j1 += nm;
00931 k1 -= nm;
00932 xr = a[j1];
00933 xi = a[j1 + 1];
00934 yr = a[k1];
00935 yi = a[k1 + 1];
00936 a[j1] = yr;
00937 a[j1 + 1] = yi;
00938 a[k1] = xr;
00939 a[k1 + 1] = xi;
00940 j1 += nm;
00941 k1 += 2 * nm;
00942 xr = a[j1];
00943 xi = a[j1 + 1];
00944 yr = a[k1];
00945 yi = a[k1 + 1];
00946 a[j1] = yr;
00947 a[j1 + 1] = yi;
00948 a[k1] = xr;
00949 a[k1 + 1] = xi;
00950 j1 += nh;
00951 k1 += 2;
00952 xr = a[j1];
00953 xi = a[j1 + 1];
00954 yr = a[k1];
00955 yi = a[k1 + 1];
00956 a[j1] = yr;
00957 a[j1 + 1] = yi;
00958 a[k1] = xr;
00959 a[k1 + 1] = xi;
00960 j1 -= nm;
00961 k1 -= 2 * nm;
00962 xr = a[j1];
00963 xi = a[j1 + 1];
00964 yr = a[k1];
00965 yi = a[k1 + 1];
00966 a[j1] = yr;
00967 a[j1 + 1] = yi;
00968 a[k1] = xr;
00969 a[k1 + 1] = xi;
00970 j1 -= nm;
00971 k1 += nm;
00972 xr = a[j1];
00973 xi = a[j1 + 1];
00974 yr = a[k1];
00975 yi = a[k1 + 1];
00976 a[j1] = yr;
00977 a[j1 + 1] = yi;
00978 a[k1] = xr;
00979 a[k1 + 1] = xi;
00980 j1 -= nm;
00981 k1 -= 2 * nm;
00982 xr = a[j1];
00983 xi = a[j1 + 1];
00984 yr = a[k1];
00985 yi = a[k1 + 1];
00986 a[j1] = yr;
00987 a[j1 + 1] = yi;
00988 a[k1] = xr;
00989 a[k1 + 1] = xi;
00990 j1 += 2;
00991 k1 += nh;
00992 xr = a[j1];
00993 xi = a[j1 + 1];
00994 yr = a[k1];
00995 yi = a[k1 + 1];
00996 a[j1] = yr;
00997 a[j1 + 1] = yi;
00998 a[k1] = xr;
00999 a[k1 + 1] = xi;
01000 j1 += nm;
01001 k1 += 2 * nm;
01002 xr = a[j1];
01003 xi = a[j1 + 1];
01004 yr = a[k1];
01005 yi = a[k1 + 1];
01006 a[j1] = yr;
01007 a[j1 + 1] = yi;
01008 a[k1] = xr;
01009 a[k1 + 1] = xi;
01010 j1 += nm;
01011 k1 -= nm;
01012 xr = a[j1];
01013 xi = a[j1 + 1];
01014 yr = a[k1];
01015 yi = a[k1 + 1];
01016 a[j1] = yr;
01017 a[j1 + 1] = yi;
01018 a[k1] = xr;
01019 a[k1 + 1] = xi;
01020 j1 += nm;
01021 k1 += 2 * nm;
01022 xr = a[j1];
01023 xi = a[j1 + 1];
01024 yr = a[k1];
01025 yi = a[k1 + 1];
01026 a[j1] = yr;
01027 a[j1 + 1] = yi;
01028 a[k1] = xr;
01029 a[k1 + 1] = xi;
01030 j1 -= nh;
01031 k1 -= 2;
01032 xr = a[j1];
01033 xi = a[j1 + 1];
01034 yr = a[k1];
01035 yi = a[k1 + 1];
01036 a[j1] = yr;
01037 a[j1 + 1] = yi;
01038 a[k1] = xr;
01039 a[k1 + 1] = xi;
01040 j1 -= nm;
01041 k1 -= 2 * nm;
01042 xr = a[j1];
01043 xi = a[j1 + 1];
01044 yr = a[k1];
01045 yi = a[k1 + 1];
01046 a[j1] = yr;
01047 a[j1 + 1] = yi;
01048 a[k1] = xr;
01049 a[k1 + 1] = xi;
01050 j1 -= nm;
01051 k1 += nm;
01052 xr = a[j1];
01053 xi = a[j1 + 1];
01054 yr = a[k1];
01055 yi = a[k1 + 1];
01056 a[j1] = yr;
01057 a[j1 + 1] = yi;
01058 a[k1] = xr;
01059 a[k1 + 1] = xi;
01060 j1 -= nm;
01061 k1 -= 2 * nm;
01062 xr = a[j1];
01063 xi = a[j1 + 1];
01064 yr = a[k1];
01065 yi = a[k1 + 1];
01066 a[j1] = yr;
01067 a[j1 + 1] = yi;
01068 a[k1] = xr;
01069 a[k1 + 1] = xi;
01070 }
01071 k1 = 4 * k + 2 * ip[m + k];
01072 j1 = k1 + 2;
01073 k1 += nh;
01074 xr = a[j1];
01075 xi = a[j1 + 1];
01076 yr = a[k1];
01077 yi = a[k1 + 1];
01078 a[j1] = yr;
01079 a[j1 + 1] = yi;
01080 a[k1] = xr;
01081 a[k1 + 1] = xi;
01082 j1 += nm;
01083 k1 += 2 * nm;
01084 xr = a[j1];
01085 xi = a[j1 + 1];
01086 yr = a[k1];
01087 yi = a[k1 + 1];
01088 a[j1] = yr;
01089 a[j1 + 1] = yi;
01090 a[k1] = xr;
01091 a[k1 + 1] = xi;
01092 j1 += nm;
01093 k1 -= nm;
01094 xr = a[j1];
01095 xi = a[j1 + 1];
01096 yr = a[k1];
01097 yi = a[k1 + 1];
01098 a[j1] = yr;
01099 a[j1 + 1] = yi;
01100 a[k1] = xr;
01101 a[k1 + 1] = xi;
01102 j1 -= 2;
01103 k1 -= nh;
01104 xr = a[j1];
01105 xi = a[j1 + 1];
01106 yr = a[k1];
01107 yi = a[k1 + 1];
01108 a[j1] = yr;
01109 a[j1 + 1] = yi;
01110 a[k1] = xr;
01111 a[k1 + 1] = xi;
01112 j1 += nh + 2;
01113 k1 += nh + 2;
01114 xr = a[j1];
01115 xi = a[j1 + 1];
01116 yr = a[k1];
01117 yi = a[k1 + 1];
01118 a[j1] = yr;
01119 a[j1 + 1] = yi;
01120 a[k1] = xr;
01121 a[k1 + 1] = xi;
01122 j1 -= nh - nm;
01123 k1 += 2 * nm - 2;
01124 xr = a[j1];
01125 xi = a[j1 + 1];
01126 yr = a[k1];
01127 yi = a[k1 + 1];
01128 a[j1] = yr;
01129 a[j1 + 1] = yi;
01130 a[k1] = xr;
01131 a[k1 + 1] = xi;
01132 }
01133 } else {
01134 for (k = 0; k < m; k++) {
01135 for (j = 0; j < k; j++) {
01136 j1 = 4 * j + ip[m + k];
01137 k1 = 4 * k + ip[m + j];
01138 xr = a[j1];
01139 xi = a[j1 + 1];
01140 yr = a[k1];
01141 yi = a[k1 + 1];
01142 a[j1] = yr;
01143 a[j1 + 1] = yi;
01144 a[k1] = xr;
01145 a[k1 + 1] = xi;
01146 j1 += nm;
01147 k1 += nm;
01148 xr = a[j1];
01149 xi = a[j1 + 1];
01150 yr = a[k1];
01151 yi = a[k1 + 1];
01152 a[j1] = yr;
01153 a[j1 + 1] = yi;
01154 a[k1] = xr;
01155 a[k1 + 1] = xi;
01156 j1 += nh;
01157 k1 += 2;
01158 xr = a[j1];
01159 xi = a[j1 + 1];
01160 yr = a[k1];
01161 yi = a[k1 + 1];
01162 a[j1] = yr;
01163 a[j1 + 1] = yi;
01164 a[k1] = xr;
01165 a[k1 + 1] = xi;
01166 j1 -= nm;
01167 k1 -= nm;
01168 xr = a[j1];
01169 xi = a[j1 + 1];
01170 yr = a[k1];
01171 yi = a[k1 + 1];
01172 a[j1] = yr;
01173 a[j1 + 1] = yi;
01174 a[k1] = xr;
01175 a[k1 + 1] = xi;
01176 j1 += 2;
01177 k1 += nh;
01178 xr = a[j1];
01179 xi = a[j1 + 1];
01180 yr = a[k1];
01181 yi = a[k1 + 1];
01182 a[j1] = yr;
01183 a[j1 + 1] = yi;
01184 a[k1] = xr;
01185 a[k1 + 1] = xi;
01186 j1 += nm;
01187 k1 += nm;
01188 xr = a[j1];
01189 xi = a[j1 + 1];
01190 yr = a[k1];
01191 yi = a[k1 + 1];
01192 a[j1] = yr;
01193 a[j1 + 1] = yi;
01194 a[k1] = xr;
01195 a[k1 + 1] = xi;
01196 j1 -= nh;
01197 k1 -= 2;
01198 xr = a[j1];
01199 xi = a[j1 + 1];
01200 yr = a[k1];
01201 yi = a[k1 + 1];
01202 a[j1] = yr;
01203 a[j1 + 1] = yi;
01204 a[k1] = xr;
01205 a[k1 + 1] = xi;
01206 j1 -= nm;
01207 k1 -= nm;
01208 xr = a[j1];
01209 xi = a[j1 + 1];
01210 yr = a[k1];
01211 yi = a[k1 + 1];
01212 a[j1] = yr;
01213 a[j1 + 1] = yi;
01214 a[k1] = xr;
01215 a[k1 + 1] = xi;
01216 }
01217 k1 = 4 * k + ip[m + k];
01218 j1 = k1 + 2;
01219 k1 += nh;
01220 xr = a[j1];
01221 xi = a[j1 + 1];
01222 yr = a[k1];
01223 yi = a[k1 + 1];
01224 a[j1] = yr;
01225 a[j1 + 1] = yi;
01226 a[k1] = xr;
01227 a[k1 + 1] = xi;
01228 j1 += nm;
01229 k1 += nm;
01230 xr = a[j1];
01231 xi = a[j1 + 1];
01232 yr = a[k1];
01233 yi = a[k1 + 1];
01234 a[j1] = yr;
01235 a[j1 + 1] = yi;
01236 a[k1] = xr;
01237 a[k1 + 1] = xi;
01238 }
01239 }
01240 }
01241
01242
01243 void bitrv2conj(int n, int *ip, double *a)
01244 {
01245 int j, j1, k, k1, l, m, nh, nm;
01246 double xr, xi, yr, yi;
01247
01248 m = 1;
01249 for (l = n >> 2; l > 8; l >>= 2) {
01250 m <<= 1;
01251 }
01252 nh = n >> 1;
01253 nm = 4 * m;
01254 if (l == 8) {
01255 for (k = 0; k < m; k++) {
01256 for (j = 0; j < k; j++) {
01257 j1 = 4 * j + 2 * ip[m + k];
01258 k1 = 4 * k + 2 * ip[m + j];
01259 xr = a[j1];
01260 xi = -a[j1 + 1];
01261 yr = a[k1];
01262 yi = -a[k1 + 1];
01263 a[j1] = yr;
01264 a[j1 + 1] = yi;
01265 a[k1] = xr;
01266 a[k1 + 1] = xi;
01267 j1 += nm;
01268 k1 += 2 * nm;
01269 xr = a[j1];
01270 xi = -a[j1 + 1];
01271 yr = a[k1];
01272 yi = -a[k1 + 1];
01273 a[j1] = yr;
01274 a[j1 + 1] = yi;
01275 a[k1] = xr;
01276 a[k1 + 1] = xi;
01277 j1 += nm;
01278 k1 -= nm;
01279 xr = a[j1];
01280 xi = -a[j1 + 1];
01281 yr = a[k1];
01282 yi = -a[k1 + 1];
01283 a[j1] = yr;
01284 a[j1 + 1] = yi;
01285 a[k1] = xr;
01286 a[k1 + 1] = xi;
01287 j1 += nm;
01288 k1 += 2 * nm;
01289 xr = a[j1];
01290 xi = -a[j1 + 1];
01291 yr = a[k1];
01292 yi = -a[k1 + 1];
01293 a[j1] = yr;
01294 a[j1 + 1] = yi;
01295 a[k1] = xr;
01296 a[k1 + 1] = xi;
01297 j1 += nh;
01298 k1 += 2;
01299 xr = a[j1];
01300 xi = -a[j1 + 1];
01301 yr = a[k1];
01302 yi = -a[k1 + 1];
01303 a[j1] = yr;
01304 a[j1 + 1] = yi;
01305 a[k1] = xr;
01306 a[k1 + 1] = xi;
01307 j1 -= nm;
01308 k1 -= 2 * nm;
01309 xr = a[j1];
01310 xi = -a[j1 + 1];
01311 yr = a[k1];
01312 yi = -a[k1 + 1];
01313 a[j1] = yr;
01314 a[j1 + 1] = yi;
01315 a[k1] = xr;
01316 a[k1 + 1] = xi;
01317 j1 -= nm;
01318 k1 += nm;
01319 xr = a[j1];
01320 xi = -a[j1 + 1];
01321 yr = a[k1];
01322 yi = -a[k1 + 1];
01323 a[j1] = yr;
01324 a[j1 + 1] = yi;
01325 a[k1] = xr;
01326 a[k1 + 1] = xi;
01327 j1 -= nm;
01328 k1 -= 2 * nm;
01329 xr = a[j1];
01330 xi = -a[j1 + 1];
01331 yr = a[k1];
01332 yi = -a[k1 + 1];
01333 a[j1] = yr;
01334 a[j1 + 1] = yi;
01335 a[k1] = xr;
01336 a[k1 + 1] = xi;
01337 j1 += 2;
01338 k1 += nh;
01339 xr = a[j1];
01340 xi = -a[j1 + 1];
01341 yr = a[k1];
01342 yi = -a[k1 + 1];
01343 a[j1] = yr;
01344 a[j1 + 1] = yi;
01345 a[k1] = xr;
01346 a[k1 + 1] = xi;
01347 j1 += nm;
01348 k1 += 2 * nm;
01349 xr = a[j1];
01350 xi = -a[j1 + 1];
01351 yr = a[k1];
01352 yi = -a[k1 + 1];
01353 a[j1] = yr;
01354 a[j1 + 1] = yi;
01355 a[k1] = xr;
01356 a[k1 + 1] = xi;
01357 j1 += nm;
01358 k1 -= nm;
01359 xr = a[j1];
01360 xi = -a[j1 + 1];
01361 yr = a[k1];
01362 yi = -a[k1 + 1];
01363 a[j1] = yr;
01364 a[j1 + 1] = yi;
01365 a[k1] = xr;
01366 a[k1 + 1] = xi;
01367 j1 += nm;
01368 k1 += 2 * nm;
01369 xr = a[j1];
01370 xi = -a[j1 + 1];
01371 yr = a[k1];
01372 yi = -a[k1 + 1];
01373 a[j1] = yr;
01374 a[j1 + 1] = yi;
01375 a[k1] = xr;
01376 a[k1 + 1] = xi;
01377 j1 -= nh;
01378 k1 -= 2;
01379 xr = a[j1];
01380 xi = -a[j1 + 1];
01381 yr = a[k1];
01382 yi = -a[k1 + 1];
01383 a[j1] = yr;
01384 a[j1 + 1] = yi;
01385 a[k1] = xr;
01386 a[k1 + 1] = xi;
01387 j1 -= nm;
01388 k1 -= 2 * nm;
01389 xr = a[j1];
01390 xi = -a[j1 + 1];
01391 yr = a[k1];
01392 yi = -a[k1 + 1];
01393 a[j1] = yr;
01394 a[j1 + 1] = yi;
01395 a[k1] = xr;
01396 a[k1 + 1] = xi;
01397 j1 -= nm;
01398 k1 += nm;
01399 xr = a[j1];
01400 xi = -a[j1 + 1];
01401 yr = a[k1];
01402 yi = -a[k1 + 1];
01403 a[j1] = yr;
01404 a[j1 + 1] = yi;
01405 a[k1] = xr;
01406 a[k1 + 1] = xi;
01407 j1 -= nm;
01408 k1 -= 2 * nm;
01409 xr = a[j1];
01410 xi = -a[j1 + 1];
01411 yr = a[k1];
01412 yi = -a[k1 + 1];
01413 a[j1] = yr;
01414 a[j1 + 1] = yi;
01415 a[k1] = xr;
01416 a[k1 + 1] = xi;
01417 }
01418 k1 = 4 * k + 2 * ip[m + k];
01419 j1 = k1 + 2;
01420 k1 += nh;
01421 a[j1 - 1] = -a[j1 - 1];
01422 xr = a[j1];
01423 xi = -a[j1 + 1];
01424 yr = a[k1];
01425 yi = -a[k1 + 1];
01426 a[j1] = yr;
01427 a[j1 + 1] = yi;
01428 a[k1] = xr;
01429 a[k1 + 1] = xi;
01430 a[k1 + 3] = -a[k1 + 3];
01431 j1 += nm;
01432 k1 += 2 * nm;
01433 xr = a[j1];
01434 xi = -a[j1 + 1];
01435 yr = a[k1];
01436 yi = -a[k1 + 1];
01437 a[j1] = yr;
01438 a[j1 + 1] = yi;
01439 a[k1] = xr;
01440 a[k1 + 1] = xi;
01441 j1 += nm;
01442 k1 -= nm;
01443 xr = a[j1];
01444 xi = -a[j1 + 1];
01445 yr = a[k1];
01446 yi = -a[k1 + 1];
01447 a[j1] = yr;
01448 a[j1 + 1] = yi;
01449 a[k1] = xr;
01450 a[k1 + 1] = xi;
01451 j1 -= 2;
01452 k1 -= nh;
01453 xr = a[j1];
01454 xi = -a[j1 + 1];
01455 yr = a[k1];
01456 yi = -a[k1 + 1];
01457 a[j1] = yr;
01458 a[j1 + 1] = yi;
01459 a[k1] = xr;
01460 a[k1 + 1] = xi;
01461 j1 += nh + 2;
01462 k1 += nh + 2;
01463 xr = a[j1];
01464 xi = -a[j1 + 1];
01465 yr = a[k1];
01466 yi = -a[k1 + 1];
01467 a[j1] = yr;
01468 a[j1 + 1] = yi;
01469 a[k1] = xr;
01470 a[k1 + 1] = xi;
01471 j1 -= nh - nm;
01472 k1 += 2 * nm - 2;
01473 a[j1 - 1] = -a[j1 - 1];
01474 xr = a[j1];
01475 xi = -a[j1 + 1];
01476 yr = a[k1];
01477 yi = -a[k1 + 1];
01478 a[j1] = yr;
01479 a[j1 + 1] = yi;
01480 a[k1] = xr;
01481 a[k1 + 1] = xi;
01482 a[k1 + 3] = -a[k1 + 3];
01483 }
01484 } else {
01485 for (k = 0; k < m; k++) {
01486 for (j = 0; j < k; j++) {
01487 j1 = 4 * j + ip[m + k];
01488 k1 = 4 * k + ip[m + j];
01489 xr = a[j1];
01490 xi = -a[j1 + 1];
01491 yr = a[k1];
01492 yi = -a[k1 + 1];
01493 a[j1] = yr;
01494 a[j1 + 1] = yi;
01495 a[k1] = xr;
01496 a[k1 + 1] = xi;
01497 j1 += nm;
01498 k1 += nm;
01499 xr = a[j1];
01500 xi = -a[j1 + 1];
01501 yr = a[k1];
01502 yi = -a[k1 + 1];
01503 a[j1] = yr;
01504 a[j1 + 1] = yi;
01505 a[k1] = xr;
01506 a[k1 + 1] = xi;
01507 j1 += nh;
01508 k1 += 2;
01509 xr = a[j1];
01510 xi = -a[j1 + 1];
01511 yr = a[k1];
01512 yi = -a[k1 + 1];
01513 a[j1] = yr;
01514 a[j1 + 1] = yi;
01515 a[k1] = xr;
01516 a[k1 + 1] = xi;
01517 j1 -= nm;
01518 k1 -= nm;
01519 xr = a[j1];
01520 xi = -a[j1 + 1];
01521 yr = a[k1];
01522 yi = -a[k1 + 1];
01523 a[j1] = yr;
01524 a[j1 + 1] = yi;
01525 a[k1] = xr;
01526 a[k1 + 1] = xi;
01527 j1 += 2;
01528 k1 += nh;
01529 xr = a[j1];
01530 xi = -a[j1 + 1];
01531 yr = a[k1];
01532 yi = -a[k1 + 1];
01533 a[j1] = yr;
01534 a[j1 + 1] = yi;
01535 a[k1] = xr;
01536 a[k1 + 1] = xi;
01537 j1 += nm;
01538 k1 += nm;
01539 xr = a[j1];
01540 xi = -a[j1 + 1];
01541 yr = a[k1];
01542 yi = -a[k1 + 1];
01543 a[j1] = yr;
01544 a[j1 + 1] = yi;
01545 a[k1] = xr;
01546 a[k1 + 1] = xi;
01547 j1 -= nh;
01548 k1 -= 2;
01549 xr = a[j1];
01550 xi = -a[j1 + 1];
01551 yr = a[k1];
01552 yi = -a[k1 + 1];
01553 a[j1] = yr;
01554 a[j1 + 1] = yi;
01555 a[k1] = xr;
01556 a[k1 + 1] = xi;
01557 j1 -= nm;
01558 k1 -= nm;
01559 xr = a[j1];
01560 xi = -a[j1 + 1];
01561 yr = a[k1];
01562 yi = -a[k1 + 1];
01563 a[j1] = yr;
01564 a[j1 + 1] = yi;
01565 a[k1] = xr;
01566 a[k1 + 1] = xi;
01567 }
01568 k1 = 4 * k + ip[m + k];
01569 j1 = k1 + 2;
01570 k1 += nh;
01571 a[j1 - 1] = -a[j1 - 1];
01572 xr = a[j1];
01573 xi = -a[j1 + 1];
01574 yr = a[k1];
01575 yi = -a[k1 + 1];
01576 a[j1] = yr;
01577 a[j1 + 1] = yi;
01578 a[k1] = xr;
01579 a[k1 + 1] = xi;
01580 a[k1 + 3] = -a[k1 + 3];
01581 j1 += nm;
01582 k1 += nm;
01583 a[j1 - 1] = -a[j1 - 1];
01584 xr = a[j1];
01585 xi = -a[j1 + 1];
01586 yr = a[k1];
01587 yi = -a[k1 + 1];
01588 a[j1] = yr;
01589 a[j1 + 1] = yi;
01590 a[k1] = xr;
01591 a[k1 + 1] = xi;
01592 a[k1 + 3] = -a[k1 + 3];
01593 }
01594 }
01595 }
01596
01597
01598 void bitrv216(double *a)
01599 {
01600 double x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i,
01601 x5r, x5i, x7r, x7i, x8r, x8i, x10r, x10i,
01602 x11r, x11i, x12r, x12i, x13r, x13i, x14r, x14i;
01603
01604 x1r = a[2];
01605 x1i = a[3];
01606 x2r = a[4];
01607 x2i = a[5];
01608 x3r = a[6];
01609 x3i = a[7];
01610 x4r = a[8];
01611 x4i = a[9];
01612 x5r = a[10];
01613 x5i = a[11];
01614 x7r = a[14];
01615 x7i = a[15];
01616 x8r = a[16];
01617 x8i = a[17];
01618 x10r = a[20];
01619 x10i = a[21];
01620 x11r = a[22];
01621 x11i = a[23];
01622 x12r = a[24];
01623 x12i = a[25];
01624 x13r = a[26];
01625 x13i = a[27];
01626 x14r = a[28];
01627 x14i = a[29];
01628 a[2] = x8r;
01629 a[3] = x8i;
01630 a[4] = x4r;
01631 a[5] = x4i;
01632 a[6] = x12r;
01633 a[7] = x12i;
01634 a[8] = x2r;
01635 a[9] = x2i;
01636 a[10] = x10r;
01637 a[11] = x10i;
01638 a[14] = x14r;
01639 a[15] = x14i;
01640 a[16] = x1r;
01641 a[17] = x1i;
01642 a[20] = x5r;
01643 a[21] = x5i;
01644 a[22] = x13r;
01645 a[23] = x13i;
01646 a[24] = x3r;
01647 a[25] = x3i;
01648 a[26] = x11r;
01649 a[27] = x11i;
01650 a[28] = x7r;
01651 a[29] = x7i;
01652 }
01653
01654
01655 void bitrv216neg(double *a)
01656 {
01657 double x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i,
01658 x5r, x5i, x6r, x6i, x7r, x7i, x8r, x8i,
01659 x9r, x9i, x10r, x10i, x11r, x11i, x12r, x12i,
01660 x13r, x13i, x14r, x14i, x15r, x15i;
01661
01662 x1r = a[2];
01663 x1i = a[3];
01664 x2r = a[4];
01665 x2i = a[5];
01666 x3r = a[6];
01667 x3i = a[7];
01668 x4r = a[8];
01669 x4i = a[9];
01670 x5r = a[10];
01671 x5i = a[11];
01672 x6r = a[12];
01673 x6i = a[13];
01674 x7r = a[14];
01675 x7i = a[15];
01676 x8r = a[16];
01677 x8i = a[17];
01678 x9r = a[18];
01679 x9i = a[19];
01680 x10r = a[20];
01681 x10i = a[21];
01682 x11r = a[22];
01683 x11i = a[23];
01684 x12r = a[24];
01685 x12i = a[25];
01686 x13r = a[26];
01687 x13i = a[27];
01688 x14r = a[28];
01689 x14i = a[29];
01690 x15r = a[30];
01691 x15i = a[31];
01692 a[2] = x15r;
01693 a[3] = x15i;
01694 a[4] = x7r;
01695 a[5] = x7i;
01696 a[6] = x11r;
01697 a[7] = x11i;
01698 a[8] = x3r;
01699 a[9] = x3i;
01700 a[10] = x13r;
01701 a[11] = x13i;
01702 a[12] = x5r;
01703 a[13] = x5i;
01704 a[14] = x9r;
01705 a[15] = x9i;
01706 a[16] = x1r;
01707 a[17] = x1i;
01708 a[18] = x14r;
01709 a[19] = x14i;
01710 a[20] = x6r;
01711 a[21] = x6i;
01712 a[22] = x10r;
01713 a[23] = x10i;
01714 a[24] = x2r;
01715 a[25] = x2i;
01716 a[26] = x12r;
01717 a[27] = x12i;
01718 a[28] = x4r;
01719 a[29] = x4i;
01720 a[30] = x8r;
01721 a[31] = x8i;
01722 }
01723
01724
01725 void bitrv208(double *a)
01726 {
01727 double x1r, x1i, x3r, x3i, x4r, x4i, x6r, x6i;
01728
01729 x1r = a[2];
01730 x1i = a[3];
01731 x3r = a[6];
01732 x3i = a[7];
01733 x4r = a[8];
01734 x4i = a[9];
01735 x6r = a[12];
01736 x6i = a[13];
01737 a[2] = x4r;
01738 a[3] = x4i;
01739 a[6] = x6r;
01740 a[7] = x6i;
01741 a[8] = x1r;
01742 a[9] = x1i;
01743 a[12] = x3r;
01744 a[13] = x3i;
01745 }
01746
01747
01748 void bitrv208neg(double *a)
01749 {
01750 double x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i,
01751 x5r, x5i, x6r, x6i, x7r, x7i;
01752
01753 x1r = a[2];
01754 x1i = a[3];
01755 x2r = a[4];
01756 x2i = a[5];
01757 x3r = a[6];
01758 x3i = a[7];
01759 x4r = a[8];
01760 x4i = a[9];
01761 x5r = a[10];
01762 x5i = a[11];
01763 x6r = a[12];
01764 x6i = a[13];
01765 x7r = a[14];
01766 x7i = a[15];
01767 a[2] = x7r;
01768 a[3] = x7i;
01769 a[4] = x3r;
01770 a[5] = x3i;
01771 a[6] = x5r;
01772 a[7] = x5i;
01773 a[8] = x1r;
01774 a[9] = x1i;
01775 a[10] = x6r;
01776 a[11] = x6i;
01777 a[12] = x2r;
01778 a[13] = x2i;
01779 a[14] = x4r;
01780 a[15] = x4i;
01781 }
01782
01783
01784 void cftf1st(int n, double *a, double *w)
01785 {
01786 int j, j0, j1, j2, j3, k, m, mh;
01787 double wn4r, csc1, csc3, wk1r, wk1i, wk3r, wk3i,
01788 wd1r, wd1i, wd3r, wd3i;
01789 double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i,
01790 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i;
01791
01792 mh = n >> 3;
01793 m = 2 * mh;
01794 j1 = m;
01795 j2 = j1 + m;
01796 j3 = j2 + m;
01797 x0r = a[0] + a[j2];
01798 x0i = a[1] + a[j2 + 1];
01799 x1r = a[0] - a[j2];
01800 x1i = a[1] - a[j2 + 1];
01801 x2r = a[j1] + a[j3];
01802 x2i = a[j1 + 1] + a[j3 + 1];
01803 x3r = a[j1] - a[j3];
01804 x3i = a[j1 + 1] - a[j3 + 1];
01805 a[0] = x0r + x2r;
01806 a[1] = x0i + x2i;
01807 a[j1] = x0r - x2r;
01808 a[j1 + 1] = x0i - x2i;
01809 a[j2] = x1r - x3i;
01810 a[j2 + 1] = x1i + x3r;
01811 a[j3] = x1r + x3i;
01812 a[j3 + 1] = x1i - x3r;
01813 wn4r = w[1];
01814 csc1 = w[2];
01815 csc3 = w[3];
01816 wd1r = 1;
01817 wd1i = 0;
01818 wd3r = 1;
01819 wd3i = 0;
01820 k = 0;
01821 for (j = 2; j < mh - 2; j += 4) {
01822 k += 4;
01823 wk1r = csc1 * (wd1r + w[k]);
01824 wk1i = csc1 * (wd1i + w[k + 1]);
01825 wk3r = csc3 * (wd3r + w[k + 2]);
01826 wk3i = csc3 * (wd3i + w[k + 3]);
01827 wd1r = w[k];
01828 wd1i = w[k + 1];
01829 wd3r = w[k + 2];
01830 wd3i = w[k + 3];
01831 j1 = j + m;
01832 j2 = j1 + m;
01833 j3 = j2 + m;
01834 x0r = a[j] + a[j2];
01835 x0i = a[j + 1] + a[j2 + 1];
01836 x1r = a[j] - a[j2];
01837 x1i = a[j + 1] - a[j2 + 1];
01838 y0r = a[j + 2] + a[j2 + 2];
01839 y0i = a[j + 3] + a[j2 + 3];
01840 y1r = a[j + 2] - a[j2 + 2];
01841 y1i = a[j + 3] - a[j2 + 3];
01842 x2r = a[j1] + a[j3];
01843 x2i = a[j1 + 1] + a[j3 + 1];
01844 x3r = a[j1] - a[j3];
01845 x3i = a[j1 + 1] - a[j3 + 1];
01846 y2r = a[j1 + 2] + a[j3 + 2];
01847 y2i = a[j1 + 3] + a[j3 + 3];
01848 y3r = a[j1 + 2] - a[j3 + 2];
01849 y3i = a[j1 + 3] - a[j3 + 3];
01850 a[j] = x0r + x2r;
01851 a[j + 1] = x0i + x2i;
01852 a[j + 2] = y0r + y2r;
01853 a[j + 3] = y0i + y2i;
01854 a[j1] = x0r - x2r;
01855 a[j1 + 1] = x0i - x2i;
01856 a[j1 + 2] = y0r - y2r;
01857 a[j1 + 3] = y0i - y2i;
01858 x0r = x1r - x3i;
01859 x0i = x1i + x3r;
01860 a[j2] = wk1r * x0r - wk1i * x0i;
01861 a[j2 + 1] = wk1r * x0i + wk1i * x0r;
01862 x0r = y1r - y3i;
01863 x0i = y1i + y3r;
01864 a[j2 + 2] = wd1r * x0r - wd1i * x0i;
01865 a[j2 + 3] = wd1r * x0i + wd1i * x0r;
01866 x0r = x1r + x3i;
01867 x0i = x1i - x3r;
01868 a[j3] = wk3r * x0r + wk3i * x0i;
01869 a[j3 + 1] = wk3r * x0i - wk3i * x0r;
01870 x0r = y1r + y3i;
01871 x0i = y1i - y3r;
01872 a[j3 + 2] = wd3r * x0r + wd3i * x0i;
01873 a[j3 + 3] = wd3r * x0i - wd3i * x0r;
01874 j0 = m - j;
01875 j1 = j0 + m;
01876 j2 = j1 + m;
01877 j3 = j2 + m;
01878 x0r = a[j0] + a[j2];
01879 x0i = a[j0 + 1] + a[j2 + 1];
01880 x1r = a[j0] - a[j2];
01881 x1i = a[j0 + 1] - a[j2 + 1];
01882 y0r = a[j0 - 2] + a[j2 - 2];
01883 y0i = a[j0 - 1] + a[j2 - 1];
01884 y1r = a[j0 - 2] - a[j2 - 2];
01885 y1i = a[j0 - 1] - a[j2 - 1];
01886 x2r = a[j1] + a[j3];
01887 x2i = a[j1 + 1] + a[j3 + 1];
01888 x3r = a[j1] - a[j3];
01889 x3i = a[j1 + 1] - a[j3 + 1];
01890 y2r = a[j1 - 2] + a[j3 - 2];
01891 y2i = a[j1 - 1] + a[j3 - 1];
01892 y3r = a[j1 - 2] - a[j3 - 2];
01893 y3i = a[j1 - 1] - a[j3 - 1];
01894 a[j0] = x0r + x2r;
01895 a[j0 + 1] = x0i + x2i;
01896 a[j0 - 2] = y0r + y2r;
01897 a[j0 - 1] = y0i + y2i;
01898 a[j1] = x0r - x2r;
01899 a[j1 + 1] = x0i - x2i;
01900 a[j1 - 2] = y0r - y2r;
01901 a[j1 - 1] = y0i - y2i;
01902 x0r = x1r - x3i;
01903 x0i = x1i + x3r;
01904 a[j2] = wk1i * x0r - wk1r * x0i;
01905 a[j2 + 1] = wk1i * x0i + wk1r * x0r;
01906 x0r = y1r - y3i;
01907 x0i = y1i + y3r;
01908 a[j2 - 2] = wd1i * x0r - wd1r * x0i;
01909 a[j2 - 1] = wd1i * x0i + wd1r * x0r;
01910 x0r = x1r + x3i;
01911 x0i = x1i - x3r;
01912 a[j3] = wk3i * x0r + wk3r * x0i;
01913 a[j3 + 1] = wk3i * x0i - wk3r * x0r;
01914 x0r = y1r + y3i;
01915 x0i = y1i - y3r;
01916 a[j3 - 2] = wd3i * x0r + wd3r * x0i;
01917 a[j3 - 1] = wd3i * x0i - wd3r * x0r;
01918 }
01919 wk1r = csc1 * (wd1r + wn4r);
01920 wk1i = csc1 * (wd1i + wn4r);
01921 wk3r = csc3 * (wd3r - wn4r);
01922 wk3i = csc3 * (wd3i - wn4r);
01923 j0 = mh;
01924 j1 = j0 + m;
01925 j2 = j1 + m;
01926 j3 = j2 + m;
01927 x0r = a[j0 - 2] + a[j2 - 2];
01928 x0i = a[j0 - 1] + a[j2 - 1];
01929 x1r = a[j0 - 2] - a[j2 - 2];
01930 x1i = a[j0 - 1] - a[j2 - 1];
01931 x2r = a[j1 - 2] + a[j3 - 2];
01932 x2i = a[j1 - 1] + a[j3 - 1];
01933 x3r = a[j1 - 2] - a[j3 - 2];
01934 x3i = a[j1 - 1] - a[j3 - 1];
01935 a[j0 - 2] = x0r + x2r;
01936 a[j0 - 1] = x0i + x2i;
01937 a[j1 - 2] = x0r - x2r;
01938 a[j1 - 1] = x0i - x2i;
01939 x0r = x1r - x3i;
01940 x0i = x1i + x3r;
01941 a[j2 - 2] = wk1r * x0r - wk1i * x0i;
01942 a[j2 - 1] = wk1r * x0i + wk1i * x0r;
01943 x0r = x1r + x3i;
01944 x0i = x1i - x3r;
01945 a[j3 - 2] = wk3r * x0r + wk3i * x0i;
01946 a[j3 - 1] = wk3r * x0i - wk3i * x0r;
01947 x0r = a[j0] + a[j2];
01948 x0i = a[j0 + 1] + a[j2 + 1];
01949 x1r = a[j0] - a[j2];
01950 x1i = a[j0 + 1] - a[j2 + 1];
01951 x2r = a[j1] + a[j3];
01952 x2i = a[j1 + 1] + a[j3 + 1];
01953 x3r = a[j1] - a[j3];
01954 x3i = a[j1 + 1] - a[j3 + 1];
01955 a[j0] = x0r + x2r;
01956 a[j0 + 1] = x0i + x2i;
01957 a[j1] = x0r - x2r;
01958 a[j1 + 1] = x0i - x2i;
01959 x0r = x1r - x3i;
01960 x0i = x1i + x3r;
01961 a[j2] = wn4r * (x0r - x0i);
01962 a[j2 + 1] = wn4r * (x0i + x0r);
01963 x0r = x1r + x3i;
01964 x0i = x1i - x3r;
01965 a[j3] = -wn4r * (x0r + x0i);
01966 a[j3 + 1] = -wn4r * (x0i - x0r);
01967 x0r = a[j0 + 2] + a[j2 + 2];
01968 x0i = a[j0 + 3] + a[j2 + 3];
01969 x1r = a[j0 + 2] - a[j2 + 2];
01970 x1i = a[j0 + 3] - a[j2 + 3];
01971 x2r = a[j1 + 2] + a[j3 + 2];
01972 x2i = a[j1 + 3] + a[j3 + 3];
01973 x3r = a[j1 + 2] - a[j3 + 2];
01974 x3i = a[j1 + 3] - a[j3 + 3];
01975 a[j0 + 2] = x0r + x2r;
01976 a[j0 + 3] = x0i + x2i;
01977 a[j1 + 2] = x0r - x2r;
01978 a[j1 + 3] = x0i - x2i;
01979 x0r = x1r - x3i;
01980 x0i = x1i + x3r;
01981 a[j2 + 2] = wk1i * x0r - wk1r * x0i;
01982 a[j2 + 3] = wk1i * x0i + wk1r * x0r;
01983 x0r = x1r + x3i;
01984 x0i = x1i - x3r;
01985 a[j3 + 2] = wk3i * x0r + wk3r * x0i;
01986 a[j3 + 3] = wk3i * x0i - wk3r * x0r;
01987 }
01988
01989
01990 void cftb1st(int n, double *a, double *w)
01991 {
01992 int j, j0, j1, j2, j3, k, m, mh;
01993 double wn4r, csc1, csc3, wk1r, wk1i, wk3r, wk3i,
01994 wd1r, wd1i, wd3r, wd3i;
01995 double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i,
01996 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i;
01997
01998 mh = n >> 3;
01999 m = 2 * mh;
02000 j1 = m;
02001 j2 = j1 + m;
02002 j3 = j2 + m;
02003 x0r = a[0] + a[j2];
02004 x0i = -a[1] - a[j2 + 1];
02005 x1r = a[0] - a[j2];
02006 x1i = -a[1] + a[j2 + 1];
02007 x2r = a[j1] + a[j3];
02008 x2i = a[j1 + 1] + a[j3 + 1];
02009 x3r = a[j1] - a[j3];
02010 x3i = a[j1 + 1] - a[j3 + 1];
02011 a[0] = x0r + x2r;
02012 a[1] = x0i - x2i;
02013 a[j1] = x0r - x2r;
02014 a[j1 + 1] = x0i + x2i;
02015 a[j2] = x1r + x3i;
02016 a[j2 + 1] = x1i + x3r;
02017 a[j3] = x1r - x3i;
02018 a[j3 + 1] = x1i - x3r;
02019 wn4r = w[1];
02020 csc1 = w[2];
02021 csc3 = w[3];
02022 wd1r = 1;
02023 wd1i = 0;
02024 wd3r = 1;
02025 wd3i = 0;
02026 k = 0;
02027 for (j = 2; j < mh - 2; j += 4) {
02028 k += 4;
02029 wk1r = csc1 * (wd1r + w[k]);
02030 wk1i = csc1 * (wd1i + w[k + 1]);
02031 wk3r = csc3 * (wd3r + w[k + 2]);
02032 wk3i = csc3 * (wd3i + w[k + 3]);
02033 wd1r = w[k];
02034 wd1i = w[k + 1];
02035 wd3r = w[k + 2];
02036 wd3i = w[k + 3];
02037 j1 = j + m;
02038 j2 = j1 + m;
02039 j3 = j2 + m;
02040 x0r = a[j] + a[j2];
02041 x0i = -a[j + 1] - a[j2 + 1];
02042 x1r = a[j] - a[j2];
02043 x1i = -a[j + 1] + a[j2 + 1];
02044 y0r = a[j + 2] + a[j2 + 2];
02045 y0i = -a[j + 3] - a[j2 + 3];
02046 y1r = a[j + 2] - a[j2 + 2];
02047 y1i = -a[j + 3] + a[j2 + 3];
02048 x2r = a[j1] + a[j3];
02049 x2i = a[j1 + 1] + a[j3 + 1];
02050 x3r = a[j1] - a[j3];
02051 x3i = a[j1 + 1] - a[j3 + 1];
02052 y2r = a[j1 + 2] + a[j3 + 2];
02053 y2i = a[j1 + 3] + a[j3 + 3];
02054 y3r = a[j1 + 2] - a[j3 + 2];
02055 y3i = a[j1 + 3] - a[j3 + 3];
02056 a[j] = x0r + x2r;
02057 a[j + 1] = x0i - x2i;
02058 a[j + 2] = y0r + y2r;
02059 a[j + 3] = y0i - y2i;
02060 a[j1] = x0r - x2r;
02061 a[j1 + 1] = x0i + x2i;
02062 a[j1 + 2] = y0r - y2r;
02063 a[j1 + 3] = y0i + y2i;
02064 x0r = x1r + x3i;
02065 x0i = x1i + x3r;
02066 a[j2] = wk1r * x0r - wk1i * x0i;
02067 a[j2 + 1] = wk1r * x0i + wk1i * x0r;
02068 x0r = y1r + y3i;
02069 x0i = y1i + y3r;
02070 a[j2 + 2] = wd1r * x0r - wd1i * x0i;
02071 a[j2 + 3] = wd1r * x0i + wd1i * x0r;
02072 x0r = x1r - x3i;
02073 x0i = x1i - x3r;
02074 a[j3] = wk3r * x0r + wk3i * x0i;
02075 a[j3 + 1] = wk3r * x0i - wk3i * x0r;
02076 x0r = y1r - y3i;
02077 x0i = y1i - y3r;
02078 a[j3 + 2] = wd3r * x0r + wd3i * x0i;
02079 a[j3 + 3] = wd3r * x0i - wd3i * x0r;
02080 j0 = m - j;
02081 j1 = j0 + m;
02082 j2 = j1 + m;
02083 j3 = j2 + m;
02084 x0r = a[j0] + a[j2];
02085 x0i = -a[j0 + 1] - a[j2 + 1];
02086 x1r = a[j0] - a[j2];
02087 x1i = -a[j0 + 1] + a[j2 + 1];
02088 y0r = a[j0 - 2] + a[j2 - 2];
02089 y0i = -a[j0 - 1] - a[j2 - 1];
02090 y1r = a[j0 - 2] - a[j2 - 2];
02091 y1i = -a[j0 - 1] + a[j2 - 1];
02092 x2r = a[j1] + a[j3];
02093 x2i = a[j1 + 1] + a[j3 + 1];
02094 x3r = a[j1] - a[j3];
02095 x3i = a[j1 + 1] - a[j3 + 1];
02096 y2r = a[j1 - 2] + a[j3 - 2];
02097 y2i = a[j1 - 1] + a[j3 - 1];
02098 y3r = a[j1 - 2] - a[j3 - 2];
02099 y3i = a[j1 - 1] - a[j3 - 1];
02100 a[j0] = x0r + x2r;
02101 a[j0 + 1] = x0i - x2i;
02102 a[j0 - 2] = y0r + y2r;
02103 a[j0 - 1] = y0i - y2i;
02104 a[j1] = x0r - x2r;
02105 a[j1 + 1] = x0i + x2i;
02106 a[j1 - 2] = y0r - y2r;
02107 a[j1 - 1] = y0i + y2i;
02108 x0r = x1r + x3i;
02109 x0i = x1i + x3r;
02110 a[j2] = wk1i * x0r - wk1r * x0i;
02111 a[j2 + 1] = wk1i * x0i + wk1r * x0r;
02112 x0r = y1r + y3i;
02113 x0i = y1i + y3r;
02114 a[j2 - 2] = wd1i * x0r - wd1r * x0i;
02115 a[j2 - 1] = wd1i * x0i + wd1r * x0r;
02116 x0r = x1r - x3i;
02117 x0i = x1i - x3r;
02118 a[j3] = wk3i * x0r + wk3r * x0i;
02119 a[j3 + 1] = wk3i * x0i - wk3r * x0r;
02120 x0r = y1r - y3i;
02121 x0i = y1i - y3r;
02122 a[j3 - 2] = wd3i * x0r + wd3r * x0i;
02123 a[j3 - 1] = wd3i * x0i - wd3r * x0r;
02124 }
02125 wk1r = csc1 * (wd1r + wn4r);
02126 wk1i = csc1 * (wd1i + wn4r);
02127 wk3r = csc3 * (wd3r - wn4r);
02128 wk3i = csc3 * (wd3i - wn4r);
02129 j0 = mh;
02130 j1 = j0 + m;
02131 j2 = j1 + m;
02132 j3 = j2 + m;
02133 x0r = a[j0 - 2] + a[j2 - 2];
02134 x0i = -a[j0 - 1] - a[j2 - 1];
02135 x1r = a[j0 - 2] - a[j2 - 2];
02136 x1i = -a[j0 - 1] + a[j2 - 1];
02137 x2r = a[j1 - 2] + a[j3 - 2];
02138 x2i = a[j1 - 1] + a[j3 - 1];
02139 x3r = a[j1 - 2] - a[j3 - 2];
02140 x3i = a[j1 - 1] - a[j3 - 1];
02141 a[j0 - 2] = x0r + x2r;
02142 a[j0 - 1] = x0i - x2i;
02143 a[j1 - 2] = x0r - x2r;
02144 a[j1 - 1] = x0i + x2i;
02145 x0r = x1r + x3i;
02146 x0i = x1i + x3r;
02147 a[j2 - 2] = wk1r * x0r - wk1i * x0i;
02148 a[j2 - 1] = wk1r * x0i + wk1i * x0r;
02149 x0r = x1r - x3i;
02150 x0i = x1i - x3r;
02151 a[j3 - 2] = wk3r * x0r + wk3i * x0i;
02152 a[j3 - 1] = wk3r * x0i - wk3i * x0r;
02153 x0r = a[j0] + a[j2];
02154 x0i = -a[j0 + 1] - a[j2 + 1];
02155 x1r = a[j0] - a[j2];
02156 x1i = -a[j0 + 1] + a[j2 + 1];
02157 x2r = a[j1] + a[j3];
02158 x2i = a[j1 + 1] + a[j3 + 1];
02159 x3r = a[j1] - a[j3];
02160 x3i = a[j1 + 1] - a[j3 + 1];
02161 a[j0] = x0r + x2r;
02162 a[j0 + 1] = x0i - x2i;
02163 a[j1] = x0r - x2r;
02164 a[j1 + 1] = x0i + x2i;
02165 x0r = x1r + x3i;
02166 x0i = x1i + x3r;
02167 a[j2] = wn4r * (x0r - x0i);
02168 a[j2 + 1] = wn4r * (x0i + x0r);
02169 x0r = x1r - x3i;
02170 x0i = x1i - x3r;
02171 a[j3] = -wn4r * (x0r + x0i);
02172 a[j3 + 1] = -wn4r * (x0i - x0r);
02173 x0r = a[j0 + 2] + a[j2 + 2];
02174 x0i = -a[j0 + 3] - a[j2 + 3];
02175 x1r = a[j0 + 2] - a[j2 + 2];
02176 x1i = -a[j0 + 3] + a[j2 + 3];
02177 x2r = a[j1 + 2] + a[j3 + 2];
02178 x2i = a[j1 + 3] + a[j3 + 3];
02179 x3r = a[j1 + 2] - a[j3 + 2];
02180 x3i = a[j1 + 3] - a[j3 + 3];
02181 a[j0 + 2] = x0r + x2r;
02182 a[j0 + 3] = x0i - x2i;
02183 a[j1 + 2] = x0r - x2r;
02184 a[j1 + 3] = x0i + x2i;
02185 x0r = x1r + x3i;
02186 x0i = x1i + x3r;
02187 a[j2 + 2] = wk1i * x0r - wk1r * x0i;
02188 a[j2 + 3] = wk1i * x0i + wk1r * x0r;
02189 x0r = x1r - x3i;
02190 x0i = x1i - x3r;
02191 a[j3 + 2] = wk3i * x0r + wk3r * x0i;
02192 a[j3 + 3] = wk3i * x0i - wk3r * x0r;
02193 }
02194
02195
02196 #ifdef USE_CDFT_THREADS
02197 struct cdft_arg_st {
02198 int n0;
02199 int n;
02200 double *a;
02201 int nw;
02202 double *w;
02203 };
02204 typedef struct cdft_arg_st cdft_arg_t;
02205
02206
02207 void cftrec4_th(int n, double *a, int nw, double *w)
02208 {
02209 void *cftrec1_th(void *p);
02210 void *cftrec2_th(void *p);
02211 int i, idiv4, m, nthread;
02212 cdft_thread_t th[4];
02213 cdft_arg_t ag[4];
02214
02215 nthread = 2;
02216 idiv4 = 0;
02217 m = n >> 1;
02218 if (n > CDFT_4THREADS_BEGIN_N) {
02219 nthread = 4;
02220 idiv4 = 1;
02221 m >>= 1;
02222 }
02223 for (i = 0; i < nthread; i++) {
02224 ag[i].n0 = n;
02225 ag[i].n = m;
02226 ag[i].a = &a[i * m];
02227 ag[i].nw = nw;
02228 ag[i].w = w;
02229 if (i != idiv4) {
02230 cdft_thread_create(&th[i], cftrec1_th, &ag[i]);
02231 } else {
02232 cdft_thread_create(&th[i], cftrec2_th, &ag[i]);
02233 }
02234 }
02235 for (i = 0; i < nthread; i++) {
02236 cdft_thread_wait(th[i]);
02237 }
02238 }
02239
02240
02241 void *cftrec1_th(void *p)
02242 {
02243 int cfttree(int n, int j, int k, double *a, int nw, double *w);
02244 void cftleaf(int n, int isplt, double *a, int nw, double *w);
02245 void cftmdl1(int n, double *a, double *w);
02246 int isplt, j, k, m, n, n0, nw;
02247 double *a, *w;
02248
02249 n0 = ((cdft_arg_t *) p)->n0;
02250 n = ((cdft_arg_t *) p)->n;
02251 a = ((cdft_arg_t *) p)->a;
02252 nw = ((cdft_arg_t *) p)->nw;
02253 w = ((cdft_arg_t *) p)->w;
02254 m = n0;
02255 while (m > 512) {
02256 m >>= 2;
02257 cftmdl1(m, &a[n - m], &w[nw - (m >> 1)]);
02258 }
02259 cftleaf(m, 1, &a[n - m], nw, w);
02260 k = 0;
02261 for (j = n - m; j > 0; j -= m) {
02262 k++;
02263 isplt = cfttree(m, j, k, a, nw, w);
02264 cftleaf(m, isplt, &a[j - m], nw, w);
02265 }
02266 return (void *) 0;
02267 }
02268
02269
02270 void *cftrec2_th(void *p)
02271 {
02272 int cfttree(int n, int j, int k, double *a, int nw, double *w);
02273 void cftleaf(int n, int isplt, double *a, int nw, double *w);
02274 void cftmdl2(int n, double *a, double *w);
02275 int isplt, j, k, m, n, n0, nw;
02276 double *a, *w;
02277
02278 n0 = ((cdft_arg_t *) p)->n0;
02279 n = ((cdft_arg_t *) p)->n;
02280 a = ((cdft_arg_t *) p)->a;
02281 nw = ((cdft_arg_t *) p)->nw;
02282 w = ((cdft_arg_t *) p)->w;
02283 k = 1;
02284 m = n0;
02285 while (m > 512) {
02286 m >>= 2;
02287 k <<= 2;
02288 cftmdl2(m, &a[n - m], &w[nw - m]);
02289 }
02290 cftleaf(m, 0, &a[n - m], nw, w);
02291 k >>= 1;
02292 for (j = n - m; j > 0; j -= m) {
02293 k++;
02294 isplt = cfttree(m, j, k, a, nw, w);
02295 cftleaf(m, isplt, &a[j - m], nw, w);
02296 }
02297 return (void *) 0;
02298 }
02299 #endif
02300
02301
02302 void cftrec4(int n, double *a, int nw, double *w)
02303 {
02304 int cfttree(int n, int j, int k, double *a, int nw, double *w);
02305 void cftleaf(int n, int isplt, double *a, int nw, double *w);
02306 void cftmdl1(int n, double *a, double *w);
02307 int isplt, j, k, m;
02308
02309 m = n;
02310 while (m > 512) {
02311 m >>= 2;
02312 cftmdl1(m, &a[n - m], &w[nw - (m >> 1)]);
02313 }
02314 cftleaf(m, 1, &a[n - m], nw, w);
02315 k = 0;
02316 for (j = n - m; j > 0; j -= m) {
02317 k++;
02318 isplt = cfttree(m, j, k, a, nw, w);
02319 cftleaf(m, isplt, &a[j - m], nw, w);
02320 }
02321 }
02322
02323
02324 int cfttree(int n, int j, int k, double *a, int nw, double *w)
02325 {
02326 void cftmdl1(int n, double *a, double *w);
02327 void cftmdl2(int n, double *a, double *w);
02328 int i, isplt, m;
02329
02330 if ((k & 3) != 0) {
02331 isplt = k & 1;
02332 if (isplt != 0) {
02333 cftmdl1(n, &a[j - n], &w[nw - (n >> 1)]);
02334 } else {
02335 cftmdl2(n, &a[j - n], &w[nw - n]);
02336 }
02337 } else {
02338 m = n;
02339 for (i = k; (i & 3) == 0; i >>= 2) {
02340 m <<= 2;
02341 }
02342 isplt = i & 1;
02343 if (isplt != 0) {
02344 while (m > 128) {
02345 cftmdl1(m, &a[j - m], &w[nw - (m >> 1)]);
02346 m >>= 2;
02347 }
02348 } else {
02349 while (m > 128) {
02350 cftmdl2(m, &a[j - m], &w[nw - m]);
02351 m >>= 2;
02352 }
02353 }
02354 }
02355 return isplt;
02356 }
02357
02358
02359 void cftleaf(int n, int isplt, double *a, int nw, double *w)
02360 {
02361 void cftmdl1(int n, double *a, double *w);
02362 void cftmdl2(int n, double *a, double *w);
02363 void cftf161(double *a, double *w);
02364 void cftf162(double *a, double *w);
02365 void cftf081(double *a, double *w);
02366 void cftf082(double *a, double *w);
02367
02368 if (n == 512) {
02369 cftmdl1(128, a, &w[nw - 64]);
02370 cftf161(a, &w[nw - 8]);
02371 cftf162(&a[32], &w[nw - 32]);
02372 cftf161(&a[64], &w[nw - 8]);
02373 cftf161(&a[96], &w[nw - 8]);
02374 cftmdl2(128, &a[128], &w[nw - 128]);
02375 cftf161(&a[128], &w[nw - 8]);
02376 cftf162(&a[160], &w[nw - 32]);
02377 cftf161(&a[192], &w[nw - 8]);
02378 cftf162(&a[224], &w[nw - 32]);
02379 cftmdl1(128, &a[256], &w[nw - 64]);
02380 cftf161(&a[256], &w[nw - 8]);
02381 cftf162(&a[288], &w[nw - 32]);
02382 cftf161(&a[320], &w[nw - 8]);
02383 cftf161(&a[352], &w[nw - 8]);
02384 if (isplt != 0) {
02385 cftmdl1(128, &a[384], &w[nw - 64]);
02386 cftf161(&a[480], &w[nw - 8]);
02387 } else {
02388 cftmdl2(128, &a[384], &w[nw - 128]);
02389 cftf162(&a[480], &w[nw - 32]);
02390 }
02391 cftf161(&a[384], &w[nw - 8]);
02392 cftf162(&a[416], &w[nw - 32]);
02393 cftf161(&a[448], &w[nw - 8]);
02394 } else {
02395 cftmdl1(64, a, &w[nw - 32]);
02396 cftf081(a, &w[nw - 8]);
02397 cftf082(&a[16], &w[nw - 8]);
02398 cftf081(&a[32], &w[nw - 8]);
02399 cftf081(&a[48], &w[nw - 8]);
02400 cftmdl2(64, &a[64], &w[nw - 64]);
02401 cftf081(&a[64], &w[nw - 8]);
02402 cftf082(&a[80], &w[nw - 8]);
02403 cftf081(&a[96], &w[nw - 8]);
02404 cftf082(&a[112], &w[nw - 8]);
02405 cftmdl1(64, &a[128], &w[nw - 32]);
02406 cftf081(&a[128], &w[nw - 8]);
02407 cftf082(&a[144], &w[nw - 8]);
02408 cftf081(&a[160], &w[nw - 8]);
02409 cftf081(&a[176], &w[nw - 8]);
02410 if (isplt != 0) {
02411 cftmdl1(64, &a[192], &w[nw - 32]);
02412 cftf081(&a[240], &w[nw - 8]);
02413 } else {
02414 cftmdl2(64, &a[192], &w[nw - 64]);
02415 cftf082(&a[240], &w[nw - 8]);
02416 }
02417 cftf081(&a[192], &w[nw - 8]);
02418 cftf082(&a[208], &w[nw - 8]);
02419 cftf081(&a[224], &w[nw - 8]);
02420 }
02421 }
02422
02423
02424 void cftmdl1(int n, double *a, double *w)
02425 {
02426 int j, j0, j1, j2, j3, k, m, mh;
02427 double wn4r, wk1r, wk1i, wk3r, wk3i;
02428 double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
02429
02430 mh = n >> 3;
02431 m = 2 * mh;
02432 j1 = m;
02433 j2 = j1 + m;
02434 j3 = j2 + m;
02435 x0r = a[0] + a[j2];
02436 x0i = a[1] + a[j2 + 1];
02437 x1r = a[0] - a[j2];
02438 x1i = a[1] - a[j2 + 1];
02439 x2r = a[j1] + a[j3];
02440 x2i = a[j1 + 1] + a[j3 + 1];
02441 x3r = a[j1] - a[j3];
02442 x3i = a[j1 + 1] - a[j3 + 1];
02443 a[0] = x0r + x2r;
02444 a[1] = x0i + x2i;
02445 a[j1] = x0r - x2r;
02446 a[j1 + 1] = x0i - x2i;
02447 a[j2] = x1r - x3i;
02448 a[j2 + 1] = x1i + x3r;
02449 a[j3] = x1r + x3i;
02450 a[j3 + 1] = x1i - x3r;
02451 wn4r = w[1];
02452 k = 0;
02453 for (j = 2; j < mh; j += 2) {
02454 k += 4;
02455 wk1r = w[k];
02456 wk1i = w[k + 1];
02457 wk3r = w[k + 2];
02458 wk3i = w[k + 3];
02459 j1 = j + m;
02460 j2 = j1 + m;
02461 j3 = j2 + m;
02462 x0r = a[j] + a[j2];
02463 x0i = a[j + 1] + a[j2 + 1];
02464 x1r = a[j] - a[j2];
02465 x1i = a[j + 1] - a[j2 + 1];
02466 x2r = a[j1] + a[j3];
02467 x2i = a[j1 + 1] + a[j3 + 1];
02468 x3r = a[j1] - a[j3];
02469 x3i = a[j1 + 1] - a[j3 + 1];
02470 a[j] = x0r + x2r;
02471 a[j + 1] = x0i + x2i;
02472 a[j1] = x0r - x2r;
02473 a[j1 + 1] = x0i - x2i;
02474 x0r = x1r - x3i;
02475 x0i = x1i + x3r;
02476 a[j2] = wk1r * x0r - wk1i * x0i;
02477 a[j2 + 1] = wk1r * x0i + wk1i * x0r;
02478 x0r = x1r + x3i;
02479 x0i = x1i - x3r;
02480 a[j3] = wk3r * x0r + wk3i * x0i;
02481 a[j3 + 1] = wk3r * x0i - wk3i * x0r;
02482 j0 = m - j;
02483 j1 = j0 + m;
02484 j2 = j1 + m;
02485 j3 = j2 + m;
02486 x0r = a[j0] + a[j2];
02487 x0i = a[j0 + 1] + a[j2 + 1];
02488 x1r = a[j0] - a[j2];
02489 x1i = a[j0 + 1] - a[j2 + 1];
02490 x2r = a[j1] + a[j3];
02491 x2i = a[j1 + 1] + a[j3 + 1];
02492 x3r = a[j1] - a[j3];
02493 x3i = a[j1 + 1] - a[j3 + 1];
02494 a[j0] = x0r + x2r;
02495 a[j0 + 1] = x0i + x2i;
02496 a[j1] = x0r - x2r;
02497 a[j1 + 1] = x0i - x2i;
02498 x0r = x1r - x3i;
02499 x0i = x1i + x3r;
02500 a[j2] = wk1i * x0r - wk1r * x0i;
02501 a[j2 + 1] = wk1i * x0i + wk1r * x0r;
02502 x0r = x1r + x3i;
02503 x0i = x1i - x3r;
02504 a[j3] = wk3i * x0r + wk3r * x0i;
02505 a[j3 + 1] = wk3i * x0i - wk3r * x0r;
02506 }
02507 j0 = mh;
02508 j1 = j0 + m;
02509 j2 = j1 + m;
02510 j3 = j2 + m;
02511 x0r = a[j0] + a[j2];
02512 x0i = a[j0 + 1] + a[j2 + 1];
02513 x1r = a[j0] - a[j2];
02514 x1i = a[j0 + 1] - a[j2 + 1];
02515 x2r = a[j1] + a[j3];
02516 x2i = a[j1 + 1] + a[j3 + 1];
02517 x3r = a[j1] - a[j3];
02518 x3i = a[j1 + 1] - a[j3 + 1];
02519 a[j0] = x0r + x2r;
02520 a[j0 + 1] = x0i + x2i;
02521 a[j1] = x0r - x2r;
02522 a[j1 + 1] = x0i - x2i;
02523 x0r = x1r - x3i;
02524 x0i = x1i + x3r;
02525 a[j2] = wn4r * (x0r - x0i);
02526 a[j2 + 1] = wn4r * (x0i + x0r);
02527 x0r = x1r + x3i;
02528 x0i = x1i - x3r;
02529 a[j3] = -wn4r * (x0r + x0i);
02530 a[j3 + 1] = -wn4r * (x0i - x0r);
02531 }
02532
02533
02534 void cftmdl2(int n, double *a, double *w)
02535 {
02536 int j, j0, j1, j2, j3, k, kr, m, mh;
02537 double wn4r, wk1r, wk1i, wk3r, wk3i, wd1r, wd1i, wd3r, wd3i;
02538 double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, y0r, y0i, y2r, y2i;
02539
02540 mh = n >> 3;
02541 m = 2 * mh;
02542 wn4r = w[1];
02543 j1 = m;
02544 j2 = j1 + m;
02545 j3 = j2 + m;
02546 x0r = a[0] - a[j2 + 1];
02547 x0i = a[1] + a[j2];
02548 x1r = a[0] + a[j2 + 1];
02549 x1i = a[1] - a[j2];
02550 x2r = a[j1] - a[j3 + 1];
02551 x2i = a[j1 + 1] + a[j3];
02552 x3r = a[j1] + a[j3 + 1];
02553 x3i = a[j1 + 1] - a[j3];
02554 y0r = wn4r * (x2r - x2i);
02555 y0i = wn4r * (x2i + x2r);
02556 a[0] = x0r + y0r;
02557 a[1] = x0i + y0i;
02558 a[j1] = x0r - y0r;
02559 a[j1 + 1] = x0i - y0i;
02560 y0r = wn4r * (x3r - x3i);
02561 y0i = wn4r * (x3i + x3r);
02562 a[j2] = x1r - y0i;
02563 a[j2 + 1] = x1i + y0r;
02564 a[j3] = x1r + y0i;
02565 a[j3 + 1] = x1i - y0r;
02566 k = 0;
02567 kr = 2 * m;
02568 for (j = 2; j < mh; j += 2) {
02569 k += 4;
02570 wk1r = w[k];
02571 wk1i = w[k + 1];
02572 wk3r = w[k + 2];
02573 wk3i = w[k + 3];
02574 kr -= 4;
02575 wd1i = w[kr];
02576 wd1r = w[kr + 1];
02577 wd3i = w[kr + 2];
02578 wd3r = w[kr + 3];
02579 j1 = j + m;
02580 j2 = j1 + m;
02581 j3 = j2 + m;
02582 x0r = a[j] - a[j2 + 1];
02583 x0i = a[j + 1] + a[j2];
02584 x1r = a[j] + a[j2 + 1];
02585 x1i = a[j + 1] - a[j2];
02586 x2r = a[j1] - a[j3 + 1];
02587 x2i = a[j1 + 1] + a[j3];
02588 x3r = a[j1] + a[j3 + 1];
02589 x3i = a[j1 + 1] - a[j3];
02590 y0r = wk1r * x0r - wk1i * x0i;
02591 y0i = wk1r * x0i + wk1i * x0r;
02592 y2r = wd1r * x2r - wd1i * x2i;
02593 y2i = wd1r * x2i + wd1i * x2r;
02594 a[j] = y0r + y2r;
02595 a[j + 1] = y0i + y2i;
02596 a[j1] = y0r - y2r;
02597 a[j1 + 1] = y0i - y2i;
02598 y0r = wk3r * x1r + wk3i * x1i;
02599 y0i = wk3r * x1i - wk3i * x1r;
02600 y2r = wd3r * x3r + wd3i * x3i;
02601 y2i = wd3r * x3i - wd3i * x3r;
02602 a[j2] = y0r + y2r;
02603 a[j2 + 1] = y0i + y2i;
02604 a[j3] = y0r - y2r;
02605 a[j3 + 1] = y0i - y2i;
02606 j0 = m - j;
02607 j1 = j0 + m;
02608 j2 = j1 + m;
02609 j3 = j2 + m;
02610 x0r = a[j0] - a[j2 + 1];
02611 x0i = a[j0 + 1] + a[j2];
02612 x1r = a[j0] + a[j2 + 1];
02613 x1i = a[j0 + 1] - a[j2];
02614 x2r = a[j1] - a[j3 + 1];
02615 x2i = a[j1 + 1] + a[j3];
02616 x3r = a[j1] + a[j3 + 1];
02617 x3i = a[j1 + 1] - a[j3];
02618 y0r = wd1i * x0r - wd1r * x0i;
02619 y0i = wd1i * x0i + wd1r * x0r;
02620 y2r = wk1i * x2r - wk1r * x2i;
02621 y2i = wk1i * x2i + wk1r * x2r;
02622 a[j0] = y0r + y2r;
02623 a[j0 + 1] = y0i + y2i;
02624 a[j1] = y0r - y2r;
02625 a[j1 + 1] = y0i - y2i;
02626 y0r = wd3i * x1r + wd3r * x1i;
02627 y0i = wd3i * x1i - wd3r * x1r;
02628 y2r = wk3i * x3r + wk3r * x3i;
02629 y2i = wk3i * x3i - wk3r * x3r;
02630 a[j2] = y0r + y2r;
02631 a[j2 + 1] = y0i + y2i;
02632 a[j3] = y0r - y2r;
02633 a[j3 + 1] = y0i - y2i;
02634 }
02635 wk1r = w[m];
02636 wk1i = w[m + 1];
02637 j0 = mh;
02638 j1 = j0 + m;
02639 j2 = j1 + m;
02640 j3 = j2 + m;
02641 x0r = a[j0] - a[j2 + 1];
02642 x0i = a[j0 + 1] + a[j2];
02643 x1r = a[j0] + a[j2 + 1];
02644 x1i = a[j0 + 1] - a[j2];
02645 x2r = a[j1] - a[j3 + 1];
02646 x2i = a[j1 + 1] + a[j3];
02647 x3r = a[j1] + a[j3 + 1];
02648 x3i = a[j1 + 1] - a[j3];
02649 y0r = wk1r * x0r - wk1i * x0i;
02650 y0i = wk1r * x0i + wk1i * x0r;
02651 y2r = wk1i * x2r - wk1r * x2i;
02652 y2i = wk1i * x2i + wk1r * x2r;
02653 a[j0] = y0r + y2r;
02654 a[j0 + 1] = y0i + y2i;
02655 a[j1] = y0r - y2r;
02656 a[j1 + 1] = y0i - y2i;
02657 y0r = wk1i * x1r - wk1r * x1i;
02658 y0i = wk1i * x1i + wk1r * x1r;
02659 y2r = wk1r * x3r - wk1i * x3i;
02660 y2i = wk1r * x3i + wk1i * x3r;
02661 a[j2] = y0r - y2r;
02662 a[j2 + 1] = y0i - y2i;
02663 a[j3] = y0r + y2r;
02664 a[j3 + 1] = y0i + y2i;
02665 }
02666
02667
02668 void cftfx41(int n, double *a, int nw, double *w)
02669 {
02670 void cftf161(double *a, double *w);
02671 void cftf162(double *a, double *w);
02672 void cftf081(double *a, double *w);
02673 void cftf082(double *a, double *w);
02674
02675 if (n == 128) {
02676 cftf161(a, &w[nw - 8]);
02677 cftf162(&a[32], &w[nw - 32]);
02678 cftf161(&a[64], &w[nw - 8]);
02679 cftf161(&a[96], &w[nw - 8]);
02680 } else {
02681 cftf081(a, &w[nw - 8]);
02682 cftf082(&a[16], &w[nw - 8]);
02683 cftf081(&a[32], &w[nw - 8]);
02684 cftf081(&a[48], &w[nw - 8]);
02685 }
02686 }
02687
02688
02689 void cftf161(double *a, double *w)
02690 {
02691 double wn4r, wk1r, wk1i,
02692 x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i,
02693 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i,
02694 y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i,
02695 y8r, y8i, y9r, y9i, y10r, y10i, y11r, y11i,
02696 y12r, y12i, y13r, y13i, y14r, y14i, y15r, y15i;
02697
02698 wn4r = w[1];
02699 wk1r = w[2];
02700 wk1i = w[3];
02701 x0r = a[0] + a[16];
02702 x0i = a[1] + a[17];
02703 x1r = a[0] - a[16];
02704 x1i = a[1] - a[17];
02705 x2r = a[8] + a[24];
02706 x2i = a[9] + a[25];
02707 x3r = a[8] - a[24];
02708 x3i = a[9] - a[25];
02709 y0r = x0r + x2r;
02710 y0i = x0i + x2i;
02711 y4r = x0r - x2r;
02712 y4i = x0i - x2i;
02713 y8r = x1r - x3i;
02714 y8i = x1i + x3r;
02715 y12r = x1r + x3i;
02716 y12i = x1i - x3r;
02717 x0r = a[2] + a[18];
02718 x0i = a[3] + a[19];
02719 x1r = a[2] - a[18];
02720 x1i = a[3] - a[19];
02721 x2r = a[10] + a[26];
02722 x2i = a[11] + a[27];
02723 x3r = a[10] - a[26];
02724 x3i = a[11] - a[27];
02725 y1r = x0r + x2r;
02726 y1i = x0i + x2i;
02727 y5r = x0r - x2r;
02728 y5i = x0i - x2i;
02729 x0r = x1r - x3i;
02730 x0i = x1i + x3r;
02731 y9r = wk1r * x0r - wk1i * x0i;
02732 y9i = wk1r * x0i + wk1i * x0r;
02733 x0r = x1r + x3i;
02734 x0i = x1i - x3r;
02735 y13r = wk1i * x0r - wk1r * x0i;
02736 y13i = wk1i * x0i + wk1r * x0r;
02737 x0r = a[4] + a[20];
02738 x0i = a[5] + a[21];
02739 x1r = a[4] - a[20];
02740 x1i = a[5] - a[21];
02741 x2r = a[12] + a[28];
02742 x2i = a[13] + a[29];
02743 x3r = a[12] - a[28];
02744 x3i = a[13] - a[29];
02745 y2r = x0r + x2r;
02746 y2i = x0i + x2i;
02747 y6r = x0r - x2r;
02748 y6i = x0i - x2i;
02749 x0r = x1r - x3i;
02750 x0i = x1i + x3r;
02751 y10r = wn4r * (x0r - x0i);
02752 y10i = wn4r * (x0i + x0r);
02753 x0r = x1r + x3i;
02754 x0i = x1i - x3r;
02755 y14r = wn4r * (x0r + x0i);
02756 y14i = wn4r * (x0i - x0r);
02757 x0r = a[6] + a[22];
02758 x0i = a[7] + a[23];
02759 x1r = a[6] - a[22];
02760 x1i = a[7] - a[23];
02761 x2r = a[14] + a[30];
02762 x2i = a[15] + a[31];
02763 x3r = a[14] - a[30];
02764 x3i = a[15] - a[31];
02765 y3r = x0r + x2r;
02766 y3i = x0i + x2i;
02767 y7r = x0r - x2r;
02768 y7i = x0i - x2i;
02769 x0r = x1r - x3i;
02770 x0i = x1i + x3r;
02771 y11r = wk1i * x0r - wk1r * x0i;
02772 y11i = wk1i * x0i + wk1r * x0r;
02773 x0r = x1r + x3i;
02774 x0i = x1i - x3r;
02775 y15r = wk1r * x0r - wk1i * x0i;
02776 y15i = wk1r * x0i + wk1i * x0r;
02777 x0r = y12r - y14r;
02778 x0i = y12i - y14i;
02779 x1r = y12r + y14r;
02780 x1i = y12i + y14i;
02781 x2r = y13r - y15r;
02782 x2i = y13i - y15i;
02783 x3r = y13r + y15r;
02784 x3i = y13i + y15i;
02785 a[24] = x0r + x2r;
02786 a[25] = x0i + x2i;
02787 a[26] = x0r - x2r;
02788 a[27] = x0i - x2i;
02789 a[28] = x1r - x3i;
02790 a[29] = x1i + x3r;
02791 a[30] = x1r + x3i;
02792 a[31] = x1i - x3r;
02793 x0r = y8r + y10r;
02794 x0i = y8i + y10i;
02795 x1r = y8r - y10r;
02796 x1i = y8i - y10i;
02797 x2r = y9r + y11r;
02798 x2i = y9i + y11i;
02799 x3r = y9r - y11r;
02800 x3i = y9i - y11i;
02801 a[16] = x0r + x2r;
02802 a[17] = x0i + x2i;
02803 a[18] = x0r - x2r;
02804 a[19] = x0i - x2i;
02805 a[20] = x1r - x3i;
02806 a[21] = x1i + x3r;
02807 a[22] = x1r + x3i;
02808 a[23] = x1i - x3r;
02809 x0r = y5r - y7i;
02810 x0i = y5i + y7r;
02811 x2r = wn4r * (x0r - x0i);
02812 x2i = wn4r * (x0i + x0r);
02813 x0r = y5r + y7i;
02814 x0i = y5i - y7r;
02815 x3r = wn4r * (x0r - x0i);
02816 x3i = wn4r * (x0i + x0r);
02817 x0r = y4r - y6i;
02818 x0i = y4i + y6r;
02819 x1r = y4r + y6i;
02820 x1i = y4i - y6r;
02821 a[8] = x0r + x2r;
02822 a[9] = x0i + x2i;
02823 a[10] = x0r - x2r;
02824 a[11] = x0i - x2i;
02825 a[12] = x1r - x3i;
02826 a[13] = x1i + x3r;
02827 a[14] = x1r + x3i;
02828 a[15] = x1i - x3r;
02829 x0r = y0r + y2r;
02830 x0i = y0i + y2i;
02831 x1r = y0r - y2r;
02832 x1i = y0i - y2i;
02833 x2r = y1r + y3r;
02834 x2i = y1i + y3i;
02835 x3r = y1r - y3r;
02836 x3i = y1i - y3i;
02837 a[0] = x0r + x2r;
02838 a[1] = x0i + x2i;
02839 a[2] = x0r - x2r;
02840 a[3] = x0i - x2i;
02841 a[4] = x1r - x3i;
02842 a[5] = x1i + x3r;
02843 a[6] = x1r + x3i;
02844 a[7] = x1i - x3r;
02845 }
02846
02847
02848 void cftf162(double *a, double *w)
02849 {
02850 double wn4r, wk1r, wk1i, wk2r, wk2i, wk3r, wk3i,
02851 x0r, x0i, x1r, x1i, x2r, x2i,
02852 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i,
02853 y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i,
02854 y8r, y8i, y9r, y9i, y10r, y10i, y11r, y11i,
02855 y12r, y12i, y13r, y13i, y14r, y14i, y15r, y15i;
02856
02857 wn4r = w[1];
02858 wk1r = w[4];
02859 wk1i = w[5];
02860 wk3r = w[6];
02861 wk3i = -w[7];
02862 wk2r = w[8];
02863 wk2i = w[9];
02864 x1r = a[0] - a[17];
02865 x1i = a[1] + a[16];
02866 x0r = a[8] - a[25];
02867 x0i = a[9] + a[24];
02868 x2r = wn4r * (x0r - x0i);
02869 x2i = wn4r * (x0i + x0r);
02870 y0r = x1r + x2r;
02871 y0i = x1i + x2i;
02872 y4r = x1r - x2r;
02873 y4i = x1i - x2i;
02874 x1r = a[0] + a[17];
02875 x1i = a[1] - a[16];
02876 x0r = a[8] + a[25];
02877 x0i = a[9] - a[24];
02878 x2r = wn4r * (x0r - x0i);
02879 x2i = wn4r * (x0i + x0r);
02880 y8r = x1r - x2i;
02881 y8i = x1i + x2r;
02882 y12r = x1r + x2i;
02883 y12i = x1i - x2r;
02884 x0r = a[2] - a[19];
02885 x0i = a[3] + a[18];
02886 x1r = wk1r * x0r - wk1i * x0i;
02887 x1i = wk1r * x0i + wk1i * x0r;
02888 x0r = a[10] - a[27];
02889 x0i = a[11] + a[26];
02890 x2r = wk3i * x0r - wk3r * x0i;
02891 x2i = wk3i * x0i + wk3r * x0r;
02892 y1r = x1r + x2r;
02893 y1i = x1i + x2i;
02894 y5r = x1r - x2r;
02895 y5i = x1i - x2i;
02896 x0r = a[2] + a[19];
02897 x0i = a[3] - a[18];
02898 x1r = wk3r * x0r - wk3i * x0i;
02899 x1i = wk3r * x0i + wk3i * x0r;
02900 x0r = a[10] + a[27];
02901 x0i = a[11] - a[26];
02902 x2r = wk1r * x0r + wk1i * x0i;
02903 x2i = wk1r * x0i - wk1i * x0r;
02904 y9r = x1r - x2r;
02905 y9i = x1i - x2i;
02906 y13r = x1r + x2r;
02907 y13i = x1i + x2i;
02908 x0r = a[4] - a[21];
02909 x0i = a[5] + a[20];
02910 x1r = wk2r * x0r - wk2i * x0i;
02911 x1i = wk2r * x0i + wk2i * x0r;
02912 x0r = a[12] - a[29];
02913 x0i = a[13] + a[28];
02914 x2r = wk2i * x0r - wk2r * x0i;
02915 x2i = wk2i * x0i + wk2r * x0r;
02916 y2r = x1r + x2r;
02917 y2i = x1i + x2i;
02918 y6r = x1r - x2r;
02919 y6i = x1i - x2i;
02920 x0r = a[4] + a[21];
02921 x0i = a[5] - a[20];
02922 x1r = wk2i * x0r - wk2r * x0i;
02923 x1i = wk2i * x0i + wk2r * x0r;
02924 x0r = a[12] + a[29];
02925 x0i = a[13] - a[28];
02926 x2r = wk2r * x0r - wk2i * x0i;
02927 x2i = wk2r * x0i + wk2i * x0r;
02928 y10r = x1r - x2r;
02929 y10i = x1i - x2i;
02930 y14r = x1r + x2r;
02931 y14i = x1i + x2i;
02932 x0r = a[6] - a[23];
02933 x0i = a[7] + a[22];
02934 x1r = wk3r * x0r - wk3i * x0i;
02935 x1i = wk3r * x0i + wk3i * x0r;
02936 x0r = a[14] - a[31];
02937 x0i = a[15] + a[30];
02938 x2r = wk1i * x0r - wk1r * x0i;
02939 x2i = wk1i * x0i + wk1r * x0r;
02940 y3r = x1r + x2r;
02941 y3i = x1i + x2i;
02942 y7r = x1r - x2r;
02943 y7i = x1i - x2i;
02944 x0r = a[6] + a[23];
02945 x0i = a[7] - a[22];
02946 x1r = wk1i * x0r + wk1r * x0i;
02947 x1i = wk1i * x0i - wk1r * x0r;
02948 x0r = a[14] + a[31];
02949 x0i = a[15] - a[30];
02950 x2r = wk3i * x0r - wk3r * x0i;
02951 x2i = wk3i * x0i + wk3r * x0r;
02952 y11r = x1r + x2r;
02953 y11i = x1i + x2i;
02954 y15r = x1r - x2r;
02955 y15i = x1i - x2i;
02956 x1r = y0r + y2r;
02957 x1i = y0i + y2i;
02958 x2r = y1r + y3r;
02959 x2i = y1i + y3i;
02960 a[0] = x1r + x2r;
02961 a[1] = x1i + x2i;
02962 a[2] = x1r - x2r;
02963 a[3] = x1i - x2i;
02964 x1r = y0r - y2r;
02965 x1i = y0i - y2i;
02966 x2r = y1r - y3r;
02967 x2i = y1i - y3i;
02968 a[4] = x1r - x2i;
02969 a[5] = x1i + x2r;
02970 a[6] = x1r + x2i;
02971 a[7] = x1i - x2r;
02972 x1r = y4r - y6i;
02973 x1i = y4i + y6r;
02974 x0r = y5r - y7i;
02975 x0i = y5i + y7r;
02976 x2r = wn4r * (x0r - x0i);
02977 x2i = wn4r * (x0i + x0r);
02978 a[8] = x1r + x2r;
02979 a[9] = x1i + x2i;
02980 a[10] = x1r - x2r;
02981 a[11] = x1i - x2i;
02982 x1r = y4r + y6i;
02983 x1i = y4i - y6r;
02984 x0r = y5r + y7i;
02985 x0i = y5i - y7r;
02986 x2r = wn4r * (x0r - x0i);
02987 x2i = wn4r * (x0i + x0r);
02988 a[12] = x1r - x2i;
02989 a[13] = x1i + x2r;
02990 a[14] = x1r + x2i;
02991 a[15] = x1i - x2r;
02992 x1r = y8r + y10r;
02993 x1i = y8i + y10i;
02994 x2r = y9r - y11r;
02995 x2i = y9i - y11i;
02996 a[16] = x1r + x2r;
02997 a[17] = x1i + x2i;
02998 a[18] = x1r - x2r;
02999 a[19] = x1i - x2i;
03000 x1r = y8r - y10r;
03001 x1i = y8i - y10i;
03002 x2r = y9r + y11r;
03003 x2i = y9i + y11i;
03004 a[20] = x1r - x2i;
03005 a[21] = x1i + x2r;
03006 a[22] = x1r + x2i;
03007 a[23] = x1i - x2r;
03008 x1r = y12r - y14i;
03009 x1i = y12i + y14r;
03010 x0r = y13r + y15i;
03011 x0i = y13i - y15r;
03012 x2r = wn4r * (x0r - x0i);
03013 x2i = wn4r * (x0i + x0r);
03014 a[24] = x1r + x2r;
03015 a[25] = x1i + x2i;
03016 a[26] = x1r - x2r;
03017 a[27] = x1i - x2i;
03018 x1r = y12r + y14i;
03019 x1i = y12i - y14r;
03020 x0r = y13r - y15i;
03021 x0i = y13i + y15r;
03022 x2r = wn4r * (x0r - x0i);
03023 x2i = wn4r * (x0i + x0r);
03024 a[28] = x1r - x2i;
03025 a[29] = x1i + x2r;
03026 a[30] = x1r + x2i;
03027 a[31] = x1i - x2r;
03028 }
03029
03030
03031 void cftf081(double *a, double *w)
03032 {
03033 double wn4r, x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i,
03034 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i,
03035 y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i;
03036
03037 wn4r = w[1];
03038 x0r = a[0] + a[8];
03039 x0i = a[1] + a[9];
03040 x1r = a[0] - a[8];
03041 x1i = a[1] - a[9];
03042 x2r = a[4] + a[12];
03043 x2i = a[5] + a[13];
03044 x3r = a[4] - a[12];
03045 x3i = a[5] - a[13];
03046 y0r = x0r + x2r;
03047 y0i = x0i + x2i;
03048 y2r = x0r - x2r;
03049 y2i = x0i - x2i;
03050 y1r = x1r - x3i;
03051 y1i = x1i + x3r;
03052 y3r = x1r + x3i;
03053 y3i = x1i - x3r;
03054 x0r = a[2] + a[10];
03055 x0i = a[3] + a[11];
03056 x1r = a[2] - a[10];
03057 x1i = a[3] - a[11];
03058 x2r = a[6] + a[14];
03059 x2i = a[7] + a[15];
03060 x3r = a[6] - a[14];
03061 x3i = a[7] - a[15];
03062 y4r = x0r + x2r;
03063 y4i = x0i + x2i;
03064 y6r = x0r - x2r;
03065 y6i = x0i - x2i;
03066 x0r = x1r - x3i;
03067 x0i = x1i + x3r;
03068 x2r = x1r + x3i;
03069 x2i = x1i - x3r;
03070 y5r = wn4r * (x0r - x0i);
03071 y5i = wn4r * (x0r + x0i);
03072 y7r = wn4r * (x2r - x2i);
03073 y7i = wn4r * (x2r + x2i);
03074 a[8] = y1r + y5r;
03075 a[9] = y1i + y5i;
03076 a[10] = y1r - y5r;
03077 a[11] = y1i - y5i;
03078 a[12] = y3r - y7i;
03079 a[13] = y3i + y7r;
03080 a[14] = y3r + y7i;
03081 a[15] = y3i - y7r;
03082 a[0] = y0r + y4r;
03083 a[1] = y0i + y4i;
03084 a[2] = y0r - y4r;
03085 a[3] = y0i - y4i;
03086 a[4] = y2r - y6i;
03087 a[5] = y2i + y6r;
03088 a[6] = y2r + y6i;
03089 a[7] = y2i - y6r;
03090 }
03091
03092
03093 void cftf082(double *a, double *w)
03094 {
03095 double wn4r, wk1r, wk1i, x0r, x0i, x1r, x1i,
03096 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i,
03097 y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i;
03098
03099 wn4r = w[1];
03100 wk1r = w[2];
03101 wk1i = w[3];
03102 y0r = a[0] - a[9];
03103 y0i = a[1] + a[8];
03104 y1r = a[0] + a[9];
03105 y1i = a[1] - a[8];
03106 x0r = a[4] - a[13];
03107 x0i = a[5] + a[12];
03108 y2r = wn4r * (x0r - x0i);
03109 y2i = wn4r * (x0i + x0r);
03110 x0r = a[4] + a[13];
03111 x0i = a[5] - a[12];
03112 y3r = wn4r * (x0r - x0i);
03113 y3i = wn4r * (x0i + x0r);
03114 x0r = a[2] - a[11];
03115 x0i = a[3] + a[10];
03116 y4r = wk1r * x0r - wk1i * x0i;
03117 y4i = wk1r * x0i + wk1i * x0r;
03118 x0r = a[2] + a[11];
03119 x0i = a[3] - a[10];
03120 y5r = wk1i * x0r - wk1r * x0i;
03121 y5i = wk1i * x0i + wk1r * x0r;
03122 x0r = a[6] - a[15];
03123 x0i = a[7] + a[14];
03124 y6r = wk1i * x0r - wk1r * x0i;
03125 y6i = wk1i * x0i + wk1r * x0r;
03126 x0r = a[6] + a[15];
03127 x0i = a[7] - a[14];
03128 y7r = wk1r * x0r - wk1i * x0i;
03129 y7i = wk1r * x0i + wk1i * x0r;
03130 x0r = y0r + y2r;
03131 x0i = y0i + y2i;
03132 x1r = y4r + y6r;
03133 x1i = y4i + y6i;
03134 a[0] = x0r + x1r;
03135 a[1] = x0i + x1i;
03136 a[2] = x0r - x1r;
03137 a[3] = x0i - x1i;
03138 x0r = y0r - y2r;
03139 x0i = y0i - y2i;
03140 x1r = y4r - y6r;
03141 x1i = y4i - y6i;
03142 a[4] = x0r - x1i;
03143 a[5] = x0i + x1r;
03144 a[6] = x0r + x1i;
03145 a[7] = x0i - x1r;
03146 x0r = y1r - y3i;
03147 x0i = y1i + y3r;
03148 x1r = y5r - y7r;
03149 x1i = y5i - y7i;
03150 a[8] = x0r + x1r;
03151 a[9] = x0i + x1i;
03152 a[10] = x0r - x1r;
03153 a[11] = x0i - x1i;
03154 x0r = y1r + y3i;
03155 x0i = y1i - y3r;
03156 x1r = y5r + y7r;
03157 x1i = y5i + y7i;
03158 a[12] = x0r - x1i;
03159 a[13] = x0i + x1r;
03160 a[14] = x0r + x1i;
03161 a[15] = x0i - x1r;
03162 }
03163
03164
03165 void cftf040(double *a)
03166 {
03167 double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
03168
03169 x0r = a[0] + a[4];
03170 x0i = a[1] + a[5];
03171 x1r = a[0] - a[4];
03172 x1i = a[1] - a[5];
03173 x2r = a[2] + a[6];
03174 x2i = a[3] + a[7];
03175 x3r = a[2] - a[6];
03176 x3i = a[3] - a[7];
03177 a[0] = x0r + x2r;
03178 a[1] = x0i + x2i;
03179 a[2] = x1r - x3i;
03180 a[3] = x1i + x3r;
03181 a[4] = x0r - x2r;
03182 a[5] = x0i - x2i;
03183 a[6] = x1r + x3i;
03184 a[7] = x1i - x3r;
03185 }
03186
03187
03188 void cftb040(double *a)
03189 {
03190 double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
03191
03192 x0r = a[0] + a[4];
03193 x0i = a[1] + a[5];
03194 x1r = a[0] - a[4];
03195 x1i = a[1] - a[5];
03196 x2r = a[2] + a[6];
03197 x2i = a[3] + a[7];
03198 x3r = a[2] - a[6];
03199 x3i = a[3] - a[7];
03200 a[0] = x0r + x2r;
03201 a[1] = x0i + x2i;
03202 a[2] = x1r + x3i;
03203 a[3] = x1i - x3r;
03204 a[4] = x0r - x2r;
03205 a[5] = x0i - x2i;
03206 a[6] = x1r - x3i;
03207 a[7] = x1i + x3r;
03208 }
03209
03210
03211 void cftx020(double *a)
03212 {
03213 double x0r, x0i;
03214
03215 x0r = a[0] - a[2];
03216 x0i = a[1] - a[3];
03217 a[0] += a[2];
03218 a[1] += a[3];
03219 a[2] = x0r;
03220 a[3] = x0i;
03221 }
03222
03223
03224 void rftfsub(int n, double *a, int nc, double *c)
03225 {
03226 int j, k, kk, ks, m;
03227 double wkr, wki, xr, xi, yr, yi;
03228
03229 m = n >> 1;
03230 ks = 2 * nc / m;
03231 kk = 0;
03232 for (j = 2; j < m; j += 2) {
03233 k = n - j;
03234 kk += ks;
03235 wkr = 0.5 - c[nc - kk];
03236 wki = c[kk];
03237 xr = a[j] - a[k];
03238 xi = a[j + 1] + a[k + 1];
03239 yr = wkr * xr - wki * xi;
03240 yi = wkr * xi + wki * xr;
03241 a[j] -= yr;
03242 a[j + 1] -= yi;
03243 a[k] += yr;
03244 a[k + 1] -= yi;
03245 }
03246 }
03247
03248
03249 void rftbsub(int n, double *a, int nc, double *c)
03250 {
03251 int j, k, kk, ks, m;
03252 double wkr, wki, xr, xi, yr, yi;
03253
03254 m = n >> 1;
03255 ks = 2 * nc / m;
03256 kk = 0;
03257 for (j = 2; j < m; j += 2) {
03258 k = n - j;
03259 kk += ks;
03260 wkr = 0.5 - c[nc - kk];
03261 wki = c[kk];
03262 xr = a[j] - a[k];
03263 xi = a[j + 1] + a[k + 1];
03264 yr = wkr * xr + wki * xi;
03265 yi = wkr * xi - wki * xr;
03266 a[j] -= yr;
03267 a[j + 1] -= yi;
03268 a[k] += yr;
03269 a[k + 1] -= yi;
03270 }
03271 }
03272
03273
03274 void dctsub(int n, double *a, int nc, double *c)
03275 {
03276 int j, k, kk, ks, m;
03277 double wkr, wki, xr;
03278
03279 m = n >> 1;
03280 ks = nc / n;
03281 kk = 0;
03282 for (j = 1; j < m; j++) {
03283 k = n - j;
03284 kk += ks;
03285 wkr = c[kk] - c[nc - kk];
03286 wki = c[kk] + c[nc - kk];
03287 xr = wki * a[j] - wkr * a[k];
03288 a[j] = wkr * a[j] + wki * a[k];
03289 a[k] = xr;
03290 }
03291 a[m] *= c[0];
03292 }
03293
03294
03295 void dstsub(int n, double *a, int nc, double *c)
03296 {
03297 int j, k, kk, ks, m;
03298 double wkr, wki, xr;
03299
03300 m = n >> 1;
03301 ks = nc / n;
03302 kk = 0;
03303 for (j = 1; j < m; j++) {
03304 k = n - j;
03305 kk += ks;
03306 wkr = c[kk] - c[nc - kk];
03307 wki = c[kk] + c[nc - kk];
03308 xr = wki * a[k] - wkr * a[j];
03309 a[k] = wkr * a[k] + wki * a[j];
03310 a[j] = xr;
03311 }
03312 a[m] *= c[0];
03313 }
03314