source: src/spectral/ooura_fft8g.c @ 54eba9d

feature/autosinkfeature/constantqfeature/pitchshiftfeature/pydocstringsfeature/timestretchsampleryinfft+
Last change on this file since 54eba9d was 392ad1c, checked in by Paul Brossier <piem@piem.org>, 3 years ago

src/spectral/ooura_fft8g.c: prefix public function with aubio_ooura_ to avoid with other apps using ooura (e.g. puredata), make internal functions static

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