lib/mdct.c

   1 /********************************************************************
   2  *                                                                  *
   3  * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE.   *
   4  * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS     *
   5  * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
   6  * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.       *
   7  *                                                                  *
   8  * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2009             *
   9  * by the Xiph.Org Foundation http://www.xiph.org/                  *
  10  *                                                                  *
  11  ********************************************************************
  12
  13  function: normalized modified discrete cosine transform
  14            power of two length transform only [64 <= n ]
  15
  16  Original algorithm adapted long ago from _The use of multirate filter
  17  banks for coding of high quality digital audio_, by T. Sporer,
  18  K. Brandenburg and B. Edler, collection of the European Signal
  19  Processing Conference (EUSIPCO), Amsterdam, June 1992, Vol.1, pp
  20  211-214
  21
  22  The below code implements an algorithm that no longer looks much like
  23  that presented in the paper, but the basic structure remains if you
  24  dig deep enough to see it.
  25
  26  This module DOES NOT INCLUDE code to generate/apply the window
  27  function.  Everybody has their own weird favorite including me... I
  28  happen to like the properties of y=sin(.5PI*sin^2(x)), but others may
  29  vehemently disagree.
  30
  31  ********************************************************************/
  32
  33 /* this can also be run as an integer transform by uncommenting a
  34    define in mdct.h; the integerization is a first pass and although
  35    it's likely stable for Vorbis, the dynamic range is constrained and
  36    roundoff isn't done (so it's noisy).  Consider it functional, but
  37    only a starting point.  There's no point on a machine with an FPU */
  38
  39 #include <stdio.h>
  40 #include <stdlib.h>
  41 #include <string.h>
  42 #include <math.h>
  43 #include "vorbis/codec.h"
  44 #include "mdct.h"
  45 #include "os.h"
  46 #include "misc.h"
  47
  48 /* build lookups for trig functions; also pre-figure scaling and
  49    some window function algebra. */
  50
  51 void mdct_init(mdct_lookup *lookup,int n){
  52   int   *bitrev=_ogg_malloc(sizeof(*bitrev)*(n/4));
  53   DATA_TYPE *T=_ogg_malloc(sizeof(*T)*(n+n/4));
  54
  55   int i;
  56   int n2=n>>1;
  57   int log2n=lookup->log2n=rint(log((float)n)/log(2.f));
  58   lookup->n=n;
  59   lookup->trig=T;
  60   lookup->bitrev=bitrev;
  61
  62 /* trig lookups... */
  63
  64   for(i=0;i<n/4;i++){
  65     T[i*2]=FLOAT_CONV(cos((M_PI/n)*(4*i)));
  66     T[i*2+1]=FLOAT_CONV(-sin((M_PI/n)*(4*i)));
  67     T[n2+i*2]=FLOAT_CONV(cos((M_PI/(2*n))*(2*i+1)));
  68     T[n2+i*2+1]=FLOAT_CONV(sin((M_PI/(2*n))*(2*i+1)));
  69   }
  70   for(i=0;i<n/8;i++){
  71     T[n+i*2]=FLOAT_CONV(cos((M_PI/n)*(4*i+2))*.5);
  72     T[n+i*2+1]=FLOAT_CONV(-sin((M_PI/n)*(4*i+2))*.5);
  73   }
  74
  75   /* bitreverse lookup... */
  76
  77   {
  78     int mask=(1<<(log2n-1))-1,i,j;
  79     int msb=1<<(log2n-2);
  80     for(i=0;i<n/8;i++){
  81       int acc=0;
  82       for(j=0;msb>>j;j++)
  83         if((msb>>j)&i)acc|=1<<j;
  84       bitrev[i*2]=((~acc)&mask)-1;
  85       bitrev[i*2+1]=acc;
  86
  87     }
  88   }
  89   lookup->scale=FLOAT_CONV(4.f/n);
  90 }
  91
  92 /* 8 point butterfly (in place, 4 register) */
  93 STIN void mdct_butterfly_8(DATA_TYPE *x){
  94   REG_TYPE r0   = x[6] + x[2];
  95   REG_TYPE r1   = x[6] - x[2];
  96   REG_TYPE r2   = x[4] + x[0];
  97   REG_TYPE r3   = x[4] - x[0];
  98
  99            x[6] = r0   + r2;
 100            x[4] = r0   - r2;
 101
 102            r0   = x[5] - x[1];
 103            r2   = x[7] - x[3];
 104            x[0] = r1   + r0;
 105            x[2] = r1   - r0;
 106
 107            r0   = x[5] + x[1];
 108            r1   = x[7] + x[3];
 109            x[3] = r2   + r3;
 110            x[1] = r2   - r3;
 111            x[7] = r1   + r0;
 112            x[5] = r1   - r0;
 113
 114 }
 115
 116 /* 16 point butterfly (in place, 4 register) */
 117 STIN void mdct_butterfly_16(DATA_TYPE *x){
 118   REG_TYPE r0     = x[1]  - x[9];
 119   REG_TYPE r1     = x[0]  - x[8];
 120
 121            x[8]  += x[0];
 122            x[9]  += x[1];
 123            x[0]   = MULT_NORM((r0   + r1) * cPI2_8);
 124            x[1]   = MULT_NORM((r0   - r1) * cPI2_8);
 125
 126            r0     = x[3]  - x[11];
 127            r1     = x[10] - x[2];
 128            x[10] += x[2];
 129            x[11] += x[3];
 130            x[2]   = r0;
 131            x[3]   = r1;
 132
 133            r0     = x[12] - x[4];
 134            r1     = x[13] - x[5];
 135            x[12] += x[4];
 136            x[13] += x[5];
 137            x[4]   = MULT_NORM((r0   - r1) * cPI2_8);
 138            x[5]   = MULT_NORM((r0   + r1) * cPI2_8);
 139
 140            r0     = x[14] - x[6];
 141            r1     = x[15] - x[7];
 142            x[14] += x[6];
 143            x[15] += x[7];
 144            x[6]  = r0;
 145            x[7]  = r1;
 146
 147            mdct_butterfly_8(x);
 148            mdct_butterfly_8(x+8);
 149 }
 150
 151 /* 32 point butterfly (in place, 4 register) */
 152 STIN void mdct_butterfly_32(DATA_TYPE *x){
 153   REG_TYPE r0     = x[30] - x[14];
 154   REG_TYPE r1     = x[31] - x[15];
 155
 156            x[30] +=         x[14];
 157            x[31] +=         x[15];
 158            x[14]  =         r0;
 159            x[15]  =         r1;
 160
 161            r0     = x[28] - x[12];
 162            r1     = x[29] - x[13];
 163            x[28] +=         x[12];
 164            x[29] +=         x[13];
 165            x[12]  = MULT_NORM( r0 * cPI1_8  -  r1 * cPI3_8 );
 166            x[13]  = MULT_NORM( r0 * cPI3_8  +  r1 * cPI1_8 );
 167
 168            r0     = x[26] - x[10];
 169            r1     = x[27] - x[11];
 170            x[26] +=         x[10];
 171            x[27] +=         x[11];
 172            x[10]  = MULT_NORM(( r0  - r1 ) * cPI2_8);
 173            x[11]  = MULT_NORM(( r0  + r1 ) * cPI2_8);
 174
 175            r0     = x[24] - x[8];
 176            r1     = x[25] - x[9];
 177            x[24] += x[8];
 178            x[25] += x[9];
 179            x[8]   = MULT_NORM( r0 * cPI3_8  -  r1 * cPI1_8 );
 180            x[9]   = MULT_NORM( r1 * cPI3_8  +  r0 * cPI1_8 );
 181
 182            r0     = x[22] - x[6];
 183            r1     = x[7]  - x[23];
 184            x[22] += x[6];
 185            x[23] += x[7];
 186            x[6]   = r1;
 187            x[7]   = r0;
 188
 189            r0     = x[4]  - x[20];
 190            r1     = x[5]  - x[21];
 191            x[20] += x[4];
 192            x[21] += x[5];
 193            x[4]   = MULT_NORM( r1 * cPI1_8  +  r0 * cPI3_8 );
 194            x[5]   = MULT_NORM( r1 * cPI3_8  -  r0 * cPI1_8 );
 195
 196            r0     = x[2]  - x[18];
 197            r1     = x[3]  - x[19];
 198            x[18] += x[2];
 199            x[19] += x[3];
 200            x[2]   = MULT_NORM(( r1  + r0 ) * cPI2_8);
 201            x[3]   = MULT_NORM(( r1  - r0 ) * cPI2_8);
 202
 203            r0     = x[0]  - x[16];
 204            r1     = x[1]  - x[17];
 205            x[16] += x[0];
 206            x[17] += x[1];
 207            x[0]   = MULT_NORM( r1 * cPI3_8  +  r0 * cPI1_8 );
 208            x[1]   = MULT_NORM( r1 * cPI1_8  -  r0 * cPI3_8 );
 209
 210            mdct_butterfly_16(x);
 211            mdct_butterfly_16(x+16);
 212
 213 }
 214
 215 /* N point first stage butterfly (in place, 2 register) */
 216 STIN void mdct_butterfly_first(DATA_TYPE *T,
 217                                         DATA_TYPE *x,
 218                                         int points){
 219
 220   DATA_TYPE *x1        = x          + points      - 8;
 221   DATA_TYPE *x2        = x          + (points>>1) - 8;
 222   REG_TYPE   r0;
 223   REG_TYPE   r1;
 224
 225   do{
 226
 227                r0      = x1[6]      -  x2[6];
 228                r1      = x1[7]      -  x2[7];
 229                x1[6]  += x2[6];
 230                x1[7]  += x2[7];
 231                x2[6]   = MULT_NORM(r1 * T[1]  +  r0 * T[0]);
 232                x2[7]   = MULT_NORM(r1 * T[0]  -  r0 * T[1]);
 233
 234                r0      = x1[4]      -  x2[4];
 235                r1      = x1[5]      -  x2[5];
 236                x1[4]  += x2[4];
 237                x1[5]  += x2[5];
 238                x2[4]   = MULT_NORM(r1 * T[5]  +  r0 * T[4]);
 239                x2[5]   = MULT_NORM(r1 * T[4]  -  r0 * T[5]);
 240
 241                r0      = x1[2]      -  x2[2];
 242                r1      = x1[3]      -  x2[3];
 243                x1[2]  += x2[2];
 244                x1[3]  += x2[3];
 245                x2[2]   = MULT_NORM(r1 * T[9]  +  r0 * T[8]);
 246                x2[3]   = MULT_NORM(r1 * T[8]  -  r0 * T[9]);
 247
 248                r0      = x1[0]      -  x2[0];
 249                r1      = x1[1]      -  x2[1];
 250                x1[0]  += x2[0];
 251                x1[1]  += x2[1];
 252                x2[0]   = MULT_NORM(r1 * T[13] +  r0 * T[12]);
 253                x2[1]   = MULT_NORM(r1 * T[12] -  r0 * T[13]);
 254
 255     x1-=8;
 256     x2-=8;
 257     T+=16;
 258
 259   }while(x2>=x);
 260 }
 261
 262 /* N/stage point generic N stage butterfly (in place, 2 register) */
 263 STIN void mdct_butterfly_generic(DATA_TYPE *T,
 264                                           DATA_TYPE *x,
 265                                           int points,
 266                                           int trigint){
 267
 268   DATA_TYPE *x1        = x          + points      - 8;
 269   DATA_TYPE *x2        = x          + (points>>1) - 8;
 270   REG_TYPE   r0;
 271   REG_TYPE   r1;
 272
 273   do{
 274
 275                r0      = x1[6]      -  x2[6];
 276                r1      = x1[7]      -  x2[7];
 277                x1[6]  += x2[6];
 278                x1[7]  += x2[7];
 279                x2[6]   = MULT_NORM(r1 * T[1]  +  r0 * T[0]);
 280                x2[7]   = MULT_NORM(r1 * T[0]  -  r0 * T[1]);
 281
 282                T+=trigint;
 283
 284                r0      = x1[4]      -  x2[4];
 285                r1      = x1[5]      -  x2[5];
 286                x1[4]  += x2[4];
 287                x1[5]  += x2[5];
 288                x2[4]   = MULT_NORM(r1 * T[1]  +  r0 * T[0]);
 289                x2[5]   = MULT_NORM(r1 * T[0]  -  r0 * T[1]);
 290
 291                T+=trigint;
 292
 293                r0      = x1[2]      -  x2[2];
 294                r1      = x1[3]      -  x2[3];
 295                x1[2]  += x2[2];
 296                x1[3]  += x2[3];
 297                x2[2]   = MULT_NORM(r1 * T[1]  +  r0 * T[0]);
 298                x2[3]   = MULT_NORM(r1 * T[0]  -  r0 * T[1]);
 299
 300                T+=trigint;
 301
 302                r0      = x1[0]      -  x2[0];
 303                r1      = x1[1]      -  x2[1];
 304                x1[0]  += x2[0];
 305                x1[1]  += x2[1];
 306                x2[0]   = MULT_NORM(r1 * T[1]  +  r0 * T[0]);
 307                x2[1]   = MULT_NORM(r1 * T[0]  -  r0 * T[1]);
 308
 309                T+=trigint;
 310     x1-=8;
 311     x2-=8;
 312
 313   }while(x2>=x);
 314 }
 315
 316 STIN void mdct_butterflies(mdct_lookup *init,
 317                              DATA_TYPE *x,
 318                              int points){
 319
 320   DATA_TYPE *T=init->trig;
 321   int stages=init->log2n-5;
 322   int i,j;
 323
 324   if(--stages>0){
 325     mdct_butterfly_first(T,x,points);
 326   }
 327
 328   for(i=1;--stages>0;i++){
 329     for(j=0;j<(1<<i);j++)
 330       mdct_butterfly_generic(T,x+(points>>i)*j,points>>i,4<<i);
 331   }
 332
 333   for(j=0;j<points;j+=32)
 334     mdct_butterfly_32(x+j);
 335
 336 }
 337
 338 void mdct_clear(mdct_lookup *l){
 339   if(l){
 340     if(l->trig)_ogg_free(l->trig);
 341     if(l->bitrev)_ogg_free(l->bitrev);
 342     memset(l,0,sizeof(*l));
 343   }
 344 }
 345
 346 STIN void mdct_bitreverse(mdct_lookup *init,
 347                             DATA_TYPE *x){
 348   int        n       = init->n;
 349   int       *bit     = init->bitrev;
 350   DATA_TYPE *w0      = x;
 351   DATA_TYPE *w1      = x = w0+(n>>1);
 352   DATA_TYPE *T       = init->trig+n;
 353
 354   do{
 355     DATA_TYPE *x0    = x+bit[0];
 356     DATA_TYPE *x1    = x+bit[1];
 357
 358     REG_TYPE  r0     = x0[1]  - x1[1];
 359     REG_TYPE  r1     = x0[0]  + x1[0];
 360     REG_TYPE  r2     = MULT_NORM(r1     * T[0]   + r0 * T[1]);
 361     REG_TYPE  r3     = MULT_NORM(r1     * T[1]   - r0 * T[0]);
 362
 363               w1    -= 4;
 364
 365               r0     = HALVE(x0[1] + x1[1]);
 366               r1     = HALVE(x0[0] - x1[0]);
 367
 368               w0[0]  = r0     + r2;
 369               w1[2]  = r0     - r2;
 370               w0[1]  = r1     + r3;
 371               w1[3]  = r3     - r1;
 372
 373               x0     = x+bit[2];
 374               x1     = x+bit[3];
 375
 376               r0     = x0[1]  - x1[1];
 377               r1     = x0[0]  + x1[0];
 378               r2     = MULT_NORM(r1     * T[2]   + r0 * T[3]);
 379               r3     = MULT_NORM(r1     * T[3]   - r0 * T[2]);
 380
 381               r0     = HALVE(x0[1] + x1[1]);
 382               r1     = HALVE(x0[0] - x1[0]);
 383
 384               w0[2]  = r0     + r2;
 385               w1[0]  = r0     - r2;
 386               w0[3]  = r1     + r3;
 387               w1[1]  = r3     - r1;
 388
 389               T     += 4;
 390               bit   += 4;
 391               w0    += 4;
 392
 393   }while(w0<w1);
 394 }
 395
 396 void mdct_backward(mdct_lookup *init, DATA_TYPE *in, DATA_TYPE *out){
 397   int n=init->n;
 398   int n2=n>>1;
 399   int n4=n>>2;
 400
 401   /* rotate */
 402
 403   DATA_TYPE *iX = in+n2-7;
 404   DATA_TYPE *oX = out+n2+n4;
 405   DATA_TYPE *T  = init->trig+n4;
 406
 407   do{
 408     oX         -= 4;
 409     oX[0]       = MULT_NORM(-iX[2] * T[3] - iX[0]  * T[2]);
 410     oX[1]       = MULT_NORM (iX[0] * T[3] - iX[2]  * T[2]);
 411     oX[2]       = MULT_NORM(-iX[6] * T[1] - iX[4]  * T[0]);
 412     oX[3]       = MULT_NORM (iX[4] * T[1] - iX[6]  * T[0]);
 413     iX         -= 8;
 414     T          += 4;
 415   }while(iX>=in);
 416
 417   iX            = in+n2-8;
 418   oX            = out+n2+n4;
 419   T             = init->trig+n4;
 420
 421   do{
 422     T          -= 4;
 423     oX[0]       =  MULT_NORM (iX[4] * T[3] + iX[6] * T[2]);
 424     oX[1]       =  MULT_NORM (iX[4] * T[2] - iX[6] * T[3]);
 425     oX[2]       =  MULT_NORM (iX[0] * T[1] + iX[2] * T[0]);
 426     oX[3]       =  MULT_NORM (iX[0] * T[0] - iX[2] * T[1]);
 427     iX         -= 8;
 428     oX         += 4;
 429   }while(iX>=in);
 430
 431   mdct_butterflies(init,out+n2,n2);
 432   mdct_bitreverse(init,out);
 433
 434   /* roatate + window */
 435
 436   {
 437     DATA_TYPE *oX1=out+n2+n4;
 438     DATA_TYPE *oX2=out+n2+n4;
 439     DATA_TYPE *iX =out;
 440     T             =init->trig+n2;
 441
 442     do{
 443       oX1-=4;
 444
 445       oX1[3]  =  MULT_NORM (iX[0] * T[1] - iX[1] * T[0]);
 446       oX2[0]  = -MULT_NORM (iX[0] * T[0] + iX[1] * T[1]);
 447
 448       oX1[2]  =  MULT_NORM (iX[2] * T[3] - iX[3] * T[2]);
 449       oX2[1]  = -MULT_NORM (iX[2] * T[2] + iX[3] * T[3]);
 450
 451       oX1[1]  =  MULT_NORM (iX[4] * T[5] - iX[5] * T[4]);
 452       oX2[2]  = -MULT_NORM (iX[4] * T[4] + iX[5] * T[5]);
 453
 454       oX1[0]  =  MULT_NORM (iX[6] * T[7] - iX[7] * T[6]);
 455       oX2[3]  = -MULT_NORM (iX[6] * T[6] + iX[7] * T[7]);
 456
 457       oX2+=4;
 458       iX    +=   8;
 459       T     +=   8;
 460     }while(iX<oX1);
 461
 462     iX=out+n2+n4;
 463     oX1=out+n4;
 464     oX2=oX1;
 465
 466     do{
 467       oX1-=4;
 468       iX-=4;
 469
 470       oX2[0] = -(oX1[3] = iX[3]);
 471       oX2[1] = -(oX1[2] = iX[2]);
 472       oX2[2] = -(oX1[1] = iX[1]);
 473       oX2[3] = -(oX1[0] = iX[0]);
 474
 475       oX2+=4;
 476     }while(oX2<iX);
 477
 478     iX=out+n2+n4;
 479     oX1=out+n2+n4;
 480     oX2=out+n2;
 481     do{
 482       oX1-=4;
 483       oX1[0]= iX[3];
 484       oX1[1]= iX[2];
 485       oX1[2]= iX[1];
 486       oX1[3]= iX[0];
 487       iX+=4;
 488     }while(oX1>oX2);
 489   }
 490 }
 491
 492 void mdct_forward(mdct_lookup *init, DATA_TYPE *in, DATA_TYPE *out){
 493   int n=init->n;
 494   int n2=n>>1;
 495   int n4=n>>2;
 496   int n8=n>>3;
 497   DATA_TYPE *w=alloca(n*sizeof(*w)); /* forward needs working space */
 498   DATA_TYPE *w2=w+n2;
 499
 500   /* rotate */
 501
 502   /* window + rotate + step 1 */
 503
 504   REG_TYPE r0;
 505   REG_TYPE r1;
 506   DATA_TYPE *x0=in+n2+n4;
 507   DATA_TYPE *x1=x0+1;
 508   DATA_TYPE *T=init->trig+n2;
 509
 510   int i=0;
 511
 512   for(i=0;i<n8;i+=2){
 513     x0 -=4;
 514     T-=2;
 515     r0= x0[2] + x1[0];
 516     r1= x0[0] + x1[2];
 517     w2[i]=   MULT_NORM(r1*T[1] + r0*T[0]);
 518     w2[i+1]= MULT_NORM(r1*T[0] - r0*T[1]);
 519     x1 +=4;
 520   }
 521
 522   x1=in+1;
 523
 524   for(;i<n2-n8;i+=2){
 525     T-=2;
 526     x0 -=4;
 527     r0= x0[2] - x1[0];
 528     r1= x0[0] - x1[2];
 529     w2[i]=   MULT_NORM(r1*T[1] + r0*T[0]);
 530     w2[i+1]= MULT_NORM(r1*T[0] - r0*T[1]);
 531     x1 +=4;
 532   }
 533
 534   x0=in+n;
 535
 536   for(;i<n2;i+=2){
 537     T-=2;
 538     x0 -=4;
 539     r0= -x0[2] - x1[0];
 540     r1= -x0[0] - x1[2];
 541     w2[i]=   MULT_NORM(r1*T[1] + r0*T[0]);
 542     w2[i+1]= MULT_NORM(r1*T[0] - r0*T[1]);
 543     x1 +=4;
 544   }
 545
 546
 547   mdct_butterflies(init,w+n2,n2);
 548   mdct_bitreverse(init,w);
 549
 550   /* roatate + window */
 551
 552   T=init->trig+n2;
 553   x0=out+n2;
 554
 555   for(i=0;i<n4;i++){
 556     x0--;
 557     out[i] =MULT_NORM((w[0]*T[0]+w[1]*T[1])*init->scale);
 558     x0[0]  =MULT_NORM((w[0]*T[1]-w[1]*T[0])*init->scale);
 559     w+=2;
 560     T+=2;
 561   }
 562 }