8.8 SPEKTRALANALYSE

 

 

EINFÜHRUNG

 

Fällt Licht auf dein Auge oder hörst du einen Ton, so entsteht ein Signal, das man als eine Funktion der Zeit y(t) auffassen kann. Für eine einzige Spektralfarbe oder einen reinen Ton handelt es sich um eine sinusförmige Schwingung mit der Amplitude A und Frequenz f, mathematisch ausgedrückt [mehr... Beginnt man mit der Zeitzählung nicht am Anfang der Sinusschwingung, so muss man noch eine Phasenverschiebung berücksichtigen.]:

y(t) = Asin(ω * t) mit ω = 2 * π * f

Ein komplexeres Signal, beispielsweise von einem konstant ausgehaltenen Ton eines Musikinstruments, ist noch immer periodisch, aber nicht mehr sinusförmig. Der berühmte Mathematiker Joseph Fourier (1768-1830) hat aber bewiesen, dass man jede periodische Funktion auch als eine Summe von Sinusschwingungen, als sogenannte Fourierreihe, auffassen kann. Er hat damit einen Grundstein von unschätzbarem Wert für den Fortschritt in der modernen Mathematik, Physik und Technik gelegt. Die Aufspaltung eines Signals in seine sinusförmigen Frequenzanteile nennt man Spektralanalyse.  

PROGRAMMIERKONZEPTE:
Sinusschwingung, Fourierreihe, Fast Fourier Transform (FFT), Spektrum, Sonogramm

 

 

SPEKTRUM EINES KLANGS, OBERTÖNE

 

Für den typischen Klangcharakter einer Stimme oder eines Musikinstruments sind die sinusförmigen Frequenzanteile wichtig, die darin enthalten sind. Für einen exakt periodischen Klang bestehen diese aus dem Grundton und den Obertönen, deren Frequenz ein ganzzahliges Vielfaches der Frequenz des Grundtons ist.  Trägst du die Amplituden der Frequenzanteile in einer Grafik auf, so erhältst du das Spektrum des Klangs. Ein Gerät, mit dem du das Spektrum bestimmen kannst, nennt man einen Spektralanalysator. TigerJython kann das Spektrum auf Grund eines Algorithmus bestimmen, der sehr berühmt ist und Fast Fourier Transform (FFT) heisst.

Um die FFT durchzuführen, übergibst du der Funktion fft(samples, n) eine Liste samples mit den zeitlich äquidistanten Abtastwerten und der Angabe n, wieviele davon vom Anfang der Liste an für die FFT zu verwenden sind.
 
Als Rückgabewerte erhältst du eine Liste mit den Amplituden der n/2 (normalisierten) Frequenzanteile zurück. Diese liegen im Abstand  r = fs / n, wobei fs die Abtastrate (sampling frequency) ist. r nennt man die Auflösung des Spektrums.

Diese n/2 Rückgabewerte im Abstand r "bevölkern" das Frequenzgebiet von 0 bis n/2*r = fs/2, oder kurz gesagt: Bei einer Abtastrate von fs liefert die FFT das Spektrum von 0 bis fs/2. Für eine Musik-CD mit der typischen Abtastrate fs = 44100 Hz entspricht dies einem Spektrum bis 22050 Hz, das den ganzen vom Menschen hörbaren Bereich umfasst.

Um den Spektrumanalysator zu testen, verwendest du vorerst einen Klang "wav/doublesine.wav" aus der Distribution von TigerJython, der zwei Sinustöne überlagert. Der Klang wurde mit einer Abtastrate von  fs = 40'000 Hz aufgenommen. Nimmst du  n = 10'000 Abtastwerte, so gibt dir die Funktion fft(samples, n)  5'000 Frequenzanteile mit der Auflösung r = 40'000 / 10'000 = 4 Hz im Bereich 0..20'000 Hz zurück, die du in einem GPanel als Spektrum mit vertikalen Linien grafisch darstellst.
 

from soundsystem import *
from gpanel import *

def showSpectrum(text):
    makeGPanel(-2000, 22000, -0.2, 1.2)
    drawGrid(0, 20000, 0, 1.0, 10, 5, "blue")
    title(text)
    lineWidth(2)
    r = fs / n # Resolution
    f = 0
    for i in range(n // 2): 
        line(f, 0, f, a[i])
        f += r

fs = 40000 # Sampling frequency
n = 10000 # Number of samples
samples = getWavMono("wav/doublesine.wav")
openMonoPlayer(samples, fs)
play()
a = fft(samples, n)
showSpectrum("Audio Spectrum")
Programmcode markieren (Ctrl+C kopieren, Ctrl+V einfügen)

 

 

MEMO

 

Wie du vermutest (und hörst),  handelt es sich um die Frequenzen 500 Hz und 1.5 kHz mit einem Amplitudenverhältnis von 1 : 1/2. Zusätzlich gibt es noch einige Störungsanteile. Die Frequenz 0 entspricht einem konstanten Signalanteil (offset).

Du hast jetzt mit TigerJython einen feudalen Spektrumanalysator vor dir, mit dem du die Grund- und Obertöne von Musikinstrumenten und menschlichen oder tierischen Lauten untersuchen kannst. In der Distribution findest du bereits den Ton einer Flöte ("wav/flute.wav")  und einer Oboe ("wav/oboe.wav"), deren Klangcharakteristik sehr unterschiedlich sind.

 
 

 

 

SPEKTRUM FÜR SELBSTDEFINIERTE FUNKTIONEN

 

Gemäss dem Satz von Fourier ist jede periodische Funktion mit der Frequenz f als Überlagerung von Sinusfunktionen mit den Frequenzen f, 2*f, 3*f, usw. darstellbar (Fourierreihe).

Mit deinem Fourieranalysator kannst du die Amplituden dieser Frequenzkomponenten experimentell bestimmen. Hier betrachtest du eine Rechteckschwingung mit der Frequenz f = 1 kHz. Die eingebaute Funktion square(A, f, t) liefert dir während der ersten Hälfte der Periode den Wert A und in der zweiten-A.  

Du wählst eine Abtastrate von fs = 40 kHz und bestimmst die Soundsamples für eine Dauer von 3 s (120'000 Werte). Dann spielst du den Soundclip ab. Für das Spektrum verwendest du aber lediglich 10'000 Werte und stellst das Spektrum dar.


 
from soundsystem import *
from gpanel import *

def showSpectrum(text):
    makeGPanel(-2000, 22000, -0.2, 1.2)
    drawGrid(0, 20000, 0, 1.0, 10, 5, "blue")
    title(text)
    lineWidth(2)
    r = fs / n # Resolution
    f = 0
    for i in range(n // 2): 
        line(f, 0, f, a[i])
        f += r  

n = 10000
fs = 40000 # Sampling frequency
f = 1000 # Signal frequency

samples = [0] * 120000  # sampled data for 3 s
t = 0
dt = 1 / fs # sampling period
for i in range(120000):
   samples[i] = square(1000, f, t)
   t += dt

openMonoPlayer(samples, 40000)
play()
a = fft(samples, n)
showSpectrum("Spectrum Square Wave")
Programmcode markieren (Ctrl+C kopieren, Ctrl+V einfügen)

 

 

MEMO

 

Das Experiment betätigt die Theorie, wonach das Spektrum einer Rechteckfunktion aus den ungeraden Vielfachen der Grundfrequenz besteht und sich die Amplituden der Spektralanteile wie 1, 1/3, 1/5, 1/7, usw. verhalten. Du kannst aber nie  experimentell herausfinden, dass die spektralen Anteile theoretisch bis ins Unendliche reichen.

 

 

SONOGRAMME

 

Die FFT ist ein perfektes Tool, um den spektralen Verlauf eines sich zeitlich ändernden Lauts aufzuzeichnen, beispielsweise eines gesprochenen Worts. Das Signal ist in diesem Fall natürlich nicht mehr periodisch, man kann aber davon ausgehen, dass es wenigstens stückweise einigermassen periodisch ist. Darum wendet man die FFT für kurze Signalblöcke an, beispielsweise für eine Blocklänge von 100 ms und wiederholt dies alle 2.5 ms. Damit ergibt sich alle 2.5 ms ein neues Spektrum, das man in einem Sonogramm als kolorierte vertikale Linie darstellen kann.

In deinem Programm beginnst du am Anfang der Abtastwerte und analysierst eine Blocklänge von 2000 Werten. Den nächsten Block beginnst du 50 Samples später, usw. In Python kannst du mit einer Slice-Operation

samples[k * 50:] mit k = 0, 1, 2,...

Es entsteht dadurch ein Sonogramm, beispielsweise für das gesprochene Wort "harris", das sich als "wav/harris.wav" in der Distribution von TigerJython befindet.


 

from soundsystem import *
from gpanel import *

def toColor(z):
    w = int(450 +  300 * z)
    c = X11Color.wavelengthToColor(w)
    return c

def drawSonogram():
    makeGPanel(0, 190, 0, 1000)
    title("Sonogramm of 'Harris'")
    lineWidth(4)
    # Analyse blocks every 50 samples
    for k in range(191):
        a = fft(samples[k * 50:], n)
        for i in range(n // 2):
            setColor(toColor(a[i]))
            point(k, i)

fs = 20000 # Sampling freq->spectrum 0..10 kHz
n = 2000 #  Size of block for analyser

samples = getWavMono("wav/harris.wav")
openMonoPlayer(samples, fs)
play()
drawSonogram()
Programmcode markieren (Ctrl+C kopieren, Ctrl+V einfügen)

 

 

MEMO

 

Das Sonogramm zeigt vertikal Frequenzen im Bereich 0..10 kHz, horizontal den zeitlichen Verlauf von 0 bis 190 * 50 / 20000 = 0.475 s.

Für die Umsetzung von Zahlen in Farben ist es zweckmässig, die Funktion X11Color.wavelengthToColor() zu verwenden, die Wellenlängen des Farbspektrums im Bereich 380...780 nm in Farben umsetzt.

Deutlich sichtbar sind die hohen Spektralanteile beim Zischlaut "s", bei dem die Grundtöne gänzlich fehlen.

 

 

LICHTSPEKTREN

 

Auch Licht kann spektral zerlegt werden, um die darin enthaltenen Wellenlängenanteile zu bestimmen. Das sichtbare Spektrum liegt im Wellenlängenbereich von ungefähr 380 nm bis 780 nm

Als Spektrumanalysator für Licht kennst du sicher das Prisma, bei dem das Licht für verschiedenen Wellenlängen gemäss dem Brechungsgesetz verschieden stark abgelenkt (gebrochen) wird.

 

In deinem Programm simulierst du den Übergang eines weissen Lichtstrahls in Glas und zeigst in einer Vergrösserung den Lichtweg für verschiedene Farben.

from gpanel import *

# K5 glass
B = 1.5220
C = 4590  # nanometer^2
# Cauchy equation for refracting index
def n(wavelength):
    return B + C / (wavelength * wavelength)

makeGPanel(-1, 1, -1, 1)
title("Refracting at the K5 glass")
bgColor("black")
setColor("white")
line(-1, 0, 1, 0)


lineWidth(4)
line(-1, 1, 0, 0)
lineWidth(1)

sineAlpha = 0.707

for i in range(51):
     wavelength = 380 + 8 * i
     setColor(X11Color.wavelengthToColor(wavelength))
     sineBeta = sineAlpha / n(wavelength)
     x = (sineBeta - 0.45) * 100 - 0.5  # magnification
     line(0, 0, x, -1)
Programmcode markieren (Ctrl+C kopieren, Ctrl+V einfügen)

 

 

MEMO

 

Um eine schöne Grafik zu erhalten, spaltest du die Farben mehr auf, als dies in Wirklichkeit der Fall ist.

 

 

AUFGABEN

 

1.


Untersuche andere Musikinstrumente oder Stimmen bezüglich ihres Obertongehalts und versuche, deine Feststellungen auf Grund des typischen Klangcharakters zu verstehen.


2.


Neben der globalen Funktion square(A, f, t) stehen dir noch die Funktionen sine(a, f, t), triangle(A, f, t), sawtooth(A, f, t) zur Verfügung. Stelle eine Vermutung über das Spektrum einer Dreieckschwingung und einer Sägezahnschwingung auf. Mit sine()  kannst du auch Überlagerungen von Sinusschwingungen untersuchen.


3.


Untersuche in Sonogrammen die Unterschiede zwischen verschiedenen Sprecherinnen und Sprechern, die das gleiche Wort sprechen.