In der Uni wird ständig mit der Fouriertransformation hantiert, aber zumeist nur in „formaler“ Weise, d.h. man löst z.B. die Wellengleichung im k-Raum oder stellt die Dirac-Funktion als Fouriertransformierte dar. Im 2. oder 3. Semester macht man auch die Spielchen, dass man eine Gaußfunktion per Hand transformiert und wieder eine Gaußfunktion herauskommt.
In der Praxis hat man ein Signal diskretisiert, wobei der Zeitabstand zwischen zwei Samplepunkten einer sinnvollen physikalischen Einheit entspricht. Nun möchte man „einfach mal kurz“ das Zeitsignal in seine Frequenzen zerlegen. Und genau hier wird der durchschnittliche Physiker durch die Implementation der FFT zur Weißglut getrieben. „Wie zur Hölle muss ich denn nun meine Frequenzachse konstruieren um physikalisch sinnvolle Frequenzen zu erhalten????“
Der Fallstrick ist eigentlich nicht die Implementation (z.B. in Python oder Matlab) selbst, sondern dass man glaubt alles - wirklich alles - über die Fouriertransformation zu wissen und daher liest man die zugehörige FFT-Dokumentation nicht richtig.
Im folgenden möchte ich diese Fallstricke anhand der Implementation der FFT im numpy-Paket in Python erläutern.
Das wichtigste:
- Die FFT spuckt ihr Ergebnis in einer „unintuitiven“ Ordnung aus („Standard Packing“)
- Die Normierung ist unsymmetrisch
- Im Exponenten ist ein Faktor 2pi versteckt
import matplotlib.pyplot as plt
import numpy as np
from math import *
N = 200
deltaT = 0.01 #ms
timeseries = np.zeros(N)
tAxis = np.linspace(0,deltaT*N,N)
damping = 15.0
frequency1 = 10.0 # 10 kHz
frequency2 = 30.0 # 30 kHz
fNyquist =0.5/deltaT
for i in range(N):
t = tAxis[i]
timeseries[i] = exp(-t*t*damping/4.0)*(1.0-exp(-t*t*damping))*sin(2.0*pi*frequency1*t)+0.05*cos(2.0*pi*frequency2*t)
plt.xkcd() #rausnehmen, falls matplotlib version <1.3.0
plt.plot(tAxis,timeseries,label = "damped sine")
plt.xlabel("t in ms")
plt.ylabel("Amplitude")
plt.legend()
plt.savefig("damped_sine.png")
Der Code zur Erzeugung dieser Grafik ist dieser:
FFTtimeseries = np.fft.fft(timeseries)
fMax = 1.0/deltaT
fAxis = np.linspace(0,fMax,N)
plt.clf()
plt.plot(fAxis,np.abs(FFTtimeseries),'-b',label="FFT")
plt.hold(True)
plt.plot(frequency1, max(np.abs(FFTtimeseries)),'ro',label="my guess!" )
plt.plot(frequency2, max(np.abs(FFTtimeseries)),'r+',label="my guess!" )
plt.plot(fNyquist, max(np.abs(FFTtimeseries)),'go',label="Nyquist frequency" )
plt.xlabel("f in kHz")
plt.ylabel("FFT")
plt.xlim([min(fAxis),max(fAxis)])
plt.savefig("wrong_axis.png")
Nun schaut man sich das an und krazt sich am Kopf: Was „zur Hölle“ sind das für 2 Peaks hinter der Nyquist-Frequenz? Ok, ich ignoriere die einfach, weil sie ja höher als die Nyquistfrequenz sind. Das ist falsch. Der Grund für die beiden hinteren Peaks ist die FFT-Implementation. Stichwort „Standard Packing“. Die Funktion np.fft.fft gibt ein Array der Länge N wieder (hier also 200 Elemente). Das 0. Element entspricht der Frequenz Null, die Elemente 1, 2, ..., N/2-1 entsprechen den DFT-Koeffizienten zur Frequenz deltaF, 2*deltaF, ..., (N/2-1)*deltaF. Wobei deltaF=1/(N*deltaT) ist. Also durch diejenige harmonische Schwingung definiert ist, die genau einmal in das vorgegebene Zeitfenster N*deltaT (hier 2 ms) passt. Die Elemente N/2,...,N-1 entsprechen dagegen den negativen Frequenzkomponenten -deltaF*(N-1), -deltaF*(N-2),... -deltaF. Um ein „erwartetes Ergebnis“ zu erhalten, muss man das Ergebnis der FFT umsortieren. Dazu ist schon die Funktion fft.fftshift implementiert.
Der Code:
FFTtimeseries = np.fft.fftshift(np.fft.fft(timeseries))
fMax = 1.0/deltaT
deltaF =1.0/(deltaT*N)
truefAxis = np.linspace(-fMax/2,fMax/2,N)
plt.clf()
plt.plot(truefAxis,np.abs(FFTtimeseries),'-b',label="FFT")
plt.hold(True)
plt.plot(frequency1, max(np.abs(FFTtimeseries)),'ro',label="my guess!" )
plt.plot(frequency2, max(np.abs(FFTtimeseries)),'r+',label="my guess!" )
plt.plot(fNyquist,max(np.abs(FFTtimeseries)) ,'go',label="Nyquist frequency" )
plt.xlabel("f in kHz")
plt.ylabel("FFT")
plt.xlim([min(truefAxis),max(truefAxis)])
plt.savefig("right_axis.png") Das Ergebnis:
Nochmal kurz die Rücktransformation austesten:
plt.clf()
plt.plot(tAxis,np.fft.ifft(np.fft.fft(timeseries)).real,'-r',label="back transformed")
plt.xlabel("t in ms")
plt.ylabel("Amplitude")
#plt.legend()
plt.savefig("back_transformed.png")
Das Ergebnis
Der zweite Punkt der Fallstrickliste soll anhand des Faltungstheorems diskutiert werden. Zur Erinnerung: Die Fouriertransformierte eines Faltungsintegrals ist einfach die Multiplikation der Fouriertransformierten der beiden einzelnen Funktionen. Außerdem ist ja bekannt, das die Faltung der Rechtecksfunktion mit sich selbst (Autokorrekation) die Zeltfunktion/Dreiecksfunktion gibt. Also versuchen wir uns mal daran:
def rect(x, L):
if np.abs(x) < L:
return 1.0
else:
return 0.0
tAxis = np.linspace(-5,5,N)
rectW = np.zeros(N)
for i in range(N):
rectW[i] = rect(tAxis[i],1.0)
deltaT =abs(tAxis[1]-tAxis[0])
plt.clf()
plt.plot(tAxis,rectW,'-b',label="rect window")
plt.xlabel("t in ms")
plt.ylabel("Amplitude")
plt.legend()
plt.ylim([0,1.25])
plt.savefig("rect.png")
iAmClever = np.fft.ifft(np.fft.fft(rectW)*np.fft.fft(rectW))
plt.clf()
plt.plot(tAxis,iAmClever,'-b',label="first try" )
plt.xlabel("t in ms")
plt.ylabel("Amplitude")
plt.legend()
plt.savefig("first_fft_conv.png")
Im obigen Code habe ich eine Funktion „rect“ definiert die uns ein Rechtecksfenster erzeugt und danach haben wir das Faltungstheorem angewendet. Die Ergebnisse sind...
FFTConvTimeseries = np.fft.ifftshift(np.fft.ifft(np.fft.fft(rectW)*np.fft.fft(rectW))).real*deltaT
plt.clf()
plt.plot(tAxis,FFTConvTimeseries,'-r',label="right axis and normalization" )
plt.xlabel("t in ms")
plt.ylabel("Amplitude")
plt.ylim([0,1.5])
plt.legend()
plt.savefig("fft_conv.png")
und man erhält das erwartete Ergebnis
