思うところあって、Python の fft と ifft の使い方の復習をしておく。
基本的な使い方
fftとifft（高速フーリエ変換および逆変換）は、numpyとscipyの
両方に含まれている。このため、双方を使うプログラムを作ってみた。
import numpy as np
from scipy import fftpack
def fft_n(nn,x):
sp=np.fft.fft(x)/nn
spa=np.sqrt(sp.real**2+sp.imag**2)
wv=np.fft.ifft(sp*nn)
return sp,spa,wv
def fft_s(nn,x):
sp=fftpack.fft(x)/nn
spa=np.sqrt(sp.real**2+sp.imag**2)
wv=fftpack.ifft(sp*nn)
return sp,spa,wv
def main():
x=np.random.random_sample(10)
nd=len(x)
nn=2
while nn<nd:
nn=nn*2
xx=np.zeros(nn)
xx[0:nd]=xx[0:nd]+x[0:nd]
for iii in [1,2]:
x=xx
if iii==1:
tstr='* numpy-fft'
sp,spa,wv=fft_n(nn,x)
if iii==2:
tstr='* scipy-fft'
sp,spa,wv=fft_s(nn,x)
print(tstr)
print('{0:6s} {1:6s} {2:6s} {3:6s} {4:6s} {5:6s}'\
.format('x','Re(Ck)','Im(Ck)','Abs(Ck)','Re(x)','Re(x)'))
for i in range(0,nn):
print('{0:6.3f} {1:6.3f} {2:6.3f} {3:6.3f} {4:6.3f} {5:6.3f}'\
.format(xx[i],sp.real[i],sp.imag[i],spa[i],wv.real[i],wv.imag[i]))
print()
#---------------
# Execute
#---------------
if __name__ == '__main__': main()
出力は以下の通り。当然のことながら同じ結果となっている。
* numpy-fft
x Re(Ck) Im(Ck) Abs(Ck) Re(x) Re(x)
0.016 0.218 0.000 0.218 0.016 0.000
0.534 -0.062 -0.114 0.130 0.534 0.000
0.018 0.030 0.019 0.035 0.018 0.000
0.420 -0.015 -0.032 0.035 0.420 0.000
0.448 0.015 -0.011 0.019 0.448 0.000
0.162 0.011 -0.014 0.018 0.162 0.000
0.815 -0.007 -0.080 0.081 0.815 0.000
0.281 -0.083 0.016 0.085 0.281 0.000
0.616 0.021 0.000 0.021 0.616 0.000
0.178 -0.083 -0.016 0.085 0.178 0.000
0.000 -0.007 0.080 0.081 -0.000 0.000
0.000 0.011 0.014 0.018 -0.000 0.000
0.000 0.015 0.011 0.019 -0.000 0.000
0.000 -0.015 0.032 0.035 0.000 0.000
0.000 0.030 -0.019 0.035 0.000 0.000
0.000 -0.062 0.114 0.130 -0.000 0.000
* scipy-fft
x Re(Ck) Im(Ck) Abs(Ck) Re(x) Re(x)
0.016 0.218 -0.000 0.218 0.016 0.000
0.534 -0.062 -0.114 0.130 0.534 0.000
0.018 0.030 0.019 0.035 0.018 0.000
0.420 -0.015 -0.032 0.035 0.420 0.000
0.448 0.015 -0.011 0.019 0.448 0.000
0.162 0.011 -0.014 0.018 0.162 -0.000
0.815 -0.007 -0.080 0.081 0.815 0.000
0.281 -0.083 0.016 0.085 0.281 0.000
0.616 0.021 -0.000 0.021 0.616 0.000
0.178 -0.083 -0.016 0.085 0.178 0.000
0.000 -0.007 0.080 0.081 -0.000 0.000
0.000 0.011 0.014 0.018 -0.000 0.000
0.000 0.015 0.011 0.019 0.000 0.000
0.000 -0.015 0.032 0.035 0.000 -0.000
0.000 0.030 -0.019 0.035 0.000 0.000
0.000 -0.062 0.114 0.130 0.000 0.000
少し応用的な使い方
以下のように４つの周期を持つサイン波を重ね合わせた波を考える。
\begin{align}
x&=1.0\cdot\sin\left(\cfrac{2\cdot\pi\cdot t}{0.1}\right)\\
&+0.5\cdot\sin\left(\cfrac{2\cdot\pi\cdot t}{1.0}\right)\\
&+1.0\cdot\sin\left(\cfrac{2\cdot\pi\cdot t}{10.0}\right)\\
&+2.0\cdot\sin\left(\cfrac{2\cdot\pi\cdot t}{100.0}\right)
\end{align}
この波形をフーリエ変換して、周期100秒のサイン波以外の成分を取り去って、フーリエ逆変換をかけてみる。
結果は以下の通り。
時刻歴の比較
フーリエスペクトルの比較
プログラム
フーリエ変換と逆変換にはnumpyのライブラリを使っている。これは、iPadで用いるPythonista 3などではscipyは使えないため、移植することを考えると、できるだけscipyを使わないほうがいいため。本来ならできるだけscipyを使いたいところなのだがしかたがない。
フーリエ変換で求まった複素フーリエ係数を扱う場合、折れ曲がり点（ナイキスト振動数）の向こう側も対象となるよういじっておく必要があることに注意。
フーリエ逆変換では、複素フーリエ係数を入力して波形の時刻歴に戻すわけであるが、ifftの戻り値は複素数なので、プロット用には実数部のみを取り出す必要があることに注意。
上で紹介したフーリエスペクトルの図ではscipy.signal.find_peaksにより、ピーク周期をいれるようにした。プログラム例では、ピークの値が50を超えるものにその周期を記載するようにしている。図中、先性データの周期は100秒であるが、フーリwスペクトルのピークは109秒となっている。これはおそらく周期100秒の波の個数が少ないためではないかと思う。（未検証）
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import find_peaks
def drawfig1(tt,xx1,xx2):
fsz=10
xmin=0 ; xmax=300; dx=50
ymin=-5; ymax=5; dy=1
iw=10
ih=3
plt.figure(figsize=(iw,ih*2),facecolor='w')
plt.rcParams['font.size']=fsz
plt.rcParams['font.family']='sans-serif'
for iii in [1,2]:
if iii==1:
tstr='Original wave'
x=tt
y=xx1[0:len(tt)]
if iii==2:
tstr='After waveform processing'
x=tt
y=xx2[0:len(tt)]
nplt=210+iii
plt.subplot(nplt)
plt.xlim([xmin,xmax])
plt.ylim([ymin,ymax])
plt.xlabel('Time (sec)')
plt.ylabel('Stress (MPa)')
plt.xticks(np.arange(xmin,xmax+dx,dx))
plt.yticks(np.arange(ymin,ymax+dy,dy))
plt.grid(color='#999999',linestyle='solid')
plt.title(tstr,loc='left',fontsize=fsz)
plt.plot(x,y,'-',color='#000080',lw=0.5)
plt.tight_layout()
fnameF='fig_wave1.jpg'
plt.savefig(fnameF, dpi=200, bbox_inches="tight", pad_inches=0.1)
def drawfig2(fk1,spa1,fk2,spa2):
fsz=10
xmin=0.01; xmax=10000; dx=1
ymin=0; ymax=500; dy=50
iw=6
ih=4
plt.figure(figsize=(iw*2,ih),facecolor='w')
plt.rcParams['font.size']=fsz
plt.rcParams['font.family']='sans-serif'
for iii in [1,2]:
if iii==1:
tstr='Original'
x=1/fk1[1:]; y=spa1[1:len(x)+1]
if iii==2:
tstr='After waveform processing'
x=1/fk2[1:]; y=spa2[1:len(x)+1]
nplt=120+iii
plt.subplot(nplt)
plt.xlim([xmin,xmax])
plt.ylim([ymin,ymax])
plt.xlabel('Period (sec)',fontsize=fsz)
plt.ylabel('Fourier Spectrunm (MPa*sec)',fontsize=fsz)
plt.xticks(np.arange(xmin,xmax+dx,dx))
plt.yticks(np.arange(ymin,ymax+dy,dy))
plt.grid(color='#999999',linestyle='solid')
plt.title(tstr,loc='left',fontsize=fsz)
plt.xscale('log')
plt.plot(x,y,color='#000080',lw=1)
id=find_peaks(y,height=50)[0]
for i in id:
ss=' {0:.2f}sec'.format(x[i])
plt.text(x[i],y[i],ss,ha='center',va='bottom',fontsize=fsz-2,rotation=90)
plt.tight_layout()
fnameF='fig_wave2.jpg'
plt.savefig(fnameF, dpi=200, bbox_inches="tight", pad_inches=0.1)
def nfft_n(x):
nd=len(x)
nn=2
while nn<nd:
nn=nn*2
xx=np.zeros(nn)
xx[0:nd]=xx[0:nd]+x[0:nd]
sp=np.fft.fft(xx)/nn # complex number
return nn,sp
def ifft_n(nn,sp):
wv=np.fft.ifft(sp*nn) # complex number
return wv.real
def main():
# preparation of test wave
dt=0.01
tt=np.arange(0,300+dt,dt)
x1=1.0*np.sin(2*np.pi*tt/0.1)
x2=0.5*np.sin(2*np.pi*tt/1.0)
x3=1.0*np.sin(2*np.pi*tt/10.0)
x4=2.0*np.sin(2*np.pi*tt/100.0)
x=x1+x2+x3+x4
xx1=x
nn,sp1=nfft_n(xx1)
spa1=np.sqrt(sp1.real[0:int(nn/2)+1]**2+sp1.imag[0:int(nn/2)+1]**2)*dt*nn
fk1=np.arange(0,nn/2+1)/nn/dt
sp=sp1
for i in range(0,len(fk1)):
if 0.02<=fk1[i]:
sp.real[i]=0; sp.real[nn-i]=0
sp.imag[i]=0; sp.imag[nn-i]=0
xx2=ifft_n(nn,sp)
nn,sp2=nfft_n(xx2)
spa2=np.sqrt(sp2.real[0:int(nn/2)+1]**2+sp2.imag[0:int(nn/2)+1]**2)*dt*nn
fk2=np.arange(0,nn/2+1)/nn/dt
drawfig1(tt,xx1,xx2) # time history
drawfig2(fk1,spa1,fk2,spa2) # Fourier spectrum
#---------------
# Execute
#---------------
if __name__ == '__main__': main()
以 上