PyLearn SIG 03: 周波数領域の特徴量

前回、横軸を時間、縦軸を信号値とするデータに対し、平均や分散などの時間領域特徴量を算出する方法について学習しました。今回は、時間領域ではなく信号の周波数領域に対して、その特徴量を算出する方法を学びます。このため、フーリエ変換と呼ばれる、時間領域の信号を周波数領域に変換する方法について学びます。その後、周波数領域特徴量として、パワーバンド、周波数領域エントロピー、エネルギーを学習します。

フーリエ変換

フーリエ変換とは、時系列信号の周波数領域を算出する方法の一つです。詳細を説明するとややこしいので、具体例から設営していきます。

1つのsin波形の場合

まず、一番シンプルな事例として、以下のsin波形について考えてみます。

$y = \sin{2\pi f t}$

$f=20$[Hz]

$t$は時間で、$f$は周波数です。今回は20Hzですから、1秒間に20回振動するsin波形であることを示しています。

In [2]:
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()
データ観測点:
[0.    0.005 0.01  0.015 0.02  0.025 0.03  0.035 0.04  0.045 0.05  0.055
 0.06  0.065 0.07  0.075 0.08  0.085 0.09  0.095 0.1   0.105 0.11  0.115
 0.12  0.125 0.13  0.135 0.14  0.145 0.15  0.155 0.16  0.165 0.17  0.175
 0.18  0.185 0.19  0.195 0.2   0.205 0.21  0.215 0.22  0.225 0.23  0.235
 0.24  0.245]

周波数が20Hzですから、1秒間に20回振動する波ということになります。周波数の逆数をとれば、一回の振動に要する時間が1/20=0.05[秒]と算出できます。確かに上の波形を見ると、0.05秒に一回振動している波であることがわかりますね。ここでわかるように、今眺めている信号は、はたしてどのくらいの周波数か?、という視点で考えることができます。

さて、このsin波はあらかじめこちらで式を立てたので、20Hzの成分が強いとすぐにわかります。ただ、センサなどで時系列信号を集めるとき、その信号はこちらが与えら数式から発生したものではないため、その周波数は未知数です。このようなとき、フーリエ変換を行うと、その信号はどのくらいの周波数か、調べることができます。実際に、今の信号を用いてフーリエ変換を行ってみます。

In [3]:
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]

In [4]:
# *** 波形の構築 ***
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()
データ観測点:
[0.    0.005 0.01  0.015 0.02  0.025 0.03  0.035 0.04  0.045 0.05  0.055
 0.06  0.065 0.07  0.075 0.08  0.085 0.09  0.095 0.1   0.105 0.11  0.115
 0.12  0.125 0.13  0.135 0.14  0.145 0.15  0.155 0.16  0.165 0.17  0.175
 0.18  0.185 0.19  0.195 0.2   0.205 0.21  0.215 0.22  0.225 0.23  0.235
 0.24  0.245]

1秒間に20回振動する成分と、1秒間に40回振動する成分が混ざった信号なので、先ほどよりもちょっと汚いかもしれません。これをフーリエ変換してみましょう。

In [5]:
# フーリエ変換
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$

In [6]:
# *** 波形の構築 ***
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()
データ観測点:
[0.    0.005 0.01  0.015 0.02  0.025 0.03  0.035 0.04  0.045 0.05  0.055
 0.06  0.065 0.07  0.075 0.08  0.085 0.09  0.095 0.1   0.105 0.11  0.115
 0.12  0.125 0.13  0.135 0.14  0.145 0.15  0.155 0.16  0.165 0.17  0.175
 0.18  0.185 0.19  0.195 0.2   0.205 0.21  0.215 0.22  0.225 0.23  0.235
 0.24  0.245]

人間の目で見ると、20Hzの信号にしか見えません。40Hzの振動は弱すぎるので、波形を見てもわからないです。ここでFFTをしてみます。

In [7]:
# フーリエ変換
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の成分が弱い信号です。

In [8]:
# *** 波形の構築 ***
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()
データ観測点:
[0.    0.005 0.01  0.015 0.02  0.025 0.03  0.035 0.04  0.045 0.05  0.055
 0.06  0.065 0.07  0.075 0.08  0.085 0.09  0.095 0.1   0.105 0.11  0.115
 0.12  0.125 0.13  0.135 0.14  0.145 0.15  0.155 0.16  0.165 0.17  0.175
 0.18  0.185 0.19  0.195 0.2   0.205 0.21  0.215 0.22  0.225 0.23  0.235
 0.24  0.245]

ノイズ付きの信号です。きれいで人工的なsin波形と比べ、少しだけ現実っぽくなってきました。この信号に対し、FFTをかけてみます。

In [9]:
# フーリエ変換
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で調べてあげるわけです。

まず、慣性センサデータをロードし、合成加速後を計算します。

In [10]:
# 信号のロード
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をしてみます。

In [11]:
# *** 信号情報取得 ***
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してみましょう。

In [12]:
for i in range(0, int(L/2)): # 信号長の半分だけ繰り返し
    print("freqの", i, "番目は", np.round(freq[i], 2), "Hzです。")
freqの 0 番目は 0.0 Hzです。
freqの 1 番目は 0.22 Hzです。
freqの 2 番目は 0.44 Hzです。
freqの 3 番目は 0.66 Hzです。
freqの 4 番目は 0.88 Hzです。
freqの 5 番目は 1.11 Hzです。
freqの 6 番目は 1.33 Hzです。
freqの 7 番目は 1.55 Hzです。
freqの 8 番目は 1.77 Hzです。
freqの 9 番目は 1.99 Hzです。
freqの 10 番目は 2.21 Hzです。
freqの 11 番目は 2.43 Hzです。
freqの 12 番目は 2.65 Hzです。
freqの 13 番目は 2.88 Hzです。
freqの 14 番目は 3.1 Hzです。
freqの 15 番目は 3.32 Hzです。
freqの 16 番目は 3.54 Hzです。
freqの 17 番目は 3.76 Hzです。
freqの 18 番目は 3.98 Hzです。
freqの 19 番目は 4.2 Hzです。
freqの 20 番目は 4.42 Hzです。
freqの 21 番目は 4.65 Hzです。
freqの 22 番目は 4.87 Hzです。
freqの 23 番目は 5.09 Hzです。
freqの 24 番目は 5.31 Hzです。
freqの 25 番目は 5.53 Hzです。
freqの 26 番目は 5.75 Hzです。
freqの 27 番目は 5.97 Hzです。
freqの 28 番目は 6.19 Hzです。
freqの 29 番目は 6.42 Hzです。
freqの 30 番目は 6.64 Hzです。
freqの 31 番目は 6.86 Hzです。
freqの 32 番目は 7.08 Hzです。
freqの 33 番目は 7.3 Hzです。
freqの 34 番目は 7.52 Hzです。
freqの 35 番目は 7.74 Hzです。
freqの 36 番目は 7.96 Hzです。
freqの 37 番目は 8.19 Hzです。
freqの 38 番目は 8.41 Hzです。
freqの 39 番目は 8.63 Hzです。
freqの 40 番目は 8.85 Hzです。
freqの 41 番目は 9.07 Hzです。
freqの 42 番目は 9.29 Hzです。
freqの 43 番目は 9.51 Hzです。
freqの 44 番目は 9.73 Hzです。
freqの 45 番目は 9.96 Hzです。
freqの 46 番目は 10.18 Hzです。
freqの 47 番目は 10.4 Hzです。
freqの 48 番目は 10.62 Hzです。
freqの 49 番目は 10.84 Hzです。
freqの 50 番目は 11.06 Hzです。
freqの 51 番目は 11.28 Hzです。
freqの 52 番目は 11.5 Hzです。
freqの 53 番目は 11.73 Hzです。
freqの 54 番目は 11.95 Hzです。
freqの 55 番目は 12.17 Hzです。
freqの 56 番目は 12.39 Hzです。
freqの 57 番目は 12.61 Hzです。
freqの 58 番目は 12.83 Hzです。
freqの 59 番目は 13.05 Hzです。
freqの 60 番目は 13.27 Hzです。
freqの 61 番目は 13.5 Hzです。
freqの 62 番目は 13.72 Hzです。
freqの 63 番目は 13.94 Hzです。
freqの 64 番目は 14.16 Hzです。
freqの 65 番目は 14.38 Hzです。
freqの 66 番目は 14.6 Hzです。
freqの 67 番目は 14.82 Hzです。
freqの 68 番目は 15.04 Hzです。
freqの 69 番目は 15.27 Hzです。
freqの 70 番目は 15.49 Hzです。
freqの 71 番目は 15.71 Hzです。
freqの 72 番目は 15.93 Hzです。
freqの 73 番目は 16.15 Hzです。
freqの 74 番目は 16.37 Hzです。
freqの 75 番目は 16.59 Hzです。
freqの 76 番目は 16.81 Hzです。
freqの 77 番目は 17.04 Hzです。
freqの 78 番目は 17.26 Hzです。
freqの 79 番目は 17.48 Hzです。
freqの 80 番目は 17.7 Hzです。
freqの 81 番目は 17.92 Hzです。
freqの 82 番目は 18.14 Hzです。
freqの 83 番目は 18.36 Hzです。
freqの 84 番目は 18.58 Hzです。
freqの 85 番目は 18.81 Hzです。
freqの 86 番目は 19.03 Hzです。
freqの 87 番目は 19.25 Hzです。
freqの 88 番目は 19.47 Hzです。
freqの 89 番目は 19.69 Hzです。
freqの 90 番目は 19.91 Hzです。
freqの 91 番目は 20.13 Hzです。
freqの 92 番目は 20.35 Hzです。
freqの 93 番目は 20.58 Hzです。
freqの 94 番目は 20.8 Hzです。
freqの 95 番目は 21.02 Hzです。
freqの 96 番目は 21.24 Hzです。
freqの 97 番目は 21.46 Hzです。
freqの 98 番目は 21.68 Hzです。
freqの 99 番目は 21.9 Hzです。
freqの 100 番目は 22.12 Hzです。
freqの 101 番目は 22.35 Hzです。
freqの 102 番目は 22.57 Hzです。
freqの 103 番目は 22.79 Hzです。
freqの 104 番目は 23.01 Hzです。
freqの 105 番目は 23.23 Hzです。
freqの 106 番目は 23.45 Hzです。
freqの 107 番目は 23.67 Hzです。
freqの 108 番目は 23.89 Hzです。
freqの 109 番目は 24.12 Hzです。
freqの 110 番目は 24.34 Hzです。
freqの 111 番目は 24.56 Hzです。
freqの 112 番目は 24.78 Hzです。
freqの 113 番目は 25.0 Hzです。
freqの 114 番目は 25.22 Hzです。
freqの 115 番目は 25.44 Hzです。
freqの 116 番目は 25.66 Hzです。
freqの 117 番目は 25.88 Hzです。
freqの 118 番目は 26.11 Hzです。
freqの 119 番目は 26.33 Hzです。
freqの 120 番目は 26.55 Hzです。
freqの 121 番目は 26.77 Hzです。
freqの 122 番目は 26.99 Hzです。
freqの 123 番目は 27.21 Hzです。
freqの 124 番目は 27.43 Hzです。
freqの 125 番目は 27.65 Hzです。
freqの 126 番目は 27.88 Hzです。
freqの 127 番目は 28.1 Hzです。
freqの 128 番目は 28.32 Hzです。
freqの 129 番目は 28.54 Hzです。
freqの 130 番目は 28.76 Hzです。
freqの 131 番目は 28.98 Hzです。
freqの 132 番目は 29.2 Hzです。
freqの 133 番目は 29.42 Hzです。
freqの 134 番目は 29.65 Hzです。
freqの 135 番目は 29.87 Hzです。
freqの 136 番目は 30.09 Hzです。
freqの 137 番目は 30.31 Hzです。
freqの 138 番目は 30.53 Hzです。
freqの 139 番目は 30.75 Hzです。
freqの 140 番目は 30.97 Hzです。
freqの 141 番目は 31.19 Hzです。
freqの 142 番目は 31.42 Hzです。
freqの 143 番目は 31.64 Hzです。
freqの 144 番目は 31.86 Hzです。
freqの 145 番目は 32.08 Hzです。
freqの 146 番目は 32.3 Hzです。
freqの 147 番目は 32.52 Hzです。
freqの 148 番目は 32.74 Hzです。
freqの 149 番目は 32.96 Hzです。
freqの 150 番目は 33.19 Hzです。
freqの 151 番目は 33.41 Hzです。
freqの 152 番目は 33.63 Hzです。
freqの 153 番目は 33.85 Hzです。
freqの 154 番目は 34.07 Hzです。
freqの 155 番目は 34.29 Hzです。
freqの 156 番目は 34.51 Hzです。
freqの 157 番目は 34.73 Hzです。
freqの 158 番目は 34.96 Hzです。
freqの 159 番目は 35.18 Hzです。
freqの 160 番目は 35.4 Hzです。
freqの 161 番目は 35.62 Hzです。
freqの 162 番目は 35.84 Hzです。
freqの 163 番目は 36.06 Hzです。
freqの 164 番目は 36.28 Hzです。
freqの 165 番目は 36.5 Hzです。
freqの 166 番目は 36.73 Hzです。
freqの 167 番目は 36.95 Hzです。
freqの 168 番目は 37.17 Hzです。
freqの 169 番目は 37.39 Hzです。
freqの 170 番目は 37.61 Hzです。
freqの 171 番目は 37.83 Hzです。
freqの 172 番目は 38.05 Hzです。
freqの 173 番目は 38.27 Hzです。
freqの 174 番目は 38.5 Hzです。
freqの 175 番目は 38.72 Hzです。
freqの 176 番目は 38.94 Hzです。
freqの 177 番目は 39.16 Hzです。
freqの 178 番目は 39.38 Hzです。
freqの 179 番目は 39.6 Hzです。
freqの 180 番目は 39.82 Hzです。
freqの 181 番目は 40.04 Hzです。
freqの 182 番目は 40.27 Hzです。
freqの 183 番目は 40.49 Hzです。
freqの 184 番目は 40.71 Hzです。
freqの 185 番目は 40.93 Hzです。
freqの 186 番目は 41.15 Hzです。
freqの 187 番目は 41.37 Hzです。
freqの 188 番目は 41.59 Hzです。
freqの 189 番目は 41.81 Hzです。
freqの 190 番目は 42.04 Hzです。
freqの 191 番目は 42.26 Hzです。
freqの 192 番目は 42.48 Hzです。
freqの 193 番目は 42.7 Hzです。
freqの 194 番目は 42.92 Hzです。
freqの 195 番目は 43.14 Hzです。
freqの 196 番目は 43.36 Hzです。
freqの 197 番目は 43.58 Hzです。
freqの 198 番目は 43.81 Hzです。
freqの 199 番目は 44.03 Hzです。
freqの 200 番目は 44.25 Hzです。
freqの 201 番目は 44.47 Hzです。
freqの 202 番目は 44.69 Hzです。
freqの 203 番目は 44.91 Hzです。
freqの 204 番目は 45.13 Hzです。
freqの 205 番目は 45.35 Hzです。
freqの 206 番目は 45.58 Hzです。
freqの 207 番目は 45.8 Hzです。
freqの 208 番目は 46.02 Hzです。
freqの 209 番目は 46.24 Hzです。
freqの 210 番目は 46.46 Hzです。
freqの 211 番目は 46.68 Hzです。
freqの 212 番目は 46.9 Hzです。
freqの 213 番目は 47.12 Hzです。
freqの 214 番目は 47.35 Hzです。
freqの 215 番目は 47.57 Hzです。
freqの 216 番目は 47.79 Hzです。
freqの 217 番目は 48.01 Hzです。
freqの 218 番目は 48.23 Hzです。
freqの 219 番目は 48.45 Hzです。
freqの 220 番目は 48.67 Hzです。
freqの 221 番目は 48.89 Hzです。
freqの 222 番目は 49.12 Hzです。
freqの 223 番目は 49.34 Hzです。
freqの 224 番目は 49.56 Hzです。
freqの 225 番目は 49.78 Hzです。

最小0Hz、最大が50Hz未満です。今回、サンプリング周波数が100Hzですから、ナイキスト周波数は50Hzになります。したがって、50Hz以上の成分をみることに、意味はありません。なぜ信号長Lの半分だけfor文を回すのか疑問に思った方もいるかもしれません。これは、信号長の半分がナイキスト周波数ギリギリのラインになるためです。これはちょっと直感的にわかりにくいので、ゆっくり考えた方がいいです。ひとまずそれを置いておくと、だいたい73番目から86番目あたりにありそうです。実際に、この範囲の周波数の強さをfor文で確認してみましょう。

In [13]:
for i in range(73, 89):
    print(np.round(freq[i], 2), "Hzの強さは", np.round(yf[i], 2), "です。")
16.15 Hzの強さは 676035460.78 です。
16.37 Hzの強さは 39554544.63 です。
16.59 Hzの強さは 1420458416.22 です。
16.81 Hzの強さは 402076768.52 です。
17.04 Hzの強さは 164021935.25 です。
17.26 Hzの強さは 1440221358.67 です。
17.48 Hzの強さは 490503011.98 です。
17.7 Hzの強さは 2716362388.43 です。
17.92 Hzの強さは 358067709.34 です。
18.14 Hzの強さは 410830147.88 です。
18.36 Hzの強さは 177688316.12 です。
18.58 Hzの強さは 303243362.05 です。
18.81 Hzの強さは 470365934.57 です。
19.03 Hzの強さは 280564671.95 です。
19.25 Hzの強さは 168230131.33 です。
19.47 Hzの強さは 88294717.07 です。

この個別の数字と周波数領域のグラフを見比べてみると、だいたい同じ値であることがわかります。また、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を自動ですると、パワーバンドを算出してくれる関数を示します。

In [14]:
# パワーバンドを求める関数
   # 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刻みで計算してみます。

In [15]:
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()
[3927731412.8775616, 2789937393.547026, 2124892031.238901, 13150698520.938198, 4461355815.556784, 6891695724.569159, 2830943498.840117, 5381923285.120466, 2285457944.565505, 3474538669.844453]

このように、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かと思います。

以下に、周波数領域エントロピーを算出する関数を記載します。

In [16]:
# 周波数領域エントロピーを求める関数
   # 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波形(複雑さ: 小)
In [17]:
# サンプリング周波数
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))
20Hz, 40Hz, 60Hzの波のエントロピー:  1.585
20Hz, 40Hzの波のエントロピー:  1.0
20Hzの波のエントロピー:  0.0

実際に出力を確認してみると、複雑な信号が一番エントロピーが高く、単純な信号であるほどエントロピーが低いことがわかります。このように、周波数領域エントロピーにより、時系列信号の複雑さを周波数領域の面から把握することが可能になります。

次に、慣性センサデータ、7チャンネル分(X軸加速度、Y軸加速度、...、合成加速度)の周波数領域エントロピーを算出してみます。

In [18]:
# 周波数領域エントロピーの算出
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()
[6.952776169695395, 6.8899938902225255, 6.694661909614187, 6.612753842650501, 7.027488561973456, 6.589424328620089, 6.699857283890942]

ほとんど同じ値になりました。あえて言えば、Y軸角速度の周波数領域エントロピーが高いようです。今回ロードしている慣性センサデータは、ある一度の観測を、7チャンネルという側面で見ているので、ほとんど同じというのもなんとなく頷けます。実際に比較すべきなのは、故障している機械の振動と、故障していない機械の振動のように、異なる観測に対してとなります。

周波数領域の特徴量: エネルギー

次に、エネルギーと呼ばれる特徴量について解説します。エネルギーは、直流成分を除いたパワースペクトルの平均値により定義されます。すなわち、

$E = \frac{1}{L/2}\sum_{i=2}^{L/2}|f_i|^2$,

により算出できます。$f_i$は$i$番目のフーリエスペクトルです。エネルギーが高いほど、信号の強度が強いことを示しています。以下に、エネルギーを算出する関数を定義します。

In [19]:
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軸分のエネルギーを算出してみます。

In [20]:
# 周波数領域エントロピーの算出
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()
[205676788.95128754, 88695212.33588338, 640176903.2983073, 51499786.36933906, 49704421.52600191, 1803797.4107932204, 647899630.1897469]

Z軸加速度とマグニチュードのみ、すごく大きな値が出ました。直流成分は抜いているので、重力加速度の強さを表しているわけではないので、注意してください。実際に、Z軸加速度の周波数領域を見てみましょう。

In [25]:
# *** 信号情報取得 ***
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しか扱えない場合などは、計算方法を工夫したり、そもそも周波数領域特徴量は使わないといった選択を行うことも必要になります。