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