source: src/spectral/ooura_fft8g.c @ 029bf4e

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

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