Implementierung eines Python-Notebooks, das für ein zu ladendes Signal eine diskrete Fourier-Analyse durchführt und anschließend mit Hilfe der Rücktransformation das Ausgangssignal rekonstruiert.
Im Projekt werden sowohl einfache, als auch komplexe Signale mit verschiedenen Implementierungen transformiert und im Anschluss rücktransformiert. Die Differenz zwischen Eingangs-und und Ausgangssignal wird visualisiert und kann durch verschiedene Tonbeispiele und Filtereinstellungen wahrgenommen werden.
Zudem wird der Code und die mathemtaische Grundlage der Fouriertransformation erklärt und demonstriert.
Die folgende Darstellung des Jupyter Notebooks ist lediglich demonstrativ und nicht funktional. Zum Testen kann das Jupyter Notebook heruntergeladen und gestartet werden.
- Repository clonen oder im Zip-Format herunterladen
- In den Repository-Ordner wechseln
- Requirements installieren:
pip install -r requirements.txt
- Im Ordner notebook den Server starten:
jupyter notebook
Schaut man sich den Schwingungsverlauf eines realen akkustischen Signals an, fällt relativ schnell auf, dass es sich hierbei nicht nur um eine einfache Sinus Schwingung handelt. Betrachtet man beispielsweise den folgenden Geigenton an der ersten rotmarkierten Stellen genauer, lässt sich eine komplexe Schwingung erkennen (vgl. Abbildung 1 Signalausschnitt beim Streichen einer Geigensaite).
Abbildung 1 Signalausschnitt beim Streichen einer Geigensaite
In der Abbildung unter dem Signalausschnitt ist das Spektrum, also die frequentielle Komposition des Signales abgebildet. Dabei lässt sich erkennen dass der Ton nicht nur aus einer einzigen Sinus/Cosinus-Schwingung besteht, sondern sich durch eine Überlagerung verschiedener Frequenzen mit unterschiedlich starken Amplituden zusammensetzt. Bezogen auf die Geige bildet sich der Ton, beim Zupfen der Saite, aus der Grundfrequenz
Um die frequentielle Zusammensetzung eines Schwingungssignals zu berechnen greift man auf die Fourier-Analysis zurück, um genau zu sein auf die Fouriertransformation.
from IPython.display import Audio
sound_file = 'files/Geigenton.wav'
Audio(sound_file)
Geigenton.mp4
Unter der Fourier-Analysis versteht man im Allgemeinen die Zerlegung und Synthese beliebig komplexer Schwingungen in einfachere Sinus/Cosinus-Schwingungen.
Der dabei wichtigeste Aspekt ist die sogenannte Fouriertransformation, wobei eine Schwingung, meistens aus dem Zeitbereich, in den Spektralbereich transformiert wird. Der Spektralbereich beschreibt die frequentielle Komposition der Schwingungsfunktion. Anhand dieser frequentiellen Signatur lässt sich das Ausgangssignal durch eine Fouriersynthese/inverse Fouriertransformation rekonstruieren
Bei Fouriertransformationen unterscheidet man zwischen der Zerlegung in eine Fourierreihe, der kontinuierlichen Fouriertransformation, der diskreten Fouriertransformation (DFT) und Fouriertransformationen für zeitdiskrete Signale (DTFT). Diese Typen unterscheiden sich hauptsächlich in der Kontinuität der Schwingung und der daraus resultierenden Spektralfunktion (vgl. Abbildung 1). Im Folgenden wird die diskrete Fouriertransformation, ihre Bedeutung, Anwendung und Theorie näher erläutert.
Abbildung 2 Fouriertransformation
Bei der diskreten Fouriertransformation liegt das Aussgangssignal als diskrete Funktion
Die diskrete Fouriertransformation ist folgendermaßen definiert:
Dabei fällt realtiv schnell auf, dass die realen Funktionswerte xj mit der komplexen Analysefunktion
Mit der eulerschen Formel
lässt sich das komplexe Ergebnis in einen Sinus und Cosinus Anteil zerlgen. Dies bedeutet, dass man für jede Frequenz k eine zugehörige Sinus und Cosinus Schwingung erhält. Der Betrag des komplexen Ausdrucks entspricht hierbei der Amplitude der Schwingungsanteile für die jeweillige Frequenz. Der Phasenwinkel des komplexen Ausdrucks entspricht der x-Achsen Verschiebung der Schwingungsanteile für die jeweillige Frequenz . Damit ist die Amplitude ein Maß dafür mit welchem Faktor die jeweillige Frequenz zur Gesamtschwingung beiträgt.
Da wir bei der diskreten Fouriertransformation nur mit diskreten Werten arbeiten können, kann nicht für jede beliebige Frequenz die dazugehörigen Schwingungsanteile errechnet werden.
Wir sind also auf die Anzahl der diskreten Werte von
Gehen wir davon aus, dass unser zu transformierendes Signal
Wollen wir also ein kontinuierliches Frequenzspektrum betrachten müsste man N gegen unendlich laufen lassen, was in der Praxis jedoch nicht möglich ist.
Im folgenden Codebeispiel wird die im oberen Abschnitt beschriebene Mathematik in ein Python-Programm übersetzt.
Als diskretes Ausgangssignal nehmen wir eine einfach Sinusschwingung mit der Frequenz
import numpy as np
import matplotlib.pyplot as plt
#Anzahl SamplePunkte
N=100
# MEssbereich
t_start=0
t_ende=2
#Samplingfrequenz
f_s=N/(t_ende-t_start)
# Erstellung der X-Achse
x = np.linspace(t_start,t_ende,N)
#Deklarieren von
y_sin = np.zeros(N)
f=2
for n in range(N):
y_sin[n]=np.sin((t_ende-t_start)*f*2*np.pi*n/N)
plt.plot(x,y_sin,'k.')
plt.ylabel('Amplitude in cm')
plt.xlabel('Zeit in s')
plt.show()
Im Anschluss wird die Funktion zur Berechnung der diskreten Fouriertransformation definiert. Dazu wird über alle Frequenzen k iteriert,auch wenn die zweite Hälfte der Frequenzen nach dem Nyquist-Abtasttheorems unwichtig ist. Für jede Frequenz k wird der dazu komlpexe Schwingungsterm c[k] berechnet.
def discrete_fourier_transformation(y):
N=y.size
c=np.zeros(N,complex)
n=np.arange(N)
#Diskrete Transformation mit For-Loops
for k in range(N):
c[k]=np.sum(y*np.exp(-2j*np.pi*k*n/N))
return c
#print(abs(c))
#Diskrete Transformation mit Matrix
#k= np.arange(N)[:,None]
#c = np.sum(y*np.exp(-2j*np.pi*k*n/N),axis=1)
Wendet man nun die diskrete Fouriertransfomation auf das oben erzeugte Sinussiganl an und schneidet die obere Hälfte aufgrund des Nyquist-Limits ab ergibt sich ein Schaubild, dass wie Erwartet zeigt, dass die Sinusfunktion mit der Frequenz
def plot_fourier_data(c,f_s):
N=c.size
plt.figure(figsize=(12,4))
plt.subplot(1,2,1)
# Es wird nur die Hälfte der Frequenzen geplottet weil die andere Hälfte nach Nyquist unzulässig ist
x = np.linspace(0,int(f_s/2),int(N/2)+1)
plt.plot(x,np.imag(c[:int(N/2)+1]), label="Imaginär")
plt.plot(x,np.real(c[:int(N/2)+1]),label="Real")
plt.ylabel('Amplitude in cm')
plt.xlabel('Frequenz in Hz')
plt.legend()
plt.subplot(1,2,2)
if N/2>50:
plt.plot(x,abs(c[:int(N/2)+1]),"r",label="Betrag")
else:
plt.plot(x,abs(c[:int(N/2)+1]),"k.",label="Betrag")
plt.ylabel('Amplitude in cm')
plt.xlabel('Frequenz in Hz')
plt.legend()
plt.show()
# Berechnung der Fouriertransformation
c_sin=discrete_fourier_transformation(y_sin)
#Ploten der Ergebnisse
plot_fourier_data(c_sin,f_s)
Im Folgenden wird aus dem Spektralbereich eine inverse Fouriertransformation durchgeführt. Dabei wird das Originalsignal aus der Fourierfunktion mit folgender Formel synthetisiert.
Durch eine Überlagerung der komplexen Schwingungsterme ergibt sich somit wieder das Originalsignal. Wenn ein konstantes Signal transformiert wird und im nachhinein synthetisiert, entsteht sogenanntes weißes Rauschen. Dieser Fall lässt sich in späteren Abschnitten erzeugen.
def discrete_fourier_synthese(c):
N=c.size
c_rück=np.zeros(N,complex)
k= np.arange(N)
for n in range(N):
c_rück[n]=1/N*np.sum(c*np.exp(2j*np.pi*k*n/N))
return c_rück
def plotcomparesignals(y,c_rück,t_start,t_ende):
N=y.size
x = np.linspace(t_start,t_ende,N)
plt.plot(x,np.real(y), label="Originales Signal")
if N>500:
plt.plot(x,np.real(c_rück), label="Synthetisiertes Signal")
else:
plt.plot(x,np.real(c_rück),"k.", label="Synthetisiertes Signal")
plt.ylabel('Amplitude a in cm')
plt.xlabel('Zeit t in s')
plt.legend()
plt.show()
c_rück_sin=discrete_fourier_synthese(c_sin)
plotcomparesignals(y_sin,c_rück_sin,t_start,t_ende)
Zur Berechnug des Fehler zwischen dem synthetisierten Signal und dem Originalsignal wird die Standardabweichung s berechnet. Es gilt:
def standardabweichung(c_rück,y):
c_rück=np.real(c_rück)
y, c_rück = np.array(y), np.array(c_rück)
return np.square(np.subtract(c_rück, y)).mean()
print("Die Standardabweichung des original Signals zum synthetisierten Signal beträgt: "+str(standardabweichung(c_rück_sin,y_sin)))
Die Standardabweichung des original Signals zum synthetisierten Signal beträgt: 1.176214079343573e-28
from __future__ import print_function
import time
import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets
def cosinus(N,f,intervall):
intervall_start,intervall_ende=intervall
#Samplingfrequenz
f_s=N/(t_ende-t_start)
# Erstellung der X-Achse
x = np.linspace(t_start,t_ende,N)
#Deklarieren von
y_cos = np.zeros(N)
for n in range(N):
y_cos[n]=np.cos((t_ende-t_start)*f*2*np.pi*n/N)
plt.plot(x,y_cos)
plt.ylabel('Amplitude in cm')
plt.xlabel('Zeit in s')
plt.show()
#DFT
c_cos=discrete_fourier_transformation(y_cos)
plot_fourier_data(c_cos,f_s)
#Intervall festlegen
#positive Frequenzen eliminieren
c_cos[0:intervall_start]=0
c_cos[intervall_ende+1:int(c_cos.size/2)]=0
#negative Frequenzen eliminieren
offset=int(c_cos.size/2)
c_cos[c_cos.size-intervall_start+1:c_cos.size]=0
c_cos[int(c_cos.size/2):c_cos.size-intervall_ende]=0
##########
c_rück_cos=discrete_fourier_synthese(c_cos)
plotcomparesignals(y_cos,c_rück_cos,t_start,t_ende)
print("Die Standardabweichung des original Signals zum synthetisierten Signal beträgt: "+str(standardabweichung(c_rück_cos,y_cos)))
# MEssbereich
t_start=0
t_ende=1
#Samplingfrequenz
Sample_Points_slider=widgets.IntSlider(
value=100,
min=2,
max=4000,
step=1,
description='Anzahl Samplepoints:',
disabled=False,
continuous_update=False,
orientation='horizontal',
readout=True,
readout_format='d'
)
Intervall_slider=intervall=widgets.IntRangeSlider(
value=[0,f_s],
min=0,
max=2056,
step=1,
description='Frequenzintervall:',
disabled=False,
continuous_update=False,
orientation='horizontal',
readout=True,
readout_format='d',
)
Frequenz_Slider=widgets.IntSlider(
value=1,
min=0,
max=50,
step=1,
description='Frequenz in Hz:',
disabled=False,
continuous_update=False,
orientsation='horizontal',
readout=True,
readout_format='d'
)
interact(cosinus,N=Sample_Points_slider,f=Frequenz_Slider,intervall=Intervall_slider)
# linking the two sliders to dynamically change the border of the slider on half of the Sample-Points
def on_value_change(change):
new_Sample_Points_slider = change['new']
Intervall_slider.max = int(new_Sample_Points_slider/2)
Sample_Points_slider.observe(on_value_change, names='value')
Die Standardabweichung des original Signals zum synthetisierten Signal beträgt: 5.455210080049842e-29
Wenn man die Frequenz des Signals auf 0 Hz stellt, fällt auf, dassbei der Rücktransformation ein spezieller Schwingungsgraph entsteht. Dieses Signal nennt man weißes Rauschen und entsteht bei einer IDFT auf einen konstanten Faktor.
from __future__ import print_function
import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets
def square(N,f,intervall):
intervall_start,intervall_ende=intervall
# Berechnung der Samplingfrequenz
f_s=N/(t_ende-t_start)
# Erstellung der X-Achse
x = np.linspace(t_start,t_ende,N)
#Deklarieren von
y_rechteck = np.zeros(N)
for n in range(N):
if np.sin((t_ende-t_start)*f*2*np.pi*n/N)>0:
y_rechteck[n]=1
else:
y_rechteck[n]=0
plt.plot(x,y_rechteck)
plt.ylabel('Amplitude in cm')
plt.xlabel('Zeit in s')
plt.show()
#DFT
c_rechteck=discrete_fourier_transformation(y_rechteck)
# Plotten der Ergebnisse
plot_fourier_data(c_rechteck,f_s)
#Intervall festlegen
#positive Frequenzen eliminieren
c_rechteck[0:intervall_start]=0
c_rechteck[intervall_ende+1:int(c_rechteck.size/2)]=0
#negative Frequenzen eliminieren
offset=int(c_rechteck.size/2)
c_rechteck[c_rechteck.size-intervall_start+1:c_rechteck.size]=0
c_rechteck[int(c_rechteck.size/2):c_rechteck.size-intervall_ende]=0
##########
#IDFT
c_rück_rechteck=discrete_fourier_synthese(c_rechteck)
# Plotten der Ergebnisse
plotcomparesignals(y_rechteck,c_rück_rechteck,t_start,t_ende)
# Fehler berechnen
print("Die Standardabweichung des original Signals zum synthetisierten Signal beträgt: "+str(standardabweichung(c_rück_rechteck,y_rechteck)))
# MEssbereich
t_start=0
t_ende=1
#Samplingfrequenz
Sample_Points_slider=widgets.IntSlider(
value=4000,
min=2,
max=4000,
step=1,
description='Anzahl Samplepoints:',
disabled=False,
continuous_update=False,
orientation='horizontal',
readout=True,
readout_format='d'
)
Intervall_slider=intervall=widgets.IntRangeSlider(
value=[0,8],
min=0,
max=2000,
step=1,
description='Frequenzintervall:',
disabled=False,
continuous_update=False,
orientation='horizontal',
readout=True,
readout_format='d',
)
Frequenz_Slider=widgets.IntSlider(
value=1,
min=0,
max=50,
step=1,
description='Frequenz in Hz:',
disabled=False,
continuous_update=False,
orientsation='horizontal',
readout=True,
readout_format='d'
)
interact(square,N=Sample_Points_slider,f=Frequenz_Slider,intervall=Intervall_slider)
plt.tight_layout()
# linking the two sliders to dynamically change the border of the slider
# linking the two sliders to dynamically change the border of the slider on half of the Sample-Points
def on_value_change(change):
new_Sample_Points_slider = change['new']
Intervall_slider.max = int(new_Sample_Points_slider/2)
Sample_Points_slider.observe(on_value_change, names='value')
Die Standardabweichung des original Signals zum synthetisierten Signal beträgt: 0.012600394253636267
Im folgenden Abschnitt werden wir die DFT auf ein komplexes, reales Tonsignal anwenden. Das vorliegende Testsignal hat folgende Daten:
- Auflösung pro Sample-Wert: 16bit
- Audiokanäle: mono
- Samplefrequenz: 16000 Hz
- Länge: 1s
In der Informatik hat sich eine Sample-Frequenz von 44100 Hz durchgesetzt. Wenn wir an dieser Stelle noch einmal das Nyquist-Abtasttheorem betrachten, fällt auch auf warum. Das menschliche Ohr kann Frequenzen zwischen 20 und 20 000 Hz erfassen. Ein mit 44100 Hz gesampeltes Signal lässt sich nach der oben beschriebenen Mathematik in 21500 Frequenzen zerlegen, was im Groben unserem Hörspektrum entspricht. Zur diskreten Fouriertransformation wird dazu der FFT (Fast Fourier Transform) Algorithmus benutzt. Der oben implementierte Algorithmus besitzt eine Komplexität von
from scipy.io import wavfile
import scipy.io
from IPython.display import Audio
import numpy as np
import matplotlib.pyplot as plt
#Abspielen der Audiodatei
sound_file = 'files/Geigenton_small_16000.wav'
display(Audio(sound_file))
Geigenton_small_16000.mp4
#Auslesen der Audiodatei
sound_file = 'files/Geigenton_small_16000.wav'
f_s, y_geige = wavfile.read(sound_file)
length=y_geige.size/f_s
x = np.linspace(0., length, y_geige.size)
plt.plot(x, y_geige)
plt.xlabel("Time t in s")
plt.ylabel("Amplitude")
plt.show()
#DFT
c_geige=discrete_fourier_transformation(y_geige)
plot_fourier_data(c_geige,f_s)
print(c_geige[1500])
#IDFT
c_rück_geige=discrete_fourier_synthese(c_geige)
plotcomparesignals(y_geige,c_rück_geige,0,length)
print("Die Standardabweichung des original Signals zum synthetisierten Signal beträgt: "+str(standardabweichung(c_rück_geige,y_geige)))
Im Vergleich dazu den FFT-Algorithmus auf das gleiche Signal. Die Ausführzeit ist um ein vielfaches geringer!
import scipy.fft
sound_file = 'files/Geigenton_small_16000.wav'
f_s, y_geige = wavfile.read(sound_file)
length=y_geige.size/f_s
x = np.linspace(0., length, y_geige.size)
plt.plot(x, y_geige)
plt.xlabel("Time t in s")
plt.ylabel("Amplitude")
plt.show()
c_geige=scipy.fft.fft(y_geige)
print(c_geige[1500])
plot_fourier_data(c_geige,f_s)
c_rück_geige=scipy.fft.ifft(c_geige)
plotcomparesignals(y_geige,c_rück_geige,0,length)
print("Die Standardabweichung des original Signals zum synthetisierten Signal beträgt: "+str(standardabweichung(c_rück_geige,y_geige)))
Die Standardabweichung des original Signals zum synthetisierten Signal beträgt: 2.0102026232882622e-16
Im Folgenden wird das komplexe Tonsignal von oben noch einmal genauer betrachtet. Wir betrachten einen kleineren Auschnitt und untersuchen welche Auswirkungen die Anzahl an Samplepoints auf die Genauigkeit des Signals hat.
Wir wissen, dass das Signal 16000 Sample-Points besitzt und betrachten im Folgenden nur
import numpy as np
import matplotlib.pyplot as plt
sound_file = 'files/Geigenton_small_16000.wav'
f_s, y_geige = wavfile.read(sound_file)
# Auslesen von 100 Sample-Points
y_geige=y_geige[0:100]
length=y_geige.size/f_s
x = np.linspace(0., length, y_geige.size)
plt.plot(x, y_geige)
plt.xlabel("Time t in s")
plt.ylabel("Amplitude")
plt.show()
c_geige=scipy.fft.fft(y_geige)
plot_fourier_data(c_geige,f_s)
c_rück_geige=scipy.fft.ifft(c_geige)
plotcomparesignals(y_geige,c_rück_geige,0,length)
print("Die Standardabweichung des original Signals zum synthetisierten Signal beträgt: "+str(standardabweichung(c_rück_geige,y_geige)))
Die Standardabweichung des original Signals zum synthetisierten Signal beträgt: 4.786527138398841e-26
Bisher haben wir alle IDFT mit allen vorliegenden Frequenzen durchgeführt. Im Folgenden wollen wir die IDFT auf einen gewissen Bereich beschränken. Dazu laden wir uns den gleichen Geigentestton wie im vorherigen Abschnitt, jedoch mit einer Samplefrequenz von
from __future__ import print_function
import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets
import scipy.fft
from scipy.io.wavfile import write, read
from IPython.display import Audio
def intervall_idft(intervall):
intervall_start,intervall_ende=intervall
x = np.linspace(0., length, y_geige.size)
plt.plot(x, y_geige)
plt.xlabel("Time t in s")
plt.ylabel("Amplitude")
plt.show()
# DFT umd Frequenzplot zu generieren
c_geige=scipy.fft.fft(y_geige)
plot_fourier_data(c_geige,f_s)
#Intervall festlegen
#positive Frequenzen eliminieren
c_geige[0:intervall_start]=0
c_geige[intervall_ende+1:int(c_geige.size/2)]=0
#negative Frequenzen eliminieren
offset=int(c_geige.size/2)
c_geige[c_geige.size-intervall_start+1:c_geige.size]=0
c_geige[int(c_geige.size/2):c_geige.size-intervall_ende]=0
##########
c_rück_geige=scipy.fft.ifft(c_geige)
plotcomparesignals(y_geige,c_rück_geige,0,length)
scaled = np.int16(np.real(c_rück_geige)/np.max(np.abs(np.real(c_rück_geige)))*32767)
write('files/own_intervall.wav', c_rück_geige.size, scaled)
print("Das synthetisierte Signale zum anhören:")
display(Audio('files/own_intervall.wav'))
print("Zum Vergleich das Originalsignal")
display(Audio('files/Geigenton_small_44100.wav'))
print("Die Standardabweichung des original Signals zum synthetisierten Signal beträgt: "+str(standardabweichung(c_rück_geige,y_geige)))
sound_file = 'files/Geigenton_small_44100.wav'
f_s, y_geige = read(sound_file)
f_s=44100
length=y_geige.size/f_s
interact(intervall_idft,intervall=widgets.IntRangeSlider(
value=[0, 44100],
min=0,
max=22050,
step=1,
description='Frequenzintervall:',
disabled=False,
continuous_update=False,
orientation='horizontal',
readout=True,
readout_format='d',
))
Das synthetisierte Signal mit eingeschänktem Frequenzintervall zum anhören:
Synthetisierter.Geigenton.mp4
Zum Vergleich das Originalsignal:
Originaler.Geigenton.mp4
Es fällt auf, dass wenn man das Samplefrequenzintervall zu klein macht, sich der Ton akkustisch wahrnehmbar verändert. Die obere Grenze des Intervalls lässt sich jedoch bis auf 8khz herunter drücken, ohne hörbare akkustische Änderungen wahrzunehmen. Dies liegt an der menschlichen Hörschwelle von ungefähr 16000 Hz (Erwachsener).
Ausführung in Visual Studio Code kann womöglich zu Fehlern führen: issue: microsoft/vscode-jupyter#1012