source: src/spectral/ooura_fft8g.c @ 1e4d90f

feature/autosinkfeature/cnnfeature/cnn_orgfeature/constantqfeature/crepefeature/crepe_orgfeature/pitchshiftfeature/pydocstringsfeature/timestretchfix/ffmpeg5pitchshiftsamplertimestretchyinfft+
Last change on this file since 1e4d90f was 5cb8abe, checked in by Paul Brossier <piem@piem.org>, 11 years ago

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