Beispielimplementierung FFT Dieser Eintrag befasst sich mit einem Beispielprogramm für die Durchführung der Fast-Fourier-Transformation. Als Programmiersprache kommt C zum Einsatz. Wer nur am Programm interessiert ist, kann die einzelnen Code-Fragemente einfach nacheinander in eine .c-Datei einfügen. GCC benötigt noch das Flag -lm, durch das die mathematischen Funktionen gelinkt werden können.
Die benötigten header sind für die Ein- und Ausgabe, die mathematischen Funktionen und die hier verwendeten komplexen Zahlen nötig.
[cpp]#include <stdio.h>
#include <math.h>
#include <complex.h>[/cpp]
Die Konstanten berücksichtigen die Eigenheiten der FFT. So muss die Anzahl der Stützstellen – hier DCOUNT – eine ganzzahlige Potenz von zwei sein. Da an einigen Stellen des Programmcodes die Basis dieser Zahl noch relevant ist, gibt es eine eigene Konstante dafür: N. Es müssen somit immer beide Zahlen so angepasst werden, dass 2^N=DCOUNT gilt. E und PI sollten selbsterklärend sein, BUFLEN wird beim Einlesen der Stützdaten genutzt. Die Konstanten werden für die Array-Größen verwendet, um auf malloc verzichten zu können.
[cpp]#define DCOUNT 32768
#define BUFLEN 80
#define PI 3.14159265
#define E 2.71828183
#define N 15[/cpp]
Der weiter unten aufgeführte Programmcode berechnet zwar die FFT eines Eingabevektors, aber die Reihenfolge der Elemente im Ergebnisvektor ist nicht so, wie man es direkt erwarten würde. Vielmehr werden die Elemente nach dem Schmetterlingsgraph verknüpft, sodass die Bitreihenfolge des Index noch umgekehrt werden muss. Angenommen es gibt einen Vektor mit vier Elementen – da vier eine Potenz von zwei ist, wäre das eine valide Eingabe. Nun lassen sich alle Indizes binär darstellen: 00, 01, 10, 11. Dabei ist es wichtig, dass genau so viele führende Nullen angegeben werden, dass alle für den jeweiligen Vektor relevanten Indizes in der Binärdarstellung gleich lang sind. Kehrt man nun die Bitreihenfolge der Zahlen um, so ergibt sich der korrekte Index im Lösungsvektor des Programmteils, der die FFT durchführt: 00, 10, 01, 11. Die folgende Funktion dreht die Bitreihenfolge eines Integers um, setzt dabei aber voraus, dass dieser kleiner als DCOUNT ist.
[cpp]int shift(int val) {
int pos;
int ret = 0;
for (pos = 0; pos < N; pos++) {
ret *= 2; ret += (val >> pos) & 1;
}
return ret;
}[/cpp]
Dabei wird ausgenutzt, dass a >> b die Binärdarstellung der Zahl a um b Stellen nach rechts verschiebt – also eine Integerdivision durch 2^b. Damit ist das letzte Bit der Zahl gleich das b-te Bit in a. Die logische AND-Verknüpfung mittels & ergibt somit 1, wenn das b-te Bit der Ursprungszahl 1 war und andernfalls 0.
Auch wenn die Stützstellen durch DCOUNT bereits berechenbar wären, werden sie hier mit eingelesen. Die Datei “daten.dat” müsste dazu zwei Spalten mit den x- und y-Werten der Stützpunkte aufführen und DCOUNT Zeilen haben.
[cpp]int main(void) {
/* Daten einlesen */
char buffer[BUFLEN];
double x[DCOUNT];
double y[DCOUNT];
int i = 0;
FILE * fp = fopen("daten.dat", "r");
if (fp == NULL) {
printf("Oeffnen der Datei fehlgeschlagen.\n");
return 1;
}
while (fgets(buffer, BUFLEN, fp) != NULL) {
sscanf(buffer, "%lf%lf", &x[i], &y[i]); i++;
}[/cpp]
Bei der FFT werden immer zwei “alte” Elemente aus dem vorherigen Schritt (hier in f1 gespeichert) zu zwei “neuen” Elementen in f2 kombiniert. Die Initialwerte von f1 entsprechen den y-Werten der Stützpunkte.
[cpp] double complex f1[DCOUNT];
double complex f2[DCOUNT];
double complex e;
for (i = 0; i < DCOUNT; i++) {
f1[i] = y[i];
}[/cpp]
Der eigentliche Teil der FFT-Berechnung erfolgt in zwei geschachtelten Schleifen.
[cpp] /* FFT */
int r, k;
int M, R;
int li, ri;
for (i = N; i > 0; i–) {
M = pow(2, i-1);
R = pow(2, N-i);
e = cpow(E, -2*PI*I/pow(2, i));
for (r=0; r < R; r++) {
for (k = 0; k < M; k++) {
li = 2*r*M+k;
ri = (2*r+1)*M+k;
f2[li] = f1[li]+f1[ri];
f2[ri] = (f1[li]-f1[ri])*cpow(e, k);
}
}
for (r = 0; r < DCOUNT; r++) {
f1[r] = f2[r];
}
}[/cpp]
Nun müssen die Indizes wie o.a. geändert und normiert werden.
[cpp] /* Koeffizienten beta sortieren und normieren*/
printf("Koeffizienten beta (die ersten 160)\n");
double a,b;
for (r = 0; r < DCOUNT; r++) {
i = shift(r);
f2[r] = f1[i] / DCOUNT;
a = creal(f2[r]);
b = cimag(f2[r]);
if (r < 160) {
printf("-b %d %f %f %f\n", r, a, b, sqrt(a*a+b*b));
}
}[/cpp]
Abschließend können die Koeffizienten A und B errechnet werden.
[cpp] /* Koeffizienten A und B */
printf("Koeffizienten A und B (die ersten 160)\n");
double complex A, B;
for (r = 0; r < 160; r++) {
A = f2[r]+f2[DCOUNT-r];
B = I*(f2[r]-f2[DCOUNT-r]);
if (r == 0) {
A = 2*f2[0];
B = 0;
}
a = creal(A);
b = cimag(A);
printf("-A %d %f %f %f\n", r, a, b, sqrt(a*a+b*b));
a = creal(B);
b = cimag(B);
printf("-B %d %f %f %f\n", r, a, b, sqrt(a*a+b*b));
}
return 0;
}[/cpp]