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}