PyLearn SIG02: 時間領域の特徴量¶
通常、信号は多数の数字の集まりによって構成されるため、情報量が多いことから、うまく解析することができません。そのため、その信号が有する特徴を表す別の数字に変換する必要があります。抽出された特徴は通常、生の信号よりも数字の数が少ないので、低次元化と見做すことができます。このような多数の数字の集まりから特徴を抽出する数学的方法として、
- 平均や分散といった統計量を使用する方法、
- なんらかの目的関数に従い低次元化写像を構成する方法（例えば、深層学習、主成分・因子分析、特異値分解、オートエンコーダなど）
この2つがあります。前者と後者では色々な違いがありますが、抑えておきたいポイントは、人間が解釈可能か、不可能か、この2つです。前者は物理的意味を伴うので、人間が解釈可能です（例えば、分散が大きい信号は、ばらつきが大きい、など）。一方、後者は変数間の相関関係の除去するような特徴抽出、元のデータをできるだけ再現できるような特徴抽出といった仕組みが導入されており、そこに「人間にとってのわかりやすさ」は一切含まれていません。機械学習と信号処理の専門家と言い張りたい場合にはぜんぶ理解することが必要ですが、大変ですので、ここでは前者のみを扱うことにします。
人間がわかりやすい信号の特徴量として、これまた2つの種類に分けることができます。
- 時間領域の特徴量
- 周波数領域の特徴量
前者は、横軸が時間、縦軸が信号の大きさとした場合に、そこから抽出された特徴を表します（つまり、そのままの信号から特徴を取り出します）。後者は、信号に対しフーリエ変換やウェーブレット変換をかけ、周波数領域を得ます。周波数領域は、いくらの周波数がどのくらい強いかを表すデータです。この中から抽出された特徴を周波数領域特徴量と言います。
これを踏まえ、ここでは時間領域特徴量について解説します。時間領域特徴量は色々な研究者が色々あれこれ提案しているので網羅は到底できません。従って、ここでは代表的なものとして、平均、分散、歪度、尖度、最小値、最大値、中央値、RMS、動的特徴量などがあります。
前準備（7チャンネル慣性センサデータの取得）¶
はじめに、特徴量の算出対象である7チャンネル慣性センサデータを計算しておきます。
# 信号のロード
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)
平均、分散、歪度、尖度、最小値、最大値、中央値¶
まず、簡単に求められる特徴量から説明します。
- 平均値: 信号の平均的な値を意味します。
- 分散: 信号のばらつきを表します。
- 歪度: 信号をヒストグラムにしたとき、左右、どちらに裾が長いか。
- 尖度: 信号をヒストグラムにしたとき、どれだけ尖っているか。
- 最小値: 信号の一番小さな値
- 最大値: 信号の一番大きな値
- 中央値: 信号の中央の値
以上、合計7個の特徴量です。SIG01で7チャンネルの慣性センサデータがありましたので、$7 \times 7 = 49$個の特徴量に変換されるわけです。それでは、この7個の特徴量を算出してみます。
# 平均値
Xacc_Mean = np.mean(sig7[:, 0]) # X軸 加速度
Yacc_Mean = np.mean(sig7[:, 1]) # Y軸 加速度
Zacc_Mean = np.mean(sig7[:, 2]) # Z軸 加速度
Xang_Mean = np.mean(sig7[:, 3]) # X軸 角速度
Yang_Mean = np.mean(sig7[:, 4]) # Y軸 角速度
Zang_Mean = np.mean(sig7[:, 5]) # Z軸 角速度
Macc_Mean = np.mean(sig7[:, 6]) # 合成加速度
print("X軸角速度の平均値: ", Xang_Mean)
という感じで、7チャンネルの平均値を求めることができます。しかし、これを分散や歪度などずらずら並べるのはあまり頭がいい感じがしません。普通は以下のようにfor文を使ってやります。
import scipy.stats as st
# 空のリストを用意
FV_time= []
# 特徴量算出
for i in range(0, len(sig7[0,:])): # range(0, 7)という意味（7チャンネル分）
FV_time.append(np.mean(sig7[:,i])) # 平均
FV_time.append(np.var(sig7[:,i])) # 分散
FV_time.append(st.skew(sig7[:,i])) # 歪度
FV_time.append(st.kurtosis(sig7[:,i])) # 尖度
FV_time.append(np.min(sig7[:,i])) # 最小
FV_time.append(np.max(sig7[:,i])) # 最大
FV_time.append(np.median(sig7[:,i])) # 中央
# X軸加速度の平均、...、X軸加速度の中央値、Y軸加速度の平均、...という順番でリストへappend
# 表示
k=0
ChannelName = ["Xacc", "Yacc", "Zacc", "Xang", "Zang", "Zang", "Mag"]
FeatureName = ["mean", "var", "sk", "ku", "min", "max", "med"]
print("特徴量の数は", len(FV_time), "個です！")
for i in range(0, len(ChannelName)):
for j in range(0, len(FeatureName)):
print("FV", k, ":", ChannelName[i], "の", FeatureName[j], "は", round(FV_time[k], 3), "です。")
k=k+1
print(" *** this channel is finished.*** ")
RMS (Root Mean Square)¶
今、X軸加速度がプラスのとき右への加速、マイナスのとき左への加速とします。このとき、ある単位時間当たりに左右に同じくらい加速していれば、その平均値は0となりますね。従って、特徴量の一つである平均値は、符号により左右の加速を打ち消し合う場合があります。これを一つの特徴として表すことも当然ありです。ただ、このような打ち消し合いがされない特徴量を欲しい場合もあるでしょう。これを実現する特徴量に、二乗平均平方根があります。英語で書くとRoot Mean Squareとなり、工学の分野ではRMSとよく呼ばれます。例えば、ある時間$i$から$i+n$のX軸加速度のRMSは、
${\rm RMS}(X_{\rm acc}, i, i+n) = \frac{1}{n}\sqrt{\sum_{t = i}^{i+n} X_{\rm acc}(t)^2},$
で定義されます。他の軸も同様の式で算出可能です。ある時間$t$の信号を二乗することで、すべてプラス方向に信号が整理されます。その後、平方根を取ることで、二乗され見かけ上値が大きくなってしまった部分がキャンセルされます。最後に、その平均値を算出しています。つまり、RMSとは、左右の向きを気にせずに、左右両者の方向においてその全体の加速の大きさを表す特徴になります。上下方向の軸であれば、上下両方においてその全体の加速の大きさを表します。
RMSを求める関数¶
# RMSを求める関数
def Calc_RMS(signal):
a = signal * signal # 二乗
sum_a = np.sum(a) # 総和
sqrt_a = np.sqrt(sum_a) # 平方根
RMS = np.mean(sqrt_a) # 平均値
return RMS
関数を利用して7チャンネルのRMSを算出¶
# RMS 取得
FV_RMS=[]
for i in range(0, len(sig7[0,:])): # 7週分
FV_RMS.append(Calc_RMS(sig7[:,i]))
# 表示
ChannelName = ["Xacc", "Yacc", "Zacc", "Xang", "Zang", "Zang", "Mag"]
for i in range(0, len(sig7[0,:])): # 7週分
print(ChannelName[i], "のRMS: ", round(FV_RMS[i], 3))
動的特徴量（回帰係数）¶
動的特徴量とは、その信号は右肩上がりか、右肩下がりかを表現する特徴量です。ある時刻$i$からある時刻$i+n$の信号において、時間$t$からX軸加速度の値$X_{\rm acc}(t)$を、
$X_{\rm acc}(t)' = a t + b$
として推定する線形回帰モデルを考えます。ここで$a$は傾きを表すパラメータ、$b$は信号のバイアス成分を表すパラメータとなります。実際の信号$X_{\rm acc}(t)$と時間$t$から推定された値$X_{\rm acc}(t)'$が一致していれば良いので、$パラメータa, b$は実測値と推定値のずれ$X_{\rm acc}(t) - X_{\rm acc}(t)'$が最小となる値が採用されます。
この$a$が、動的特徴量となり、
- $a=0$: その信号は、右肩上がりでも右肩下がりでもない。
- $a>0$: その信号は、右肩上がりである。
- $a<0$: その信号は、右肩下がりである。
という意味を持ちます。$a$の絶対値が大きいほど、上がり/下がりが急激であることを意味しています。例えば、以下はX軸加速度の時刻50から時刻80の信号をプロットしたものです。
import matplotlib.pyplot as plt
plt.plot(sig7[50:80, 0])
plt.xlabel("t")
plt.ylabel("Acceleration")
plt.title("Xacc")
plt.show()
これを見ると、信号が右肩下がりであることがわかります。ということは、動的特徴量はマイナスの値となるはずです。本当にそうなるか、確認してみます。
from sklearn.linear_model import LinearRegression
# 入力変数（信号値に対応する連番の数字）
X=np.arange(80-50)
X=np.reshape(X, (len(X), 1))
# 出力変数（X軸加速度の50〜80番目の信号値）
Y=sig7[50:80, 0]
# 線形回帰モデルのパラメータフィッティング
ModelLR = LinearRegression() # 線形回帰モデルのセット
ModelLR.fit(X, Y) # パラメータ獲得（Xは時間、Yは信号値）
print('回帰係数（動的特徴量）a = ', round(ModelLR.coef_[0], 3)) # a
print('切片 b = ', round(ModelLR.intercept_, 3)) # b
動的特徴量は-18となりました（切片bは参考のため出しました。気にしなくてokです）。マイナスなので、動的特徴量は「右肩下がりの信号だ」と主張しています。これがわかったところで、動的特徴量を算出する関数を定義して、実際に7チャンネル分の特徴量を抽出してみます。
動的特徴量を算出する関数¶
# 動的特徴量を求める関数
def Calc_LinearCoef(signal):
# 入力変数（信号値に対応する連番の数字）
X=np.arange(len(signal))
X=np.reshape(X, (len(X), 1))
# 出力変数
Y=signal
# 線形回帰モデルのパラメータフィッティング
ModelLR = LinearRegression() # 線形回帰モデルのセット
ModelLR.fit(X, Y) # パラメータ獲得（Xは時間、Yは信号値）
# 動的特徴量を返り値に
return ModelLR.coef_[0]
関数を利用して動的特徴量を7チャンネル分取得¶
# 動的特徴量 取得
FV_LinearCoef=[]
for i in range(0, len(sig7[0,:])): # 7週分
FV_LinearCoef.append(Calc_LinearCoef(sig7[:,i]))
# 表示
ChannelName = ["Xacc", "Yacc", "Zacc", "Xang", "Zang", "Zang", "Mag"]
for i in range(0, len(sig7[0,:])): # 7週分
print(ChannelName[i], "の動的特徴量: ", round(FV_LinearCoef[i], 3))
このように、7チャンネル慣性センサデータの動的特徴量を算出できました。 今回、約5秒分の信号に対して、動的特徴量を算出しました。動的特徴量はほぼ0付近ですので、約5秒間を総合的に解釈して、右肩上がり、右肩下がりの傾向がないということになります。先ほど、X軸加速度で動的特徴量がマイナスになりましたが、あれは50番目から79番目、サンプリング周波数が100Hzですから、測定開始してから0.5秒〜0.79秒の間、0.3秒間分の信号を対象にしたからです。つまり、0.3秒間という小さい視点で見れば、X軸加速度は右肩下がりの部分もありますが、約5秒分を一気に対象にすると、そうでもない、ということになります。
まとめ¶
今回、時間領域の特徴量として、平均、分散、歪度、尖度、最小値、最大値、中央値、RMS、動的特徴量、合計9個の特徴量について説明しました。7チャンネル慣性センサデータの1チャンネル毎に9個の特徴量が得られるので、1回の観測に対し、合計63個の特徴量が抽出されたことになります。この63個の特徴を元に、いろいろな解析をしていくわけです。
人工知能あるいは人間は、複雑な情報を複雑なまま理解することは苦手です。今回の信号は、453点分、7チャンネルですから、$453 \times 7 = 3171$個の数字で構成されています。このような莫大な数字の集まりを、今回の特徴量変換により、63個にまで圧縮することができました。3171次元のベクトルが63次元まで低次元化されたとみなすことも可能です。すなわち、情報の複雑さは激減しました。また、一つ一つの数字も物理的な意味を保持しています（X軸加速度は右肩上がりだとか、Z軸加速度のばらつきは大きいだとか）。そのため、信号が有する物理的な状況をたいへん考えやすくなります。
これらの特徴量を人工知能に入力し、センサが装着された機械が壊れているだとか、センサが装着された人間はかなりいい腕立て伏せをしていただとか、判断させることができるかもしれないわけです。また、不調の機械と良好な機械に慣性センサを接続し、特徴量を算出します。そして、その特徴量を不調-良好の2群に振り分け、t検定などで有意差を算出してはどうでしょうか。もしX軸加速度の分散に有意差が生じていれば、その特徴が、不調-良好の識別にたいへん有効、と判断することもできます。以上のように、信号を特徴量に変換することで、いろいろな可能性が見えてきます。