Ing Mathematik: FFT und Wavelets

Zurück zu Näherungsweises Lösen gewöhnlicher Dgl. | Hoch zu Gesamtinhaltsverzeichnis |

Fourierreihen

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

$f(x)\approx {\frac {a_{0}}{2}}+\sum _{k=1}^{\infty }\left(a_{k}\cos \left(kx\right)+b_{k}\sin \left(kx\right)\right)$

$a_{k}={\frac {1}{\pi }}\int _{-\pi }^{\pi }f(x)\cdot \cos \left(kx\right)\mathrm {d} x\quad {\text{für }}k\geq 0$

$b_{k}={\frac {1}{\pi }}\int _{-\pi }^{\pi }f(x)\cdot \sin \left(kx\right)\mathrm {d} x\quad {\text{für }}k\geq 1$

bzw. in komplexer Schreibweise

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

$c_{k}={\frac {1}{2\pi }}\int _{-\pi }^{\pi }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  Fourierreihe

Fourier-Transformation

${\mathcal {F}}{\{f(t)\}}=\int _{-\infty }^{\infty }\mathrm {e} ^{-j\omega t}f(t)\mathrm {d} t$

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

Siehe auch  Fourier-Transformation

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}\quad$; $k=0,\dotsc ,N-1$

Matrix-Vektor-Form: ${\hat {a}}=\underbrace {W} _{\text{Fouriermatrix}}\cdot a\quad$ mit $\quad 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}\quad$ für $\quad k=0,\dotsc ,N-1$

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

Siehe auch  Diskrete Fourier-Transformation

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  Schnelle Fourier-Transformation

Short-time-Fourier-Transformation, STFT

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

Typische Fensterfunktionen:  Hamming-Fenster,  von-Hann-Fenster,  Gauß-Fenster.

Siehe auch  Fensterfunktion,  Kurzzeit-Fourier-Transformation.

Z-Transformation

Z-Transformation: ${\mathcal {Z}}\{f_{n}\}=\sum _{n=0}^{\infty }f_{n}{\frac {1}{z^{n}}}$; $\quad z\in \mathbb {C}$

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

Siehe auch  Z-Transformation

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{cases}1;&\quad 0\leq x<{\frac {1}{2}}\\-1;&\quad {\frac {1}{2}}\leq x\leq 1\\0;&\quad {\text{sonst}}\end{cases}}$

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

${\mathcal {W}}_{\psi }x(a,b)={\frac {1}{\sqrt {a}}}\int _{-\infty }^{\infty }{\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  Wavelet,  Wavelet-Transformation,  Haar-Wavelet, 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

Zurück zu Näherungsweises Lösen gewöhnlicher Dgl. | Hoch zu Gesamtinhaltsverzeichnis |