Pythonで高速フーリエ変換(FFT)の練習-2 信号を時間軸と周波数軸で表現する

Pythonで高速フーリエ変換(FFT)の練習-2 信号を時間軸と周波数軸で表現する

Pythonで高速フリーエ変換(FFT)を行う方法をモモノキ&ナノネと一緒に学習していきます。

モモノキ&ナノネと一緒に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としよう。あとは前回と同じ要領でグラフまで作成してみて。

In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
In [2]:
# 簡単な信号の作成
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)
Out[2]:
[<matplotlib.lines.Line2D at 0xad33106c>]

うまくできたかも。

OK、x軸が時間になったね。次はそのデータでFFTを実行してみよう。

In [3]:
F = np.fft.fft(f) # 高速フーリエ変換(FFT)
In [4]:
# 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、周波数軸に変更する必要あり
Out[4]:
[<matplotlib.lines.Line2D at 0xad24b60c>]

FFTを実行して前回みたいにグラフまで描いてみたけど、周波数の検出点がちょっとズレてる感じがする。

うん、そのままではダメでx軸は周波数(Hz)に変更する必要があるよ。周波数(Hz)の取る範囲は、1/サンプリング周期(秒)で求まるよ。この場合は1/0.01だから100Hzまで。

つまり解析できる周波数は、サンプリング周期(周波数)に依存することになるね。今回のサンプリング周期は0.01秒に設定しているから、サンプリング周波数としては100Hzになるよ。

でも、前回の話で少し出てきた『なんとか定理』があるから、解析できる周波数はその半分まで?

うん、そうなんだ。標本化定理(サンプリング定理)と呼ばれるよ。解析できる周波数は最大でもサンプリング周波数の半分まで。その周波数の値をナイキスト定数と呼ぶから覚えておいてね。

OK、半分の値がナイキスト定数ね。

周波数軸の値を作るには、Numpyのlinspace(開始,終了,分割数)が便利だから使ってみて。

In [5]:
# 周波数軸のデータ作成
fq = np.linspace(0, 1.0/dt, N) # 周波数軸 linspace(開始,終了,分割数)

周波数軸のデータを作った。グラフにプロットしてみる。

In [6]:
# 周波数軸に変更してグラフを再表示
plt.xlabel('freqency(Hz)', fontsize=14)
plt.ylabel('signal amplitude', fontsize=14)
plt.plot(fq, F_abs_amp)
Out[6]:
[<matplotlib.lines.Line2D at 0xad17d1ec>]

周波数軸に変えたら、ピークも10のところにキチンと出た!

前回と一緒で後半の鏡像ピークが紛らわしいから、グラフの表示範囲はナイキスト定数までにしておく。

In [7]:
# グラフ表示
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]) # ナイキスト定数まで表示
Out[7]:
[<matplotlib.lines.Line2D at 0xad14b70c>]

うまくいったね。今回の練習はここまでにして、まとめたコードを下に載せておくね。

信号をもっと複雑にしたりするのは、前回練習通りだから一度試してみてね。

In [8]:
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]) # ナイキスト定数まで表示
Out[8]:
[<matplotlib.lines.Line2D at 0xad12556c>]

次回は逆FFTをやってみるよ。FFT(周波数軸)の結果を元の信号(時間軸)に戻す練習だよ。





スポンサーリンク

0 件のコメント :

コメントを投稿