PyLearn SIG 03: 周波数領域の特徴量¶
前回、横軸を時間、縦軸を信号値とするデータに対し、平均や分散などの時間領域特徴量を算出する方法について学習しました。今回は、時間領域ではなく信号の周波数領域に対して、その特徴量を算出する方法を学びます。このため、フーリエ変換と呼ばれる、時間領域の信号を周波数領域に変換する方法について学びます。その後、周波数領域特徴量として、パワーバンド、周波数領域エントロピー、エネルギーを学習します。
import numpy as np
import math as mt
import matplotlib.pyplot as plt
# *** 波形の構築 ***
f = 20 # 周波数
sf = 200 # サンプリング周波数
L = 50 # 信号長
t = np.arange(L) * (1/sf) # 時間
y = np.sin(2 * mt.pi * f * t) # sin波形取得
# *** 観測点の表示 ***
print("データ観測点:")
print(t)
# *** 波形の表示 ***
plt.plot(t, y) # 縦軸-横軸
plt.xlabel("time [s]")
plt.ylabel("value")
plt.show()
周波数が20Hzですから、1秒間に20回振動する波ということになります。周波数の逆数をとれば、一回の振動に要する時間が1/20=0.05[秒]と算出できます。確かに上の波形を見ると、0.05秒に一回振動している波であることがわかりますね。ここでわかるように、今眺めている信号は、はたしてどのくらいの周波数か？、という視点で考えることができます。
さて、このsin波はあらかじめこちらで式を立てたので、20Hzの成分が強いとすぐにわかります。ただ、センサなどで時系列信号を集めるとき、その信号はこちらが与えら数式から発生したものではないため、その周波数は未知数です。このようなとき、フーリエ変換を行うと、その信号はどのくらいの周波数か、調べることができます。実際に、今の信号を用いてフーリエ変換を行ってみます。
from scipy.fftpack import fft
# フーリエ変換
freq = np.linspace(0, sf, L) # 周波数領域の横軸
yf = np.abs(fft(y)) # 振幅スペクトル
yf = yf * yf # パワースペクトル
# 周波数領域の描画
plt.plot(freq[0:int(L/2)], yf[0:int(L/2)])
plt.xlabel("Hz")
plt.ylabel("Power")
plt.show()
横軸が周波数、縦軸がその成分の強さとする図が描画されました。このようなデータを周波数領域、グラフの名前をパワースペクトルと言います。横軸が20のところだけ値が高くなりました。つまり、20Hzの成分が強い信号だということを表しています。このように、FFTを行うことで、いくらの周波数帯が強い信号か、把握することが可能となります。サンプリング周波数の1/2の周波数をナイキスト周波数と呼びます。ナイキスト周波数以上の周波数帯はFFTをしてもうまく調べられないので、注意してください（これを標本化定理と呼びます。）。具体例を挙げると、
- サンプリング周波数が50Hzでは、ナイキスト周波数が25Hzなので、20Hzの信号をFFTすると解釈可能
- サンプリング周波数が50Hzでは、ナイキスト周波数が25Hzなので、30Hzの信号をFFTしても解釈不可能
- サンプリング周波数が50Hzでは、ナイキスト周波数が25Hzなので、25Hzの信号をFFTしても解釈不可能
あなたが慣性センサのサンプリング周波数を100Hzにしたならば、0Hzから50Hz未満の振動しかFFTで解析することができなくなるので、注意が必要です。
余談ですが、フーリエ変換は通常、高速フーリエ変換と呼ばれる方法で算出します。英語で書くと、Fast Fourier Translationとなりますので、略してFFTと呼ぶことが普通です。エンジニアの人たちがFFT、FFTということがよくありますので、それを聞いたときは、「あ、時系列信号にフーリエ変換をかけて、どういう周波数の信号か解析しようとしているんだな」と考えるようにして下さい。
2つのsin関数の合成の場合¶
次に、sin関数を足し算で結合させた以下の波形について考えてみます。
$y = \sin{2\pi f_1 t} + \sin{2\pi f_2 t}$
$f_1=20$[Hz], $f_2=40$[Hz]
# *** 波形の構築 ***
f1 = 20 # 周波数
f2 = 40
# *** 信号の情報 ***
sf = 200 # サンプリング周波数
L = 50 # 信号長
# 波形の算出
t = np.arange(L) * (1/sf)
y2 = np.sin(2 * mt.pi * f1 * t) + np.sin(2 * mt.pi * f2 * t)
# *** 観測点の表示 ***
print("データ観測点:")
print(t)
# *** 波形の表示 ***
plt.plot(t, y2) # 縦軸-横軸
plt.xlabel("time [s]")
plt.ylabel("value")
plt.show()
1秒間に20回振動する成分と、1秒間に40回振動する成分が混ざった信号なので、先ほどよりもちょっと汚いかもしれません。これをフーリエ変換してみましょう。
# フーリエ変換
freq = np.linspace(0, sf, L) # 周波数領域の横軸
yf = np.abs(fft(y2)) # 振幅スペクトル
yf = yf * yf # パワースペクトル
# 周波数領域の描画
plt.plot(freq[0:int(L/2)], yf[0:int(L/2)])
plt.xlabel("Hz")
plt.ylabel("Power")
plt.show()
このように、FFTを行うことで、20Hzと40Hzの成分が強い信号であることがわかります。
重みがある場合¶
以下の信号を見てください。sin関数の前に重み$a_1, a_2$がついています。20Hzの成分は強くて、40Hzの成分が小さいことを見て取れます。
$y = a_1 \sin{2\pi f_1 t} + a_2 \sin{2\pi f_2 t}$
$f_1=20$[Hz], $f_2=40$[Hz]
$a_1=300, a_2=50$
# *** 波形の構築 ***
f1 = 20 # 周波数
f2 = 40
a1 = 300 # 重み
a2 = 50
# *** 信号の情報 ***
sf = 200 # サンプリング周波数
L = 50 # 信号長
# 波形の算出
t = np.arange(L) * (1/sf)
y3 = a1 * np.sin(2 * mt.pi * f1 * t) + a2 * np.sin(2 * mt.pi * f2 * t)
# *** 観測点の表示 ***
print("データ観測点:")
print(t)
# *** 波形の表示 ***
plt.plot(t, y3) # 縦軸-横軸
plt.xlabel("time [s]")
plt.ylabel("value")
plt.show()
人間の目で見ると、20Hzの信号にしか見えません。40Hzの振動は弱すぎるので、波形を見てもわからないです。ここでFFTをしてみます。
# フーリエ変換
freq = np.linspace(0, sf, L) # 周波数領域の横軸
yf = np.abs(fft(y3)) # 振幅スペクトル
yf = yf * yf # パワースペクトル
# 周波数領域の描画
plt.plot(freq[0:int(L/2)], yf[0:int(L/2)])
plt.xlabel("Hz")
plt.ylabel("Power")
plt.show()
20Hzの振動にしか見えなかった波形ですが、FFTをすると、40Hzの成分もあるよ、ということがわかりました。このように、目ではわからない微小な周波数も把握することができます。
ノイズがある場合¶
さて、現実的な時系列信号は色々なノイズを含みます。ノイズがあるときのFFTも見てみましょう。対象とする信号は、以下の通りです。
$y = a_1 \sin{2\pi f_1 t} + a_2 \sin{2\pi f_2 t}$
$f_1=20$[Hz], $f_2=100$[Hz]
$a_1=300, a_2=50$
一つ上とまったく同じ信号です。ですので、20Hzの成分が強く、40Hzの成分が弱い信号です。
# *** 波形の構築 ***
f1 = 20 # 周波数
f2 = 40
a1 = 300 # 重み
a2 = 100
# *** 信号の情報 ***
sf = 200 # サンプリング周波数
L = 50 # 信号長
# 波形の算出
t = np.arange(L) * (1/sf)
y4 = a1 * np.sin(2 * mt.pi * f1 * t) + a2 * np.sin(2 * mt.pi * f2 * t)
# ノイズ生成 & 付与
np.random.seed(1) # 再現可能な乱数
noise = np.random.normal(0, 100, L) # ノイズ生成（平均0, 標準偏差100の正規分布に従う乱数）
y4 = y4 + noise # 綺麗な波形にノイズを付与
# *** 観測点の表示 ***
print("データ観測点:")
print(t)
# *** 波形の表示 ***
plt.plot(t, y4) # 縦軸-横軸
plt.xlabel("time [s]")
plt.ylabel("value")
plt.show()
ノイズ付きの信号です。きれいで人工的なsin波形と比べ、少しだけ現実っぽくなってきました。この信号に対し、FFTをかけてみます。
# フーリエ変換
freq = np.linspace(0, sf, L) # 周波数領域の横軸
yf = np.abs(fft(y4)) # 振幅スペクトル
yf = yf * yf # パワースペクトル
# 周波数領域の描画
plt.plot(freq[0:int(L/2)], yf[0:int(L/2)])
plt.xlabel("Hz")
plt.ylabel("Power")
plt.show()
FFTの図も歪んできました。ノイズがあるとこのような結果になります。ただし、20Hzの成分があるということは一目瞭然ですし、40Hzの成分があることも見て取れますね。なお、乱数を大きくすると、40Hzの成分は観測不可能になりそうですね。実際にそうなるか、自分で試してみてください。正規乱数の標準偏差を大きくすれば、乱数が強くなります。
フーリエ変換のまとめ¶
- 周波数が$n$ Hzである信号とは、1秒間に$n$回振動する波という意味である。
- サンプリング周波数$f$ Hzはデータを1秒回に$f$回記録するという意味で、異なるので注意すること。
- フーリエ変換とは、時系列信号がどのような周波数の成分を有するか調べる分析手法である。（1秒回に何回振動する成分が強いか？）
- 複数の周波数成分が混ざっていても分析可能である。
- フーリエ変換を略すとFFTと呼ぶ。
- サンプリング周波数を$f$ Hzとした場合、FFTで解析可能な周波数は$\frac{f}{2}$ Hz未満の範囲である。これを標本化定理と呼ぶ。
- $\frac{f}{2}$ Hz をナイキスト周波数と呼ぶ。
備考: pythonのfft関数に時系列信号を入力して得られる結果を、フーリエスペクトルと言います。これは複素数表記により位相情報を含む値ですので、情報量は多いのですが、ちょっとわかりにくいです。そのため、絶対値をつけるか二乗を行うことで、複素数情報を消失されます。フーリエスペクトルに絶対値をつけたものを振幅スペクトル、振幅スペクトルを二乗したものをパワースペクトルと呼びます。今回行なったプログラムでは、fft関数をabs関数により絶対値に変換し、2乗しているため、パワースペクトルを求めていることになります。グラフを書くときは、パワースペクトルの場合は、Powerを縦軸に、振幅スペクトルの場合はAmplitudeを縦軸に記載することが一般的です。パワースペクトルも振幅スペクトルも、物理的な解釈の仕方はほとんど一緒です。
センサデータのフーリエ変換¶
さて、これまで、「40Hzのsin波でFFTをしたら、40Hzの成分が強いことがわかった」みたいなことをやってきました。これまではあくまで基礎の学習でしたので、40Hzと答えがわかっている信号に、40Hzだと当ててきたわけです。しかし、実際の分析で使用する慣性センサデータは、当然のことながら、いくつの周波数の成分が強いかわかりません。そのため、FFTで調べてあげるわけです。
まず、慣性センサデータをロードし、合成加速後を計算します。
# 信号のロード
import numpy as np
signal = np.loadtxt("dat/1.csv", delimiter=",", skiprows=0, usecols=range(2,8))
# 合成加速度の算出
M = signal[:,0] * signal[:,0] + signal[:,1] * signal[:,1] + (signal[:,2]) * (signal[:,2])
M = np.sqrt(M)
M_new = np.reshape(M, (len(M), 1))
# 7チャンネル慣性センサデータの取得
sig7 = np.concatenate([signal, M_new], 1)
このセンサデータはサンプリング周波数100Hzで取得しています。これを踏まえてX軸加速度に対してFFTをしてみます。
# *** 信号情報取得 ***
L = len(sig7[:, 0]) # 信号長
sf = 100 # サンプリング周波数
# フーリエ変換
freq = np.linspace(0, sf, L) # 周波数領域の横軸
yf = np.abs(fft(sig7[:, 0])) # 振幅スペクトル
yf = yf * yf # パワースペクトル
fig = plt.figure(figsize=(20, 8)) # 描画領域の横、縦のサイズ
# 周波数領域の描画
plt.subplot(2,1,1)
plt.plot(freq[0:int(L/2)], yf[0:int(L/2)])
plt.xlabel("Hz", size = 20)
plt.ylabel("Power", size = 20)
plt.title("Frequency Domain of Xacc", size = 20)
# 時間領域の描画
plt.subplot(2,1,2)
plt.plot(sig7[:, 0])
plt.xlabel("time: 0.01[sec]", size=20)
plt.ylabel("Accelecration: 0.1[mG]", size=20)
plt.title("Time Domain of Xacc", size = 20)
plt.tight_layout() # グラフが被るのを防止
plt.show()
上が周波数領域、下が時間領域です（どちらもX軸加速度の信号）。人工的に作った信号と違い、かなり複雑な感じになりました。自然な信号はなかなか複雑ということです。傾向を見ると、
- 0Hz
- 16〜19Hz付近
- 26〜32Hz付近
- 38〜41Hz付近
- 48〜50Hz付近
の成分が強いことがわかりますね。ここで、0Hzを信号の直流成分といい、信号の振動ではなく、高さを表します。1秒間に0回振動する波を思い浮かべてください。それは横軸に平行な直線になることは想像がつくでしょうか。ですので、FFTなどで周波数解析を行うとき、0Hz成分の強さに目を向けることはあまりありません。
これがわかったところで、0Hzを無視して16〜19Hzに目を向けてみます。この高さを取得してみます。これを知るには、そもそも配列の何番目に16〜19Hzの成分の強さが入っているんだ？ということを知らねばなりません。このために、周波数領域の横軸であるfreqの値をfor文で一つずつprintしてみましょう。
for i in range(0, int(L/2)): # 信号長の半分だけ繰り返し
print("freqの", i, "番目は", np.round(freq[i], 2), "Hzです。")
最小0Hz、最大が50Hz未満です。今回、サンプリング周波数が100Hzですから、ナイキスト周波数は50Hzになります。したがって、50Hz以上の成分をみることに、意味はありません。なぜ信号長Lの半分だけfor文を回すのか疑問に思った方もいるかもしれません。これは、信号長の半分がナイキスト周波数ギリギリのラインになるためです。これはちょっと直感的にわかりにくいので、ゆっくり考えた方がいいです。ひとまずそれを置いておくと、だいたい73番目から86番目あたりにありそうです。実際に、この範囲の周波数の強さをfor文で確認してみましょう。
for i in range(73, 89):
print(np.round(freq[i], 2), "Hzの強さは", np.round(yf[i], 2), "です。")
この個別の数字と周波数領域のグラフを見比べてみると、だいたい同じ値であることがわかります。また、17.7Hzの強さが27億を超えています（他は10億を超えていません）。すなわち、X軸加速度の信号は、1秒間に17.7回振動する波の成分が最も強いということがわかりました。
周波数領域の特徴量: パワーバンド¶
周波数領域のことがわかったところで、ここから特徴量を取り出す方法について説明します。この表し方には色々な方法がありますが、特徴量は人間が解釈しやすいような物理的意味を伴った方があれこれやりやすいので、例えば以下のようにしてみます。
- 0〜5 Hzのパワーの総和（ただし、直流成分は除く）
- 5〜10 Hzのパワーの総和
...
- 45〜50 Hzのパワーの総和（ただし、ナイキスト周波数は除く）
この特徴量を、パワーバンドと呼びます。この名称は、以下の論文で使用されているものを参考にしました。
Morris, D., Saponas, T. S., Guillory, A., & Kelner, I. (2014, April). RecoFit: using a wearable sensor to find, recognize, and count repetitive exercises. In Proceedings of the SIGCHI Conference on Human Factors in Computing Systems (pp. 3225-3234). ACM.
パワーバンドとして今回、刻み幅を5Hzごとにしましたが、10Hzごとなどでも良いと思います。もし、機械の故障を検出したいとして、機械が故障すると、10Hzの成分が弱まり、20Hzの成分が強まるのではないか、という仮説があるのであれば、それらが互いに独立するような刻み幅にすべきでしょう。そういった仮説などがない場合は、取り敢えずてきとうに決めましょう。
以下に、時系列信号を入力したら、FFTを自動ですると、パワーバンドを算出してくれる関数を示します。
# パワーバンドを求める関数
# signal: 時系列信号（1軸分のみ）
# sf: サンプリング周波数
# min_f: パワーバンドを算出する周波数の最小値
# max_f: パワーバンドを算出する周波数の最小値
def Calc_PowerBand(signal, sf, min_f, max_f):
import numpy as np # <- どこかでインポートしていればなくて良い。
# パワースペクトル算出
L = len(signal) # 信号長
freq = np.linspace(0, sf, L) # 周波数
yf = np.abs(fft(signal)) # 振幅スペクトル
yf = yf * yf # パワースペクトル
# 指定された周波数の探索
n1, n2 = [], []
for i in range(1, int(L/2)): # <- 直流成分は探索範囲から除外 & ナイキスト周波数まで
n1.append(np.abs(freq[i] - min_f))
n2.append(np.abs(freq[i] - max_f))
min_id = np.argmin(n1) # min_fと最も近い周波数が格納された配列番号取得
max_id = np.argmin(n2)# max_fと最も近い周波数が格納された配列番号取得
# 指定された周波数のパワーバンド算出
PB = np.sum(yf[min_id:max_id])
return PB
関数ができたところで、実際にX軸加速度のパワーバンドを5Hz刻みで計算してみます。
FV_PB = []
axis_name = []
sf = 100 # サンプリング周波数
# 5Hz刻みでパワーバンドを計算
for i in range(0, 10):
FV_PB.append(Calc_PowerBand(sig7[:,0], sf, i*5, (i+1)*5))
axis_name.append(str(i*5) + "-" + str((i+1)*5) + "Hz")
# printで表示
print(FV_PB)
# 棒グラフで表示
plt.bar(axis_name, FV_PB)
plt.xticks(rotation=90) # 横軸ラベル回転
plt.title("Power Bands of Xacc Signal")
plt.show()
このように、X軸加速度の信号を、周波数領域に変換し、そこから10個の特徴量に変換することができました。先ほど描画したX軸加速度のパワースペクトルとこの結果を見比べてみてください。結果が一致していることがわかりますね。これにより、X軸加速度の信号に対して、10個の特徴量を算出することができました。同じように、先ほど定義した関数で、Y軸やZ軸の加速度、角速度、合成加速度で同じように特徴量を算出することが可能です。これにより、1回の観測データに対し、合計70個の周波数領域特徴量を算出することができます。
周波数領域の特徴量: 周波数領域エントロピー¶
もう一つ、周波数領域エントロピーと呼ばれる特徴量も存在します。エントロピーとは、状態の複雑さを表す指標です。エントロピーが0に近いとき状態は極めて単純なこと、値が大きいとき状態が極めて複雑であることを表します。周波数領域エントロピーとは、周波数領域、すなわち、パワースペクトルの複雑さを表します。パワースペクトルとは「1回に何回振動するか」に関する総合的な情報を持つデータです。これが「単純である」とは、30Hzただ一つの周波数により構成された信号を表します。一方、「複雑である」とは、20Hz、24Hz、32Hz、41Hz…、というように、複数の周波数の成分が強く混ざり合った信号を意味します。この複雑さを表す特徴量として、周波数領域エントロピーが生まれました。周波数領域は以下の計算により算出できます。
- $i$番目のフーリエスペクトル$f_i$に対して、そのパワーの総和を算出する。なお、信号長は$L$とし、直流成分($j=1$)は計算対象から除くものとする。総和を$L/2$までとしている理由は、$L/2$がナイキスト周波数ギリギリのラインであるためである（これ以上見ても意味がない）。:
$ a = \sum^{L/2}_{j=2} |f_j|^2$
- $i$番目のパワースペクトル$|f_i|^2$を$a$によって規格化する。これを規格化パワースペクトル$p(i)$とする:
$p(i) = \frac{|f_i|^2}{a}$
- $i$番目の規格化されたパワースペクトル$p(i)$に対して、直流成分を除くすべての$i$に対し、そのエントロピーを算出する。これを周波数領域エントロピー$E$とする。
$E = -\sum^{L/2}_{j=2} p(i) \log p(i)$
以上、3つの手続きで周波数領域エントロピーが算出できます。規格化して、エントロピーで複雑さ・乱雑さを測るイメージです。よくわからない人は、この特徴量が大きいとき・小さいとき、それぞれどういう意味かわかればひとまずはokかと思います。
以下に、周波数領域エントロピーを算出する関数を記載します。
# 周波数領域エントロピーを求める関数
# signal: 時系列信号（1軸分のみ）
# sf: サンプリング周波数
def Calc_FreqEnt(signal, sf):
import numpy as np # <- どこかでインポートしていればなくて良い。
# パワースペクトル算出
L = len(signal) # 信号長
freq = np.linspace(0, sf, L) # 周波数
yf = np.abs(fft(signal)) # 振幅スペクトル
yf = yf * yf # パワースペクトル
# 手続き1
a = 0
for i in range(1, int(L/2)): # <- 直流成分を抜く & ナイキスト周波数まで
a = a + yf[i]
# 手続き2 & 3
E = 0
for i in range(1, int(L/2)): # <- 直流成分を抜く & ナイキスト周波数まで
E = E + (yf[i]/a) * np.log2(yf[i]/a)
E = -E
return E
関数ができたところで、以下の3つのsin波形に対して、周波数領域エントロピーを算出してみます。
- 20Hz, 40Hz, 60Hzの合成により得られるsin波形（複雑さ: 大）
- 20Hz, 40Hzの合成により得られるsin波形（複雑さ: 中）
- 20Hzのsin波形（複雑さ: 小）
# サンプリング周波数
sf = 200
# 3つの周波数の場合
f1 = 20
f2 = 40
f3 = 60
y = np.sin(2 * mt.pi * f1 * t) + np.sin(2 * mt.pi * f2 * t) + np.sin(2 * mt.pi * f3 * t)
print("20Hz, 40Hz, 60Hzの波のエントロピー: ", round(Calc_FreqEnt(y, sf), 3))
# 2つの周波数の場合
f1 = 20
f2 = 40
y = np.sin(2 * mt.pi * f1 * t) + np.sin(2 * mt.pi * f2 * t)
print("20Hz, 40Hzの波のエントロピー: ", round(Calc_FreqEnt(y, sf), 3))
# 1つの周波数の場合
f1 = 20
y = np.sin(2 * mt.pi * f1 * t)
print("20Hzの波のエントロピー: ", round(Calc_FreqEnt(y, sf), 3))
実際に出力を確認してみると、複雑な信号が一番エントロピーが高く、単純な信号であるほどエントロピーが低いことがわかります。このように、周波数領域エントロピーにより、時系列信号の複雑さを周波数領域の面から把握することが可能になります。
次に、慣性センサデータ、7チャンネル分（X軸加速度、Y軸加速度、...、合成加速度）の周波数領域エントロピーを算出してみます。
# 周波数領域エントロピーの算出
sf = 100 # サンプリング周波数
FV_FE = []
for i in range(0, len(sig7[0, :])):
FV_FE.append(Calc_FreqEnt(sig7[:,i], sf))
# printで表示
print(FV_FE)
# 棒グラフで表示
axis_name_FE = ["Xacc", "Yacc", "Zacc", "Xang", "Yang", "Zang", "Mag"]
plt.bar(axis_name_FE, FV_FE)
plt.xticks(rotation=90) # 横軸ラベル回転
plt.title("Frequency Domain Entropy")
plt.ylim((6, 7.2)) # y軸のレンジ
plt.show()
ほとんど同じ値になりました。あえて言えば、Y軸角速度の周波数領域エントロピーが高いようです。今回ロードしている慣性センサデータは、ある一度の観測を、7チャンネルという側面で見ているので、ほとんど同じというのもなんとなく頷けます。実際に比較すべきなのは、故障している機械の振動と、故障していない機械の振動のように、異なる観測に対してとなります。
周波数領域の特徴量: エネルギー¶
次に、エネルギーと呼ばれる特徴量について解説します。エネルギーは、直流成分を除いたパワースペクトルの平均値により定義されます。すなわち、
$E = \frac{1}{L/2}\sum_{i=2}^{L/2}|f_i|^2$,
により算出できます。$f_i$は$i$番目のフーリエスペクトルです。エネルギーが高いほど、信号の強度が強いことを示しています。以下に、エネルギーを算出する関数を定義します。
def Calc_Energy(signal, sf):
import numpy as np # <- どこかでインポートしていればなくて良い。
# パワースペクトル算出
L = len(signal) # 信号長
freq = np.linspace(0, sf, L) # 周波数
yf = np.abs(fft(signal)) # 振幅スペクトル
yf = yf * yf # パワースペクトル
# エネルギー算出
En = np.mean(yf[1:int(L/2)]) # 直流成分を除く & ナイキスト周波数まで
return En
次に、7軸分のエネルギーを算出してみます。
# 周波数領域エントロピーの算出
sf = 100 # サンプリング周波数
FV_En = []
for i in range(0, len(sig7[0, :])):
FV_En.append(Calc_Energy(sig7[:,i], sf))
# printで表示
print(FV_En)
# 棒グラフで表示
axis_name_En = ["Xacc", "Yacc", "Zacc", "Xang", "Yang", "Zang", "Mag"]
plt.bar(axis_name_En, FV_En)
plt.xticks(rotation=90) # 横軸ラベル回転
plt.title("Energy")
plt.show()
Z軸加速度とマグニチュードのみ、すごく大きな値が出ました。直流成分は抜いているので、重力加速度の強さを表しているわけではないので、注意してください。実際に、Z軸加速度の周波数領域を見てみましょう。
# *** 信号情報取得 ***
L = len(sig7[:, 2]) # 信号長
sf = 100 # サンプリング周波数
# フーリエ変換
freq = np.linspace(0, sf, L) # 周波数領域の横軸
yf = np.abs(fft(sig7[:, 2])) # <- Z軸加速度を対象とする
yf = yf * yf # パワースペクトル
fig = plt.figure(figsize=(20, 8)) # 描画領域の横、縦のサイズ
# 周波数領域の描画
plt.subplot(2,1,1)
plt.plot(freq[1:int(L/2)], yf[1:int(L/2)]) # <- 直流成分は重力加速度が入るので見ない
plt.xlabel("Hz", size = 20)
plt.ylabel("Power", size = 20)
plt.title("Frequency Domain of Zacc", size = 20)
# 時間領域の描画
plt.subplot(2,1,2)
plt.plot(sig7[:, 0])
plt.xlabel("time: 0.01[sec]", size=20)
plt.ylabel("Accelecration: 0.1[mG]", size=20)
plt.title("Time Domain of Zacc", size = 20)
plt.tight_layout() # グラフが被るのを防止
plt.show()
このように、だいたい30Hz以降に強い成分があることがわかります。これが、エネルギーが大きい理由の一つでしょう。
まとめ¶
今回、周波数領域特徴量として、以下を取り扱いました。
- パワーバンド
- 周波数領域エントロピー
- エネルギー
これにより、信号の波に関する成分に対して、様々な特徴を抽出することができます。時間領域特徴量などではみられない特徴となりますので、検討すると良いと思います。なお、フーリエ変換を行う関係上、時間領域特徴量よりも計算時間は大きいです。軟弱なCPUしか扱えない場合などは、計算方法を工夫したり、そもそも周波数領域特徴量は使わないといった選択を行うことも必要になります。