jorbis/src/com/jcraft/jorbis/Mdct.java

   1 /* JOrbis
   2  * Copyright (C) 2000 ymnk, JCraft,Inc.
   3  *
   4  * Written by: 2000 ymnk<ymnk@jcraft.com>
   5  *
   6  * Many thanks to
   7  *   Monty <monty@xiph.org> and
   8  *   The XIPHOPHORUS Company http://www.xiph.org/ .
   9  * JOrbis has been based on their awesome works, Vorbis codec.
  10  *
  11  * This program is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Library General Public License
  13  * as published by the Free Software Foundation; either version 2 of
  14  * the License, or (at your option) any later version.
  15
  16  * This program is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  19  * GNU Library General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Library General Public
  22  * License along with this program; if not, write to the Free Software
  23  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  24  */
  25
  26 package com.jcraft.jorbis;
  27
  28 class Mdct{
  29
  30   static private final float cPI3_8=0.38268343236508977175f;
  31   static private final float cPI2_8=0.70710678118654752441f;
  32   static private final float cPI1_8=0.92387953251128675613f;
  33
  34   int n;
  35   int log2n;
  36
  37   float[] trig;
  38   int[] bitrev;
  39
  40   float scale;
  41
  42   void init(int n){
  43     bitrev=new int[n/4];
  44     trig=new float[n+n/4];
  45
  46     int n2=n>>>1;
  47     log2n=(int)Math.rint(Math.log(n)/Math.log(2));
  48     this.n=n;
  49
  50
  51     int AE=0;
  52     int AO=1;
  53     int BE=AE+n/2;
  54     int BO=BE+1;
  55     int CE=BE+n/2;
  56     int CO=CE+1;
  57     // trig lookups...
  58     for(int i=0;i<n/4;i++){
  59       trig[AE+i*2]=(float)Math.cos((Math.PI/n)*(4*i));
  60       trig[AO+i*2]=(float)-Math.sin((Math.PI/n)*(4*i));
  61       trig[BE+i*2]=(float)Math.cos((Math.PI/(2*n))*(2*i+1));
  62       trig[BO+i*2]=(float)Math.sin((Math.PI/(2*n))*(2*i+1));
  63     }
  64     for(int i=0;i<n/8;i++){
  65       trig[CE+i*2]=(float)Math.cos((Math.PI/n)*(4*i+2));
  66       trig[CO+i*2]=(float)-Math.sin((Math.PI/n)*(4*i+2));
  67     }
  68
  69     {
  70       int mask=(1<<(log2n-1))-1;
  71       int msb=1<<(log2n-2);
  72       for(int i=0;i<n/8;i++){
  73         int acc=0;
  74         for(int j=0;msb>>>j!=0;j++)
  75           if(((msb>>>j)&i)!=0)acc|=1<<j;
  76         bitrev[i*2]=((~acc)&mask);
  77 //      bitrev[i*2]=((~acc)&mask)-1;
  78         bitrev[i*2+1]=acc;
  79       }
  80     }
  81     scale=4.f/n;
  82   }
  83
  84   void clear(){
  85   }
  86
  87   void forward(float[] in, float[] out){
  88   }
  89
  90   float[] _x=new float[1024];
  91   float[] _w=new float[1024];
  92
  93   synchronized void backward(float[] in, float[] out){
  94     if(_x.length<n/2){_x=new float[n/2];}
  95     if(_w.length<n/2){_w=new float[n/2];}
  96     float[] x=_x;
  97     float[] w=_w;
  98     int n2=n>>>1;
  99     int n4=n>>>2;
 100     int n8=n>>>3;
 101
 102     // rotate + step 1
 103     {
 104       int inO=1;
 105       int xO=0;
 106       int A=n2;
 107
 108       int i;
 109       for(i=0;i<n8;i++){
 110         A-=2;
 111         x[xO++]=-in[inO+2]*trig[A+1] - in[inO]*trig[A];
 112         x[xO++]= in[inO]*trig[A+1] - in[inO+2]*trig[A];
 113         inO+=4;
 114       }
 115
 116       inO=n2-4;
 117
 118       for(i=0;i<n8;i++){
 119         A-=2;
 120         x[xO++]=in[inO]*trig[A+1] + in[inO+2]*trig[A];
 121         x[xO++]=in[inO]*trig[A] - in[inO+2]*trig[A+1];
 122         inO-=4;
 123       }
 124     }
 125
 126     float[] xxx=mdct_kernel(x,w,n,n2,n4,n8);
 127     int xx=0;
 128
 129     // step 8
 130
 131     {
 132       int B=n2;
 133       int o1=n4,o2=o1-1;
 134       int o3=n4+n2,o4=o3-1;
 135
 136       for(int i=0;i<n4;i++){
 137         float temp1= (xxx[xx] * trig[B+1] - xxx[xx+1] * trig[B]);
 138         float temp2=-(xxx[xx] * trig[B] + xxx[xx+1] * trig[B+1]);
 139
 140         out[o1]=-temp1;
 141         out[o2]= temp1;
 142         out[o3]= temp2;
 143         out[o4]= temp2;
 144
 145         o1++;
 146         o2--;
 147         o3++;
 148         o4--;
 149         xx+=2;
 150         B+=2;
 151       }
 152     }
 153   }
 154   private float[] mdct_kernel(float[] x, float[] w,
 155                                int n, int n2, int n4, int n8){
 156     // step 2
 157
 158     int xA=n4;
 159     int xB=0;
 160     int w2=n4;
 161     int A=n2;
 162
 163     for(int i=0;i<n4;){
 164       float x0=x[xA] - x[xB];
 165       float x1;
 166       w[w2+i]=x[xA++]+x[xB++];
 167
 168       x1=x[xA]-x[xB];
 169       A-=4;
 170
 171       w[i++]=   x0 * trig[A] + x1 * trig[A+1];
 172       w[i]=     x1 * trig[A] - x0 * trig[A+1];
 173
 174       w[w2+i]=x[xA++]+x[xB++];
 175       i++;
 176     }
 177
 178     // step 3
 179
 180     {
 181       for(int i=0;i<log2n-3;i++){
 182         int k0=n>>>(i+2);
 183         int k1=1<<(i+3);
 184         int wbase=n2-2;
 185
 186         A=0;
 187         float[] temp;
 188
 189         for(int r=0;r<(k0>>>2);r++){
 190           int w1=wbase;
 191           w2=w1-(k0>>1);
 192           float AEv= trig[A],wA;
 193           float AOv= trig[A+1],wB;
 194           wbase-=2;
 195
 196           k0++;
 197           for(int s=0;s<(2<<i);s++){
 198             wB     =w[w1]   -w[w2];
 199             x[w1]  =w[w1]   +w[w2];
 200
 201             wA     =w[++w1] -w[++w2];
 202             x[w1]  =w[w1]   +w[w2];
 203
 204             x[w2]  =wA*AEv  - wB*AOv;
 205             x[w2-1]=wB*AEv  + wA*AOv;
 206
 207             w1-=k0;
 208             w2-=k0;
 209           }
 210           k0--;
 211           A+=k1;
 212         }
 213
 214         temp=w;
 215         w=x;
 216         x=temp;
 217       }
 218     }
 219
 220     // step 4, 5, 6, 7
 221     {
 222       int C=n;
 223       int bit=0;
 224       int x1=0;
 225       int x2=n2-1;
 226
 227       for(int i=0;i<n8;i++){
 228         int t1=bitrev[bit++];
 229         int t2=bitrev[bit++];
 230
 231         float wA=w[t1]-w[t2+1];
 232         float wB=w[t1-1]+w[t2];
 233         float wC=w[t1]+w[t2+1];
 234         float wD=w[t1-1]-w[t2];
 235
 236         float wACE=wA* trig[C];
 237         float wBCE=wB* trig[C++];
 238         float wACO=wA* trig[C];
 239         float wBCO=wB* trig[C++];
 240
 241         x[x1++]=( wC+wACO+wBCE)*.5f;
 242         x[x2--]=(-wD+wBCO-wACE)*.5f;
 243         x[x1++]=( wD+wBCO-wACE)*.5f;
 244         x[x2--]=( wC-wACO-wBCE)*.5f;
 245       }
 246     }
 247     return(x);
 248   }
 249 }