Main Page | Modules | Class Hierarchy | Class List | Directories | File List | Class Members | File Members | Related Pages

fftsg.c

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

Generated on Tue Dec 20 10:14:57 2005 for vlc-0.8.4a by  doxygen 1.4.2