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.

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.

zeta_5
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.

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 24 = 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.

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.

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.

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.

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

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.