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.