source: src/spectral/ooura_fft8g.c @ d389e23

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

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