Frequenzanalyse mit fft

Wenn du dir nicht sicher bist, in welchem der anderen Foren du die Frage stellen sollst, dann bist du hier im Forum für allgemeine Fragen sicher richtig.
Antworten
pypet
User
Beiträge: 6
Registriert: Dienstag 7. Juli 2015, 11:27

Hallo,
der u.g. Code stellt aus einem Ausschnitt eines Audiosignals (z.B. 1/20 Sekunde) die entsprechenden Frequenzanteile grafisch dar.

Code: Alles auswählen

import numpy as np
import matplotlib.pyplot as plt

from scipy.io.wavfile import read
(fs, at) = read('440_880_1320_4_3_2.wav')
# at: Amplitude in der zeitlichen Domaene
np.int_(at) # Umwandlung von int16 in int

Y = np.fft.fft(at)
N = len(Y)/2+1

t = np.arange(at.size)/float(fs)
dt = t[1] - t[0]
af = 1.0/dt # scan frequency

X = np.linspace(0, af/2, N, endpoint=True)

hann = np.hanning(len(at))
Yhann = np.fft.fft(hann*at)

plt.figure(figsize=(7,3))
plt.subplot(121)
plt.plot(t,at)
plt.title('Time Domain Signal')
plt.ylim(np.min(at)*3, np.max(at)*3)
plt.xlabel('Time ($sec$)')
plt.ylabel('Amplitude ($Unit$)')

plt.subplot(122)
plt.plot(X, 2.0*np.abs(Yhann[:N])/N)
plt.title('Frequency Domain Signal')
plt.xlabel('Frequency ($Hz$)')
plt.ylabel('Amplitude ($Unit$)')

plt.show()
Der Code orientiert sich an dieser Seite:
Paul Balzer: Die FFT mit Python einfach erklärt

Als Anfänger kann ich nicht sagen, daß ich den Code vollständig verstanden habe, doch er funktioniert im Wesentlichen, wie die folgende Probe zeigt:

Es wurde ein wav-file erstellt, welches ein Mischsignal enthält mit den drei Frequenzen (Sinus) 440, 880 und 1320 Hz sowie den entsprechenden Amplituden 0.4, 0.3 und 0.2 (Audacity). Download hier.

Wie kann ich jetzt die Werte der Maxima (Frequenz und Amplitude) numerisch ermitteln? Wenn man einen entsprechenden Ausschnitt aus der Grafik anschaut, habe ich den Verdacht, daß die ermittelten Frequenzen nicht genau mit den Ausgangsmaterial übereinstimmt. Wie kann das verbessert werden?

Was bedeutet [:N] in folgendem Ausdruck?
plt.plot(X, 2.0*np.abs(Yhann[:N])/N)

Wie kann die Darstellung im Frequenzbereich auf max. 20000Hz beschränkt werden?

Viele Grüße
Pypet
pypet
User
Beiträge: 6
Registriert: Dienstag 7. Juli 2015, 11:27

Hallo,

meine gestellten Fragen konnte ich inzwischen selbst beantworten und ich habe meinen Code auf das wav-Beispile abgestimmt.

Die tolle Methode von MagBen, lokale Maxima in einem Histogramm zu ermitteln und auszudrucken wurde gleich angewandt. Ich wollte ja prüfen, ob das Ergebnis der fft stimmt, bzw. wie genau die Methode für das definierte Beispiel ist.

Mein Verdacht der Ungenauigkeit hat sich zunächst bestätigt:

Code: Alles auswählen

Frequenz      Amplitude
443.697478992 129927.32287
887.394957983  97429.6595396
1331.09243697  64935.3416517
Ich kontrollierte daraufhin den wav-File und siehe da, es fehlte der erste Messpunkt (0/0). (Die automatischen Suche von Audacity nach Nulldurchgängen sollte kontrolliert werden.)
Mit einem entsprechend neu erstellten wav-File entsprach das Ergebnis den Erwartungen:

Code: Alles auswählen

Frequenz  Amplitude
440.0     129934.693614
880.0      97451.1102525
1320.0     64967.4107177
Die Frequenzgenauigkeit beträgt im Beispiel 100% und auch die Amplitudengenauigkeit ist sehr hoch:

Die doppelte Amplitude des 1320 Hz-Tons sollte diejenige des 440 Hz-Tones ergeben.
Ergebnis: 129 934,821435 statt 129 934.693614

Verbesserungen am Code werden gerne entgegengenommen.

Viele Grüße
pypet

Code: Alles auswählen

import numpy as np
import matplotlib.pyplot as plt

from scipy.io.wavfile import read
(fs, at) = read('440_880_1320_4_3_2_.wav')
# at: Amplitude in der zeitlichen Domaene
np.int_(at) # Umwandlung von int16 in int

Y = np.fft.fft(at)
# N = len(Y)/2+1
N = 0.1 * len(Y)/2+1

t = np.arange(at.size)/float(fs)
# t = np.arange(at.size)/(float(fs) * 0.1)
dt = t[1] - t[0]
af = 0.1/dt # scan frequency

X = np.linspace(0, af/2, N, endpoint=True)

hann = np.hanning(len(at))
Yhann = np.fft.fft(hann*at)
Yh = 2.0*np.abs(Yhann[:N])/N

# Maxima suchen und ausdrucken
y2 = Yh[1:-1]
x2 = X[1:-1]
sw = 100 # Schwellwert
maxima = (y2 > sw) & (y2 > Yh[:-2]) & (y2 > Yh[2:])

for x_max, y_max in zip(x2[maxima], y2[maxima]):
    print x_max, y_max

plt.figure(figsize=(7,3))
plt.subplot(121)
plt.plot(t,at)
plt.title('Time Domain Signal')
plt.ylim(np.min(at)*3, np.max(at)*3)
plt.xlabel('Time ($sec$)')
plt.ylabel('Amplitude ($Unit$)')

plt.subplot(122)
plt.plot(X, 0.1 * Yh)
plt.title('Frequency Domain Signal')
plt.xlabel('Frequency ($Hz$)')
plt.ylabel('Amplitude ($Unit$)')

plt.show()
Antworten