Fourier-Reihen, Teil 6 – Diskrete Fourier-Transformation (DFT)

In Teil 3 haben wir gesehen, dass wir ein periodisches Signal s mit Periodendauer T als Summe rotierender Zeiger

$\displaystyle s(t) = \sum_{-\infty}^{+\infty}\underline{S}_k \cdot e^{\underline{i}k\omega_1 t}$

schreiben können (zumindest wenn s »schön« ist). Dabei ist die Grundfrequenz $f_1 = 1/T$ und die Grundkreisfrequenz $\omega_1 = \tau/T$ mit $\tau = 2\pi$ .

Wir haben auch gesehen, dass wir die Fourier-Koeffizienten $\underline{S}_k$ über die Mittelwerte

$\displaystyle\underline{S}_k = \frac{1}{T} \int_0^T s(t) \cdot e^{-\underline{i}k\omega_1 t} \, \mathrm{d}t$

erhalten. Dabei müssen wir über eine ganze Periode integrieren, egal wo wir anfangen: 0 bis T, $-T/2$ bis $+T/2$ , $-T/4$ bis $+3T/4$ , …

Wenn wir den Verlauf des Signals s tatsächlich als mathematischen Funktionsterm kennen, sind diese Integrale prinzipiell berechenbar – auch wenn es manchmal kompliziert werden kann. Aber was, wenn wir den Funktionsterm des Signals nicht kennen, z.B. weil wir es gemessen haben? – In beiden Fällen können wir die Integrale zumindest näherungsweise numerisch berechnen.

Periodendauer des Signals bekannt

Zunächst wollen wir annehmen, dass wir zumindest die Periodendauer T des Signals tatsächlich kennen. Dann messen wir das Signal N-mal innerhalb der Periodendauer, sprich wir kennen seine Werte zu den Zeiten

$t_n = n \cdot \Delta t = n \cdot \dfrac{T}{N}$ für $n = 0, 1, 2, \dotsc, N - 1$

(s. Abb. 1). Man beachte, dass der Wert am Ende der Periode bei $t = T = 2\,\text{s}$ nicht dabei ist. Wenn N klein ist, kann das den Mittelwert beeinflussen.

Der Kehrwert des Messintervalls $\Delta t = T/N$ ist die Abtastrate bzw. Sampling-Frequenz

$f_s = \dfrac{1}{\Delta t} = \dfrac{N}{T}$ .

sampled_sig — Abb. 1: Ein periodisches Signal, von dem eine Periode hervorgehoben ist. Die roten Punkte sind die N gemessenen Werte in der Periodendauer T. Ihr zeitlicher Abstand ist $\Delta t = T/N$ . In diesem Beispiel ist N = 10.

Abb. 1: Ein periodisches Signal, von dem eine Periode hervorgehoben ist. Die roten Punkte sind die N gemessenen Werte in der Periodendauer T. Ihr zeitlicher Abstand ist $\Delta t = T/N$ . In diesem Beispiel ist N = 10.

Hintransformation

Mit den N gemessenen Signalwerten $s_n = s(t_n)$ versuchen wir jetzt, das Integral für die Fourier-Koeffizienten $\underline{S}_k$ näherungsweise zu berechnen, indem wir es durch eine Summe ersetzen:

$\displaystyle\begin{aligned}\underline{S}_k &\approx \frac{1}{T} \sum_{n = 0}^{N-1} s(t_n) \cdot e^{-\underline{i}k\omega_1 t_n} \, \Delta t\\ &= \frac{1}{T} \sum_{n = 0}^{N-1} s_n \cdot e^{-\underline{i}k \cdot \frac{\tau}{T} \cdot n \cdot \frac{T}{N}} \, \frac{T}{N}\\ &= \frac{1}{N} \sum_{n = 0}^{N-1} s_n \cdot e^{-\underline{i}\tau k n / N} \, ,\end{aligned}$

wobei wir $\omega_1 = \tau / T$ und $t_n = n \cdot \Delta t = n\cdot T / N$ benutzt haben. Je größer die Zahl der Messwerte $s_n$ ist, desto besser ist die Näherung; Integrale sind ja der Grenzwerte so einer Summe für $N \to \infty$ .

Für die Fourier-Reihe gab es unendlich viele Fourier-Koeffizienten – auch wenn fast alle gleich 0 sein konnten. Wie sieht es hier aus? Berechnen wir z.B.

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

Der Faktor $e^{-\underline{i}\tau n}$ ist aber gleich 1, weil

$\displaystyle e^{-\underline{i}\tau n} = \left(e^{-\underline{i}\tau}\right)^n = 1^n = 1$

ist; wir drehen unseren Zeiger ja nur n-mal im Kreis, was seine Richtung nicht ändert. Daher gilt

$\boxed{\underline{S}_{k+N} = \underline{S}_k}$ .

Anstatt unendlich vieler Fourier-Koeffizienten haben wir jetzt nur noch N verschiedene Koeffizienten für N Messwerte $s_n$ .

Rücktransformation

Wie können wir aus den $\underline{S}_k$ wieder unsere Messwerte $s_n$ erhalten? Wir können nicht so einfach in die Formel für $s(t)$ einsetzen, weil dort unendlich viele Koeffizienten benötigt werden, und wir haben nur noch N Stück.

Die N-te Einheitswurzel $e^{-\underline{i}\tau/N}$ wollen wir $\underline{\zeta}_N$ (»zeta-N«) nennen. Dann ist $e^{-\underline{i}\tau k n / N} = \left(e^{-\underline{i}\tau/N}\right)^{k n} = \underline{\zeta}_N^{k \cdot n}$ , und wir können unsere Formel für $\underline{S}_k$ schreiben als:

$\displaystyle N \cdot \underline{S}_k = s_0 \cdot \underline{\zeta}_N^{k \cdot 0} + s_1 \cdot \underline{\zeta}_N^{k \cdot 1} + \dotsb + s_n \cdot \underline{\zeta}_n^{k \cdot n} + \dotsb + s_{N - 1} \cdot \underline{\zeta}_n^{k \cdot (N - 1)}$ ,

weil wir ja von $n = 0$ bis $n = N - 1$ summieren. Von dieser Gleichung haben wir insgesamt N Stück, nämlich von $k = 0$ bis $k = N - 1$ . Multiplizieren wir alle diese Gleichungen mit $\underline{\zeta}_N^{-k \cdot n} = e^{\underline{i}\tau k n / N}$ , erhalten wir

$\displaystyle\begin{aligned}N \cdot \underline{S}_k \cdot e^{\underline{i}\tau k n / N} &= s_0 \cdot \underline{\zeta}_N^{k \cdot (0 - n)} + s_1 \cdot \underline{\zeta}_N^{k \cdot (1 - n)} + \dotsb \\ &\quad + s_n \cdot \underline{\zeta}_n^{k \cdot (n - n)} + \dotsb + s_{N - 1} \cdot \underline{\zeta}_n^{k \cdot (N - 1 - n)} \, .\end{aligned}$

Davon haben wir wieder N Stück. Und die addieren wir jetzt alle von $k = 0$ bis $k = N - 1$ :

$\displaystyle\begin{aligned}N \cdot \sum_{k=0}^{N - 1}\underline{S}_k \cdot e^{\underline{i}\tau k n / N} = {}&s_0 \cdot \left(\underline{\zeta}_N^{0 \cdot (0 - n)} + \underline{\zeta}_N^{1 \cdot (0 - n)} + \dotsb + \underline{\zeta}_N^{(N - 1) \cdot (0 - n)}\right) + {}\\ &s_1 \cdot \left(\underline{\zeta}_N^{0 \cdot (1 - n)} + \underline{\zeta}_N^{1 \cdot (1 - n)} + \dotsb + \underline{\zeta}_N^{(N - 1) \cdot (1 - n)}\right) + {}\\ &\dotsb +{}\\ &s_n \cdot \left(\underline{\zeta}_N^{0 \cdot (n - n)} + \underline{\zeta}_N^{1 \cdot (n - n)} + \dotsb + \underline{\zeta}_N^{(N - 1) \cdot (n - n)}\right) + {}\\ &\dotsb +{}\\ &s_{N - 1} \cdot \left(\underline{\zeta}_N^{0 \cdot ((N - 1) - n)} + \underline{\zeta}_N^{1 \cdot ((N - 1) - n)} + \dotsb + \underline{\zeta}_N^{(N - 1) \cdot ((N - 1) - n)}\right) \, .\end{aligned}$

Das sieht furchteinflößend aus, hat aber eine gewisse Struktur. Spannend ist, was jeweils in den großen runden Klammern passiert.

Dazu müssen wir uns wieder mit den Einheitswurzeln $\underline{\zeta}_N$ beschäftigen. Das sind komplexe Zahlen, die alle die Länge 1 haben; und wenn man sie zur N-ten Potenz nimmt, kommt genau 1 heraus. Für $N = 2$ wären es die beiden Einheitswurzeln $1$ und $-1$ , weil $1^2 = 1$ und $(-1)^2 = 1$ sind.

Schauen wir sie uns für $N = 5$ in Abb. 2(links) an. Die Einheitswurzeln sind Pfeile, die alle denselben Winkel zueinander haben, in diesem Fall $\tau/5 = 72^\circ$ . Wenn wir sie alle addieren, würde sich ein geschlossenes Polygon ergeben, und die Summe wäre gleich 0.

Abb. 2: Die 5-ten Einheitswurzeln $\underline{\zeta}_5^0$ bis $\underline{\zeta_5}^4$ (links). Multiplizieren wir alle Exponenten mit derselben ganzen Zahl, z.B. 2, erhalten wir dieselben Pfeile (rechts).

Abb. 2: Die 5-ten Einheitswurzeln $\underline{\zeta}_5^0$ bis $\underline{\zeta_5}^4$ (links). Multiplizieren wir alle Exponenten mit derselben ganzen Zahl, z.B. 2, erhalten wir dieselben Pfeile (rechts).

In den Klammern oben stehen aber leider nicht direkt die Einheitswurzeln, sondern sie werden alle mit derselben ganzen Zahl im Exponenten multipliziert. Wenn diese Zahl z.B. 2 ist, erhalten wir die Pfeile in Abb. 2(rechts). Das sind aber exakt dieselben Pfeile, nur anders nummeriert. D.h., auch so eine Summe muss gleich 0 sein.

Die vielen Klammern oben sind also alle gleich 0, mit einer Ausnahme: die Zeile, in der $s_n$ steht. Dort steht im Exponenten nämlich $(n - n) = 0$ , es werden also nur N Einser addiert. Letztlich reduziert sich dieser komplizierte Ausdruck also auf

$\displaystyle N \cdot \sum_{k=0}^{N - 1}\underline{S}_k \cdot e^{\underline{i}\tau k n / N} = s_n \cdot N$ ,

wo wir das N einfach kürzen können.

Dieses Argument ist ähnlich zu dem, das wir zur Berechnung der Fourier-Koeffizienten in Teil 3 verwendet haben. Es ist etwas komplizierter, weil die Zeiger hier jetzt nicht kontinuierlich rotieren, sondern wir nur »Schnappschüsse« zu den Zeiten $t_n = n \cdot T / N$ haben.

Diskrete Fourier-Transformation (DFT)

Für $n,k = 0, 1, 2, \dotsc, N - 1$ haben wir also das Transformations-Paar

$\displaystyle \boxed{s(n \cdot \Delta t) = s_n = \sum_{k = 0}^{N-1} \underline{S}_k \cdot e^{\underline{i}\tau k n / N}}$

und

$\displaystyle \boxed{\underline{S}(k \cdot \Delta f) = \underline{S}_k = \frac{1}{N} \sum_{n = 0}^{N-1} s_n \cdot e^{-\underline{i}\tau k n / N}}$

mit dem Zeitintervall $\Delta t = T/N = 1/f_s$ und dem Frequenzintervall $\Delta f = 1/T = f_1 = f_s / N$ .

Diese beiden Formeln gelten übrigens unabhängig davon, ob wir unsere Zahlen aus einer Messung oder sonst wie erhalten haben. Ich nehme N (komplexe) Zahlen, wende eine der beiden Formeln an und erhalte die anderen N Zahlen – und kann es auch wieder rückgängig machen.

Die komplexen Amplituden – die ich lieber mag – bekommen wir wieder mit

$\displaystyle\underline{A}_k = \begin{cases}\phantom{2}\underline{i} \cdot \underline{S}_0 & \text{wenn }k = 0\text{ ist,}\\ 2\underline{i} \cdot \underline{S}_k & \text{wenn }k > 0\text{ ist.}\end{cases}$

Abtasttheorem

Nehmen wir als Beispiel das Signal s mit

$s(t) = -1 + 3\sin\left(\frac{2\pi}{T} \, t + \pi\right) + 2\sin\left(5 \, \frac{2\pi}{T} \, t - \frac{\pi}{2}\right)$ ,

wie es in Abb. 1 gezeigt ist. Die Periodendauer soll $T = 2\,\text{s}$ sein. Die Grundfrequenz ist daher gleich $f_1 = 0.5\,\text{Hz}$ , der Oberton hat die Frequenz $f_5 = 5f_1 = 2.5\,\text{Hz}$ . Der Mittelwert beträgt -1 und die Amplituden von Grund- und Oberton sind 3 bzw. 2. Die entsprechenden Phasen sind $\pi$ bzw. $-\pi/2$ .

Die Berechnung der DFT überlassen wir dem Computer. MATLAB ist eine proprietäre Standardsoftware für solche Rechnungen. Der folgende Code wird aber auch vom freien Octave verstanden. Der fft-Befehl berechnet die Diskrete Fourier-Transformation mit dem Fast Fourier Transform-Algorithmus; dieser ist besonders schnell, wenn die Zahl der Messwerte (Samples) eine Zweierpotenz ist.

% Periodendauer (in s)
T = 2.0;

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

% Abtastzeiten (in s)
tn = (0:(N - 1))*T/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))/T;

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

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

plot(fk, abs(Ak), 'o')
xlabel('Frequenz in Hz')
ylabel('Amplitude')

Das Ergebnis ist das Amplitudenspektrum in Abb. 3. Es beginnt vielversprechend: Der Mittelwert hat Amplitude 1, der Grundton bei 0.5 Hz hat Amplitude 3 und der Oberton bei 2.5 Hz hat Amplitude 2, so wie es sein soll.

DFT_nur_A — Abb. 3: Das Amplitudenspektrum des Signals aus Abb. 1, berechnet mit einer DFT aus 2⁴ = 16 Werten.

Aber was sollen die Amplituden bei 5.5 Hz und 7.5 Hz? Wenn wir den Gleichanteil bei 0 Hz ignorieren, scheint das Spektrum symmetrisch um 4 Hz herum zu sein. Warum ist das so?

Es gilt:

$\displaystyle\begin{aligned}\underline{S}_{N - k} &= \frac{1}{N} \sum_{n=0}^{N - 1} s_n \cdot e^{-\underline{i}\tau(N - k) n / N}\\ &= \frac{1}{N} \sum_{n=0}^{N - 1} s_n \cdot \underbrace{e^{-\underline{i}\tau n}}_1 \cdot e^{+\underline{i}\tau k n / N}\\ &= \frac{1}{N} \sum_{n=0}^{N - 1} s_n \cdot e^{+\underline{i}\tau k n / N} \, .\end{aligned}$

Weil die $s_n$ reelle Zahlen sind, ändert sich für $\underline{S}_{N-k}$ verglichen mit $\underline{S}_k$ nur das Vorzeichen des Imaginärteils. Für reelle Signale s gilt daher

$\boxed{\underline{S}_{N - k} = \underline{S}_k^*}$ .

Die Beträge einer komplexen Zahl und ihrer konjugiert komplexen sind gleich, daher ist das Amplitudenspektrum symmetrisch.

D.h., ab $k = N - k$ erhalten wir keine neue Information mehr; also ab $k = \frac{N}{2}$ . Die Frequenz, bei der das passiert, ist die Nyquist-Frequenz

$f_\text{Nyquist} = \dfrac{N}{2} \cdot \Delta f = \dfrac{N}{2} \cdot \dfrac{1}{T} = \dfrac{1}{2} \cdot \dfrac{N}{T} = \dfrac{1}{2} \cdot f_s$ .

(Falls N keine gerade Zahl ist, müssten wir immer $\left\lfloor\frac{N}{2}\right\rfloor$ nehmen, die größte ganze Zahl, die nicht größer als $N/2$ ist.)

Im Beispiel oben ist die Sampling-Frequenz $f_s = 2^4/2\,\text{s} = 8\,\text{Hz}$ ; die Nyquist-Frequenz, um die herum das Spektrum symmetrisch ist, beträgt $f_\text{Nyquist} = f_s/2 = 4\,\text{Hz}$ . Eigentlich bräuchte man das Amplitudenspektrum nur bis zur Nyquist-Frequenz zu zeichnen.

Weil $f_s = 2 \cdot f_\text{Nyquist}$ ist, lautet das Abtasttheorem:

Die Sampling-Frequenz muss mindestens doppelt so groß sein wie die größte im Signal vorkommende Frequenz.

Das Problem ist, dass man das bei vielen Signalen nicht erfüllen kann. Wenn wir z.B. die Sägezahnschwingung aus Teil 2 nehmen, dann nehmen die Amplituden indirekt proportional mit der Frequenz ab, werden also nie exakt gleich 0. Es gibt daher oft keine größte Frequenz des Signals.

Aliasing

Was passiert, wenn unsere Abtastrate zu klein ist (Undersampling); wenn unser Signal also Frequenzen über $f_\text{Nyquist}$ beinhaltet? Lassen wir die Abtastrate wie oben bei $f_s = 8\,\text{Hz}$ , aber addieren zu unserem Signal einen Sinus mit Frequenz 7 Hz:

$s_2(t) = s(t) + 4\sin\left(14 \, \frac{2\pi}{T} \, t\right)$ .

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

Das entsprechende Amplitudenspektrum zeigt Abb. 4. Tatsächlich haben wir die richtige Amplitude bei 7 Hz, aber es taucht dieselbe Amplitude fälschlicherweise auch bei 1 Hz auf.

DFT_nur_A_Alias — Abb. 4: Amplitudenspektrum des Signals $s_2$ . Das 7 Hz-Signal sieht man fälschlicherweise auch bei 1 Hz.

Abb. 4: Amplitudenspektrum des Signals $s_2$ . Das 7 Hz-Signal sieht man fälschlicherweise auch bei 1 Hz.

Der Grund ist wieder, dass $\underline{S}_k^* = \underline{S}_{N - k}$ ist. Dadurch werden nicht nur die Amplituden kleiner Frequenzen »nach oben« gespiegelt, sondern auch umgekehrt. Diesen Effekt nennt man Aliasing.

Wie es anschaulich zu dieser Rückfaltung von Amplituden kommen kann, zeigt Abb. 5. Die grüne Kurve ist unser zusätzliches 7 Hz-Signal und die blaue Kurve ist das entsprechende 1 Hz-Signal mit einer Phasenverschiebung von $\pi$ (komplex konjugieren bedeutet, das Vorzeichen des Imaginärteils zu ändern, und $-\sin(x) = \sin(x + \pi)$ !).

Ursache_Aliasing — Abb. 5: Die Funktionen $4\sin\left(2\pi \cdot 7\,\text{Hz} \cdot t\right)$ (grün) und $4\sin\left(2\pi \cdot 1\,\text{Hz} \cdot t + \pi\right)$ (blau). Die roten Punkte sind $N = 2^4 = 16$ Messwerte in einer 2 s-Periode der Grundschwingung.

Abb. 5: Die Funktionen $4\sin\left(2\pi \cdot 7\,\text{Hz} \cdot t\right)$ (grün) und $4\sin\left(2\pi \cdot 1\,\text{Hz} \cdot t + \pi\right)$ (blau). Die roten Punkte sind $N = 2^4 = 16$ Messwerte in einer 2 s-Periode der Grundschwingung.

Die roten Punkte sind unsere $N = 2^4 = 16$ Messwerte. Wie man sieht, können wir aufgrund dieser Messwerte nicht zwischen den beiden Kurven unterscheiden. Hätten wir mehr Punkte gemessen, z.B. doppelt so viele, könnten wir das sehr wohl!

Sollte die tiefere Frequenz von Haus aus eine Amplitude größer 0 haben, müsste man diese Amplitude zu der nach unten gefalteten addieren (unter Berücksichtigung der relativen Phasen von $\underline{S}_k$ und $\underline{S}_{N-k}^*$ ).

Ergänzung 2018-08-14: Wenn Frequenzen jenseits der Sampling-Frequenz im Signal vorhanden sind, bekommt man sogar mehrere Kopien. Z.B., wenn wir zu unserem Ausgangssignal eine Anteil mit 11 Hz addieren:

$s_3(t) = s(t) + 4\sin\left(22 \, \frac{2\pi}{T} \, t\right)$ .

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

Der Fourier-Index dafür ist $k = 11\,\text{Hz}/0.5\,\text{Hz} = 22$ . Das entsprechende Amplitudenspektrum zeigt Abb. 4-b. Die 11 Hz kommen in dem Spektrum prinzipiell nicht vor. Aber: $22 = 6 + 16$ , d.h.

$\underline{S}_{22} = \underline{S}_{6 + 16} = \underline{S}_{6}$ und $\underline{S}_{6} = \underline{S}_{16 - 6}^* =\underline{S}_{10}^*$ ,

weil $\underline{S}_{k + N} = \underline{S}_k$ und $\underline{S}_{N - k} = \underline{S}_k^*$ sind. Daher sehen wir das 11 Hz-Signal fälschlicherweise bei $6 \cdot 0.5\,\text{Hz} = 3\,\text{Hz}$ und bei $10 \cdot 0.5\,\text{Hz} = 5\,\text{Hz}$ .

DFT_nur_A_Alias_22 — Abb. 4-b: Amplitudenspektrum des Signals $s_3$ . Das 11 Hz-Signal sieht man fälschlicherweise bei 3 Hz und 5 Hz.

Abb. 4-b: Amplitudenspektrum des Signals $s_3$ . Das 11 Hz-Signal sieht man fälschlicherweise bei 3 Hz und 5 Hz.

Schauen wir uns noch die Sägezahnfunktion aus Abb. 6 in Teil 2 an. Dazu müssen wir nur die Berechnung der Messwerte durch

% das Signal
sn = pi/2 - pi/2*tn; % Sägezahn

ersetzen.

Abb. 6 vergleicht die Amplitudenspektren für verschiedene Anzahlen von Messwerten. Für N = 16 gibt es durch das Aliasing einige Abweichungen vom erwarteten Verhalten. Der Gleichrichtwert ist nicht gleich 0, weil der letzte Punkt am Ende des Signals fehlt – der wäre bei $t = N \cdot \Delta t$ und das gehört nicht mehr zu unseren Messwerten. Für N = 1024 ist das Aliasing stark reduziert; wir bekommen fast die exakten Werte heraus. Im unteren Spektrum könnten wir natürlich weiter als bis 4 Hz gehen, nicht aber im oberen.

DFT_Saegezahn — Abb. 6: Das DFT-Amplitudenspektrum einer Sägezahnfunktion aus N = 16 Messwerten (oben) und N = 1024 Messwerten (unten). Die grün strichlierte Linie zeigt den erwarteten Verlauf $0.5\,\text{Hz} / f$ .

Abb. 6: Das DFT-Amplitudenspektrum einer Sägezahnfunktion aus N = 16 Messwerten (oben) und N = 1024 Messwerten (unten). Die grün strichlierte Linie zeigt den erwarteten Verlauf $0.5\,\text{Hz} / f$ .

Phasenspektrum

Gehen wir wieder zum Signal s aus Abb. 1 zurück, dessen Amplitudenspektrum in Abb. 2 gezeigt wurde. Um das Phasenspektrum zu berechnen müssen wir ab der Berechnung der komplexen Amplituden den folgenden Code einfügen:

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

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

% Phasen der komplexen Amplituden
phik = angle(Ak)/pi; % Phasen als Vielfache von pi

subplot(2,1,1);
plot(fk, absAk, 'o')
ylabel('Amplitude')

subplot(2,1,2);
plot(fk, phik, '+')
xlabel('Frequenz in Hz')
ylabel('Phase / pi')

Dieser Code zeichnet Amplituden- und Phasenspektrum übereinander. In Abb. 7 zeige ich jetzt nur das Phasenspektrum. Und das sieht schlimm aus.

DFT_nur_phi — Abb. 7: Das Phasenspektrum des Signals aus Abb. 1 bis zur Nyquist-Frequenz von 4 Hz.

Der Gleichrichtwert bei 0 Hz hat die Phase $-0.5\pi$ , weil der Mittelwert negativ ist. Passt! Die Grundschwingung bei 0.5 Hz hat die Phase $\pi$ und die Oberschwingung bei 2.5 Hz hat die Phase $-0.5\pi$ . Auch die beiden passen. Alle anderen Phasen sollten aber 0 sein! Stattdessen sehen wir irgendwas.

Warum ist das so? Obwohl Abb. 2 die Amplituden bei diesen Frequenzen als 0 zeigt, sind sie nicht wirklich exakt 0. Die komplexe Amplitude bei 1 Hz ist z.B.

$\underline{A}_2 \approx 3.5746 \times 10^{-16} - 7.0481 \times 10^{-16} \, \underline{i}$ ,

und ihre berechnete Phase ist $\arg(\underline{A}_2) \approx -0.3506\pi \approx -1.1014 \approx -63^\circ$ . Bei einer Amplitude, deren Betrag praktisch 0 ist, ist das allerdings nur ein numerisches Artefakt.

Die Lösung besteht darin, den größten Absolutwert der Amplituden zu suchen (außer dem Mittelwert, der kann irgendwas sein). Für jede Amplitude, deren Betrag viel kleiner als dieses Maximum ist, setzen wir die Phase manuell auf 0. Das erledigt der folgende Code:

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

% maximaler Betrag der komplexen Amplituden
% Mittelwert spielt nicht mit
maxAbsAk = max(absAk(2:end));

% Phasen der komplexen Amplituden
phik = angle(Ak)/pi; % Phasen als Vielfache von pi
phik(absAk < maxAbsAk/1000) = 0;

Für jede Amplitude, die 1000-mal kleiner als das Maximum ist, wird die Phase auf 0 gesetzt. Welchen genauen Faktor man nimmt, hängt von der Anwendung ab. Das korrigierte Phasenspektrum zeigt die Abb. 8.

DFT_nur_phi_korr — Abb. 8: Das korrigierte Phasenspektrum.

Periodendauer des Signals unbekannt

Messtechnisch ist natürlich der Fall interessanter, dass wir die Periodendauer des Signals nicht im Vorhinein kennen. Nachdem der Beitrag aber schon so lang geworden ist, verschiebe ich diese Diskussion auf den nächsten Teil 6b.

Kommentar verfassen Antwort abbrechen

Gib hier deinen Kommentar ein ...

Trage deine Daten unten ein oder klicke ein Icon um dich einzuloggen:

E-Mail (Adresse wird niemals veröffentlicht)

Name

Website

Du kommentierst mit Deinem WordPress.com-Konto. ( Abmelden / Ändern )

Du kommentierst mit Deinem Twitter-Konto. ( Abmelden / Ändern )

Du kommentierst mit Deinem Facebook-Konto. ( Abmelden / Ändern )

Abbrechen

Verbinde mit %s

Diese Seite verwendet Akismet, um Spam zu reduzieren. Erfahre, wie deine Kommentardaten verarbeitet werden..