モモノキ＆ナノネと一緒にPythonでFFTの使い方を覚えよう（2）
信号を時間軸と周波数軸でグラフに表現してみよう。
前回の続きでPythonで高速フリー変換（FFT）の練習をするよ。
ナノネ、周波数は知ってる？
たしか波の数だったと思うけど、高周波、低周波とかあるやつだよね。
うん、周波数は1秒間当たりの波の数で単位はヘルツ（Hz）だよ。ヘルツはよく耳にするよね。
電気だと地域によって50Hz、60Hzとかある。
周波数Hzは1秒間に繰り返される波の数だから、 50Hzは波（山と谷で1組）が1秒間に50回繰り返すことになるね。
そこまでは大丈夫。
前回FFTの練習用に作った簡易信号は時間を特に意識していなかったけど、今回はキチンと時間を考慮した信号を作ってFFTを実践してみよう。FFTは信号の時間軸を周波数軸に変換する意味合いがあるんだよ。
時間をキチンと考慮した信号を作るには、どうすればいいの？
データ間隔の時間を定義すればOKだよ。データの間隔はサンプリング周期ともいうよ。
データ数を128個（2の6乗）、サンプリング周期は0.01秒（10ms）としておこう。データ全体の時間は1.28秒（0.01×128）になるね。グラフで使う時間軸はNumpyのarrange(開始,終了,間隔)で簡単に作れるよ。
それと、周波数は10Hz（波数が1秒間に10個)に設定して、振幅は1としよう。あとは前回と同じ要領でグラフまで作成してみて。
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
# 簡単な信号の作成
N = 128 # サンプル数
dt = 0.01 # サンプリング周期(sec):100ms =>サンプリング周波数100Hz
freq = 10 # 周波数(10Hz) =>正弦波の周期0.1sec
amp = 1 # 振幅
t = np.arange(0, N*dt, dt) # 時間軸
f = amp * np.sin(2*np.pi*freq*t) # 信号（周波数10、振幅1の正弦波）
# グラフ表示
plt.xlabel('time(sec)', fontsize=14)
plt.ylabel('signal amplitude', fontsize=14)
plt.plot(t, f)
うまくできたかも。
OK、x軸が時間になったね。次はそのデータでFFTを実行してみよう。
F = np.fft.fft(f) # 高速フーリエ変換(FFT)
# 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) # ->NG、周波数軸に変更する必要あり
FFTを実行して前回みたいにグラフまで描いてみたけど、周波数の検出点がちょっとズレてる感じがする。
うん、そのままではダメでx軸は周波数（Hz）に変更する必要があるよ。周波数（Hz）の取る範囲は、1/サンプリング周期（秒）で求まるよ。この場合は1/0.01だから100Hzまで。
つまり解析できる周波数は、サンプリング周期（周波数）に依存することになるね。今回のサンプリング周期は0.01秒に設定しているから、サンプリング周波数としては100Hzになるよ。
でも、前回の話で少し出てきた『なんとか定理』があるから、解析できる周波数はその半分まで？
うん、そうなんだ。標本化定理（サンプリング定理）と呼ばれるよ。解析できる周波数は最大でもサンプリング周波数の半分まで。その周波数の値をナイキスト定数と呼ぶから覚えておいてね。
OK、半分の値がナイキスト定数ね。
周波数軸の値を作るには、Numpyのlinspace（開始,終了,分割数）が便利だから使ってみて。
# 周波数軸のデータ作成
fq = np.linspace(0, 1.0/dt, N) # 周波数軸 linspace(開始,終了,分割数)
周波数軸のデータを作った。グラフにプロットしてみる。
# 周波数軸に変更してグラフを再表示
plt.xlabel('freqency(Hz)', fontsize=14)
plt.ylabel('signal amplitude', fontsize=14)
plt.plot(fq, F_abs_amp)
周波数軸に変えたら、ピークも10のところにキチンと出た！
前回と一緒で後半の鏡像ピークが紛らわしいから、グラフの表示範囲はナイキスト定数までにしておく。
# グラフ表示
plt.xlabel('freqency(Hz)', fontsize=14)
plt.ylabel('signal amplitude', fontsize=14)
plt.plot(fq[:int(N/2)+1], F_abs_amp[:int(N/2)+1]) # ナイキスト定数まで表示
うまくいったね。今回の練習はここまでにして、まとめたコードを下に載せておくね。
信号をもっと複雑にしたりするのは、前回練習通りだから一度試してみてね。
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
# 簡単な信号の作成
N = 128 # サンプル数
dt = 0.01 # サンプリング周期(sec):100ms =>サンプリング周波数100Hz
freq = 10 # 周波数(10Hz) =>正弦波の周期0.1sec
amp = 1 # 振幅
t = np.arange(0, N*dt, dt) # 時間軸
f = amp * np.sin(2*np.pi*freq*t) # 信号（周波数10、振幅1の正弦波）
# 高速フーリエ変換(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倍不要
# 周波数軸のデータ作成
fq = np.linspace(0, 1.0/dt, N) # 周波数軸 linspace(開始,終了,分割数)
# グラフ表示
fig = plt.figure(figsize=(12, 4))
# 信号のグラフ（時間軸）
ax2 = fig.add_subplot(121)
plt.xlabel('time(sec)', fontsize=14)
plt.ylabel('amplitude', fontsize=14)
plt.plot(t, f)
# FFTのグラフ（周波数軸）
ax2 = fig.add_subplot(122)
plt.xlabel('freqency(Hz)', fontsize=14)
plt.ylabel('amplitude', fontsize=14)
plt.plot(fq[:int(N/2)+1], F_abs_amp[:int(N/2)+1]) # ナイキスト定数まで表示
次回は逆FFTをやってみるよ。FFT（周波数軸）の結果を元の信号（時間軸）に戻す練習だよ。
0 件のコメント :
コメントを投稿