Fast Fourier Transform
Was haben Computertomographie, Kristallographie, Frequenzspektralanalyse, sprachverstehende Computersysteme, Wetterprognose und Multiplikation von Polynomen n-ten Grades gemeinsam? All diese Fachgebiete - und viele andere - wären heute ohne die schnelle Fourier-Transformation kaum mehr vorstellbar.
Spätestens seit dem Physikunterricht wissen wir, daß wir jede Schwingung als aus mehr oder weniger, verschieden frequenten Sinusschwingungen zusammengesetzt betrachten dürfen. Welchen Werte hätte die Erkenntnis aber, gäbe es nicht auch eine Methode, diese Frequenzanteile - und damit das Frequenzspektrum - zu errechnen: Die Fourier-Transformation.
Eine Funktion wird in eine Funktion der Frequenz transformiert. Dabei ist A(f) die Funktion der Frequenz, a(t) eine Funktion der Zeit, e die Basis des natürlichen Logarithmus, i die imaginäre Einheit die Frequenz und t die Zeit. Mit der inversen Fouriertransformation kann man eine Frequenzfunktion wieder in die Zeitfunktion zurücktransformieren:
Natürlich kann man jede beliebige Funktion bzw. Reihe transformieren, so wird zum Beispiel in der Computertomographie mittels Fourier- und inverser Fouriertransformation ein Schnittbild aus vielen Röntgenprojektionen rekonstruiert [3].
Hiiilfe!
Feider laden die beiden Formeln nicht besonders zum Nachprogrammieren ein, doch wir haben Glück, es gibt auch eine computergerechte Schreibweise, die diskrete Fouriertransformation. Sie transformiert genau N Meßwerte in ein Spektrum von N Koeffizienten:
bzw. invers:
Wenn wir nun noch berücksichtigen, daß e(ix) = cos(x) + i*sin(x), wenn wir das Rechnen mit komplexen Zahlen beherrschen und viel, viel Zeit haben, können wir uns bereits ans Werk machen. Feider dauert die Berechnung eines Satzes von 300 Meßwerten mit einem ST 12 Minuten. Nun benötigen wir aber z. B. in der Computertomographie viele, viele dieser Transformationen (üblicherweise 900je Schnittbild) und wir wollen unserem Patienten nicht zumuten, mehr als 7 Tage für ein Bild im Röntgenstrahl zu schmachten, daher suchen wir eine ganz, ganz schnelle Routine.
Andere haben schon vor uns gesucht und herausgefunden, daß die Koeffizienten der diskreten Fourier-Transformation der Lösung eines Polynoms (N-1 )-ten Grades für die komplexen Wurzeln der imaginären Einheit entsprechend ?). Für jedes N gibt es genau N komplexe Wurzeln (z) der imaginären Einheit, für die gilt: zN= 1. Wir müssen also ,nur noch ein simples Polynom4 lösen.
Tragen Sie es mit Fassung, die mathematischen Grundlagen sollen uns nicht weiter kümmern. Was für uns zählt, sind Ergebnisse.
Das Programm
Das vorgestellte Programm berechnet fortlaufend Sinus- oder Rechteckfolgen unterschiedlicher Frequenz oder liest die Werte wahlweise aus der parallelen Schnittstelle (Sampler!, [5]), berechnet daraus die Fourierkoeffizienten und aus diesen wieder die Ausgangswerte (inverse Fouriertransformation). Jeder Zwischenschritt wird graphisch angezeigt. Die Rekonstruktion ist allerdings so genau, daß man die Ausgangswerte (Punkte) kaum je zu Gesicht bekommt. Der verwendete Algorithmus (frei nach [2]) arbeitet rekursiv, N kann daher nur 2-er-Potenzen annehmen.
Um die Ausgabe einerseits schnell und flackerfrei, andererseits aber GEM-kompatibel zu halten, wird das Bild in einem virtuellen Bildschirm aufgebaut und erst dann auf den logischen Bildschirm kopiert. Die Vorzüge dieser Technik werden erst so richtig bei einer geschwindigkeitsoptimierten Programmversion erkennbar. Übrigens verwenden die Assemblerroutinen kein Clipping, wer experimentieren möchte sollte daher die Koordinaten vor dem Aufruf prüfen lassen. Die Funktion dft() ist nur zur Veranschaulichung im Listing enthalten, niemand würde jemals ernsthaft davon Gebrauch machen.
Das beschriebene Programm ist zwar schnell und genau, von Echtzeitspektralanalyse kann aber noch keine Rede sein. Diese erhält man, wenn man statt mit Double-Long mit Fixkomma-Long- oder Fixkomma-Integer-Werten rechnet und nur mehr das reine Spektrum anzeigt. Die aufwendigen Double-Rechnungen können so durch wenige Maschinenbefehle ersetzt werden. Auch wäre ein iterativer Algorithmus etwas schneller. Die lange Wartezeit am Anfang könnte durch Laden einer Tabelle der komplexen Wurzeln vermieden werden. Natürlich gibt es noch eine Menge weiterer Optimierungsmöglichkeiten, doch diese - Sie erraten es schon -würden wieder einmal den Rahmen des Artikels sprengen. (Dieser Satz steht immer dann, wenn es interessant wird oder wenn ein Autor selbst nicht mehr weiter weiß! - Ich versichere Ihnen, ersteres ist der Fall... [siehe Monats-Disk])
Die Interpretation
Die Fourier-Transformation liefert komplexe Zahlen. Der Betrag der Koeffizienten wird daher durch die pythagoräische Addition von Real- und Imaginärteil gebildet: Betrag = sqrt(Realteil2 + Imaginärteil2). Da die obere Hälfte der Werte der gespiegelten unteren Hälfte entspricht und daher redundant ist, betrachten wir nur die untere. Es bleiben uns also N/2 Werte. Bei einer Abtastrate von F Hz sind nur Frequenzen von 0 bis F/2 Hz darstellbar. Jeder der N/2 Koeffizienten einer Transformation mit N Werten und einer Abtastrate von F entspricht daher einem Frequenzbereich von (F/2)/(N/2)=F/N [Hz]. Bei Signalen, die nicht auf den Nullwert abgeglichen sind, wird dem Gleichspannungsanteil entsprechend zusätzlich eine ,Schwingung4 mit Frequenz 0 angezeigt.
Ein Beispiel
Werfen wir nun einen Blick auf das Bild: Es zeigt 512 Meßwerte einer Rechteckschwingung, eingelesen mit einer Abtastrate von etwa 50060 Hz, und das berechnete Frequenzspektrum. Einem Balken des Spektrums entsprechen hier 50060/5J2 = 97,9 Hz. Der stärkste Ausschlag befindet sich an Position 11, Nach unserer Rechnung sollte die Rechteckschwingung daher eine Grundfrequenz von / / *50060/512=1075 Hz haben. Tatsächlich - wer hätte das gedacht? - waren es 1066 Hz. Bei einer Auflösung von F/N müssen wir natürlich mit einer Toleranz von +/- F/(2N) rechnen. Wenn wir noch einen Blick in unsere Schulbücher werfen, können wir uns auch die anderen Schwingungen erklären. Eine Rechteckfunktion ist (nach Fourier) definiert durch
was wir also noch sehen, sind Schwingungen am 3fachen, 5fa-chen, 7fachen, etc. der Grundfrequenz mit einem Drittel, Fünftel, Siebtel, etc. ihrer Intensität (Amplitude).
Hardware-Projekt
Jetzt wollen wir uns einen Computertomographen basteln. Zuerst stellen wir eine preiswerte Röntgenröhre her: Öffnen Sie dazu Ihren Monitor und schaben Sie mit einem Messer die Graphitschicht auf einer Seite der Bildröhre ab (vorher Netzstecker ziehen!)...
Literatur:
[1] Henri Nussbaumer, Fast Fourier Transform and Convolution Algorithms, Springer Verlag, 2. Aufl., 1982.
[2] Robert Sedgewick, Algorithms, Addison Wesley Publishing Company, 2. Aufl., 1988.
[3] Vinash C. Kalk, Malcolm Slaney, Principles of Computerized Tomographie Imaging, IEEE Press.
[4] Christian Nieher, ATARI ST, Programmieren in Maschinensprache, Sybex- Verlag.
[5] Martin Backschat, Sample mir’s noch einmal, Sam!, Soundsampler im Selbstbau, ST-Computer, 7/88, S. 22 ff.
fft.prg
=
tcstart.o
fft_c.c
fft_s.s
tcfltlib.lib ; floating point library
tcstdlib.lib ; Standard library
tcextlib.lib ; extended library
tctoslib.lib ; TOS library
tcgemlib.lib ; AES and VDI library
; FFT.S
; (c) MAXON Computer 1993
EXPORT set_scrn
EXPORT scrn_line
EXPORT scrn_bar
EXPORT cntrcs
;--------------------------------------------------
; void set_scrn(void *adr, int wid)
;
; setzt Adresse adr für scr_line()
; und Subscreen-Breite für
; scrn_bar() und scrn_line()
set_scrn: move d0,linwid
move.l a0,scrnad
rts
;--------------------------------------------------
; void scrn_bar(void *addr, int x, int hgt)
;
; zeichnet Balken in subscreen
; addr: Basisadresse (linke untere Ecke)
; x: x-Koordinate (x Pixel nach rechts)
; hgt: Höhe des zu zeichnenden Balkens - 1
scrn_bar: move d0,d2
lsr #3,d0
lea 0(a0,d0.w),a0
andi.b #7,d2
move.b #%10000000,d0
lsr.b d2,d0
move linwid,d2
neg d2
barnext: or.b d0,(a0)
lea 0(a0,d2.w),a0
dbf d1,barnext
rts
;--------------------------------------------------
; void scrn_line(int x1, int y1, int x2, int y2)
;
; zeichnet Linie in subscreen
; nach: SYBEX; NIEBER; "ATARI ST,
; Programmieren in Maschinensprache"
scrn_line:
movem.l d3-d7,-(a7)
move 24(a7),d3
sf minus
move #$8000,d4
emp d1,d3
bhi noexch
exg d0,d2
exg d1,d3
noexch: sub d0,d2
bpi notneg
neg d2
st minus
notneg: sub d1,d3
move d0,d6
andi #7,d6
eori #7,d6
move d0, d7
lsr #3,d7
move dl,d5
mulu linwid,d5
add d5,d7
ext.l d7
add.l scrnad,d7
movea.l d7,a0
emp d2, d3
bhi county
bne countx
move #$ffff,d3
bra x_loop
countx: move #16,d5
lsl.l d5,d3
divu d2, d3
x_loop: bset d6,(a0)
tst.b minus
bne x_sub
x_add: subq #l,d6
bpi countx_1
move #7,d6
addq.l #l,a0
bra countx_1
x_sub: addq #l,d6
cmpi #8,d6
bne countx_1
clr d6
subq.l #l,a0
countx_1: add d3,d4
bc.c countx_2
adda linwid,a0
countx_2 : dbf d2,x_loop
bra lexit
county: move #16,d5
lsl.l d5,d2
divu d3,d2
y_loop: bset d6,(a0)
adda linwid,a0
add d2,d4
bcc county_2
tst.b minus
bne ysub
y_add: subq #l,d6
bpi county_2
move #7,d6
addq.l #l,a0
bra county_2
y_sub: addq #l,d6
cmpi #8,d6
bne county_2
clr d6
subq.l #1,a0
county_2: dbf d3,y_loop
lexit: movem.l (a7)+,d3-d7
rts
;--------------------------------------------------
; void cntrcs(signed char v[], int num)
;
; schreibt (num-1) Werte aus der parallelen
; Schnittstelle in das signed char array v
; (Wertebereich -128 bis +127)
cntrcs: movem.l d0/a0,cntpar
pea _cntrcs
move #38,-(a7) ; Supexec
trap #14
addq.l #6,a7
rts
_cntrcs: ori #$700,sr
movem.l cntpar, d0/a0
movea.l #$ff8800,a1
move.b #7,(a1)
move.b (a1),d1
andi.b #$7f,d1
move.b dl,2(a1)
_cntlp: move.b #$e,(a1) ; strobe high
move.b (a1),d1
ori.b #$20,d1
move.b dl,2(a1)
move.b #$f,(a1)
move.b (a1),d1 ; Scannerwert
sub.b #128,d1 ; - 128
move.b dl,(a0)+ ; ins Array
move.b #$e,(a1) ; strobe low
move.b (a1),d1
andi.b #$df,d1
move.b dl,2(a1)
dbf d0,_cntlp
andi #$fbff,sr
rts
BSS
minus: ds.w 1 ; Flag für rechts-links
scrnad: ds.l 1 ; Subscreen-Adresse
linwid: ds.w 1 ; Zeilenbreite Subscreen
cntpar: ds.l 2 ; Parameterzwischenspeicher
/* FFT.C
(c) MAXON Computer 1993
*/
#include <aes.h>
#include <float.h>
#include <math.h>
#include <stddef.h>
#include <tos.h>
#include <vdi.h>
#define ESCAPE 0x011b
#define FALSE (0)
#define TRDE (!FALSE)
#define ADDR(a) (int)((long)(a) » 16), \(int)((long)(a) & 0xffff)
#define SETPXY(p, a, b, c, d, e, f, g, h) \
p[0] = (a), p[l] = (b),\
p[2] = (c), p[3] = (d),\
p[4] = (e), p[5] = (f),\
p[6] = (g), p[7] = (h)
typedef struct { double r, /* Realteil */
double i; /* Xmaginärteil */
} complex;
#define fft(c, n) _fft(c, n, 0)
void ifft(complex [], int);
void _fft(complex [], int, int);
void _ifft(complex [], int, int);
void wfft(int);
void wifft(int);
void set_scrn(void *, int);
void scrn_line(int, int, int, int);
void scrn_bar(void *, int, int);
void cntrcs(signed char [], int);
int open_vwork(void) ;
void close_vwork(void);
void process(void);
/*----------------------------------------------*/
/* EXP beeinflußt die Fensterbreite und Breite */
/* der Transformation (N), daher auch die */
/* Geschwindigkeit; vernünftige Werte 7, 8, 9 */
#define EXP 9 /* Zweierexponent für N */
/*----------------------------------------------*/
/* globale Variablen und Konstanten, die */
/* in den Routinen für FFT und iFFT */
/* verwendet werden */
#define N (1 << EXP) /* N = 2^EXP */
complex w[N], /* komplexe Wurzeln für FFT */
wi[N], /* komplexe Wurzeln für iFFT */
t[N]; /* buf für FFT u. iFFT */
/*----------------------------------------------*/
/*----------------------------------------------*/
/* Macros für komplexe Operationen */
/* */
/* VAL: Betrag einer komplexen Zahl */
/* ASS: Zuweisung einer komplexen Zahl */
/* CMX: Bildung einer komplexen Zahl */
/* ADD: Addition zweier komplexer Zahlen */
/* SUB: Subtraktion zweier komplexer Zahlen */
/* MUL: Multiplikation zweier komplexer Zahlen */
#define VAL(c) (sqrt((c).r * (c).r\
+ (c).i * (c),i))
#define ASS(c, x) ((c).r = (x).r),\
((c).i = (x).i)
#define CMX(c, rp, ip) ((c).r = (rp)),\
((c).i = (ip))
#define ADD(c, x, y) ((c).r = (x).r + (y).r),\
((c).i = (x).i + (y).i)
#define SUB(c, x, y) ((c).r = (x).r - (y).r),\
((c).i = (x).i - (y).i)
#define MUL(c, x, y) ((c).r = (x).r * (y).r\
- (x).i * (y).i), \
((c).i = (x).r * (y).i\
+ (x).i * (y).r)
/*-----------------------------------------------*/
/*-----------------------------------------------*/
/* Programmspezif. Konstanten und Variablen */
#define FRQH 70 /* Höhe des Spektrums */
#define WORKW N /* Fensterbreite */
#define WORKH (280 + FRQH) /* Fensterhöhe */
complex a[N]; /* in & out für FFT & iFFT */
signed char vector[N];
/* log. und virtueller Bildschirm */
MFDB screen = {NULL, 0, 0, 0, 1, 0, 0, 0, 0},
sub = {NULL, WORKW, WORKH, WORKW >> 4,
1, 1, 0, 0, 0};
const char *string[] =
{"Fourier-Transformation",
"[l][Out of memory!][ Abort ]",
" [1][Illegal resolution!] [ Abort ]"};
int work_in[12], work_out[57], handle,
phys_handle, gl_hchar, gl_wchar, gl_hbox,
gl_wbox, gl_apid;
int open_vwork(void)
{ int i;
if ((gl_apid = appl_init()) != -1)
{ for (i = 1; i < 10; work_in[i++] = 1);
work_in[10] =2;
phys_handle = grafhandle(&gl_wchar,
&gl_hchar, &gl_wbox, &gl_hbox);
work_in[0] = handle = phys_handle;
v_opnvwk(work_in, &handle, work_out);
return TRUE;
}
else
return FALSE;
}
void close_vwork(void)
{ v_clsvwk(handle);
appl_exit();
}
int main(void)
{ if (open_vwork() == FALSE)
return -1;
process();
close_vwork();
return 0;
}
void process(void)
{ void *base;
int winx, winy, winw, winh, window,
scrx, scry, scrw, scrh, workx, worky, i,
d, key, message[8], event,
copy[8], delete[8], invert[8];
/* Speicher für subscreen */
if ((sub.fd_addr = Malloc((WORKW >> 3)
* WORKH)) == NULL)
{ form_alert(l, string[1]);
return;
}
wind_get(0, WF_WORKXYWH, &scrx, &scry, &scrw, &scrh);
workx = scrx + (scrw - WORKW) / 2;
worky = scry + (scry + scrh - WORKH) / 2;
SETPXY(delete, 0, 0, WORKW - 1, WORKH - 1,
0, 0, WORKW - 1, WORKH - 1);
SETPXY(copy, 0, 0, WORKW - 1, WORKH - 1,
workx, worky, workx + WORKW - 1,
worky + WORKH - 1);
SETPXY(invert, 0, WORKH - FRQH, WORKW - 1,
WORKH - 1, 0, WORKH - FRQH, WORKW - 1,
WORKH - 1);
/* Assemblerroutinen vorbereiten */
set_scrn(sub.fd_addr, WORKW » 3);
base = (void *)((long)sub.fd_addr
+ (WORKW >> 3) * (WORKH - 3));
if ((window = wind_create(NAME, scrx, scry, scrw, scrh)) >= 0)
{ wind_set(window, WF_NAME, ADDR(string[0]));
wind_calc(WC_BORDER, NAME,
workx, worky, WORKW, WORKH,
&winx, &winy, &winw, &winh);
if (winx < 0 || winw > work_out[0]
|| winy < 0 || winh > work_out[1])
{ graf_mouse(ARROW, NULL);
Mfree(sub.fd_addr);
wind_delete(window);
form_alert(1, string[2]);
return;
}
/*-----------------------------------------------*/
/* Initialisierung der komplexen Wurzeln */
/* der imag. Einheit (i) */
wfft(N);
wifft(N);
graf_mouse(M_OFF, NULL);
wind_open(window, winx, winy, winw, winh);
vro_cpyfm(handle, ALL_WHITE,copy, &sub, &screen);
while (1)
{ event = evnt_multi (MU_TIMER|MU__KEYBD,
0, 0, 0,
/* 10 mal 0 */ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
message, 5, 0,
&d, &d, &.d, &d, &key, &d);
if (event & MU_KEYBD)
if (key == ESCAPE)
break;
/*-----------------------------------------------*/
/* subscreen löschen */
vro_cpyfm(handle, ALL_WHITE, delete, &sub, &sub);
/*-----------------------------------------------*/
/* Ausgangskurve berechnen (Imaginärteile =0) */
/* Sinus */
/* { static k = 2;
**
** for (i = 0; i < N; i++)
** CMX(a[i], 80 *
** sin(2 * M_PI * i * k / N), 0);
** k += 4;
** if (k > N/2)
** k = 2;
** }
*/
/* oder Rechteck */
/* { static int k = 1;
** for (i = 0; i < N; i++)
** CMX(a[i],
** (i >> k) & 1 ? -80 : 80, 0);
**
** if (k++ > 8)
** k = 1;
** }
*/
/* oder vom Sampler abholen */
cntrcs(vector, N - 1);
for (i = 0; i < N; i++)
CMX(a[i], vector[i], 0);
/*-----------------------------------------------*/
/* Ausgangskurve zeichnen (dots) */
for (i = 0; i < N; i++)
scrn_line(i,
(WORKH - FRQH) / 2 - 1 - a[i].r, i,
(WORKH - FRQH) / 2 - 1 - a[i].r);
/*-----------------------------------------------*/
/* Spektrum berechnen und zeichnen */
fft(a, N);
for (i = 0; i < N / 2; i++)
scrn_bar(base, i * 2,
((int)VAL(a[i])) / N) & Oxff);
/* Betrag skaliert mit 1/N; Pseudoclipping */
/*----------------------------------------------*/
/* Rücktransformieren und Kurve zeichnen */
/* (dot joiner) */
ifft(a, N);
{ int x1 = 0, x2 = 1, y1, y2;
y1 = (WORKH - FRQH) / 2 - 1 - (int)a[0].r;
for (i = 1; i < N; i++)
{ y2 = (WORKH - FRQH) / 2 - 1 - (int)a[i].r;
scrn_line(x1++, y1, x2++, y2);
y1 = y2;
}
}
/*----------------------------------------------*/
/* Teil des subscreen invertieren (Spektrum) */
/* subscreen auf den Bildschirm kopieren */
vro_cpyfm(handle, D_INVERT, invert, &sub, &sub);
vro_cpyfm(handle, S_ONLY, copy, &sub, ascreen);
} /* while */
wind_close(window);
wind_delete(window);
graf_mouse(M_ON, NULL);
}
Mfree(sub.fd_addr);
return;
}
/*----------------------------------------------*/
/* berechnet in w[] die Wurzeln der imaginären */
/* Einheit für die FFT */
void wfft(int n)
{ double d;
int j;
for (j = 0; j < n; j++)
{ d = (2 * M_PI * j) / n;
CMX(w[j], cos(d), sin(d));
} /* Winkel in RAD! */
}
/*-----------------------------------------------*/
/* Fourier-Transformation; rekursiv; in-place */
/* p: Eingangs- und zugl. Ausgangsarray */
/* n: Transformation über n Werte */
/* k: Anfangs immer 0; wird erst in der */
/* Rekursion verändert; */
void _fft(complex p[], int n, int k)
{ int i, j, ndiv2 = n >> 1;
complex a, b;
if (n == 2)
{ ASS (b, p[k] ) ;
ASS (a, p[k + 1] ) ;
ADD(p[k], b, a);
SUB(p[k + 1], b, a);
}
else
{ for (i = 0; i < ndiv2; i++)
{ j = k + 2 * i;
ASS(t[i], p[j]);
ASS(t[i + ndiv2], p[j + 1]);
}
for (i = 0; i < n; i++)
ASS(p[k + i], t[i]);
_fft(p, ndiv2, k) ;
_fft(p, ndiv2, k + ndiv2);
j = N / n;
for (i = 0; i < ndiv2; i++)
{ MUL(b, w[i * j], p[k + ndiv2 + i]);
ADD(t[i], p[k + i], b);
SUB(t[i + ndiv2], p[k + i], b);
}
for (i = 0; i < n; i++)
ASS (p[k + i], t[i]);
}
}
/*-----------------------------------------------*/
/* berechnet in wi[] die Wurzeln der imag. */
/* Einheit für die iFFT */
void wifft(int n)
{ double d;
int j;
for (j = 0; j < n; j++)
{ d = (-2 * M_PI * j) / n;
CMX(wi[j], cos(d), sin(d));
} /* Winkel in RAD! */
}
/*-----------------------------------------------*/
/* Parameter wie _fft() */
/* Aufruf der iFFT; Ergebnis muß anschließend */
/* durch n dividiert werden; siehe Formel 4 */
void ifft(complex p[3, int n)
{ int i;
_ifft(p, n, 0);
for (i =0; i < n; i++)
p[i].r /= n, p[i].i /= n;
}
/* weitgehend identisch mit .. fft(); */
/* => einfach kopieren und ändern! */
void _ifft(complex p[], int n, int k)
{ int i, j, ndiv2 = n >> 1;
complex a, b;
if (n == 2)
{ ASS(b, p[k]);
ASS(a, p[k + 1]);
ADD(p[k], b, a);
SUB(p[k + 1], b, a);
}
else
{ for (i = 0; i < ndiv2; i++)
{ j = k + 2 * i;
ASS(t[i], p[j]);
ASS(t[i + ndiv2], p[j + 1]);
}
for (i = 0; i < n; i++)
ASS(p[k + i] , t [i]);
_ifft(p, ndiv2, k);
_ifft(p, ndiv2, k + ndiv2);
j = N / n;
for (i = 0; i < ndiv2; i++)
{ MUL(b, wi[i * j], p[k + ndiv2 + i]);
ADD(t[i], p[k + i], b);
SUB(t[i + ndiv2], p[k + i], b);
}
for (i = 0; i < n; i++)
ASS(p[k + i], t [i]);
}
}
/***********************************************/
/* NICHT EINTIPPEN! */
/***********************************************/
/* diskrete Fourier-Transform. (Formel 3) */
/* (leicht verständliche Funktionsweise */
/* aber unzumutbare Arbeitsgeschwindigkeit) */
void dft(complex a[], complex A[], int n)
{ int j, k;
double pi2ijk_div_N;
complex e, m, summe;
for (j = 0; j < n; j++)
{ CMX(summe, 0, 0);
for (k = 0; k < n; k++)
{ pi2ijk_div_N = 2 * M_PI * j * k / n;
CMX(e, cos(pi2ijk_div_N), sin(pi2ijk_div_N));
MUL(m, a[k], e);
ADD(summe, summe, m);
}
ASS(A[j], summe);
}
}
Melchior Franz