[Scientific / technical calculation by Python] 1D-3D discrete fast Fourier transform, scipy

Introduction

The Fourier transform is an indispensable technique for carrying out scientific and technical research. ** In this article, scipy's scipy.fftpack module [ The method of high-speed discrete Fourier transform using 1] is summarized through a simple example. ** ** For the usual method of discrete Fourier transform, please refer to the literature as appropriate because there is an elementary text [2,3].

** In this paper, the use of scipy's fast Fourier transform library is brought to the fore in consideration of the usage scene in the workplace [4]. If you practice while checking the operation of Example (1), you will soon be able to use it in the field. ** **

Contents

(1) Example of one-dimensional discrete Fourier transform (2) Example of 3D Fourier transform


Example (1): One-dimensional discrete Fourier transform

In this example, the Fourier transform $ g (\ omega) $ of the Gaussian function $ f (t) = e ^ {(-(t-5) ^ 2)} $ is obtained. Only the following two methods are used.

  1. scipy.fftpack. ** fft ** Fourier transform
  2. scipy.fftpack. ** ift ** inverse Fourier transform

Let the total number of data N be 40 and the sampling width T be from $ t = 0 $ to $ t = 10 $ ($ T = 10 $).

Therefore time sampling

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

Will be.

Along with that, the discrete points of the frequency $ \ omega $

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

Will be.

The flow of the program is as follows.

  1. Sample the Gaussian function $ f (t) = e ^ {(-(t-5) ^ 2)} $ uniformly at 40 points from t = 0 to 10.
  2. Use ** fft to transform it into a complex discrete Fourier transform to get $ g (\ omega_n) $. ** **
  3. Visualize it (** Figure 1 **). It is divided into a real part and an imaginary part. 4.Power spectrum|g(\omega_n)|^2Is calculated and illustrated(Figure 2)。
  4. ** Fourier transform $ g (\ omega_n) $ again (inverse Fourier transform) to get the function $ ff (t) $. ** **
  5. Compare $ ff (t) $ with $ f (t) $ (** Figure 3 **). Strictly reproduce the given data points.

** Fast Fourier transform of the desired function is possible by changing the functional form of $ f (t) $. ** **

FFT1D.py


""" 
One-dimensional discrete Fourier transform
"""
import numpy as np
from scipy.fftpack import fft, ifft
import matplotlib.pyplot as plt


#Creating time-series sample data
N = 40                      #The number of data
T=10            #Sampling width
del_t= T/N   #Sampling interval

del_w=2*np.pi/T #Discrete Fourier Transform Frequency Interval
#

#Discrete point generation
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) # #Function that gives sample data

#
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()


#Power spectrum display

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()


#Inverse discrete Fourier transform

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()

Result (1) One-dimensional high-speed discrete Fourier transform

t1.png

Figure 1. Discrete Fourier Transform $ g (\ omega_n) $

t2.png

Figure 2.Power spectrum|g(\omega_n)|^2

t3.png

Figure 3. Comparison between the inverse complex Fourier transform (line) and the sampling data (green circle). Since I was thinking of a real function, the imaginary part of the inverse complex Fourier transform is zero.


Example (2) Three-dimensional high-speed discrete Fourier transform

In this example, the 3D Gaussian function $ f (t_x, t_y, t_z) = e ^ {-[(t_x-5) ^ 2 + (t_y-5) ^ 2 + (t_z-5) ^ 2)]} $ Is performing a high-speed complex discrete Fourier transform of. Compared to the example (1), which was a one-dimensional problem, the number of subscripts increases, but the content is not significantly difficult.

As with the one-dimensional method, there are only the following two methods used for multidimensional (two or more).

  1. scipy.fftpack. ** fftn ** n-dimensional Fourier transform
  2. scipy.fftpack. ** iftn ** n-dimensional inverse Fourier transform

** The dimension of the Fourier transform does not need to be very nervous on the programmer side because scipy arbitrarily determines the dimension of the array when performing the Fourier transform **

This example is intended only for operation check, and is not visualized. But if needed, the desired plot could easily be made.

The flow of the program is as follows.

  1. Gaussian function $ f (t_x, t_y, t_z) = e ^ {-[(t_x-5) ^ 2 + (t_y-5) ^ 2 + (t_z-5) ^ 2)] Dimensional time "$ tx, ty, tz $ is uniformly sampled at 4 points from t = 0 to 10 respectively (Nx, Ny, Nz = 4, 4, 4).
  2. Use ** fftn to perform a complex discrete Fourier transform on it to get $ g (\ omega_ {x_n}, \ omega_ {y_m}, \ omega_ {z_k}) $. ** **
  3. ** Inverse Fourier transform using ifftn to get the function $ f (t_x, t_y, t_z) $. ** **
  4. Try to output the numerical value.
""" 
Multidimensional discrete Fourier transform
August 12, 2017
"""
import numpy as np
from scipy.fftpack import fftn, ifftn #n-dimensional discrete Fourier transform / inverse Fourier transform
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


#Creating time-series sample data
Nx = 4        #Number of data in Nx direction
Ny = 4          #Number of data in Ny direction
Nz = 4          #Number of data in the Nz direction
Tx=10            #Sampling width
Ty = 10
Tz = 10

del_t_x= Tx/Nx   #Sampling interval in the tx direction
del_t_y= Ty/Ny   #Sampling interval in the ty direction
del_t_z= Tz/Nz   #Sampling interval in the ty direction


del_w_x=2*np.pi/Tx #Interval of frequency of x component of discrete Fourier transform
del_w_y=2*np.pi/Ty #Spacing of the y component of the discrete Fourier transform
del_w_z=2*np.pi/Tz #Interval of frequency of z component of discrete Fourier transform

#
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))  #Mesh generation
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))  #Mesh generation


#
f = np.exp(-((t_x-5)**2+(t_y-5)**2+(t_z-5)**2)) #Function that gives sample data

#
g = fftn(f)  #Multidimensional fast Fourier transform:3D in this case
g_real=np.real(g) #Real part
g_imag=np.imag(g) #Imaginary part


#Inverse multidimensional discrete Fourier transform
ff = ifftn(g)
ff_real=np.real(ff)#Real part
ff_imag=np.imag(ff) #Imaginary part


Result (2): 3D fast discrete Fourier transform

Let's output a little numerical value.


print(g[1,2,1])

out: (-0.500000003576+0.862688207649j)

It means that the [1,2,1] component of the Fourier transformed g is -0.5 + 0.86i (complex number).

Next, let's check the inverse Fourier transform.


print(ff[1,2,2])

out: (0.00193045413623+0j)

This is a real number. I can calculate without any problem.


References

  1. scipy's scipy.fftpack module web page
  2. Kenichi Kanatani, [“Applied Mathematics Classes You Can Understand”](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.

You can get information from many books and websites about the "not fast" ordinary discrete Fourier transform. References 2 and 3 also refer to the Fast Fourier Transform. It is very easy to read due to the careful consideration of the reader. 4 is a foreign book, but the explanation with the essence is worth reading.

Recommended Posts

[Scientific / technical calculation by Python] 1D-3D discrete fast Fourier transform, scipy
[Scientific / technical calculation by Python] 3rd order spline interpolation, scipy
[Scientific / technical calculation by Python] Sum calculation, numerical calculation
[Scientific / technical calculation by Python] Numerical integration, trapezoidal / Simpson law, numerical calculation, scipy
[Scientific / technical calculation by Python] Inverse matrix calculation, numpy
[Scientific / technical calculation by Python] Histogram, visualization, matplotlib
[Scientific / technical calculation by Python] Lagrange interpolation, numerical calculation
[Scientific / technical calculation by Python] Fitting by nonlinear function, equation of state, scipy
[Scientific / technical calculation by Python] Solving (generalized) eigenvalue problem using numpy / scipy, using library
[Scientific / technical calculation by Python] Logarithmic graph, visualization, matplotlib
[Scientific / technical calculation by Python] Polar graph, visualization, matplotlib
[Scientific / technical calculation by Python] Basic operation of arrays, numpy
[Scientific / technical calculation by Python] Monte Carlo integration, numerical calculation, numpy
[Scientific / technical calculation by Python] List of usage of (special) functions used in physics by using scipy
[Scientific / technical calculation by Python] Solving simultaneous linear equations, numerical calculation, numpy
[Scientific / technical calculation by Python] Vector field visualization example, electrostatic field, matplotlib
[Scientific / technical calculation by Python] 2D random walk (drunk walk problem), numerical calculation
[Scientific / technical calculation by Python] Calculation of matrix product by @ operator, python3.5 or later, numpy
[Scientific / technical calculation by Python] Solving ordinary differential equations, mathematical formulas, sympy
[Scientific / technical calculation by Python] Drawing animation of parabolic motion with locus, matplotlib
[Scientific / technical calculation by Python] Solving second-order ordinary differential equations by Numerov method, numerical calculation
[Scientific / technical calculation by Python] Numerical calculation to find the value of derivative (differential)
[Scientific / technical calculation by Python] Analytical solution to find the solution of equation sympy
[Scientific / technical calculation by Python] Plot, visualize, matplotlib 2D data with error bars
[Scientific / technical calculation by Python] Drawing of 3D curved surface, surface, wireframe, visualization, matplotlib
[Scientific / technical calculation by Python] Solving one-dimensional Newton equation by the 4th-order Runge-Kutta method
[Scientific / technical calculation by Python] Plot, visualization, matplotlib of 2D data read from file
[Scientific / technical calculation by Python] Drawing, visualization, matplotlib of 2D (color) contour lines, etc.
[Python] One-liner FFT (Fast Fourier Transform) and other madness
[Scientific / technical calculation by Python] Numerical solution of second-order ordinary differential equations, initial value problem, numerical calculation
[Scientific / technical calculation by Python] List of matrices that appear in Hinpan in numerical linear algebra
[Scientific and technical calculation by Python] Drawing of fractal figures [Sierpinski triangle, Bernsley fern, fractal tree]
[Scientific / technical calculation by Python] Wave "beat" and group velocity, wave superposition, visualization, high school physics
[Scientific / technical calculation by Python] Numerical solution of one-dimensional harmonic oscillator problem by velocity Verlet method
[Scientific / technical calculation by Python] Numerical solution of eigenvalue problem of matrix by power method, numerical linear algebra
[Scientific / technical calculation by Python] Monte Carlo simulation of thermodynamics of 2D Ising spin system by Metropolis method
Scientific / technical calculation by Python] Drawing and visualization of 3D isosurface and its cross section using mayavi
[Scientific / technical calculation by Python] Derivation of analytical solutions for quadratic and cubic equations, mathematical formulas, sympy
Signal processing in Python (1): Fourier transform
Properties of the discrete Fourier transform
[Scientific / technical calculation by Python] Solving the boundary value problem of ordinary differential equations in matrix format, numerical calculation
[Scientific / technical calculation by Python] Solving one-dimensional Schrodinger equation in stationary state by shooting method (1), well-type potential, quantum mechanics