apps/codecs/lib/mdct2.c

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