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

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

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



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


この記事の目次



背景・目的


とにかく何か分析めいたものがやりたくて、特に理由はないけれどフーリエ変換をやりました。

波形の分析の初歩なので、まとめておいて損はないと思いまして…



動作環境


Windows 7
・winpython 64bit 3.4.4



プログラム

ソースコード


###############################################################################
# 高速フーリエ変換( FFT )を計算するプログラム
###############################################################################
# インポート
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)
 
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]]))
 
# グラフ表示
fig, (axL, axR) = plt.subplots(ncols=2, figsize=(10,4))
 
# 信号のグラフ(時間軸)
axL.plot(t, f)
axL.set_title('Input Wave')
axL.set_xlabel('time[sec]')
axL.set_ylabel('amplitude')
axL.grid(True)
 
# FFTのグラフ(周波数軸)
axR.plot(fq[:int(N/2)+1], absf_amp[:int(N/2)+1])
axR.set_title('Fast ')
axR.set_xlabel('freqency[Hz]')
axR.set_ylabel('amplitude')
axR.grid(True)
 
# グラフの出力
file_dir  = 'C:\\Users\\UserName\\Desktop\\'
file_name = 'fft'
fig.savefig(file_dir + file_name + '0.0.jpg', bbox_unches="tight")



結果

numpy には高速フーリエ変換のメソッドがあります。
numpy は便利ですね。

図1 の Input Wave は分析した波形で、周波数 10Hz, 振幅 1 の正弦波と周波数 15Hz, 振幅 1 の正弦波の合成波でした。

図1 の Fast Fourier Transform は分析結果で、周波数が 10Hz と 15Hz でところで波形がとんがっています。
さらに、振幅がそれぞれ 1 を示しています。

見事に入力波が周波数 10Hz, 振幅 1 の正弦波と周波数 15Hz, 振幅 1 の正弦波の合成波である事が分かりました。

図1 変換結果



In [1]: runfile('C:/WinPython-64bit-3.4.4.6Qt5/settings/.spyder-py3/fft0.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



コメント

入力波形がどんな単純波形の合成で出来ているのかが分かる様になりました。

自動車の振動の測定結果や観測した地震波がどんな単純波形の合成で出来ているのかがわかる様になるため、対策すれば良い単純な波形がわかる様になります。

ただし、振動の測定結果や観測地震波などは波形が複雑なため、構成される単純波形がものすごい数になる事でしょう…



以上