DiskretFourier FH Augsburg Best Ever

5/17/2018 DiskretFourier FH Augsburg Best Ever - slidepdf.com

http://slidepdf.com/reader/full/diskretfourier-fh-augsburg-best-ever 1/8

Prof. Dr. Elmar Müller-Horsche1FH Augsburg Diskrete FourierTransformationen

Bedeutung der diskreten und schnellen Fouriertransformationen (DFT und FFT)

FFT’s sind die am häufigsten gebrauchten mathematischenAlgorithmen des Alltags. Myriaden solcher Transformationen werdenzum Erstellen komprimierter Musik (mp3) oder komprimierter Bilder

(jpeg) durchgeführt. Bei der Schwingungsanalyse gehören sie zumAlltag 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 manschreiben:

*

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 Achseentsteht.Beim Potenzieren solcher Zahlen werden nur die Phasen vervielfacht.z2 hat z.B. die Phase 2x, (z*)3 die Phase –3x. Die Punkte drehen sichalso auf dem Einheitskreis beim Potenzieren.

Komplexe Fourierkoeffizienten

0 0.5 1 1.5 2

x/2Pi

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

0

1

2

0

0

( ) ( 2 ) ( cos sin )2

cos1( ) 0,1.. 0

sin

n n

n

n

n

a f x f x n a nx b nx

a nx f x dx n b

b nx

(2)

Aus ihnen kann man folgendermaßen komplexe Fourierkoeffizienten ckonstruieren (nicht um die Studenten zu ärgern, sondern weil dieFormeln viel ergonomischer werden und man sich viel Schreib- undRechenarbeit erspart!). In einem komplexen c sind dabei die beidenreellen Zahlen a und b zusammengefasst:

2(2)

0

2(1)

0

1: ( ) (cos sin )

2 2

1( ) 0,1,2...2

wgn n

n

wg

jnx

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 dieIntegrale in (2) oder (3) analytisch nicht mehr ausgewertet werdenkönnen, oder wenn von vorneherein die Funktionswerte nur punktweise vorliegen, also z.B. bei elektronischer Datenerfassung. DasIntegrationsintervall 0 .. 2 wird wie gewohnt in N gleich breite Streifenunterteilt, 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 zumüssen, nennen wir das Ergebnis der numerischen Integrationwiederum c, behalten aber im Hinterkopf, dass es sich um eineNä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 k N

n

k k k

N jnx jnx

n k

k

N jn k N

k

k k

k

x k y f x k N N

c f x e x f x e dx

N N yc ee y

(4)

Die blaue Formel in (4) stellt genau die Diskrete Fouriertransformation(DFT) dar. Mithilfe des sog. „Schnörkelfaktors“ (twiddle factor im

Englischen)2

: j

N w e kann die DFT in Matrixschreibweise dargestelltwerden:

2

4

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 N N N

nnk

c y

w w wc y

c yw w w N

c yw w w

c y N

F w

F

Fkü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 Nz.B. würde gemäß 2 1 N k j k w e für alle Spalten k lauter Einsenenthalten genau wie Zeilenindex 0, Zeilenindex N+1 ergäbe dieselbenZahlen wie Zeilenindex 1 usw. Entsprechend würde sich imSpaltenvektor 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 = 89

10

11

12

13

14

N -1 = 15

In Polardiagrammen dargestellt ist oben die Matrix F mit N=16. Zeile 0sowie Spalte 0 bestehen immer aus lauter Einsen. Der Schnörkelfaktor w befindet sich immer an Indexposition 1,1 (rot). Mit ihm ergibt sichZeile 1 durch Rotation um jeweils den roten Winkel (2 /16 = 22,5° imBeispiel). Die gesamte Zeile 2 ergibt sich durch Winkelverdopplung vonZeile 1, Zeile 3 durch Winkelverdreifachung von Zeile 1 usw. Manerkennt: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 diereelle Achse gespiegelt. Dem entspricht die mathematische Operationkonjugiert komplex, wie eingangs erwähnt.Die gesamte Matrix ist außerdem symmetrisch (transponierte MatrixFT= F).In Formeln:

Spiegelsymmetriebzgl. mittlerer Zeile

Die c’s sind damitauch „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 k

k k

nk k n

k

F F

c F y F y

N N

F y c N

(6)




Eigenschaften des Vektors c der komplexen Fourier-Koeffizienten

Bsp.: abfallender SägezahnAmplitude = 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

c N

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 diekomplexen Komponenten von c auch wieder als Polardiagramme

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

c0 = a0/2 DC-AnteilcN/2 Oberwellenrest

Alle Informationen stecken also bereits in c0 ... cN/2. Der Rest desVektors ist das konjugiert Komplexe der oberen Hälfte.

Der Index n gibt aber die Frequenz der zugehörigen Oberwelle an. Diefett gedruckte Erkenntnis spiegelt damit schon ein wichtiges Theoremin der Signaltechnik wieder: Wenn ich N äquidistante Punkte einesperiodischen Signals kenne, so kann ich nur Amplituden und Phasenbis 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 Mikrofonsignalmindestens 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

Amplitudenspektrumdes obigenSä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 zuverwenden:

0

1

0

1

2 2

( ) ( cos sin )2

cos( )2

arctan

n n

n

n n

n

nn n n n

n

a f x a nx b nx

a A nx

b A a ba

(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... 12

n n n n

N N

N A c c n

A c cDC-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 denc’s und man kann aus ihnen sofort die wichtigen Amplituden- undPhasenspektren 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 relativeinfach im Wesentlichen durch konjugiert komplex Bildung jedesElementes von F:

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

1 1(5)

0 0

2 2 21 1 ( )

0 0

* ( *) N N wg

n i i k ni ik nk

i i

n i i k i N N j 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 angewendetwerden mit dem Resultat:

geometrische Summe:

1

0

1

1

N N i

i

aa

a

2( )2 2 ( )1 ( )

2 2( ) ( )0

1 10

1 1

0

i j n k N j n k N N j n k

N nk

j n k j n k i N N

nk

e ee

e e

N 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 / 2

1

1 1

:

( *)

*

N N wg j n k n k N

k n n

n n

N N wg j n k j N n k j k

N N N n n N

n

N c y y y y

N N N

y 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

wg j k n k j N n k j n k j x n j k N N N

wg j n k j N n k j x n j x n N N

n n n n

wg j 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 cos2 2

k n j x n

n k n k

N wg j k j k k

N k

ea nx b nx

j

N N e e k k x N




Trigonometrische Interpolation

Setzt man die Vereinfachungen in (11) ein, so erhält man die Formelfür die sog. trigonometrische Interpolation:

1

20 / 2

1

1(8) 2

0 / 2

1

cos sin cos2

cos( ) cos2

N

k n k n k N k

n

N

wg

n k n N k

n

N y c a nx b nx c x

N A 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 zueiner Gesamtfunktion g(x) überlagern, dass g(x) exakt durch allevorgegebenen Wertepaare (xk,yk) durchgeht. Solche Funktionen nenntman 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-Phasen

exakte Phasen (immer 90°

nur Sinuswellen)

12

0 / 2

1

(12)

( ) : cos( ) cos2

( )

N

n n N

n

wg

k k

N g x A A nx A x

g x y

Man beachte den Unterschied zur Fourierreihe: diese konvergiertgegen die zu analysierende Funktion f(x) (meistens) im gesamtenBereich 0 x 2 . Bricht man die Fourierreihe zur Fouriersumme ab, soerhält man eine Näherung, die zwar nahe bei f(x) liegt, mit dieser aber

nur nicht vorhersehbare Schnittpunkte gemeinsam hat. Dietrigonometrische Interpolation g(x) dagegen geht durch alleStü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 denexakten Fourierkoeffizienten (2) überein, wenn man Beispieleheranzieht, 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ützs tellen

nach der Oberwelle n=8 abgebrochene

Fourierreihe

trigonometrische Interpolation N=16




Divide et Impera: Fast Fourier Transform (FFT)

Teile und herrsche, dann fällt das herrschen leichter, meinten schondie Römer. Gauß erkannte 1805 die Möglichkeit der Arbeitsersparnisbei der Berechnung der DFT, Cooley und Tukey veröffentlichten dann

1965 den legendären Algorithmus, der ihren Namen trägt. Siespezialisierten sich auf den in der Praxis häufig auftretenden undschon anfangs erwähnten Spezialfall N=2

m. Die Berechnung der

Matrixmultiplikation1

c y N

F in (5) erfordert grob N 2 komplexe

Multiplikationen und Additionen, also je über 1Mio. Operationen beiN=1024. Man kann die F(N=1024)-Multiplikation aber auf dieF(N=512)-Muliplikation zurückführen, wie die folgende Rechnung zeigt.Dies spart ungeheuer! Im weiteren wird bei der Matrix F und beimSchnörkelfaktor w in Klammern dahinter geschrieben, für welcheEinteilung 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 kund der ungeraden in (4):

Beispiel: N=8 2 2 22 11 2 (2 1)

2 2 1

0 0

2 222 1 2 12 2

2 2 1

0 0

2 1 2 1

2 2 1

0 0

2 2

N N j n k j n m j n m N N N

n k m m

k m

N N j n m j n m j n N N N

m m

m m

N N n m n n m

m m

m 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 nn

w N w N sowie 2( ) ( )

N n

nw N w N zu

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

0 12

N n einschränken und schreiben:

2 1 2 1

2 2 1

0 02

2 ( ) 2 N N n m n mn n

m m N n m m

Ncw N y w N w N y Nc

N/2-DFT N/2-DFT

2.Stufe: gemischte Addition

(13)

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

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

werden addiert

0

2

..

N-2

F(N/2)

1

3

..

N-1

F(N/2)

Nc0

Nc1

..

Nc N/2-1

Nc N/2

Nc N/2+1 ..

Nc N-1

1. Stufe 2. Stufe

w =1

w1

...

w N/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




Mehrstufige FFT’s : Libellen und Schmetterlinge

Libelle1 1 1 11 11 1 1 1

1 1

(4)

a ab b j jc c

j jd d

F

Schmetterling1 1

(2)1 1

a a

b bF

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 nachdemselben Schema in F(N/4)-Multiplikationen aufzuspalten usw. Zubeachten ist dabei, dass bei jeder neuen Stufe die Eingangswerte yk umgeordnet werden müssen. Im 3-stufigen Fall wäre die Reihenfolgevon 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ührtenDiagrammstil 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üssendie zugehörigen Schnörkelfaktoren 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östeSchema:

Die dicken schwarzen Striche bedeuten wieder ·(-1), die dicken roten·(-j). Die Knoten wurden der Übersichtlichkeit halber weggelassen, beiden Potenzen der Schnörkelfaktoren stehen aus demselben Grund nur die Phasen in den Kästchen. Man sieht, dass bei kompletter Zerlegungrunter bis F(2) die Eingangswerte letztlich in sog. „bit reverse order“angeordnet werden müssen, die Binärzähler für die Indizes also vonlinks nach rechts arbeiten müssen. Der Deutlichkeit halber ist auf der linken Seite deshalb auch noch das Bitmuster der Indizes mitaufgefü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)

DiskretFourier FH Augsburg Best Ever

Documents

Transcript of DiskretFourier FH Augsburg Best Ever