モモノキ＆ナノネと一緒にPythonでFFTの使い方を覚えよう（1）
簡単な信号を作って高速フーリエ変換（FFT）に挑戦してみよう
今回はPythonで高速フリー変換（FFT）を使ってみよう。
モモノキ、高速フーリエ変換（FFT）で何ができるの？
簡単に言うと周期を持った信号の解析ができるよ。周期的な信号は正弦波の足し合わせで全て表現できるんだけど、FFTで解析すると、どんなの信号が組み合わせで構成されているかが一発で分かるんだ。
解析で各信号成分の周期や振幅などが分かると、特定の波形だけ抽出したり取り除いたりすることもできる。信号の加工処理にも応用されるよ。
高速フーリエ変換はFFTと呼ばれて、Fast Fourier Transformの略称。フリーエ変換にもいくつか方式があって、FFTは離散フーリエ変換を使っているよ。FFTは離散フーリエ変換を高速に処理できるアルゴリズムなんだ。もっと詳しくは書籍などで調べてみてね。
フーリエ変換ってなんか難しそうだけど、FFTは簡単に使えるの？
フーリエ変換は複素数を扱うから計算式の理解はとても難しいけど、PythonでFFTを使って信号を解析するだけなら比較的簡単にできるよ。
難しい式が分からなくても信号解析ができるならFFTを試してみる。
ナノネ、サインカーブ（正弦波）のグラフ描ける？範囲は適当でもいいから。
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from scipy import signal
x = np.arange(0, 2*np.pi, 0.1)
y = np.sin(x)
plt.plot(x, y)
サインカーブのグラフは、なんとか描ける。こんなんでいい？
# グラフのx軸をπで表示
x = np.arange(0, 2*np.pi, 0.1)
y = np.sin(x)
plt.xticks(np.arange(0, 2.1*np.pi, 0.5*np.pi),
['0','1/2$\pi$',u'$\pi$','3/2$\pi$',u'2$\pi$'], fontsize=14) # πtex表記
plt.plot(x, y)
追加でx軸のラベルをπ表記に変えてみた。もうちょっとグラフをオシャレにする？
オシャレにしなくても大丈夫。そのグラフは使わないから。ただの練習。
なんか損した感じ。
次にFFTに使う簡単な信号データを作ってみよう。
FFTで使うデータ数は、2のn乗個（nは整数）にする必要があるから注意ね。FFTはデータ数を制限することで計算を簡略化して高速処理を可能にしているよ。
まずは少ないデータ数で簡単な信号を作って、FFTを一度体験してみよう。（周波数とかはまた別途ね）
データ数はとりあず32個（2の5乗）でいいよ。さっきサインカーブを描いた要領で、32個のデータの中に正弦波を1サイクル入れた信号を作ってみて。
N = 32 # データ数
n = np.arange(32)
signal = np.sin(2*np.pi*n/N)
# グラフ表示
plt.figure(figsize=(8, 4))
plt.xlabel('n')
plt.ylabel('Signal')
plt.plot(signal)
できたよ。32個のデータの中に1サイクル入れたよ。
FFTを実行すれば周期が確認できるんだけど、1サイクルだと面白くないよね。周期を少し増やして3サイクル分をデータの中に入れてみよう。
信号の周期を変えたい場合は、sin()の中を周期で定数倍すればOKだよね。
N = 32 # データ数
n = np.arange(N)
freq = 3 # 周期
f = np.sin(freq * 2 * np.pi * (n/N)) # freq=3周期分
# グラフ表示
plt.figure(figsize=(8, 4))
plt.xlabel('n')
plt.ylabel('Signal')
plt.plot(f)
できたけど、データ32個しか使ってないから波形がカクカクだけどいいの？
大丈夫、FFTの練習には十分だよ。
今回、高速フリーエ変換はNumpyに実装されているfft.fft()を使うよ。Numpyは最初にnpという名前でインポート済みだから、np.fft.fft(データ)とすればOK。実行するとFFT処理結果が返ってくるから変数で受け取ってね。
F = np.fft.fft(f) # 高速フーリエ変換(FFT)
実行したけど、これでFFT完了？
うん、完了。FFTの処理結果はnumpy.ndarray形式で値は複素数で取得されるよ。
print(type(F), F.dtype)
print(F) # FFT結果
周期はどうややって確認するの？複素数を32個列挙されても分からないよ。
周期を確認するにはFFTで取得した複素数を絶対値に変換する必要があるよ。複素数を絶対値に変換する方法は、np.abs(FFT結果)とすればOK。変換した値でグラフにプロットしてみて。信号に入れた周期(x軸の値)のところでピークが現れるはずだよ。
F_abs = np.abs(F)
plt.plot(F_abs)
3くらいのところにピークが出た！信号に周期3のサイン波をいれたからFFTの結果は合ってそうだね。でも、右側30弱のところにもピークが出現してるよ？こんな大きな周期は入れてないはずだけど。。。
周期を確認するとき、後半側に出現するピークは無視してね。難しい話になるけど複素共役に関係して鏡像とういう形で出現するピークなんだ。
あと大事なポイントがあって、周期はデータ数の半分までしか確認できなから注意ね。周波数についてはまたあとで取り扱う予定だけど、詳しくは標本化定理（サンプリング定理）というのも調べてみてね。
周期を確認できるのは、データ数の半分までね。メモっておく。
plt.plot(F_abs[:int(N/2)+1])
後半側に出現する鏡像ピークが紛らわしいから、グラフ表示範囲は半分までに変更した。
周期を確認するだけなら、このままでも十分なんだけど。今、y軸の値は意味のない数値になっているから、元の信号の振幅に合うように値を変換してみよう。
元の振幅に揃えるには、ピーク強度の値をデータ数で割ればOKだよ。ただし、対となる鏡像の値も足しわせる必要があるのでデータ数で割ったあと2倍してあげてね。（データ数の半分で割っても同じだよ）
F_abs_amp = F_abs / N * 2 # 交流成分はデータ数で割って2倍する
F_abs_amp[0] = F_abs_amp[0] / 2 # 直流成分（今回は扱わないけど）は2倍不要
plt.plot(F_abs_amp[:int(N/2)+1])
周期3のままで、振幅は1に変わった。元の信号の振幅と一緒。
信号の振幅を変更して、FFTで正しい値が得られるかも確認してみよう。振幅を変更するにはsinの値を定数倍すればOKだね。
N = 32 # データ数
n = np.arange(N)
freq = 3 # 周期
amp = 4 # 振幅
f = amp * np.sin(freq * 2 * np.pi * (n/N)) # 周期freq=3、振幅amp=4
# グラフ表示
plt.figure(figsize=(8, 4))
plt.xlabel('n')
plt.ylabel('Signal')
plt.plot(f)
周期は3のままで、振幅を4に変えて信号を作成。
# 高速フーリエ変換(FFT)
F = np.fft.fft(f)
# FFT結果（複素数）を絶対値に変換
F_abs = np.abs(F)
# 振幅を元に信号に揃える
F_abs_amp = F_abs / N * 2 # 交流成分はデータ数で割って2倍する
F_abs_amp[0] = F_abs_amp[0] / 2 # 直流成分（今回は扱わないけど）は2倍不要
# グラフ表示（データ数の半分の周期を表示）
plt.plot(F_abs_amp[:int(N/2)+1])
FFTの結果を確認。振幅4と正しく変換できた！
その調子。次はもう少し信号を複雑にしてみよう。異なる2つ周期を足し合わせた信号を作ってでFFTを実行してみて。信号は周期2と6のサインカーブを足し合わせた波形で試してみよう。sin関数の値を単純に足すだけでOKだよ。
N = 32 # データ数
n = np.arange(N)
f1 = 2 # 周期
f2 = 6
f = np.sin(f1 * 2 * np.pi * (n/N)) + np.sin(f2 * 2 * np.pi * (n/N))
# グラフ表示
plt.figure(figsize=(8, 4))
plt.xlabel('n')
plt.ylabel('Signal')
plt.plot(f)
周期2と6のサインカーブを足した信号を作った。振幅はどっちも1。作った信号をプロット、んっ。さっきよりさらにカクカク感が増してバルタン星人みないな波形だけど大丈夫？
大丈夫だと思うけど、ちょっとカッコ悪いから少しデータ数増やして滑らかにしておこう。データ数64個（2の6乗）で。
N = 64 # データ数
n = np.arange(N)
f1 = 2 # 周期
f2 = 6
f = np.sin(f1 * 2 * np.pi * (n/N)) + np.sin(f2 * 2 * np.pi * (n/N))
# グラフ表示
plt.figure(figsize=(8, 4))
plt.xlabel('n')
plt.ylabel('Signal')
plt.plot(f)
さっきよりは丸まった感じ。FFTを実行してみる。
# 高速フーリエ変換(FFT)
F = np.fft.fft(f)
# FFT結果（複素数）を絶対値に変換
F_abs = np.abs(F)
# 振幅を元に信号に揃える
F_abs_amp = F_abs / N * 2 # 交流成分はデータ数で割って2倍する
F_abs_amp[0] = F_abs_amp[0] / 2 # 直流成分（今回は扱わないけど）は2倍不要
# グラフ表示（データ数の半分の周期を表示）
plt.plot(F_abs_amp[:int(N/2)+1])
成功、信号に入れた周期2と6のところにちゃんとピークがでた。振幅も1でOK。
少し慣れたから、周期と振幅を変えて試してみる。
N = 64 # データ数
n = np.arange(N)
f1 = 4 # 周期①
f2 = 10 # 周期②
a1 = 1.5 # 振幅①
a2 = 3 # 振幅②
f = a1 * np.sin(f1 * 2 * np.pi * (n/N)) + a2 * np.sin(f2 * 2 * np.pi * (n/N))
# グラフ表示
plt.figure(figsize=(8, 4))
plt.xlabel('n')
plt.ylabel('Signal')
plt.plot(f)
振幅1.5で周期4と振幅3で周期3の正弦波を足した信号を作ってみた。FFTでうまく分析でるかな？
# 高速フーリエ変換(FFT)
F = np.fft.fft(f)
# FFT結果（複素数）を絶対値に変換
F_abs = np.abs(F)
# 振幅を元に信号に揃える
F_abs_amp = F_abs / N * 2 # 交流成分はデータ数で割って2倍する
F_abs_amp[0] = F_abs_amp[0] / 2 # 直流成分（今回は扱わないけど）は2倍不要
# グラフ表示（データ数の半分の周期を表示）
plt.plot(F_abs_amp[:int(N/2)+1])
振幅と周期どっちも正解！FFTけっこう賢いね。
OKだね。今回はここまで、続きはまたにしよう。次回は周波数にキチンと変換してFFTを練習するよ。
In16で周期10のはずですがセリフが3になってますよ
返信削除