前回の『Pythonで初めてプログラムを書いてみた』でSIN曲線を作ったので今回は周波数分析の計算をしてみます。
仕事などでFFTするときはエクセルでデータを入手することが多いのでそのまま、エクセルのFFTマクロで計算していた。理由は出力画面を簡単に作成できるからです。
Pythonもライブラリを使うと簡単に出力グラフが作画できるので、今回はCSVファイルでデータを入手する前提でPythonのプログラムを使ってFFTを試してみました。
CSVデータとして、「sin」を16倍サンプリング、10周期の160個のデータを取り込んで、FFTでは最初の128個のデータで計算をしています。
使うライブラリを最初に宣言しておきます。FFTは色々なライブラリにありますが、numpyの中にあるFFTを利用します。
import csv
import matplotlib.pyplot as plt
import math
import numpy as np
別に関数定義をする必要はないと思いますが、関数として呼び出しています。最初にCSVファイルを読み込みますが、CSVファイルは呼び込むと文字列として扱われてしまうので数値変換をします。
def main():
file = './wavedata.csv' #ファイルのパスを指定
with open(file,'r',newline='') as f: #ファイルをオープン
reader = csv.reader(f) #ファイルオブジェクトをcsv.readerオブジェクトに変換
for row in reader:
y = [float(s) for s in row] #文字→数値変換
f.close()
なくても問題ないですが、入力データ数をFFT処理できるようにデータ数内で最大の2のn乗に自動で制限します。
N = 2**int(math.log2(len(row))) #入力を2のn乗に制限
XX = range(0,N,1) #X軸データ作成
YY = y[0:N] #Y軸データ作成
ライブラリを利用してFFT演算と、出力の振幅を計算しています。実際に周波数分析に必要なデータはN個ではなくて、N/2個なので半分の配列だけ取り出しています。「N/2」演算すると自動的に「float」になってしまうので型変換しています。
FF = np.fft.fft(YY) #FFTを計算(複素数出力)
absFF = np.abs(FF) #絶対値を計算
X0 = XX[0:int(N/2)]
absF0 = absFF[0:int(N/2)]
以下の部分でFFTの出力をグラフにプロットしています。
plt.clf()
plt.subplot(2,1,1)
plt.plot(XX,absFF,label = "FFT(all)")
plt.legend(loc = "upper right")
plt.subplot(2,1,2)
plt.plot(X0,absF0,label = "FFT(Frequency analysis)")
plt.legend(loc = "upper right")
plt.show()
if __name__ == '__main__':
main()
プログラムを動作させた結果として、図の上はすべてのデータ、下の図は半分のデータを表しています。FFTの振幅は真ん中で左右対称になるので周波数分析を行うには左半分だけで十分です。
FFTは複素入力ですが、一般的に分析するときには実数入力がほとんどなのでnumpyでは「np.fft.fft」以外に「np.fft.rfft」という実数専用のFFTもあるみたいです。若干、計算結果は異なるようですが、実用上は問題ないとのこと．．．
np.fft.fft input(64) → output(64) → 32
np.fft.rfft input(32) → output(32) → 32
処理時間は半分にはならないみたいですが高速化されるようです。
ここではFFTの使い方をまとめましたが、実際に使うときの実用的なノウハウや方法は『数学嫌いなエンジニアが楽しむFFT』のページに書きました。興味のある方は、ぜひクリックして見てください。
