Python: fft と ifft の使い方復習

思うところあって、Python の fft と ifft の使い方の復習をしておく。

# 基本的な使い方
fftとifft（高速フーリエ変換および逆変換）は、numpyとscipyの
両方に含まれている。このため、双方を使うプログラムを作ってみた。

```python
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
```

# 少し応用的な使い方
以下のように４つの周期を持つサイン波を重ね合わせた波を考える。

```math
\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秒のサイン波以外の成分を取り去って、フーリエ逆変換をかけてみる。


結果は以下の通り。

### 時刻歴の比較
![fig_wave1.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/129300/0fa353a4-d0ae-5a25-f5ea-f1dd09b7463f.jpeg)

### フーリエスペクトルの比較
![fig_wave2.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/129300/523db6d2-fc82-f7e8-8194-319760ef10e7.jpeg)


## プログラム
フーリエ変換と逆変換にはnumpyのライブラリを使っている。これは、iPadで用いるPythonista 3などではscipyは使えないため、移植することを考えると、できるだけscipyを使わないほうがいいため。本来ならできるだけscipyを使いたいところなのだがしかたがない。

フーリエ変換で求まった複素フーリエ係数を扱う場合、折れ曲がり点（ナイキスト振動数）の向こう側も対象となるよういじっておく必要があることに注意。

フーリエ逆変換では、複素フーリエ係数を入力して波形の時刻歴に戻すわけであるが、ifftの戻り値は複素数なので、プロット用には実数部のみを取り出す必要があることに注意。

上で紹介したフーリエスペクトルの図では```scipy.signal.find_peaks```により、ピーク周期をいれるようにした。プログラム例では、ピークの値が50を超えるものにその周期を記載するようにしている。図中、先性データの周期は100秒であるが、フーリwスペクトルのピークは109秒となっている。これはおそらく周期100秒の波の個数が少ないためではないかと思う。（未検証）


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

以　上