#import "fft.h"

void rfa42(komplex *Daten, int n);

@implementation FFT

/* instance methods */

- initPoints:(int)points
{
	self = [super init];
	punkte = points;
	daten = (float *)malloc(sizeof(float) * points);
	return self;
}

- (float *)performFFTWithStereo:(short *)inData;
{
	int i;
	float tr,ti;
	komplex *ergebnis = (komplex *)daten;
	
	for(i=0;i<punkte;i++)
		daten[i] = inData[i*4] + inData[i*4 + 1];				// test, war * 4

	rfa42(ergebnis, punkte);

	for(i=0;i<punkte/2;i++)
	{
		tr = ergebnis[i].real;
		ti = ergebnis[i].imag;
	  daten[i] = sqrt(tr*tr + ti*ti) / punkte / 180 ; // 4096;	//TEST
	}

	return daten;
}

- (float *)performFFTWith:(short *)inData;
{
	int i;
	float tr,ti;
	komplex *ergebnis = (komplex *)daten;
	
	for(i=0;i<punkte;i++)
		daten[i] = inData[i];

	rfa42(ergebnis, punkte);

	for(i=0;i<punkte/2;i++)
	{
		tr = ergebnis[i].real;
		ti = ergebnis[i].imag;
	  daten[i] = sqrt(tr*tr + ti*ti) / punkte / 128; // 4096;	//TEST
	}

	return daten;
}

@end

void rfa42(komplex *Daten, int n)
{
	float sr,siwr,siwi,di,drwr,drwi,theta;
	int i,ldn,nhalbe,j,j2,j3,m,step,weite;
	komplex t,w,w2,w3,d0,d1,d2,d3;
	komplex *cd = Daten;
	
	for(i=1, ldn=0; 2*i <= n; ldn++)
	  i=i*2;
  n=i;
  if(n<2)
		return;
  nhalbe=n/2;
  for(j=0, i=0; i<nhalbe; i++)
  {
    if(i<j)
    { /* Austausch */
      t=cd[j];
      cd[j]=cd[i];
      cd[i]=t;
  	}
    m=nhalbe >> 1;
    while ((j & m) != 0)
		{
      j=j - m;
      m=m >> 1;
    }
		j+=m;
  }

  weite=1;
  if ((ldn % 2) == 0)
	{
    for(i=0; i<(nhalbe / 2); i++)
		{
      j=i << 1;
      t=cd[j];
      cd[j].real+=cd[j+1].real;
			cd[j].imag+=cd[j+1].imag;
      cd[j+1].real=t.real-cd[j+1].real;
			cd[j+1].imag=t.imag-cd[j+1].imag;
    }
    weite=2;
  }
  while (weite < nhalbe)
  {
    step=4*weite;
    for(m=0; m<weite; m++)
    {
      /* Exponentialfaktoren berechnen */
      theta=-M_PI*m/weite/2;
      w.real=cos(theta);
			w.imag=sin(theta);
      w2.real=w.real*w.real-w.imag*w.imag;
			w2.imag=w.real*w.imag*2;
      w3.real=w2.real*w.real-w2.imag*w.imag;
			w3.imag=w2.real*w.imag+w2.imag*w.real;

      i=m;
      do
			{
        j=i+weite;
				j2=j+weite;
				j3=j2+weite;
				
				d0=cd[i];
				d1.real=cd[j].real*w2.real-cd[j].imag*w2.imag;
				d1.imag=cd[j].real*w2.imag+cd[j].imag*w2.real;

				d2.real=cd[j2].real*w.real-cd[j2].imag*w.imag;
				d2.imag=cd[j2].real*w.imag+cd[j2].imag*w.real;

				d3.real=cd[j3].real*w3.real-cd[j3].imag*w3.imag;
				d3.imag=cd[j3].real*w3.imag+cd[j3].imag*w3.real;

				cd[i].real=d0.real+d1.real+d2.real+d3.real;
				cd[i].imag=d0.imag+d1.imag+d2.imag+d3.imag;

				cd[j].real=d0.real-d1.real+d2.imag-d3.imag;
				cd[j].imag=d0.imag-d1.imag-d2.real+d3.real;
				
				cd[j2].real=d0.real+d1.real-d2.real-d3.real;
				cd[j2].imag=d0.imag+d1.imag-d2.imag-d3.imag;

				cd[j3].real=d0.real-d1.real-d2.imag+d3.imag;
				cd[j3].imag=d0.imag-d1.imag+d2.real-d3.real;

				i+=step;
			}
      while(i<nhalbe);
    }
    weite=step;
  }

  t=cd[0];
	cd[0].real=t.real+t.imag;
	cd[0].imag=t.real-t.imag;

  for (i=1; i<(nhalbe / 2); i++)
  {
    theta=M_PI*i/nhalbe;
		w.real=cos(theta);
		w.imag=sin(theta);

    sr=(cd[i].real+cd[nhalbe-i].real)/2;
		drwr=(sr-cd[nhalbe-i].real)*w.real;
		drwi=(sr-cd[nhalbe-i].real)*w.imag;

    di=(cd[i].imag-cd[nhalbe-i].imag)/2;
		siwr=(di+cd[nhalbe-i].imag)*w.real;
		siwi=(di+cd[nhalbe-i].imag)*w.imag;
    
		cd[nhalbe-i].real=sr-siwr+drwi;
		cd[nhalbe-i].imag=-di-siwi-drwr;
		
    cd[i].real=sr+siwr-drwi;
		cd[i].imag=di-siwi-drwr;
	}
  cd[nhalbe / 2].imag=-cd[nhalbe / 2].imag;
}


