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;
}
}
}