summaryrefslogtreecommitdiff
path: root/apps/codecs/lib/mdct2.c
diff options
context:
space:
mode:
Diffstat (limited to 'apps/codecs/lib/mdct2.c')
-rw-r--r--apps/codecs/lib/mdct2.c513
1 files changed, 0 insertions, 513 deletions
diff --git a/apps/codecs/lib/mdct2.c b/apps/codecs/lib/mdct2.c
deleted file mode 100644
index ba8b5ca6be..0000000000
--- a/apps/codecs/lib/mdct2.c
+++ /dev/null
@@ -1,513 +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#ifdef ROCKBOX
41#include <codecs/lib/codeclib.h>
42#endif /* ROCKBOX */
43
44#if defined(CPU_ARM)
45
46extern void mdct_butterfly_32(int32_t *x);
47extern void mdct_butterfly_generic_loop(int32_t *x1, int32_t *x2,
48 const int32_t *T0, int step,
49 const int32_t *Ttop);
50
51static inline void mdct_butterfly_generic(int32_t *x,int points, int step){
52 mdct_butterfly_generic_loop(x + points, x + (points>>1), sincos_lookup0, step, sincos_lookup0+1024);
53}
54
55#else
56
57/* 8 point butterfly (in place) */
58static inline void mdct_butterfly_8(int32_t *x){
59 register int32_t r0 = x[4] + x[0];
60 register int32_t r1 = x[4] - x[0];
61 register int32_t r2 = x[5] + x[1];
62 register int32_t r3 = x[5] - x[1];
63 register int32_t r4 = x[6] + x[2];
64 register int32_t r5 = x[6] - x[2];
65 register int32_t r6 = x[7] + x[3];
66 register 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}
77
78/* 16 point butterfly (in place, 4 register) */
79static inline void mdct_butterfly_16(int32_t *x){
80
81 register int32_t r0, r1;
82
83 r0 = x[ 0] - x[ 8]; x[ 8] += x[ 0];
84 r1 = x[ 1] - x[ 9]; x[ 9] += x[ 1];
85 x[ 0] = MULT31((r0 + r1) , cPI2_8);
86 x[ 1] = MULT31((r1 - r0) , cPI2_8);
87
88 r0 = x[10] - x[ 2]; x[10] += x[ 2];
89 r1 = x[ 3] - x[11]; x[11] += x[ 3];
90 x[ 2] = r1; x[ 3] = r0;
91
92 r0 = x[12] - x[ 4]; x[12] += x[ 4];
93 r1 = x[13] - x[ 5]; x[13] += x[ 5];
94 x[ 4] = MULT31((r0 - r1) , cPI2_8);
95 x[ 5] = MULT31((r0 + r1) , cPI2_8);
96
97 r0 = x[14] - x[ 6]; x[14] += x[ 6];
98 r1 = x[15] - x[ 7]; x[15] += x[ 7];
99 x[ 6] = r0; x[ 7] = r1;
100
101 mdct_butterfly_8(x);
102 mdct_butterfly_8(x+8);
103}
104
105/* 32 point butterfly (in place, 4 register) */
106static inline void mdct_butterfly_32(int32_t *x){
107
108 register int32_t r0, r1;
109
110 r0 = x[30] - x[14]; x[30] += x[14];
111 r1 = x[31] - x[15]; x[31] += x[15];
112 x[14] = r0; x[15] = r1;
113
114 r0 = x[28] - x[12]; x[28] += x[12];
115 r1 = x[29] - x[13]; x[29] += x[13];
116 XNPROD31( r0, r1, cPI1_8, cPI3_8, &x[12], &x[13] );
117
118 r0 = x[26] - x[10]; x[26] += x[10];
119 r1 = x[27] - x[11]; x[27] += x[11];
120 x[10] = MULT31((r0 - r1) , cPI2_8);
121 x[11] = MULT31((r0 + r1) , cPI2_8);
122
123 r0 = x[24] - x[ 8]; x[24] += x[ 8];
124 r1 = x[25] - x[ 9]; x[25] += x[ 9];
125 XNPROD31( r0, r1, cPI3_8, cPI1_8, &x[ 8], &x[ 9] );
126
127 r0 = x[22] - x[ 6]; x[22] += x[ 6];
128 r1 = x[ 7] - x[23]; x[23] += x[ 7];
129 x[ 6] = r1; x[ 7] = r0;
130
131 r0 = x[ 4] - x[20]; x[20] += x[ 4];
132 r1 = x[ 5] - x[21]; x[21] += x[ 5];
133 XPROD31 ( r0, r1, cPI3_8, cPI1_8, &x[ 4], &x[ 5] );
134
135 r0 = x[ 2] - x[18]; x[18] += x[ 2];
136 r1 = x[ 3] - x[19]; x[19] += x[ 3];
137 x[ 2] = MULT31((r1 + r0) , cPI2_8);
138 x[ 3] = MULT31((r1 - r0) , cPI2_8);
139
140 r0 = x[ 0] - x[16]; x[16] += x[ 0];
141 r1 = x[ 1] - x[17]; x[17] += x[ 1];
142 XPROD31 ( r0, r1, cPI1_8, cPI3_8, &x[ 0], &x[ 1] );
143
144 mdct_butterfly_16(x);
145 mdct_butterfly_16(x+16);
146}
147
148/* N/stage point generic N stage butterfly (in place, 4 register) */
149void mdct_butterfly_generic(int32_t *x,int points, int step)
150 ICODE_ATTR_TREMOR_MDCT;
151void mdct_butterfly_generic(int32_t *x,int points, int step){
152 const int32_t *T = sincos_lookup0;
153 int32_t *x1 = x + points - 8;
154 int32_t *x2 = x + (points>>1) - 8;
155 register int32_t r0;
156 register int32_t r1;
157 register int32_t r2;
158 register int32_t r3;
159
160 do{
161 r0 = x1[6] - x2[6]; x1[6] += x2[6];
162 r1 = x2[7] - x1[7]; x1[7] += x2[7];
163 r2 = x1[4] - x2[4]; x1[4] += x2[4];
164 r3 = x2[5] - x1[5]; x1[5] += x2[5];
165 XPROD31( r1, r0, T[0], T[1], &x2[6], &x2[7] ); T+=step;
166 XPROD31( r3, r2, T[0], T[1], &x2[4], &x2[5] ); T+=step;
167
168 r0 = x1[2] - x2[2]; x1[2] += x2[2];
169 r1 = x2[3] - x1[3]; x1[3] += x2[3];
170 r2 = x1[0] - x2[0]; x1[0] += x2[0];
171 r3 = x2[1] - x1[1]; x1[1] += x2[1];
172 XPROD31( r1, r0, T[0], T[1], &x2[2], &x2[3] ); T+=step;
173 XPROD31( r3, r2, T[0], T[1], &x2[0], &x2[1] ); T+=step;
174
175 x1-=8; x2-=8;
176 }while(T<sincos_lookup0+1024);
177 do{
178 r0 = x1[6] - x2[6]; x1[6] += x2[6];
179 r1 = x1[7] - x2[7]; x1[7] += x2[7];
180 r2 = x1[4] - x2[4]; x1[4] += x2[4];
181 r3 = x1[5] - x2[5]; x1[5] += x2[5];
182 XNPROD31( r0, r1, T[0], T[1], &x2[6], &x2[7] ); T-=step;
183 XNPROD31( r2, r3, T[0], T[1], &x2[4], &x2[5] ); T-=step;
184
185 r0 = x1[2] - x2[2]; x1[2] += x2[2];
186 r1 = x1[3] - x2[3]; x1[3] += x2[3];
187 r2 = x1[0] - x2[0]; x1[0] += x2[0];
188 r3 = x1[1] - x2[1]; x1[1] += x2[1];
189 XNPROD31( r0, r1, T[0], T[1], &x2[2], &x2[3] ); T-=step;
190 XNPROD31( r2, r3, T[0], T[1], &x2[0], &x2[1] ); T-=step;
191
192 x1-=8; x2-=8;
193 }while(T>sincos_lookup0);
194 do{
195 r0 = x2[6] - x1[6]; x1[6] += x2[6];
196 r1 = x2[7] - x1[7]; x1[7] += x2[7];
197 r2 = x2[4] - x1[4]; x1[4] += x2[4];
198 r3 = x2[5] - x1[5]; x1[5] += x2[5];
199 XPROD31( r0, r1, T[0], T[1], &x2[6], &x2[7] ); T+=step;
200 XPROD31( r2, r3, T[0], T[1], &x2[4], &x2[5] ); T+=step;
201
202 r0 = x2[2] - x1[2]; x1[2] += x2[2];
203 r1 = x2[3] - x1[3]; x1[3] += x2[3];
204 r2 = x2[0] - x1[0]; x1[0] += x2[0];
205 r3 = x2[1] - x1[1]; x1[1] += x2[1];
206 XPROD31( r0, r1, T[0], T[1], &x2[2], &x2[3] ); T+=step;
207 XPROD31( r2, r3, T[0], T[1], &x2[0], &x2[1] ); T+=step;
208
209 x1-=8; x2-=8;
210 }while(T<sincos_lookup0+1024);
211 do{
212 r0 = x1[6] - x2[6]; x1[6] += x2[6];
213 r1 = x2[7] - x1[7]; x1[7] += x2[7];
214 r2 = x1[4] - x2[4]; x1[4] += x2[4];
215 r3 = x2[5] - x1[5]; x1[5] += x2[5];
216 XNPROD31( r1, r0, T[0], T[1], &x2[6], &x2[7] ); T-=step;
217 XNPROD31( r3, r2, T[0], T[1], &x2[4], &x2[5] ); T-=step;
218
219 r0 = x1[2] - x2[2]; x1[2] += x2[2];
220 r1 = x2[3] - x1[3]; x1[3] += x2[3];
221 r2 = x1[0] - x2[0]; x1[0] += x2[0];
222 r3 = x2[1] - x1[1]; x1[1] += x2[1];
223 XNPROD31( r1, r0, T[0], T[1], &x2[2], &x2[3] ); T-=step;
224 XNPROD31( r3, r2, T[0], T[1], &x2[0], &x2[1] ); T-=step;
225
226 x1-=8; x2-=8;
227 }while(T>sincos_lookup0);
228}
229
230#endif /* CPU_ARM */
231
232static inline void mdct_butterflies(int32_t *x,int points,int shift) {
233
234 int stages=8-shift;
235 int i,j;
236
237 for(i=0;--stages>0;i++){
238 for(j=0;j<(1<<i);j++)
239 mdct_butterfly_generic(x+(points>>i)*j,points>>i,4<<(i+shift));
240 }
241
242 for(j=0;j<points;j+=32)
243 mdct_butterfly_32(x+j);
244}
245
246static const unsigned char bitrev[] ICONST_ATTR =
247{
248 0, 32, 16, 48, 8, 40, 24, 56, 4, 36, 20, 52, 12, 44, 28, 60,
249 2, 34, 18, 50, 10, 42, 26, 58, 6, 38, 22, 54, 14, 46, 30, 62,
250 1, 33, 17, 49, 9, 41, 25, 57, 5, 37, 21, 53, 13, 45, 29, 61,
251 3, 35, 19, 51, 11, 43, 27, 59, 7, 39, 23, 55, 15, 47, 31, 63
252};
253
254static inline int bitrev12(int x){
255 return bitrev[x>>6]|((bitrev[x&0x03f])<<6);
256}
257
258static inline void mdct_bitreverse(int32_t *x,int n,int step,int shift) {
259
260 int bit = 0;
261 int32_t *w0 = x;
262 int32_t *w1 = x = w0+(n>>1);
263 const int32_t *T = (step>=4)?(sincos_lookup0+(step>>1)):sincos_lookup1;
264 const int32_t *Ttop = T+1024;
265 register int32_t r2;
266
267 do{
268 register int32_t r3 = bitrev12(bit++);
269 int32_t *x0 = x + ((r3 ^ 0xfff)>>shift) -1;
270 int32_t *x1 = x + (r3>>shift);
271
272 register int32_t r0 = x0[0] + x1[0];
273 register int32_t r1 = x1[1] - x0[1];
274
275 XPROD32( r0, r1, T[1], T[0], r2, r3 ); T+=step;
276
277 w1 -= 4;
278
279 r0 = (x0[1] + x1[1])>>1;
280 r1 = (x0[0] - x1[0])>>1;
281 w0[0] = r0 + r2;
282 w0[1] = r1 + r3;
283 w1[2] = r0 - r2;
284 w1[3] = r3 - r1;
285
286 r3 = bitrev12(bit++);
287 x0 = x + ((r3 ^ 0xfff)>>shift) -1;
288 x1 = x + (r3>>shift);
289
290 r0 = x0[0] + x1[0];
291 r1 = x1[1] - x0[1];
292
293 XPROD32( r0, r1, T[1], T[0], r2, r3 ); T+=step;
294
295 r0 = (x0[1] + x1[1])>>1;
296 r1 = (x0[0] - x1[0])>>1;
297 w0[2] = r0 + r2;
298 w0[3] = r1 + r3;
299 w1[0] = r0 - r2;
300 w1[1] = r3 - r1;
301
302 w0 += 4;
303 }while(T<Ttop);
304 do{
305 register int32_t r3 = bitrev12(bit++);
306 int32_t *x0 = x + ((r3 ^ 0xfff)>>shift) -1;
307 int32_t *x1 = x + (r3>>shift);
308
309 register int32_t r0 = x0[0] + x1[0];
310 register int32_t r1 = x1[1] - x0[1];
311
312 T-=step; XPROD32( r0, r1, T[0], T[1], r2, r3 );
313
314 w1 -= 4;
315
316 r0 = (x0[1] + x1[1])>>1;
317 r1 = (x0[0] - x1[0])>>1;
318 w0[0] = r0 + r2;
319 w0[1] = r1 + r3;
320 w1[2] = r0 - r2;
321 w1[3] = r3 - r1;
322
323 r3 = bitrev12(bit++);
324 x0 = x + ((r3 ^ 0xfff)>>shift) -1;
325 x1 = x + (r3>>shift);
326
327 r0 = x0[0] + x1[0];
328 r1 = x1[1] - x0[1];
329
330 T-=step; XPROD32( r0, r1, T[0], T[1], r2, r3 );
331
332 r0 = (x0[1] + x1[1])>>1;
333 r1 = (x0[0] - x1[0])>>1;
334 w0[2] = r0 + r2;
335 w0[3] = r1 + r3;
336 w1[0] = r0 - r2;
337 w1[1] = r3 - r1;
338
339 w0 += 4;
340 }while(w0<w1);
341}
342
343
344void mdct_backward(int n, int32_t *in, int32_t *out)
345 ICODE_ATTR_TREMOR_MDCT;
346void mdct_backward(int n, int32_t *in, int32_t *out) {
347 int n2=n>>1;
348 int n4=n>>2;
349 int32_t *iX;
350 int32_t *oX;
351 const int32_t *T;
352 const int32_t *V;
353 int shift;
354 int step;
355 for (shift=6;!(n&(1<<shift));shift++);
356 shift=13-shift;
357 step=2<<shift;
358
359 /* rotate */
360
361 iX = in+n2-7;
362 oX = out+n2+n4;
363 T = sincos_lookup0;
364
365 do{
366 oX-=4;
367 XPROD31( iX[4], iX[6], T[0], T[1], &oX[2], &oX[3] ); T+=step;
368 XPROD31( iX[0], iX[2], T[0], T[1], &oX[0], &oX[1] ); T+=step;
369 iX-=8;
370 }while(iX>=in+n4);
371 do{
372 oX-=4;
373 XPROD31( iX[4], iX[6], T[1], T[0], &oX[2], &oX[3] ); T-=step;
374 XPROD31( iX[0], iX[2], T[1], T[0], &oX[0], &oX[1] ); T-=step;
375 iX-=8;
376 }while(iX>=in);
377
378 iX = in+n2-8;
379 oX = out+n2+n4;
380 T = sincos_lookup0;
381
382 do{
383 T+=step; XNPROD31( iX[6], iX[4], T[0], T[1], &oX[0], &oX[1] );
384 T+=step; XNPROD31( iX[2], iX[0], T[0], T[1], &oX[2], &oX[3] );
385 iX-=8;
386 oX+=4;
387 }while(iX>=in+n4);
388 do{
389 T-=step; XNPROD31( iX[6], iX[4], T[1], T[0], &oX[0], &oX[1] );
390 T-=step; XNPROD31( iX[2], iX[0], T[1], T[0], &oX[2], &oX[3] );
391 iX-=8;
392 oX+=4;
393 }while(iX>=in);
394
395 mdct_butterflies(out+n2,n2,shift);
396 mdct_bitreverse(out,n,step,shift);
397 /* rotate + window */
398
399 step>>=2;
400 {
401 int32_t *oX1=out+n2+n4;
402 int32_t *oX2=out+n2+n4;
403 int32_t *iX =out;
404
405 switch(step) {
406 default: {
407 T=(step>=4)?(sincos_lookup0+(step>>1)):sincos_lookup1;
408 do{
409 oX1-=4;
410 XPROD31( iX[0], -iX[1], T[0], T[1], &oX1[3], &oX2[0] ); T+=step;
411 XPROD31( iX[2], -iX[3], T[0], T[1], &oX1[2], &oX2[1] ); T+=step;
412 XPROD31( iX[4], -iX[5], T[0], T[1], &oX1[1], &oX2[2] ); T+=step;
413 XPROD31( iX[6], -iX[7], T[0], T[1], &oX1[0], &oX2[3] ); T+=step;
414 oX2+=4;
415 iX+=8;
416 }while(iX<oX1);
417 break;
418 }
419
420 case 1: {
421 /* linear interpolation between table values: offset=0.5, step=1 */
422 register int32_t t0,t1,v0,v1;
423 T = sincos_lookup0;
424 V = sincos_lookup1;
425 t0 = (*T++)>>1;
426 t1 = (*T++)>>1;
427 do{
428 oX1-=4;
429
430 t0 += (v0 = (*V++)>>1);
431 t1 += (v1 = (*V++)>>1);
432 XPROD31( iX[0], -iX[1], t0, t1, &oX1[3], &oX2[0] );
433 v0 += (t0 = (*T++)>>1);
434 v1 += (t1 = (*T++)>>1);
435 XPROD31( iX[2], -iX[3], v0, v1, &oX1[2], &oX2[1] );
436 t0 += (v0 = (*V++)>>1);
437 t1 += (v1 = (*V++)>>1);
438 XPROD31( iX[4], -iX[5], t0, t1, &oX1[1], &oX2[2] );
439 v0 += (t0 = (*T++)>>1);
440 v1 += (t1 = (*T++)>>1);
441 XPROD31( iX[6], -iX[7], v0, v1, &oX1[0], &oX2[3] );
442
443 oX2+=4;
444 iX+=8;
445 }while(iX<oX1);
446 break;
447 }
448
449 case 0: {
450 /* linear interpolation between table values: offset=0.25, step=0.5 */
451 register int32_t t0,t1,v0,v1,q0,q1;
452 T = sincos_lookup0;
453 V = sincos_lookup1;
454 t0 = *T++;
455 t1 = *T++;
456 do{
457 oX1-=4;
458
459 v0 = *V++;
460 v1 = *V++;
461 t0 += (q0 = (v0-t0)>>2);
462 t1 += (q1 = (v1-t1)>>2);
463 XPROD31( iX[0], -iX[1], t0, t1, &oX1[3], &oX2[0] );
464 t0 = v0-q0;
465 t1 = v1-q1;
466 XPROD31( iX[2], -iX[3], t0, t1, &oX1[2], &oX2[1] );
467
468 t0 = *T++;
469 t1 = *T++;
470 v0 += (q0 = (t0-v0)>>2);
471 v1 += (q1 = (t1-v1)>>2);
472 XPROD31( iX[4], -iX[5], v0, v1, &oX1[1], &oX2[2] );
473 v0 = t0-q0;
474 v1 = t1-q1;
475 XPROD31( iX[6], -iX[7], v0, v1, &oX1[0], &oX2[3] );
476
477 oX2+=4;
478 iX+=8;
479 }while(iX<oX1);
480 break;
481 }
482 }
483
484 iX=out+n2+n4;
485 oX1=out+n4;
486 oX2=oX1;
487
488 do{
489 oX1-=4;
490 iX-=4;
491
492 oX2[0] = -(oX1[3] = iX[3]);
493 oX2[1] = -(oX1[2] = iX[2]);
494 oX2[2] = -(oX1[1] = iX[1]);
495 oX2[3] = -(oX1[0] = iX[0]);
496
497 oX2+=4;
498 }while(oX2<iX);
499
500 iX=out+n2+n4;
501 oX1=out+n2+n4;
502 oX2=out+n2;
503
504 do{
505 oX1-=4;
506 oX1[0]= iX[3];
507 oX1[1]= iX[2];
508 oX1[2]= iX[1];
509 oX1[3]= iX[0];
510 iX+=4;
511 }while(oX1>oX2);
512 }
513}
generated by cgit v1.2.3 (git 2.39.1) at 2024-09-12 18:36:13 +0000