Fourier-Reihen, Teil 6b – DFT gemessener Signale

Im letzten Teil haben wir die Fourier-Koeffizienten eines Signals s numerisch berechnet, unter der Voraussetzung, die Periodendauer des Signals zu kennen.

Wenn wir ein Signal messen, kennen wir dessen Periodendauer normalerweise nicht. Wir messen einfach während der Messdauer T_m mit der Sampling-Frequenz (Abtastrate) f_s die momentanen Werte s(t). Wie beeinflusst das die Fourier-Koeffizienten?

Abb. 1 zeigt nochmals unser Signal

s(t) = -1 + 3\sin(2\pi \cdot 0.5\,\text{Hz} \cdot t + \pi) + 2\sin(5 \cdot 2\pi \cdot 0.5\,\text{Hz} \cdot t - \tfrac{\pi}{2})

aus dem letzten Teil.

sampled_meas_sig
Abb. 1: Das Signal aus Teil 6. Innerhalb der Messdauer von 3.5 s ist das Signal dicker gezeichnet. Der hellblaue Verlauf ist die tatsächliche Periodizität, der hellrote Verlauf die scheinbare Periodizität. Die roten Punkte sind die 32 Messwerte.

Weil wir die Periodendauer T = 2\,\text{s} jetzt nicht kennen,  haben wir das Signal einfach eine beliebige Messdauer T_m, z.B. T_m = 3.5\,\text{s}, lang gemessen. Dabei haben wir N = 2^5 = 32 Werte aufgenommen (rote Punkte). Die Sampling-Frequenz ist also

\displaystyle f_s = \frac{1}{\Delta t} = \frac{N}{T_m} = \frac{32}{3.5\,\text{s}} \approx 9.1\,\text{Hz} ,

und hängt neben der Anzahl der Messwerte auch von der Messdauer ab.

scheinbare Periodizität mit der Messdauer

Der hellblaue Verlauf zeigt wie unser Signal tatsächlich periodisch weitergeht. Aufgrund unserer Messung glauben wir aber, dass unser gemessenes Signal sich mit der Messdauer periodisch wiederholt (hellrote Kurve). Warum ist das so?

Aus unseren Messwerten s_n berechnen wir die Fourier-Koeffizienten \underline{S}_k. Wie wir schon gesehen haben, sind letztere periodisch: \underline{S}_{k+N} = \underline{S}_k. Wenn wir daraus unser Signal zurückrechnen, erhalten wir ebenfalls eine Periodizität:

\displaystyle\begin{aligned}s_{n + N} &= \sum_{k = 0}^{N - 1} \underline{S}_k \cdot e^{\underline{i} \tau k (n + N) / N}\\ &= \sum_{k = 0}^{N - 1} \underline{S}_k \cdot e^{\underline{i} \tau k n / N} \cdot e^{\underline{i} \tau k}\\ &= \sum_{k = 0}^{N - 1} \underline{S}_k \cdot e^{\underline{i} \tau k n / N} = s_n \, ,\end{aligned}

weil e^{\underline{i} \tau k} = (e^{\underline{i} \tau})^k = 1^k = 1 ist. Wenn wir aber N zu unserem Zeitindex addieren, sind wir zeitlich um N \cdot \Delta t = N \cdot T_m / N = T_m weitergegangen:

s(t_n  + T_m) = s(t_n) .

Diese Periodizität ist nicht real, sondern ein Artefakt, weil unsere Messdauer kein ganzzahliges Vielfaches der tatsächlichen Periodendauer ist. In den meisten Fällen sorgt sie für einen Sprung im Signal am Ende der scheinbaren Periode.

Amplitudenspektrum

Mit dem folgenden MATLAB/Octave-Code können wir jetzt unser Amplitudenspektrum berechnen:

% tatsächliche Periodendauer (in s)
T = 2.0;

% Messdauer (in s)
Tm = 3.5;

% Anzahl der Samples
% FFT ist am effizientesten mit 2er-Potenzen
N = 2^5;

% Abtastzeiten (in s)
tn = (0:(N - 1))*Tm/N;

% das Signal
sn = -1 + 3*sin(2*pi/T*tn + pi) + ...
     2*sin(5*2*pi/T*tn + pi/2);

% die Frequenzen der Fourier-Koeffizienten (in Hz)
fk = (0:(N - 1))/Tm;

% Fourier-Koeffizienten berechnen
Sk = fft(sn)/N;

% komplexe Amplituden
Ak = 2i*Sk;
Ak(1) = Ak(1)/2; % MATLAB beginnt bei 1 zu zählen

% Absolutwerte der komplexen Amplituden
absAk = abs(Ak);

% bis 4 Hz zeichnen
kmax = floor(4*Tm) + 1;

plot(fk(1:kmax), absAk(1:kmax), 'o')
xlabel('Frequenz in Hz')
ylabel('Amplitude')

Das Ergebnis zeigt Abb. 2, und es ist – hmmm… Zunächst einmal sieht man große Amplituden nahe bei 0.5 Hz und 2.5 Hz, was gut ist. Allerdings sehen wir sie nicht bei exakt 0.5 Hz und exakt 2.5 Hz, und die Größe passt auch nicht. Zusätzlich gibt es bei vielen anderen Frequenzen auch relativ große Amplituden.

DFT_Amp_3_5s_32
Abb. 2: Amplitudenspektrum des Signals aus Abb. 1 mit einer Messdauer von 3.5 s und 32 Messwerten.

Warum passen die Frequenzen nicht? Weil unser gemessenes Signal jetzt periodisch mit der Messdauer T_m ist; daher ist die Frequenzauflösung

\Delta f = 1/T_m .

Für unseren Wert T_m = 3.5\,\text{s} von oben ergibt sich \Delta f = 1/3.5\,\text{s} \approx 0.286\,\text{Hz}. Alle anderen Frequenzen sind ganzzahlige Vielfache davon – für 0.5 Hz und 2.5 Hz geht sich das aber nicht aus.

OK, wenn wir \Delta f kleiner machen, indem wir die Messdauer T_m verlängern, können wir prinzipiell näher an unsere tatsächlichen Frequenzen herankommen. Vielleicht hilft das ja. Wenn wir das tun, müssen wir aber gleichzeitig die Anzahl der Samples erhöhen, weil die Sampling-Frequenz f_s = N/T_m sonst zu klein wird. Nehmen wir an, wir messen etwa 4-mal so lang, z.B. T_m = 13.1\,\text{s} und vervierfachen auch die Anzahl der Messwerte. Das Ergebnis zeigt Abb. 3.

DFT_Amp_13_1s_128
Abb. 3: Amplitudenspektrum des Signals aus Abb. 1 mit einer Messdauer von 13.1 s und 128 Messwerten.

Durch das kleinere Frequenzintervall lässt sich die Lage der beiden gesuchten Amplituden zwar besser eingrenzen, wirklich gut funktioniert es aber immer noch nicht. Dabei ist die 2.5 Hz-Amplitude besser getroffen als die 0.5 Hz-Amplitude. Der Grund ist, dass 33 \cdot 1/2.5\,\text{Hz} = 33 \cdot 0.4\,\text{s} = 13.2\,\text{s} ist, was fast mit unserer Messdauer übereinstimmt.

Sollte die Messdauer zufällig in der Nähe eines ganzzahligen Vielfachen der tatsächlichen Periodendauer sein, bei uns T = 2\,\text{s}, würden alle Amplituden besser herauskommen. Das zeigt Abb. 4, wo die Messdauer von 13.1 s auf 14.3 s erhöht wurde, und 14.3 s ist ungefähr 7-mal 2 s.

DFT_Amp_14_3s_128
Abb. 4: Amplitudenspektrum des Signals aus Abb. 1 mit einer Messdauer von 14.3 s und 128 Messwerten.

Eine weitere kleine Erhöhung der Messdauer auf 15 s verschlechtert das Ergebnis für beide Amplituden wieder (s. Abb. 5), weil 15 s weder ein ganzzahliges Vielfaches von 2 s noch von 0.4 s ist.

DFT_Amp_15s_128
Abb. 5: Amplitudenspektrum des Signals aus Abb. 1 mit einer Messdauer von 15 s und 128 Messwerten.

Exakt wird es erst wieder, wenn unsere Messdauer genau ein ganzzahliges Vielfaches der Periodendauer ist (s. Abb.5b). In allen anderen Fällen scheint es so, als würden die Amplituden in die benachbarten Frequenzen ausfließen.

DFT_Amp_16s_128.png
Abb. 5b: Amplitudenspektrum des Signals aus Abb. 1 mit einer Messdauer von 16 s und 128 Messwerten.

Leck-Effekt (Leakage)

Dieses »Ausfließen« der Amplituden zu benachbarten Frequenzen nennt man Leck-Effekt bzw. Leakage. Er hat nichts mit der diskreten Abtastung des Signals zu tun, sondern kommt durch die Sprünge am Rand des Messintervalls zustande.

Wie wir schon in Teil 1 bzw. Teil 2 gesehen haben, benötigt so ein Sprung viele kleine, schnell rotierende Zeiger, die sich zur richtigen Zeit zu einer Art »Peitschenschlag« addieren. Dies ist noch einmal in Abb. 6 anhand eines Cosinus mit Periodendauer 4/7 s gezeigt, der alle 2 s (entsprechend der Messdauer) wiederholt wird. Weil 2 s kein ganzzahliges Vielfaches von 4/7 s ist, kommt es zu einem deutlichen Sprung.

AnimCutCos
Abb. 6: Eine cos-Funktion mit Periodendauer 4/7 s, die alle 2 s periodisch wiederholt wird, wodurch es zu einem Sprung kommt. Die Animation ist um den Faktor 4 verlangsamt, um den Sprung besser sehen zu können.

Die Zeigersumme läuft die meiste Zeit entlang eines Kreises mit Radius 2, wie wir es für eine reine cos-Funktion mit Amplitude 2 erwarten würden. Allerdings haben wir nicht einen Zeiger, der mit der Eigenfrequenz f_0 = \frac{7}{4}\,\text{Hz} = 1.75\,\text{Hz} umläuft, sondern zwei fast gleich große, deren Frequenzen leicht unter- bzw. oberhalb der Eigenfrequenz liegen. In der Mitte der Messdauer, z.B. bei 1 s, ist der Summenzeiger praktisch die Summe dieser beiden Zeiger. Alle anderen Zeiger sind zu diesem Zeitpunkt »eingerollt«.

Wenn wir jedoch in die Nähe des Sprungs kommen, rollen sich die höherfrequenten Zeiger aus und erzeugen uns den Sprung. Wie wir schon bei der Sägezahnfunktion in Teil 2 gesehen haben, wird ihre Amplitude mit 1/f kleiner. Das zeigt auch das Amplitudenspektrum in Abb. 7. Wir haben etwas unter- und etwas oberhalb der Eigenfrequenz zwei große Amplituden. Je weiter wir zu Frequenzen höher als f_0 gehen, desto mehr nehmen die Amplituden wie f/f^2 = 1/f ab.

SpektCutCos
Abb. 7: Das Amplituden- und Phasenspektrum der Funktion aus Abb. 6. Die »Eigenfrequenz« des Cosinus ist f_0 = \frac{7}{4}\,\text{Hz} = 1.75\,\text{Hz}.

Aber warum haben wir auch Amplituden unterhalb der Eigenfrequenz? Wie gesagt, in der Mitte des Messintervalls sind praktisch nur die beiden größten Zeiger relevant. Zum Rand des Intervalls hin aber, wenn sich die hochfrequenten Zeiger auszurollen beginnen, sind die beiden größten Zeiger fast antiparallel, ihre Summe also fast 0. Deshalb werden weitere Zeiger benötigt, um die richtige Stärke des Summensignals auch am Rand des Messintervalls zu erhalten. Das sieht man noch besser, wenn man die Eigenfrequenz erhöht, z.B. auf f_0 = \frac{19}{4}\,\text{Hz} wie in Abb. 8.

AnimCutCos_19_4__243
Abb. 8: Eine cos-Funktion mit Periodendauer 4/19 s, die alle 2 s periodisch wiederholt wird. Gezeigt ist ein Zeitpunkt kurz vor dem Sprung.

Das entsprechende Spektrum zeigt Abb. 9. Außer dass f_0 jetzt größer ist, gibt es keinen Unterschied zum Spektrum aus Abb. 7.

SpektCutCos_19_4
Abb. 9: Das Amplituden- und Phasenspektrum der Funktion aus Abb. 8. Die »Eigenfrequenz« des Cosinus ist f_0 = \frac{19}{4}\,\text{Hz} = 4.75\,\text{Hz}.

Wären alle Amplituden unterhalb der beiden größten gleich 0, ergäbe sich zum selben Zeitpunkt wie in Abb. 8 die Situation in Abb. 10. In der Mitte des Messintervalls sind die hoch- und niederfrequenten Anteile eingerollt und spielen keine Rolle – am Rand des Messintervalls aber schon. Wie Abb. 10 auch zeigt, schaffen die hochfrequenten Anteile den Sprung ganz alleine.

AnimCutCos_19_4_0__243
Abb. 10: Wie Abb. 8, nur wurden alle Amplituden unterhalb der beiden größten auf 0 gesetzt.

Je kleiner die Sprünge am Rand des Messintervalls sind – je näher also die Messdauer einem ganzzahligen Vielfachen von 1/f_0 ist – desto kleiner sind die Amplituden der Frequenzen weit weg von der Eigenfrequenz.

Ein Sprung ist allerdings nicht das einzige Problem. In Abb. 11 ist der Cosinus einfach durch einen Sinus ersetzt worden. Dadurch gibt es am Rand der Messdauer keinen Sprung, aber eine Spitze. Weil die Eigenfrequenz des Sinus kein ganzzahliges Vielfaches von 1/T_m ist, müssen auch hier zwei relativ große Zeiger zusammenarbeiten, um näherungsweise einen Kreis vom Radius 2 zu erzeugen. Für die Spitze brauchen wir aber wieder viele kleine Zeiger. Im Prinzip haben wir das schon bei dem gleichgerichteten Sinus in Teil 2 gesehen.

AnimCutSin
Abb. 11: Eine sin-Funktion mit Periodendauer 4/7 s, die alle 2 s periodisch wiederholt wird, wodurch es jetzt zu einer Spitze kommt. Die Animation ist um den Faktor 4 verlangsamt, um den Knick besser sehen zu können.

Das Spektrum zeigt Abb. 12. Beim abgeschnittenen Cosinus haben sich die Zeiger für den Sprung entlang der reellen Achse addiert. Weil die Amplituden dort praktisch mit 1/f abnehmen, ist der Realteil im Wesentlichen die harmonische Reihe \sum_{k=1}^\infty\frac{1}{k}, die nicht konvergiert. Da uns aber nur der Imaginärteil interessiert, war das egal. Für den Knick in Abb. 11 müssen sich die Zeiger der hohen Frequenzen aber entlang der imaginären Achse ausrollen. Weil die Amplituden hier praktisch mit 1/f^2 abnehmen, ist der Imaginärteil im Wesentlichen die Reihe \sum_{k=1}^\infty\frac{1}{k^2}, die zu einem endlichen Wert konvergiert.

SpektCutSin
Abb. 12: Das Amplituden- und Phasenspektrum der Funktion aus Abb. 11. Die »Eigenfrequenz« des Sinus ist f_0 = \frac{7}{4}\,\text{Hz} = 1.75\,\text{Hz}.

Die Abb. 6 und 11 sind quasi die beiden Extremfälle. Meistens liegt man dazwischen, und hat einfach einen kleineren Sprung als von Minus die Amplitude auf Plus die Amplitude. Im Wesentlichen treten dadurch in den Zählern der strichlierten Einhüllenden in den Abb. 7, 9 und 12 Mischungen von f und f_0 auf, was qualitativ keine großen Auswirkungen hat.

Was wenn, wir kein reines Cosinus-Signal haben, sondern die Überlagerung der Signale aus den Abb. 6 und 8? Die komplexen Amplituden der beiden Signale mit denselben Frequenzen addieren sich, aber als Pfeile! Wir müssen daher auch ihre Phasen (Richtungen) beachten, nicht einfach nur ihre Beträge. Das Spektrum des Summensignals zeigt Abb. 13. Die strichlierte Kurve ist die Summe der Einhüllenden der beiden Einzelsignale. Offensichtlich addieren sich die Beträge der Amplituden der Einzelsignale nicht so einfach. Für sehr große bzw. sehr kleine Frequenzen passt es, aber nicht zwischen den Eigenfrequenzen der beiden Einzelsignale. (Würden die beiden Eigenfrequenzen weiter auseinander liegen, würde es auch dazwischen halbwegs passen.)

SpektCutCos_7_4_19_4
Abb. 13: Amplituden- und Phasenspektrum der Summe Funktionen aus den Abb. 6 und 8. Die strichlierte Kurve ist die Summe der Einhüllenden aus den Abb. 7 und 9.

Fensterfunktion (window function)

Die Signalsprünge/-spitzen an den Rändern des Messintervalls führen also zu unerwünschten Verbreiterungen im Spektrum. Es gibt eine radikale Methode, um dieses Problem zu beseitigen: Man multipliziert das Signal mit einer Funktion, die in der Mitte der Messdauer gleich 1 ist und an den Rändern langsam gegen 0 geht. Eine solche »Fensterfunktion« wäre zB das Hann-Fenster

\displaystyle \frac{1}{2}\left(1 - \cos\left(\frac{2\pi}{T_m} \cdot t\right)\right) .

Multiplizieren wir unseren abgeschnittenen Cosinus aus Abb. 6 mit dieser Fensterfunktion, ergibt sich die Funktion in Abb. 14. Zugegeben, das Signal sieht jetzt anders aus, aber die Eigenfrequenz des Cosinus ist noch da. Und um die geht es in erster Linie. Weil das Hann-Fenster dieselbe Periodizität wie das Messintervall hat, kommen auch keine zusätzlichen Frequenzen ins Spektrum.

AnimCutCosWindow
Abb. 14: Der abgeschnittene Cosinus mit Periodendauer 4/7 s aus Abb. 6, multipliziert mit einem Hann-Fenster. Die Animation ist für einen besseren Vergleich ebenfalls um den Faktor 4 verlangsamt.

Das Spektrum dieser Funktion zeigt Abb. 15. Die Einhüllende, auf der alle Amplituden liegen, ist jetzt etwas komplizierter (strichlierte Kurve). Im Vergleich zur Einhüllenden ohne Fenster (punktierte Kurve) gehen die Amplituden auf beiden Seiten der Eigenfrequenz jetzt viel schneller gegen 0. Zusammen mit einer noch größeren Frequenzauflösung (längeren Messdauer) sorgt das dafür, dass wir die Frequenz genauer bestimmen können.

SpektCutCosWindow
Abb. 15: Das Amplituden- und Phasenspektrum der Funktion aus Abb. 14. Die »Eigenfrequenz« des Cosinus ist f_0 = \frac{7}{4}\,\text{Hz} = 1.75\,\text{Hz}. Die rot strichlierte Kurve ist die Einhüllende der Amplituden. Die dunkelroten Punkte zeigen die Einhüllende aus Abb. 7.

Die Amplituden selber haben mit der ursprünglichen Funktion allerdings nicht mehr viel zu tun. Für die meisten Anwendungen geht es aber um eine genaue Bestimmung der Frequenzen.

Es gibt zig andere Fensterfunktionen  – mit speziellen Vor- und Nachteilen –  auf die ich hier nicht eingehen möchte.

Phasenspektrum

Nachdem die tatsächlich im Signal vorhandenen Frequenzen nicht direkt in der DFT vorkommen, dürfen wir uns vom Phasenspektrum nicht viel erwarten. Beim abgeschnittenen Cosinus in Abb. 7 hätten wir bei 1.75 Hz eine Phase von \pi/2 = \tau/4, beim abgeschnittenen Sinus in Abb. 9 entsprechend die Phase 0. Stattdessen sieht man einen Phasen»sprung«: Unterhalb von 1.75 Hz sind alle Phasen gleich \tau/2 (bzw. \tau/4), darüber sind alle 0 (bzw. -\tau/4).

Ein weiteres Problem für die Phase ist, dass wir unsere Messung irgendwann gestartet haben können, und nicht dann, wenn für das Signal t = 0 war. Deshalb haben wir es mit einer Zeitverschiebung \Delta t gemessen, die ebenfalls einen Einfluss auf die Phase hat (\varphi = -\omega \cdot \Delta t).

Diskussion

Wenn wir ein Signal, über das wir fast nichts wissen, einfach eine Zeit lang messen, erhalten wir sein Frequenz- und Amplitudenspektrum nur näherungsweise. Solange es nur darum geht, die im Signal vorhandenen Frequenz zu bestimmen, können verschiedene Fensterfunktionen ihre Dienste leisten. Die exakten Amplituden und speziell die Phasen erhält man so nicht.

Die folgenden Punkte sollte man jedoch beachten:

  • Die gewünschte Frequenzauflösung \Delta f = 1/T_m legt die Mindestmessdauer T_m \geq 1/\Delta f fest.
  • Man kann die Messdauer T_m nicht einfach vergrößern ohne auch die Anzahl N der Messwerte zu erhöhen, weil die Samplingfrequenz f_s = N/T_m mindestens doppelt so groß sein muss, wie die höchste interessante Frequenz im Signal.

Abschließend ist aber zu sagen, auch wenn wir nicht das exakte Spektrum des Signals erhalten, so können wir die gemessenen Signalwerte (und nur die) doch exakt aus dem Spektrum zurückrechnen.