PyLearn SIG01: 7チャンネル慣性センサデータ¶
物体の動きを測定するデバイスとして、慣性センサというものがあります。多くの慣性センサは、3軸加速度、3軸角速度を測定することができます。人間に慣性センサを装着したならば、人間の動きを測定できますし、機械に装着したならば、機械の動きを測定することができます。ここでは、慣性センサデータの基本について説明します。
センサデータのロード¶
センサデータを以下においてあるので、ダウンロードしてください。 http://int-info.com/PyLearn/dat.zip
これは、あるモータの隣に慣性センサを接続し、モータを回転させたときに得られたデータとなります。 こちらを開いてください。列の意味は以下の通りとなります。
- 0列目: 測定情報（ags -> 加速度、角速度）
- 1列目: コンピュータ時間（ミリ秒）
- 2列目: X軸加速度
- 3列目: Y軸加速度
- 4列目: Z軸加速度
- 5列目: X軸角速度
- 6列目: Y軸角速度
- 7列目: Z軸角速度
こちらをロードして、pythonの変数に代入してみます。skiprows=0は0行分スキップして読み込むという意味があります。たまにヘッダ（0列目は時間、1列目はX軸加速度のような文字情報）が入っている場合がありますが、その場合はskiprowsを1にして1行目を呼び飛ばします。usecols=range(2,8)は2列目〜7列目を読み込むという意味になります。0列目から開始なので、0と1列目は読み込まない、ということになります。
import numpy as np
signal = np.loadtxt("dat/1.csv", delimiter=",", skiprows=0, usecols=range(2,8))
次に、代入された配列サイズとその値を見てみます。
print("サイズ", np.shape(signal))
print(signal)
plot関数による可視化¶
453行、6列のデータということがわかります。この6列は、3軸加速度＋3軸角速度それぞれのデータをいう意味があります。中身を表示してみましたが、数字が多すぎてよくわかりません。そのため、plotにより可視化してみます。
import matplotlib.pyplot as plt
fig = plt.figure(figsize=(20, 8)) # 描画領域の横、縦のサイズ
# 3軸加速度
plt.subplot(2, 1, 1)
plt.plot(signal[:, 0], label="Xacc")
plt.plot(signal[:, 1], label="Yacc")
plt.plot(signal[:, 2], label="Zacc")
plt.xlabel("time: 0.01[sec]", size=20)
plt.ylabel("Accelecration: 0.1[mG]", size=20)
plt.legend(fontsize=20) # 凡例表示
# 3軸角速度
plt.subplot(2, 1, 2)
plt.plot(signal[:, 3], label="Xang")
plt.plot(signal[:, 4], label="Yang")
plt.plot(signal[:, 5], label="Zang")
plt.xlabel("time: 0.01[sec]", size=20)
plt.ylabel("gyro: 0.01[dps]", size=20)
plt.legend(fontsize=20) # 凡例表示
plt.show()
こんな感じで可視化できました。よく見ると、Z軸加速度だけ10000という値を示していますね。これを計算してみると、$10000 \times 0.1 \times 10^{-3} = 1$[G]となります。1Gは、地球の重力加速度です。あらゆる物体は、重力により地球の中心方向に加速していますので、このような値が現れるわけです。
合成加速度¶
次に、合成加速度について説明します。合成加速度とは、軸成分を考慮しない物体の動きを意味します。ある時間$t$における合成加速度$M(t)$は、
$M(t) = \sqrt{X_{\rm acc}(t)^2+Y_{\rm acc}(t)^2+Z_{\rm acc}(t)^2}$,
で定義されます。2乗でマイナス方向の動きを全てプラスにして、各軸の和を算出しているので、センサ装着位置の全体的な動きを意味しています。プラスマイナス方向や軸成分の情報を消失させているので、合成加速度では細かな動きはわかりませんが、全体的な動きをスカラー値という単純な形で知ることが可能になります。これを求めるコードは以下の通りです。
# 合成加速度の算出
M = signal[:,0] * signal[:,0] + signal[:,1] * signal[:,1] + (signal[:,2]) * (signal[:,2])
M = np.sqrt(M) - 10000
これで合成加速度を算出できました。今回は重力加速度を除去した上で、合成加速度を求めてみました（いろんな流儀がありますが、単なる線形変換ですので、正直なところなんでもいいような気もします）。続いて、3軸加速度と3軸角速度、合計6列のデータに対し、その7列目に合成加速度を添加してみます。
# 配列サイズを変更（Mは1次元配列なので、2次元配列に変換）
print(np.shape(M)) #プログラム的にはここは不要
M_new = np.reshape(M, (len(M), 1))
print(np.shape(M_new)) #プログラム的にはここは不要
# 2次元配列 signal と 2次元配列 Mを結合
signal_7channel = np.concatenate([signal, M_new], 1)
ちょっとややこしいですが、Mは最初、1次元の配列（行成分しかない453行のデータ）です。これを、2次元の配列（453行1列のデータ）に変換します。np.reshapeの第一引数が変換対象の行列、第二引数が変換後のサイズです。(len(M), 1)におけるMは、Mの長さ分の行数を意味し、1は列数を意味します。これによって、453行の合成加速度Mを、453行1列の合成加速度M_newに変換できました。np.shapeにて行列のサイズをprintで見てみると、どういうことかわかると思います。
次に、np.concatenateにより、453行6列のsignalと453行1列のM_newを結合しています。concatenateにある1とは、列成分を結合させるという意味があります。ここに0を指定した場合、行成分で結合しようとするので、エラーが出ます。
これにより、6軸だったセンサデータが7軸になりました。これを7チャンネル慣性センサデータなどと呼んだりします。
サンプリング周波数 / サンプリング周期¶
信号を扱うとき、測定データのサンプリング周波数、あるいは、サンプリング周期が明らかでないと、まともな分析が行えません。サンプリング周波数とは1秒間に何点データを記録するか、サンプリング周期は何秒ごとにデータを記録するか、という意味があります。サンプリング周波数が200Hzの場合は1秒間に200点データを記録します。1秒間に200点なので、1/200=0.05となり、サンプリング周期は0.005秒となります。サンプリング周波数とサンプリング周期は各々逆数の関係があります。
今回、データ数が453行分ありました。これは一体何秒分のデータでしょうか。これは、サンプリング周波数が明らかでないとわかりません。サンプリング周波数が1000Hzであれば0.453秒分のデータになりますし、100Hzであれば4.53秒分のデータ、200Hzの場合は2.265秒分のデータとなります。訳がわかりません。なので、慣性センサでデータ測定を行う場合は、必ずサンプリング周波数/サンプリング周期のどちらかを明確にした上で、記録してください。これは後述する周波数領域特徴量を扱う場合にも、必須となります。
なお、今回のデータはサンプリング周波数を100Hzと設定してデータを取得したので、4.53秒分のデータというのが正解です。サンプリング周波数を高めるごとに、精密な測定が可能になります。一方、データサイズが増えるので、通信であるとか、解析であるとかに時間がかかるようになります。いい感じのサンプリング周波数の値を探索したり、熟練の勘や経験で決めてください。大前の個人的な体感ですが、瞬間的な衝撃が重要になる場合は、200Hz以上がいいです。サンプリング周波数を小さくすると、データの記録間隔が長く開くので、衝撃の取り漏らしが多くなります。
SF = input("サンプリング周波数はいくつ？ ") # サンプリング周波数を手入力
SF = int(SF) # string型をint型へキャスト
print("サンプリング周波数は", SF, "Hzです！")
print("データを", 1/SF, "秒おきに記録しています！")
print("この信号長は", str(len(signal_7channel[:, 0])), "です！")
print("すなわち、", str((1/SF)*len(signal_7channel[:, 0])), "秒分のデータです！")