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



Links

Copyright-Bestimmungen: siehe Über diese Seite
Classic Computer Magazines
[ Join Now | Ring Hub | Random | << Prev | Next >> ]