source: src/spectral/ooura_fft8g.c @ e2010b3

feature/crepe_org
Last change on this file since e2010b3 was 8d38841, checked in by Paul Brossier <piem@piem.org>, 8 years ago

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