[PYTHON] FFT (Fast Fourier Transform): formules et exemples de mise en œuvre

Principes et formules de la FFT

Puisqu'il n'y avait pas d'article portant sur [Formule → Implémentation], je vais l'épeler. Cette fois, je vais vous présenter la FFT de type Cooley - Tukey la plus simple appelée radix-2.

Premièrement, la DFT (transformée de Fourier discrète) est exprimée par l'équation suivante. $F(t) = \displaystyle \sum_{x=0}^{N-1}f(x) e^{-i \frac{2πtx}{N}} \qquad t=0,1,...,N-1\tag{1}$ Ici, $ N $ est la taille d'entrée et est limité à la puissance de ** 2. ** ** Définissez $ W_N = e ^ {-i \ frac {2π} {N}} $ et réécrivez comme suit. $F(t) = \displaystyle \sum_{x=0}^{N-1}f(x) W_{N}^{tx}\qquad t=0,1,...,N-1\tag{2}$ La DFT effectuée uniquement sur l'élément pair de l'entrée est exprimée par la formule suivante. $F_{even}(t) = \displaystyle \sum_{k=0}^{\frac{N}{2}-1}f(2k) W_\frac{N} {2}^{tk}\tag{3}$ De même, la DFT effectuée uniquement sur l'élément impair est exprimée par la formule suivante. $F_{odd}(t) = \displaystyle \sum_{k=0}^{\frac{N}{2}-1}f(2k+1) W_\frac{N}{2}^{tk}\tag{4}$ En utilisant $ (3) (4) $, $ (2) $ peut être réécrit comme suit. $F(t) = F_{even}(t) +W_N^t F_{odd}(t)\tag{5}$ Jusqu'à présent, vous venez de transformer le DFT en parties impaires et impaires. Maintenant, voici le point important de FFT. En fait, $ F_ {pair} (t), F_ {impair} (t), W_N ^ t $ ont la périodicité suivante.

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}

En utilisant $ (6) (7) (8) $ pour $ (5) $, le DFT peut enfin être réécrit comme suit.

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}

Puisque $ F_ {pair} (t) $ et $ F_ {impair} (t) $ sont communs, $ F (t) $ et $ F (t + \ frac {N} {2}) $ sont calculés en même temps. C'est pourquoi vous le voulez. C'est le mécanisme pour réduire la quantité de calcul.

Sur la base de ce qui précède, passons à la mise en œuvre.

Programme FFT

Écrivons une fonction FFT en Python. Voici l'idée de base. Divisez les données d'entrée en éléments pairs et en éléments impairs, effectuez une FFT (récursive) comme $ (3) (4) $, puis utilisez $ (9) (10) $ à 0 ~ (N / 2) - Calculez le 1er et N / 2 ~ N-1th. Par conséquent, FFT est une fonction récursive.

def FFT(f):
    #f:Données d'entrée de taille N
    N = len(f)
    
    if N == 1:    #Lorsque N vaut 1, les données d'entrée sont renvoyées telles quelles
        return f[0]

    f_even = f[0:N:2]    #Even-th élément de f
    f_odd = f[1:N:2]    #L'élément impair de f
    F_even = FFT(f_even)    #(3)FFT sur élément pair
    F_odd = FFT(f_odd)    #(4)FFT sur élément pair

    W_N = np.exp(-1j * (2 * np.pi * np.arange(0, N // 2)) / N)    #t vaut 0~N/2-Un tableau qui calcule W jusqu'au premier

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

    return F

Contrôle de fonctionnement

Vérifions si nous pouvons obtenir le même résultat que la FFT de numpy. Créez une onde composite de deux ondes sinusoïdales avec les paramètres suivants et utilisez-la comme entrée.

import numpy as np
import matplotlib.pyplot as plt

N = 64 #Le nombre de données
n = np.arange(N)
f1 = 2 #Cycle 1
f2 = 5 #Cycle 2
A1 = 3 #Amphithéâtre 1
A2 = 5 #Amphithéâtre 2
INPUT_f = A1 * np.sin(f1 * 2 * np.pi * (n/N)) + A2 * np.sin(f2 * 2 * np.pi * (n/N)) #Synthèse de deux ondes de péché

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

Exécutez FFT. Comme la sortie de FFT est un nombre complexe, voici une visualisation de la composante d'amplitude extraite en la multipliant par une valeur absolue.

OUTPUT_F_numpy = np.fft.fft(INPUT_f)    #version numpy
OUTPUT_F_original = FFT(INPUT_f)    #Version faite maison

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

plt.subplot(1,2,1)    #Affichage des composants d'oscillation(version numpy)
plt.title('numpy_ver')
plt.plot(np.abs(OUTPUT_F_numpy))

plt.subplot(1,2,2)    #Affichage des composants d'oscillation(Version faite maison)
plt.title('original')
plt.plot(np.abs(OUTPUT_F_original))

J'ai obtenu le même résultat en toute sécurité.