bpmdsp/fftsg.c

00001 /*
00002 Fast Fourier/Cosine/Sine Transform
00003     dimension   :one
00004     data length :power of 2
00005     decimation  :frequency
00006     radix       :split-radix
00007     data        :inplace
00008     table       :use
00009 functions
00010     cdft: Complex Discrete Fourier Transform
00011     rdft: Real Discrete Fourier Transform
00012     ddct: Discrete Cosine Transform
00013     ddst: Discrete Sine Transform
00014     dfct: Cosine Transform of RDFT (Real Symmetric DFT)
00015     dfst: Sine Transform of RDFT (Real Anti-symmetric DFT)
00016 function prototypes
00017     void cdft(int, int, double *, int *, double *);
00018     void rdft(int, int, double *, int *, double *);
00019     void ddct(int, int, double *, int *, double *);
00020     void ddst(int, int, double *, int *, double *);
00021     void dfct(int, double *, double *, int *, double *);
00022     void dfst(int, double *, double *, int *, double *);
00023 macro definitions
00024     USE_CDFT_PTHREADS : default=not defined
00025         CDFT_THREADS_BEGIN_N  : must be >= 512, default=8192
00026         CDFT_4THREADS_BEGIN_N : must be >= 512, default=65536
00027     USE_CDFT_WINTHREADS : default=not defined
00028         CDFT_THREADS_BEGIN_N  : must be >= 512, default=32768
00029         CDFT_4THREADS_BEGIN_N : must be >= 512, default=524288
00030 
00031 
00032 -------- Complex DFT (Discrete Fourier Transform) --------
00033     [definition]
00034         <case1>
00035             X[k] = sum_j=0^n-1 x[j]*exp(2*pi*i*j*k/n), 0<=k<n
00036         <case2>
00037             X[k] = sum_j=0^n-1 x[j]*exp(-2*pi*i*j*k/n), 0<=k<n
00038         (notes: sum_j=0^n-1 is a summation from j=0 to n-1)
00039     [usage]
00040         <case1>
00041             ip[0] = 0; // first time only
00042             cdft(2*n, 1, a, ip, w);
00043         <case2>
00044             ip[0] = 0; // first time only
00045             cdft(2*n, -1, a, ip, w);
00046     [parameters]
00047         2*n            :data length (int)
00048                         n >= 1, n = power of 2
00049         a[0...2*n-1]   :input/output data (double *)
00050                         input data
00051                             a[2*j] = Re(x[j]), 
00052                             a[2*j+1] = Im(x[j]), 0<=j<n
00053                         output data
00054                             a[2*k] = Re(X[k]), 
00055                             a[2*k+1] = Im(X[k]), 0<=k<n
00056         ip[0...*]      :work area for bit reversal (int *)
00057                         length of ip >= 2+sqrt(n)
00058                         strictly, 
00059                         length of ip >= 
00060                             2+(1<<(int)(log(n+0.5)/log(2))/2).
00061                         ip[0],ip[1] are pointers of the cos/sin table.
00062         w[0...n/2-1]   :cos/sin table (double *)
00063                         w[],ip[] are initialized if ip[0] == 0.
00064     [remark]
00065         Inverse of 
00066             cdft(2*n, -1, a, ip, w);
00067         is 
00068             cdft(2*n, 1, a, ip, w);
00069             for (j = 0; j <= 2 * n - 1; j++) {
00070                 a[j] *= 1.0 / n;
00071             }
00072         .
00073 
00074 
00075 -------- Real DFT / Inverse of Real DFT --------
00076     [definition]
00077         <case1> RDFT
00078             R[k] = sum_j=0^n-1 a[j]*cos(2*pi*j*k/n), 0<=k<=n/2
00079             I[k] = sum_j=0^n-1 a[j]*sin(2*pi*j*k/n), 0<k<n/2
00080         <case2> IRDFT (excluding scale)
00081             a[k] = (R[0] + R[n/2]*cos(pi*k))/2 + 
00082                    sum_j=1^n/2-1 R[j]*cos(2*pi*j*k/n) + 
00083                    sum_j=1^n/2-1 I[j]*sin(2*pi*j*k/n), 0<=k<n
00084     [usage]
00085         <case1>
00086             ip[0] = 0; // first time only
00087             rdft(n, 1, a, ip, w);
00088         <case2>
00089             ip[0] = 0; // first time only
00090             rdft(n, -1, a, ip, w);
00091     [parameters]
00092         n              :data length (int)
00093                         n >= 2, n = power of 2
00094         a[0...n-1]     :input/output data (double *)
00095                         <case1>
00096                             output data
00097                                 a[2*k] = R[k], 0<=k<n/2
00098                                 a[2*k+1] = I[k], 0<k<n/2
00099                                 a[1] = R[n/2]
00100                         <case2>
00101                             input data
00102                                 a[2*j] = R[j], 0<=j<n/2
00103                                 a[2*j+1] = I[j], 0<j<n/2
00104                                 a[1] = R[n/2]
00105         ip[0...*]      :work area for bit reversal (int *)
00106                         length of ip >= 2+sqrt(n/2)
00107                         strictly, 
00108                         length of ip >= 
00109                             2+(1<<(int)(log(n/2+0.5)/log(2))/2).
00110                         ip[0],ip[1] are pointers of the cos/sin table.
00111         w[0...n/2-1]   :cos/sin table (double *)
00112                         w[],ip[] are initialized if ip[0] == 0.
00113     [remark]
00114         Inverse of 
00115             rdft(n, 1, a, ip, w);
00116         is 
00117             rdft(n, -1, a, ip, w);
00118             for (j = 0; j <= n - 1; j++) {
00119                 a[j] *= 2.0 / n;
00120             }
00121         .
00122 
00123 
00124 -------- DCT (Discrete Cosine Transform) / Inverse of DCT --------
00125     [definition]
00126         <case1> IDCT (excluding scale)
00127             C[k] = sum_j=0^n-1 a[j]*cos(pi*j*(k+1/2)/n), 0<=k<n
00128         <case2> DCT
00129             C[k] = sum_j=0^n-1 a[j]*cos(pi*(j+1/2)*k/n), 0<=k<n
00130     [usage]
00131         <case1>
00132             ip[0] = 0; // first time only
00133             ddct(n, 1, a, ip, w);
00134         <case2>
00135             ip[0] = 0; // first time only
00136             ddct(n, -1, a, ip, w);
00137     [parameters]
00138         n              :data length (int)
00139                         n >= 2, n = power of 2
00140         a[0...n-1]     :input/output data (double *)
00141                         output data
00142                             a[k] = C[k], 0<=k<n
00143         ip[0...*]      :work area for bit reversal (int *)
00144                         length of ip >= 2+sqrt(n/2)
00145                         strictly, 
00146                         length of ip >= 
00147                             2+(1<<(int)(log(n/2+0.5)/log(2))/2).
00148                         ip[0],ip[1] are pointers of the cos/sin table.
00149         w[0...n*5/4-1] :cos/sin table (double *)
00150                         w[],ip[] are initialized if ip[0] == 0.
00151     [remark]
00152         Inverse of 
00153             ddct(n, -1, a, ip, w);
00154         is 
00155             a[0] *= 0.5;
00156             ddct(n, 1, a, ip, w);
00157             for (j = 0; j <= n - 1; j++) {
00158                 a[j] *= 2.0 / n;
00159             }
00160         .
00161 
00162 
00163 -------- DST (Discrete Sine Transform) / Inverse of DST --------
00164     [definition]
00165         <case1> IDST (excluding scale)
00166             S[k] = sum_j=1^n A[j]*sin(pi*j*(k+1/2)/n), 0<=k<n
00167         <case2> DST
00168             S[k] = sum_j=0^n-1 a[j]*sin(pi*(j+1/2)*k/n), 0<k<=n
00169     [usage]
00170         <case1>
00171             ip[0] = 0; // first time only
00172             ddst(n, 1, a, ip, w);
00173         <case2>
00174             ip[0] = 0; // first time only
00175             ddst(n, -1, a, ip, w);
00176     [parameters]
00177         n              :data length (int)
00178                         n >= 2, n = power of 2
00179         a[0...n-1]     :input/output data (double *)
00180                         <case1>
00181                             input data
00182                                 a[j] = A[j], 0<j<n
00183                                 a[0] = A[n]
00184                             output data
00185                                 a[k] = S[k], 0<=k<n
00186                         <case2>
00187                             output data
00188                                 a[k] = S[k], 0<k<n
00189                                 a[0] = S[n]
00190         ip[0...*]      :work area for bit reversal (int *)
00191                         length of ip >= 2+sqrt(n/2)
00192                         strictly, 
00193                         length of ip >= 
00194                             2+(1<<(int)(log(n/2+0.5)/log(2))/2).
00195                         ip[0],ip[1] are pointers of the cos/sin table.
00196         w[0...n*5/4-1] :cos/sin table (double *)
00197                         w[],ip[] are initialized if ip[0] == 0.
00198     [remark]
00199         Inverse of 
00200             ddst(n, -1, a, ip, w);
00201         is 
00202             a[0] *= 0.5;
00203             ddst(n, 1, a, ip, w);
00204             for (j = 0; j <= n - 1; j++) {
00205                 a[j] *= 2.0 / n;
00206             }
00207         .
00208 
00209 
00210 -------- Cosine Transform of RDFT (Real Symmetric DFT) --------
00211     [definition]
00212         C[k] = sum_j=0^n a[j]*cos(pi*j*k/n), 0<=k<=n
00213     [usage]
00214         ip[0] = 0; // first time only
00215         dfct(n, a, t, ip, w);
00216     [parameters]
00217         n              :data length - 1 (int)
00218                         n >= 2, n = power of 2
00219         a[0...n]       :input/output data (double *)
00220                         output data
00221                             a[k] = C[k], 0<=k<=n
00222         t[0...n/2]     :work area (double *)
00223         ip[0...*]      :work area for bit reversal (int *)
00224                         length of ip >= 2+sqrt(n/4)
00225                         strictly, 
00226                         length of ip >= 
00227                             2+(1<<(int)(log(n/4+0.5)/log(2))/2).
00228                         ip[0],ip[1] are pointers of the cos/sin table.
00229         w[0...n*5/8-1] :cos/sin table (double *)
00230                         w[],ip[] are initialized if ip[0] == 0.
00231     [remark]
00232         Inverse of 
00233             a[0] *= 0.5;
00234             a[n] *= 0.5;
00235             dfct(n, a, t, ip, w);
00236         is 
00237             a[0] *= 0.5;
00238             a[n] *= 0.5;
00239             dfct(n, a, t, ip, w);
00240             for (j = 0; j <= n; j++) {
00241                 a[j] *= 2.0 / n;
00242             }
00243         .
00244 
00245 
00246 -------- Sine Transform of RDFT (Real Anti-symmetric DFT) --------
00247     [definition]
00248         S[k] = sum_j=1^n-1 a[j]*sin(pi*j*k/n), 0<k<n
00249     [usage]
00250         ip[0] = 0; // first time only
00251         dfst(n, a, t, ip, w);
00252     [parameters]
00253         n              :data length + 1 (int)
00254                         n >= 2, n = power of 2
00255         a[0...n-1]     :input/output data (double *)
00256                         output data
00257                             a[k] = S[k], 0<k<n
00258                         (a[0] is used for work area)
00259         t[0...n/2-1]   :work area (double *)
00260         ip[0...*]      :work area for bit reversal (int *)
00261                         length of ip >= 2+sqrt(n/4)
00262                         strictly, 
00263                         length of ip >= 
00264                             2+(1<<(int)(log(n/4+0.5)/log(2))/2).
00265                         ip[0],ip[1] are pointers of the cos/sin table.
00266         w[0...n*5/8-1] :cos/sin table (double *)
00267                         w[],ip[] are initialized if ip[0] == 0.
00268     [remark]
00269         Inverse of 
00270             dfst(n, a, t, ip, w);
00271         is 
00272             dfst(n, a, t, ip, w);
00273             for (j = 1; j <= n - 1; j++) {
00274                 a[j] *= 2.0 / n;
00275             }
00276         .
00277 
00278 
00279 Appendix :
00280     The cos/sin table is recalculated when the larger table required.
00281     w[] and ip[] are compatible with all routines.
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 /* -------- initializing routines -------- */
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 /* -------- child routines -------- */
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 /* USE_CDFT_PTHREADS */
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 /* USE_CDFT_WINTHREADS */
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 /* USE_CDFT_THREADS */
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 /* USE_CDFT_THREADS */
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 /* USE_CDFT_THREADS */
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 /* USE_CDFT_THREADS */
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 /* USE_CDFT_THREADS */
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 

Generated on Thu Jan 10 10:18:04 2008 for libbpm by  doxygen 1.5.1