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