anx.framework/Support/OggUtils/csvorbis/Mdct.cs

/* csvorbis
 * Copyright (C) 2000 ymnk, JCraft,Inc.
 *  
 * Written by: 2000 ymnk<ymnk@jcraft.com>
 * Ported to C# from JOrbis by: Mark Crichton <crichton@gimp.org> 
 *   
 * Thanks go to the JOrbis team, for licencing the code under the
 * LGPL, making my job a lot easier.
 * 
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public License
 * as published by the Free Software Foundation; either version 2 of
 * the License, or (at your option) any later version.
   
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Library General Public License for more details.
 * 
 * You should have received a copy of the GNU Library General Public
 * License along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */


using System;
using System.Runtime.CompilerServices;
using csogg;

namespace csvorbis 
{
	class Mdct
	{

		//static private float cPI3_8=0.38268343236508977175f;
		//static private float cPI2_8=0.70710678118654752441f;
		//static private float cPI1_8=0.92387953251128675613f;

		int n;
		int log2n;
  
		float[] trig;
		int[] bitrev;

		float scale;

		internal void init(int n)
		{
			bitrev=new int[n/4];
			trig=new float[n+n/4];

			int n2=(int)((uint)n >> 1);
			log2n=(int)Math.Round(Math.Log(n)/Math.Log(2));
			this.n=n;


			int AE=0;
			int AO=1;
			int BE=AE+n/2;
			int BO=BE+1;
			int CE=BE+n/2;
			int CO=CE+1;
			// trig lookups...
			for(int i=0;i<n/4;i++)
			{
				trig[AE+i*2]=(float)Math.Cos((Math.PI/n)*(4*i));
				trig[AO+i*2]=(float)-Math.Sin((Math.PI/n)*(4*i));
				trig[BE+i*2]=(float)Math.Cos((Math.PI/(2*n))*(2*i+1));
				trig[BO+i*2]=(float)Math.Sin((Math.PI/(2*n))*(2*i+1));
			}
			for(int i=0;i<n/8;i++)
			{
				trig[CE+i*2]=(float)Math.Cos((Math.PI/n)*(4*i+2));
				trig[CO+i*2]=(float)-Math.Sin((Math.PI/n)*(4*i+2));
			}

		{
			int mask=(1<<(log2n-1))-1;
			int msb=1<<(log2n-2);
			for(int i=0;i<n/8;i++)
			{
				int acc=0;
				for(int j=0; (((uint)msb) >> j) != 0; j++)
					if(((((uint)msb>>j))&i) != 0) 
						acc |= 1 << j;
				bitrev[i*2]=((~acc)&mask);
				//	bitrev[i*2]=((~acc)&mask)-1;
				bitrev[i*2+1]=acc;
			}
		}
			scale=4.0f/n;
		}

		internal void clear()
		{
		}

		internal void forward(float[] fin, float[] fout)
		{
		}

		float[] _x=new float[1024];
		float[] _w=new float[1024];

		
		[MethodImpl(MethodImplOptions.Synchronized)]
		internal void backward(float[] fin, float[] fout)
		{
				if(_x.Length < n/2){_x=new float[n/2];}
				if(_w.Length < n/2){_w=new float[n/2];}
				float[] x=_x;
				float[] w=_w;
				int n2=(int)((uint)n >> 1);
				int n4=(int)((uint)n >> 2);
				int n8=(int)((uint)n >> 3);

				// rotate + step 1
			{
				int inO=1;
				int xO=0;
				int A=n2;

				int i;
				for(i=0;i<n8;i++)
				{
					A-=2;
					x[xO++]=-fin[inO+2]*trig[A+1] - fin[inO]*trig[A];
					x[xO++]= fin[inO]*trig[A+1] - fin[inO+2]*trig[A];
					inO+=4;
				}

				inO=n2-4;

				for(i=0;i<n8;i++)
				{
					A-=2;
					x[xO++]=fin[inO]*trig[A+1] + fin[inO+2]*trig[A];
					x[xO++]=fin[inO]*trig[A] - fin[inO+2]*trig[A+1];
					inO-=4;
				}
			}

				float[] xxx=mdct_kernel(x,w,n,n2,n4,n8);
				int xx=0;

				// step 8

			{
				int B=n2;
				int o1=n4,o2=o1-1;
				int o3=n4+n2,o4=o3-1;
    
				for(int i=0;i<n4;i++)
				{
					float temp1= (xxx[xx] * trig[B+1] - xxx[xx+1] * trig[B]);
					float temp2=-(xxx[xx] * trig[B] + xxx[xx+1] * trig[B+1]);
    
					fout[o1]=-temp1;
					fout[o2]= temp1;
					fout[o3]= temp2;
					fout[o4]= temp2;

					o1++;
					o2--;
					o3++;
					o4--;
					xx+=2;
					B+=2;
				}
			}
		}
		internal float[] mdct_kernel(float[] x, float[] w,
			int n, int n2, int n4, int n8)
		{
			// step 2

			int xA=n4;
			int xB=0;
			int w2=n4;
			int A=n2;

			for(int i=0;i<n4;)
			{
				float x0=x[xA] - x[xB];
				float x1;
				w[w2+i]=x[xA++]+x[xB++];

				x1=x[xA]-x[xB];
				A-=4;

				w[i++]=   x0 * trig[A] + x1 * trig[A+1];
				w[i]=     x1 * trig[A] - x0 * trig[A+1];

				w[w2+i]=x[xA++]+x[xB++];
				i++;
			}

			// step 3

		{
			for(int i=0;i<log2n-3;i++)
			{
				int k0=(int)((uint)n >> (i+2));
				int k1=1 << (i+3);
				int wbase=n2-2;

				A=0;
				float[] temp;

				for(int r=0; r<((uint)k0>>2); r++)
				{
					int w1=wbase;
					w2=w1-(k0>>1);
					float AEv= trig[A],wA;
					float AOv= trig[A+1],wB;
					wbase-=2;
		      
					k0++;
					for(int s=0;s<(2<<i);s++)
					{
						wB     =w[w1]   -w[w2];
						x[w1]  =w[w1]   +w[w2];

						wA     =w[++w1] -w[++w2];
						x[w1]  =w[w1]   +w[w2];
	    
						x[w2]  =wA*AEv  - wB*AOv;
						x[w2-1]=wB*AEv  + wA*AOv;

						w1-=k0;
						w2-=k0;
					}
					k0--;
					A+=k1;
				}

				temp=w;
				w=x;
				x=temp;
			}
		}

			// step 4, 5, 6, 7
		{
			int C=n;
			int bit=0;
			int x1=0;
			int x2=n2-1;

			for(int i=0;i<n8;i++)
			{
				int t1=bitrev[bit++];
				int t2=bitrev[bit++];

				float wA=w[t1]-w[t2+1];
				float wB=w[t1-1]+w[t2];
				float wC=w[t1]+w[t2+1];
				float wD=w[t1-1]-w[t2];

				float wACE=wA* trig[C];
				float wBCE=wB* trig[C++];
				float wACO=wA* trig[C];
				float wBCO=wB* trig[C++];
      
				x[x1++]=( wC+wACO+wBCE)*.5f;
				x[x2--]=(-wD+wBCO-wACE)*.5f;
				x[x1++]=( wD+wBCO-wACE)*.5f; 
				x[x2--]=( wC-wACO-wBCE)*.5f;
			}
		}
			return(x);
		}
	}
}