高速フーリエ変換（波数領域と空間領域）

モジュールのインポート

import numpy as np  
import matplotlib.pyplot as plt

import scienceplots
#Warning : As of version 2.0.0, you need to add import scienceplots before setting the style (plt.style.use('science')).
plt.style.use(['science', 'notebook'])

# rcPramsを使ってTeXのフォントをデフォルトにする
#plt.rcParams['font.family'] = 'Times New Roman' # font familyの設定
plt.rcParams['font.family'] = 'serif' # font familyの設定
plt.rcParams['mathtext.fontset'] = 'cm' # math fontの設定

波数領域と空間領域

pi2 = 2.0*np.pi
eps = 1.0e-10
n = 32
span = 20.0
diameter = 4.0
print('データ数 n:',n)

データ数 n: 32

n_x = n
x_span = span
d_x = x_span/n_x
x_min = -x_span/2.0
x_max = -x_min-d_x
x = np.linspace(x_min, x_max, n_x)
print('xmin =',x[0],', xmax =',x[n_x-1],', d_x =',d_x)

xmin = -10.0 , xmax = 9.375 , d_x = 0.625

d_kx = pi2/x_span
kx_min = -d_kx*n_x/2
kx_max = -kx_min-d_kx
kx = np.linspace(kx_min, kx_max, n_x)
print('kx_min =',kx[0],', kx_max =',kx[n_x-1],', d_kx =',d_kx)

kx_min = -5.026548245743669 , kx_max = 4.71238898038469 , d_kx = 0.3141592653589793

xr = diameter # 方形関数の範囲
a = 1.25 # 方形関数の振幅
f = np.where((x>=-xr) & (x<xr), a, 0.0) # 方形関数値の計算
fs = np.fft.ifftshift(f) # データ逆シフト
gs = d_x*np.fft.fft(fs) # FFT
# gs = d_x*np.fft.ifft(fs) # iFFT
g = np.fft.fftshift(gs) # データのシフト
g = g/pi2

nn = 100*n
xx_min = x_min
xx_max = -xx_min
xx = np.linspace(xx_min, xx_max, nn)
ff = np.where((xx>=-xr) & (xx<xr), a, 0.0) # 方形関数値の計算
kkx_min = kx_min
kkx_max = -kkx_min
kkx = np.linspace(kkx_min, kkx_max, nn)
kkx = kkx + (kkx==0)*eps # ゼロのときだけepsを加える
u = xr*kkx
gg = 2.0*a*xr*np.sin(u)/u # 方形関数のフーリエ変換
gg = gg/pi2

fig = plt.figure() # グラフ領域の作成
plt.plot(xx,ff,'r')
plt.stem(x,f,'o')
plt.title('Space domain')
plt.xlabel(r'$x$')
plt.ylabel(r'$f(x)$')
fig.savefig('p2_ex10_1.pdf')

fig = plt.figure() # グラフ領域の作成
plt.plot(kkx,gg,'r')
plt.stem(kx,g.real,'o')
plt.title('wavenumber domain')
plt.xlabel(r'$k_x$')
plt.ylabel(r'$g(k_x)$')
fig.savefig('p2_ex10_2.pdf')