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