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 *bitA=bitrev;
  59   int *bitB=bitrev+1;
  60
  61   int i;
  62   int log2n=lookup->log2n=rint(log(n)/log(2));
  63   lookup->n=n;
  64   lookup->trig=trig;
  65   lookup->bitrev=bitrev;
  66
  67   /* trig lookups... */
  68
  69   for(i=0;i<n/4;i++){
  70     AE[i*2]=cos((M_PI/n)*(4*i));
  71     AO[i*2]=-sin((M_PI/n)*(4*i));
  72     BE[i*2]=cos((M_PI/(2*n))*(2*i+1));
  73     BO[i*2]=sin((M_PI/(2*n))*(2*i+1));
  74   }
  75   for(i=0;i<n/8;i++){
  76     CE[i*2]=cos((M_PI/n)*(4*i+2));
  77     CO[i*2]=-sin((M_PI/n)*(4*i+2));
  78   }
  79
  80   /* bitreverse lookup... */
  81
  82   {
  83     int mask=(1<<(log2n-1))-1,i,j;
  84     int msb=1<<(log2n-2);
  85     for(i=0;i<n/8;i++){
  86       int acc=0;
  87       for(j=0;msb>>j;j++)
  88         if((msb>>j)&i)acc|=1<<j;
  89       bitA[i*2]=((~acc)&mask);
  90       bitB[i*2+1]=acc;
  91     }
  92   }
  93 }
  94
  95 void mdct_clear(mdct_lookup *l){
  96   if(l){
  97     if(l->trig)free(l->trig);
  98     if(l->bitrev)free(l->bitrev);
  99     memset(l,0,sizeof(mdct_lookup));
 100   }
 101 }
 102
 103 static inline double *_mdct_kernel(double *x, double *w,
 104                                 int n, int n2, int n4, int n8,
 105                                 mdct_lookup *init){
 106   int i;
 107   /* step 2 */
 108
 109   {
 110     double *xA=x+n4;
 111     double *xB=x;
 112     double *w2=w+n4;
 113     double *A=init->trig+n2;
 114
 115     for(i=0;i<n4;){
 116       double x0=*xA - *xB;
 117       double x1;
 118       w2[i]=    *xA++ + *xB++;
 119
 120
 121       x1=       *xA - *xB;
 122       A-=4;
 123
 124       w[i++]=   x0 * A[0] + x1 * A[1];
 125       w[i]=     x1 * A[0] - x0 * A[1];
 126
 127       w2[i++]=  *xA++ + *xB++;
 128
 129     }
 130   }
 131
 132   /* step 3 */
 133
 134   {
 135     int r,s;
 136     for(i=0;i<init->log2n-3;i++){
 137       int k0=n>>(i+2);
 138       int k1=1<<(i+3);
 139       int wbase=n2-2;
 140       double *A=init->trig;
 141       double *temp;
 142
 143       for(r=0;r<(k0>>2);r++){
 144         int w1=wbase;
 145         int w2=w1-(k0>>1);
 146         double AEv= A[0],wA;
 147         double AOv= A[1],wB;
 148         wbase-=2;
 149
 150         k0++;
 151         for(s=0;s<(2<<i);s++){
 152           wB     =w[w1]   -w[w2];
 153           x[w1]  =w[w1]   +w[w2];
 154
 155           wA     =w[++w1] -w[++w2];
 156           x[w1]  =w[w1]   +w[w2];
 157
 158           x[w2]  =wA*AEv  - wB*AOv;
 159           x[w2-1]=wB*AEv  + wA*AOv;
 160
 161           w1-=k0;
 162           w2-=k0;
 163         }
 164         k0--;
 165
 166         A+=k1;
 167       }
 168
 169       temp=w;
 170       w=x;
 171       x=temp;
 172     }
 173   }
 174
 175   /* step 4, 5, 6, 7 */
 176   {
 177     double *C=init->trig+n;
 178     int *bit=init->bitrev;
 179     double *x1=x;
 180     double *x2=x+n2-1;
 181     for(i=0;i<n8;i++){
 182       int t1=*bit++;
 183       int t2=*bit++;
 184
 185       double wA=w[t1]-w[t2+1];
 186       double wB=w[t1-1]+w[t2];
 187       double wC=w[t1]+w[t2+1];
 188       double wD=w[t1-1]-w[t2];
 189
 190       double wACE=wA* *C;
 191       double wBCE=wB* *C++;
 192       double wACO=wA* *C;
 193       double wBCO=wB* *C++;
 194
 195       *x1++=( wC+wACO+wBCE)*.5;
 196       *x2--=(-wD+wBCO-wACE)*.5;
 197       *x1++=( wD+wBCO-wACE)*.5;
 198       *x2--=( wC-wACO-wBCE)*.5;
 199     }
 200   }
 201   return(x);
 202 }
 203
 204 void mdct_forward(mdct_lookup *init, double *in, double *out, double *window){
 205   int n=init->n;
 206   double x[n/2];
 207   double w[n/2];
 208   double *xx;
 209   int n2=n>>1;
 210   int n4=n>>2;
 211   int n8=n>>3;
 212   int i;
 213
 214   /* window + rotate + step 1 */
 215   {
 216     double tempA,tempB;
 217     int in1=n2+n4-4;
 218     int in2=in1+5;
 219     double *A=init->trig+n2;
 220
 221     i=0;
 222
 223     for(i=0;i<n8;i+=2){
 224       A-=2;
 225       tempA= in[in1+2]*window[in1+2] + in[in2]*window[in2];
 226       tempB= in[in1]*window[in1] + in[in2+2]*window[in2+2];
 227       in1 -=4;in2 +=4;
 228       x[i]=   tempB*A[1] + tempA*A[0];
 229       x[i+1]= tempB*A[0] - tempA*A[1];
 230     }
 231
 232     in2=1;
 233
 234     for(;i<n2-n8;i+=2){
 235       A-=2;
 236       tempA= in[in1+2]*window[in1+2] - in[in2]*window[in2];
 237       tempB= in[in1]*window[in1] - in[in2+2]*window[in2+2];
 238       in1 -=4;in2 +=4;
 239       x[i]=   tempB*A[1] + tempA*A[0];
 240       x[i+1]= tempB*A[0] - tempA*A[1];
 241     }
 242
 243     in1=n-4;
 244
 245     for(;i<n2;i+=2){
 246       A-=2;
 247       tempA= -in[in1+2]*window[in1+2] - in[in2]*window[in2];
 248       tempB= -in[in1]*window[in1] - in[in2+2]*window[in2+2];
 249       in1 -=4;in2 +=4;
 250       x[i]=   tempB*A[1] + tempA*A[0];
 251       x[i+1]= tempB*A[0] - tempA*A[1];
 252     }
 253   }
 254
 255   xx=_mdct_kernel(x,w,n,n2,n4,n8,init);
 256
 257   /* step 8 */
 258
 259   {
 260     double *B=init->trig+n2;
 261     double *out2=out+n2;
 262     double scale=4./n;
 263     for(i=0;i<n4;i++){
 264       out[i]   =(xx[0]*B[0]+xx[1]*B[1])*scale;
 265       *(--out2)=(xx[0]*B[1]-xx[1]*B[0])*scale;
 266
 267       xx+=2;
 268       B+=2;
 269     }
 270   }
 271 }
 272
 273 void mdct_backward(mdct_lookup *init, double *in, double *out, double *window){
 274   int n=init->n;
 275   double x[n/2];
 276   double w[n/2];
 277   double *xx;
 278   int n2=n>>1;
 279   int n4=n>>2;
 280   int n8=n>>3;
 281   int i;
 282
 283   /* window + rotate + step 1 */
 284   {
 285     double *inO=in+1;
 286     double  *xO= x;
 287     double  *A=init->trig+n2;
 288
 289     for(i=0;i<n8;i++){
 290       A-=2;
 291       *xO++=-*(inO+2)*A[1] - *inO*A[0];
 292       *xO++= *inO*A[1] - *(inO+2)*A[0];
 293       inO+=4;
 294     }
 295
 296     inO=in+n2-4;
 297
 298     for(i=0;i<n8;i++){
 299       A-=2;
 300       *xO++=*inO*A[1] + *(inO+2)*A[0];
 301       *xO++=*inO*A[0] - *(inO+2)*A[1];
 302       inO-=4;
 303     }
 304
 305   }
 306
 307   xx=_mdct_kernel(x,w,n,n2,n4,n8,init);
 308
 309   /* step 8 */
 310
 311   {
 312     double *B=init->trig+n2;
 313     int o1=n4,o2=o1-1;
 314     int o3=n4+n2,o4=o3-1;
 315
 316     for(i=0;i<n4;i++){
 317       double temp1= (*xx * B[1] - *(xx+1) * B[0]);
 318       double temp2=-(*xx * B[0] + *(xx+1) * B[1]);
 319
 320       out[o1]=-temp1*window[o1];
 321       out[o2]= temp1*window[o2];
 322       out[o3]= temp2*window[o3];
 323       out[o4]= temp2*window[o4];
 324
 325       o1++;
 326       o2--;
 327       o3++;
 328       o4--;
 329       xx+=2;
 330       B+=2;
 331     }
 332   }
 333 }
 334
 335
 336
 337
 338