/[pcsx2_0.9.7]/trunk/3rdparty/liba52/imdct.c
ViewVC logotype

Annotation of /trunk/3rdparty/liba52/imdct.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 10 - (hide annotations) (download)
Mon Sep 6 11:40:06 2010 UTC (9 years, 5 months ago) by william
File MIME type: text/plain
File size: 11818 byte(s)
exported r3113 from ./upstream/trunk
1 william 10 /*
2     * imdct.c
3     * Copyright (C) 2000-2002 Michel Lespinasse <walken@zoy.org>
4     * Copyright (C) 1999-2000 Aaron Holtzman <aholtzma@ess.engr.uvic.ca>
5     *
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
8     *
9     * This file is part of a52dec, a free ATSC A-52 stream decoder.
10     * See http://liba52.sourceforge.net/ for updates.
11     *
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.
16     *
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.
21     *
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
25     */
26     #pragma warning(disable:4244)
27     #include "config.h"
28    
29     #include <math.h>
30     #include <stdio.h>
31     #ifdef LIBA52_DJBFFT
32     #include <fftc4.h>
33     #endif
34     #ifndef M_PI
35     #define M_PI 3.1415926535897932384626433832795029
36     #endif
37     #include "inttypes.h"
38    
39     #include "a52.h"
40     #include "a52_internal.h"
41     #include "mm_accel.h"
42    
43     typedef struct complex_s {
44     sample_t real;
45     sample_t imag;
46     } complex_t;
47    
48     static uint8_t fftorder[] = {
49     0,128, 64,192, 32,160,224, 96, 16,144, 80,208,240,112, 48,176,
50     8,136, 72,200, 40,168,232,104,248,120, 56,184, 24,152,216, 88,
51     4,132, 68,196, 36,164,228,100, 20,148, 84,212,244,116, 52,180,
52     252,124, 60,188, 28,156,220, 92, 12,140, 76,204,236,108, 44,172,
53     2,130, 66,194, 34,162,226, 98, 18,146, 82,210,242,114, 50,178,
54     10,138, 74,202, 42,170,234,106,250,122, 58,186, 26,154,218, 90,
55     254,126, 62,190, 30,158,222, 94, 14,142, 78,206,238,110, 46,174,
56     6,134, 70,198, 38,166,230,102,246,118, 54,182, 22,150,214, 86
57     };
58    
59     /* Root values for IFFT */
60     static sample_t roots16[3];
61     static sample_t roots32[7];
62     static sample_t roots64[15];
63     static sample_t roots128[31];
64    
65     /* Twiddle factors for IMDCT */
66     static complex_t pre1[128];
67     static complex_t post1[64];
68     static complex_t pre2[64];
69     static complex_t post2[32];
70    
71     static sample_t a52_imdct_window[256];
72    
73     static void (* ifft128) (complex_t * buf);
74     static void (* ifft64) (complex_t * buf);
75    
76     static inline void ifft2 (complex_t * buf)
77     {
78     double r, i;
79    
80     r = buf[0].real;
81     i = buf[0].imag;
82     buf[0].real += buf[1].real;
83     buf[0].imag += buf[1].imag;
84     buf[1].real = r - buf[1].real;
85     buf[1].imag = i - buf[1].imag;
86     }
87    
88     static inline void ifft4 (complex_t * buf)
89     {
90     double tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, tmp8;
91    
92     tmp1 = buf[0].real + buf[1].real;
93     tmp2 = buf[3].real + buf[2].real;
94     tmp3 = buf[0].imag + buf[1].imag;
95     tmp4 = buf[2].imag + buf[3].imag;
96     tmp5 = buf[0].real - buf[1].real;
97     tmp6 = buf[0].imag - buf[1].imag;
98     tmp7 = buf[2].imag - buf[3].imag;
99     tmp8 = buf[3].real - buf[2].real;
100    
101     buf[0].real = tmp1 + tmp2;
102     buf[0].imag = tmp3 + tmp4;
103     buf[2].real = tmp1 - tmp2;
104     buf[2].imag = tmp3 - tmp4;
105     buf[1].real = tmp5 + tmp7;
106     buf[1].imag = tmp6 + tmp8;
107     buf[3].real = tmp5 - tmp7;
108     buf[3].imag = tmp6 - tmp8;
109     }
110    
111     /* the basic split-radix ifft butterfly */
112    
113     #define BUTTERFLY(a0,a1,a2,a3,wr,wi) do { \
114     tmp5 = a2.real * wr + a2.imag * wi; \
115     tmp6 = a2.imag * wr - a2.real * wi; \
116     tmp7 = a3.real * wr - a3.imag * wi; \
117     tmp8 = a3.imag * wr + a3.real * wi; \
118     tmp1 = tmp5 + tmp7; \
119     tmp2 = tmp6 + tmp8; \
120     tmp3 = tmp6 - tmp8; \
121     tmp4 = tmp7 - tmp5; \
122     a2.real = a0.real - tmp1; \
123     a2.imag = a0.imag - tmp2; \
124     a3.real = a1.real - tmp3; \
125     a3.imag = a1.imag - tmp4; \
126     a0.real += tmp1; \
127     a0.imag += tmp2; \
128     a1.real += tmp3; \
129     a1.imag += tmp4; \
130     } while (0)
131    
132     /* split-radix ifft butterfly, specialized for wr=1 wi=0 */
133    
134     #define BUTTERFLY_ZERO(a0,a1,a2,a3) do { \
135     tmp1 = a2.real + a3.real; \
136     tmp2 = a2.imag + a3.imag; \
137     tmp3 = a2.imag - a3.imag; \
138     tmp4 = a3.real - a2.real; \
139     a2.real = a0.real - tmp1; \
140     a2.imag = a0.imag - tmp2; \
141     a3.real = a1.real - tmp3; \
142     a3.imag = a1.imag - tmp4; \
143     a0.real += tmp1; \
144     a0.imag += tmp2; \
145     a1.real += tmp3; \
146     a1.imag += tmp4; \
147     } while (0)
148    
149     /* split-radix ifft butterfly, specialized for wr=wi */
150    
151     #define BUTTERFLY_HALF(a0,a1,a2,a3,w) do { \
152     tmp5 = (a2.real + a2.imag) * w; \
153     tmp6 = (a2.imag - a2.real) * w; \
154     tmp7 = (a3.real - a3.imag) * w; \
155     tmp8 = (a3.imag + a3.real) * w; \
156     tmp1 = tmp5 + tmp7; \
157     tmp2 = tmp6 + tmp8; \
158     tmp3 = tmp6 - tmp8; \
159     tmp4 = tmp7 - tmp5; \
160     a2.real = a0.real - tmp1; \
161     a2.imag = a0.imag - tmp2; \
162     a3.real = a1.real - tmp3; \
163     a3.imag = a1.imag - tmp4; \
164     a0.real += tmp1; \
165     a0.imag += tmp2; \
166     a1.real += tmp3; \
167     a1.imag += tmp4; \
168     } while (0)
169    
170     static inline void ifft8 (complex_t * buf)
171     {
172     double tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, tmp8;
173    
174     ifft4 (buf);
175     ifft2 (buf + 4);
176     ifft2 (buf + 6);
177     BUTTERFLY_ZERO (buf[0], buf[2], buf[4], buf[6]);
178     BUTTERFLY_HALF (buf[1], buf[3], buf[5], buf[7], roots16[1]);
179     }
180    
181     static void ifft_pass (complex_t * buf, sample_t * weight, int n)
182     {
183     complex_t * buf1;
184     complex_t * buf2;
185     complex_t * buf3;
186     double tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, tmp8;
187     int i;
188    
189     buf++;
190     buf1 = buf + n;
191     buf2 = buf + 2 * n;
192     buf3 = buf + 3 * n;
193    
194     BUTTERFLY_ZERO (buf[-1], buf1[-1], buf2[-1], buf3[-1]);
195    
196     i = n - 1;
197    
198     do {
199     BUTTERFLY (buf[0], buf1[0], buf2[0], buf3[0], weight[n], weight[2*i]);
200     buf++;
201     buf1++;
202     buf2++;
203     buf3++;
204     weight++;
205     } while (--i);
206     }
207    
208     static void ifft16 (complex_t * buf)
209     {
210     ifft8 (buf);
211     ifft4 (buf + 8);
212     ifft4 (buf + 12);
213     ifft_pass (buf, roots16 - 4, 4);
214     }
215    
216     static void ifft32 (complex_t * buf)
217     {
218     ifft16 (buf);
219     ifft8 (buf + 16);
220     ifft8 (buf + 24);
221     ifft_pass (buf, roots32 - 8, 8);
222     }
223    
224     static void ifft64_c (complex_t * buf)
225     {
226     ifft32 (buf);
227     ifft16 (buf + 32);
228     ifft16 (buf + 48);
229     ifft_pass (buf, roots64 - 16, 16);
230     }
231    
232     static void ifft128_c (complex_t * buf)
233     {
234     ifft32 (buf);
235     ifft16 (buf + 32);
236     ifft16 (buf + 48);
237     ifft_pass (buf, roots64 - 16, 16);
238    
239     ifft32 (buf + 64);
240     ifft32 (buf + 96);
241     ifft_pass (buf, roots128 - 32, 32);
242     }
243    
244     void a52_imdct_512 (sample_t * data, sample_t * delay, sample_t bias)
245     {
246     int i, k;
247     sample_t t_r, t_i, a_r, a_i, b_r, b_i, w_1, w_2;
248     const sample_t * window = a52_imdct_window;
249     complex_t buf[128];
250    
251     for (i = 0; i < 128; i++) {
252     k = fftorder[i];
253     t_r = pre1[i].real;
254     t_i = pre1[i].imag;
255    
256     buf[i].real = t_i * data[255-k] + t_r * data[k];
257     buf[i].imag = t_r * data[255-k] - t_i * data[k];
258     }
259    
260     ifft128 (buf);
261    
262     /* Post IFFT complex multiply plus IFFT complex conjugate*/
263     /* Window and convert to real valued signal */
264     for (i = 0; i < 64; i++) {
265     /* y[n] = z[n] * (xcos1[n] + j * xsin1[n]) ; */
266     t_r = post1[i].real;
267     t_i = post1[i].imag;
268    
269     a_r = t_r * buf[i].real + t_i * buf[i].imag;
270     a_i = t_i * buf[i].real - t_r * buf[i].imag;
271     b_r = t_i * buf[127-i].real + t_r * buf[127-i].imag;
272     b_i = t_r * buf[127-i].real - t_i * buf[127-i].imag;
273    
274     w_1 = window[2*i];
275     w_2 = window[255-2*i];
276     data[2*i] = delay[2*i] * w_2 - a_r * w_1 + bias;
277     data[255-2*i] = delay[2*i] * w_1 + a_r * w_2 + bias;
278     delay[2*i] = a_i;
279    
280     w_1 = window[2*i+1];
281     w_2 = window[254-2*i];
282     data[2*i+1] = delay[2*i+1] * w_2 + b_r * w_1 + bias;
283     data[254-2*i] = delay[2*i+1] * w_1 - b_r * w_2 + bias;
284     delay[2*i+1] = b_i;
285     }
286     }
287    
288     void a52_imdct_256(sample_t * data, sample_t * delay, sample_t bias)
289     {
290     int i, k;
291     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;
292     const sample_t * window = a52_imdct_window;
293     complex_t buf1[64], buf2[64];
294    
295     /* Pre IFFT complex multiply plus IFFT cmplx conjugate */
296     for (i = 0; i < 64; i++) {
297     k = fftorder[i];
298     t_r = pre2[i].real;
299     t_i = pre2[i].imag;
300    
301     buf1[i].real = t_i * data[254-k] + t_r * data[k];
302     buf1[i].imag = t_r * data[254-k] - t_i * data[k];
303    
304     buf2[i].real = t_i * data[255-k] + t_r * data[k+1];
305     buf2[i].imag = t_r * data[255-k] - t_i * data[k+1];
306     }
307    
308     ifft64 (buf1);
309     ifft64 (buf2);
310    
311     /* Post IFFT complex multiply */
312     /* Window and convert to real valued signal */
313     for (i = 0; i < 32; i++) {
314     /* y1[n] = z1[n] * (xcos2[n] + j * xs in2[n]) ; */
315     t_r = post2[i].real;
316     t_i = post2[i].imag;
317    
318     a_r = t_r * buf1[i].real + t_i * buf1[i].imag;
319     a_i = t_i * buf1[i].real - t_r * buf1[i].imag;
320     b_r = t_i * buf1[63-i].real + t_r * buf1[63-i].imag;
321     b_i = t_r * buf1[63-i].real - t_i * buf1[63-i].imag;
322    
323     c_r = t_r * buf2[i].real + t_i * buf2[i].imag;
324     c_i = t_i * buf2[i].real - t_r * buf2[i].imag;
325     d_r = t_i * buf2[63-i].real + t_r * buf2[63-i].imag;
326     d_i = t_r * buf2[63-i].real - t_i * buf2[63-i].imag;
327    
328     w_1 = window[2*i];
329     w_2 = window[255-2*i];
330     data[2*i] = delay[2*i] * w_2 - a_r * w_1 + bias;
331     data[255-2*i] = delay[2*i] * w_1 + a_r * w_2 + bias;
332     delay[2*i] = c_i;
333    
334     w_1 = window[128+2*i];
335     w_2 = window[127-2*i];
336     data[128+2*i] = delay[127-2*i] * w_2 + a_i * w_1 + bias;
337     data[127-2*i] = delay[127-2*i] * w_1 - a_i * w_2 + bias;
338     delay[127-2*i] = c_r;
339    
340     w_1 = window[2*i+1];
341     w_2 = window[254-2*i];
342     data[2*i+1] = delay[2*i+1] * w_2 - b_i * w_1 + bias;
343     data[254-2*i] = delay[2*i+1] * w_1 + b_i * w_2 + bias;
344     delay[2*i+1] = d_r;
345    
346     w_1 = window[129+2*i];
347     w_2 = window[126-2*i];
348     data[129+2*i] = delay[126-2*i] * w_2 + b_r * w_1 + bias;
349     data[126-2*i] = delay[126-2*i] * w_1 - b_r * w_2 + bias;
350     delay[126-2*i] = d_i;
351     }
352     }
353    
354     static double besselI0 (double x)
355     {
356     double bessel = 1;
357     int i = 100;
358    
359     do
360     bessel = bessel * x / (i * i) + 1;
361     while (--i);
362     return bessel;
363     }
364    
365     void a52_imdct_init (uint32_t mm_accel)
366     {
367     int i, k;
368     double sum;
369    
370     /* compute imdct window - kaiser-bessel derived window, alpha = 5.0 */
371     sum = 0;
372     for (i = 0; i < 256; i++) {
373     sum += besselI0 (i * (256 - i) * (5 * M_PI / 256) * (5 * M_PI / 256));
374     a52_imdct_window[i] = sum;
375     }
376     sum++;
377     for (i = 0; i < 256; i++)
378     a52_imdct_window[i] = sqrt (a52_imdct_window[i] / sum);
379    
380     for (i = 0; i < 3; i++)
381     roots16[i] = cos ((M_PI / 8) * (i + 1));
382    
383     for (i = 0; i < 7; i++)
384     roots32[i] = cos ((M_PI / 16) * (i + 1));
385    
386     for (i = 0; i < 15; i++)
387     roots64[i] = cos ((M_PI / 32) * (i + 1));
388    
389     for (i = 0; i < 31; i++)
390     roots128[i] = cos ((M_PI / 64) * (i + 1));
391    
392     for (i = 0; i < 64; i++) {
393     k = fftorder[i] / 2 + 64;
394     pre1[i].real = cos ((M_PI / 256) * (k - 0.25));
395     pre1[i].imag = sin ((M_PI / 256) * (k - 0.25));
396     }
397    
398     for (i = 64; i < 128; i++) {
399     k = fftorder[i] / 2 + 64;
400     pre1[i].real = -cos ((M_PI / 256) * (k - 0.25));
401     pre1[i].imag = -sin ((M_PI / 256) * (k - 0.25));
402     }
403    
404     for (i = 0; i < 64; i++) {
405     post1[i].real = cos ((M_PI / 256) * (i + 0.5));
406     post1[i].imag = sin ((M_PI / 256) * (i + 0.5));
407     }
408    
409     for (i = 0; i < 64; i++) {
410     k = fftorder[i] / 4;
411     pre2[i].real = cos ((M_PI / 128) * (k - 0.25));
412     pre2[i].imag = sin ((M_PI / 128) * (k - 0.25));
413     }
414    
415     for (i = 0; i < 32; i++) {
416     post2[i].real = cos ((M_PI / 128) * (i + 0.5));
417     post2[i].imag = sin ((M_PI / 128) * (i + 0.5));
418     }
419    
420     #ifdef LIBA52_DJBFFT
421     if (mm_accel & MM_ACCEL_DJBFFT) {
422     fprintf (stderr, "Using djbfft for IMDCT transform\n");
423     ifft128 = (void (*) (complex_t *)) fftc4_un128;
424     ifft64 = (void (*) (complex_t *)) fftc4_un64;
425     } else
426     #endif
427     {
428     //fprintf (stderr, "No accelerated IMDCT transform found\n");
429     ifft128 = ifft128_c;
430     ifft64 = ifft64_c;
431     }
432     }

  ViewVC Help
Powered by ViewVC 1.1.22