lib/mdct.c

   1 /********************************************************************
   2  *                                                                  *
   3  * THIS FILE IS PART OF THE Ogg Vorbis SOFTWARE CODEC SOURCE CODE.  *
   4  * USE, DISTRIBUTION AND REPRODUCTION OF THIS SOURCE IS GOVERNED BY *
   5  * THE GNU PUBLIC LICENSE 2, WHICH IS INCLUDED WITH THIS SOURCE.    *
   6  * PLEASE READ THESE TERMS DISTRIBUTING.                            *
   7  *                                                                  *
   8  * THE OggSQUISH SOURCE CODE IS (C) COPYRIGHT 1994-1999             *
   9  * by 1999 Monty <monty@xiph.org> and The XIPHOPHORUS Company       *
  10  * http://www.xiph.org/                                             *
  11  *                                                                  *
  12  ********************************************************************
  13
  14  function: modified discrete cosine transform
  15            power of two length transform only [16 <= n ]
  16
  17  author: Monty <xiphmont@mit.edu>
  18  modifications by: Monty
  19  last modification date: Oct 10 1999
  20
  21  Algorithm adapted from _The use of multirate filter banks for coding
  22  of high quality digital audio_, by T. Sporer, K. Brandenburg and
  23  B. Edler, collection of the European Signal Processing Conference
  24  (EUSIPCO), Amsterdam, June 1992, Vol.1, pp 211-214
  25
  26  Note that the below code won't make much sense without the paper;
  27  The presented algorithm was already fairly polished, and the code
  28  once followed it closely.  The current code both corrects several
  29  typos in the paper and goes beyond the presented optimizations
  30  (steps 4 through 6 are, for example, entirely eliminated).
  31
  32  This module DOES NOT INCLUDE code to generate the window function.
  33  Everybody has their own weird favorite including me... I happen to
  34  like the properties of y=sin(2PI*sin^2(x)), but others may vehemently
  35  disagree.
  36
  37  ********************************************************************/
  38
  39 #include <stdio.h>
  40 #include <stdlib.h>
  41 #include <string.h>
  42 #include <math.h>
  43 #include "mdct.h"
  44
  45 /* build lookups for trig functions; also pre-figure scaling and
  46    some window function algebra. */
  47
  48 void mdct_init(mdct_lookup *lookup,int n){
  49   int    *bitrev=malloc(sizeof(int)*(n/4));
  50   double *trig=malloc(sizeof(double)*(n+n/4));
  51   double *AE=trig;
  52   double *AO=trig+1;
  53   double *BE=AE+n/2;
  54   double *BO=BE+1;
  55   double *CE=BE+n/2;
  56   double *CO=CE+1;
  57
  58   int i;
  59   int log2n=lookup->log2n=rint(log(n)/log(2));
  60   lookup->n=n;
  61   lookup->trig=trig;
  62   lookup->bitrev=bitrev;
  63
  64   /* trig lookups... */
  65
  66   for(i=0;i<n/4;i++){
  67     AE[i*2]=cos((M_PI/n)*(4*i));
  68     AO[i*2]=-sin((M_PI/n)*(4*i));
  69     BE[i*2]=cos((M_PI/(2*n))*(2*i+1));
  70     BO[i*2]=sin((M_PI/(2*n))*(2*i+1));
  71   }
  72   for(i=0;i<n/8;i++){
  73     CE[i*2]=cos((M_PI/n)*(4*i+2));
  74     CO[i*2]=-sin((M_PI/n)*(4*i+2));
  75   }
  76
  77   /* bitreverse lookup... */
  78
  79   {
  80     int mask=(1<<(log2n-1))-1,i,j;
  81     int msb=1<<(log2n-2);
  82     for(i=0;i<n/8;i++){
  83       int acc=0;
  84       for(j=0;msb>>j;j++)
  85         if((msb>>j)&i)acc|=1<<j;
  86       bitrev[i*2]=((~acc)&mask);
  87       bitrev[i*2+1]=acc;
  88     }
  89   }
  90 }
  91
  92 void mdct_clear(mdct_lookup *l){
  93   if(l){
  94     if(l->trig)free(l->trig);
  95     if(l->bitrev)free(l->bitrev);
  96     memset(l,0,sizeof(mdct_lookup));
  97   }
  98 }
  99
 100 static inline double *_mdct_kernel(double *x, double *w,
 101                                 int n, int n2, int n4, int n8,
 102                                 mdct_lookup *init){
 103   int i;
 104   /* step 2 */
 105
 106   {
 107     double *xA=x+n4;
 108     double *xB=x;
 109     double *w2=w+n4;
 110     double *A=init->trig+n2;
 111
 112     for(i=0;i<n4;){
 113       double x0=*xA - *xB;
 114       double x1;
 115       w2[i]=    *xA++ + *xB++;
 116
 117
 118       x1=       *xA - *xB;
 119       A-=4;
 120
 121       w[i++]=   x0 * A[0] + x1 * A[1];
 122       w[i]=     x1 * A[0] - x0 * A[1];
 123
 124       w2[i++]=  *xA++ + *xB++;
 125
 126     }
 127   }
 128
 129   /* step 3 */
 130
 131   {
 132     int r,s;
 133     for(i=0;i<init->log2n-3;i++){
 134       int k0=n>>(i+2);
 135       int k1=1<<(i+3);
 136       int wbase=n2-2;
 137       double *A=init->trig;
 138       double *temp;
 139
 140       for(r=0;r<(k0>>2);r++){
 141         int w1=wbase;
 142         int w2=w1-(k0>>1);
 143         double AEv= A[0],wA;
 144         double AOv= A[1],wB;
 145         wbase-=2;
 146
 147         k0++;
 148         for(s=0;s<(2<<i);s++){
 149           wB     =w[w1]   -w[w2];
 150           x[w1]  =w[w1]   +w[w2];
 151
 152           wA     =w[++w1] -w[++w2];
 153           x[w1]  =w[w1]   +w[w2];
 154
 155           x[w2]  =wA*AEv  - wB*AOv;
 156           x[w2-1]=wB*AEv  + wA*AOv;
 157
 158           w1-=k0;
 159           w2-=k0;
 160         }
 161         k0--;
 162
 163         A+=k1;
 164       }
 165
 166       temp=w;
 167       w=x;
 168       x=temp;
 169     }
 170   }
 171
 172   /* step 4, 5, 6, 7 */
 173   {
 174     double *C=init->trig+n;
 175     int *bit=init->bitrev;
 176     double *x1=x;
 177     double *x2=x+n2-1;
 178     for(i=0;i<n8;i++){
 179       int t1=*bit++;
 180       int t2=*bit++;
 181
 182       double wA=w[t1]-w[t2+1];
 183       double wB=w[t1-1]+w[t2];
 184       double wC=w[t1]+w[t2+1];
 185       double wD=w[t1-1]-w[t2];
 186
 187       double wACE=wA* *C;
 188       double wBCE=wB* *C++;
 189       double wACO=wA* *C;
 190       double wBCO=wB* *C++;
 191
 192       *x1++=( wC+wACO+wBCE)*.5;
 193       *x2--=(-wD+wBCO-wACE)*.5;
 194       *x1++=( wD+wBCO-wACE)*.5;
 195       *x2--=( wC-wACO-wBCE)*.5;
 196     }
 197   }
 198   return(x);
 199 }
 200
 201 void mdct_forward(mdct_lookup *init, double *in, double *out, double *window){
 202   int n=init->n;
 203   double x[n/2];
 204   double w[n/2];
 205   double *xx;
 206   int n2=n>>1;
 207   int n4=n>>2;
 208   int n8=n>>3;
 209   int i;
 210
 211   /* window + rotate + step 1 */
 212   {
 213     double tempA,tempB;
 214     int in1=n2+n4-4;
 215     int in2=in1+5;
 216     double *A=init->trig+n2;
 217
 218     i=0;
 219
 220     for(i=0;i<n8;i+=2){
 221       A-=2;
 222       tempA= in[in1+2]*window[in1+2] + in[in2]*window[in2];
 223       tempB= in[in1]*window[in1] + in[in2+2]*window[in2+2];
 224       in1 -=4;in2 +=4;
 225       x[i]=   tempB*A[1] + tempA*A[0];
 226       x[i+1]= tempB*A[0] - tempA*A[1];
 227     }
 228
 229     in2=1;
 230
 231     for(;i<n2-n8;i+=2){
 232       A-=2;
 233       tempA= in[in1+2]*window[in1+2] - in[in2]*window[in2];
 234       tempB= in[in1]*window[in1] - in[in2+2]*window[in2+2];
 235       in1 -=4;in2 +=4;
 236       x[i]=   tempB*A[1] + tempA*A[0];
 237       x[i+1]= tempB*A[0] - tempA*A[1];
 238     }
 239
 240     in1=n-4;
 241
 242     for(;i<n2;i+=2){
 243       A-=2;
 244       tempA= -in[in1+2]*window[in1+2] - in[in2]*window[in2];
 245       tempB= -in[in1]*window[in1] - in[in2+2]*window[in2+2];
 246       in1 -=4;in2 +=4;
 247       x[i]=   tempB*A[1] + tempA*A[0];
 248       x[i+1]= tempB*A[0] - tempA*A[1];
 249     }
 250   }
 251
 252   xx=_mdct_kernel(x,w,n,n2,n4,n8,init);
 253
 254   /* step 8 */
 255
 256   {
 257     double *B=init->trig+n2;
 258     double *out2=out+n2;
 259     double scale=4./n;
 260     for(i=0;i<n4;i++){
 261       out[i]   =(xx[0]*B[0]+xx[1]*B[1])*scale;
 262       *(--out2)=(xx[0]*B[1]-xx[1]*B[0])*scale;
 263
 264       xx+=2;
 265       B+=2;
 266     }
 267   }
 268 }
 269
 270 void mdct_backward(mdct_lookup *init, double *in, double *out, double *window){
 271   int n=init->n;
 272   double x[n/2];
 273   double w[n/2];
 274   double *xx;
 275   int n2=n>>1;
 276   int n4=n>>2;
 277   int n8=n>>3;
 278   int i;
 279
 280   /* window + rotate + step 1 */
 281   {
 282     double *inO=in+1;
 283     double  *xO= x;
 284     double  *A=init->trig+n2;
 285
 286     for(i=0;i<n8;i++){
 287       A-=2;
 288       *xO++=-*(inO+2)*A[1] - *inO*A[0];
 289       *xO++= *inO*A[1] - *(inO+2)*A[0];
 290       inO+=4;
 291     }
 292
 293     inO=in+n2-4;
 294
 295     for(i=0;i<n8;i++){
 296       A-=2;
 297       *xO++=*inO*A[1] + *(inO+2)*A[0];
 298       *xO++=*inO*A[0] - *(inO+2)*A[1];
 299       inO-=4;
 300     }
 301
 302   }
 303
 304   xx=_mdct_kernel(x,w,n,n2,n4,n8,init);
 305
 306   /* step 8 */
 307
 308   {
 309     double *B=init->trig+n2;
 310     int o1=n4,o2=o1-1;
 311     int o3=n4+n2,o4=o3-1;
 312
 313     for(i=0;i<n4;i++){
 314       double temp1= (*xx * B[1] - *(xx+1) * B[0]);
 315       double temp2=-(*xx * B[0] + *(xx+1) * B[1]);
 316
 317       out[o1]=-temp1*window[o1];
 318       out[o2]= temp1*window[o2];
 319       out[o3]= temp2*window[o3];
 320       out[o4]= temp2*window[o4];
 321
 322       o1++;
 323       o2--;
 324       o3++;
 325       o4--;
 326       xx+=2;
 327       B+=2;
 328     }
 329   }
 330 }
 331
 332
 333
 334
 335