Prof. Dr. Elmar Müller -Horsche FH Augsburg 1 Diskrete...

Prof. Dr. Elmar Müller-Horsche 1 FH Augsburg Diskrete FourierTransformationen

Bedeutung der diskreten und schnellen Fouriertransformationen (DFT und FFT)

FFT’s sind die am häufigsten gebrauchten mathematischen

Algorithmen des Alltags. Myriaden solcher Transformationen werden zum Erstellen komprimierter Musik (mp3) oder komprimierter Bilder (jpeg) durchgeführt. Bei der Schwingungsanalyse gehören sie zum Alltag des Messingenieurs („Spektralanalyse“).

Komplexe Zahlen (Wiederholung)

werden als Punkte in der komplexen Zahlenebene dargestellt. Hier sind vor allem komplexe Zahlen mit Betrag 1 interessant. Sie liegen auf dem Einheitskreis. Ihre Phasen nennen wir jetzt mal x. Das ist der Winkel zur Rechtsachse hin (= reelle Achse). Nach Euler kann man schreiben:

*

cos sin

cos sin

jx

jx

z x j x e

z x j x e−

= + =

= − = (1)

z* ist dabei die zu z konjugiert komplexe Zahl, die durch Vorzeichen-

wechsel beim Imaginärteil oder durch Spiegelung an der reellen Achse entsteht. Beim Potenzieren solcher Zahlen werden nur die Phasen vervielfacht. z2 hat z.B. die Phase 2x, (z*)3 die Phase –3x. Die Punkte drehen sich also auf dem Einheitskreis beim Potenzieren.

Komplexe Fourierkoeffizienten

0 0.5 1 1.5 2

x/2Pi

In den trigonometrischen Reihen reeller 2π-periodischer Funktionen kommen nur reelle a- und b-Koeffizienten vor:

0

1

2

00

( ) ( 2 ) ( cos sin )2

cos1( ) 0,1.. 0

sin

n nn

n

n

af x f x n a nx b nx

a nxf x dx n b

b nx

π

π

π

∞

== ± ⋅ = + +

= = ∞ =

∑

∫ (2)

Aus ihnen kann man folgendermaßen komplexe Fourierkoeffizienten c

konstruieren (nicht um die Studenten zu ärgern, sondern weil die Formeln viel ergonomischer werden und man sich viel Schreib- und Rechenarbeit erspart!). In einem komplexen c sind dabei die beiden reellen Zahlen a und b zusammengefasst:

2(2)

0

2(1)

0

1: ( ) (cos sin )

2 2

1( ) 0,1,2...

2

wgn n

n

wgjnx

a j bc f x nx j nx dx

f x e dx n

π

π

π

π−

− ⋅= = ⋅ − =

= ⋅ = ∞

∫

∫

(3)

x 2x

z

z2

z*

z*3

-3x


Diskrete Fouriertransformation

wird die numerische Integration von Integral (3) mit Hilfe der Treppenstufenmethode genannt. Dies wird notwendig, wenn die Integrale in (2) oder (3) analytisch nicht mehr ausgewertet werden können, oder wenn von vorneherein die Funktionswerte nur punktweise vorliegen, also z.B. bei elektronischer Datenerfassung. Das Integrationsintervall 0 .. 2π wird wie gewohnt in N gleich breite Streifen unterteilt, die Funktionswerte am linken Streifenrand sollen yk heißen, k = 0 .. N-1. N hat oft den Wert 512, 1024 usw., lässt sich also als 2er-Potenz schreiben: N = 2m. Um nicht neue Buchstaben einführen zu müssen, nennen wir das Ergebnis der numerischen Integration wiederum c, behalten aber im Hinterkopf, dass es sich um eine Näherung der komplexen Fourierkoeffizienten handelt.

21

0 0

21

0

21

0

2: : ( ) 0,1.. 1

1 1( ) ( )

2 1

2 2

1

2

k

N j n kN

n

k k k

Njnx jnx

n kk

N jn kN

kkk

k

x k y f x k NN

c f x e x f x e dx

N Nyc ee y

π

π π

π

π π

ππ

−− −

=

− − ⋅

=

− − ⋅

=

= ⋅ = = −

= ⋅ ⋅∆ ≅ ⋅

= ⋅ ⋅ = ⋅

∑ ∫

∑ ∑

(4)

Die blaue Formel in (4) stellt genau die Diskrete Fouriertransformation (DFT) dar. Mithilfe des sog. „Drehfaktors“

2

:j

Nw eπ− ⋅

= kann die DFT in Matrixschreibweise dargestellt werden:

24

1 1 1 1

1 1

1 1 1 1

1 1

Beispiel N=4

j

j j

j j

w e jπ− ⋅

− −

− −

− −

= = −

=

F

0 01 2 1

1 1

2 4 2 22 2

1 2 2 ( 1) ( 1)1 1

.1 1 1 1

.11

.1

. .. . . ..

1 .

1

N

N

N N N NN N

nnk

c y

w w wc y

c yw w wN

c yw w w

c yN

F w

−

−

− − − ⋅ −− −

= ⋅

= ⋅ ⋅

=

F

F� ��

kürzer:

Matrix hat dabei die komplexen Komponenten:

, 0,1.. 1k n k N⋅ = −

(5)

Eigentlich müsste in (5) der Index n bei den c’s bis unendlich laufen (s. Formel 3). Beim Versuch, die Matrix F nach unten zu verlängern, stellt man aber fest, dass sich die Zeilen wiederholen. Zeilenindex N z.B. würde gemäß 2 1N k j kw e π⋅ − ⋅ ⋅= = für alle Spalten k lauter Einsen enthalten genau wie Zeilenindex 0, Zeilenindex N+1 ergäbe dieselben Zahlen wie Zeilenindex 1 usw. Entsprechend würde sich im Spaltenvektor c die Zahlenfolge ab cN wiederholen. Man erkennt: Aus N reellen Funktionswerten yk (k=0,1 .. N-1) lassen sich nur N verschiedene komplexe Fourierkoeffizienten cn (n=0,1 .. N-1) mittels DFT bestimmen


Eigenschaften der DFT-Matrix F (N geradzahlig)

Zeile 0

1

2

3

4

5

6

7

N/2 = 8

9

10

11

12

13

14

N -1 = 15

In Polardiagrammen dargestellt ist oben die Matrix F mit N=16. Zeile 0

sowie Spalte 0 bestehen immer aus lauter Einsen. Der Drehfaktor w befindet sich immer an Indexposition 1,1 (rot). Mit ihm ergibt sich Zeile 1 durch Rotation um jeweils den roten Winkel (2π/16 = 22,5° im Beispiel). Die gesamte Zeile 2 ergibt sich durch Winkelverdopplung von Zeile 1, Zeile 3 durch Winkelverdreifachung von Zeile 1 usw. Man erkennt: In Zeile cN/2 stehen 1, -1 im Wechsel. Der untere Rest (Zeile 9 .. 15) der Matrix ergibt sich durch Spiegelung der Zeilen 7 .. 1 an dieser Zeile 8. Die einzelnen Polardiagramme werden dabei auch um die reelle Achse gespiegelt. Dem entspricht die mathematische Operation konjugiert komplex, wie eingangs erwähnt. Die gesamte Matrix ist außerdem symmetrisch (transponierte Matrix FT= F). In Formeln:

Spiegelsymmetrie bzgl. mittlerer Zeile

Die c’s sind damit

auch „symmetrisch“ bezüglich cN/2 !

k

( )

*

( ) ( ) ( )

(y und N reell!)

*

1 1* *

1

N n k nk

N n N n k k N n k kk k

nk k nk

F F

c F y F yN N

F y cN

−

− − −

=

= ⋅ = ⋅ =

= ⋅ =

∑ ∑

∑

(6)


Eigenschaften des Vektors c der komplexen Fourier-Koeffizienten

Bsp.: abfallender Sägezahn

Amplitude = 16, N = 16

16

8

0

8

16

16 1

14 1 5, 03

12 1 2, 41

10 1 1, 50

8 1

6 1 0, 67

4 1 0, 41

2 1 0, 2

0 1

2 1 0, 2

4 1 0, 41

6 1 0, 67

8 1

10 1 1, 50

12 1 2, 41

14 1 5, 03

1

j

j

j

j

j

j

j

j

j

j

j

j

j

j

cN

−−−

−−−−

− +− +− +− +

− +− +− +

= ⋅ ⋅ =

F�

Wie in der 2. Zeile von (6) gezeigt, ist auch c�

wie F „spiegel-symmetrisch“ bezüglich der mittleren Komponente cN/2 (sofern man die komplexen Komponenten von c

�

auch wieder als Polardiagramme darstellen würde) , also c(N-n)* = cn. Verwendet wurde dabei die komplexe Rechenregel: (u·v·w)* = u*·v*·w* (u, v, w komplex) sowie r* = r für reelles r. Außerdem müssen c0 und cN/2 reell werden, da die entsprechenden Zeilen in F nur aus reellen Zahlen bestehen ( 0 2,1; 1k N kF F= = ± ). Sie heißen:

c0 = a0/2 DC-Anteil cN/2 Oberwellenrest Alle Informationen stecken also bereits in c0 ... cN/2. Der Rest des Vektors ist das konjugiert Komplexe der oberen Hälfte. Der Index n gibt aber die Frequenz der zugehörigen Oberwelle an. Die fett gedruckte Erkenntnis spiegelt damit schon ein wichtiges Theorem in der Signaltechnik wieder: Wenn ich N äquidistante Punkte eines periodischen Signals kenne, so kann ich nur Amplituden und Phasen bis zur Oberwelle mit N/2-facher Frequenz der Grundwelle bestimmen („Nyquist-Theorem“). Anwendung: ich will Musik bis 20kHz analysieren (= höchster hörbarer Ton), dann muss ich das Mikrofonsignal mindestens mit 40 kHz abtasten, also alle 25µs das Signal digi-talisieren (ist bei hochwertigen Audio-wav-files auch ungefähr der Fall).

(7)

Gesamtamplituden und Phasen

Amplitudenspektrum des obigen Sägezahns

0 50

5

10

15

Phasenspektrum

0 2 4 6 80

30

60

90

Bei den reellen Fourierkoeffizienten gibt es die Möglichkeit, Gesamtamplituden An und Phasen ϕn anstelle von an und bn zu verwenden:

0

1

0

1

2 2

( ) ( cos sin )2

cos( )2

arctan

n nn

n nn

nn n n n

n

af x a nx b nx

aA nx

bA a b

a

ϕ

ϕ

∞

=∞

=

= + +

= + −

= + =

∑

∑ (8)

Die letzte Zeile lässt sich aber auch leicht direkt aus den cn bestimmen, wie ein Blick auf die Definition (3) zeigt. Zusammengefasst mit (7) ergibt sich schließlich:

0 0 2 2

2 arg 1... 12n n n n

N N

NA c c n

A c c

ϕ= = − = −

= =DC-Anteil A Oberwellenrest

(9)

Für den Praktiker sind diese Größen meist interessanter als die a- und b-Koeffizienten. Sie ergeben sich gemäß (9) direkt und einfach aus den c’s und man kann aus ihnen sofort die wichtigen Amplituden- und Phasenspektren ableiten.


Inverse diskrete Fouriertransformation (IDFT): wie bekomme ich yk wieder aus cn?

Beispiel N=4

1 1 1 11 11 1 1 11 1

1 1 1 1 4 0 0 01 1 0 4 0 01 1 1 1 0 0 4 01 1 0 0 0 4

j j

j j

j j

j j

− −− −

− −

− −− −− −

Mithilfe der Inversen Matrix F-1 natürlich. Und die ergibt sich relativ einfach im Wesentlichen durch konjugiert komplex Bildung jedes Elementes von F:

F*·F = N·1 (10) wobei F* aus den Elementen Fnk* besteht. Beweis:

( )1 1(5)

0 0

2 2 21 1 ( )

0 0

* ( *)N Nwg

n i i kni iknk

i i

n i i k iN Nj j j n k

N N N

i i

F F w w

e e eπ π π

− −⋅ ⋅

= =⋅ ⋅

− −⋅ − ⋅ ⋅ −

= =

⋅ = ⋅ = ⋅ =

= ⋅ =

∑ ∑

∑ ∑

F * F

Auf der Diagonalen ist n=k und die Summe besteht aus lauter Einsen. Für n≠k muss die Formel für geometrische Summen angewendet werden mit dem Resultat:

geometrische Summe:

1

0

1

1

NNi

i

aa

a

−

=

−=−∑

( )

( ) {

2( )2 2 ( )1 ( )

2 2( ) ( )0

1 10

1 1

0

i j n k Nj n kN Nj n k

Nnk

j n k j n ki N N

nk

e ee

e eN für n k

für n k

ππ π

π π

⋅ − ⋅ ⋅ −− ⋅ −

⋅ − ⋅ −=

− −⋅ = = = = − −

=⋅ = ≠

∑F * F

F * Fzusammengefasst :

Mit diesen Resultaten wollen wir Gleichung (5) von links mit F* multiplizieren und erhalten:

21 1(5)

0 0

1 2 2 2(6) 2 ( )2

0 / 21

1 1

:

( *)

*

N Nwg j n kn k N

k n nn n

NNwg j n k j N n k j k

N N Nn n N

n

Nc y y y y

N N Ny c oder

y w c e c

c e c e c e c

π

π π π

− − ⋅ ⋅ ⋅⋅

= =

−⋅ ⋅ ⋅ ⋅ ⋅ − ⋅ ⋅ ⋅ ⋅

=

⋅ = ⋅ ⋅ ⋅ = ⋅ ⋅ ⋅ = ⋅ =

= ⋅

= ⋅ = ⋅ =

= + ⋅ + ⋅ + ⋅

∑ ∑

∑

F * F * F F * F

F *

� � � � �

� �

(11)

In der letzten Zeile wurde die Summe in 2 Hälften aufgespalten (z.B. N=16: 1...7 und 15...9) und die Symmetrie (6) der cn wurde ausgenutzt (z.B. c3 = c13* oder c13 = c3*). (11) kann weiter vereinfacht werden:

22 2 (4)2( ) 2

2 2 (4)( )

(3)

1

* *

2 2

2

k

k k

k k

k k k

wgj k n kj N n k j n k j x nj kNN N

wgj n k j N n k j x n j x nN Nn n n n

wgj x n j x nn n n n

j x n j x n j x

n n

e e e e e

e c e c c e c e

a jb a jbe e

e e ea b

ππ πππ

π π

⋅ − ⋅ ⋅⋅ ⋅ − ⋅ − ⋅ ⋅ ⋅ − ⋅ ⋅⋅

⋅ ⋅ ⋅ ⋅ ⋅ − ⋅ ⋅ ⋅ − ⋅ ⋅

⋅ ⋅ − ⋅ ⋅

⋅ ⋅ − ⋅ ⋅ ⋅

= = ⋅ = ⋅

⋅ + ⋅ = ⋅ + ⋅ =− += ⋅ + ⋅ =

+= +

2 (4)2

cos sin2

2( 1) cos cos cos

2 2

kn j x n

n k n k

N wgj kj k kN

k

ea nx b nx

j

N Ne e k k x

N

ππ π π

⋅ − ⋅ ⋅

⋅ ⋅ ⋅ ⋅ ⋅

− = +

= = − = = ⋅ ⋅ = ⋅


Trigonometrische Interpolation

Setzt man die Vereinfachungen in (11) ein, so erhält man die Formel

für die sog. trigonometrische Interpolation:

( )1

2

0 / 21

1(8) 2

0 / 21

cos sin cos2

cos( ) cos2

N

k n k n k N kn

Nwg

n k n N kn

Ny c a nx b nx c x

NA A nx A xϕ

−

=

−

=

= + + + ⋅ ⋅

= + − + ⋅ ⋅

∑

∑

(12)

Mithilfe der aus der DFT gewonnenen Koeffizienten cn oder (an,bn) oder (An,ϕn) kann man also Grundwelle cos x und Oberwellen cos nx so zu einer Gesamtfunktion g(x) überlagern, dass g(x) exakt durch alle vorgegebenen Wertepaare (xk,yk) durchgeht. Solche Funktionen nennt man Interpolationsfunktionen.

Beispiel: Sägezahn (N=16)

0 5 10 150

5

10

15

genäherte DFT-Amplituden (N=16)exakte Fourierkoeffizienten

0 5 10 150

30

60

90

genäherte DFT-Phasenexakte Phasen (immer 90° nur Sinuswellen)

12

0 / 21

(12)

( ) : cos( ) cos2

( )

N

n n Nn

wg

k k

Ng x A A nx A x

g x y

ϕ−

=

= + − + ⋅

=

∑

Man beachte den Unterschied zur Fourierreihe: diese konvergiert

gegen die zu analysierende Funktion f(x) (meistens) im gesamten Bereich 0≤x≤2π. Bricht man die Fourierreihe zur Fouriersumme ab, so erhält man eine Näherung, die zwar nahe bei f(x) liegt, mit dieser aber nur nicht vorhersehbare Schnittpunkte gemeinsam hat. Die trigonometrische Interpolation g(x) dagegen geht durch alle Stützstellen (xk,yk) exakt durch, die den Dateninput der DFT bildeten. Für große N spielen diese Unterschiede allerdings keine Rolle mehr, die Amplituden und Phasen aus der DFT stimmen praktisch mit den exakten Fourierkoeffizienten (2) überein, wenn man Beispiele heranzieht, bei denen man (2) auch analytisch exakt bestimmen kann.

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 120

16

12

8

4

0

4

8

12

16

20

x/2Pi

zu analysierende Funktion f(x)Stützstellennach der Oberwelle n=8 abgebrochene Fourierreihetrigonometrische Interpolation N=16


Divide et Impera: Fast Fourier Transform (FFT) Teile und herrsche, dann fällt das herrschen leichter, meinten schon

die Römer. Gauß erkannte 1805 die Möglichkeit der Arbeitsersparnis bei der Berechnung der DFT, Cooley und Tukey veröffentlichten dann 1965 den legendären Algorithmus, der ihren Namen trägt. Sie spezialisierten sich auf den in der Praxis häufig auftretenden und schon anfangs erwähnten Spezialfall N=2m. Die Berechnung der

Matrixmultiplikation 1

c yN

= ⋅ ⋅F� ��

in (5) erfordert grob N2 komplexe

Multiplikationen und Additionen, also je über 1Mio. Operationen bei N=1024. Man kann die F(N=1024)-Multiplikation aber auf die F(N=512)-Multiplikation zurückführen, wie die folgende Rechnung zeigt. Dies spart ungeheuer! Im Weiteren wird bei der Matrix F und beim Schnörkelfaktor w in Klammern dahinter geschrieben, für welche Einteilung N die Größen gelten. Also F(1024) ist eine 1024x1024-

Matrix,2

1024(1024)j

w eπ− ⋅

= . Der Trick besteht lediglich im Zusammenfassen der geraden Indizes k und der ungeraden in (4):

Beispiel: N=8

( ) ( ) ( )

2 2 22 11 2 (2 1)

2 2 10 0

2 222 1 2 12 2

2 2 10 0

2 1 2 1

2 2 10 0

2 2

NN j n k j n m j n mN N N

n k m mk m

N Nj n m j n mj nN NN

m mm m

N Nn m n n m

m mm m

Nc e y e y e y

e y e e y

w N y w N w N y

π π π

π ππ

−− − ⋅ − ⋅ − ⋅ ++

= =

− −− ⋅ − ⋅−

+= =

− −⋅ ⋅+

= =

= ⋅ = ⋅ + ⋅ =

= ⋅ + ⋅ ⋅ =

= ⋅ + ⋅ ⋅

∑ ∑

∑ ∑

∑ ∑

In der letzten Zeile stehen aber im Wesentlichen zwei N/2-DFT’s, jedenfalls solange n<N/2. Für n≥N/2 ist die Periodizität

2( 2) ( 2)N

nnw N w N

+= sowie 2( ) ( )

Nn n

w N w N+

= − zu berücksichtigen (s. Torten links). Man kann dann n auf den Bereich

0 12Nn≤ ≤ − einschränken und schreiben:

( ) { } ( )

2 1 2 1

2 2 10 02

2 ( ) 2N Nn m n mn n

m mNn m m

Ncw N y w N w N yNc

− −⋅ ⋅

++ = =

+ = ⋅ ⋅ ⋅−

∑ ∑N/2-DFT N/2-DFT

2.Stufe: gemischte Addition

��

��

(13)

Damit ist das Ziel erreicht. Man erhält eine 2-stufige DFT. In der 1. Stufe wird mit den geraden y2m und den ungeraden y2m+1 je eine N/2-DFT durchgeführt. In der 2. Stufe erfolgt die gemischte Addition gemäß (13). Am besten macht man sich das Verfahren in einem Diagramm deutlich. In der 2.Stufe bedeuten:

Dünner Strich = Weiterleitung der Zahl nach rechts Dicker Strich = Multiplikation mit –1 und Weiterleitung Kasten = Multiplikation mit Zahl im Kasten und weiter Alle von links in einen Knoten einlaufende Striche

werden addiert

y0

y2

..

yN-2

F(N/2)

y1

y3

..

yN-1

F(N/2)

Nc0

Nc1

..

NcN/2-1

NcN/2

NcN/2+1

..

NcN-1

1. Stufe 2. Stufe

w0=1

w1

...

wN/2-1

w=w(N)

w(N/2) = w(4)

w(4)7=w(4)3

w(8)7= - w(8)3

w(N) = w(8) w(8)3

Prof. Dr. Elmar Müller-Horsche

8 FH Augsburg Diskrete FourierTransformationen

Divide et Impera: Nochmal mit Bildern erklärt Die Matrixgleichung für die 16-Punkt-DFT:

16 c y F (input: Signalwerte yk, output: Spektralwerte cn)

ist im unteren Bild so dargestellt, dass die auf der rechten Seite auftretenden Produktpaare Fnk·yk schön untereinander stehen. Dazu wurde der Spaltenvektor y in einen Zeilenvektor nach links gekippt

und als Spaltenüberschrift für F missbraucht. Das Bild ist jetzt so zu lesen: um z.B. den Spektralwert 16·c6 zu erhalten (er enthält Informationen zur Oberwelle mit 6-facher Grundfrequenz), müssen die roten Wertepaare ausmultipliziert und addiert werden!

X

+

X

+

X

+

X

+

X

+

X

+

X

+

X

+

X

+

X

+

X

+

X

+

X

+

X

+

X

+

X



Sortieren nach geraden und ungeraden Indizes Die Umordnung der Spalten incl. Spaltenüberschriften ergibt dasselbe

Ergebnis, da Summen problemlos umgeordnet werden können. In der linken Hälfte ergibt sich schon eine wesentliche Vereinfachung: die linke obere Hälfte der Drehfaktoren entspricht der linken unteren Hälfte! Sie bestehen jeweils aus der kleineren 8x8 Matrix F(N=8), die

mit dem Drehfaktor ( 45 )je

gebildet wird. Damit können für die linke Hälfte die Ergebnisse von 16·c0 … 16·c7 hergenommen werden um 16·c8 … 16·c15 zu bestimmen! Mit den geraden y wird dazu eine 8-Punkt-DFT durchgeführt. Leider sieht die rechte Seite noch nicht so schön aus. Das ändert sich aber, wenn man in jeder Zeile einen komplexen Faktor auf dem Einheitskreis ausklammert. Dem entsprechen reine Drehungen, die „Uhren“ in den Zeilen der rechten Hälfte werden um jeweils einen bestimmten Betrag zurückgedreht…



Ausklammern von Multiplikatoren … und schon entstehen auch auf der rechten Seite F(N=8)-Blöcke,

also auch die ungeraden y werden einer 8-Punkt-DFT unterzogen! Allerdings müssen die entstehenden Werte noch mit dem angezeigten Multiplikator multipliziert werden, der für jede Zeile unterschiedlich ist. Die Multiplikatoren in der unteren Hälfte sind gerade das Negative der oberen Hälfte. Wenn man unten anstelle einer Addition eine Subtraktion einführt, werden die Multiplikatoren gleich…



Oben plus, unten minus Damit hat man den Rechenaufwand grob halbiert:

anstelle einer Multiplikation

(16 16) (16)xF y (grob 256 komplexe Multiplikationen und Additionen)

hat man 2 Multiplikationen der Art

(8 8) (8)x

gerade

ungerade

F y (grob 64 komplexe Multiplikationen und Additionen)

durchzuführen plus Multiplikation mit den Multiplikatoren und Addition / Subtraktion der linken und rechten Anteile



Schaltbild Führt man ähnlich wie in der Elektrotechnik Strichdiagramme ein mit

Input links und Output rechts, so entsteht letztendlich das untere Bild. Mit den geraden und ungeraden y wird jeweils eine 8-Punkt-DFT berechnet (oben / unten). Die unteren Ergebnisse müssen noch mit den Multiplikatoren verdreht werden. Anschließend erfolgt die Addition (dünner Strich) oder Subtraktion (dicker Strich). Um z.B. 16·c10 zu erhalten, muss man also den 3. Wert vom unteren Block vom 3. Wert des oberen Blockes subtrahieren


Mehrstufige FFT’s : Libellen und Schmetterlinge Libelle

1 1 1 11 11 1 1 11 1

(4)a ab bj jc c

j jd d

F − −− −

− −

⋅ = ⋅ =

Schmetterling

( ) ( ) ( )1 1(2)

1 1a ab b

⋅⋅ = − =F

Schema 13 halbiert ungefähr den Rechenaufwand bei großen DFT’s (z.B. N=1024). Es liegt natürlich nahe, die F(N/2)-Multiplikation nach demselben Schema in F(N/4)-Multiplikationen aufzuspalten usw. Zu beachten ist dabei, dass bei jeder neuen Stufe die Eingangswerte yk umgeordnet werden müssen. Im 3-stufigen Fall wäre die Reihenfolge von oben nach unten z.B.: y0, y4, y8 .. y2, y6, y10 .. y1, y5, y9 .. y3, y7, y11 .. Beim betrachteten Spezialfall N=2m gelangt man so nach m-1 Stufen zur F(4) und nach m Stufen zur F(2). Mit dem bei (13) eingeführten Diagrammstil kann man F(4) als „Libelle“ (dragonfly) und F(2) als „Schmetterling“ (butterfly) darstellen.

Bei jeder neuen Stufe, die auf der linken Seite dazu kommt, müssen die zugehörigen Drehfaktoren w bestimmt werden. Ihre Winkel verdoppeln sich jeweils von rechts nach links. So bekommt man für N=16 letztendlich folgendes völlig in Schmetterlinge aufgelöste Schema:

Die dicken schwarzen Striche bedeuten wieder ·(-1), die dicken roten ·(-j). Die Knoten wurden der Übersichtlichkeit halber weggelassen, bei den Potenzen der Dreh faktoren stehen aus demselben Grund nur die Phasen in den Kästchen. Man sieht, dass bei kompletter Zerlegung runter bis F(2) die Eingangswerte letztlich in sog. „bit reverse order“ angeordnet werden müssen, die Binärzähler für die Indizes also von links nach rechts arbeiten müssen. Der Deutlichkeit halber ist auf der linken Seite deshalb auch noch das Bitmuster der Indizes mit aufgeführt. Man hat ausgerechnet, dass bei völliger Zerlegung in Schmetterlinge die 1024-Punkte-DFT nur noch 0,4% der ursprünglichen Multiplikationen und nur noch 1% der ursprünglichen Additionen benötigt!

- 45°

-135°

- 45°

-135°

- 22,5°

- 45°

- 67,5°

-112,5°

- 135°

-157,5°

w(16) = e - j·22,5° w(8) = e - j·45°

w(4) = e - j·90°=-j

0000 y0

1000 y8

0100 y4

1100 y12

0010 y2

1010 y10

0110 y6

1110 y14

0001 y1

1001 y9

0101 y5

1101 y13

0011 y3

1011 y11

0111 y7

1111 y15

16·c0

16·c1

16·c2

16·c3

16·c4

16·c5

16·c6

16·c7

16·c8

16·c9

16·c10

16·c11

16·c12

16·c13

16·c14

16·c15

a

b

a

b

c

d

·1 ·(-1) · j ·(-j)

Prof. Dr. Elmar Müller -Horsche FH Augsburg 1 Diskrete...

Documents

Transcript of Prof. Dr. Elmar Müller -Horsche FH Augsburg 1 Diskrete...