高速フーリエ変換

モジュールのインポート

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'])

フーリエ変換（時間領域と周波数領域）

eps = 1.0e-10
n = 32
print('データ数 n:',n)

データ数 n: 32

tspan = 20.0
dt = tspan/n
tmin = -tspan/2.0
tmax = -tmin-dt
t = np.linspace(tmin,tmax,n)
print('tの最小値 tmin =',t[0],', 最大値 tmax =',t[n-1],', データ間隔 dt =',dt)

tの最小値 tmin = -10.0 , 最大値 tmax = 9.375 , データ間隔 dt = 0.625

df = 1.0/tspan
fmin = -df*n/2
fmax = -fmin-df
f = np.linspace(fmin,fmax,n)
print('fの最小値 fmin =',f[0],', 最大値 fmax =',f[n-1],', データ間隔 df =',df)

fの最小値 fmin = -0.8 , 最大値 fmax = 0.75 , データ間隔 df = 0.05

tr = 4.0 # 方形関数の範囲
a = 1.25 # 方形関数の振幅
z = np.where((t>=-tr) & (t<tr), a, 0.0) # 方形関数値の計算
zs = np.fft.ifftshift(z) # データ逆シフト
#Zs = dt*np.fft.fft(zs) # FFT
Zs = dt*np.fft.ifft(zs)*n # iFFT
Z = np.fft.fftshift(Zs) # データのシフト

nn = 100*n
ttmin = tmin
ttmax = -ttmin
tt = np.linspace(ttmin,ttmax,nn)
zz = np.where((tt>=-tr) & (tt<tr), a, 0.0) # 方形関数値の計算
ffmin = fmin
ffmax = -ffmin
ff = np.linspace(ffmin,ffmax,nn)
ff = ff + (ff==0)*eps # ゼロのときだけeps を加える
fc = 2.0*np.pi*tr*ff
ZZ = 2.0*a*tr*np.sin(fc)/fc # 方形関数のフーリエ変換

Z

array([ 0.78125   +0.j,  0.22788222+0.j, -0.66231177+0.j, -0.63108889+0.j,
        0.32360435+0.j,  0.88158444+0.j,  0.18330714+0.j, -0.8913217 +0.j,
       -0.78125   +0.j,  0.58052112+0.j,  1.37919286+0.j,  0.16244468+0.j,
       -1.88610435+0.j, -1.70735761+0.j,  2.22481177+0.j,  7.62733574+0.j,
       10.15625   +0.j,  7.62733574+0.j,  2.22481177+0.j, -1.70735761+0.j,
       -1.88610435+0.j,  0.16244468+0.j,  1.37919286+0.j,  0.58052112+0.j,
       -0.78125   +0.j, -0.8913217 +0.j,  0.18330714+0.j,  0.88158444+0.j,
        0.32360435+0.j, -0.63108889+0.j, -0.66231177+0.j,  0.22788222+0.j])

Z.real

array([ 0.78125   ,  0.22788222, -0.66231177, -0.63108889,  0.32360435,
        0.88158444,  0.18330714, -0.8913217 , -0.78125   ,  0.58052112,
        1.37919286,  0.16244468, -1.88610435, -1.70735761,  2.22481177,
        7.62733574, 10.15625   ,  7.62733574,  2.22481177, -1.70735761,
       -1.88610435,  0.16244468,  1.37919286,  0.58052112, -0.78125   ,
       -0.8913217 ,  0.18330714,  0.88158444,  0.32360435, -0.63108889,
       -0.66231177,  0.22788222])

fig = plt.figure() # グラフ領域の作成
plt.plot(tt,zz,'r')
plt.stem(t,z,'o')
plt.title('Time domain')
plt.xlabel(r'$t$')
plt.ylabel(r'$h(t)$')
fig.savefig('p1_ex10_1.pdf')

fig = plt.figure() # グラフ領域の作成
plt.plot(ff,ZZ,'r')
plt.stem(f,Z.real,'o')
plt.title('Frequency domain')
plt.xlabel(r'$f$')
plt.ylabel(r'$H(f)$')
fig.savefig('p1_ex10_2.pdf')