解析エンジニアの自動化 blog

コツコツと自動化した方法を残す blog

Python で逆高速フーリエ変換(IFFT)をする



こんにちは。
仕事の自動化にやりがいと達成感を感じるガッくんです。


この記事の目次



背景・目的


フーリエ変換をやったので、逆フーリエ変換をやります。



動作環境


Windows 7
・winpython 64bit 3.4.4



プログラム

ソースコード


###############################################################################
# 逆高速フーリエ変換( IFFT )を計算するプログラム
###############################################################################
# インポート
import matplotlib.pyplot as plt
import numpy as np
 
# 簡単な信号の作成
N     = 2**6 # サンプル数
dt    = 0.025 # サンプリング周期( sec )
freq1 = 10 # 周波数( Hz )
ampl1 = 1 # 振幅
freq2 = 15 # 周波数( Hz )
ampl2 = 1 # 振幅
print("■■■ 入力条件 ■■■")
print("サンプル数 : " + str(N))
print("サンプリング周期 : " + str(dt) + ' sec')
print("サンプリング周波数 : " + str(1/dt) + ' Hz')
print("入力波1")
print("  周波数 : " + str(freq1) + ' Hz')
print("  振 幅 : " + str(ampl1))
print("入力波2")
print("  周波数 : " + str(freq2) + ' Hz')
print("  振 幅 : " + str(ampl2))
 
# 時間軸
t = np.arange(0, N*dt, dt)
print("サンプリング時間 : " + str(N*dt-dt) + ' sec')
print("サンプリング時刻 : " + str(t))
 
# {周波数 freq1, 振幅 amp1 の正弦入力波} + {周波数 freq2, 振幅 amp2 の正弦入力波}
f = ampl1*np.sin(2*np.pi*freq1*t) + ampl2*np.sin(2*np.pi*freq2*t)
# 周波数 freq1, 振幅 amp1 の正弦入力波
#f = amp1 * np.sin(2*np.pi*freq1*t)
 
# 高速フーリエ変換( FFT )
F = np.fft.fft(f)
 
# FFT の複素数結果を絶対に変換
absf = np.abs(F)
 
# 振幅をもとの信号に揃える
absf_amp = absf / N * 2
absf_amp[0] = absf_amp[0] / 2
 
# 周波数軸のデータ作成
fq = np.linspace(0, 1.0/dt, N) # 周波数軸 linspace(開始, 終了, 分割数)
 
idx = np.argmax(f)
print("\n■■■ 入力信号特性 ■■■")
print("入力信号の最大idx : " + str(idx))
print("入力信号の最大振幅 : " + str(f[idx]))
print("入力信号の最大時刻 : " + str(t[idx]))
 
idx = np.array(absf_amp[:int(N/2)+1]) # コピー
idx = idx.argsort()[::-1] # 降順に並べ替えた時のインデックスを取得する
print("\n■■■ フーリエ変換信号特性 ■■■")
print("フーリエ変換信号の最大idx : " + str(idx[0]))
print("フーリエ変換信号の最大振幅 : " + str(absf_amp[idx[0]]))
print("フーリエ変換信号の最大周波数 : " + str(fq[idx[0]]))
 
print("\nフーリエ変換信号の次点idx : " + str(idx[1]))
print("フーリエ変換信号の次点振幅 : " + str(absf_amp[idx[1]]))
print("フーリエ変換信号の次点周波数 : " + str(fq[idx[1]]))
 
# 逆フーリエ変換( IFFT )
F_ifft = np.fft.ifft(F)
 
# 実数部の取得
F_ifft_real = F_ifft.real
 
# スライス
F_ifft_real[:10]
 
idx = np.argmax(F_ifft_real)
print("\n■■■ 逆フーリエ変換信号特性 ■■■")
print("逆フーリエ変換信号の最大idx : " + str(idx))
print("逆フーリエ変換信号の最大振幅 : " + str(F_ifft_real[idx]))
print("逆フーリエ変換信号の最大時刻 : " + str(t[idx]))
 
# グラフ表示
fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(10, 12))
 
# 信号のグラフ(時間軸)
axes[0, 0].plot(t, f)
axes[0, 0].set_title('Input Wave')
axes[0, 0].set_xlabel('time[sec]')
axes[0, 0].set_ylabel('amplitude')
axes[0, 0].grid(True)
 
# FFTのグラフ(周波数軸)
axes[0, 1].plot(fq[:int(N/2)+1], absf_amp[:int(N/2)+1])
axes[0, 1].set_title('Fast Fourier Transform')
axes[0, 1].set_xlabel('freqency[Hz]')
axes[0, 1].set_ylabel('amplitude')
axes[0, 1].grid(True)
 
# IFFTのグラフ(時間軸)
axes[1, 0].plot(t, F_ifft_real, c="g")
axes[1, 0].set_title('Inverse Fast Fourier Transform')
axes[1, 0].set_xlabel('time[sec]')
axes[1, 0].set_ylabel('amplitude')
axes[1, 0].grid(True)
 
# Input Wave と IFFT Wave の重ね合わせ
axes[1, 1].plot(t, f, c="g")
axes[1, 1].plot(t, F_ifft_real, c="b")
axes[1, 1].set_title('Input Wave and IFFT')
axes[1, 1].set_xlabel('time[sec]')
axes[1, 1].set_ylabel('amplitude')
axes[1, 1].grid(True)
 
# グラフの出力
file_dir  = 'C:\\Users\\UserName\\Desktop\\'
file_name = 'ifft'
fig.savefig(file_dir + file_name + '0.0.jpg', bbox_unches="tight")



結果

図1 の 4 つのグラフの説明をします。

Input Wave は入力波のグラフです。

Fast Fourier Transform は Input Wave を高速フーリエ変換したグラフです。

Inverse Fast Fourier Transform は Fast Fourier Transform を逆高速フーリエ変換したグラフです。

Input Wave and IFFT は Input Wave と Inverse Fast Fourier Transform のグラフを重ね合わせたグラフです。

図1 変換結果


Input Wave and IFFT の図から Fast Fourier Transform の波形を逆高速フーリエ変換して求めた波形が入力した Input Wave とピッタリ重なっています。

そのため、正しく逆高速フーリエ変換が出来た事がわかります。

Python の 出力結果を示します。

In [1]: runfile('C:/WinPython-64bit-3.4.4.6Qt5/settings/.spyder-py3/fft-0.0.py', wdir='C:/WinPython-64bit-3.4.4.6Qt5/settings/.spyder-py3')
■■■ 入力条件 ■■■
サンプル数 : 64
サンプリング周期 : 0.025 sec
サンプリング周波数 : 40.0 Hz
入力波1
  周波数 : 10 Hz
  振 幅 : 1
入力波2
  周波数 : 15 Hz
  振 幅 : 1
サンプリング時間 : 1.5750000000000002 sec
サンプリング時刻 : [ 0. 0.025 0.05 ..., 1.525 1.55 1.575]
 
■■■ 入力信号特性 ■■■
入力信号の最大idx : 57
入力信号の最大振幅 : 1.70710678119
入力信号の最大時刻 : 1.425
 
■■■ フーリエ変換信号特性 ■■■
フーリエ変換信号の最大idx : 24
フーリエ変換信号の最大振幅 : 1.0
フーリエ変換信号の最大周波数 : 15.2380952381
 
フーリエ変換信号の次点idx : 16
フーリエ変換信号の次点振幅 : 1.0
フーリエ変換信号の次点周波数 : 10.1587301587
 
■■■ 逆フーリエ変換信号特性 ■■■
逆フーリエ変換信号の最大idx : 57
逆フーリエ変換信号の最大振幅 : 1.70710678119
逆フーリエ変換信号の最大時刻 : 1.425



コメント

Input Wave を FFT した後、ノイズを除去して、 IFFT する事があります。

Input Wave は横軸が時間のグラフで FFT したグラフは横軸が周波数です。

横軸が周波数になるとノイズなどが見やすくなることが多いので、波形をキレイにするときは、『入力波形 → FFT → ノイズ除去 → IFFT → キレイな波形』という手法は王道です。



以上