package HTools.Math {
	
	import HTools.Math.Complex;
	import flash.utils.*;
	import flash.errors.*;
	
	public class FFT {
		
		
		
		public var values:Vector.<Complex>;
		
		private var exponentials:Vector.<Vector.<Complex> >;
		private var posMap:Vector.<uint>;
		private var invPosMap:Vector.<uint>;
		private var points:uint;
		private var logPoints:uint;
		
		private function buildExponentials(points:uint) {
			
			var i:uint;
			
			this.points=points;
			
			exponentials=new Vector.<Vector.<Complex> >();
			
			var sqrtPoints:Number = Math.sqrt(points);
			// calculate binary log
			
			var tempPoints=points-1;
			while (tempPoints != 0)	{
				tempPoints >>= 1;
				logPoints++;
			}
			
			var _2_l:int = 2;
			for (var l:uint = 1; l <= logPoints; l++)	{
				var subV:Vector.<Complex>=new Vector.<Complex>();
	
				for (i = 0; i < points; i++ )	{
						var re:Number =  Math.cos (2. * Math.PI * i / _2_l);
						var im:Number = -Math.sin (2. * Math.PI * i / _2_l);
						subV.push( new Complex(re, im) );
				}
				exponentials.push(subV);
				_2_l *= 2;
			}
			
			
			// set up bit reverse mapping
			posMap=new Vector.<uint>();
			invPosMap=new Vector.<uint>()
			invPosMap.length=points;
			var rev:uint = 0;
			var halfPoints:uint = points/2;
			for (i = 0; i < points - 1; i++){
				posMap[i] = rev;
				invPosMap[rev]=i;
				var mask:uint = halfPoints;
				// add 1 backwards

				while (rev >= mask) {
					rev -= mask; // turn off this bit
					mask >>= 1;
					if(rev==0 && mask == 0) {
						throw ArgumentError("Points has to be a multiple of 2!");
					}
				}
				rev += mask;
			}
			posMap[points-1] = points-1;
		}
		
		public function loadSamples(b:ByteArray) {
			
			values=new Vector.<Complex>();
			values.length=points;
			
			for(var i:uint;i<points;++i) {
				values[posMap[points-i-1]]=new Complex(b.readFloat(),0);
				CONFIG::DOUBLE_READ {
					b.readFloat();
				}
			}
			
		}
		
		public function unloadSamples(b:ByteArray) {
			
			for(var i:uint=points;i>0;i--) {
				var sample=values[posMap[i-1]].real;
				b.writeFloat(sample);
				CONFIG::DOUBLE_READ {
					b.writeFloat(sample)
				}
			}
			
		}
		
		public function FFT(points:uint) {
			buildExponentials(points);
		}
		
		public function newLoadSamples(b:ByteArray):void {
			
			values=new Vector.<Complex>();
			values.length=points;
			
			var nd2:uint=points/2;
			
			var j:uint=nd2;
			var i:uint;
			var t:Complex;
			
			for(i=0;i<points;i++) {
				values[i]=new Complex(b.readFloat());
				CONFIG::DOUBLE_READ {
					b.readFloat();
				}
			}
			
			//bitreversal sorting
			for(i=1;i<=points-2;++i) {
				if(i<j){
					//swap points
					t=values[j];
					values[j]=values[i];
					values[i]=t;
				}
				var k:uint=nd2;
				while(k<=j) {
					j-=k;
					k/=2;
				}
				j+=k;
			}
		}
		
		public function newTransform():void {
			var m:uint=Math.floor(Math.log(points)/Math.log(2));
			var nm1:uint=points-1;
			
			var t:Complex;
			var j:uint,i:uint;
			
			for(var l:uint=1;l<=m;++l) {//each stage
				var le:uint=Math.floor(Math.pow(2,l));
				var le2:uint=le/2;
				
				var u:Complex=new Complex(1,0);
				var s:Complex=new Complex(Math.cos(Math.PI/le2),-Math.sin(Math.PI/le2));
				
				for(j=1;j<=le2;++j) {//each sub dft
					var jm1:uint=j-1;
					
					for(i=jm1;i<=nm1;i+=le) {
						var ip:uint=i+le2;
						
						t=values[i].mulClone(values[ip]);
						values[ip]=values[i].subClone(t);
						values[i].addEq(t);
						
					}//end of i
					
					u.mulEq(s);
					
				}//end of j
				
			}//end of l
			
		}//end of method
		
		public function transform() {
			
			/*
				Fast fouirer transform class based on code copyright to Reliable Software 2006 (unknown license)
			*/
			
			var step:uint = 1;
			for (var level:uint = 1; level <= logPoints; level++)	{
				var increm:uint = step * 2;
				
				for (var j:uint = 0; j < step; j++) {
					// U = exp ( - 2 PI j / 2 ^ level )
					var U:Complex = exponentials[level-1][j];
					
					for (var i:uint = j; i < points; i += increm)	{
						// in-place butterfly
						// Xnew[i]      = X[i] + U * X[i+step]
						// Xnew[i+step] = X[i] - U * X[i+step]
						// X'=X+ U * Y
						// Y'=X- U * Y
						var T:Complex = U.clone();
						T.mulEq(values[i+step]);
						values[i+step]=values[i].clone();
						values[i+step].subEq(T);
						values[i].addEq(T);
					}
				}
				step *= 2;
			}
		}
		
		public function untransform():void {
			var i:uint;
			for(i=0;i<points;++i) {
				values[i].im=-values[i].im;
			}
			newTransform();
			for(i=0;i<points;++i) {
				values[i].real=values[i].real/points;
				values[i].im=values[i].im/-points;
			}
		}
		
		public function oldUntransform() {
			
			var step:uint = 1;
			for (var level:uint = 1; level <= logPoints; level++)	{
				var increm:uint = step * 2;
				
				for (var j:uint = 0; j < step; j++) {
					// U = exp ( - 2 PI j / 2 ^ level )
					var U:Complex = exponentials[level-1][j];
					
					for (var i:uint = j; i < points; i += increm)	{
						// in-place butterfly
						// Xnew[i]      = X[i] + U * X[i+step]
						// Xnew[i+step] = X[i] - U * X[i+step]
						// X'=X+ U * Y
						// Y'=X- U * Y
						
						var newX:Complex=values[i].clone();
						newX.addEq(values[i+step]);
						newX.intDivEq(2);
						//I can now clober the real copies
						values[i+step].intDivEq(2);
						values[i].subEq(values[i+step]);
						values[i].divEq(U);
						var newY:Complex=values[i];
						
						//now, I can overwrite the old data with these fresh objects.
						values[i]=newX;
						values[i+step]=newY;
					}
				}
				step *= 2;
			}
		}
		
		public function pointToFreq(point:uint):Number {
			return 44100/points*point;
		}
		
		public function freqToPoint(freq:Number) {
			return points * freq / 44100;
		}
	}
}