Ing Mathematik: FFT und Wavelets

Fourierreihen

Jean Baptiste Joseph Fourier (französischer Mathematiker und Physiker, 1768 - 1830)

$f(x)\approx \frac{a_0}{2}+\sum _{k=1}^{\mathrm{\infty }}(a_k\cos \left(kx\right)+b_k\sin \left(kx\right))$

$a_k=\frac{1}{\pi }\underset{-\pi }{\overset{\pi }{\int }}f(x)\cdot \cos \left(kx\right)\mathrm{d}x\text{für }k\ge 0$

$b_k=\frac{1}{\pi }\underset{-\pi }{\overset{\pi }{\int }}f(x)\cdot \sin \left(kx\right)\mathrm{d}x\text{für }k\ge 1$

bzw. in komplexer Schreibweise

$f(x)\approx \sum _{k=-\mathrm{\infty }}^{\mathrm{\infty }}{c}_k{\mathrm{e}}^{jkx}$

$c_k=\frac{1}{2\pi }\underset{-\pi }{\overset{\pi }{\int }}f(x)\cdot {\mathrm{e}}^{-jkx}\mathrm{d}x$

Beispiel: Gegeben sei eine Rechteckfunktion. Sie soll mit der Fourierreihe angenähert werden.

• 
• 
• 

Dazugehöriger Python-Code (rudimentär):

from scipy import integrate
import numpy as np
import matplotlib.pyplot as plt

def ak(k):
    def y(x):
        if x <= -1 or x >= 2:
            return 0
        else:
            return 10 * np.cos(x*k)

    I, err = integrate.quad(y, -np.pi, np.pi)
    return 1./np.pi * I

def bk(k):
    def y(x):
        if x <= -1 or x >= 2:
            return 0
        else:
            return 10 * np.sin(x*k)

    I, err = integrate.quad(y, -np.pi, np.pi)
    return 1./np.pi * I

n = 2

a = np.zeros([n+1])
b = np.zeros([n+1])

for i in np.arange(0, n):
    a[i] = ak(i)
    b[i] = bk(i)

fx = a[0]/2.

x = np.arange(-1, 2.01, 0.01)

for i in np.arange(1, n):
    fx += a[i]*np.cos(i*x) + b[i]*np.sin(i*x)

plt.plot([-1,2], [10,10], color="green")
plt.plot([-1,-1], [0,10], color="green")
plt.plot([2,2], [0,10], color="green")
plt.plot(x, fx, label="n = " + str(n))
plt.legend()
plt.show()

Siehe auch Vorlage:W

Fourier-Transformation

$ℱ\{f(t)\}=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\mathrm{e}}^{-j\omega t}f(t)\mathrm{d}t$

${ℱ}^{-1}\{F(\omega )\}=\frac{1}{2\pi }\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\mathrm{e}}^{j\omega t}F(\omega )\mathrm{d}\omega$.

Siehe auch Vorlage:W

DFT - Diskrete Fourier Transformation

Diskrete Fourier-Transformation: ${\hat{a}}_k=\sum _{j=0}^{N-1}{\mathrm{e}}^{-2\pi \mathrm{i}\cdot \frac{jk}{N}}\cdot a_j$; $k=0,\dots ,N-1$

Matrix-Vektor-Form: $\hat{a}=\underset{\text{Fouriermatrix}}{\underbrace{W}}\cdot a$ mit $W[k,j]={\mathrm{e}}^{-2\pi \mathrm{i}\cdot \frac{jk}{N}}$

Inverse Diskrete Fourier-Transformation, iDFT: $a_k=\frac{1}{N}\sum _{j=0}^{N-1}{\mathrm{e}}^{2\pi \mathrm{i}\cdot \frac{jk}{N}}\cdot {\hat{a}}_j$ für $k=0,\dots ,N-1$

iDFT in Matrix-Vektor-Form: $a=W^{-1}\cdot \hat{a}$ mit $W^{-1}[k,j]=\frac{1}{N}{\mathrm{e}}^{2\pi \mathrm{i}\cdot \frac{jk}{N}}$

Siehe auch Vorlage:W

FFT - Fast Fourier Transformation

Eine effizientere Methode ist die FFT. Sie verwendet die Divide-and-Conquer-Methode. Zu diesem Verfahren gibt es vorgefertigte SciPy-Funktionen. Siehe z.B. [1]

Python-Code: FFT und Rücktransformation

import numpy as np
from scipy.fftpack import fft, ifft

y = np.array([1., 0.5, 2.5, 0.5, 3.5, 2.5]) 

u = fft(y)
v = ifft(u)

print("Messwerte = ", y)
print("FFT = ", u)
print("inverse FFT = ", v)

Ausgabe:

Messwerte =  [1.  0.5 2.5 0.5 3.5 2.5]
FFT =  [10.5-0.j         -1. +2.59807621j -3. +0.8660254j   3.5-0.j
 -3. -0.8660254j  -1. -2.59807621j]
inverse FFT =  [1. +0.j 0.5+0.j 2.5+0.j 0.5+0.j 3.5-0.j 2.5+0.j]

Python-Code: Rekonstruktion (Siehe auch Steinkamp: Der Python-Kurs für Ingenieure und Naturwissenschaftler; Rheinwerk, 2021, Seite 328ff oder Woyand: Python für Ingenieure und Naturwissenschaftler. 4. Aufl., Hanser, 2021, Seite 294ff)

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft, ifft, fftfreq

f = 50              # Frequenz
fgrenz = 1.5 * f    # Grenzfrequenz
n = 1000            # Anzahl Abtastungen
t = np.linspace(0, 1/f, n)
ut = np.sin(np.pi*f*t) + 3*np.random.randn(t.size)

u = fft(ut)
fk = fftfreq(ut.size, 1/(f*n))
Fg = u *(np.abs(fk) < fgrenz)
v = ifft(Fg)
z = np.abs(u[0:n])

# Grafikausgabe
fig, ax = plt.subplots(3,1)
ax[0].set_xlabel("Zeit")
ax[0].set_ylabel("u(t)")
ax[0].set_title("Messwerte")
ax[0].plot(t, ut, lw=1)
ax[1].set_xlabel("Zeit")
ax[1].set_ylabel("u(t)")
ax[1].set_title("Rekonstruktion")
ax[1].plot(t, v.real, lw=2)
ax[2].set_xlabel("Frequenz")
ax[2].set_ylabel("|z(t)|")
ax[2].set_title("Frequenzspektrum")
ax[2].plot(t, z, lw=1)
fig.tight_layout()
plt.show()

Ausgabe:

Siehe auch Vorlage:W

Short-time-Fourier-Transformation, STFT

Kontinuierliche STFT: $X(\tau ,\omega )=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}x(t)\underset{\text{Fensterfunktion}}{\underbrace{w(t-\tau )}}{\mathrm{e}}^{-\mathrm{j}\omega t}\mathrm{d}t$

Typische Fensterfunktionen: Vorlage:W, Vorlage:W, Vorlage:W.

Siehe auch Vorlage:W, Vorlage:W.

Z-Transformation

Z-Transformation: $𝒵\{f_n\}=\sum _{n=0}^{\mathrm{\infty }}{f}_n\frac{1}{z^n}$; $z\in ℂ$

Inverse Z-Transformation: ${𝒵}^{-1}\{F(z)\}=\frac{1}{2\pi j}{\oint }_C{z}^{n-1}F(z)\mathrm{d}z$

Siehe auch Vorlage:W

Wavelets

• 
  Ingrid Chantal Daubechies (belgische Physikerin und Mathematikerin, geb. 1954)
• 
  Yves François Meyer (französischer Mathematiker, geb. 1939)
• 
  Alfréd Haar (ungarischer Mathematiker, 1885-1933)

Wavelet ... kleine Welle

Haar-Wavelet

$\psi =\{\begin{array}{cc}1;& 0\le x<\frac{1}{2}\\ -1;& \frac{1}{2}\le x\le 1\\ 0;& \text{sonst}\end{array}$

Python-Code:

import numpy as np
import matplotlib.pyplot as plt

def psi(x):
    u = np.zeros(len(x))
    k = 0
    
    for i in x:
        if 0<=i<=1/2:
            u[k] = 1
        elif 1/2<=i<=1:
            u[k] = -1
        else:
            u[k] = 0
        k += 1
    return u

x = np.arange(-0.5, 1.5, 0.001)
y = psi(x)

plt.plot(x,y)
plt.grid()
plt.show()

Mexikanischer Hut

$\psi =(1-x^2)\exp (\frac{-x^2}{2})$

Python-Code:

import numpy as np
import matplotlib.pyplot as plt

def psi(x):
    p = (1-x**2) * np.exp(-x**2 / 2)
    return p

x = np.arange(-np.pi, np.pi, 0.01)
y = psi(x)

plt.plot(x,y)
plt.grid()
plt.show()

Morlet-Wavelet

Daubechies-Wavelet

Wavelet-Transformation

${𝒲}_{\psi }x(a,b)=\frac{1}{\sqrt{a}}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\overline{\psi \left(\frac{t-b}{a}\right)}x(t)dt.$

• $x(t)$: zu transformierende Funktion
• $\psi (t)$: Wavelet-Funktion
• $b$: Translationsparameter
• $a$: Skalierungsparameter

Siehe auch Vorlage:W, Vorlage:W, Vorlage:W, PyWavelets

Gedruckte Bücher (auszugsweise)

• Burg, Haf, Wille, Meister: Höhere Mathematik für Ingenieure, Band I: Analysis. 9. Auflage, Vieweg+Teubner, 2011, ISBN 978-3-8348-1218-6
• Burg, Haf, Wille, Meister: Höhere Mathematik für Ingenieure, Band III: Gewöhnliche Differentialgleichungen, Distributionen, Integraltransformationen. 5. Auflage, Vieweg+Teubner, 2009, ISBN 978-3-8348-0565-2
• Hanke-Bourgeois: Grundlagen der Numerischen Mathematik und des Wissenschaftlichen Rechnens. 3. Aufl., Vieweg+Teubner, 2009, ISBN 978-3-8348-0708-3
• Preuß: Funktionaltransformationen. 2. Aufl., Hanser, 2014, ISBN 978-3-446-41787-8
• Thuselt, Gennrich: Praktische Mathematik mit MATLAB, Scilab und Octave. Springer, 2013, ISBN 978-3-642-25824-4

Vorlage:Navigation zurückhochvor buch