lib/mdct.c

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