Fourier-Analysis undder Algorithmus des 20. Jahrhunderts

Dipl.-Inform. Wolfgang Globke

Institut fur Algebra und GeometrieArbeitsgruppe Differentialgeometrie

Universitat Karlsruhe

1 / 54

Die schwingende Saite

Schlagt man die Saite einer Gitarre (oder eines anderenSaiteninstruments) an, so schwingt die Saite um ihreRuhelage.

Sie erzeugt dadurch einen Ton, dessen Tonhohe (Frequenz)von der Wellenlange der Schallwellen abhangt.

Die Wellenlange wiederum ist durch die Lange der Saitegegeben.

Verkurzt man die Wellenlange (z.B. durch Abgreifen in derMitte der Saite), so erhoht sich die Frequenz.

2 / 54

Die schwingende Saite

Lasst man die Saite frei schwingen, so hort man den Grundton.

0 1 2 3 4 5 6

-2

-1

1

2

Im Bild ist der Einfachheit wegen die Saitenlange π angenommenworden.Sind die Enden der Saite fest eingespannt, so nimmt dieSchwingung (fortgesetzt auf [0, 2π]) die Form sin(x) an.

3 / 54

Grund- und Obertone

Verkurzt man die Saite um den Faktor 2, 3, 4,. . ., so hort man den2., 3., 4.,. . . Oberton.

0 1 2 3 4 5 6

-2

-1

1

2

0 1 2 3 4 5 6

-2

-1

1

2

Der n-te Oberton nimmt die Form sin(nx) an.Sind die Enden der Saite nicht fest eingespannt, konnen auchSchwingungen der Form cos(nx) (oder Uberlagerungen) auftreten.

4 / 54

Superpositionsprinzip

Aus dem physikalischen Superpositionsprinzip ergibt sich, dasssich jede denkbare Schwingung als Uberlagerung der Grund-und Oberschwingungen darstellen lasst.

Genauer gesagt: Die Schwingung f nimmt an der Stelle x derSaite die Form

f (x) = β1 sin(x) + β2 sin(2x) + β3 sin(3x) + . . .

an.Hierbei ist βk die Amplitude, mit der die k-teOberschwingung in die Uberlagerung eingeht.

Wir stellen fest:Superposition bedeutet nichts anderes, als dass f eineLinearkombination der Grund- und Oberschwingungen ist.

5 / 54

Signale

Die Schwingungen einer Saite sind ein Spezialfall dessen, was imtechnischen Sprachgebrauch als Signal bezeichnet wird.Andere Beispiele fur Signale konnen etwa elektromagnetischeWellen (Funksignale) oder digitale Bilder sein.

6 / 54

Signale

Nach dem Superpositionsprinzip ist ein Signal eineLinearkombination gewisser

”elementarer Signale“

(z.B. die Grund- und Oberschwingungen).

7 / 54

Signale

Auch gewisse unstetige Signale lassen sich durch Superposition vonSchwingungen beliebig gut approximieren.Lassen wir statt Linearkombinationen auch

”unendliche

Linearkombinationen“ (also Reihen) zu, so konnen wir auch dieseSignale aus unseren elementaren Signalen zusammensetzen.

8 / 54

Signale

Sagezahnfunktion approximiert durch

− 2

π

2∑k=1

sin(πkx)

k.

0 1 2 3 4 5 6

-1

1

9 / 54

Signale

Sagezahnfunktion approximiert durch

− 2

π

7∑k=1

sin(πkx)

k.

0 1 2 3 4 5 6

-1

1

10 / 54

Signale

Sagezahnfunktion approximiert durch

− 2

π

15∑k=1

sin(πkx)

k.

0 1 2 3 4 5 6

-1

1

11 / 54

Signale

Wir erhalten die Sagezahnfunktion im Grenzwert:

− 2

π

∞∑k=1

sin(πkx)

k.

12 / 54

Signale

Kastenfunktion approximiert durch

4

π

2∑k=1

sin(π(k + 12 )x)

2k + 1.

0 1 2 3 4 5 6

-1

1

13 / 54

Signale

Kastenfunktion approximiert durch

4

π

7∑k=1

sin(π(k + 12 )x)

2k + 1.

0 1 2 3 4 5 6

-1

1

14 / 54

Signale

Kastenfunktion approximiert durch

4

π

15∑k=1

sin(π(k + 12 )x)

2k + 1.

0 1 2 3 4 5 6

-1

1

15 / 54

Signale

Wir erhalten die Kastenfunktion im Grenzwert:

4

π

∞∑k=1

sin(π(k + 12 )x)

2k + 1.

16 / 54

Analyse

Es stellt sich nun die Frage, wie man fur eine gegebenes Signal fbestimmen kann, wie es sich aus den elementaren Signalenzusammensetzt.

Wir wissen, dass man Funktionen als Linearkombinationen derelementaren Funktionen schreiben kann.

Wir suchen also einen geeigneten Vektorraum V , der alleSignale enthalt, die uns interessieren.

Wurden unsere Signale z.B. im Rn leben, so konnten wir siedurch einfache Rechnung in ihre Komponenten bzgl. derelementaren Signale zerlegen.

Aber man sieht schnell, dass V unendliche Dimension habenmuss.

17 / 54

Signale endlicher Energie

Der geeignete Vektorraum V ist der Raum L2[0, 2π] derSignale mit endlicher Energie auf dem Intervall [0, 2π]:

L2[0, 2π] =

f : [0, 2π]→ R

∣∣∣ ∫ 2π

0f (x)2dx <∞

/∼ .

Die Aquivalenzrelation”modulo ∼“ bedeutet hier nur, dass wir

Signale, die sich an hochstens abzahlbar vielen Stellenunterscheiden, als identisch betrachten.L2[0, 2π] enthalt unter anderem folgende wichtige Funktionen:

Die stetigen Funktionen auf [0, 2π].

Die stuckweise stetigen Funktionen auf [0, 2π].

18 / 54

Signale endlicher Energie

Auf L2[0, 2π] kann man das folgende L2-Skalarprodukt definieren:

〈f , g〉 =1

π

∫ 2π

0f (x)g(x)dx .

Man sieht, dass die Bedingung endlicher Energie gerade derEndlichkeit der L2-Norm entspricht:

‖f ‖2 =1

π

∫ 2π

0f (x)2dx <∞.

19 / 54

Eine Orthonormalbasis

Man erhalt nun den folgenden Satz:Die Menge der elementaren Funktionen 1√

2, cos(n·), sin(n·)

∣∣∣ n ∈ N

ist eine Orthonormalbasis von L2[0, 2π](allerdings nur in dem Sinne, dass wir auch unendliche Reihen alsLinearkombinationen zulassen).

20 / 54

Skalarprodukte in Aktion

Durch Berechnen des L2-Skalarproduktes konnen wir nun fur jedesSignal f endlicher Energie seine Komponente in Richtung derorthogonalen Funktionen ermitteln:

α0 = 〈f , 1/√

2〉, αk = 〈f , cos(k ·)〉, βk = 〈f , sin(k·)〉.

Jeder Wert αk bzw. βk stellt die Intensitat (Amplitude) dar,mit der das Signal cos(k ·) bzw. sin(k·) zum Signal f beitragt.

21 / 54

Fourier-Reihen

Somit konnen wir f durch diese ONB darstellen:

f (x) = α0 · 1√2

+∞∑

k=1

(αkcos(kx) + βksin(kx)

).

(”Gleichheit“ gilt hierbei nur modulo ∼.)

Diese Darstellung ist als Fourier-Reihe von f bekannt,die αk , βk ∈ R heißen Fourier-Koeffizienten.

22 / 54

Unendliche Signale

Es gibt auch auf R = (−∞,+∞) Signale mit endlicherEnergie, d.h. der Vektorraum L2(R) existiert.

Aber cos und sin sind nicht Elemente von L2(R).

Es gibt keine Fourier-Reihendarstellung fur Signale f : R→ C

endlicher Energie.

Die Fourier-Reihe wird hier durch das Fourier-Integral ersetzt,

f (ω) = 〈f , eiω·〉 =1√2π

∫R

f (t)e−iωtdt,

das man als”Grenzfall a→∞“ fur die Fourier-Reihen im

Intervall [−a, a] erhalt.

Statt abzahlbar vieler Fourier-Koeffizienten hat man hier einKontinuum von Frequenzen ω, die mit Intensitat f (ω) zu fbeitragen.

23 / 54

Komplexe Fourier-Reihen

Um die Darstellung zu vereinfachen, kann man die Sinus- undCosinusterme in der Fourier-Reihe auf dem Intervall [0, 2π] durchkomplexe e-Funktionen ersetzen:

f (t) =∞∑

k=−∞γkeikt ,

wobei die γk jetzt durch das hermitesche Skalarprodukt

γk = 〈f , eik·〉 =1

2π

∫ 2π

0f (t)e−iktdt

gegeben ist.Der Zusammenhang mit den ursprunglichen Koeffizienten ist

αk = γk + γ−k , βk = i(γk − γ−k).

24 / 54

Digitale Signale

0 0,8 1,6 2,4 3,2 4 4,8 5,6

-2,4

-1,6

-0,8

0,8

1,6

2,4

tt0 t1 tn!1

f

25 / 54

Digitale Signale

Um aus einem kontinuierlichen (analogen) Signal f eindigitales Signal zu machen, wird f an n aquidistant verteiltenZeitpunkten t0, . . . , tn−1 abgetastet (genauer: tj = j

n2π).

Das analoge Signal f wird dann durch ein stuckweisekonstantes (digitales) Signal (f0, . . . , fn−1) mit den Wertenfj := f (tj) approximiert.

26 / 54

Digitale Signale: Abtasttheorem

In der Praxis hat man oft bandbegrenzte Signale, d.h. es gehen inder Fourier-Reihe nur Oberschwingungen bis zu einer Ordnungk0 <∞ ein:

f (t) =

k0∑k=−k0

γkeikt

In diesem Fall besagt das Abtasttheorem, dass n = 2k0 + 1Abtaststellen ausreichen, um das analoge Signal f vollstandig zurekonstruieren.Dies sieht man leicht, da man 2k0 + 1 lineare Gleichungen

fj =

k0∑k=−k0

γkeiktj

benotigt, um 2k0 + 1 Unbekannte γk eindeutig zu bestimmen.

27 / 54

Diskrete Fourier-Transformation

Wir approximieren die Fourier-Koeffizienten mit Hilfe derAbtastwerte:

γk ≈ ck :=1

n

n−1∑j=0

fje−2πik j

n .

Wir werden sehen, dass diese Werte (c0, . . . , cn−1) ausreichen,um das digitale Signal (f0, . . . , fn−1) zu rekonstruieren. Dahernennen wir (c0, . . . , cn−1) die diskrete Fourier-Transformiertevon (f0, . . . , fn−1).

28 / 54

Diskrete Fourier-Transformation

Die diskrete Fourier-Transformation ist eine lineare Abbildung

Fn : Cn → Cn, f 7→ Ωn · f ,

mit f =

f0...

fn−1

und

Ωn =1

n

1 1 1 · · · 11 ω1 ω2 · · · ωn−1

1 ω2 ω4 · · · ω2(n−1)

......

.... . .

...

1 ωn−1 ω2(n−1) · · · ω(n−1)2

∈ Cn×n,

wobei wir abkurzend ω := e−2πin schreiben.

29 / 54

Diskrete Fourier-Transformation

Die Matrix Ωn ist invertierbar und es ist Ω−1n = nΩn.

Es folgtf = Ω−1

n Ωnf = Ω−1n c ,

wobei c der Vektor mit den diskreten Fourier-Koeffizienten ist.

Daraus folgt, dass das digitale Signal f durch die diskretenFourier-Koeffizienten festgelegt ist.

Bei der Berechnung der diskreten Fourier-Transformationdurch die Formel f 7→ Ωf werden O(n2) Additionen undMultiplikationen durchgefuhrt.

Bei großen n ist dieser Aufwand unakzeptabel.Ein verbessertes Verfahren muss her!

30 / 54

Flinke Fourier-Transformation (FFT)

Es gibt einen Algorithmus, der die Berechnung der diskretenFourier-Transformation erheblich beschleunigt:Die schnelle Fourier-Transformation (FFT).

Wir nehmen nun an, n sei eine Zweierpotenz n = 2p.

Dann konnen wir die spezielle Struktur der Matrix Ωn

ausnutzen:

Ωn =1

n

1 1 1 · · · 11 ω1 ω2 · · · ωn−1

1 ω2 ω4 · · · ω2(n−1)

......

.... . .

...

1 ωn−1 ω2(n−1) · · · ω(n−1)2

.

31 / 54

Flinke Fourier-Transformation (FFT)

Mit dieser Matrix schreibt sich der k-te diskreteFourier-Koeffizient

ck =1

n

n−1∑j=0

fjωjk .

Diese Summe kann man in zwei Teilsummen aufteilen: In derersten sammelt man alle Summanden mit geraden Indizes, inder zweiten alle Summanden mit ungeraden Indizes:

ck =1

n

n2−1∑

j=0

f2jω2jk +

1

n

n2−1∑

j=0

f2j+1ω(2j+1)k

=(1

n

n2−1∑

j=0

f2jω2jk)

+ ωk(1

n

n2−1∑

j=0

f2j+1ω2jk)

Die Formeln in den Klammern beschreiben die diskreteFourier-Transformation F n

2der Ordnung n

2 = 2p−1.

32 / 54

Flinke Fourier-Transformation (FFT)

Wir konnen also Fn(f ) berechnen, indem wir rekursiv F n2(f+)

und F n2(f−) berechnen. Hierbei sind

f+ =

f0f2...

fn−2

, f− =

f1f3...

fn−1

∈ C n2

die diskreten Signale, die jeweils die Koeffizienten mit geradenbzw. ungeraden Indizes von f enthalten.

Außerdem giltωn = 1, ω

n2 = −1.

nach Definition von ω.

33 / 54

Flinke Fourier-Transformation (FFT)

Damit schreiben wir die Koeffizienten der diskretenFourier-Transformierten c = Ωnf als

ck = (Ω n2f+)k + ωk · (Ω n

2f−)k

c n2

+k = (Ω n2f+)k − ωk · (Ω n

2f−)k

(1)

fur k = 0, . . . , n2 − 1.

In ubersichtlicher Matrixschreibweise:

c = Ωnf =

(E n

2W n

2

E n2−W n

2

)·(

Ω n2

0

0 Ω n2

)·(

f+f−

)=

(E n

2W n

2

E n2−W n

2

)·(

Ω n2

0

0 Ω n2

)· Pn · f ,

(2)

wobei W n2

= diag(1, ω, . . . , ωn2 ) ∈ C n

2× n

2

und Pn ∈ Cn×n ist eine Permutationsmatrix, die gerade undungerade Indizes von f entsprechend sortiert.

34 / 54

Flinke Fourier-Transformation (FFT)

Fur n = 22 = 4 z.B. ist

P4 =

1 0 0 00 0 1 00 1 0 00 0 0 1

·

f0f1f2f3

=

f0f2f1f3

.

Da n eine Zweierpotenz ist, lasst sich dieses Prinzip (1)rekursiv auch auf Ω n

2anwenden, bis schließlich bei Ω1 die

Rekursion endet.

Mit Hilfe der Matrixschreibweise (2) lasst sich die Rekursionin eine iterative Prozedur auflosen.

35 / 54

FFT Aufwandsanalyse

Um den Aufwand T (n) fur Fn(f ) zu ermitteln, betrachten wirwieder Gleichung (1):

ck = (Ω n2f+)k + ωk · (Ω n

2f−)k

c n2

+k = (Ω n2f+)k − ωk · (Ω n

2f−)k

fur k = 0, . . . , n2 − 1.

Hier geht zweimal der Aufwand T (n2 ) fur F n

2ein, und

außerdem n2 Multiplikationen mit den ωk :

T (n) = n2 + 2 · T (n

2 ).Mit Induktion konnen wir T (n

2 ) ≤ n4 log2(n

2 ) annehmen.Es folgt:

T (n) =n

2+ 2 · T

(n

2

)≤ n

2+ 2 ·

(n

4log2

(n

2

))=

n

2

(log2(2) + log2

(n

2

))=

n

2log2(n).

36 / 54

FFT Aufwandsanalyse

Der Aufwand fur die FFT liegt also in O(n log(n)).Dies ist eine enorme Verbesserung gegenuber O(n2):

0 80 160 240 320 400 480 560 640 720 800 880 960 1040 1120 1200

800

1600

2400

3200

O(n log(n))

O(n2)

37 / 54

Algorithmus des 20. Jahrhunderts

Die FFT wurde Mitte der 1960er Jahre entwickelt.

Durch sie wurde erstmals die effektive Verarbeitung (Echtzeit)großer Datenmengen moglich.

Sie wird uberall eingesetzt, wo in irgendeiner Weise digitaleSignale verarbeitet werden mussen (und daruber hinaus).

Wegen dieser Bedeutung bezeichnet man die FFT auch alsden Algorithmus des 20. Jahrhunderts.

38 / 54

Anwendung: Faltung von Signalen

Betrachten wir ein LTI-System H, dass fur ein digitalesEingabesignal f ein Ausgabesignal H(f ) liefert.

Hf H(f)

LTI steht dabei fur

Linear, d.h. H(αf1 + f2) = αH(f1) + H(f2).

Time Independent, d.h. fur die gleiche Eingabe f hangt dieAusgabe H(f ) nicht vom Zeitpunkt des Empfangs ab.

Außerdem verlangen wir noch, dass H stetig ist.

39 / 54

Anwendung: Faltung von Signalen

Jedes digitale Signal f = (f0, . . . , fn−1) lasst sich schreiben als

f =n−1∑k=0

fkδk ,

wobei δk = (0, . . . , 0, 1︸︷︷︸k-te Stelle

, 0, . . . , 0) sei.

Das LTI-System H ist vollstandig durch die Impulsantwort

H(δ0) =: h = (h0, . . . , hn−1)

bestimmt (es ist dann H(δk)j = hj−k).

40 / 54

Anwendung: Faltung von Signalen

Fur ein beliebiges Signal f gilt wegen der Linearitat von H:

H(f )j =n−1∑k=0

fkH(δk)j =n−1∑k=0

fkhj−k =: (f ∗ h)j .

Diesen Ausdruck bezeichnen wir als die diskrete Faltung f ∗ h vonf mit dem Faltungskern h.Faltungen lassen sich sehr effizient mit der FFT berechnen, denn esgilt der Faltungssatz

Fn(f ∗ h) = Fn(f ) Fn(h),

wobei das Symbol die koeffizientenweise Multiplikation derVektoren Fn(f ) und Fn(h) bedeutet.

41 / 54

Anwendung: Faltung von Signalen

Anstatt die Faltung zu berechnen, kann man also zunachst ftransformieren, dann mit Fn(h) bestimmen und dannzurucktransformieren:

Fn

!h

!Fn(h)

F!1n

f f ! h

Fn(f) Fn(f)! Fn(h)

42 / 54

Anwendung: Langzahlarithmetik

Es seien a, b ∈ Z ganze n-stellige Zahlen, die zu groß sind, umdirekt im Register des Prozessor multipliziert zu werden.Zur Basis β konnen wir a und b schreiben als

a =n−1∑k=0

akβk , b =

n−1∑k=0

bkβk .

Das Produkt c = a · b ist dann

c =2n−2∑k=0

ckβk ,

wobei ck durch die Faltung

ck =n−1∑j=0

ak−jbj

gegeben ist.43 / 54

Anwendung: Langzahlarithmetik

Angenommen, wir konnen auf unserem Rechner modulo mrechnen, und es ist m > ck fur alle k .

Dann konnen wir die Faltungen mit Hilfe der FFT berechnen,wenn wir ω durch eine Einheitswurzel in Z/mZ ersetzen.

Dieses Vorgehen ist als Schonhage-Strassen-Algorithmusbekannt.

44 / 54

Anwendung: JPEG

Das JPEG-Verfahren (Joint Photographics Expert Group) ist einverlustbehaftetes Kompressionsverfahren fur Bilder.

Wir fassen das Bild als Signal f (x , y) auf, wobei x und y dieKoordinaten der einzelnen Pixel beschreiben, und f (x , y) derFarbwert am Pixel (x , y) ist.

JPEG verwendet eine Variante der FFT,die diskrete Cosinus-Transformation (DCT).

Das Bild wird in 8× 8-Pixelblocke unterteilt, die jeweilsCosinus-transformiert werden:

f (u, v) = cu,v

7∑x=0

7∑y=0

f (x , y) cos

((2x + 1)uπ

16

)cos

((2y + 1)vπ

16

),

wobei die Konstanten cu,v nur der Normierung dienen.

45 / 54

Anwendung: JPEG

Dabei stellen die Werte f (u, v) die Intensitat dar, mit der dieWellenlangen u in x-Richtung bzw. v in y -Richtung in den8× 8-Block eingehen.

(u, v) = (0, 0) entspricht der niedrigsten Frequenz,(u, v) = (7, 7) der hochsten.

Mit Hilfe von empirisch ermittelten Tabellen werden dietransformierten Werte f (u, v) quantisiert, d.h. auf denWertebereich 0, . . . , 255 abgebildet:

Q(f (u, v)) = fQ(u, v) ∈ 0, . . . , 255.

Dabei werden die hohen Frequenzanteile unterdruckt.

46 / 54

Anwendung: JPEG

Die quantisierten Frequenzanteile fQ werden nun sequentiellgeordnet, beginnend mit den niedrigsten Frequenanteilen:

fQ(u, v)

u

v

Der so erhaltene Datenstrom wird nun mit derHuffman-Codierung komprimiert. Diese Daten werdenschließlich abgespeichert.Um das Bild zu rekonstruieren, werden alle Schritte ruckwartsdurchgefuhrt.Da bei der Quantisierung Information verloren gegangen ist,wird das Bild in der Regel nicht exakt rekonstruiert.

47 / 54

Anwendung: MP3

Das MP3-Format fur Audiodaten verwendet ein ahnlichesVerfahren wie die JPEG-Komprimierung, angewandt aufAudiodaten.

48 / 54

Wavelets

Die diskrete Fourier-Transformation eignet sich vor allem furdas Verarbeiten periodischer Signale.

Kunstliche erzeugte Signale sind oft periodisch, aber siewerden von unregelmaßigen Storungen uberlagert(ein kurzer Knall oder das Knacken auf einer Schallplatte).

Um solche kurzen, unregelmaßigen Storungen darzustellen, isteine Wavelet-Basis fur unseren Signalraum L2(R) besondersgut geeignet.

49 / 54

Wavelets

Es gibt viele Arten von Wavelets. Ihnen allen ist gemeinsam, dasses ein Mutter-Wavelet ψ ∈ L2(R) gibt, und die Wavelet-Basis Bentsteht durch Skalieren und Verschieben von ψ:

B = ψs,k | s ∈ N0, k ∈ Z,

wobei das Wavelet ψs,k definiert ist durch

ψs,k(x) = 2−s/2ψ(2−sx − k).

Der Parameter s, die Skala von ψs,k , ist ein Maß fur die Feinheitder Auflosung, die durch ψs,k moglich ist.

50 / 54

Wavelets

Beispiele fur Wavelets sind die Haar-Wavelets oder die MexicanHat-Wavelets:

-3 -2,5 -2 -1,5 -1 -0,5 0 0,5 1 1,5 2 2,5 3

-1

-0,5

0,5

1

!Mex

!Haar

51 / 54

Wavelets

Da Wavelets sehr stark an einer Stelle lokalisiert sind, eignen siesich gut, um lokale Storungen in einem Signal f herauszufiltern:

= +f = f ! "f ,!#!f !

52 / 54

Flinke Wavelet Transformation (FWT)

Die Skalenstruktur der Wavelets fuhrt auf einen Algorithmusum die diskrete Wavelet-Transformation zu berechnen.

Dieser Algorithmus FWT ist sogar noch schneller als die FFT;er liegt in O(n).

Da bei diesem Verfahren Koeffizienten fur die verschiedenenSkalen eines Signals schrittweise bestimmt werden, sprichtman von Multi-Skalen-Analyse.

53 / 54

D. BensonMusic: A Mathematical Offering (CUP)

O. ForsterAlgorithmische Zahlentheorie (Vieweg)

C. Gasquet, P. WitomskiFourier Analysis and Applications (Springer)

G. GoosVorlesungen uber Informatik, Band 4 (Springer)

T.W. KornerFourier Analysis (CUP)

A.K. Louis, P. Maaß, A. RiederWavelets (Teubner)

Fastest Fourier Transform in the West (FFTW)http://www.fftw.org/

54 / 54