3 * Copyright (C) 2000-2002 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 * Modified for use with MPlayer, changes contained in liba52_changes.diff.
13 * detailed changelog at http://svn.mplayerhq.hu/mplayer/trunk/
16 * a52dec is free software; you can redistribute it and/or modify
17 * it under the terms of the GNU General Public License as published by
18 * the Free Software Foundation; either version 2 of the License, or
19 * (at your option) any later version.
21 * a52dec is distributed in the hope that it will be useful,
22 * but WITHOUT ANY WARRANTY; without even the implied warranty of
23 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 * GNU General Public License for more details.
26 * You should have received a copy of the GNU General Public License
27 * along with this program; if not, write to the Free Software
28 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
30 * SSE optimizations from Michael Niedermayer (michaelni@gmx.at)
31 * 3DNOW optimizations from Nick Kurshev <nickols_k@mail.ru>
32 * michael did port them from libac3 (untested, perhaps totally broken)
33 * AltiVec optimizations from Romain Dolbeau (romain@dolbeau.org)
44 #define M_PI 3.1415926535897932384626433832795029
49 #include "a52_internal.h"
53 void (*a52_imdct_512
) (sample_t
* data
, sample_t
* delay
, sample_t bias
);
55 #ifdef RUNTIME_CPUDETECT
59 typedef struct complex_s
{
64 static const int pm128
[128] attribute_used
__attribute__((aligned(16))) =
66 0, 16, 32, 48, 64, 80, 96, 112, 8, 40, 72, 104, 24, 56, 88, 120,
67 4, 20, 36, 52, 68, 84, 100, 116, 12, 28, 44, 60, 76, 92, 108, 124,
68 2, 18, 34, 50, 66, 82, 98, 114, 10, 42, 74, 106, 26, 58, 90, 122,
69 6, 22, 38, 54, 70, 86, 102, 118, 14, 46, 78, 110, 30, 62, 94, 126,
70 1, 17, 33, 49, 65, 81, 97, 113, 9, 41, 73, 105, 25, 57, 89, 121,
71 5, 21, 37, 53, 69, 85, 101, 117, 13, 29, 45, 61, 77, 93, 109, 125,
72 3, 19, 35, 51, 67, 83, 99, 115, 11, 43, 75, 107, 27, 59, 91, 123,
73 7, 23, 39, 55, 71, 87, 103, 119, 15, 31, 47, 63, 79, 95, 111, 127
76 static uint8_t attribute_used bit_reverse_512
[] = {
77 0x00, 0x40, 0x20, 0x60, 0x10, 0x50, 0x30, 0x70,
78 0x08, 0x48, 0x28, 0x68, 0x18, 0x58, 0x38, 0x78,
79 0x04, 0x44, 0x24, 0x64, 0x14, 0x54, 0x34, 0x74,
80 0x0c, 0x4c, 0x2c, 0x6c, 0x1c, 0x5c, 0x3c, 0x7c,
81 0x02, 0x42, 0x22, 0x62, 0x12, 0x52, 0x32, 0x72,
82 0x0a, 0x4a, 0x2a, 0x6a, 0x1a, 0x5a, 0x3a, 0x7a,
83 0x06, 0x46, 0x26, 0x66, 0x16, 0x56, 0x36, 0x76,
84 0x0e, 0x4e, 0x2e, 0x6e, 0x1e, 0x5e, 0x3e, 0x7e,
85 0x01, 0x41, 0x21, 0x61, 0x11, 0x51, 0x31, 0x71,
86 0x09, 0x49, 0x29, 0x69, 0x19, 0x59, 0x39, 0x79,
87 0x05, 0x45, 0x25, 0x65, 0x15, 0x55, 0x35, 0x75,
88 0x0d, 0x4d, 0x2d, 0x6d, 0x1d, 0x5d, 0x3d, 0x7d,
89 0x03, 0x43, 0x23, 0x63, 0x13, 0x53, 0x33, 0x73,
90 0x0b, 0x4b, 0x2b, 0x6b, 0x1b, 0x5b, 0x3b, 0x7b,
91 0x07, 0x47, 0x27, 0x67, 0x17, 0x57, 0x37, 0x77,
92 0x0f, 0x4f, 0x2f, 0x6f, 0x1f, 0x5f, 0x3f, 0x7f};
94 static uint8_t fftorder
[] = {
95 0,128, 64,192, 32,160,224, 96, 16,144, 80,208,240,112, 48,176,
96 8,136, 72,200, 40,168,232,104,248,120, 56,184, 24,152,216, 88,
97 4,132, 68,196, 36,164,228,100, 20,148, 84,212,244,116, 52,180,
98 252,124, 60,188, 28,156,220, 92, 12,140, 76,204,236,108, 44,172,
99 2,130, 66,194, 34,162,226, 98, 18,146, 82,210,242,114, 50,178,
100 10,138, 74,202, 42,170,234,106,250,122, 58,186, 26,154,218, 90,
101 254,126, 62,190, 30,158,222, 94, 14,142, 78,206,238,110, 46,174,
102 6,134, 70,198, 38,166,230,102,246,118, 54,182, 22,150,214, 86
105 static complex_t
__attribute__((aligned(16))) buf
[128];
107 /* Twiddle factor LUT */
108 static complex_t
__attribute__((aligned(16))) w_1
[1];
109 static complex_t
__attribute__((aligned(16))) w_2
[2];
110 static complex_t
__attribute__((aligned(16))) w_4
[4];
111 static complex_t
__attribute__((aligned(16))) w_8
[8];
112 static complex_t
__attribute__((aligned(16))) w_16
[16];
113 static complex_t
__attribute__((aligned(16))) w_32
[32];
114 static complex_t
__attribute__((aligned(16))) w_64
[64];
115 static complex_t
__attribute__((aligned(16))) * w
[7] = {w_1
, w_2
, w_4
, w_8
, w_16
, w_32
, w_64
};
117 /* Twiddle factors for IMDCT */
118 static sample_t
__attribute__((aligned(16))) xcos1
[128];
119 static sample_t
__attribute__((aligned(16))) xsin1
[128];
121 #if defined(ARCH_X86) || defined(ARCH_X86_64)
122 // NOTE: SSE needs 16byte alignment or it will segfault
124 static float __attribute__((aligned(16))) sseSinCos1c
[256];
125 static float __attribute__((aligned(16))) sseSinCos1d
[256];
126 static float attribute_used
__attribute__((aligned(16))) ps111_1
[4]={1,1,1,-1};
127 //static float __attribute__((aligned(16))) sseW0[4];
128 static float __attribute__((aligned(16))) sseW1
[8];
129 static float __attribute__((aligned(16))) sseW2
[16];
130 static float __attribute__((aligned(16))) sseW3
[32];
131 static float __attribute__((aligned(16))) sseW4
[64];
132 static float __attribute__((aligned(16))) sseW5
[128];
133 static float __attribute__((aligned(16))) sseW6
[256];
134 static float __attribute__((aligned(16))) *sseW
[7]=
135 {NULL
/*sseW0*/,sseW1
,sseW2
,sseW3
,sseW4
,sseW5
,sseW6
};
136 static float __attribute__((aligned(16))) sseWindow
[512];
139 /* Root values for IFFT */
140 static sample_t roots16
[3];
141 static sample_t roots32
[7];
142 static sample_t roots64
[15];
143 static sample_t roots128
[31];
145 /* Twiddle factors for IMDCT */
146 static complex_t pre1
[128];
147 static complex_t post1
[64];
148 static complex_t pre2
[64];
149 static complex_t post2
[32];
151 static sample_t a52_imdct_window
[256];
153 static void (* ifft128
) (complex_t
* buf
);
154 static void (* ifft64
) (complex_t
* buf
);
156 static inline void ifft2 (complex_t
* buf
)
162 buf
[0].real
+= buf
[1].real
;
163 buf
[0].imag
+= buf
[1].imag
;
164 buf
[1].real
= r
- buf
[1].real
;
165 buf
[1].imag
= i
- buf
[1].imag
;
168 static inline void ifft4 (complex_t
* buf
)
170 double tmp1
, tmp2
, tmp3
, tmp4
, tmp5
, tmp6
, tmp7
, tmp8
;
172 tmp1
= buf
[0].real
+ buf
[1].real
;
173 tmp2
= buf
[3].real
+ buf
[2].real
;
174 tmp3
= buf
[0].imag
+ buf
[1].imag
;
175 tmp4
= buf
[2].imag
+ buf
[3].imag
;
176 tmp5
= buf
[0].real
- buf
[1].real
;
177 tmp6
= buf
[0].imag
- buf
[1].imag
;
178 tmp7
= buf
[2].imag
- buf
[3].imag
;
179 tmp8
= buf
[3].real
- buf
[2].real
;
181 buf
[0].real
= tmp1
+ tmp2
;
182 buf
[0].imag
= tmp3
+ tmp4
;
183 buf
[2].real
= tmp1
- tmp2
;
184 buf
[2].imag
= tmp3
- tmp4
;
185 buf
[1].real
= tmp5
+ tmp7
;
186 buf
[1].imag
= tmp6
+ tmp8
;
187 buf
[3].real
= tmp5
- tmp7
;
188 buf
[3].imag
= tmp6
- tmp8
;
191 /* the basic split-radix ifft butterfly */
193 #define BUTTERFLY(a0,a1,a2,a3,wr,wi) do { \
194 tmp5 = a2.real * wr + a2.imag * wi; \
195 tmp6 = a2.imag * wr - a2.real * wi; \
196 tmp7 = a3.real * wr - a3.imag * wi; \
197 tmp8 = a3.imag * wr + a3.real * wi; \
198 tmp1 = tmp5 + tmp7; \
199 tmp2 = tmp6 + tmp8; \
200 tmp3 = tmp6 - tmp8; \
201 tmp4 = tmp7 - tmp5; \
202 a2.real = a0.real - tmp1; \
203 a2.imag = a0.imag - tmp2; \
204 a3.real = a1.real - tmp3; \
205 a3.imag = a1.imag - tmp4; \
212 /* split-radix ifft butterfly, specialized for wr=1 wi=0 */
214 #define BUTTERFLY_ZERO(a0,a1,a2,a3) do { \
215 tmp1 = a2.real + a3.real; \
216 tmp2 = a2.imag + a3.imag; \
217 tmp3 = a2.imag - a3.imag; \
218 tmp4 = a3.real - a2.real; \
219 a2.real = a0.real - tmp1; \
220 a2.imag = a0.imag - tmp2; \
221 a3.real = a1.real - tmp3; \
222 a3.imag = a1.imag - tmp4; \
229 /* split-radix ifft butterfly, specialized for wr=wi */
231 #define BUTTERFLY_HALF(a0,a1,a2,a3,w) do { \
232 tmp5 = (a2.real + a2.imag) * w; \
233 tmp6 = (a2.imag - a2.real) * w; \
234 tmp7 = (a3.real - a3.imag) * w; \
235 tmp8 = (a3.imag + a3.real) * w; \
236 tmp1 = tmp5 + tmp7; \
237 tmp2 = tmp6 + tmp8; \
238 tmp3 = tmp6 - tmp8; \
239 tmp4 = tmp7 - tmp5; \
240 a2.real = a0.real - tmp1; \
241 a2.imag = a0.imag - tmp2; \
242 a3.real = a1.real - tmp3; \
243 a3.imag = a1.imag - tmp4; \
250 static inline void ifft8 (complex_t
* buf
)
252 double tmp1
, tmp2
, tmp3
, tmp4
, tmp5
, tmp6
, tmp7
, tmp8
;
257 BUTTERFLY_ZERO (buf
[0], buf
[2], buf
[4], buf
[6]);
258 BUTTERFLY_HALF (buf
[1], buf
[3], buf
[5], buf
[7], roots16
[1]);
261 static void ifft_pass (complex_t
* buf
, sample_t
* weight
, int n
)
266 double tmp1
, tmp2
, tmp3
, tmp4
, tmp5
, tmp6
, tmp7
, tmp8
;
274 BUTTERFLY_ZERO (buf
[-1], buf1
[-1], buf2
[-1], buf3
[-1]);
279 BUTTERFLY (buf
[0], buf1
[0], buf2
[0], buf3
[0], weight
[n
], weight
[2*i
]);
288 static void ifft16 (complex_t
* buf
)
293 ifft_pass (buf
, roots16
- 4, 4);
296 static void ifft32 (complex_t
* buf
)
301 ifft_pass (buf
, roots32
- 8, 8);
304 static void ifft64_c (complex_t
* buf
)
309 ifft_pass (buf
, roots64
- 16, 16);
312 static void ifft128_c (complex_t
* buf
)
317 ifft_pass (buf
, roots64
- 16, 16);
321 ifft_pass (buf
, roots128
- 32, 32);
324 void imdct_do_512 (sample_t
* data
, sample_t
* delay
, sample_t bias
)
327 sample_t t_r
, t_i
, a_r
, a_i
, b_r
, b_i
, w_1
, w_2
;
328 const sample_t
* window
= a52_imdct_window
;
331 for (i
= 0; i
< 128; i
++) {
336 buf
[i
].real
= t_i
* data
[255-k
] + t_r
* data
[k
];
337 buf
[i
].imag
= t_r
* data
[255-k
] - t_i
* data
[k
];
342 /* Post IFFT complex multiply plus IFFT complex conjugate*/
343 /* Window and convert to real valued signal */
344 for (i
= 0; i
< 64; i
++) {
345 /* y[n] = z[n] * (xcos1[n] + j * xsin1[n]) ; */
349 a_r
= t_r
* buf
[i
].real
+ t_i
* buf
[i
].imag
;
350 a_i
= t_i
* buf
[i
].real
- t_r
* buf
[i
].imag
;
351 b_r
= t_i
* buf
[127-i
].real
+ t_r
* buf
[127-i
].imag
;
352 b_i
= t_r
* buf
[127-i
].real
- t_i
* buf
[127-i
].imag
;
355 w_2
= window
[255-2*i
];
356 data
[2*i
] = delay
[2*i
] * w_2
- a_r
* w_1
+ bias
;
357 data
[255-2*i
] = delay
[2*i
] * w_1
+ a_r
* w_2
+ bias
;
361 w_2
= window
[254-2*i
];
362 data
[2*i
+1] = delay
[2*i
+1] * w_2
+ b_r
* w_1
+ bias
;
363 data
[254-2*i
] = delay
[2*i
+1] * w_1
- b_r
* w_2
+ bias
;
370 #ifdef HAVE_ALTIVEC_H
374 // used to build registers permutation vectors (vcprm)
375 // the 's' are for words in the _s_econd vector
376 #define WORD_0 0x00,0x01,0x02,0x03
377 #define WORD_1 0x04,0x05,0x06,0x07
378 #define WORD_2 0x08,0x09,0x0a,0x0b
379 #define WORD_3 0x0c,0x0d,0x0e,0x0f
380 #define WORD_s0 0x10,0x11,0x12,0x13
381 #define WORD_s1 0x14,0x15,0x16,0x17
382 #define WORD_s2 0x18,0x19,0x1a,0x1b
383 #define WORD_s3 0x1c,0x1d,0x1e,0x1f
385 #define vcprm(a,b,c,d) (const vector unsigned char){WORD_ ## a, WORD_ ## b, WORD_ ## c, WORD_ ## d}
386 #define vcii(a,b,c,d) (const vector float){FLOAT_ ## a, FLOAT_ ## b, FLOAT_ ## c, FLOAT_ ## d}
388 #define FOUROF(a) {a,a,a,a}
390 // vcprmle is used to keep the same index as in the SSE version.
391 // it's the same as vcprm, with the index inversed
392 // ('le' is Little Endian)
393 #define vcprmle(a,b,c,d) vcprm(d,c,b,a)
395 // used to build inverse/identity vectors (vcii)
396 // n is _n_egative, p is _p_ositive
402 imdct_do_512_altivec(sample_t data
[],sample_t delay
[], sample_t bias
)
418 sample_t
*window_ptr
;
420 /* 512 IMDCT with source and dest data in 'data' */
422 /* Pre IFFT complex multiply plus IFFT cmplx conjugate & reordering*/
423 for( i
=0; i
< 128; i
++) {
424 /* z[i] = (X[256-2*i-1] + j * X[2*i]) * (xcos1[i] + j * xsin1[i]) ; */
425 int j
= bit_reverse_512
[i
];
426 buf
[i
].real
= (data
[256-2*j
-1] * xcos1
[j
]) - (data
[2*j
] * xsin1
[j
]);
427 buf
[i
].imag
= -1.0 * ((data
[2*j
] * xcos1
[j
]) + (data
[256-2*j
-1] * xsin1
[j
]));
431 for(i
= 0; i
< 128; i
+= 2) {
433 tmp_a_r
= buf
[i
].real
;
434 tmp_a_i
= buf
[i
].imag
;
435 tmp_b_r
= buf
[i
+1].real
;
436 tmp_b_i
= buf
[i
+1].imag
;
437 buf
[i
].real
= tmp_a_r
+ tmp_b_r
;
438 buf
[i
].imag
= tmp_a_i
+ tmp_b_i
;
439 buf
[i
+1].real
= tmp_a_r
- tmp_b_r
;
440 buf
[i
+1].imag
= tmp_a_i
- tmp_b_i
;
442 vector
float temp
, bufv
;
444 bufv
= vec_ld(i
<< 3, (float*)buf
);
445 temp
= vec_perm(bufv
, bufv
, vcprm(2,3,0,1));
446 bufv
= vec_madd(bufv
, vcii(p
,p
,n
,n
), temp
);
447 vec_st(bufv
, i
<< 3, (float*)buf
);
452 // Note w[1]={{1,0}, {0,-1}}
453 for(i
= 0; i
< 128; i
+= 4) {
455 tmp_a_r
= buf
[i
].real
;
456 tmp_a_i
= buf
[i
].imag
;
457 tmp_b_r
= buf
[i
+2].real
;
458 tmp_b_i
= buf
[i
+2].imag
;
459 buf
[i
].real
= tmp_a_r
+ tmp_b_r
;
460 buf
[i
].imag
= tmp_a_i
+ tmp_b_i
;
461 buf
[i
+2].real
= tmp_a_r
- tmp_b_r
;
462 buf
[i
+2].imag
= tmp_a_i
- tmp_b_i
;
463 tmp_a_r
= buf
[i
+1].real
;
464 tmp_a_i
= buf
[i
+1].imag
;
465 /* WARNING: im <-> re here ! */
466 tmp_b_r
= buf
[i
+3].imag
;
467 tmp_b_i
= buf
[i
+3].real
;
468 buf
[i
+1].real
= tmp_a_r
+ tmp_b_r
;
469 buf
[i
+1].imag
= tmp_a_i
- tmp_b_i
;
470 buf
[i
+3].real
= tmp_a_r
- tmp_b_r
;
471 buf
[i
+3].imag
= tmp_a_i
+ tmp_b_i
;
473 vector
float buf01
, buf23
, temp1
, temp2
;
475 buf01
= vec_ld((i
+ 0) << 3, (float*)buf
);
476 buf23
= vec_ld((i
+ 2) << 3, (float*)buf
);
477 buf23
= vec_perm(buf23
,buf23
,vcprm(0,1,3,2));
479 temp1
= vec_madd(buf23
, vcii(p
,p
,p
,n
), buf01
);
480 temp2
= vec_madd(buf23
, vcii(n
,n
,n
,p
), buf01
);
482 vec_st(temp1
, (i
+ 0) << 3, (float*)buf
);
483 vec_st(temp2
, (i
+ 2) << 3, (float*)buf
);
488 for(i
= 0; i
< 128; i
+= 8) {
490 tmp_a_r
= buf
[i
].real
;
491 tmp_a_i
= buf
[i
].imag
;
492 tmp_b_r
= buf
[i
+4].real
;
493 tmp_b_i
= buf
[i
+4].imag
;
494 buf
[i
].real
= tmp_a_r
+ tmp_b_r
;
495 buf
[i
].imag
= tmp_a_i
+ tmp_b_i
;
496 buf
[i
+4].real
= tmp_a_r
- tmp_b_r
;
497 buf
[i
+4].imag
= tmp_a_i
- tmp_b_i
;
498 tmp_a_r
= buf
[1+i
].real
;
499 tmp_a_i
= buf
[1+i
].imag
;
500 tmp_b_r
= (buf
[i
+5].real
+ buf
[i
+5].imag
) * w
[2][1].real
;
501 tmp_b_i
= (buf
[i
+5].imag
- buf
[i
+5].real
) * w
[2][1].real
;
502 buf
[1+i
].real
= tmp_a_r
+ tmp_b_r
;
503 buf
[1+i
].imag
= tmp_a_i
+ tmp_b_i
;
504 buf
[i
+5].real
= tmp_a_r
- tmp_b_r
;
505 buf
[i
+5].imag
= tmp_a_i
- tmp_b_i
;
506 tmp_a_r
= buf
[i
+2].real
;
507 tmp_a_i
= buf
[i
+2].imag
;
508 /* WARNING re <-> im & sign */
509 tmp_b_r
= buf
[i
+6].imag
;
510 tmp_b_i
= - buf
[i
+6].real
;
511 buf
[i
+2].real
= tmp_a_r
+ tmp_b_r
;
512 buf
[i
+2].imag
= tmp_a_i
+ tmp_b_i
;
513 buf
[i
+6].real
= tmp_a_r
- tmp_b_r
;
514 buf
[i
+6].imag
= tmp_a_i
- tmp_b_i
;
515 tmp_a_r
= buf
[i
+3].real
;
516 tmp_a_i
= buf
[i
+3].imag
;
517 tmp_b_r
= (buf
[i
+7].real
- buf
[i
+7].imag
) * w
[2][3].imag
;
518 tmp_b_i
= (buf
[i
+7].imag
+ buf
[i
+7].real
) * w
[2][3].imag
;
519 buf
[i
+3].real
= tmp_a_r
+ tmp_b_r
;
520 buf
[i
+3].imag
= tmp_a_i
+ tmp_b_i
;
521 buf
[i
+7].real
= tmp_a_r
- tmp_b_r
;
522 buf
[i
+7].imag
= tmp_a_i
- tmp_b_i
;
524 vector
float buf01
, buf23
, buf45
, buf67
;
526 buf01
= vec_ld((i
+ 0) << 3, (float*)buf
);
527 buf23
= vec_ld((i
+ 2) << 3, (float*)buf
);
529 tmp_b_r
= (buf
[i
+5].real
+ buf
[i
+5].imag
) * w
[2][1].real
;
530 tmp_b_i
= (buf
[i
+5].imag
- buf
[i
+5].real
) * w
[2][1].real
;
531 buf
[i
+5].real
= tmp_b_r
;
532 buf
[i
+5].imag
= tmp_b_i
;
533 tmp_b_r
= (buf
[i
+7].real
- buf
[i
+7].imag
) * w
[2][3].imag
;
534 tmp_b_i
= (buf
[i
+7].imag
+ buf
[i
+7].real
) * w
[2][3].imag
;
535 buf
[i
+7].real
= tmp_b_r
;
536 buf
[i
+7].imag
= tmp_b_i
;
538 buf23
= vec_ld((i
+ 2) << 3, (float*)buf
);
539 buf45
= vec_ld((i
+ 4) << 3, (float*)buf
);
540 buf67
= vec_ld((i
+ 6) << 3, (float*)buf
);
541 buf67
= vec_perm(buf67
, buf67
, vcprm(1,0,2,3));
543 vec_st(vec_add(buf01
, buf45
), (i
+ 0) << 3, (float*)buf
);
544 vec_st(vec_madd(buf67
, vcii(p
,n
,p
,p
), buf23
), (i
+ 2) << 3, (float*)buf
);
545 vec_st(vec_sub(buf01
, buf45
), (i
+ 4) << 3, (float*)buf
);
546 vec_st(vec_nmsub(buf67
, vcii(p
,n
,p
,p
), buf23
), (i
+ 6) << 3, (float*)buf
);
550 /* 4-7. iterations */
551 for (m
=3; m
< 7; m
++) {
554 two_m_plus_one
= two_m
<<1;
556 for(i
= 0; i
< 128; i
+= two_m_plus_one
) {
557 for(k
= 0; k
< two_m
; k
+=2) {
561 tmp_a_r
= buf
[p
].real
;
562 tmp_a_i
= buf
[p
].imag
;
564 buf
[q
].real
* w
[m
][k
].real
-
565 buf
[q
].imag
* w
[m
][k
].imag
;
567 buf
[q
].imag
* w
[m
][k
].real
+
568 buf
[q
].real
* w
[m
][k
].imag
;
569 buf
[p
].real
= tmp_a_r
+ tmp_b_r
;
570 buf
[p
].imag
= tmp_a_i
+ tmp_b_i
;
571 buf
[q
].real
= tmp_a_r
- tmp_b_r
;
572 buf
[q
].imag
= tmp_a_i
- tmp_b_i
;
574 tmp_a_r
= buf
[(p
+ 1)].real
;
575 tmp_a_i
= buf
[(p
+ 1)].imag
;
577 buf
[(q
+ 1)].real
* w
[m
][(k
+ 1)].real
-
578 buf
[(q
+ 1)].imag
* w
[m
][(k
+ 1)].imag
;
580 buf
[(q
+ 1)].imag
* w
[m
][(k
+ 1)].real
+
581 buf
[(q
+ 1)].real
* w
[m
][(k
+ 1)].imag
;
582 buf
[(p
+ 1)].real
= tmp_a_r
+ tmp_b_r
;
583 buf
[(p
+ 1)].imag
= tmp_a_i
+ tmp_b_i
;
584 buf
[(q
+ 1)].real
= tmp_a_r
- tmp_b_r
;
585 buf
[(q
+ 1)].imag
= tmp_a_i
- tmp_b_i
;
589 vector
float vecp
, vecq
, vecw
, temp1
, temp2
, temp3
, temp4
;
590 const vector
float vczero
= (const vector
float)FOUROF(0.);
591 // first compute buf[q] and buf[q+1]
592 vecq
= vec_ld(q
<< 3, (float*)buf
);
593 vecw
= vec_ld(0, (float*)&(w
[m
][k
]));
594 temp1
= vec_madd(vecq
, vecw
, vczero
);
595 temp2
= vec_perm(vecq
, vecq
, vcprm(1,0,3,2));
596 temp2
= vec_madd(temp2
, vecw
, vczero
);
597 temp3
= vec_perm(temp1
, temp2
, vcprm(0,s0
,2,s2
));
598 temp4
= vec_perm(temp1
, temp2
, vcprm(1,s1
,3,s3
));
599 vecq
= vec_madd(temp4
, vcii(n
,p
,n
,p
), temp3
);
600 // then butterfly with buf[p] and buf[p+1]
601 vecp
= vec_ld(p
<< 3, (float*)buf
);
603 temp1
= vec_add(vecp
, vecq
);
604 temp2
= vec_sub(vecp
, vecq
);
606 vec_st(temp1
, p
<< 3, (float*)buf
);
607 vec_st(temp2
, q
<< 3, (float*)buf
);
613 /* Post IFFT complex multiply plus IFFT complex conjugate*/
614 for( i
=0; i
< 128; i
+=4) {
615 /* y[n] = z[n] * (xcos1[n] + j * xsin1[n]) ; */
617 tmp_a_r
= buf
[(i
+ 0)].real
;
618 tmp_a_i
= -1.0 * buf
[(i
+ 0)].imag
;
620 (tmp_a_r
* xcos1
[(i
+ 0)]) - (tmp_a_i
* xsin1
[(i
+ 0)]);
622 (tmp_a_r
* xsin1
[(i
+ 0)]) + (tmp_a_i
* xcos1
[(i
+ 0)]);
624 tmp_a_r
= buf
[(i
+ 1)].real
;
625 tmp_a_i
= -1.0 * buf
[(i
+ 1)].imag
;
627 (tmp_a_r
* xcos1
[(i
+ 1)]) - (tmp_a_i
* xsin1
[(i
+ 1)]);
629 (tmp_a_r
* xsin1
[(i
+ 1)]) + (tmp_a_i
* xcos1
[(i
+ 1)]);
631 tmp_a_r
= buf
[(i
+ 2)].real
;
632 tmp_a_i
= -1.0 * buf
[(i
+ 2)].imag
;
634 (tmp_a_r
* xcos1
[(i
+ 2)]) - (tmp_a_i
* xsin1
[(i
+ 2)]);
636 (tmp_a_r
* xsin1
[(i
+ 2)]) + (tmp_a_i
* xcos1
[(i
+ 2)]);
638 tmp_a_r
= buf
[(i
+ 3)].real
;
639 tmp_a_i
= -1.0 * buf
[(i
+ 3)].imag
;
641 (tmp_a_r
* xcos1
[(i
+ 3)]) - (tmp_a_i
* xsin1
[(i
+ 3)]);
643 (tmp_a_r
* xsin1
[(i
+ 3)]) + (tmp_a_i
* xcos1
[(i
+ 3)]);
645 vector
float bufv_0
, bufv_2
, cosv
, sinv
, temp1
, temp2
;
646 vector
float temp0022
, temp1133
, tempCS01
;
647 const vector
float vczero
= (const vector
float)FOUROF(0.);
649 bufv_0
= vec_ld((i
+ 0) << 3, (float*)buf
);
650 bufv_2
= vec_ld((i
+ 2) << 3, (float*)buf
);
652 cosv
= vec_ld(i
<< 2, xcos1
);
653 sinv
= vec_ld(i
<< 2, xsin1
);
655 temp0022
= vec_perm(bufv_0
, bufv_0
, vcprm(0,0,2,2));
656 temp1133
= vec_perm(bufv_0
, bufv_0
, vcprm(1,1,3,3));
657 tempCS01
= vec_perm(cosv
, sinv
, vcprm(0,s0
,1,s1
));
658 temp1
= vec_madd(temp0022
, tempCS01
, vczero
);
659 tempCS01
= vec_perm(cosv
, sinv
, vcprm(s0
,0,s1
,1));
660 temp2
= vec_madd(temp1133
, tempCS01
, vczero
);
661 bufv_0
= vec_madd(temp2
, vcii(p
,n
,p
,n
), temp1
);
663 vec_st(bufv_0
, (i
+ 0) << 3, (float*)buf
);
665 /* idem with bufv_2 and high-order cosv/sinv */
667 temp0022
= vec_perm(bufv_2
, bufv_2
, vcprm(0,0,2,2));
668 temp1133
= vec_perm(bufv_2
, bufv_2
, vcprm(1,1,3,3));
669 tempCS01
= vec_perm(cosv
, sinv
, vcprm(2,s2
,3,s3
));
670 temp1
= vec_madd(temp0022
, tempCS01
, vczero
);
671 tempCS01
= vec_perm(cosv
, sinv
, vcprm(s2
,2,s3
,3));
672 temp2
= vec_madd(temp1133
, tempCS01
, vczero
);
673 bufv_2
= vec_madd(temp2
, vcii(p
,n
,p
,n
), temp1
);
675 vec_st(bufv_2
, (i
+ 2) << 3, (float*)buf
);
682 window_ptr
= a52_imdct_window
;
684 /* Window and convert to real valued signal */
685 for(i
=0; i
< 64; i
++) {
686 *data_ptr
++ = -buf
[64+i
].imag
* *window_ptr
++ + *delay_ptr
++ + bias
;
687 *data_ptr
++ = buf
[64-i
-1].real
* *window_ptr
++ + *delay_ptr
++ + bias
;
690 for(i
=0; i
< 64; i
++) {
691 *data_ptr
++ = -buf
[i
].real
* *window_ptr
++ + *delay_ptr
++ + bias
;
692 *data_ptr
++ = buf
[128-i
-1].imag
* *window_ptr
++ + *delay_ptr
++ + bias
;
695 /* The trailing edge of the window goes into the delay line */
698 for(i
=0; i
< 64; i
++) {
699 *delay_ptr
++ = -buf
[64+i
].real
* *--window_ptr
;
700 *delay_ptr
++ = buf
[64-i
-1].imag
* *--window_ptr
;
703 for(i
=0; i
<64; i
++) {
704 *delay_ptr
++ = buf
[i
].imag
* *--window_ptr
;
705 *delay_ptr
++ = -buf
[128-i
-1].real
* *--window_ptr
;
711 // Stuff below this line is borrowed from libac3
713 #if defined(ARCH_X86) || defined(ARCH_X86_64)
717 #include "srfftp_3dnow.h"
719 const i_cmplx_t x_plus_minus_3dnow
__attribute__ ((aligned (8))) = {{ 0x00000000UL
, 0x80000000UL
}};
720 const i_cmplx_t x_minus_plus_3dnow
__attribute__ ((aligned (8))) = {{ 0x80000000UL
, 0x00000000UL
}};
721 const complex_t HSQRT2_3DNOW
__attribute__ ((aligned (8))) = { 0.707106781188, 0.707106781188 };
724 #include "imdct_3dnow.h"
726 #include "imdct_3dnow.h"
729 imdct_do_512_sse(sample_t data
[],sample_t delay
[], sample_t bias
)
736 long two_m_plus_one_shl3
;
737 complex_t
*buf_offset
;
746 sample_t
*window_ptr
;
748 /* 512 IMDCT with source and dest data in 'data' */
749 /* see the c version (dct_do_512()), its allmost identical, just in C */
751 /* Pre IFFT complex multiply plus IFFT cmplx conjugate */
752 /* Bit reversed shuffling */
754 "xor %%"REG_S
", %%"REG_S
" \n\t"
755 "lea "MANGLE(bit_reverse_512
)", %%"REG_a
"\n\t"
756 "mov $1008, %%"REG_D
" \n\t"
757 "push %%"REG_BP
" \n\t" //use ebp without telling gcc
760 "movlps (%0, %%"REG_S
"), %%xmm0 \n\t" // XXXI
761 "movhps 8(%0, %%"REG_D
"), %%xmm0 \n\t" // RXXI
762 "movlps 8(%0, %%"REG_S
"), %%xmm1 \n\t" // XXXi
763 "movhps (%0, %%"REG_D
"), %%xmm1 \n\t" // rXXi
764 "shufps $0x33, %%xmm1, %%xmm0 \n\t" // irIR
765 "movaps "MANGLE(sseSinCos1c
)"(%%"REG_S
"), %%xmm2\n\t"
766 "mulps %%xmm0, %%xmm2 \n\t"
767 "shufps $0xB1, %%xmm0, %%xmm0 \n\t" // riRI
768 "mulps "MANGLE(sseSinCos1d
)"(%%"REG_S
"), %%xmm0\n\t"
769 "subps %%xmm0, %%xmm2 \n\t"
770 "movzb (%%"REG_a
"), %%"REG_d
" \n\t"
771 "movzb 1(%%"REG_a
"), %%"REG_BP
" \n\t"
772 "movlps %%xmm2, (%1, %%"REG_d
", 8) \n\t"
773 "movhps %%xmm2, (%1, %%"REG_BP
", 8) \n\t"
774 "add $16, %%"REG_S
" \n\t"
775 "add $2, %%"REG_a
" \n\t" // avoid complex addressing for P4 crap
776 "sub $16, %%"REG_D
" \n\t"
778 "pop %%"REG_BP
" \n\t"//no we didnt touch ebp *g*
779 :: "b" (data
), "c" (buf
)
780 : "%"REG_S
, "%"REG_D
, "%"REG_a
, "%"REG_d
785 /* unoptimized variant
786 for (m=1; m < 7; m++) {
792 two_m_plus_one = (1 << (m+1));
794 for(i = 0; i < 128; i += two_m_plus_one) {
795 for(k = 0; k < two_m; k++) {
798 tmp_a_r = buf[p].real;
799 tmp_a_i = buf[p].imag;
800 tmp_b_r = buf[q].real * w[m][k].real - buf[q].imag * w[m][k].imag;
801 tmp_b_i = buf[q].imag * w[m][k].real + buf[q].real * w[m][k].imag;
802 buf[p].real = tmp_a_r + tmp_b_r;
803 buf[p].imag = tmp_a_i + tmp_b_i;
804 buf[q].real = tmp_a_r - tmp_b_r;
805 buf[q].imag = tmp_a_i - tmp_b_i;
812 // Note w[0][0]={1,0}
814 "xorps %%xmm1, %%xmm1 \n\t"
815 "xorps %%xmm2, %%xmm2 \n\t"
816 "mov %0, %%"REG_S
" \n\t"
819 "movlps (%%"REG_S
"), %%xmm0\n\t" //buf[p]
820 "movlps 8(%%"REG_S
"), %%xmm1\n\t" //buf[q]
821 "movhps (%%"REG_S
"), %%xmm0\n\t" //buf[p]
822 "movhps 8(%%"REG_S
"), %%xmm2\n\t" //buf[q]
823 "addps %%xmm1, %%xmm0 \n\t"
824 "subps %%xmm2, %%xmm0 \n\t"
825 "movaps %%xmm0, (%%"REG_S
")\n\t"
826 "add $16, %%"REG_S
" \n\t"
827 "cmp %1, %%"REG_S
" \n\t"
829 :: "g" (buf
), "r" (buf
+ 128)
834 // Note w[1]={{1,0}, {0,-1}}
836 "movaps "MANGLE(ps111_1
)", %%xmm7\n\t" // 1,1,1,-1
837 "mov %0, %%"REG_S
" \n\t"
840 "movaps 16(%%"REG_S
"), %%xmm2 \n\t" //r2,i2,r3,i3
841 "shufps $0xB4, %%xmm2, %%xmm2 \n\t" //r2,i2,i3,r3
842 "mulps %%xmm7, %%xmm2 \n\t" //r2,i2,i3,-r3
843 "movaps (%%"REG_S
"), %%xmm0 \n\t" //r0,i0,r1,i1
844 "movaps (%%"REG_S
"), %%xmm1 \n\t" //r0,i0,r1,i1
845 "addps %%xmm2, %%xmm0 \n\t"
846 "subps %%xmm2, %%xmm1 \n\t"
847 "movaps %%xmm0, (%%"REG_S
") \n\t"
848 "movaps %%xmm1, 16(%%"REG_S
") \n\t"
849 "add $32, %%"REG_S
" \n\t"
850 "cmp %1, %%"REG_S
" \n\t"
852 :: "g" (buf
), "r" (buf
+ 128)
858 Note sseW2+0={1,1,sqrt(2),sqrt(2))
859 Note sseW2+16={0,0,sqrt(2),-sqrt(2))
860 Note sseW2+32={0,0,-sqrt(2),-sqrt(2))
861 Note sseW2+48={1,-1,sqrt(2),-sqrt(2))
864 "movaps 48+"MANGLE(sseW2
)", %%xmm6\n\t"
865 "movaps 16+"MANGLE(sseW2
)", %%xmm7\n\t"
866 "xorps %%xmm5, %%xmm5 \n\t"
867 "xorps %%xmm2, %%xmm2 \n\t"
868 "mov %0, %%"REG_S
" \n\t"
871 "movaps 32(%%"REG_S
"), %%xmm2 \n\t" //r4,i4,r5,i5
872 "movaps 48(%%"REG_S
"), %%xmm3 \n\t" //r6,i6,r7,i7
873 "movaps "MANGLE(sseW2
)", %%xmm4 \n\t" //r4,i4,r5,i5
874 "movaps 32+"MANGLE(sseW2
)", %%xmm5\n\t" //r6,i6,r7,i7
875 "mulps %%xmm2, %%xmm4 \n\t"
876 "mulps %%xmm3, %%xmm5 \n\t"
877 "shufps $0xB1, %%xmm2, %%xmm2 \n\t" //i4,r4,i5,r5
878 "shufps $0xB1, %%xmm3, %%xmm3 \n\t" //i6,r6,i7,r7
879 "mulps %%xmm6, %%xmm3 \n\t"
880 "mulps %%xmm7, %%xmm2 \n\t"
881 "movaps (%%"REG_S
"), %%xmm0 \n\t" //r0,i0,r1,i1
882 "movaps 16(%%"REG_S
"), %%xmm1 \n\t" //r2,i2,r3,i3
883 "addps %%xmm4, %%xmm2 \n\t"
884 "addps %%xmm5, %%xmm3 \n\t"
885 "movaps %%xmm2, %%xmm4 \n\t"
886 "movaps %%xmm3, %%xmm5 \n\t"
887 "addps %%xmm0, %%xmm2 \n\t"
888 "addps %%xmm1, %%xmm3 \n\t"
889 "subps %%xmm4, %%xmm0 \n\t"
890 "subps %%xmm5, %%xmm1 \n\t"
891 "movaps %%xmm2, (%%"REG_S
") \n\t"
892 "movaps %%xmm3, 16(%%"REG_S
") \n\t"
893 "movaps %%xmm0, 32(%%"REG_S
") \n\t"
894 "movaps %%xmm1, 48(%%"REG_S
") \n\t"
895 "add $64, %%"REG_S
" \n\t"
896 "cmp %1, %%"REG_S
" \n\t"
898 :: "g" (buf
), "r" (buf
+ 128)
902 /* 4-7. iterations */
903 for (m
=3; m
< 7; m
++) {
905 two_m_plus_one
= two_m
<<1;
906 two_m_plus_one_shl3
= (two_m_plus_one
<<3);
907 buf_offset
= buf
+128;
909 "mov %0, %%"REG_S
" \n\t"
912 "xor %%"REG_D
", %%"REG_D
" \n\t" // k
913 "lea (%%"REG_S
", %3), %%"REG_d
" \n\t"
915 "movaps (%%"REG_d
", %%"REG_D
"), %%xmm1 \n\t"
916 "movaps (%4, %%"REG_D
", 2), %%xmm2 \n\t"
917 "mulps %%xmm1, %%xmm2 \n\t"
918 "shufps $0xB1, %%xmm1, %%xmm1 \n\t"
919 "mulps 16(%4, %%"REG_D
", 2), %%xmm1 \n\t"
920 "movaps (%%"REG_S
", %%"REG_D
"), %%xmm0 \n\t"
921 "addps %%xmm2, %%xmm1 \n\t"
922 "movaps %%xmm1, %%xmm2 \n\t"
923 "addps %%xmm0, %%xmm1 \n\t"
924 "subps %%xmm2, %%xmm0 \n\t"
925 "movaps %%xmm1, (%%"REG_S
", %%"REG_D
") \n\t"
926 "movaps %%xmm0, (%%"REG_d
", %%"REG_D
") \n\t"
927 "add $16, %%"REG_D
" \n\t"
928 "cmp %3, %%"REG_D
" \n\t" //FIXME (opt) count against 0
930 "add %2, %%"REG_S
" \n\t"
931 "cmp %1, %%"REG_S
" \n\t"
933 :: "g" (buf
), "m" (buf_offset
), "m" (two_m_plus_one_shl3
), "r" (two_m
<<3),
935 : "%"REG_S
, "%"REG_D
, "%"REG_d
939 /* Post IFFT complex multiply plus IFFT complex conjugate*/
941 "mov $-1024, %%"REG_S
" \n\t"
944 "movaps (%0, %%"REG_S
"), %%xmm0 \n\t"
945 "movaps (%0, %%"REG_S
"), %%xmm1 \n\t"
946 "shufps $0xB1, %%xmm0, %%xmm0 \n\t"
947 "mulps 1024+"MANGLE(sseSinCos1c
)"(%%"REG_S
"), %%xmm1\n\t"
948 "mulps 1024+"MANGLE(sseSinCos1d
)"(%%"REG_S
"), %%xmm0\n\t"
949 "addps %%xmm1, %%xmm0 \n\t"
950 "movaps %%xmm0, (%0, %%"REG_S
") \n\t"
951 "add $16, %%"REG_S
" \n\t"
960 window_ptr
= a52_imdct_window
;
962 /* Window and convert to real valued signal */
964 "xor %%"REG_D
", %%"REG_D
" \n\t" // 0
965 "xor %%"REG_S
", %%"REG_S
" \n\t" // 0
966 "movss %3, %%xmm2 \n\t" // bias
967 "shufps $0x00, %%xmm2, %%xmm2 \n\t" // bias, bias, ...
970 "movlps (%0, %%"REG_S
"), %%xmm0 \n\t" // ? ? A ?
971 "movlps 8(%0, %%"REG_S
"), %%xmm1 \n\t" // ? ? C ?
972 "movhps -16(%0, %%"REG_D
"), %%xmm1 \n\t" // ? D C ?
973 "movhps -8(%0, %%"REG_D
"), %%xmm0 \n\t" // ? B A ?
974 "shufps $0x99, %%xmm1, %%xmm0 \n\t" // D C B A
975 "mulps "MANGLE(sseWindow
)"(%%"REG_S
"), %%xmm0\n\t"
976 "addps (%2, %%"REG_S
"), %%xmm0 \n\t"
977 "addps %%xmm2, %%xmm0 \n\t"
978 "movaps %%xmm0, (%1, %%"REG_S
") \n\t"
979 "add $16, %%"REG_S
" \n\t"
980 "sub $16, %%"REG_D
" \n\t"
981 "cmp $512, %%"REG_S
" \n\t"
983 :: "r" (buf
+64), "r" (data_ptr
), "r" (delay_ptr
), "m" (bias
)
991 "mov $1024, %%"REG_D
" \n\t" // 512
992 "xor %%"REG_S
", %%"REG_S
" \n\t" // 0
993 "movss %3, %%xmm2 \n\t" // bias
994 "shufps $0x00, %%xmm2, %%xmm2 \n\t" // bias, bias, ...
997 "movlps (%0, %%"REG_S
"), %%xmm0 \n\t" // ? ? ? A
998 "movlps 8(%0, %%"REG_S
"), %%xmm1 \n\t" // ? ? ? C
999 "movhps -16(%0, %%"REG_D
"), %%xmm1 \n\t" // D ? ? C
1000 "movhps -8(%0, %%"REG_D
"), %%xmm0 \n\t" // B ? ? A
1001 "shufps $0xCC, %%xmm1, %%xmm0 \n\t" // D C B A
1002 "mulps 512+"MANGLE(sseWindow
)"(%%"REG_S
"), %%xmm0\n\t"
1003 "addps (%2, %%"REG_S
"), %%xmm0 \n\t"
1004 "addps %%xmm2, %%xmm0 \n\t"
1005 "movaps %%xmm0, (%1, %%"REG_S
") \n\t"
1006 "add $16, %%"REG_S
" \n\t"
1007 "sub $16, %%"REG_D
" \n\t"
1008 "cmp $512, %%"REG_S
" \n\t"
1010 :: "r" (buf
), "r" (data_ptr
), "r" (delay_ptr
), "m" (bias
)
1011 : "%"REG_S
, "%"REG_D
1016 /* The trailing edge of the window goes into the delay line */
1020 "xor %%"REG_D
", %%"REG_D
" \n\t" // 0
1021 "xor %%"REG_S
", %%"REG_S
" \n\t" // 0
1024 "movlps (%0, %%"REG_S
"), %%xmm0 \n\t" // ? ? ? A
1025 "movlps 8(%0, %%"REG_S
"), %%xmm1 \n\t" // ? ? ? C
1026 "movhps -16(%0, %%"REG_D
"), %%xmm1 \n\t" // D ? ? C
1027 "movhps -8(%0, %%"REG_D
"), %%xmm0 \n\t" // B ? ? A
1028 "shufps $0xCC, %%xmm1, %%xmm0 \n\t" // D C B A
1029 "mulps 1024+"MANGLE(sseWindow
)"(%%"REG_S
"), %%xmm0\n\t"
1030 "movaps %%xmm0, (%1, %%"REG_S
") \n\t"
1031 "add $16, %%"REG_S
" \n\t"
1032 "sub $16, %%"REG_D
" \n\t"
1033 "cmp $512, %%"REG_S
" \n\t"
1035 :: "r" (buf
+64), "r" (delay_ptr
)
1036 : "%"REG_S
, "%"REG_D
1042 "mov $1024, %%"REG_D
" \n\t" // 1024
1043 "xor %%"REG_S
", %%"REG_S
" \n\t" // 0
1046 "movlps (%0, %%"REG_S
"), %%xmm0 \n\t" // ? ? A ?
1047 "movlps 8(%0, %%"REG_S
"), %%xmm1 \n\t" // ? ? C ?
1048 "movhps -16(%0, %%"REG_D
"), %%xmm1 \n\t" // ? D C ?
1049 "movhps -8(%0, %%"REG_D
"), %%xmm0 \n\t" // ? B A ?
1050 "shufps $0x99, %%xmm1, %%xmm0 \n\t" // D C B A
1051 "mulps 1536+"MANGLE(sseWindow
)"(%%"REG_S
"), %%xmm0\n\t"
1052 "movaps %%xmm0, (%1, %%"REG_S
") \n\t"
1053 "add $16, %%"REG_S
" \n\t"
1054 "sub $16, %%"REG_D
" \n\t"
1055 "cmp $512, %%"REG_S
" \n\t"
1057 :: "r" (buf
), "r" (delay_ptr
)
1058 : "%"REG_S
, "%"REG_D
1061 #endif // ARCH_X86 || ARCH_X86_64
1063 void a52_imdct_256(sample_t
* data
, sample_t
* delay
, sample_t bias
)
1066 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
;
1067 const sample_t
* window
= a52_imdct_window
;
1068 complex_t buf1
[64], buf2
[64];
1070 /* Pre IFFT complex multiply plus IFFT cmplx conjugate */
1071 for (i
= 0; i
< 64; i
++) {
1076 buf1
[i
].real
= t_i
* data
[254-k
] + t_r
* data
[k
];
1077 buf1
[i
].imag
= t_r
* data
[254-k
] - t_i
* data
[k
];
1079 buf2
[i
].real
= t_i
* data
[255-k
] + t_r
* data
[k
+1];
1080 buf2
[i
].imag
= t_r
* data
[255-k
] - t_i
* data
[k
+1];
1086 /* Post IFFT complex multiply */
1087 /* Window and convert to real valued signal */
1088 for (i
= 0; i
< 32; i
++) {
1089 /* y1[n] = z1[n] * (xcos2[n] + j * xs in2[n]) ; */
1090 t_r
= post2
[i
].real
;
1091 t_i
= post2
[i
].imag
;
1093 a_r
= t_r
* buf1
[i
].real
+ t_i
* buf1
[i
].imag
;
1094 a_i
= t_i
* buf1
[i
].real
- t_r
* buf1
[i
].imag
;
1095 b_r
= t_i
* buf1
[63-i
].real
+ t_r
* buf1
[63-i
].imag
;
1096 b_i
= t_r
* buf1
[63-i
].real
- t_i
* buf1
[63-i
].imag
;
1098 c_r
= t_r
* buf2
[i
].real
+ t_i
* buf2
[i
].imag
;
1099 c_i
= t_i
* buf2
[i
].real
- t_r
* buf2
[i
].imag
;
1100 d_r
= t_i
* buf2
[63-i
].real
+ t_r
* buf2
[63-i
].imag
;
1101 d_i
= t_r
* buf2
[63-i
].real
- t_i
* buf2
[63-i
].imag
;
1104 w_2
= window
[255-2*i
];
1105 data
[2*i
] = delay
[2*i
] * w_2
- a_r
* w_1
+ bias
;
1106 data
[255-2*i
] = delay
[2*i
] * w_1
+ a_r
* w_2
+ bias
;
1109 w_1
= window
[128+2*i
];
1110 w_2
= window
[127-2*i
];
1111 data
[128+2*i
] = delay
[127-2*i
] * w_2
+ a_i
* w_1
+ bias
;
1112 data
[127-2*i
] = delay
[127-2*i
] * w_1
- a_i
* w_2
+ bias
;
1113 delay
[127-2*i
] = c_r
;
1115 w_1
= window
[2*i
+1];
1116 w_2
= window
[254-2*i
];
1117 data
[2*i
+1] = delay
[2*i
+1] * w_2
- b_i
* w_1
+ bias
;
1118 data
[254-2*i
] = delay
[2*i
+1] * w_1
+ b_i
* w_2
+ bias
;
1121 w_1
= window
[129+2*i
];
1122 w_2
= window
[126-2*i
];
1123 data
[129+2*i
] = delay
[126-2*i
] * w_2
+ b_r
* w_1
+ bias
;
1124 data
[126-2*i
] = delay
[126-2*i
] * w_1
- b_r
* w_2
+ bias
;
1125 delay
[126-2*i
] = d_i
;
1129 static double besselI0 (double x
)
1135 bessel
= bessel
* x
/ (i
* i
) + 1;
1140 void a52_imdct_init (uint32_t mm_accel
)
1145 /* compute imdct window - kaiser-bessel derived window, alpha = 5.0 */
1147 for (i
= 0; i
< 256; i
++) {
1148 sum
+= besselI0 (i
* (256 - i
) * (5 * M_PI
/ 256) * (5 * M_PI
/ 256));
1149 a52_imdct_window
[i
] = sum
;
1152 for (i
= 0; i
< 256; i
++)
1153 a52_imdct_window
[i
] = sqrt (a52_imdct_window
[i
] / sum
);
1155 for (i
= 0; i
< 3; i
++)
1156 roots16
[i
] = cos ((M_PI
/ 8) * (i
+ 1));
1158 for (i
= 0; i
< 7; i
++)
1159 roots32
[i
] = cos ((M_PI
/ 16) * (i
+ 1));
1161 for (i
= 0; i
< 15; i
++)
1162 roots64
[i
] = cos ((M_PI
/ 32) * (i
+ 1));
1164 for (i
= 0; i
< 31; i
++)
1165 roots128
[i
] = cos ((M_PI
/ 64) * (i
+ 1));
1167 for (i
= 0; i
< 64; i
++) {
1168 k
= fftorder
[i
] / 2 + 64;
1169 pre1
[i
].real
= cos ((M_PI
/ 256) * (k
- 0.25));
1170 pre1
[i
].imag
= sin ((M_PI
/ 256) * (k
- 0.25));
1173 for (i
= 64; i
< 128; i
++) {
1174 k
= fftorder
[i
] / 2 + 64;
1175 pre1
[i
].real
= -cos ((M_PI
/ 256) * (k
- 0.25));
1176 pre1
[i
].imag
= -sin ((M_PI
/ 256) * (k
- 0.25));
1179 for (i
= 0; i
< 64; i
++) {
1180 post1
[i
].real
= cos ((M_PI
/ 256) * (i
+ 0.5));
1181 post1
[i
].imag
= sin ((M_PI
/ 256) * (i
+ 0.5));
1184 for (i
= 0; i
< 64; i
++) {
1185 k
= fftorder
[i
] / 4;
1186 pre2
[i
].real
= cos ((M_PI
/ 128) * (k
- 0.25));
1187 pre2
[i
].imag
= sin ((M_PI
/ 128) * (k
- 0.25));
1190 for (i
= 0; i
< 32; i
++) {
1191 post2
[i
].real
= cos ((M_PI
/ 128) * (i
+ 0.5));
1192 post2
[i
].imag
= sin ((M_PI
/ 128) * (i
+ 0.5));
1194 for (i
= 0; i
< 128; i
++) {
1195 xcos1
[i
] = -cos ((M_PI
/ 2048) * (8 * i
+ 1));
1196 xsin1
[i
] = -sin ((M_PI
/ 2048) * (8 * i
+ 1));
1198 for (i
= 0; i
< 7; i
++) {
1200 for (k
= 0; k
< j
; k
++) {
1201 w
[i
][k
].real
= cos (-M_PI
* k
/ j
);
1202 w
[i
][k
].imag
= sin (-M_PI
* k
/ j
);
1205 #if defined(ARCH_X86) || defined(ARCH_X86_64)
1206 for (i
= 0; i
< 128; i
++) {
1207 sseSinCos1c
[2*i
+0]= xcos1
[i
];
1208 sseSinCos1c
[2*i
+1]= -xcos1
[i
];
1209 sseSinCos1d
[2*i
+0]= xsin1
[i
];
1210 sseSinCos1d
[2*i
+1]= xsin1
[i
];
1212 for (i
= 1; i
< 7; i
++) {
1214 for (k
= 0; k
< j
; k
+=2) {
1216 sseW
[i
][4*k
+ 0] = w
[i
][k
+0].real
;
1217 sseW
[i
][4*k
+ 1] = w
[i
][k
+0].real
;
1218 sseW
[i
][4*k
+ 2] = w
[i
][k
+1].real
;
1219 sseW
[i
][4*k
+ 3] = w
[i
][k
+1].real
;
1221 sseW
[i
][4*k
+ 4] = -w
[i
][k
+0].imag
;
1222 sseW
[i
][4*k
+ 5] = w
[i
][k
+0].imag
;
1223 sseW
[i
][4*k
+ 6] = -w
[i
][k
+1].imag
;
1224 sseW
[i
][4*k
+ 7] = w
[i
][k
+1].imag
;
1226 //we multiply more or less uninitalized numbers so we need to use exactly 0.0
1229 // sseW[i][4*k + 0]= sseW[i][4*k + 1]= 1.0;
1230 sseW
[i
][4*k
+ 4]= sseW
[i
][4*k
+ 5]= 0.0;
1235 sseW
[i
][4*k
+ 0]= sseW
[i
][4*k
+ 1]= 0.0;
1236 // sseW[i][4*k + 4]= -(sseW[i][4*k + 5]= -1.0);
1241 for(i
=0; i
<128; i
++)
1243 sseWindow
[2*i
+0]= -a52_imdct_window
[2*i
+0];
1244 sseWindow
[2*i
+1]= a52_imdct_window
[2*i
+1];
1249 sseWindow
[256 + 2*i
+0]= -a52_imdct_window
[254 - 2*i
+1];
1250 sseWindow
[256 + 2*i
+1]= a52_imdct_window
[254 - 2*i
+0];
1251 sseWindow
[384 + 2*i
+0]= a52_imdct_window
[126 - 2*i
+1];
1252 sseWindow
[384 + 2*i
+1]= -a52_imdct_window
[126 - 2*i
+0];
1255 a52_imdct_512
= imdct_do_512
;
1256 ifft128
= ifft128_c
;
1259 #if defined(ARCH_X86) || defined(ARCH_X86_64)
1260 if(mm_accel
& MM_ACCEL_X86_SSE
)
1262 fprintf (stderr
, "Using SSE optimized IMDCT transform\n");
1263 a52_imdct_512
= imdct_do_512_sse
;
1266 if(mm_accel
& MM_ACCEL_X86_3DNOWEXT
)
1268 fprintf (stderr
, "Using 3DNowEx optimized IMDCT transform\n");
1269 a52_imdct_512
= imdct_do_512_3dnowex
;
1272 if(mm_accel
& MM_ACCEL_X86_3DNOW
)
1274 fprintf (stderr
, "Using 3DNow optimized IMDCT transform\n");
1275 a52_imdct_512
= imdct_do_512_3dnow
;
1278 #endif // ARCH_X86 || ARCH_X86_64
1280 if (mm_accel
& MM_ACCEL_PPC_ALTIVEC
)
1282 fprintf(stderr
, "Using AltiVec optimized IMDCT transform\n");
1283 a52_imdct_512
= imdct_do_512_altivec
;
1288 #ifdef LIBA52_DJBFFT
1289 if (mm_accel
& MM_ACCEL_DJBFFT
) {
1290 fprintf (stderr
, "Using djbfft for IMDCT transform\n");
1291 ifft128
= (void (*) (complex_t
*)) fftc4_un128
;
1292 ifft64
= (void (*) (complex_t
*)) fftc4_un64
;
1296 fprintf (stderr
, "No accelerated IMDCT transform found\n");