PyLearnSIGHT 01: 視線データの特徴量¶
0. 信号・信号長・サンプリング周波数¶
ここでは、視線計測の解析方法の基礎について言及します。視線は通常、なんらかの画面を見たときのXY座標のデータとして記録されます。でも、1度の計測でこの2つの数字が出てくるわけではなくて、XY座標のデータは時間によって変化するわけですから、計測時間分だけ、その数字が記録されることになります。このように、時間に従って変化するデータを、信号、信号の長さを信号長と呼びます。X軸の信号、Y軸の信号、2つの信号が得られるわけです。
具体的な信号長を知るには、「1秒間に何回データを記録するのか」という情報がなければいけません。例えば、10秒間のデータ測定をした状況を想定します。このとき、1秒間に10回データを記録するなら、信号長は100になります。でも、1秒間に50回データを記録するなら、信号長は500になります。このように、測定時間が同じでも、1秒間にデータを記録する回数がわからなければ、信号長がどれくらいになるかわかりません。
したがいまして、「 1秒間にデータを記録する回数」というのは、信号の計測においてすごく大事な情報です。これを、サンプリング周波数と言います。単位はHzになります。例えば、サンプリング周波数が100Hzであれば、1秒間に100個のデータを記録することになります。20秒間の計測なら、信号長が2000だ！、などがわかるわけです。逆に、信号長2000の信号があれば、20秒間の計測だったんだな、ということもわかることになります。
でも、サンプリング周波数がわからないと、ただただデータがあるだけで、何秒の計測だったのかがわからなくなってしまいます。したがって、データを記録・解析する際は、サンプリング周波数が幾つなのか、確実に把握するようにしてください。学会発表や論文執筆の際には、絶対に書かなければいけない情報なので、注意しましょう。
1. データのロード¶
今、jupyterを起動しているディレクトリに、データを保存する「dat_eye」という下位ディレクトリがあり、そこに「csv」拡張子の視線計測データが保存されているとします。
この状況で、視線データをロードしてみます。……の前に、ちょっとおまじない。コンピュータの初学者はスルーしていいです。
# Q: これは何？
# A: 測定データの「文字コード」を、UTF8に変換するコードです。
# 測定環境によって文字コードがバラバラになったり、
# Excelが勝手に文字コードを変えてきたりするので、文字コード統一が必要になります。
# 謎のエラーは大体文字コードが原因なので、やらないとめんどくさいことが起こるかも？
import os
import glob as gl
filename = gl.glob("dat_eye/*.csv")
# brew install nkf（ターミナルで、必須）
for i in range(len(filename)):
os.system("nkf -w8 --overwrite "+ filename[i])
おまじないが済んだところで、データのロードをしてみます。np.loadtxt関数がファイルの読み取り、delimiterはデータが区切られている記号（csvの場合は、カンマ）、skiprows引数は上何行を読み取らないか、dtypeは読みった後のデータ型です。
何行スキップすればいいのかなどは、自分でcsvを開いて、チェックしましょう。また、普通、信号は実数（float型）で扱う必要がありますが、視線計測データには時間など、floatでは扱えないデータも含まれているので、一旦str型で読み取っておかなければいけません。dtypeをfloatにすると、実数以外のデータがあるので読み取れないよとエラーが出てしまいます。
import numpy as np
# skiprowsで、上の行をロードをスキップ
filename = "dat_eye/20191118-155920.csv"
dat = np.loadtxt(filename, delimiter=",", skiprows=25, dtype="str")
# print(dat) # データを見たい場合
# print(np.shape(dat)) # データのサイズを見たい場合
ロードが終わったら、きちんと想定通りのデータが入っているか、printで確認してみましょう。
続いて、データを見て、XY軸のデータだけを取り出してみます。このために、numpyの配列を定義しておきます。英語で、rowは行（下に増えていくデータ）、colは列（右に増えていくデータ）です。
# 何行何列のデータを格納したいか、算出
size_rows = len(dat[:, 0]) #行サイズ（信号長だけ必要）
size_cols = 2 #列サイズ（XY座標なので、2列でok）
# 行サイズが信号長で、列サイズが2の零行列
npdat = np.zeros([size_rows, size_cols])
# ロードしたデータから、XY座標の信号を、npdatに代入
for i in range(size_rows):
npdat[i, 0] = float(dat[i,9]) # CX（9列目）
npdat[i, 1] = float(dat[i,10]) # CY（10列目）
これで、npdatの0列目にX座標の信号、1列目にY座標の信号を、代入することができました。信号長と計測時間を表示してみます。信号長は行の長さ、測定時間は信号長をサンプリング周波数で割った値になります。今回の測定器はサンプリング周波数が60Hzらしいので、その値を採用します。
sf = 60 # サンプリング周波数
print("信号長: ", size_rows)
print("測定時間: ", np.round(size_rows/sf, 2), " sec")
2. エラー値の平均値補完¶
8分弱の計測であることがわかりました。次に、X, Y座標で、どんな信号が得られたのか可視化してみます。
import matplotlib.pyplot as plt
# X軸
plt.subplot(2, 1, 1)
plt.plot(npdat[:, 0])
plt.ylabel("x val")
# Y軸
plt.subplot(2, 1, 2)
plt.plot(npdat[:, 1])
plt.ylabel("y val")
plt.xlabel("time") # サンプリング周波数60で割ったら、秒になる
なんだかおかしなデータがたくさんあります。どうやら、変なところを向くと、値が極端に大きくなるようです。どうやら、仕様として、変なところを見ると、値が0か999になるらしいです。このような不適合な数をカウントしてみます。
x_999 = np.sum(npdat[:,0] == 999)
print("X座標999の数: ", x_999)
x_000 = np.sum(npdat[:,0] == 0)
print("X座標0の数: ", x_000)
y_999 = np.sum(npdat[:,1] == 999)
print("Y座標999の数: ", y_999)
y_000 = np.sum(npdat[:,1] == 0)
print("Y座標0の数: ", y_000)
数百程度、よくない値があることがわかりました。このようなデータは、その信号の平均値で置き換えましょう。
# まず、999を0に変換
for i in range(size_rows):
if npdat[i, 0] == 999:
npdat[i, 0] = 0
if npdat[i, 1] == 999:
npdat[i, 1] = 0
# 次に、各信号の平均値を算出
temp_x = npdat[npdat[:, 0] != 0, 0] # 信号xの0を除去
temp_y = npdat[npdat[:, 1] != 0, 1] # 信号yの0を除去
mean_x = np.mean(temp_x) # 平均
mean_y = np.mean(temp_y)
# 信号の値が0のところを、平均値で補完
for i in range(size_rows):
if npdat[i, 0] == 0:
npdat[i, 0] = mean_x
if npdat[i, 1] == 0:
npdat[i, 1] = mean_y
# 可視化
# X軸
plt.subplot(2, 1, 1)
plt.plot(npdat[:, 0])
plt.ylabel("x val")
# Y軸
plt.subplot(2, 1, 2)
plt.plot(npdat[:, 1])
plt.ylabel("y val")
plt.xlabel("time") # サンプリング周波数60で割ったら、秒になる
おかしなデータがなくなっているように見えます。本当に消えたのかチェックしてみます。
x_999 = np.sum(npdat[:,0] == 999)
print("X座標999の数: ", x_999)
x_000 = np.sum(npdat[:,0] == 0)
print("X座標0の数: ", x_000)
y_999 = np.sum(npdat[:,1] == 999)
print("Y座標999の数: ", y_999)
y_000 = np.sum(npdat[:,1] == 0)
print("Y座標0の数: ", y_000)
確かに、999と0がなくなっていることを確認できました。
この信号は、散布図として可視化することで、画面上のどこをみているのかを表示することもできます。ただし、時間情報は消失するので注意してください。
plt.scatter(npdat[:, 0], npdat[:, 1])
plt.xlabel("x")
plt.ylabel("y")
十字をクロスするような視線移動が多いことがわかります。斜めでは、左上がやや多いけど、右上が少ない、などの情報もわかりますね。
3. 周波数領域¶
周波数領域とは、信号に対する振動の強さを解析する方法です。詳しい説明は PyLearn SIG 03: 周波数領域の特徴量 を参照ください。今回はサンプリング周波数が60Hzなので、ナイキスト周波数は30Hz、つまり30Hzまでの振動解析が可能です。理論は難しいですけど、結果を見るだけなら、結構簡単です。結果の見方がだいたいわかったら、実装してみます。
実装も、ちょっとややこしいです。信号処理の勉強をたくさんせねばならないので、慣れないうちはコピペで乗り切ってください。ただし、結果の解釈はできなくてはいけません。
from scipy.fftpack import fft
# パワースペクトル算出関数
# 入力: 1次元信号、サンプリング周波数
# 出力: パワースペクトルの値、その横軸になる周波数
def fft_sig(signal, sf):
L = len(signal) # 信号長
freq = np.linspace(0, sf, L) # 周波数領域の横軸
# フーリエ変換
freq = np.linspace(0, sf, L) # 周波数領域の横軸
yf = np.abs(fft(signal)) # 振幅スペクトル
yf = yf * yf # パワースペクトル
return [yf, freq]
# 作った関数を使用して、パワースペクトルを得る
ps_x = fft_sig(npdat[:, 0], sf)
ps_y = fft_sig(npdat[:, 1], sf)
これで、信号の周波数領域が計算できました。といっても意味不明だと思うので、結果を可視化してみます。
fig = plt.figure(figsize=(15, 6)) # 描画領域の横、縦のサイズ
plt.subplots_adjust(wspace=0.4, hspace=0.6) # グラフ配置の間隔調整
NyqFreq = int(len(ps_x[0])/2) # ナイキスト周波数の配列番号
# 信号Xのパワースペクトル
plt.subplot(2,1,1)
plt.plot(ps_x[1][0:NyqFreq], ps_x[0][0:NyqFreq])
plt.xlabel("Hz", size = 12)
plt.ylabel("Power", size = 12)
plt.title("Frequency Domain of X", size = 12)
plt.xlim([0, 0.5]) # 横軸: 0Hz 〜 0.5Hzまでを可視化（30Hzまで見れる）
plt.ylim([0, np.power(10, 12)]) # 縦軸: 0 〜 10^12まで
# 信号Xのパワースペクトル
plt.subplot(2,1,2)
plt.plot(ps_y[1][0:NyqFreq], ps_y[0][0:NyqFreq])
plt.xlabel("Hz", size = 12)
plt.ylabel("Power", size = 12)
plt.title("Frequency Domain of Y", size = 12)
plt.xlim([0, 0.5]) # 横軸: 0Hz 〜 0.2Hzまでを可視化（30Hzまで見れる）
plt.ylim([0, np.power(10, 12)]) # 縦軸: 0 〜 10^12まで
plt.show()
X信号のについては、0.025Hzくらいの成分が強いでしょうか？ 0.025Hzの逆数は、40秒です。すなわち、40秒に1回振動するような成分（そのくらいのペースで視線を左右に動かす）が強いことがわかります。また、0.050Hzもちょっと強いので、20秒に1回の振動もあるかもしれません。一方、Y信号の方は、あまり強くありません。視線の縦移動に対する規則はないのかもしれません。
注意点: 本来、周波数解析は定常波（長時間同じ振動をする波）に対して行うものです。視線の動きは非定常波なので、信号全体にfftを適用するのは不適切かもしれません。周波数から得られる特徴量を使用する場合は、相談が必要です。
4. 特徴量の算出¶
視線の動きに対する、時間領域（横軸が時間）、周波数領域（横軸が周波数）がわかったところで、特徴量を算出してみます。特徴量の算出とは、莫大な数字で構成される信号から、物理的に意味のある情報を単純化した上で取り出すことを意味します。なぜ特徴量を求めなければいけないか？それは、生の信号は数字が多すぎて、分析ができなくなるためです（今回のような8分弱の信号であれば、1信号で2万個以上、X, Y信号の2つがあるので、合計4万個の数字の集まりになります）。したがって、特徴量に変換しなければいけません。
今回は、時間領域の特徴量として、
- 平均: 信号の平均値
- 分散: 信号のばらつき
- 歪度: ヒストグラム化したときの、左右の裾の偏り
- 尖度: ヒストグラム化したときの、尖り具合
- 絶対平均: 信号の符号をすべてプラスにしたときの平均値
- 最大: 信号の最大値
- 最小: 信号の最小値
を算出してみます（メモ: マイナスをはかないのが標準なら、絶対平均入らないかも？）。
import scipy.stats as st
# 特徴量
print("*** X座標の信号に対する時間領域特徴量 ***")
print("X平均: ", np.mean(npdat[:, 0]))
print("X分散: ", np.var(npdat[:, 0]))
print("X歪度: ", st.skew(npdat[:, 0]))
print("X尖度: ", st.kurtosis(npdat[:, 0]))
print("X絶対平均: ", np.mean(np.abs(npdat[:, 0])))
print("X最大: ", np.max(npdat[:, 0]))
print("X最小: ", np.min(npdat[:, 0]))
print("X中央値: ", np.median(npdat[:, 0]))
print("")
print("*** Y座標の信号に対する時間領域特徴量 ***")
print("Y平均: ", np.mean(npdat[:, 1]))
print("Y分散: ", np.var(npdat[:, 1]))
print("Y歪度: ", st.skew(npdat[:, 1]))
print("Y尖度: ", st.kurtosis(npdat[:, 1]))
print("Y絶対平均: ", np.mean(np.abs(npdat[:, 1])))
print("Y最大: ", np.max(npdat[:, 1]))
print("Y最小: ", np.min(npdat[:, 1]))
print("Y中央値: ", np.median(npdat[:, 1]))
それぞれ、print関数の中に使用している関数で算出可能です。小数点以下が多くてちょっとみにくいですが、こんな感じで求められます。
また、周波数領域の特徴量として、
- パワーバンド（特定周波数帯のパワースペクトルの総和）: その成分の強さ
- 周波数領域エントロピー: 周波数領域全体の曖昧さ
- エネルギー: 周波数全体のパワーの強さ
があります。詳しい説明は PyLearn SIG 03: 周波数領域の特徴量 を参照ください。ただし、初学者には難易度が高いので、時間領域の特徴量だけの方が良いかもしれません。
# パワーバンドを求める関数
# signal: 1次元信号
# sf: サンプリング周波数
# min_f: パワーバンドを算出する周波数の最小値
# max_f: パワーバンドを算出する周波数の最小値
def Calc_PowerBand(signal, sf, min_f, max_f):
# パワースペクトル算出
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
# 周波数領域エントロピーを求める関数
# 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
# エネルギー
# signal: 1次元信号
# sf: サンプリング周波数
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
上から順に、パワーバンド、周波数領域エントロピー、エネルギーを求める関数です。時間領域特徴量と違って、初めから関数が用意されていないので、自分で作る必要があります。
# 特徴量
print("*** X座標の信号に対する時間領域特徴量 ***")
print("Xパワーバンド_0.000〜0.025Hz: ", Calc_PowerBand(npdat[:,0], sf, 0, 0.025))
print("Xパワーバンド_0.025〜0.050Hz: ", Calc_PowerBand(npdat[:,0], sf, 0.025, 0.05))
print("X周波数領域エントロピー: ", Calc_FreqEnt(npdat[:,0], sf))
print("Xエネルギー: ", Calc_Energy(npdat[:,0], sf))
print("")
print("*** Y座標の信号に対する時間領域特徴量 ***")
print("Yパワーバンド_0.000〜0.025Hz: ", Calc_PowerBand(npdat[:,1], sf, 0, 0.025))
print("Yパワーバンド_0.025〜0.050Hz: ", Calc_PowerBand(npdat[:,1], sf, 0.025, 0.05))
print("Y周波数領域エントロピー: ", Calc_FreqEnt(npdat[:,1], sf))
print("Yエネルギー: ", Calc_Energy(npdat[:,1], sf))
周波数領域の特徴量が求まりました。値が大きくてわかりにくいです。まあ、大事なのは相対的にどっちが大きいかです。
例えば、X信号の0〜0.025Hzと0.025〜0.050Hzのパワーバンドを比較すると、前者の方が値が大きいです。つまり、視線の横振動は、0〜0.025Hzの方が活発である、ということがわかります。
またX信号とY信号の0〜0.025Hzのパワーバンドを比較すると、X信号の方が3〜4倍くらい大きいです。つまり、この成分の振動は、縦よりも横の方がはるかに強いことがわかります。
5. 特徴量を求めた後は？¶
基本的にやることは、以下の4つです。順番は関係ありません。興味があることをやるといいでしょう。
(0) 便利な関数を作っておく。
分析をする際、取得した信号一つ一つの特徴量を計算するコードを、毎回書くのは大変です。そのため、例えば、ファイル名を入力したら、すべての特徴量を一気に計算して、出力してくれる関数があればどうでしょうか。一度作るだけで、一行コードを書くだけで特徴量を計算することができます。これは、今後の分析をする上で、大変便利です。
(1) 新たな特徴量がないか検討する。
今回説明した特徴量は、信号でよく使われるものをピックアップしたにすぎません。視線解析に有効な特徴量は、もっとたくさんあるでしょう。それを先行研究から調べて、自分で実装すると、分析の幅が広がります（視線移動のユークリッド距離やマンハッタン距離などがあるようです）。また、既存研究にとらわれず、自分で新たな特徴量を開発するというのも良いことです。
(2) 複数の人間から視線計測を実施し、特徴量を比較する。
もし、男女差を調べたいなら、男性と女性の特徴量を比較します。例えば、男性の方が女性よりも「X信号の分散が大きい」ということがわかれば、男性は左右に視線をばらつかせる傾向がある、など、これまでわからなかった新たな事実を知ることができます。他には、ある製品に興味があるグループと、興味がないグループの、製品広告を見る際の視線の動き、車酔いしやすいグループ、しにくいグループなど、いろんな比較があると思います。生の信号を比較しても意味がさっぱりわかりませんが、特徴量を比較することで、いろいろなことがわかります。
(3) 特徴量を入力し、何らかの状態判定を実現する人工知能を構築する。
視線の動きから性別を明らかにしたい、車酔いしやすいしにくいを判定する人工知能を作る場合、生の信号をそのまま入力させるのではなく、だいたい、信号の特徴量を入力して、状態を予測するというアプローチが取られます（生の信号から判定するアプローチもありますが、その場合は、必要なデータ数が爆発的に増えます）。
Appendix 関数化のヒント（上記(0)の課題）¶
以下のコードは、ファイル名を入れたら、X信号、Y信号の平均値2つを出力してくれる関数です。
def FeatureVector(filename):
# データロード
dat = np.loadtxt(filename, delimiter=",", skiprows=25, dtype="str")
# 何行何列のデータを格納したいか、算出
size_rows = len(dat[:, 0]) #行サイズ（信号長だけ必要）
size_cols = 2 #列サイズ（XY座標なので、2列でok）
# 行サイズが信号長で、列サイズが2の零行列
npdat = np.zeros([size_rows, size_cols])
# ロードしたデータから、XY座標の信号を、npdatに代入
for i in range(size_rows):
npdat[i, 0] = float(dat[i,9]) # CX（9列目）
npdat[i, 1] = float(dat[i,10]) # CY（10列目）
# まず、999を0に変換
for i in range(size_rows):
if npdat[i, 0] == 999:
npdat[i, 0] = 0
if npdat[i, 1] == 999:
npdat[i, 1] = 0
# 次に、各信号の平均値を算出
temp_x = npdat[npdat[:, 0] != 0, 0] # 信号xの0を除去
temp_y = npdat[npdat[:, 1] != 0, 1] # 信号yの0を除去
mean_x = np.mean(temp_x) # 平均
mean_y = np.mean(temp_y)
# 信号の値が0のところを、平均値で補完
for i in range(size_rows):
if npdat[i, 0] == 0:
npdat[i, 0] = mean_x
if npdat[i, 1] == 0:
npdat[i, 1] = mean_y
# *** ここまで、エラー値の補完処理 ***
# （特徴量変換）X信号とY信号の平均値を算出
FV_mean_x = np.mean(npdat[:, 0])
FV_mean_y = np.mean(npdat[:, 1])
# → この下にどんどん特徴量算出のコードを追加していく
# 出力
return [FV_mean_x, FV_mean_y]
関数化できたら、以下のように数行で特徴量を算出することができます。これは大変便利ですので、ぜひ、検討しましょう。毎回組んでると大変です。
fname = "dat_eye/20191118-155920.csv"
a = FeatureVector(fname)
print("X平均", a[0])
print("Y平均", a[1])