Im Folgenden möchte ich kurz eine mögliche Lösung kommentiert vorstellen, durch die von gegebenen Eingangsdaten die diskrete Fourier-Transformierte berechnet werden kann. 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.
Zunächst müssen die benötigten header eingefügt werden. Hier sind das diejenigen für die Ein- und Ausgabe, mathematische Funktionen und komplexe Zahlen. Diese Implementierung ist zwar geringfügig langsamer, als sie es bei einer eigenen Lösung mit je zwei reellen Zahlen statt nur einer komplexen sein könnte, aber so ist das Ergebnis besser lesbar.
[cpp]#include <stdio.h>
#include <math.h>
#include <complex.h>[/cpp]
Die DFT setzt voraus, dass die Eingangsdaten einerseits gleichverteilte Stützstellen aufweisen und andererseits die Anzahl der Stützstellen eine ganzzahlige Potenz von 2 ist. Um diesem Umstand Rechnung zu tragen werden ein paar Konstanten eingeführt. Dabei ist DCOUNT die Anzahl der Datensätze, BUFLEN wird beim Einlesen der Daten verwendet und E enthält einen Näherungswert für die eulersche Zahl.
[cpp]#define DCOUNT 32768
#define BUFLEN 80
#define E 2.71828183[/cpp]
Dass DCOUNT als Konstante eingeführt wird, ist lediglich der Bequemlichkeit geschuldet. Da C keine nativen dynamischen Arrays kennt, wäre der Gebrauch von malloc nötig, was zwar nicht schwierig ist, aber die Übersicht im Code weiter reduziert und den Fokus von der DFT ablenkt. Aus diesem Grund ist auch der Dateiname mit den Daten direkt im Programmcode eingetragen. Das Einlesen einer Datei, in der in jeder Zeile je ein x- und ein y-Wert eingetragen ist, erfolgt dann so:
[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]
Nach diesem Abschnitt liegen die Stützstellen im Vektor x und die Funktionswerte an den Stützstellen im Vektor y vor. Der Vektor x ließe sich auch berechnen, da die Intervallbreite 2*PI beträgt und die Stützstellen ja gleichverteilt sein müssen. Das Herzstück, die eigentliche DFT besteht aus zwei geschachtelten Summen, wobei der Vektor w normalerweise mit omega bezeichnet wird.
[cpp] complex double w[DCOUNT];
complex double b[DCOUNT];
complex double a;
int k;
int m;
for (m = 0; m < DCOUNT; m++) {
w[m] = cpow(E, I*x[m]);
}
for (m = 0; m < DCOUNT; m++) {
b[m] = 0;
for (k = 0; k < DCOUNT; k++) {
a= y[k]*cpow(w[k], -m);
b[m] += a;
}
b[m] /= DCOUNT;
}[/cpp]
Daraus lassen sich dann noch die Koeffizienten berechnen, die hier auf der Konsole ausgegeben werden.
[cpp] /* Koeffizienten beta berechnen */
for (i = 0; i < DCOUNT; i++) {
a = b[i];
printf("-b %d %f %f %f\n", i, creal(b[i]), cimag(b[i]), cabs(a));
}
/* Koeffizienten B berechnen */
double Br[DCOUNT/2];
double Bi[DCOUNT/2];
for (i = 0; i < DCOUNT/2; i++) {
Bi[i] = creal(b[i])-creal(b[DCOUNT-i]);
Br[i] = cimag(b[DCOUNT-i])-cimag(b[i]);
printf("-B %d %f %f %f\n", i, Br[i], Bi[i], sqrt(Br[i]*Br[i]+Bi[i]*Bi[i]));
}
/* Koeffizienten A berechnen */
double Ar[DCOUNT/2];
double Ai[DCOUNT/2];
for (i = 0; i < DCOUNT/2; i++) {
Ai[i] = cimag(b[i])+cimag(b[DCOUNT-i]);
Ar[i] = creal(b[DCOUNT-i])+creal(b[i]);
printf("-A %d %f %f %f\n", i, Ar[i], Ai[i], sqrt(Ar[i]*Ar[i]+Ai[i]*Ai[i]));
}
return 0;
}[/cpp]