001/* JOrbis
002 * Copyright (C) 2000 ymnk, JCraft,Inc.
003 *  
004 * Written by: 2000 ymnk<ymnk@jcaft.com>
005 *   
006 * Many thanks to 
007 *   Monty <monty@xiph.org> and 
008 *   The XIPHOPHORUS Company http://www.xiph.org/ .
009 * JOrbis has been based on their awesome works, Vorbis codec.
010 *   
011 * This program is free software; you can redistribute it and/or
012 * modify it under the terms of the GNU Library General Public License
013 * as published by the Free Software Foundation; either version 2 of
014 * the License, or (at your option) any later version.
015   
016 * This program is distributed in the hope that it will be useful,
017 * but WITHOUT ANY WARRANTY; without even the implied warranty of
018 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
019 * GNU Library General Public License for more details.
020 * 
021 * You should have received a copy of the GNU Library General Public
022 * License along with this program; if not, write to the Free Software
023 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
024 */
025
026package com.jcraft.jorbis;
027
028class Mdct{
029
030  static private final float cPI3_8=0.38268343236508977175f;
031  static private final float cPI2_8=0.70710678118654752441f;
032  static private final float cPI1_8=0.92387953251128675613f;
033
034  int n;
035  int log2n;
036  
037  float[] trig;
038  int[] bitrev;
039
040  float scale;
041
042  void init(int n){
043    bitrev=new int[n/4];
044    trig=new float[n+n/4];
045
046    int n2=n>>>1;
047    log2n=(int)Math.rint(Math.log(n)/Math.log(2));
048    this.n=n;
049
050
051    int AE=0;
052    int AO=1;
053    int BE=AE+n/2;
054    int BO=BE+1;
055    int CE=BE+n/2;
056    int CO=CE+1;
057    // trig lookups...
058    for(int i=0;i<n/4;i++){
059      trig[AE+i*2]=(float)Math.cos((Math.PI/n)*(4*i));
060      trig[AO+i*2]=(float)-Math.sin((Math.PI/n)*(4*i));
061      trig[BE+i*2]=(float)Math.cos((Math.PI/(2*n))*(2*i+1));
062      trig[BO+i*2]=(float)Math.sin((Math.PI/(2*n))*(2*i+1));
063    }
064    for(int i=0;i<n/8;i++){
065      trig[CE+i*2]=(float)Math.cos((Math.PI/n)*(4*i+2));
066      trig[CO+i*2]=(float)-Math.sin((Math.PI/n)*(4*i+2));
067    }
068
069    {
070      int mask=(1<<(log2n-1))-1;
071      int msb=1<<(log2n-2);
072      for(int i=0;i<n/8;i++){
073        int acc=0;
074        for(int j=0;msb>>>j!=0;j++)
075          if(((msb>>>j)&i)!=0)acc|=1<<j;
076        bitrev[i*2]=((~acc)&mask);
077//      bitrev[i*2]=((~acc)&mask)-1;
078        bitrev[i*2+1]=acc;
079      }
080    }
081    scale=4.f/n;
082  }
083
084  void clear(){
085  }
086
087  void forward(float[] in, float[] out){
088  }
089
090  float[] _x=new float[1024];
091  float[] _w=new float[1024];
092
093  synchronized void backward(float[] in, float[] out){
094    if(_x.length<n/2){_x=new float[n/2];}
095    if(_w.length<n/2){_w=new float[n/2];}
096    float[] x=_x;
097    float[] w=_w;
098    int n2=n>>>1;
099    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}