This is fft.m in view mode; [Download] [Up]
#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;
}
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.