Mathematik für Ingenieure mit Python

Mathematische Grundlagen
www.grund-wissen.de/mathematik
WinPython
Alles was man für Mathematik mit Python unter Windows braucht.
winpython.sourceforge.net
Download WinPython
Python-Code
vollständiger Python-Code: fourier.py

Fourier-Analyse mit Numpy

Die Beispiel Funktion

import numpy as np

a     = 6
b     = 5.5
alpha = 0.1  
beta  = 1.5 

def fkt(t):
    "The example-function to be transformed"
    return a*np.exp(-alpha*t) - b*np.exp(-beta*t)

Die Fourier-Transformierte analytisch berechnet

def Fkt(f):
    "The analytical solution of the Fourier-transformation of the example-function."
    om = 2*np.pi*f
    return a/(alpha+1j*om)- b/(beta+1j*om)

Die Fourier-Transformierte aus der FFT berechnen

def fourier_transform(t, fkt):
    """
    Calculates the Fourier-Transformation of fkt with FFT. 
    The output Fkt of the Fourier-Transformation of fkt 
        is correctly normed (multiply with dt)
        is correctly ordered (swap first and second half of output-array).
    Calculates the frequencies from t.
    """
    n = t.size
    dt = t[1]-t[0]
    
    f=np.linspace(-1/(2*dt), 1/(2*dt), t.size, endpoint=False)
    Fkt = np.fft.fft(fkt)*dt
    
    # swap first and second half of output-array
    Fkt2 = np.fft.fftshift(Fkt)
    
    return f, Fkt2

Die Fourier-Koeffizienten der Fourier-Reihe berechnen

def fourier_coefficient(t, fkt, f):
    "Calculate only one fourier coefficient of fkt(t) for frequency f."
    phi = -2*np.pi*f*t
    return (fkt * (np.cos(phi)+1j*np.sin(phi))).sum()*(t[-1]-t[0])/t.size

Rechnung

n = 1e5
t_max = 300
t = np.linspace(0., t_max, n+1)[:-1]
s = fkt(t)
f, S1 = fourier_transform(t, s)

S2 = Fkt(f)

f3 = np.array([.1, .15, .2, -0.03])
S3 = np.array([fourier_coefficient(t,s,fi) for fi in f3])

Plotten

import pylab as plt
plt.figure()
plt.title("Time Domain")
plt.figtext(.6, .5, 
            r"$s(t) = a e^{-\alpha t}-b e^{-\beta t}$" "\n\n" +
            (r"$a = %G$" % a) + "\n" +
            (r"$b = %G$" % b) + "\n" +
            (r"$\alpha = %G$" % alpha) + "\n" +
            (r"$\beta = %G$" % beta),
            fontsize=20)
plt.plot(t,s)
plt.xlim(0,20)
plt.grid(True)
time domain function
plt.figure()
plt.title("Fourier-Analysis")
plt.figtext(.13, .25,
          r"$S(f)=\frac{a}{\alpha+j2\pi f}-\frac{b}{\beta + j 2\pi f}$",
          fontsize=24)
plt.plot(f,S2.real,              label="analytic real")
plt.plot(f,S2.imag,              label="analytic imag")
plt.plot(f,np.abs(S2),           label="analytic abs")
plt.plot(f,S1.real,   ".",       label="from FFT real")
plt.plot(f,S1.imag,   ".",       label="from FFT imag")
plt.plot(f,np.abs(S1), ".",      label="from FFT abs")
plt.plot(f3,S3.real,  "o", ms=7, label="coefficients real")
plt.plot(f3,S3.imag,  "o", ms=7, label="coefficients imag")
plt.plot(f3,np.abs(S3),"o", ms=7,label="coefficients abs")
plt.xlim(-.3,.3)
plt.grid(True)
plt.legend()

plt.show()
frequency domain Fourier transformed function