source: src/spectral/ooura_fft8g.c @ 6bbdcff

feature/autosinkfeature/cnnfeature/cnn_orgfeature/constantqfeature/crepefeature/crepe_orgfeature/pitchshiftfeature/pydocstringsfeature/timestretchfix/ffmpeg5pitchshiftsamplertimestretchyinfft+
Last change on this file since 6bbdcff was cfaa3c4, checked in by Paul Brossier <piem@piem.org>, 13 years ago

src/spectral/ooura_fft8g.c: add ooura
  • Property mode set to 100644
File size: 46.8 KB
Line 
1/*
2Fast Fourier/Cosine/Sine Transform
3    dimension   :one
4    data length :power of 2
5    decimation  :frequency
6    radix       :8, 4, 2
7    data        :inplace
8    table       :use
9functions
10    cdft: Complex Discrete Fourier Transform
11    rdft: Real Discrete Fourier Transform
12    ddct: Discrete Cosine Transform
13    ddst: Discrete Sine Transform
14    dfct: Cosine Transform of RDFT (Real Symmetric DFT)
15    dfst: Sine Transform of RDFT (Real Anti-symmetric DFT)
16function prototypes
17    void cdft(int, int, double *, int *, double *);
18    void rdft(int, int, double *, int *, double *);
19    void ddct(int, int, double *, int *, double *);
20    void ddst(int, int, double *, int *, double *);
21    void dfct(int, double *, double *, int *, double *);
22    void dfst(int, double *, double *, int *, double *);
23
24
25-------- Complex DFT (Discrete Fourier Transform) --------
26    [definition]
27        <case1>
28            X[k] = sum_j=0^n-1 x[j]*exp(2*pi*i*j*k/n), 0<=k<n
29        <case2>
30            X[k] = sum_j=0^n-1 x[j]*exp(-2*pi*i*j*k/n), 0<=k<n
31        (notes: sum_j=0^n-1 is a summation from j=0 to n-1)
32    [usage]
33        <case1>
34            ip[0] = 0; // first time only
35            cdft(2*n, 1, a, ip, w);
36        <case2>
37            ip[0] = 0; // first time only
38            cdft(2*n, -1, a, ip, w);
39    [parameters]
40        2*n            :data length (int)
41                        n >= 1, n = power of 2
42        a[0...2*n-1]   :input/output data (double *)
43                        input data
44                            a[2*j] = Re(x[j]),
45                            a[2*j+1] = Im(x[j]), 0<=j<n
46                        output data
47                            a[2*k] = Re(X[k]),
48                            a[2*k+1] = Im(X[k]), 0<=k<n
49        ip[0...*]      :work area for bit reversal (int *)
50                        length of ip >= 2+sqrt(n)
51                        strictly,
52                        length of ip >=
53                            2+(1<<(int)(log(n+0.5)/log(2))/2).
54                        ip[0],ip[1] are pointers of the cos/sin table.
55        w[0...n/2-1]   :cos/sin table (double *)
56                        w[],ip[] are initialized if ip[0] == 0.
57    [remark]
58        Inverse of
59            cdft(2*n, -1, a, ip, w);
60        is
61            cdft(2*n, 1, a, ip, w);
62            for (j = 0; j <= 2 * n - 1; j++) {
63                a[j] *= 1.0 / n;
64            }
65        .
66
67
68-------- Real DFT / Inverse of Real DFT --------
69    [definition]
70        <case1> RDFT
71            R[k] = sum_j=0^n-1 a[j]*cos(2*pi*j*k/n), 0<=k<=n/2
72            I[k] = sum_j=0^n-1 a[j]*sin(2*pi*j*k/n), 0<k<n/2
73        <case2> IRDFT (excluding scale)
74            a[k] = (R[0] + R[n/2]*cos(pi*k))/2 +
75                   sum_j=1^n/2-1 R[j]*cos(2*pi*j*k/n) +
76                   sum_j=1^n/2-1 I[j]*sin(2*pi*j*k/n), 0<=k<n
77    [usage]
78        <case1>
79            ip[0] = 0; // first time only
80            rdft(n, 1, a, ip, w);
81        <case2>
82            ip[0] = 0; // first time only
83            rdft(n, -1, a, ip, w);
84    [parameters]
85        n              :data length (int)
86                        n >= 2, n = power of 2
87        a[0...n-1]     :input/output data (double *)
88                        <case1>
89                            output data
90                                a[2*k] = R[k], 0<=k<n/2
91                                a[2*k+1] = I[k], 0<k<n/2
92                                a[1] = R[n/2]
93                        <case2>
94                            input data
95                                a[2*j] = R[j], 0<=j<n/2
96                                a[2*j+1] = I[j], 0<j<n/2
97                                a[1] = R[n/2]
98        ip[0...*]      :work area for bit reversal (int *)
99                        length of ip >= 2+sqrt(n/2)
100                        strictly,
101                        length of ip >=
102                            2+(1<<(int)(log(n/2+0.5)/log(2))/2).
103                        ip[0],ip[1] are pointers of the cos/sin table.
104        w[0...n/2-1]   :cos/sin table (double *)
105                        w[],ip[] are initialized if ip[0] == 0.
106    [remark]
107        Inverse of
108            rdft(n, 1, a, ip, w);
109        is
110            rdft(n, -1, a, ip, w);
111            for (j = 0; j <= n - 1; j++) {
112                a[j] *= 2.0 / n;
113            }
114        .
115
116
117-------- DCT (Discrete Cosine Transform) / Inverse of DCT --------
118    [definition]
119        <case1> IDCT (excluding scale)
120            C[k] = sum_j=0^n-1 a[j]*cos(pi*j*(k+1/2)/n), 0<=k<n
121        <case2> DCT
122            C[k] = sum_j=0^n-1 a[j]*cos(pi*(j+1/2)*k/n), 0<=k<n
123    [usage]
124        <case1>
125            ip[0] = 0; // first time only
126            ddct(n, 1, a, ip, w);
127        <case2>
128            ip[0] = 0; // first time only
129            ddct(n, -1, a, ip, w);
130    [parameters]
131        n              :data length (int)
132                        n >= 2, n = power of 2
133        a[0...n-1]     :input/output data (double *)
134                        output data
135                            a[k] = C[k], 0<=k<n
136        ip[0...*]      :work area for bit reversal (int *)
137                        length of ip >= 2+sqrt(n/2)
138                        strictly,
139                        length of ip >=
140                            2+(1<<(int)(log(n/2+0.5)/log(2))/2).
141                        ip[0],ip[1] are pointers of the cos/sin table.
142        w[0...n*5/4-1] :cos/sin table (double *)
143                        w[],ip[] are initialized if ip[0] == 0.
144    [remark]
145        Inverse of
146            ddct(n, -1, a, ip, w);
147        is
148            a[0] *= 0.5;
149            ddct(n, 1, a, ip, w);
150            for (j = 0; j <= n - 1; j++) {
151                a[j] *= 2.0 / n;
152            }
153        .
154
155
156-------- DST (Discrete Sine Transform) / Inverse of DST --------
157    [definition]
158        <case1> IDST (excluding scale)
159            S[k] = sum_j=1^n A[j]*sin(pi*j*(k+1/2)/n), 0<=k<n
160        <case2> DST
161            S[k] = sum_j=0^n-1 a[j]*sin(pi*(j+1/2)*k/n), 0<k<=n
162    [usage]
163        <case1>
164            ip[0] = 0; // first time only
165            ddst(n, 1, a, ip, w);
166        <case2>
167            ip[0] = 0; // first time only
168            ddst(n, -1, a, ip, w);
169    [parameters]
170        n              :data length (int)
171                        n >= 2, n = power of 2
172        a[0...n-1]     :input/output data (double *)
173                        <case1>
174                            input data
175                                a[j] = A[j], 0<j<n
176                                a[0] = A[n]
177                            output data
178                                a[k] = S[k], 0<=k<n
179                        <case2>
180                            output data
181                                a[k] = S[k], 0<k<n
182                                a[0] = S[n]
183        ip[0...*]      :work area for bit reversal (int *)
184                        length of ip >= 2+sqrt(n/2)
185                        strictly,
186                        length of ip >=
187                            2+(1<<(int)(log(n/2+0.5)/log(2))/2).
188                        ip[0],ip[1] are pointers of the cos/sin table.
189        w[0...n*5/4-1] :cos/sin table (double *)
190                        w[],ip[] are initialized if ip[0] == 0.
191    [remark]
192        Inverse of
193            ddst(n, -1, a, ip, w);
194        is
195            a[0] *= 0.5;
196            ddst(n, 1, a, ip, w);
197            for (j = 0; j <= n - 1; j++) {
198                a[j] *= 2.0 / n;
199            }
200        .
201
202
203-------- Cosine Transform of RDFT (Real Symmetric DFT) --------
204    [definition]
205        C[k] = sum_j=0^n a[j]*cos(pi*j*k/n), 0<=k<=n
206    [usage]
207        ip[0] = 0; // first time only
208        dfct(n, a, t, ip, w);
209    [parameters]
210        n              :data length - 1 (int)
211                        n >= 2, n = power of 2
212        a[0...n]       :input/output data (double *)
213                        output data
214                            a[k] = C[k], 0<=k<=n
215        t[0...n/2]     :work area (double *)
216        ip[0...*]      :work area for bit reversal (int *)
217                        length of ip >= 2+sqrt(n/4)
218                        strictly,
219                        length of ip >=
220                            2+(1<<(int)(log(n/4+0.5)/log(2))/2).
221                        ip[0],ip[1] are pointers of the cos/sin table.
222        w[0...n*5/8-1] :cos/sin table (double *)
223                        w[],ip[] are initialized if ip[0] == 0.
224    [remark]
225        Inverse of
226            a[0] *= 0.5;
227            a[n] *= 0.5;
228            dfct(n, a, t, ip, w);
229        is
230            a[0] *= 0.5;
231            a[n] *= 0.5;
232            dfct(n, a, t, ip, w);
233            for (j = 0; j <= n; j++) {
234                a[j] *= 2.0 / n;
235            }
236        .
237
238
239-------- Sine Transform of RDFT (Real Anti-symmetric DFT) --------
240    [definition]
241        S[k] = sum_j=1^n-1 a[j]*sin(pi*j*k/n), 0<k<n
242    [usage]
243        ip[0] = 0; // first time only
244        dfst(n, a, t, ip, w);
245    [parameters]
246        n              :data length + 1 (int)
247                        n >= 2, n = power of 2
248        a[0...n-1]     :input/output data (double *)
249                        output data
250                            a[k] = S[k], 0<k<n
251                        (a[0] is used for work area)
252        t[0...n/2-1]   :work area (double *)
253        ip[0...*]      :work area for bit reversal (int *)
254                        length of ip >= 2+sqrt(n/4)
255                        strictly,
256                        length of ip >=
257                            2+(1<<(int)(log(n/4+0.5)/log(2))/2).
258                        ip[0],ip[1] are pointers of the cos/sin table.
259        w[0...n*5/8-1] :cos/sin table (double *)
260                        w[],ip[] are initialized if ip[0] == 0.
261    [remark]
262        Inverse of
263            dfst(n, a, t, ip, w);
264        is
265            dfst(n, a, t, ip, w);
266            for (j = 1; j <= n - 1; j++) {
267                a[j] *= 2.0 / n;
268            }
269        .
270
271
272Appendix :
273    The cos/sin table is recalculated when the larger table required.
274    w[] and ip[] are compatible with all routines.
275*/
276
277
278void cdft(int n, int isgn, double *a, int *ip, double *w)
279{
280    void makewt(int nw, int *ip, double *w);
281    void bitrv2(int n, int *ip, double *a);
282    void bitrv2conj(int n, int *ip, double *a);
283    void cftfsub(int n, double *a, double *w);
284    void cftbsub(int n, double *a, double *w);
285   
286    if (n > (ip[0] << 2)) {
287        makewt(n >> 2, ip, w);
288    }
289    if (n > 4) {
290        if (isgn >= 0) {
291            bitrv2(n, ip + 2, a);
292            cftfsub(n, a, w);
293        } else {
294            bitrv2conj(n, ip + 2, a);
295            cftbsub(n, a, w);
296        }
297    } else if (n == 4) {
298        cftfsub(n, a, w);
299    }
300}
301
302
303void rdft(int n, int isgn, double *a, int *ip, double *w)
304{
305    void makewt(int nw, int *ip, double *w);
306    void makect(int nc, int *ip, double *c);
307    void bitrv2(int n, int *ip, double *a);
308    void cftfsub(int n, double *a, double *w);
309    void cftbsub(int n, double *a, double *w);
310    void rftfsub(int n, double *a, int nc, double *c);
311    void rftbsub(int n, double *a, int nc, double *c);
312    int nw, nc;
313    double xi;
314   
315    nw = ip[0];
316    if (n > (nw << 2)) {
317        nw = n >> 2;
318        makewt(nw, ip, w);
319    }
320    nc = ip[1];
321    if (n > (nc << 2)) {
322        nc = n >> 2;
323        makect(nc, ip, w + nw);
324    }
325    if (isgn >= 0) {
326        if (n > 4) {
327            bitrv2(n, ip + 2, a);
328            cftfsub(n, a, w);
329            rftfsub(n, a, nc, w + nw);
330        } else if (n == 4) {
331            cftfsub(n, a, w);
332        }
333        xi = a[0] - a[1];
334        a[0] += a[1];
335        a[1] = xi;
336    } else {
337        a[1] = 0.5 * (a[0] - a[1]);
338        a[0] -= a[1];
339        if (n > 4) {
340            rftbsub(n, a, nc, w + nw);
341            bitrv2(n, ip + 2, a);
342            cftbsub(n, a, w);
343        } else if (n == 4) {
344            cftfsub(n, a, w);
345        }
346    }
347}
348
349
350void ddct(int n, int isgn, double *a, int *ip, double *w)
351{
352    void makewt(int nw, int *ip, double *w);
353    void makect(int nc, int *ip, double *c);
354    void bitrv2(int n, int *ip, double *a);
355    void cftfsub(int n, double *a, double *w);
356    void cftbsub(int n, double *a, double *w);
357    void rftfsub(int n, double *a, int nc, double *c);
358    void rftbsub(int n, double *a, int nc, double *c);
359    void dctsub(int n, double *a, int nc, double *c);
360    int j, nw, nc;
361    double xr;
362   
363    nw = ip[0];
364    if (n > (nw << 2)) {
365        nw = n >> 2;
366        makewt(nw, ip, w);
367    }
368    nc = ip[1];
369    if (n > nc) {
370        nc = n;
371        makect(nc, ip, w + nw);
372    }
373    if (isgn < 0) {
374        xr = a[n - 1];
375        for (j = n - 2; j >= 2; j -= 2) {
376            a[j + 1] = a[j] - a[j - 1];
377            a[j] += a[j - 1];
378        }
379        a[1] = a[0] - xr;
380        a[0] += xr;
381        if (n > 4) {
382            rftbsub(n, a, nc, w + nw);
383            bitrv2(n, ip + 2, a);
384            cftbsub(n, a, w);
385        } else if (n == 4) {
386            cftfsub(n, a, w);
387        }
388    }
389    dctsub(n, a, nc, w + nw);
390    if (isgn >= 0) {
391        if (n > 4) {
392            bitrv2(n, ip + 2, a);
393            cftfsub(n, a, w);
394            rftfsub(n, a, nc, w + nw);
395        } else if (n == 4) {
396            cftfsub(n, a, w);
397        }
398        xr = a[0] - a[1];
399        a[0] += a[1];
400        for (j = 2; j < n; j += 2) {
401            a[j - 1] = a[j] - a[j + 1];
402            a[j] += a[j + 1];
403        }
404        a[n - 1] = xr;
405    }
406}
407
408
409void ddst(int n, int isgn, double *a, int *ip, double *w)
410{
411    void makewt(int nw, int *ip, double *w);
412    void makect(int nc, int *ip, double *c);
413    void bitrv2(int n, int *ip, double *a);
414    void cftfsub(int n, double *a, double *w);
415    void cftbsub(int n, double *a, double *w);
416    void rftfsub(int n, double *a, int nc, double *c);
417    void rftbsub(int n, double *a, int nc, double *c);
418    void dstsub(int n, double *a, int nc, double *c);
419    int j, nw, nc;
420    double xr;
421   
422    nw = ip[0];
423    if (n > (nw << 2)) {
424        nw = n >> 2;
425        makewt(nw, ip, w);
426    }
427    nc = ip[1];
428    if (n > nc) {
429        nc = n;
430        makect(nc, ip, w + nw);
431    }
432    if (isgn < 0) {
433        xr = a[n - 1];
434        for (j = n - 2; j >= 2; j -= 2) {
435            a[j + 1] = -a[j] - a[j - 1];
436            a[j] -= a[j - 1];
437        }
438        a[1] = a[0] + xr;
439        a[0] -= xr;
440        if (n > 4) {
441            rftbsub(n, a, nc, w + nw);
442            bitrv2(n, ip + 2, a);
443            cftbsub(n, a, w);
444        } else if (n == 4) {
445            cftfsub(n, a, w);
446        }
447    }
448    dstsub(n, a, nc, w + nw);
449    if (isgn >= 0) {
450        if (n > 4) {
451            bitrv2(n, ip + 2, a);
452            cftfsub(n, a, w);
453            rftfsub(n, a, nc, w + nw);
454        } else if (n == 4) {
455            cftfsub(n, a, w);
456        }
457        xr = a[0] - a[1];
458        a[0] += a[1];
459        for (j = 2; j < n; j += 2) {
460            a[j - 1] = -a[j] - a[j + 1];
461            a[j] -= a[j + 1];
462        }
463        a[n - 1] = -xr;
464    }
465}
466
467
468void dfct(int n, double *a, double *t, int *ip, double *w)
469{
470    void makewt(int nw, int *ip, double *w);
471    void makect(int nc, int *ip, double *c);
472    void bitrv2(int n, int *ip, double *a);
473    void cftfsub(int n, double *a, double *w);
474    void rftfsub(int n, double *a, int nc, double *c);
475    void dctsub(int n, double *a, int nc, double *c);
476    int j, k, l, m, mh, nw, nc;
477    double xr, xi, yr, yi;
478   
479    nw = ip[0];
480    if (n > (nw << 3)) {
481        nw = n >> 3;
482        makewt(nw, ip, w);
483    }
484    nc = ip[1];
485    if (n > (nc << 1)) {
486        nc = n >> 1;
487        makect(nc, ip, w + nw);
488    }
489    m = n >> 1;
490    yi = a[m];
491    xi = a[0] + a[n];
492    a[0] -= a[n];
493    t[0] = xi - yi;
494    t[m] = xi + yi;
495    if (n > 2) {
496        mh = m >> 1;
497        for (j = 1; j < mh; j++) {
498            k = m - j;
499            xr = a[j] - a[n - j];
500            xi = a[j] + a[n - j];
501            yr = a[k] - a[n - k];
502            yi = a[k] + a[n - k];
503            a[j] = xr;
504            a[k] = yr;
505            t[j] = xi - yi;
506            t[k] = xi + yi;
507        }
508        t[mh] = a[mh] + a[n - mh];
509        a[mh] -= a[n - mh];
510        dctsub(m, a, nc, w + nw);
511        if (m > 4) {
512            bitrv2(m, ip + 2, a);
513            cftfsub(m, a, w);
514            rftfsub(m, a, nc, w + nw);
515        } else if (m == 4) {
516            cftfsub(m, a, w);
517        }
518        a[n - 1] = a[0] - a[1];
519        a[1] = a[0] + a[1];
520        for (j = m - 2; j >= 2; j -= 2) {
521            a[2 * j + 1] = a[j] + a[j + 1];
522            a[2 * j - 1] = a[j] - a[j + 1];
523        }
524        l = 2;
525        m = mh;
526        while (m >= 2) {
527            dctsub(m, t, nc, w + nw);
528            if (m > 4) {
529                bitrv2(m, ip + 2, t);
530                cftfsub(m, t, w);
531                rftfsub(m, t, nc, w + nw);
532            } else if (m == 4) {
533                cftfsub(m, t, w);
534            }
535            a[n - l] = t[0] - t[1];
536            a[l] = t[0] + t[1];
537            k = 0;
538            for (j = 2; j < m; j += 2) {
539                k += l << 2;
540                a[k - l] = t[j] - t[j + 1];
541                a[k + l] = t[j] + t[j + 1];
542            }
543            l <<= 1;
544            mh = m >> 1;
545            for (j = 0; j < mh; j++) {
546                k = m - j;
547                t[j] = t[m + k] - t[m + j];
548                t[k] = t[m + k] + t[m + j];
549            }
550            t[mh] = t[m + mh];
551            m = mh;
552        }
553        a[l] = t[0];
554        a[n] = t[2] - t[1];
555        a[0] = t[2] + t[1];
556    } else {
557        a[1] = a[0];
558        a[2] = t[0];
559        a[0] = t[1];
560    }
561}
562
563
564void dfst(int n, double *a, double *t, int *ip, double *w)
565{
566    void makewt(int nw, int *ip, double *w);
567    void makect(int nc, int *ip, double *c);
568    void bitrv2(int n, int *ip, double *a);
569    void cftfsub(int n, double *a, double *w);
570    void rftfsub(int n, double *a, int nc, double *c);
571    void dstsub(int n, double *a, int nc, double *c);
572    int j, k, l, m, mh, nw, nc;
573    double xr, xi, yr, yi;
574   
575    nw = ip[0];
576    if (n > (nw << 3)) {
577        nw = n >> 3;
578        makewt(nw, ip, w);
579    }
580    nc = ip[1];
581    if (n > (nc << 1)) {
582        nc = n >> 1;
583        makect(nc, ip, w + nw);
584    }
585    if (n > 2) {
586        m = n >> 1;
587        mh = m >> 1;
588        for (j = 1; j < mh; j++) {
589            k = m - j;
590            xr = a[j] + a[n - j];
591            xi = a[j] - a[n - j];
592            yr = a[k] + a[n - k];
593            yi = a[k] - a[n - k];
594            a[j] = xr;
595            a[k] = yr;
596            t[j] = xi + yi;
597            t[k] = xi - yi;
598        }
599        t[0] = a[mh] - a[n - mh];
600        a[mh] += a[n - mh];
601        a[0] = a[m];
602        dstsub(m, a, nc, w + nw);
603        if (m > 4) {
604            bitrv2(m, ip + 2, a);
605            cftfsub(m, a, w);
606            rftfsub(m, a, nc, w + nw);
607        } else if (m == 4) {
608            cftfsub(m, a, w);
609        }
610        a[n - 1] = a[1] - a[0];
611        a[1] = a[0] + a[1];
612        for (j = m - 2; j >= 2; j -= 2) {
613            a[2 * j + 1] = a[j] - a[j + 1];
614            a[2 * j - 1] = -a[j] - a[j + 1];
615        }
616        l = 2;
617        m = mh;
618        while (m >= 2) {
619            dstsub(m, t, nc, w + nw);
620            if (m > 4) {
621                bitrv2(m, ip + 2, t);
622                cftfsub(m, t, w);
623                rftfsub(m, t, nc, w + nw);
624            } else if (m == 4) {
625                cftfsub(m, t, w);
626            }
627            a[n - l] = t[1] - t[0];
628            a[l] = t[0] + t[1];
629            k = 0;
630            for (j = 2; j < m; j += 2) {
631                k += l << 2;
632                a[k - l] = -t[j] - t[j + 1];
633                a[k + l] = t[j] - t[j + 1];
634            }
635            l <<= 1;
636            mh = m >> 1;
637            for (j = 1; j < mh; j++) {
638                k = m - j;
639                t[j] = t[m + k] + t[m + j];
640                t[k] = t[m + k] - t[m + j];
641            }
642            t[0] = t[m + mh];
643            m = mh;
644        }
645        a[l] = t[0];
646    }
647    a[0] = 0;
648}
649
650
651/* -------- initializing routines -------- */
652
653
654#include <math.h>
655
656void makewt(int nw, int *ip, double *w)
657{
658    void bitrv2(int n, int *ip, double *a);
659    int j, nwh;
660    double delta, x, y;
661   
662    ip[0] = nw;
663    ip[1] = 1;
664    if (nw > 2) {
665        nwh = nw >> 1;
666        delta = atan(1.0) / nwh;
667        w[0] = 1;
668        w[1] = 0;
669        w[nwh] = cos(delta * nwh);
670        w[nwh + 1] = w[nwh];
671        if (nwh > 2) {
672            for (j = 2; j < nwh; j += 2) {
673                x = cos(delta * j);
674                y = sin(delta * j);
675                w[j] = x;
676                w[j + 1] = y;
677                w[nw - j] = y;
678                w[nw - j + 1] = x;
679            }
680            for (j = nwh - 2; j >= 2; j -= 2) {
681                x = w[2 * j];
682                y = w[2 * j + 1];
683                w[nwh + j] = x;
684                w[nwh + j + 1] = y;
685            }
686            bitrv2(nw, ip + 2, w);
687        }
688    }
689}
690
691
692void makect(int nc, int *ip, double *c)
693{
694    int j, nch;
695    double delta;
696   
697    ip[1] = nc;
698    if (nc > 1) {
699        nch = nc >> 1;
700        delta = atan(1.0) / nch;
701        c[0] = cos(delta * nch);
702        c[nch] = 0.5 * c[0];
703        for (j = 1; j < nch; j++) {
704            c[j] = 0.5 * cos(delta * j);
705            c[nc - j] = 0.5 * sin(delta * j);
706        }
707    }
708}
709
710
711/* -------- child routines -------- */
712
713
714void bitrv2(int n, int *ip, double *a)
715{
716    int j, j1, k, k1, l, m, m2;
717    double xr, xi, yr, yi;
718   
719    ip[0] = 0;
720    l = n;
721    m = 1;
722    while ((m << 3) < l) {
723        l >>= 1;
724        for (j = 0; j < m; j++) {
725            ip[m + j] = ip[j] + l;
726        }
727        m <<= 1;
728    }
729    m2 = 2 * m;
730    if ((m << 3) == l) {
731        for (k = 0; k < m; k++) {
732            for (j = 0; j < k; j++) {
733                j1 = 2 * j + ip[k];
734                k1 = 2 * k + ip[j];
735                xr = a[j1];
736                xi = a[j1 + 1];
737                yr = a[k1];
738                yi = a[k1 + 1];
739                a[j1] = yr;
740                a[j1 + 1] = yi;
741                a[k1] = xr;
742                a[k1 + 1] = xi;
743                j1 += m2;
744                k1 += 2 * m2;
745                xr = a[j1];
746                xi = a[j1 + 1];
747                yr = a[k1];
748                yi = a[k1 + 1];
749                a[j1] = yr;
750                a[j1 + 1] = yi;
751                a[k1] = xr;
752                a[k1 + 1] = xi;
753                j1 += m2;
754                k1 -= m2;
755                xr = a[j1];
756                xi = a[j1 + 1];
757                yr = a[k1];
758                yi = a[k1 + 1];
759                a[j1] = yr;
760                a[j1 + 1] = yi;
761                a[k1] = xr;
762                a[k1 + 1] = xi;
763                j1 += m2;
764                k1 += 2 * m2;
765                xr = a[j1];
766                xi = a[j1 + 1];
767                yr = a[k1];
768                yi = a[k1 + 1];
769                a[j1] = yr;
770                a[j1 + 1] = yi;
771                a[k1] = xr;
772                a[k1 + 1] = xi;
773            }
774            j1 = 2 * k + m2 + ip[k];
775            k1 = j1 + m2;
776            xr = a[j1];
777            xi = a[j1 + 1];
778            yr = a[k1];
779            yi = a[k1 + 1];
780            a[j1] = yr;
781            a[j1 + 1] = yi;
782            a[k1] = xr;
783            a[k1 + 1] = xi;
784        }
785    } else {
786        for (k = 1; k < m; k++) {
787            for (j = 0; j < k; j++) {
788                j1 = 2 * j + ip[k];
789                k1 = 2 * k + ip[j];
790                xr = a[j1];
791                xi = a[j1 + 1];
792                yr = a[k1];
793                yi = a[k1 + 1];
794                a[j1] = yr;
795                a[j1 + 1] = yi;
796                a[k1] = xr;
797                a[k1 + 1] = xi;
798                j1 += m2;
799                k1 += m2;
800                xr = a[j1];
801                xi = a[j1 + 1];
802                yr = a[k1];
803                yi = a[k1 + 1];
804                a[j1] = yr;
805                a[j1 + 1] = yi;
806                a[k1] = xr;
807                a[k1 + 1] = xi;
808            }
809        }
810    }
811}
812
813
814void bitrv2conj(int n, int *ip, double *a)
815{
816    int j, j1, k, k1, l, m, m2;
817    double xr, xi, yr, yi;
818   
819    ip[0] = 0;
820    l = n;
821    m = 1;
822    while ((m << 3) < l) {
823        l >>= 1;
824        for (j = 0; j < m; j++) {
825            ip[m + j] = ip[j] + l;
826        }
827        m <<= 1;
828    }
829    m2 = 2 * m;
830    if ((m << 3) == l) {
831        for (k = 0; k < m; k++) {
832            for (j = 0; j < k; j++) {
833                j1 = 2 * j + ip[k];
834                k1 = 2 * k + ip[j];
835                xr = a[j1];
836                xi = -a[j1 + 1];
837                yr = a[k1];
838                yi = -a[k1 + 1];
839                a[j1] = yr;
840                a[j1 + 1] = yi;
841                a[k1] = xr;
842                a[k1 + 1] = xi;
843                j1 += m2;
844                k1 += 2 * m2;
845                xr = a[j1];
846                xi = -a[j1 + 1];
847                yr = a[k1];
848                yi = -a[k1 + 1];
849                a[j1] = yr;
850                a[j1 + 1] = yi;
851                a[k1] = xr;
852                a[k1 + 1] = xi;
853                j1 += m2;
854                k1 -= m2;
855                xr = a[j1];
856                xi = -a[j1 + 1];
857                yr = a[k1];
858                yi = -a[k1 + 1];
859                a[j1] = yr;
860                a[j1 + 1] = yi;
861                a[k1] = xr;
862                a[k1 + 1] = xi;
863                j1 += m2;
864                k1 += 2 * m2;
865                xr = a[j1];
866                xi = -a[j1 + 1];
867                yr = a[k1];
868                yi = -a[k1 + 1];
869                a[j1] = yr;
870                a[j1 + 1] = yi;
871                a[k1] = xr;
872                a[k1 + 1] = xi;
873            }
874            k1 = 2 * k + ip[k];
875            a[k1 + 1] = -a[k1 + 1];
876            j1 = k1 + m2;
877            k1 = j1 + m2;
878            xr = a[j1];
879            xi = -a[j1 + 1];
880            yr = a[k1];
881            yi = -a[k1 + 1];
882            a[j1] = yr;
883            a[j1 + 1] = yi;
884            a[k1] = xr;
885            a[k1 + 1] = xi;
886            k1 += m2;
887            a[k1 + 1] = -a[k1 + 1];
888        }
889    } else {
890        a[1] = -a[1];
891        a[m2 + 1] = -a[m2 + 1];
892        for (k = 1; k < m; k++) {
893            for (j = 0; j < k; j++) {
894                j1 = 2 * j + ip[k];
895                k1 = 2 * k + ip[j];
896                xr = a[j1];
897                xi = -a[j1 + 1];
898                yr = a[k1];
899                yi = -a[k1 + 1];
900                a[j1] = yr;
901                a[j1 + 1] = yi;
902                a[k1] = xr;
903                a[k1 + 1] = xi;
904                j1 += m2;
905                k1 += m2;
906                xr = a[j1];
907                xi = -a[j1 + 1];
908                yr = a[k1];
909                yi = -a[k1 + 1];
910                a[j1] = yr;
911                a[j1 + 1] = yi;
912                a[k1] = xr;
913                a[k1 + 1] = xi;
914            }
915            k1 = 2 * k + ip[k];
916            a[k1 + 1] = -a[k1 + 1];
917            a[k1 + m2 + 1] = -a[k1 + m2 + 1];
918        }
919    }
920}
921
922
923void cftfsub(int n, double *a, double *w)
924{
925    void cft1st(int n, double *a, double *w);
926    void cftmdl(int n, int l, double *a, double *w);
927    int j, j1, j2, j3, l;
928    double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
929   
930    l = 2;
931    if (n >= 16) {
932        cft1st(n, a, w);
933        l = 16;
934        while ((l << 3) <= n) {
935            cftmdl(n, l, a, w);
936            l <<= 3;
937        }
938    }
939    if ((l << 1) < n) {
940        for (j = 0; j < l; j += 2) {
941            j1 = j + l;
942            j2 = j1 + l;
943            j3 = j2 + l;
944            x0r = a[j] + a[j1];
945            x0i = a[j + 1] + a[j1 + 1];
946            x1r = a[j] - a[j1];
947            x1i = a[j + 1] - a[j1 + 1];
948            x2r = a[j2] + a[j3];
949            x2i = a[j2 + 1] + a[j3 + 1];
950            x3r = a[j2] - a[j3];
951            x3i = a[j2 + 1] - a[j3 + 1];
952            a[j] = x0r + x2r;
953            a[j + 1] = x0i + x2i;
954            a[j2] = x0r - x2r;
955            a[j2 + 1] = x0i - x2i;
956            a[j1] = x1r - x3i;
957            a[j1 + 1] = x1i + x3r;
958            a[j3] = x1r + x3i;
959            a[j3 + 1] = x1i - x3r;
960        }
961    } else if ((l << 1) == n) {
962        for (j = 0; j < l; j += 2) {
963            j1 = j + l;
964            x0r = a[j] - a[j1];
965            x0i = a[j + 1] - a[j1 + 1];
966            a[j] += a[j1];
967            a[j + 1] += a[j1 + 1];
968            a[j1] = x0r;
969            a[j1 + 1] = x0i;
970        }
971    }
972}
973
974
975void cftbsub(int n, double *a, double *w)
976{
977    void cft1st(int n, double *a, double *w);
978    void cftmdl(int n, int l, double *a, double *w);
979    int j, j1, j2, j3, j4, j5, j6, j7, l;
980    double wn4r, x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, 
981        y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, 
982        y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i;
983   
984    l = 2;
985    if (n > 16) {
986        cft1st(n, a, w);
987        l = 16;
988        while ((l << 3) < n) {
989            cftmdl(n, l, a, w);
990            l <<= 3;
991        }
992    }
993    if ((l << 2) < n) {
994        wn4r = w[2];
995        for (j = 0; j < l; j += 2) {
996            j1 = j + l;
997            j2 = j1 + l;
998            j3 = j2 + l;
999            j4 = j3 + l;
1000            j5 = j4 + l;
1001            j6 = j5 + l;
1002            j7 = j6 + l;
1003            x0r = a[j] + a[j1];
1004            x0i = -a[j + 1] - a[j1 + 1];
1005            x1r = a[j] - a[j1];
1006            x1i = -a[j + 1] + a[j1 + 1];
1007            x2r = a[j2] + a[j3];
1008            x2i = a[j2 + 1] + a[j3 + 1];
1009            x3r = a[j2] - a[j3];
1010            x3i = a[j2 + 1] - a[j3 + 1];
1011            y0r = x0r + x2r;
1012            y0i = x0i - x2i;
1013            y2r = x0r - x2r;
1014            y2i = x0i + x2i;
1015            y1r = x1r - x3i;
1016            y1i = x1i - x3r;
1017            y3r = x1r + x3i;
1018            y3i = x1i + x3r;
1019            x0r = a[j4] + a[j5];
1020            x0i = a[j4 + 1] + a[j5 + 1];
1021            x1r = a[j4] - a[j5];
1022            x1i = a[j4 + 1] - a[j5 + 1];
1023            x2r = a[j6] + a[j7];
1024            x2i = a[j6 + 1] + a[j7 + 1];
1025            x3r = a[j6] - a[j7];
1026            x3i = a[j6 + 1] - a[j7 + 1];
1027            y4r = x0r + x2r;
1028            y4i = x0i + x2i;
1029            y6r = x0r - x2r;
1030            y6i = x0i - x2i;
1031            x0r = x1r - x3i;
1032            x0i = x1i + x3r;
1033            x2r = x1r + x3i;
1034            x2i = x1i - x3r;
1035            y5r = wn4r * (x0r - x0i);
1036            y5i = wn4r * (x0r + x0i);
1037            y7r = wn4r * (x2r - x2i);
1038            y7i = wn4r * (x2r + x2i);
1039            a[j1] = y1r + y5r;
1040            a[j1 + 1] = y1i - y5i;
1041            a[j5] = y1r - y5r;
1042            a[j5 + 1] = y1i + y5i;
1043            a[j3] = y3r - y7i;
1044            a[j3 + 1] = y3i - y7r;
1045            a[j7] = y3r + y7i;
1046            a[j7 + 1] = y3i + y7r;
1047            a[j] = y0r + y4r;
1048            a[j + 1] = y0i - y4i;
1049            a[j4] = y0r - y4r;
1050            a[j4 + 1] = y0i + y4i;
1051            a[j2] = y2r - y6i;
1052            a[j2 + 1] = y2i - y6r;
1053            a[j6] = y2r + y6i;
1054            a[j6 + 1] = y2i + y6r;
1055        }
1056    } else if ((l << 2) == n) {
1057        for (j = 0; j < l; j += 2) {
1058            j1 = j + l;
1059            j2 = j1 + l;
1060            j3 = j2 + l;
1061            x0r = a[j] + a[j1];
1062            x0i = -a[j + 1] - a[j1 + 1];
1063            x1r = a[j] - a[j1];
1064            x1i = -a[j + 1] + a[j1 + 1];
1065            x2r = a[j2] + a[j3];
1066            x2i = a[j2 + 1] + a[j3 + 1];
1067            x3r = a[j2] - a[j3];
1068            x3i = a[j2 + 1] - a[j3 + 1];
1069            a[j] = x0r + x2r;
1070            a[j + 1] = x0i - x2i;
1071            a[j2] = x0r - x2r;
1072            a[j2 + 1] = x0i + x2i;
1073            a[j1] = x1r - x3i;
1074            a[j1 + 1] = x1i - x3r;
1075            a[j3] = x1r + x3i;
1076            a[j3 + 1] = x1i + x3r;
1077        }
1078    } else {
1079        for (j = 0; j < l; j += 2) {
1080            j1 = j + l;
1081            x0r = a[j] - a[j1];
1082            x0i = -a[j + 1] + a[j1 + 1];
1083            a[j] += a[j1];
1084            a[j + 1] = -a[j + 1] - a[j1 + 1];
1085            a[j1] = x0r;
1086            a[j1 + 1] = x0i;
1087        }
1088    }
1089}
1090
1091
1092void cft1st(int n, double *a, double *w)
1093{
1094    int j, k1;
1095    double wn4r, wtmp, wk1r, wk1i, wk2r, wk2i, wk3r, wk3i, 
1096        wk4r, wk4i, wk5r, wk5i, wk6r, wk6i, wk7r, wk7i;
1097    double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, 
1098        y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, 
1099        y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i;
1100   
1101    wn4r = w[2];
1102    x0r = a[0] + a[2];
1103    x0i = a[1] + a[3];
1104    x1r = a[0] - a[2];
1105    x1i = a[1] - a[3];
1106    x2r = a[4] + a[6];
1107    x2i = a[5] + a[7];
1108    x3r = a[4] - a[6];
1109    x3i = a[5] - a[7];
1110    y0r = x0r + x2r;
1111    y0i = x0i + x2i;
1112    y2r = x0r - x2r;
1113    y2i = x0i - x2i;
1114    y1r = x1r - x3i;
1115    y1i = x1i + x3r;
1116    y3r = x1r + x3i;
1117    y3i = x1i - x3r;
1118    x0r = a[8] + a[10];
1119    x0i = a[9] + a[11];
1120    x1r = a[8] - a[10];
1121    x1i = a[9] - a[11];
1122    x2r = a[12] + a[14];
1123    x2i = a[13] + a[15];
1124    x3r = a[12] - a[14];
1125    x3i = a[13] - a[15];
1126    y4r = x0r + x2r;
1127    y4i = x0i + x2i;
1128    y6r = x0r - x2r;
1129    y6i = x0i - x2i;
1130    x0r = x1r - x3i;
1131    x0i = x1i + x3r;
1132    x2r = x1r + x3i;
1133    x2i = x1i - x3r;
1134    y5r = wn4r * (x0r - x0i);
1135    y5i = wn4r * (x0r + x0i);
1136    y7r = wn4r * (x2r - x2i);
1137    y7i = wn4r * (x2r + x2i);
1138    a[2] = y1r + y5r;
1139    a[3] = y1i + y5i;
1140    a[10] = y1r - y5r;
1141    a[11] = y1i - y5i;
1142    a[6] = y3r - y7i;
1143    a[7] = y3i + y7r;
1144    a[14] = y3r + y7i;
1145    a[15] = y3i - y7r;
1146    a[0] = y0r + y4r;
1147    a[1] = y0i + y4i;
1148    a[8] = y0r - y4r;
1149    a[9] = y0i - y4i;
1150    a[4] = y2r - y6i;
1151    a[5] = y2i + y6r;
1152    a[12] = y2r + y6i;
1153    a[13] = y2i - y6r;
1154    if (n > 16) {
1155        wk1r = w[4];
1156        wk1i = w[5];
1157        x0r = a[16] + a[18];
1158        x0i = a[17] + a[19];
1159        x1r = a[16] - a[18];
1160        x1i = a[17] - a[19];
1161        x2r = a[20] + a[22];
1162        x2i = a[21] + a[23];
1163        x3r = a[20] - a[22];
1164        x3i = a[21] - a[23];
1165        y0r = x0r + x2r;
1166        y0i = x0i + x2i;
1167        y2r = x0r - x2r;
1168        y2i = x0i - x2i;
1169        y1r = x1r - x3i;
1170        y1i = x1i + x3r;
1171        y3r = x1r + x3i;
1172        y3i = x1i - x3r;
1173        x0r = a[24] + a[26];
1174        x0i = a[25] + a[27];
1175        x1r = a[24] - a[26];
1176        x1i = a[25] - a[27];
1177        x2r = a[28] + a[30];
1178        x2i = a[29] + a[31];
1179        x3r = a[28] - a[30];
1180        x3i = a[29] - a[31];
1181        y4r = x0r + x2r;
1182        y4i = x0i + x2i;
1183        y6r = x0r - x2r;
1184        y6i = x0i - x2i;
1185        x0r = x1r - x3i;
1186        x0i = x1i + x3r;
1187        x2r = x1r + x3i;
1188        x2i = x3r - x1i;
1189        y5r = wk1i * x0r - wk1r * x0i;
1190        y5i = wk1i * x0i + wk1r * x0r;
1191        y7r = wk1r * x2r + wk1i * x2i;
1192        y7i = wk1r * x2i - wk1i * x2r;
1193        x0r = wk1r * y1r - wk1i * y1i;
1194        x0i = wk1r * y1i + wk1i * y1r;
1195        a[18] = x0r + y5r;
1196        a[19] = x0i + y5i;
1197        a[26] = y5i - x0i;
1198        a[27] = x0r - y5r;
1199        x0r = wk1i * y3r - wk1r * y3i;
1200        x0i = wk1i * y3i + wk1r * y3r;
1201        a[22] = x0r - y7r;
1202        a[23] = x0i + y7i;
1203        a[30] = y7i - x0i;
1204        a[31] = x0r + y7r;
1205        a[16] = y0r + y4r;
1206        a[17] = y0i + y4i;
1207        a[24] = y4i - y0i;
1208        a[25] = y0r - y4r;
1209        x0r = y2r - y6i;
1210        x0i = y2i + y6r;
1211        a[20] = wn4r * (x0r - x0i);
1212        a[21] = wn4r * (x0i + x0r);
1213        x0r = y6r - y2i;
1214        x0i = y2r + y6i;
1215        a[28] = wn4r * (x0r - x0i);
1216        a[29] = wn4r * (x0i + x0r);
1217        k1 = 4;
1218        for (j = 32; j < n; j += 16) {
1219            k1 += 4;
1220            wk1r = w[k1];
1221            wk1i = w[k1 + 1];
1222            wk2r = w[k1 + 2];
1223            wk2i = w[k1 + 3];
1224            wtmp = 2 * wk2i;
1225            wk3r = wk1r - wtmp * wk1i;
1226            wk3i = wtmp * wk1r - wk1i;
1227            wk4r = 1 - wtmp * wk2i;
1228            wk4i = wtmp * wk2r;
1229            wtmp = 2 * wk4i;
1230            wk5r = wk3r - wtmp * wk1i;
1231            wk5i = wtmp * wk1r - wk3i;
1232            wk6r = wk2r - wtmp * wk2i;
1233            wk6i = wtmp * wk2r - wk2i;
1234            wk7r = wk1r - wtmp * wk3i;
1235            wk7i = wtmp * wk3r - wk1i;
1236            x0r = a[j] + a[j + 2];
1237            x0i = a[j + 1] + a[j + 3];
1238            x1r = a[j] - a[j + 2];
1239            x1i = a[j + 1] - a[j + 3];
1240            x2r = a[j + 4] + a[j + 6];
1241            x2i = a[j + 5] + a[j + 7];
1242            x3r = a[j + 4] - a[j + 6];
1243            x3i = a[j + 5] - a[j + 7];
1244            y0r = x0r + x2r;
1245            y0i = x0i + x2i;
1246            y2r = x0r - x2r;
1247            y2i = x0i - x2i;
1248            y1r = x1r - x3i;
1249            y1i = x1i + x3r;
1250            y3r = x1r + x3i;
1251            y3i = x1i - x3r;
1252            x0r = a[j + 8] + a[j + 10];
1253            x0i = a[j + 9] + a[j + 11];
1254            x1r = a[j + 8] - a[j + 10];
1255            x1i = a[j + 9] - a[j + 11];
1256            x2r = a[j + 12] + a[j + 14];
1257            x2i = a[j + 13] + a[j + 15];
1258            x3r = a[j + 12] - a[j + 14];
1259            x3i = a[j + 13] - a[j + 15];
1260            y4r = x0r + x2r;
1261            y4i = x0i + x2i;
1262            y6r = x0r - x2r;
1263            y6i = x0i - x2i;
1264            x0r = x1r - x3i;
1265            x0i = x1i + x3r;
1266            x2r = x1r + x3i;
1267            x2i = x1i - x3r;
1268            y5r = wn4r * (x0r - x0i);
1269            y5i = wn4r * (x0r + x0i);
1270            y7r = wn4r * (x2r - x2i);
1271            y7i = wn4r * (x2r + x2i);
1272            x0r = y1r + y5r;
1273            x0i = y1i + y5i;
1274            a[j + 2] = wk1r * x0r - wk1i * x0i;
1275            a[j + 3] = wk1r * x0i + wk1i * x0r;
1276            x0r = y1r - y5r;
1277            x0i = y1i - y5i;
1278            a[j + 10] = wk5r * x0r - wk5i * x0i;
1279            a[j + 11] = wk5r * x0i + wk5i * x0r;
1280            x0r = y3r - y7i;
1281            x0i = y3i + y7r;
1282            a[j + 6] = wk3r * x0r - wk3i * x0i;
1283            a[j + 7] = wk3r * x0i + wk3i * x0r;
1284            x0r = y3r + y7i;
1285            x0i = y3i - y7r;
1286            a[j + 14] = wk7r * x0r - wk7i * x0i;
1287            a[j + 15] = wk7r * x0i + wk7i * x0r;
1288            a[j] = y0r + y4r;
1289            a[j + 1] = y0i + y4i;
1290            x0r = y0r - y4r;
1291            x0i = y0i - y4i;
1292            a[j + 8] = wk4r * x0r - wk4i * x0i;
1293            a[j + 9] = wk4r * x0i + wk4i * x0r;
1294            x0r = y2r - y6i;
1295            x0i = y2i + y6r;
1296            a[j + 4] = wk2r * x0r - wk2i * x0i;
1297            a[j + 5] = wk2r * x0i + wk2i * x0r;
1298            x0r = y2r + y6i;
1299            x0i = y2i - y6r;
1300            a[j + 12] = wk6r * x0r - wk6i * x0i;
1301            a[j + 13] = wk6r * x0i + wk6i * x0r;
1302        }
1303    }
1304}
1305
1306
1307void cftmdl(int n, int l, double *a, double *w)
1308{
1309    int j, j1, j2, j3, j4, j5, j6, j7, k, k1, m;
1310    double wn4r, wtmp, wk1r, wk1i, wk2r, wk2i, wk3r, wk3i, 
1311        wk4r, wk4i, wk5r, wk5i, wk6r, wk6i, wk7r, wk7i;
1312    double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, 
1313        y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, 
1314        y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i;
1315   
1316    m = l << 3;
1317    wn4r = w[2];
1318    for (j = 0; j < l; j += 2) {
1319        j1 = j + l;
1320        j2 = j1 + l;
1321        j3 = j2 + l;
1322        j4 = j3 + l;
1323        j5 = j4 + l;
1324        j6 = j5 + l;
1325        j7 = j6 + l;
1326        x0r = a[j] + a[j1];
1327        x0i = a[j + 1] + a[j1 + 1];
1328        x1r = a[j] - a[j1];
1329        x1i = a[j + 1] - a[j1 + 1];
1330        x2r = a[j2] + a[j3];
1331        x2i = a[j2 + 1] + a[j3 + 1];
1332        x3r = a[j2] - a[j3];
1333        x3i = a[j2 + 1] - a[j3 + 1];
1334        y0r = x0r + x2r;
1335        y0i = x0i + x2i;
1336        y2r = x0r - x2r;
1337        y2i = x0i - x2i;
1338        y1r = x1r - x3i;
1339        y1i = x1i + x3r;
1340        y3r = x1r + x3i;
1341        y3i = x1i - x3r;
1342        x0r = a[j4] + a[j5];
1343        x0i = a[j4 + 1] + a[j5 + 1];
1344        x1r = a[j4] - a[j5];
1345        x1i = a[j4 + 1] - a[j5 + 1];
1346        x2r = a[j6] + a[j7];
1347        x2i = a[j6 + 1] + a[j7 + 1];
1348        x3r = a[j6] - a[j7];
1349        x3i = a[j6 + 1] - a[j7 + 1];
1350        y4r = x0r + x2r;
1351        y4i = x0i + x2i;
1352        y6r = x0r - x2r;
1353        y6i = x0i - x2i;
1354        x0r = x1r - x3i;
1355        x0i = x1i + x3r;
1356        x2r = x1r + x3i;
1357        x2i = x1i - x3r;
1358        y5r = wn4r * (x0r - x0i);
1359        y5i = wn4r * (x0r + x0i);
1360        y7r = wn4r * (x2r - x2i);
1361        y7i = wn4r * (x2r + x2i);
1362        a[j1] = y1r + y5r;
1363        a[j1 + 1] = y1i + y5i;
1364        a[j5] = y1r - y5r;
1365        a[j5 + 1] = y1i - y5i;
1366        a[j3] = y3r - y7i;
1367        a[j3 + 1] = y3i + y7r;
1368        a[j7] = y3r + y7i;
1369        a[j7 + 1] = y3i - y7r;
1370        a[j] = y0r + y4r;
1371        a[j + 1] = y0i + y4i;
1372        a[j4] = y0r - y4r;
1373        a[j4 + 1] = y0i - y4i;
1374        a[j2] = y2r - y6i;
1375        a[j2 + 1] = y2i + y6r;
1376        a[j6] = y2r + y6i;
1377        a[j6 + 1] = y2i - y6r;
1378    }
1379    if (m < n) {
1380        wk1r = w[4];
1381        wk1i = w[5];
1382        for (j = m; j < l + m; j += 2) {
1383            j1 = j + l;
1384            j2 = j1 + l;
1385            j3 = j2 + l;
1386            j4 = j3 + l;
1387            j5 = j4 + l;
1388            j6 = j5 + l;
1389            j7 = j6 + l;
1390            x0r = a[j] + a[j1];
1391            x0i = a[j + 1] + a[j1 + 1];
1392            x1r = a[j] - a[j1];
1393            x1i = a[j + 1] - a[j1 + 1];
1394            x2r = a[j2] + a[j3];
1395            x2i = a[j2 + 1] + a[j3 + 1];
1396            x3r = a[j2] - a[j3];
1397            x3i = a[j2 + 1] - a[j3 + 1];
1398            y0r = x0r + x2r;
1399            y0i = x0i + x2i;
1400            y2r = x0r - x2r;
1401            y2i = x0i - x2i;
1402            y1r = x1r - x3i;
1403            y1i = x1i + x3r;
1404            y3r = x1r + x3i;
1405            y3i = x1i - x3r;
1406            x0r = a[j4] + a[j5];
1407            x0i = a[j4 + 1] + a[j5 + 1];
1408            x1r = a[j4] - a[j5];
1409            x1i = a[j4 + 1] - a[j5 + 1];
1410            x2r = a[j6] + a[j7];
1411            x2i = a[j6 + 1] + a[j7 + 1];
1412            x3r = a[j6] - a[j7];
1413            x3i = a[j6 + 1] - a[j7 + 1];
1414            y4r = x0r + x2r;
1415            y4i = x0i + x2i;
1416            y6r = x0r - x2r;
1417            y6i = x0i - x2i;
1418            x0r = x1r - x3i;
1419            x0i = x1i + x3r;
1420            x2r = x1r + x3i;
1421            x2i = x3r - x1i;
1422            y5r = wk1i * x0r - wk1r * x0i;
1423            y5i = wk1i * x0i + wk1r * x0r;
1424            y7r = wk1r * x2r + wk1i * x2i;
1425            y7i = wk1r * x2i - wk1i * x2r;
1426            x0r = wk1r * y1r - wk1i * y1i;
1427            x0i = wk1r * y1i + wk1i * y1r;
1428            a[j1] = x0r + y5r;
1429            a[j1 + 1] = x0i + y5i;
1430            a[j5] = y5i - x0i;
1431            a[j5 + 1] = x0r - y5r;
1432            x0r = wk1i * y3r - wk1r * y3i;
1433            x0i = wk1i * y3i + wk1r * y3r;
1434            a[j3] = x0r - y7r;
1435            a[j3 + 1] = x0i + y7i;
1436            a[j7] = y7i - x0i;
1437            a[j7 + 1] = x0r + y7r;
1438            a[j] = y0r + y4r;
1439            a[j + 1] = y0i + y4i;
1440            a[j4] = y4i - y0i;
1441            a[j4 + 1] = y0r - y4r;
1442            x0r = y2r - y6i;
1443            x0i = y2i + y6r;
1444            a[j2] = wn4r * (x0r - x0i);
1445            a[j2 + 1] = wn4r * (x0i + x0r);
1446            x0r = y6r - y2i;
1447            x0i = y2r + y6i;
1448            a[j6] = wn4r * (x0r - x0i);
1449            a[j6 + 1] = wn4r * (x0i + x0r);
1450        }
1451        k1 = 4;
1452        for (k = 2 * m; k < n; k += m) {
1453            k1 += 4;
1454            wk1r = w[k1];
1455            wk1i = w[k1 + 1];
1456            wk2r = w[k1 + 2];
1457            wk2i = w[k1 + 3];
1458            wtmp = 2 * wk2i;
1459            wk3r = wk1r - wtmp * wk1i;
1460            wk3i = wtmp * wk1r - wk1i;
1461            wk4r = 1 - wtmp * wk2i;
1462            wk4i = wtmp * wk2r;
1463            wtmp = 2 * wk4i;
1464            wk5r = wk3r - wtmp * wk1i;
1465            wk5i = wtmp * wk1r - wk3i;
1466            wk6r = wk2r - wtmp * wk2i;
1467            wk6i = wtmp * wk2r - wk2i;
1468            wk7r = wk1r - wtmp * wk3i;
1469            wk7i = wtmp * wk3r - wk1i;
1470            for (j = k; j < l + k; j += 2) {
1471                j1 = j + l;
1472                j2 = j1 + l;
1473                j3 = j2 + l;
1474                j4 = j3 + l;
1475                j5 = j4 + l;
1476                j6 = j5 + l;
1477                j7 = j6 + l;
1478                x0r = a[j] + a[j1];
1479                x0i = a[j + 1] + a[j1 + 1];
1480                x1r = a[j] - a[j1];
1481                x1i = a[j + 1] - a[j1 + 1];
1482                x2r = a[j2] + a[j3];
1483                x2i = a[j2 + 1] + a[j3 + 1];
1484                x3r = a[j2] - a[j3];
1485                x3i = a[j2 + 1] - a[j3 + 1];
1486                y0r = x0r + x2r;
1487                y0i = x0i + x2i;
1488                y2r = x0r - x2r;
1489                y2i = x0i - x2i;
1490                y1r = x1r - x3i;
1491                y1i = x1i + x3r;
1492                y3r = x1r + x3i;
1493                y3i = x1i - x3r;
1494                x0r = a[j4] + a[j5];
1495                x0i = a[j4 + 1] + a[j5 + 1];
1496                x1r = a[j4] - a[j5];
1497                x1i = a[j4 + 1] - a[j5 + 1];
1498                x2r = a[j6] + a[j7];
1499                x2i = a[j6 + 1] + a[j7 + 1];
1500                x3r = a[j6] - a[j7];
1501                x3i = a[j6 + 1] - a[j7 + 1];
1502                y4r = x0r + x2r;
1503                y4i = x0i + x2i;
1504                y6r = x0r - x2r;
1505                y6i = x0i - x2i;
1506                x0r = x1r - x3i;
1507                x0i = x1i + x3r;
1508                x2r = x1r + x3i;
1509                x2i = x1i - x3r;
1510                y5r = wn4r * (x0r - x0i);
1511                y5i = wn4r * (x0r + x0i);
1512                y7r = wn4r * (x2r - x2i);
1513                y7i = wn4r * (x2r + x2i);
1514                x0r = y1r + y5r;
1515                x0i = y1i + y5i;
1516                a[j1] = wk1r * x0r - wk1i * x0i;
1517                a[j1 + 1] = wk1r * x0i + wk1i * x0r;
1518                x0r = y1r - y5r;
1519                x0i = y1i - y5i;
1520                a[j5] = wk5r * x0r - wk5i * x0i;
1521                a[j5 + 1] = wk5r * x0i + wk5i * x0r;
1522                x0r = y3r - y7i;
1523                x0i = y3i + y7r;
1524                a[j3] = wk3r * x0r - wk3i * x0i;
1525                a[j3 + 1] = wk3r * x0i + wk3i * x0r;
1526                x0r = y3r + y7i;
1527                x0i = y3i - y7r;
1528                a[j7] = wk7r * x0r - wk7i * x0i;
1529                a[j7 + 1] = wk7r * x0i + wk7i * x0r;
1530                a[j] = y0r + y4r;
1531                a[j + 1] = y0i + y4i;
1532                x0r = y0r - y4r;
1533                x0i = y0i - y4i;
1534                a[j4] = wk4r * x0r - wk4i * x0i;
1535                a[j4 + 1] = wk4r * x0i + wk4i * x0r;
1536                x0r = y2r - y6i;
1537                x0i = y2i + y6r;
1538                a[j2] = wk2r * x0r - wk2i * x0i;
1539                a[j2 + 1] = wk2r * x0i + wk2i * x0r;
1540                x0r = y2r + y6i;
1541                x0i = y2i - y6r;
1542                a[j6] = wk6r * x0r - wk6i * x0i;
1543                a[j6 + 1] = wk6r * x0i + wk6i * x0r;
1544            }
1545        }
1546    }
1547}
1548
1549
1550void rftfsub(int n, double *a, int nc, double *c)
1551{
1552    int j, k, kk, ks, m;
1553    double wkr, wki, xr, xi, yr, yi;
1554   
1555    m = n >> 1;
1556    ks = 2 * nc / m;
1557    kk = 0;
1558    for (j = 2; j < m; j += 2) {
1559        k = n - j;
1560        kk += ks;
1561        wkr = 0.5 - c[nc - kk];
1562        wki = c[kk];
1563        xr = a[j] - a[k];
1564        xi = a[j + 1] + a[k + 1];
1565        yr = wkr * xr - wki * xi;
1566        yi = wkr * xi + wki * xr;
1567        a[j] -= yr;
1568        a[j + 1] -= yi;
1569        a[k] += yr;
1570        a[k + 1] -= yi;
1571    }
1572}
1573
1574
1575void rftbsub(int n, double *a, int nc, double *c)
1576{
1577    int j, k, kk, ks, m;
1578    double wkr, wki, xr, xi, yr, yi;
1579   
1580    a[1] = -a[1];
1581    m = n >> 1;
1582    ks = 2 * nc / m;
1583    kk = 0;
1584    for (j = 2; j < m; j += 2) {
1585        k = n - j;
1586        kk += ks;
1587        wkr = 0.5 - c[nc - kk];
1588        wki = c[kk];
1589        xr = a[j] - a[k];
1590        xi = a[j + 1] + a[k + 1];
1591        yr = wkr * xr + wki * xi;
1592        yi = wkr * xi - wki * xr;
1593        a[j] -= yr;
1594        a[j + 1] = yi - a[j + 1];
1595        a[k] += yr;
1596        a[k + 1] = yi - a[k + 1];
1597    }
1598    a[m + 1] = -a[m + 1];
1599}
1600
1601
1602void dctsub(int n, double *a, int nc, double *c)
1603{
1604    int j, k, kk, ks, m;
1605    double wkr, wki, xr;
1606   
1607    m = n >> 1;
1608    ks = nc / n;
1609    kk = 0;
1610    for (j = 1; j < m; j++) {
1611        k = n - j;
1612        kk += ks;
1613        wkr = c[kk] - c[nc - kk];
1614        wki = c[kk] + c[nc - kk];
1615        xr = wki * a[j] - wkr * a[k];
1616        a[j] = wkr * a[j] + wki * a[k];
1617        a[k] = xr;
1618    }
1619    a[m] *= c[0];
1620}
1621
1622
1623void dstsub(int n, double *a, int nc, double *c)
1624{
1625    int j, k, kk, ks, m;
1626    double wkr, wki, xr;
1627   
1628    m = n >> 1;
1629    ks = nc / n;
1630    kk = 0;
1631    for (j = 1; j < m; j++) {
1632        k = n - j;
1633        kk += ks;
1634        wkr = c[kk] - c[nc - kk];
1635        wki = c[kk] + c[nc - kk];
1636        xr = wki * a[k] - wkr * a[j];
1637        a[k] = wkr * a[k] + wki * a[j];
1638        a[j] = xr;
1639    }
1640    a[m] *= c[0];
1641}
1642
Note: See TracBrowser for help on using the repository browser.