[PYTHON] FFT (Fast Fourier Transform): Formeln und Beispiele für die Implementierung

FFT-Prinzipien und -Formeln

Da es keinen Artikel gab, der sich mit [Formel → Implementierung] befasste, werde ich ihn buchstabieren. Dieses Mal werde ich die einfachste Cooley - Tukey Typ FFT namens Radix-2 vorstellen.

Zunächst wird DFT (diskrete Fourier-Transformation) durch die folgende Gleichung ausgedrückt. $F(t) = \displaystyle \sum_{x=0}^{N-1}f(x) e^{-i \frac{2πtx}{N}} \qquad t=0,1,...,N-1\tag{1}$ Hier ist $ N $ die Eingabegröße und auf die Potenz von ** 2 begrenzt. ** ** ** Setze $ W_N = e ^ {-i \ frac {2π} {N}} $ und schreibe wie folgt um. $F(t) = \displaystyle \sum_{x=0}^{N-1}f(x) W_{N}^{tx}\qquad t=0,1,...,N-1\tag{2}$ Die DFT, die nur für das gerade Element der Eingabe ausgeführt wird, wird durch die folgende Formel ausgedrückt. $F_{even}(t) = \displaystyle \sum_{k=0}^{\frac{N}{2}-1}f(2k) W_\frac{N} {2}^{tk}\tag{3}$ In ähnlicher Weise wird die DFT, die nur für das ungerade-te Element ausgeführt wird, durch die folgende Formel ausgedrückt. $F_{odd}(t) = \displaystyle \sum_{k=0}^{\frac{N}{2}-1}f(2k+1) W_\frac{N}{2}^{tk}\tag{4}$ Mit $ (3) (4) $ kann $ (2) $ wie folgt umgeschrieben werden. $F(t) = F_{even}(t) +W_N^t F_{odd}(t)\tag{5}$ Bisher haben Sie die DFT in ungerade und ungerade Teile unterteilt und transformiert. Hier ist der wichtige Punkt der FFT. Tatsächlich haben $ F_ {gerade} (t), F_ {ungerade} (t), W_N ^ t $ die folgende Periodizität.

F_{even}(t) = F_{even}(t+\frac{N}{2})\tag{6} F_{odd}(t) = F_{odd}(t+\frac{N}{2})\tag{7} W_N^t = -W_N^{t+\frac{N}{2}}\tag{8}

Durch Verwendung von $ (6) (7) (8) $ für $ (5) $ kann die DFT schließlich wie folgt umgeschrieben werden.

F(t) = F_{even}(t) +W_N^t F_{odd}(t)\qquad t=0,1,...,\frac{N}{2}-1\tag{9} F(t+\frac{N}{2}) = F_{even}(t) - W_N^tF_{odd}(t)\qquad t=0,1,...,\frac{N}{2}-1\tag{10}

Da $ F_ {gerade} (t) $ und $ F_ {ungerade} (t) $ gemeinsam sind, werden $ F (t) $ und $ F (t + \ frac {N} {2}) $ gleichzeitig berechnet. Deshalb willst du es. Dies ist der Mechanismus zum Reduzieren des Rechenaufwands.

Fahren wir auf der Grundlage der obigen Ausführungen mit der Implementierung fort.

FFT-Programm

Schreiben wir eine FFT-Funktion in Python. Hier ist die Grundidee. Teilen Sie die Eingabedaten in gerade und ungerade Elemente, führen Sie eine FFT (rekursiv) wie $ (3) (4) $ durch und verwenden Sie dann $ (9) (10) $ bis 0 ~ (N / 2) - Berechnen Sie die 1. und N / 2 ~ N-1. Daher ist FFT eine rekursive Funktion.

def FFT(f):
    #f:Eingabedaten der Größe N.
    N = len(f)
    
    if N == 1:    #Wenn N 1 ist, werden die Eingabedaten so zurückgegeben, wie sie sind
        return f[0]

    f_even = f[0:N:2]    #Gerades Element von f
    f_odd = f[1:N:2]    #Das ungerade Element von f
    F_even = FFT(f_even)    #(3)FFT am geraden Element
    F_odd = FFT(f_odd)    #(4)FFT am geraden Element

    W_N = np.exp(-1j * (2 * np.pi * np.arange(0, N // 2)) / N)    #t ist 0~N/2-Ein Array, das W bis zum ersten berechnet

    F = np.zeros(N, dtype ='complex')    #FFT-Ausgang
    F[0:N//2] = F_even + W_N * F_odd    #(9)Berechnung(t:0~N/2-1)
    F[N//2:N] = F_even - W_N * F_odd    #(10)Berechnung(t:N/2~N-1)

    return F

Funktionsprüfung

Lassen Sie uns prüfen, ob wir die gleiche Ausgabe wie die FFT von numpy erhalten können. Erstellen Sie eine zusammengesetzte Welle aus zwei Sinuswellen mit den folgenden Parametern und verwenden Sie diese als Eingabe.

import numpy as np
import matplotlib.pyplot as plt

N = 64 #Die Anzahl der Daten
n = np.arange(N)
f1 = 2 #Zyklus 1
f2 = 5 #Zyklus 2
A1 = 3 #Amphitheater 1
A2 = 5 #Amphitheater 2
INPUT_f = A1 * np.sin(f1 * 2 * np.pi * (n/N)) + A2 * np.sin(f2 * 2 * np.pi * (n/N)) #Synthese zweier Sinuswellen

#Anzeige
plt.figure(figsize=(10, 5))
plt.style.use('seaborn')
plt.plot(INPUT_f)

Führen Sie FFT aus. Da die Ausgabe von FFT eine komplexe Zahl ist, ist hier eine Visualisierung der Amplitudenkomponente, die durch Multiplizieren mit einem absoluten Wert extrahiert wurde.

OUTPUT_F_numpy = np.fft.fft(INPUT_f)    #numpy Version
OUTPUT_F_original = FFT(INPUT_f)    #Selbstgemachte Version

plt.figure(figsize=(10, 5))

plt.subplot(1,2,1)    #Anzeige der Oszillationskomponenten(numpy Version)
plt.title('numpy_ver')
plt.plot(np.abs(OUTPUT_F_numpy))

plt.subplot(1,2,2)    #Anzeige der Oszillationskomponenten(Selbstgemachte Version)
plt.title('original')
plt.plot(np.abs(OUTPUT_F_original))

Ich habe das gleiche Ergebnis sicher erhalten.