[Calcul scientifique / technique par Python] Transformation de Fourier à grande vitesse discrète en 3D unidimensionnelle, scipy

introduction

La transformation de Fourier est une technique indispensable pour mener des recherches scientifiques et techniques. ** Dans cet article, scipy module scipy.fftpack [ La méthode de transformation de Fourier discrète à grande vitesse utilisant 1] est résumée par un exemple simple. ** ** Pour la méthode de transformation de Fourier discrète ordinaire, veuillez vous référer à la littérature le cas échéant car il existe un texte élémentaire [2,3].

** Dans cet article, l'utilisation de la bibliothèque de conversion de Fourier à grande vitesse de scipy est mise en avant en tenant compte de la scène d'utilisation sur le lieu de travail [4]. Si vous vous entraînez en vérifiant le fonctionnement de l'exemple (1), vous pourrez bientôt l'utiliser sur le terrain. ** **

Contenu

(1) Exemple de transformée de Fourier discrète unidimensionnelle (2) Exemple de transformée de Fourier 3D

Exemple (1): Transformée de Fourier discrète unidimensionnelle

Dans cet exemple, la transformée de Fourier $ g (\ omega) $ de la fonction gaussienne $ f (t) = e ^ {(- (t-5) ^ 2)} $ est obtenue. Seules les deux méthodes suivantes sont utilisées.

1. scipy.fftpack. ** fft ** Transformée de Fourier
2. scipy.fftpack. ** ift ** Transformée de Fourier inverse

Soit le nombre total de données N égal à 40 et la largeur d'échantillonnage T comprise entre $ t = 0 $ et $ t = 10 $ ($ T = 10 $).

Par conséquent, l'échantillonnage temporel

t_n \ = \frac{nT}{N} (n =0, 1, ..., N-1)

Sera.

Parallèlement à cela, les points discrets de la fréquence $ \ omega $

\omega_k = \frac{2k\pi}{T} (k=1,2,3,...,N)

Sera.

Le déroulement du programme est le suivant.

1. La fonction gaussienne $ f (t) = e ^ {(- (t-5) ^ 2)} $ est uniformément échantillonnée à 40 points de t = 0 à 10.
2. Utilisez ** fft pour le transformer en une transformée de Fourier discrète complexe pour obtenir $ g (\ omega_n) $. ** **
3. Visualisez-le (** Figure 1 **). Il est divisé en une partie réelle et une partie imaginaire. 4.Spectre de puissance|g(\omega_n)|^2Est calculé et illustré(Figure 2)。
4. ** Convertissez à nouveau $ g (\ omega_n) $ par transformée de Fourier (transformée de Fourier inverse) pour obtenir la fonction $ ff (t) $. ** **
5. Comparez $ ff (t) $ avec $ f (t) $ (** Figure 3 **). Reproduisez strictement les points de données donnés.

** La transformation de Fourier à grande vitesse de la fonction souhaitée est possible en changeant la forme fonctionnelle de $ f (t) $. ** **

FFT1D.py

""" 
Transformée de Fourier discrète unidimensionnelle
"""
import numpy as np
from scipy.fftpack import fft, ifft
import matplotlib.pyplot as plt


#Création d'exemples de données de séries chronologiques
N = 40                      #Le nombre de données
T=10            #Largeur d'échantillonnage
del_t= T/N   #Intervalle d'échantillonnage

del_w=2*np.pi/T #Intervalle de fréquence de la transformée de Fourier discrète
#

#Génération de points discrets
t = np.arange(0,T-del_t,del_t)
w=np.arange(2*np.pi/T, 2*np.pi*N/T, del_w)

#

f = np.exp(-(t-5)**2) # #Fonction qui donne des exemples de données

#
g = fft(f)

g_real=np.real(g)
g_imag=np.imag(g)

#plot
plt.xlabel('w', fontsize=24)
plt.ylabel('g(w)', fontsize=24)

plt.plot(w,g_real,marker='o', markersize=4,label='Fourier transform: Real part')
plt.plot(w,g_imag,marker='o',markersize=4, label='Fourier transform: Imaginary part')

plt.legend(loc='best')

plt.show()


#Affichage du spectre de puissance

plt.plot(w,np.abs(g)**2,marker='o',markersize=4,label='|g(w)|^2')

plt.xlabel('w', fontsize=24)
plt.ylabel('Power spectrum |g(w)|^2', fontsize=24)
plt.legend(loc='best')

plt.show()


#Transformée de Fourier discrète inverse

ff = ifft(g)
plt.plot(t, np.real(ff), label='Fourier inverse transform: Real part')
plt.plot(t, np.imag(ff), label='Fourier inverse transform: Imaginary part')

plt.plot(t, f,'o',markersize=4,label='Raw data')


plt.xlabel('t', fontsize=24)
plt.ylabel('f(t)', fontsize=24)
plt.legend(loc='best')

plt.show()

Résultat (1) Transformée de Fourier discrète à grande vitesse unidimensionnelle

Figure 1. Transformée de Fourier discrète $ g (\ omega_n) $

Figure 2.Spectre de puissance|g(\omega_n)|^2

Figure 3. Comparaison entre la transformée de Fourier complexe inverse (ligne) et les données échantillonnées (cercle vert). Puisque je pensais à une fonction réelle, la partie imaginaire de la transformée de Fourier complexe inverse est nulle.

Exemple (2) Transformée de Fourier discrète à grande vitesse tridimensionnelle

Dans cet exemple, la fonction gaussienne tridimensionnelle $ f (t_x, t_y, t_z) = e ^ {- [(t_x-5) ^ 2 + (t_y-5) ^ 2 + (t_z-5) ^ 2)]} $ Effectue une transformée de Fourier discrète complexe à haute vitesse. Par rapport à l'exemple (1), qui était un problème unidimensionnel, il y a plus d'indices, etc., mais le contenu n'est pas significativement difficile.

Comme pour la méthode unidimensionnelle, seules les deux méthodes suivantes sont utilisées pour le multidimensionnel (deux ou plus).

1. scipy.fftpack. ** fftn ** Transformée de Fourier à n dimensions
2. scipy.fftpack. ** iftn ** Transformée de Fourier inverse à n dimensions

** La dimension de la transformation de Fourier n'a pas besoin d'être très nerveuse du côté du programmeur car scipy la déterminera à partir de la taille du tableau lors de l'exécution de la transformation de Fourier **

Cet exemple est destiné uniquement au contrôle de fonctionnement et n'est pas visualisé. Mais si nécessaire, l'intrigue souhaitée pourrait facilement être réalisée.

Le déroulement du programme est le suivant.

1. Fonction gaussienne $ f (t_x, t_y, t_z) = e ^ {- [(t_x-5) ^ 2 + (t_y-5) ^ 2 + (t_z-5) ^ 2)]} , "3 Les temps dimensionnels " tx, ty, tz $ sont uniformément échantillonnés en 4 points de t = 0 à 10 respectivement (Nx, Ny, Nz = 4, 4, 4).
2. Utilisez ** fftn pour le transformer en une transformée de Fourier discrète complexe pour obtenir $ g (\ omega_ {x_n}, \ omega_ {y_m}, \ omega_ {z_k}) $. ** **
3. ** Transformée de Fourier inverse en utilisant ifftn pour obtenir la fonction $ f (t_x, t_y, t_z) $. ** **
4. Essayez de sortir la valeur numérique.

""" 
Transformée de Fourier discrète multidimensionnelle
12 août 2017
"""
import numpy as np
from scipy.fftpack import fftn, ifftn #Transformée de Fourier discrète à n dimensions / inverse de Fourier
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


#Création d'exemples de données de séries chronologiques
Nx = 4        #Nombre de données dans la direction Nx
Ny = 4          #Nombre de données dans la direction Ny
Nz = 4          #Nombre de données dans la direction Nz
Tx=10            #Largeur d'échantillonnage
Ty = 10
Tz = 10

del_t_x= Tx/Nx   #Intervalle d'échantillonnage dans la direction tx
del_t_y= Ty/Ny   #Intervalle d'échantillonnage dans le sens ty
del_t_z= Tz/Nz   #Intervalle d'échantillonnage dans le sens ty


del_w_x=2*np.pi/Tx #Espacement des fréquences des composantes x de la transformée de Fourier discrète
del_w_y=2*np.pi/Ty #Espacement des fréquences des composantes y de la transformée de Fourier discrète
del_w_z=2*np.pi/Tz #Espacement de la fréquence de composante z de la transformée de Fourier discrète

#
t_x, t_y, t_z= np.meshgrid(np.arange(0,Tx-del_t_x,del_t_x),np.arange(0,Ty-del_t_y,del_t_y), np.arange(0,Tz-del_t_z,del_t_z))  #Génération de maillage
w_x, w_y, w_z= np.meshgrid(np.arange(2*np.pi/Tx, 2*np.pi*Nx/Tx, del_w_x),np.arange(2*np.pi/Ty, 2*np.pi*Ny/Ty, del_w_y), np.arange(2*np.pi/Tz, 2*np.pi*Nz/Tz, del_w_z))  #Génération de maillage


#
f = np.exp(-((t_x-5)**2+(t_y-5)**2+(t_z-5)**2)) #Fonction qui donne des exemples de données

#
g = fftn(f)  #Transformée de Fourier rapide multidimensionnelle:3D dans ce cas
g_real=np.real(g) #Partie réelle
g_imag=np.imag(g) #Partie imaginaire


#Transformée de Fourier discrète multidimensionnelle inverse
ff = ifftn(g)
ff_real=np.real(ff)#Partie réelle
ff_imag=np.imag(ff) #Partie imaginaire

Résultat (2): Transformée de Fourier discrète rapide 3D

Sortons une petite valeur numérique.

print(g[1,2,1])

out: (-0.500000003576+0.862688207649j)

Cela signifie que la composante [1,2,1] de la transformée de Fourier g est -0,5 + 0,86i (nombre complexe).

Ensuite, vérifions la transformée de Fourier inverse.

print(ff[1,2,2])

out: (0.00193045413623+0j)

C'est un vrai nombre. Je peux calculer sans aucun problème.

Les références

1. Page Web de scipy module scipy.fftpack
2. Kenichi Kanaya, ["Cours de mathématiques appliquées que vous pouvez comprendre"](https://www.amazon.co.jp/%E3%81%93%E3%82%8C%E3%81%AA%E3%82 % 89% E5% 88% 86% E3% 81% 8B% E3% 82% 8B% E5% BF% 9C% E7% 94% A8% E6% 95% B0% E5% AD% A6% E6% 95% 99 % E5% AE% A4% E2% 80% 95% E6% 9C% 80% E5% B0% 8F% E4% BA% 8C% E4% B9% 97% E6% B3% 95% E3% 81% 8B% E3 % 82% 89% E3% 82% A6% E3% 82% A7% E3% 83% BC% E3% 83% 96% E3% 83% AC% E3% 83% 83% E3% 83% 88% E3% 81 % BE% E3% 81% A7-% E9% 87% 91% E8% B0% B7-% E5% 81% A5% E4% B8% 80 / dp / 4320017382), Kyoritsu Publishing, 2003.
3. Strang, ["Linear Algebra Introduction"](https://www.amazon.co.jp/%E4%B8%96%E7%95%8C%E6%A8%99%E6%BA%96MIT%E6 % 95% 99% E7% A7% 91% E6% 9B% B8-% E3% 82% B9% E3% 83% 88% E3% 83% A9% E3% 83% B3% E3% 82% B0-% E7 % B7% 9A% E5% BD% A2% E4% BB% A3% E6% 95% B0% E3% 82% A4% E3% 83% B3% E3% 83% 88% E3% 83% AD% E3% 83 % 80% E3% 82% AF% E3% 82% B7% E3% 83% A7% E3% 83% B3-% E3% 82% AE% E3% 83% AB% E3% 83% 90% E3% 83% BC% E3% 83% 88 / dp / 4764904055), Modern Science Co., Ltd., 2015.
4. Mark Newman, "Computational Physics", Createspace Independent Publishing Platform, 2012.

Vous pouvez obtenir des informations dans de nombreux livres et sites Web sur les transformées de Fourier discrètes ordinaires «pas rapides». Les références 2 et 3 font également référence aux transformées de Fourier rapides. Il est très facile à lire grâce à un examen attentif du lecteur. 4 est un livre étranger, mais l'explication avec l'essence vaut la peine d'être lue.