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