3 * Copyright (C) 2000-2003 Michel Lespinasse <walken@zoy.org>
4 * Copyright (C) 1999-2000 Aaron Holtzman <aholtzma@ess.engr.uvic.ca>
6 * The ifft algorithms in this file have been largely inspired by Dan
7 * Bernstein's work, djbfft, available at http://cr.yp.to/djbfft.html
9 * This file is part of a52dec, a free ATSC A-52 stream decoder.
10 * See http://liba52.sourceforge.net/ for updates.
12 * a52dec is free software; you can redistribute it and/or modify
13 * it under the terms of the GNU General Public License as published by
14 * the Free Software Foundation; either version 2 of the License, or
15 * (at your option) any later version.
17 * a52dec is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU General Public License for more details.
22 * You should have received a copy of the GNU General Public License
23 * along with this program; if not, write to the Free Software
24 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
27 #include "config-a52.h"
36 #define M_PI 3.1415926535897932384626433832795029
41 #include "a52_internal.h"
44 typedef struct complex_s
{
49 static const uint8_t fftorder
[] = {
50 0,128, 64,192, 32,160,224, 96, 16,144, 80,208,240,112, 48,176,
51 8,136, 72,200, 40,168,232,104,248,120, 56,184, 24,152,216, 88,
52 4,132, 68,196, 36,164,228,100, 20,148, 84,212,244,116, 52,180,
53 252,124, 60,188, 28,156,220, 92, 12,140, 76,204,236,108, 44,172,
54 2,130, 66,194, 34,162,226, 98, 18,146, 82,210,242,114, 50,178,
55 10,138, 74,202, 42,170,234,106,250,122, 58,186, 26,154,218, 90,
56 254,126, 62,190, 30,158,222, 94, 14,142, 78,206,238,110, 46,174,
57 6,134, 70,198, 38,166,230,102,246,118, 54,182, 22,150,214, 86
60 /* Root values for IFFT */
61 //static sample_t roots16[3];
62 //static sample_t roots32[7];
63 //static sample_t roots64[15];
64 //static sample_t roots128[31];
66 /* Twiddle factors for IMDCT */
67 //static complex_t pre1[128];
68 //static complex_t post1[64];
69 //static complex_t pre2[64];
70 //static complex_t post2[32];
72 //static sample_t a52_imdct_window[256];
73 #include "imdct_lookups.h"
75 static void (* ifft128
) (complex_t
* buf
);
76 static void (* ifft64
) (complex_t
* buf
);
78 static inline void ifft2 (complex_t
* buf
)
84 buf
[0].real
+= buf
[1].real
;
85 buf
[0].imag
+= buf
[1].imag
;
86 buf
[1].real
= r
- buf
[1].real
;
87 buf
[1].imag
= i
- buf
[1].imag
;
90 static inline void ifft4 (complex_t
* buf
)
92 sample_t tmp1
, tmp2
, tmp3
, tmp4
, tmp5
, tmp6
, tmp7
, tmp8
;
94 tmp1
= buf
[0].real
+ buf
[1].real
;
95 tmp2
= buf
[3].real
+ buf
[2].real
;
96 tmp3
= buf
[0].imag
+ buf
[1].imag
;
97 tmp4
= buf
[2].imag
+ buf
[3].imag
;
98 tmp5
= buf
[0].real
- buf
[1].real
;
99 tmp6
= buf
[0].imag
- buf
[1].imag
;
100 tmp7
= buf
[2].imag
- buf
[3].imag
;
101 tmp8
= buf
[3].real
- buf
[2].real
;
103 buf
[0].real
= tmp1
+ tmp2
;
104 buf
[0].imag
= tmp3
+ tmp4
;
105 buf
[2].real
= tmp1
- tmp2
;
106 buf
[2].imag
= tmp3
- tmp4
;
107 buf
[1].real
= tmp5
+ tmp7
;
108 buf
[1].imag
= tmp6
+ tmp8
;
109 buf
[3].real
= tmp5
- tmp7
;
110 buf
[3].imag
= tmp6
- tmp8
;
113 /* basic radix-2 ifft butterfly */
115 #define BUTTERFLY_0(t0,t1,W0,W1,d0,d1) do { \
116 t0 = MUL (W1, d1) + MUL (W0, d0); \
117 t1 = MUL (W0, d1) - MUL (W1, d0); \
120 /* radix-2 ifft butterfly with bias */
122 #define BUTTERFLY_B(t0,t1,W0,W1,d0,d1) do { \
123 t0 = BIAS (MUL (d1, W1) + MUL (d0, W0)); \
124 t1 = BIAS (MUL (d1, W0) - MUL (d0, W1)); \
127 /* the basic split-radix ifft butterfly */
129 #define BUTTERFLY(a0,a1,a2,a3,wr,wi) do { \
130 BUTTERFLY_0 (tmp5, tmp6, wr, wi, a2.real, a2.imag); \
131 BUTTERFLY_0 (tmp8, tmp7, wr, wi, a3.imag, a3.real); \
132 tmp1 = tmp5 + tmp7; \
133 tmp2 = tmp6 + tmp8; \
134 tmp3 = tmp6 - tmp8; \
135 tmp4 = tmp7 - tmp5; \
136 a2.real = a0.real - tmp1; \
137 a2.imag = a0.imag - tmp2; \
138 a3.real = a1.real - tmp3; \
139 a3.imag = a1.imag - tmp4; \
146 /* split-radix ifft butterfly, specialized for wr=1 wi=0 */
148 #define BUTTERFLY_ZERO(a0,a1,a2,a3) do { \
149 tmp1 = a2.real + a3.real; \
150 tmp2 = a2.imag + a3.imag; \
151 tmp3 = a2.imag - a3.imag; \
152 tmp4 = a3.real - a2.real; \
153 a2.real = a0.real - tmp1; \
154 a2.imag = a0.imag - tmp2; \
155 a3.real = a1.real - tmp3; \
156 a3.imag = a1.imag - tmp4; \
163 /* split-radix ifft butterfly, specialized for wr=wi */
165 #define BUTTERFLY_HALF(a0,a1,a2,a3,w) do { \
166 tmp5 = MUL (a2.real + a2.imag, w); \
167 tmp6 = MUL (a2.imag - a2.real, w); \
168 tmp7 = MUL (a3.real - a3.imag, w); \
169 tmp8 = MUL (a3.imag + a3.real, w); \
170 tmp1 = tmp5 + tmp7; \
171 tmp2 = tmp6 + tmp8; \
172 tmp3 = tmp6 - tmp8; \
173 tmp4 = tmp7 - tmp5; \
174 a2.real = a0.real - tmp1; \
175 a2.imag = a0.imag - tmp2; \
176 a3.real = a1.real - tmp3; \
177 a3.imag = a1.imag - tmp4; \
184 static inline void ifft8 (complex_t
* buf
)
186 sample_t tmp1
, tmp2
, tmp3
, tmp4
, tmp5
, tmp6
, tmp7
, tmp8
;
191 BUTTERFLY_ZERO (buf
[0], buf
[2], buf
[4], buf
[6]);
192 BUTTERFLY_HALF (buf
[1], buf
[3], buf
[5], buf
[7], roots16
[1]);
195 static void ifft_pass (complex_t
* buf
, const sample_t
* weight
, int n
)
200 sample_t tmp1
, tmp2
, tmp3
, tmp4
, tmp5
, tmp6
, tmp7
, tmp8
;
208 BUTTERFLY_ZERO (buf
[-1], buf1
[-1], buf2
[-1], buf3
[-1]);
213 BUTTERFLY (buf
[0], buf1
[0], buf2
[0], buf3
[0],
214 weight
[0], weight
[2*i
-n
]);
223 static void ifft16 (complex_t
* buf
)
228 ifft_pass (buf
, roots16
, 4);
231 static void ifft32 (complex_t
* buf
)
236 ifft_pass (buf
, roots32
, 8);
239 static void ifft64_c (complex_t
* buf
)
244 ifft_pass (buf
, roots64
, 16);
247 static void ifft128_c (complex_t
* buf
)
252 ifft_pass (buf
, roots64
, 16);
256 ifft_pass (buf
, roots128
, 32);
259 void a52_imdct_512 (sample_t
* data
, sample_t
* delay
, sample_t bias
)
262 sample_t t_r
, t_i
, a_r
, a_i
, b_r
, b_i
, w_1
, w_2
;
263 const sample_t
* window
= a52_imdct_window
;
266 for (i
= 0; i
< 128; i
++) {
270 BUTTERFLY_0 (buf
[i
].real
, buf
[i
].imag
, t_r
, t_i
, data
[k
], data
[255-k
]);
275 /* Post IFFT complex multiply plus IFFT complex conjugate*/
276 /* Window and convert to real valued signal */
277 for (i
= 0; i
< 64; i
++) {
278 /* y[n] = z[n] * (xcos1[n] + j * xsin1[n]) ; */
281 BUTTERFLY_0 (a_r
, a_i
, t_i
, t_r
, buf
[i
].imag
, buf
[i
].real
);
282 BUTTERFLY_0 (b_r
, b_i
, t_r
, t_i
, buf
[127-i
].imag
, buf
[127-i
].real
);
285 w_2
= window
[255-2*i
];
286 BUTTERFLY_B (data
[255-2*i
], data
[2*i
], w_2
, w_1
, a_r
, delay
[2*i
]);
290 w_2
= window
[254-2*i
];
291 BUTTERFLY_B (data
[2*i
+1], data
[254-2*i
], w_1
, w_2
, b_r
, delay
[2*i
+1]);
296 void a52_imdct_256 (sample_t
* data
, sample_t
* delay
, sample_t bias
)
299 sample_t t_r
, t_i
, a_r
, a_i
, b_r
, b_i
, c_r
, c_i
, d_r
, d_i
, w_1
, w_2
;
300 const sample_t
* window
= a52_imdct_window
;
301 complex_t buf1
[64], buf2
[64];
303 /* Pre IFFT complex multiply plus IFFT cmplx conjugate */
304 for (i
= 0; i
< 64; i
++) {
308 BUTTERFLY_0 (buf1
[i
].real
, buf1
[i
].imag
, t_r
, t_i
, data
[k
], data
[254-k
]);
309 BUTTERFLY_0 (buf2
[i
].real
, buf2
[i
].imag
, t_r
, t_i
, data
[k
+1], data
[255-k
]);
315 /* Post IFFT complex multiply */
316 /* Window and convert to real valued signal */
317 for (i
= 0; i
< 32; i
++) {
318 /* y1[n] = z1[n] * (xcos2[n] + j * xs in2[n]) ; */
321 BUTTERFLY_0 (a_r
, a_i
, t_i
, t_r
, buf1
[i
].imag
, buf1
[i
].real
);
322 BUTTERFLY_0 (b_r
, b_i
, t_r
, t_i
, buf1
[63-i
].imag
, buf1
[63-i
].real
);
323 BUTTERFLY_0 (c_r
, c_i
, t_i
, t_r
, buf2
[i
].imag
, buf2
[i
].real
);
324 BUTTERFLY_0 (d_r
, d_i
, t_r
, t_i
, buf2
[63-i
].imag
, buf2
[63-i
].real
);
327 w_2
= window
[255-2*i
];
328 BUTTERFLY_B (data
[255-2*i
], data
[2*i
], w_2
, w_1
, a_r
, delay
[2*i
]);
331 w_1
= window
[128+2*i
];
332 w_2
= window
[127-2*i
];
333 BUTTERFLY_B (data
[128+2*i
], data
[127-2*i
], w_1
, w_2
, a_i
, delay
[127-2*i
]);
334 delay
[127-2*i
] = c_r
;
337 w_2
= window
[254-2*i
];
338 BUTTERFLY_B (data
[254-2*i
], data
[2*i
+1], w_2
, w_1
, b_i
, delay
[2*i
+1]);
341 w_1
= window
[129+2*i
];
342 w_2
= window
[126-2*i
];
343 BUTTERFLY_B (data
[129+2*i
], data
[126-2*i
], w_1
, w_2
, b_r
, delay
[126-2*i
]);
344 delay
[126-2*i
] = d_i
;
349 static double besselI0 (double x)
355 bessel = bessel * x / (i * i) + 1;
361 void a52_imdct_init (uint32_t mm_accel
)
366 double local_imdct_window[256];*/
368 /* compute imdct window - kaiser-bessel derived window, alpha = 5.0 */
370 for (i = 0; i < 256; i++) {
371 sum += besselI0 (i * (256 - i) * (5 * M_PI / 256) * (5 * M_PI / 256));
372 local_imdct_window[i] = sum;
376 /* for (i = 0; i < 256; i++)
377 a52_imdct_window[i] = SAMPLE (sqrt (local_imdct_window[i] / sum));
379 printf("static sample_t a52_imdct_window[256]={");
380 for (i=0;i<256;i++) {
381 if ((i % 16)==0) { printf("\n"); }
382 printf("%d,",a52_imdct_window[i]);
387 /* for (i = 0; i < 3; i++)
388 roots16[i] = SAMPLE (cos ((M_PI / 8) * (i + 1)));
390 printf("static sample_t roots16[3]={%d,%d,%d};\n\n",roots16[0],roots16[1],roots16[2]);
392 for (i = 0; i < 7; i++)
393 roots32[i] = SAMPLE (cos ((M_PI / 16) * (i + 1)));
395 printf("static sample_t roots32[7]={");
396 for (i=0;i<7;i++) { printf("%d%s",roots32[i],(i < 6 ? "," : "")); }
399 for (i = 0; i < 15; i++)
400 roots64[i] = SAMPLE (cos ((M_PI / 32) * (i + 1)));
402 printf("static sample_t roots64[15]={");
403 for (i=0;i<15;i++) { printf("%d%s",roots64[i],(i < 14 ? "," : "")); }
406 for (i = 0; i < 31; i++)
407 roots128[i] = SAMPLE (cos ((M_PI / 64) * (i + 1)));
409 printf("static sample_t roots128[31]={");
410 for (i=0;i<31;i++) { printf("%d%s",roots128[i],(i < 30 ? "," : "")); }
414 for (i = 0; i < 64; i++) {
415 k = fftorder[i] / 2 + 64;
416 pre1[i].real = SAMPLE (cos ((M_PI / 256) * (k - 0.25)));
417 pre1[i].imag = SAMPLE (sin ((M_PI / 256) * (k - 0.25)));
420 for (i = 64; i < 128; i++) {
421 k = fftorder[i] / 2 + 64;
422 pre1[i].real = SAMPLE (-cos ((M_PI / 256) * (k - 0.25)));
423 pre1[i].imag = SAMPLE (-sin ((M_PI / 256) * (k - 0.25)));
426 printf("static complex_t pre1[128]={");
427 for (i=0;i<128;i++) { printf("{%d,%d}%s",pre1[i].real,pre1[i].imag,(i < 127 ? "," : "")); }
431 for (i = 0; i < 64; i++) {
432 post1[i].real = SAMPLE (cos ((M_PI / 256) * (i + 0.5)));
433 post1[i].imag = SAMPLE (sin ((M_PI / 256) * (i + 0.5)));
436 printf("static complex_t post1[64]={");
437 for (i=0;i<64;i++) { printf("{%d,%d}%s",post1[i].real,post1[i].imag,(i < 63 ? "," : "")); }
442 for (i = 0; i < 64; i++) {
444 pre2[i].real = SAMPLE (cos ((M_PI / 128) * (k - 0.25)));
445 pre2[i].imag = SAMPLE (sin ((M_PI / 128) * (k - 0.25)));
448 printf("static complex_t pre2[64]={");
449 for (i=0;i<64;i++) { printf("{%d,%d}%s",pre2[i].real,pre2[i].imag,(i < 63 ? "," : "")); }
452 for (i = 0; i < 32; i++) {
453 post2[i].real = SAMPLE (cos ((M_PI / 128) * (i + 0.5)));
454 post2[i].imag = SAMPLE (sin ((M_PI / 128) * (i + 0.5)));
457 printf("static complex_t post2[32]={");
458 for (i=0;i<32;i++) { printf("{%d,%d}%s",post2[i].real,post2[i].imag,(i < 31 ? "," : "")); }
463 if (mm_accel
& MM_ACCEL_DJBFFT
) {
464 #ifndef LIBA52_DOUBLE
465 ifft128
= (void (*) (complex_t
*)) fftc4_un128
;
466 ifft64
= (void (*) (complex_t
*)) fftc4_un64
;
468 ifft128
= (void (*) (complex_t
*)) fftc8_un128
;
469 ifft64
= (void (*) (complex_t
*)) fftc8_un64
;