5.4 Diskrete Signale

Die diskrete Fourier-Transformation überführt eine periodische doppelt unendliche Zahlenfolge in eine endliche Folge von Zahlen.

Liegt nur ein diskretes Signal enlicher Länge vor, so kann dieses stets zu einer periodischen Folge fortgesetzt werden. Die diskrete Fourier-Transformation ist also auch auf Signale endlicher Länge anwendbar.

5.4.1Definition¶

Aus der Fourier-Transformierten kann man stets das ursprüngliche Signal zurückgewinnen.

Die diskrete Fourier-Transformation kann als Matrix-Vektor-Produkt formuliert werden: Ist $F\in\bbC^{N\times N}$ die durch

definierte Matrix und ist $a:=(a_0,a_1,\ldots,a_{N-1})$ eine Periode der Folge $(a_k)_{k\in\bbZ}$, so gilt

und umgekehrt

5.4.2Zusammenhang mit periodischen Funktionen¶

Sei $(\hat{a}_0,\hat{a}_1,\ldots,\hat{a}_{N-1})\in\bbC^N$ die Fourier-Transformierte von $(a_k)_{k\in\bbZ}$ und sei $T>0$. Die durch

definierte Funktion $f:\bbR\to\bbC$ hat folgende Eigenschaften:

• $f$ ist periodisch mit Periode $T$.

• Für jedes $k\in\bbZ$ gilt

  $$f\bigl(k\,\tfrac{T}{N}\bigr)=a_k,$$

  (8)

  d.h. $f$ interpoliert die Folge $(a_k)_{k\in\bbZ}$ auf den Stützstellen $0,\frac{T}{N},2\,\frac{T}{N},\ldots,(N-1)\,\frac{T}{N}$.

• Die Fourier-Transformierte von $f$ ist durch

  $$(\cF\,f)_k=\begin{cases}\hat{a}_k,&\text{für }k=0,\ldots,N-1,\\0,&\text{sonst}\end{cases}$$

  (9)

  gegeben.

In diesem Sinne ist die diskrete Fourier-Transformation also eine diskretisierte Version der Fourier-Transformation periodischer Funktionen.

5.4.3Reellwertige Signale¶

Analog zur Fourier-Transformation für periodische Funktionen können für reellwertige $N$-periodische Signale $(a_k)_{k\in\bbZ}$ komplexe Zahlen in bei der Fourier-Transformation vermieden werden, wenn Cosinus- und Sinus-Anteile separat betrachtet werden.

Für reellwertige Signale gilt außerdem

d.h. die Information über das ursprüngliche Signal ist vollständig in der Hälfte der Fourier-Komponenten enthalten. Genauer: Die Fourier-Transformierte eines $N$-periodischen reellwertigen Signals ist durch $N$ reelle Zahlen beschrieben (IDV 520).

5.4.4Mehrdimensionale Fourier-Transformation¶

Führen hier nur die Fourier-Transformation für zweidimensionale Signale ein. Für höhere Dimensionen erfolgt die Berechnung völlig analog.

Sei $A\in\bbC^{m\times n}$ ein zweidimensionales Signal mit einträgen $a_{m,n}$ (welches man sich auch periodisch auf $\bbZ\times\bbZ$ fortgesetzt vorstellen kann). Dann ist die zweidimensionale Fourier-Transformierte gegeben durch

für $k=0,\ldots,M-1$ und $l=0,\ldots,N-1$.

Die inverse Transformation ist

5.4.5Anwendungen¶

Die diskrete Fourier-Transformation ist von grundlegender Bedeutung in vielen Bereichen von Naturwissenschaft und Technik. Die Anwendung reiche von einfachen Aufgaben der Signalverarbeitung (Filter) bis zur Kompression von Audio- und Videosignalen und Schwingungsanalysen aller Art.

5.4.6Beispiel: Audiosignale¶

Wollen ein Audiosignal (Klaviermusik) im Frequenzbereich untersuchen.

5.4.6.1Digitale Audiosignale¶

Digitale Audiosignale sind folgen reeller Zahlen, die gewissen äquidistanten Zeitpunkten eine Auslenkung (des Luftdruck bzw. der Lautsprechermembran) zuordnen. Dabei werden üblicherweise einige 10000 Auslenkungswerte pro Sekunde gespeichert. Diese Hohe Anzahl ist nötig um auch hochfrequente Schwingungen hinreichend gut wiedergeben zu können.

from IPython.display import Audio
import matplotlib.pyplot as plt
import numpy as np
import scipy

# WAV-Datei mit einem Track (mono!) laden
samples_per_sec, data = scipy.io.wavfile.read('audio.wav')

# Daten plotten
t = np.arange(0, len(data)) / samples_per_sec  # Zeitpunkte
fig, ax = plt.subplots(figsize=(10, 3))
ax.plot(t, data, lw=0.1)
ax.set_xlim(t.min(), t.max())
ax.set_xlabel('Zeit (Sekunden)')
ax.set_ylabel('Auslenkung (beliebige Einheit)')
ax.set_title('Signal im Zeitbereich')
plt.show()

# Ausschnitt plotten
t_min = 1.4
t_max = 1.6
mask = np.logical_and(t > t_min, t < t_max)
fig, ax = plt.subplots(figsize=(10, 3))
ax.plot(t[mask], data[mask], lw=1)
ax.set_xlim(t_min, t_max)
ax.set_xlabel('Zeit (Sekunden)')
ax.set_ylabel('Auslenkung (beliebige Einheit)')
ax.set_title('Ausschnitt des Signals im Zeitbereich')
plt.show()

# abspielen
display(Audio(data, rate=samples_per_sec, autoplay=False))

5.4.6.2Spektrum¶

Die diskrete Fourier-Transformation kann mit numpy.fft.fft (Komponenten) und numpy.fft.fftfreq (zugehörige Frequenzen) leicht berechnet werden.

# diskrete F-Trafo mittels schneller F-Trafo
spec = 1 / len(data) * np.fft.fft(data)

# Frequenzen zu den Fourier-Komponenten
freq = np.fft.fftfreq(len(data), 1 / samples_per_sec)

# Fourier-Komponenten plotten
fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(10, 8), sharex=True)
ax1.plot(spec.real, lw=0.2)
ax2.plot(spec.imag, lw=0.2)
ax3.plot(np.abs(spec), lw=0.2)
ax3.set_xlim(0, len(spec))
ax3.set_xlabel('Index der Fourier-Komponenten')
ax1.set_title('Realteil der Fourier-Komponenten')
ax2.set_title('Imaginärteil der Fourier-Komponenten')
ax3.set_title('Betrag der Fourier-Komponenten')
plt.show()

Für die Darstellung als Spektrum müssen die Kompmentenindizes durch die zugehörigen Frequenzen ersetzt werden. Da Audiosignale reelle Zahlenfolgen sind, ist die vollständige Information über das Signal bereits in er Hälfte der Fourier-Komponenten enthalten.

Beachte: Da bei der diskreten Fourier-Transformation die Signale als periodische doppelt unendliche Folgen interpretiert werden, ist es egal, ob die Indizes (und damit die Frequenzen) von 0 bis $N-1$ oder von $-\frac{N-1}{2}$ bis $\frac{N-1}{2}$ laufen. Letztere Variante liefert statt hohen Frequenzen negative Frequenzen, was jedoch aus Sicht der diskreten Fourier-Transformation völlig äquivalent ist.

# Frequenzen plotten
fig, ax = plt.subplots(figsize=(10, 3))
ax.plot(freq)
ax.set_xlabel('Index der Fourier-Komponenten')
ax.set_ylabel('Frequenz (Hertz)')
ax.set_xlim(0, len(freq))
ax.set_title('Frequenzen zu den Fourier-Komponenten')
plt.show()

# halbes Amplitudenspektrum plotten (logarithmische Frequenzachse)
freq_pos_mask = freq > 0
fig, ax = plt.subplots(figsize=(10, 3))
ax.semilogx(freq[freq_pos_mask], np.abs(spec[freq_pos_mask]), lw=0.2)
ax.set_xlim(freq[freq_pos_mask].min(), freq[freq_pos_mask].max())
ax.set_xlabel('Frequenz (Hertz, logarithmisch)')
ax.set_ylabel('Auslenkung (beliebige Einheit)')
ax.set_title('Amplitudenspektrum (nur Hälfte der Frequenzen)')
plt.show()

5.4.6.3Interpretation¶

Im Spektrum kann man sehr deutlich Peaks bei den in der Musik üblichen Frequenzen erkennen (vgl. Frequenzen der gleichstufigen Stimmung).

# Frequenzbereich eines Klaviers (ohne Obertöne)
piano_mask = np.logical_and(freq > 16, freq < 4187)

# Spektrum mit Markierung der zu hörenden Noten
fig, ax = plt.subplots(figsize=(10, 3))
ax.semilogx(freq[piano_mask], np.abs(spec[piano_mask]), lw=1)
ax.set_xlim(freq[piano_mask].min(), freq[piano_mask].max())
ax.set_xticks([16.3516, 110, 130.813, 164.814, 195.998, 220, 261.626, 4186.01])
ax.set_xticklabels(['C0', 'A2', 'C3', 'E3', 'G3', 'A3', 'C4', 'C8'])
ax.get_xaxis().set_tick_params(which='minor', size=0)
ax.grid()
ax.set_xlabel('Frequenz (Klaviertaste, logarithmisch)')
ax.set_ylabel('Auslenkung (beliebige Einheit)')
ax.set_title('Spektrum und Noten')
plt.show()

Die nicht markierten Peaks entstehen durch die Obertöne, die zu jedem Grundton erklingen.

Der Peak bei E3 ist höher als bei anderen Noten, da diese Note im Audiosignal öfter vorkommt als die anderen.

5.4.6.4Kurzzeit-Fourier-Transformation¶

Die Kurzzeit-Fourier-Transformation (kurz: STFT für “short-time Fourier transform”; auch gefensterte Fourier-Transformation genannt) zerlegt das zu transformierende Signal in kleine, überlappende Teilstücke. Für jedes Teilstück wird dann eine diskrete Fourier-Transformation durchgeführt, sodass eine Folge zeitabhängiger Spektren entsteht. Das Zerlegen erfolgt dabei nicht durch Abschneiden, sondern durch Multiplikation mit einer glatten Funktion, z.B. einer Gauß’schen Glockenkurve. So werden zu starke Artefakte im Frequenzbereich durch die Sprünge an den Schnittstellen vermieden.

Spektrogramme lassen sich relativ leicht mit scipy.signal.ShortTimeFFT.spectrogram erzeugen.

# Fensterbreite und -versatz
win_width = 10000  # Breit in Samples
win_std = 5000  # Standardabweichung in Samples
hop = 200  # Versatz in Samples

# Spektrogramm erzeugen
win = scipy.signal.windows.gaussian(win_width, std=win_std, sym=True)
SFT = scipy.signal.ShortTimeFFT(win, hop=hop, fs=samples_per_sec)#, scale_to='magnitude')
s = SFT.spectrogram(data)

# Spektrogramm logarithmisch skalieren (für bessere Darstellung)
s_scaled =  np.log10(np.fmax(s, 5e7))

# Indizes in Sekunden und Hertz umrechnen
t_min, t_max, freq_min, freq_max = SFT.extent(len(data))
freq = np.linspace(freq_min, freq_max, s.shape[0])
freq_mask = freq < 500  # nur unteren Teil des Spektrogramms plotten

# plotten
fig, ax = plt.subplots(figsize=(10, 5))

# Spektrogramm
img = ax.imshow(
    s_scaled[freq_mask, :],
    origin='lower', aspect='auto',
    extent=(t_min, t_max, freq[freq_mask].min(), freq[freq_mask].max()),
    cmap='jet'
)
fig.colorbar(img, label='Amplitude (beliebige Einheiten, logarithmisch)')

# Linien für Klaviertasten
for f in [110, 130.813, 164.814, 195.998, 220, 261.626]:
    ax.plot([t_min, t_max], [f, f], '--k', lw=1)

ax.set_xlabel('Zeit (Sekunden)')
ax.set_ylabel('Frequenz (Hertz)')
plt.show()

Die farbkodierte Darstellung der Amplituden in Abhängigkeit von Zeit und Freqenz wird als Spektrogramm bezeichnet. Im Spektrogramm erkennt man neben den dominanten Frequenzen auch deren zeitliches Auftreten.