summaryrefslogtreecommitdiff
path: root/apps/codecs/libwma/mdct2.c
diff options
context:
space:
mode:
Diffstat (limited to 'apps/codecs/libwma/mdct2.c')
-rw-r--r--apps/codecs/libwma/mdct2.c522
1 files changed, 0 insertions, 522 deletions
diff --git a/apps/codecs/libwma/mdct2.c b/apps/codecs/libwma/mdct2.c
deleted file mode 100644
index c492403b41..0000000000
--- a/apps/codecs/libwma/mdct2.c
+++ /dev/null
@@ -1,522 +0,0 @@
1/********************************************************************
2 * *
3 * THIS FILE IS PART OF THE OggVorbis 'TREMOR' CODEC SOURCE CODE. *
4 * *
5 * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS *
6 * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
7 * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. *
8 * *
9 * THE OggVorbis 'TREMOR' SOURCE CODE IS (C) COPYRIGHT 1994-2002 *
10 * BY THE Xiph.Org FOUNDATION http://www.xiph.org/ *
11 * *
12 ********************************************************************
13
14 function: normalized modified discrete cosine transform
15 power of two length transform only [64 <= n ]
16
17
18 Original algorithm adapted long ago from _The use of multirate filter
19 banks for coding of high quality digital audio_, by T. Sporer,
20 K. Brandenburg and B. Edler, collection of the European Signal
21 Processing Conference (EUSIPCO), Amsterdam, June 1992, Vol.1, pp
22 211-214
23
24 The below code implements an algorithm that no longer looks much like
25 that presented in the paper, but the basic structure remains if you
26 dig deep enough to see it.
27
28 This module DOES NOT INCLUDE code to generate/apply the window
29 function. Everybody has their own weird favorite including me... I
30 happen to like the properties of y=sin(.5PI*sin^2(x)), but others may
31 vehemently disagree.
32
33 ********************************************************************/
34
35/*Tremor IMDCT adapted for use with libwmai*/
36
37
38#include "mdct2.h"
39#include "mdct_lookup.h"
40#include <codecs/lib/codeclib.h>
41
42#if defined(CPU_ARM) && CONFIG_CPU != S3C2440
43/* C code is faster on S3C2440 */
44
45extern void mdct_butterfly_32(ogg_int32_t *x);
46extern void mdct_butterfly_generic_loop(ogg_int32_t *x1, ogg_int32_t *x2,
47 LOOKUP_T *T0, int step,
48 LOOKUP_T *Ttop);
49
50static inline void mdct_butterfly_generic(ogg_int32_t *x,int points, int step){
51 mdct_butterfly_generic_loop(x + points, x + (points>>1),
52 sincos_lookup0, step, sincos_lookup0+1024);
53}
54
55#else
56
57/* 8 point butterfly (in place) */
58static inline void mdct_butterfly_8(ogg_int32_t *x){
59 register ogg_int32_t r0 = x[4] + x[0];
60 register ogg_int32_t r1 = x[4] - x[0];
61 register ogg_int32_t r2 = x[5] + x[1];
62 register ogg_int32_t r3 = x[5] - x[1];
63 register ogg_int32_t r4 = x[6] + x[2];
64 register ogg_int32_t r5 = x[6] - x[2];
65 register ogg_int32_t r6 = x[7] + x[3];
66 register ogg_int32_t r7 = x[7] - x[3];
67
68 x[0] = r5 + r3;
69 x[1] = r7 - r1;
70 x[2] = r5 - r3;
71 x[3] = r7 + r1;
72 x[4] = r4 - r0;
73 x[5] = r6 - r2;
74 x[6] = r4 + r0;
75 x[7] = r6 + r2;
76 MB();
77}
78
79/* 16 point butterfly (in place, 4 register) */
80static inline void mdct_butterfly_16(ogg_int32_t *x){
81
82 register ogg_int32_t r0, r1;
83
84 r0 = x[ 0] - x[ 8]; x[ 8] += x[ 0];
85 r1 = x[ 1] - x[ 9]; x[ 9] += x[ 1];
86 x[ 0] = MULT31((r0 + r1) , cPI2_8);
87 x[ 1] = MULT31((r1 - r0) , cPI2_8);
88 MB();
89
90 r0 = x[10] - x[ 2]; x[10] += x[ 2];
91 r1 = x[ 3] - x[11]; x[11] += x[ 3];
92 x[ 2] = r1; x[ 3] = r0;
93 MB();
94
95 r0 = x[12] - x[ 4]; x[12] += x[ 4];
96 r1 = x[13] - x[ 5]; x[13] += x[ 5];
97 x[ 4] = MULT31((r0 - r1) , cPI2_8);
98 x[ 5] = MULT31((r0 + r1) , cPI2_8);
99 MB();
100
101 r0 = x[14] - x[ 6]; x[14] += x[ 6];
102 r1 = x[15] - x[ 7]; x[15] += x[ 7];
103 x[ 6] = r0; x[ 7] = r1;
104 MB();
105
106 mdct_butterfly_8(x);
107 mdct_butterfly_8(x+8);
108}
109
110/* 32 point butterfly (in place, 4 register) */
111static inline void mdct_butterfly_32(ogg_int32_t *x){
112
113 register ogg_int32_t r0, r1;
114
115 r0 = x[30] - x[14]; x[30] += x[14];
116 r1 = x[31] - x[15]; x[31] += x[15];
117 x[14] = r0; x[15] = r1;
118 MB();
119
120 r0 = x[28] - x[12]; x[28] += x[12];
121 r1 = x[29] - x[13]; x[29] += x[13];
122 XNPROD31( r0, r1, cPI1_8, cPI3_8, &x[12], &x[13] );
123 MB();
124
125 r0 = x[26] - x[10]; x[26] += x[10];
126 r1 = x[27] - x[11]; x[27] += x[11];
127 x[10] = MULT31((r0 - r1) , cPI2_8);
128 x[11] = MULT31((r0 + r1) , cPI2_8);
129 MB();
130
131 r0 = x[24] - x[ 8]; x[24] += x[ 8];
132 r1 = x[25] - x[ 9]; x[25] += x[ 9];
133 XNPROD31( r0, r1, cPI3_8, cPI1_8, &x[ 8], &x[ 9] );
134 MB();
135
136 r0 = x[22] - x[ 6]; x[22] += x[ 6];
137 r1 = x[ 7] - x[23]; x[23] += x[ 7];
138 x[ 6] = r1; x[ 7] = r0;
139 MB();
140
141 r0 = x[ 4] - x[20]; x[20] += x[ 4];
142 r1 = x[ 5] - x[21]; x[21] += x[ 5];
143 XPROD31 ( r0, r1, cPI3_8, cPI1_8, &x[ 4], &x[ 5] );
144 MB();
145
146 r0 = x[ 2] - x[18]; x[18] += x[ 2];
147 r1 = x[ 3] - x[19]; x[19] += x[ 3];
148 x[ 2] = MULT31((r1 + r0) , cPI2_8);
149 x[ 3] = MULT31((r1 - r0) , cPI2_8);
150 MB();
151
152 r0 = x[ 0] - x[16]; x[16] += x[ 0];
153 r1 = x[ 1] - x[17]; x[17] += x[ 1];
154 XPROD31 ( r0, r1, cPI1_8, cPI3_8, &x[ 0], &x[ 1] );
155 MB();
156
157 mdct_butterfly_16(x);
158 mdct_butterfly_16(x+16);
159}
160
161/* N/stage point generic N stage butterfly (in place, 4 register) */
162void mdct_butterfly_generic(ogg_int32_t *x,int points, int step)
163 ICODE_ATTR_TREMOR_MDCT;
164void mdct_butterfly_generic(ogg_int32_t *x,int points, int step){
165 LOOKUP_T *T = sincos_lookup0;
166 ogg_int32_t *x1 = x + points - 8;
167 ogg_int32_t *x2 = x + (points>>1) - 8;
168 register ogg_int32_t r0;
169 register ogg_int32_t r1;
170 register ogg_int32_t r2;
171 register ogg_int32_t r3;
172
173 do{
174 r0 = x1[6] - x2[6]; x1[6] += x2[6];
175 r1 = x2[7] - x1[7]; x1[7] += x2[7];
176 r2 = x1[4] - x2[4]; x1[4] += x2[4];
177 r3 = x2[5] - x1[5]; x1[5] += x2[5];
178 XPROD31( r1, r0, T[0], T[1], &x2[6], &x2[7] ); T+=step;
179 XPROD31( r3, r2, T[0], T[1], &x2[4], &x2[5] ); T+=step;
180
181 r0 = x1[2] - x2[2]; x1[2] += x2[2];
182 r1 = x2[3] - x1[3]; x1[3] += x2[3];
183 r2 = x1[0] - x2[0]; x1[0] += x2[0];
184 r3 = x2[1] - x1[1]; x1[1] += x2[1];
185 XPROD31( r1, r0, T[0], T[1], &x2[2], &x2[3] ); T+=step;
186 XPROD31( r3, r2, T[0], T[1], &x2[0], &x2[1] ); T+=step;
187
188 x1-=8; x2-=8;
189 }while(T<sincos_lookup0+1024);
190 do{
191 r0 = x1[6] - x2[6]; x1[6] += x2[6];
192 r1 = x1[7] - x2[7]; x1[7] += x2[7];
193 r2 = x1[4] - x2[4]; x1[4] += x2[4];
194 r3 = x1[5] - x2[5]; x1[5] += x2[5];
195 XNPROD31( r0, r1, T[0], T[1], &x2[6], &x2[7] ); T-=step;
196 XNPROD31( r2, r3, T[0], T[1], &x2[4], &x2[5] ); T-=step;
197
198 r0 = x1[2] - x2[2]; x1[2] += x2[2];
199 r1 = x1[3] - x2[3]; x1[3] += x2[3];
200 r2 = x1[0] - x2[0]; x1[0] += x2[0];
201 r3 = x1[1] - x2[1]; x1[1] += x2[1];
202 XNPROD31( r0, r1, T[0], T[1], &x2[2], &x2[3] ); T-=step;
203 XNPROD31( r2, r3, T[0], T[1], &x2[0], &x2[1] ); T-=step;
204
205 x1-=8; x2-=8;
206 }while(T>sincos_lookup0);
207 do{
208 r0 = x2[6] - x1[6]; x1[6] += x2[6];
209 r1 = x2[7] - x1[7]; x1[7] += x2[7];
210 r2 = x2[4] - x1[4]; x1[4] += x2[4];
211 r3 = x2[5] - x1[5]; x1[5] += x2[5];
212 XPROD31( r0, r1, T[0], T[1], &x2[6], &x2[7] ); T+=step;
213 XPROD31( r2, r3, T[0], T[1], &x2[4], &x2[5] ); T+=step;
214
215 r0 = x2[2] - x1[2]; x1[2] += x2[2];
216 r1 = x2[3] - x1[3]; x1[3] += x2[3];
217 r2 = x2[0] - x1[0]; x1[0] += x2[0];
218 r3 = x2[1] - x1[1]; x1[1] += x2[1];
219 XPROD31( r0, r1, T[0], T[1], &x2[2], &x2[3] ); T+=step;
220 XPROD31( r2, r3, T[0], T[1], &x2[0], &x2[1] ); T+=step;
221
222 x1-=8; x2-=8;
223 }while(T<sincos_lookup0+1024);
224 do{
225 r0 = x1[6] - x2[6]; x1[6] += x2[6];
226 r1 = x2[7] - x1[7]; x1[7] += x2[7];
227 r2 = x1[4] - x2[4]; x1[4] += x2[4];
228 r3 = x2[5] - x1[5]; x1[5] += x2[5];
229 XNPROD31( r1, r0, T[0], T[1], &x2[6], &x2[7] ); T-=step;
230 XNPROD31( r3, r2, T[0], T[1], &x2[4], &x2[5] ); T-=step;
231
232 r0 = x1[2] - x2[2]; x1[2] += x2[2];
233 r1 = x2[3] - x1[3]; x1[3] += x2[3];
234 r2 = x1[0] - x2[0]; x1[0] += x2[0];
235 r3 = x2[1] - x1[1]; x1[1] += x2[1];
236 XNPROD31( r1, r0, T[0], T[1], &x2[2], &x2[3] ); T-=step;
237 XNPROD31( r3, r2, T[0], T[1], &x2[0], &x2[1] ); T-=step;
238
239 x1-=8; x2-=8;
240 }while(T>sincos_lookup0);
241}
242
243#endif /* CPU_ARM */
244
245static inline void mdct_butterflies(ogg_int32_t *x,int points,int shift) {
246
247 int stages=8-shift;
248 int i,j;
249
250 for(i=0;--stages>0;i++){
251 for(j=0;j<(1<<i);j++)
252 mdct_butterfly_generic(x+(points>>i)*j,points>>i,4<<(i+shift));
253 }
254
255 for(j=0;j<points;j+=32)
256 mdct_butterfly_32(x+j);
257}
258
259
260static const unsigned char bitrev[16] ICONST_ATTR =
261 {0,8,4,12,2,10,6,14,1,9,5,13,3,11,7,15};
262
263static inline int bitrev12(int x){
264 return bitrev[x>>8]|(bitrev[(x&0x0f0)>>4]<<4)|(((int)bitrev[x&0x00f])<<8);
265}
266
267static inline void mdct_bitreverse(ogg_int32_t *x,int n,int step,int shift) {
268
269 int bit = 0;
270 ogg_int32_t *w0 = x;
271 ogg_int32_t *w1 = x = w0+(n>>1);
272 LOOKUP_T *T = (step>=4)?(sincos_lookup0+(step>>1)):sincos_lookup1;
273 LOOKUP_T *Ttop = T+1024;
274 register ogg_int32_t r2;
275
276 do{
277 register ogg_int32_t r3 = bitrev12(bit++);
278 ogg_int32_t *x0 = x + ((r3 ^ 0xfff)>>shift) -1;
279 ogg_int32_t *x1 = x + (r3>>shift);
280
281 register ogg_int32_t r0 = x0[0] + x1[0];
282 register ogg_int32_t r1 = x1[1] - x0[1];
283
284 XPROD32( r0, r1, T[1], T[0], r2, r3 ); T+=step;
285
286 w1 -= 4;
287
288 r0 = (x0[1] + x1[1])>>1;
289 r1 = (x0[0] - x1[0])>>1;
290 w0[0] = r0 + r2;
291 w0[1] = r1 + r3;
292 w1[2] = r0 - r2;
293 w1[3] = r3 - r1;
294
295 r3 = bitrev12(bit++);
296 x0 = x + ((r3 ^ 0xfff)>>shift) -1;
297 x1 = x + (r3>>shift);
298
299 r0 = x0[0] + x1[0];
300 r1 = x1[1] - x0[1];
301
302 XPROD32( r0, r1, T[1], T[0], r2, r3 ); T+=step;
303
304 r0 = (x0[1] + x1[1])>>1;
305 r1 = (x0[0] - x1[0])>>1;
306 w0[2] = r0 + r2;
307 w0[3] = r1 + r3;
308 w1[0] = r0 - r2;
309 w1[1] = r3 - r1;
310
311 w0 += 4;
312 }while(T<Ttop);
313 do{
314 register ogg_int32_t r3 = bitrev12(bit++);
315 ogg_int32_t *x0 = x + ((r3 ^ 0xfff)>>shift) -1;
316 ogg_int32_t *x1 = x + (r3>>shift);
317
318 register ogg_int32_t r0 = x0[0] + x1[0];
319 register ogg_int32_t r1 = x1[1] - x0[1];
320
321 T-=step; XPROD32( r0, r1, T[0], T[1], r2, r3 );
322
323 w1 -= 4;
324
325 r0 = (x0[1] + x1[1])>>1;
326 r1 = (x0[0] - x1[0])>>1;
327 w0[0] = r0 + r2;
328 w0[1] = r1 + r3;
329 w1[2] = r0 - r2;
330 w1[3] = r3 - r1;
331
332 r3 = bitrev12(bit++);
333 x0 = x + ((r3 ^ 0xfff)>>shift) -1;
334 x1 = x + (r3>>shift);
335
336 r0 = x0[0] + x1[0];
337 r1 = x1[1] - x0[1];
338
339 T-=step; XPROD32( r0, r1, T[0], T[1], r2, r3 );
340
341 r0 = (x0[1] + x1[1])>>1;
342 r1 = (x0[0] - x1[0])>>1;
343 w0[2] = r0 + r2;
344 w0[3] = r1 + r3;
345 w1[0] = r0 - r2;
346 w1[1] = r3 - r1;
347
348 w0 += 4;
349 }while(w0<w1);
350}
351
352
353void mdct_backward(int n, ogg_int32_t *in, ogg_int32_t *out)
354 ICODE_ATTR_TREMOR_MDCT;
355void mdct_backward(int n, ogg_int32_t *in, ogg_int32_t *out) {
356 int n2=n>>1;
357 int n4=n>>2;
358 ogg_int32_t *iX;
359 ogg_int32_t *oX;
360 LOOKUP_T *T;
361 LOOKUP_T *V;
362 int shift;
363 int step;
364 for (shift=6;!(n&(1<<shift));shift++);
365 shift=13-shift;
366 step=2<<shift;
367
368 /* rotate */
369
370 iX = in+n2-7;
371 oX = out+n2+n4;
372 T = sincos_lookup0;
373
374 do{
375 oX-=4;
376 XPROD31( iX[4], iX[6], T[0], T[1], &oX[2], &oX[3] ); T+=step;
377 XPROD31( iX[0], iX[2], T[0], T[1], &oX[0], &oX[1] ); T+=step;
378 iX-=8;
379 }while(iX>=in+n4);
380 do{
381 oX-=4;
382 XPROD31( iX[4], iX[6], T[1], T[0], &oX[2], &oX[3] ); T-=step;
383 XPROD31( iX[0], iX[2], T[1], T[0], &oX[0], &oX[1] ); T-=step;
384 iX-=8;
385 }while(iX>=in);
386
387 iX = in+n2-8;
388 oX = out+n2+n4;
389 T = sincos_lookup0;
390
391 do{
392 T+=step; XNPROD31( iX[6], iX[4], T[0], T[1], &oX[0], &oX[1] );
393 T+=step; XNPROD31( iX[2], iX[0], T[0], T[1], &oX[2], &oX[3] );
394 iX-=8;
395 oX+=4;
396 }while(iX>=in+n4);
397 do{
398 T-=step; XNPROD31( iX[6], iX[4], T[1], T[0], &oX[0], &oX[1] );
399 T-=step; XNPROD31( iX[2], iX[0], T[1], T[0], &oX[2], &oX[3] );
400 iX-=8;
401 oX+=4;
402 }while(iX>=in);
403
404 mdct_butterflies(out+n2,n2,shift);
405 mdct_bitreverse(out,n,step,shift);
406 /* rotate + window */
407
408 step>>=2;
409 {
410 ogg_int32_t *oX1=out+n2+n4;
411 ogg_int32_t *oX2=out+n2+n4;
412 ogg_int32_t *iX =out;
413
414 switch(step) {
415 default: {
416 T=(step>=4)?(sincos_lookup0+(step>>1)):sincos_lookup1;
417 do{
418 oX1-=4;
419 XPROD31( iX[0], -iX[1], T[0], T[1], &oX1[3], &oX2[0] ); T+=step;
420 XPROD31( iX[2], -iX[3], T[0], T[1], &oX1[2], &oX2[1] ); T+=step;
421 XPROD31( iX[4], -iX[5], T[0], T[1], &oX1[1], &oX2[2] ); T+=step;
422 XPROD31( iX[6], -iX[7], T[0], T[1], &oX1[0], &oX2[3] ); T+=step;
423 oX2+=4;
424 iX+=8;
425 }while(iX<oX1);
426 break;
427 }
428
429 case 1: {
430 /* linear interpolation between table values: offset=0.5, step=1 */
431 register ogg_int32_t t0,t1,v0,v1;
432 T = sincos_lookup0;
433 V = sincos_lookup1;
434 t0 = (*T++)>>1;
435 t1 = (*T++)>>1;
436 do{
437 oX1-=4;
438
439 t0 += (v0 = (*V++)>>1);
440 t1 += (v1 = (*V++)>>1);
441 XPROD31( iX[0], -iX[1], t0, t1, &oX1[3], &oX2[0] );
442 v0 += (t0 = (*T++)>>1);
443 v1 += (t1 = (*T++)>>1);
444 XPROD31( iX[2], -iX[3], v0, v1, &oX1[2], &oX2[1] );
445 t0 += (v0 = (*V++)>>1);
446 t1 += (v1 = (*V++)>>1);
447 XPROD31( iX[4], -iX[5], t0, t1, &oX1[1], &oX2[2] );
448 v0 += (t0 = (*T++)>>1);
449 v1 += (t1 = (*T++)>>1);
450 XPROD31( iX[6], -iX[7], v0, v1, &oX1[0], &oX2[3] );
451
452 oX2+=4;
453 iX+=8;
454 }while(iX<oX1);
455 break;
456 }
457
458 case 0: {
459 /* linear interpolation between table values: offset=0.25, step=0.5 */
460 register ogg_int32_t t0,t1,v0,v1,q0,q1;
461 T = sincos_lookup0;
462 V = sincos_lookup1;
463 t0 = *T++;
464 t1 = *T++;
465 do{
466 oX1-=4;
467
468 v0 = *V++;
469 v1 = *V++;
470 t0 += (q0 = (v0-t0)>>2);
471 t1 += (q1 = (v1-t1)>>2);
472 XPROD31( iX[0], -iX[1], t0, t1, &oX1[3], &oX2[0] );
473 t0 = v0-q0;
474 t1 = v1-q1;
475 XPROD31( iX[2], -iX[3], t0, t1, &oX1[2], &oX2[1] );
476
477 t0 = *T++;
478 t1 = *T++;
479 v0 += (q0 = (t0-v0)>>2);
480 v1 += (q1 = (t1-v1)>>2);
481 XPROD31( iX[4], -iX[5], v0, v1, &oX1[1], &oX2[2] );
482 v0 = t0-q0;
483 v1 = t1-q1;
484 XPROD31( iX[6], -iX[7], v0, v1, &oX1[0], &oX2[3] );
485
486 oX2+=4;
487 iX+=8;
488 }while(iX<oX1);
489 break;
490 }
491 }
492
493 iX=out+n2+n4;
494 oX1=out+n4;
495 oX2=oX1;
496
497 do{
498 oX1-=4;
499 iX-=4;
500
501 oX2[0] = -(oX1[3] = iX[3]);
502 oX2[1] = -(oX1[2] = iX[2]);
503 oX2[2] = -(oX1[1] = iX[1]);
504 oX2[3] = -(oX1[0] = iX[0]);
505
506 oX2+=4;
507 }while(oX2<iX);
508
509 iX=out+n2+n4;
510 oX1=out+n2+n4;
511 oX2=out+n2;
512
513 do{
514 oX1-=4;
515 oX1[0]= iX[3];
516 oX1[1]= iX[2];
517 oX1[2]= iX[1];
518 oX1[3]= iX[0];
519 iX+=4;
520 }while(oX1>oX2);
521 }
522}