======================
(野村さんへのメモ)
・図2-6-4はネットから落としたので、そのまま使わずトレースして下さい
・「★補足資料」の部分は囲み記事としてフォントを落として添えて下さい
・リスト2-6-3は長くて無理なら以下にずっと置きます
	http://nagasm.org/ASL/mbed3/fig4/CQ_nagasm_list_2-6-3.txt
======================

■第4部　マイコンでトライ！ 生体情報の信号処理の基礎■

第2章　mbedによる筋電情報システムの例

2-6　FFTと筋電情報パターン認識

本節では、前節までに製作したNucleoF401REによる筋電センサ・システムを進化させて、単に筋電情報をセンシングする(計測する)だけでなく、リアルタイムに筋電(生体)情報を信号処理することで、何らかの「意味」を見出す(パターン認識)ようなシステムを目指します。筋電センサで歴史的にもっとも古いパターン認識の応用は「義手の制御」、つまり「筋電義手」ですが、調べてみると多くの文献がありますので、ここでは省略します。

ここでのトピックとして大きく3つ、まずは2チャンネル筋電センサの情報を2次元パターン認識する「リサジュー図形」によるアプローチ、次に定番のFFTによるスペクトル分析、そしてニューラルネットワークによる認識応用について検討します。

●リサジュー図形によるパターン認識

まずは図2-6-1を見てみましょう。これは前節2-5の最後、リスト2-5-5のNucleoF401REプログラムで筋電センサ・システムについてはまったくそのままで、XBee経由で送られてくる情報を表示するホストPC側のMax6のパッチを改訂したスクリーンショットです。画面中央の2段の筋電信号グラフはそのままですが、右下に正方形のスクリーンが新しく登場しました。この右下スクリーン内にぐちゃぐちゃ動いている「軌跡」に注目して下さい。これはx座標とy座標に2チャンネルの筋電情報を入力して「lineto」コマンドで、つまり「その座標を終点とする直線を引く」→「始点をこの座標に移動する」という処理を連続している軌跡です。

(図2-6-1　リサジュー図形によるホスト側のMax6パッチ画面例(1))

これは、本システムの2チャンネルの筋電情報のように、相互に何らかの関係がある2系統の情報を視覚化する「リサジュー図形」というもので、例えばサインとコサインだと円になるように、両チャンネルの情報の関係が浮き出てきます。ここから何らかの「ジェスチャー」・「ニュアンス」・「コントロール」・「コミュニケーション」のような情報を抽出できる可能性をがあります。

図2-6-1は、第2章の本システムで製作した筋電センサの電極を腕に付けて２チャンネル出力を「USBオーディオ治具」を介して記録した筋電信号データの1本目です。画面の筋電波形のようにちょっとレベルが小さめなので、右下の2次元グラフ内の軌跡がこじんまりとしています。さらに両チャンネルの変化がほどよく同期しているので、右上がりの傾きに沿っているのがわかります。2チャンネルの入力が完全に同一だと、右上がり45°の直線(y = x)の上に乗りますが、本システムでは、NucleoF401REからホストに送るシリアル情報を7.5msec間隔で交互に(各チャンネルは15msecサンプリング)送っているので、厳密には必ず同期からちょっとずれています。

図2-6-2は筋電信号データの2本目で、基本的には右上がりの傾向がありますが、ときに大きく離脱して、つまり各チャンネル間の信号に大きな違いがある瞬間が確認できます。

(図2-6-2　リサジュー図形によるホスト側のMax6パッチ画面例(2))

図2-6-3は筋電信号データの3本目で、前半は上チャンネルの筋電信号が良好に出ているのに、下チャンネルが低レベルで沈黙している区間です。この場合、右下の2次元グラフ内の軌跡としては、下の方に「濃い」(なんども行き来した)エリアがあり、左右方向の変動(上チャンネル)があるのに、上下方向の変動(下チャンネル)があまり無かった、という状況が推測できます。このように、リサジュー図形を眺めると、片方の電極の入力が不調だった、というような「意味」を抽出することができます。

(図2-6-3　リサジュー図形によるホスト側のMax6パッチ画面例(3))

ここまでは、NucleoF401REの筋電システムからXBeeでホストPCに刻々の筋電情報を送って、ホストPC上でその解釈を行ったわけですが、せっかく余力のあるNucleoF401REなので、なんとか簡単な「パターン認識」みたいなものを行って、その結果をホストに送りたくなります。ここまでの例では、XBeeで7バイトフォーマットでシリアル伝送するトラフィックの上限から、数msec〜10msecあたりのサンプリングでしたが、このシリアル転送を止めてNucleoF401RE内部で処理すれば、データ解釈のスピードは1桁ないし2桁、もっと高速にできる感触があります。

ただし、時間的に変化する生体情報をこのように「静止画」にしてしまうので、コマ送りアニメのように、「静止画にして分析する」サンプリングのスピードを上げる必要があります。また、図2-6-1の右下グラフのように、2次元プロット画像データをそのままホストに送るのでは、画素数を考えれば判るようにデータ量が増えて重くなります。上の例では「lineto」で描画した「軌跡」ですが、2チャンネルの情報を座標としてプロットすると、図2-6-4のように点が分布するグラフとなります。この雲状の分布のパターンとか点の密度などを分析する、つまり、リアルタイムな「統計分析」ということになります。

(図2-6-4　リサジュー図形のいろいろなプロットパターンの例)

◆Maxプログラムによるシミュレーション検討

この2チャンネルの筋電情報によるリサジュー図形の様子をどのようにNucleoF401REで処理するかについて、まずホストPCの環境「Max6」でシミュレーション実験しました。図2-6-5から図2-6-7までの3つのスクリーンショットでは、USBオーディオ治具を使って記録した筋電センサ信号をステレオ再生して、「10msecサンプリング、検波+20段移動平均による平滑」信号を求めて、これを正方形ウインドウのリサジュー図形のxとyの座標にしてプロットしています。全部で30秒間のデータ時間を8分割して、3750msecごとに「左から右」、さらに上段から下段へと次々に切り替えて表示しました。スクリーン上方にある2段の波形ウインドウは信号の全区間ではない(一部)ことに注意して下さい。

(図2-6-5　筋電記録データのリサジュー図形シミュレーション画面例(1))

(図2-6-6　筋電記録データのリサジュー図形シミュレーション画面例(2))

(図2-6-7　筋電記録データのリサジュー図形シミュレーション画面例(3))

これを見ると、図2-6-5の「記録(1)」では、8つの正方形ウインドウの区間の多くでシルエットが45度よりも寝ていて、上段の筋電レベル(x方向)が下段の筋電レベル(y方向)よりも大きい(優勢)こと、しかし全体としては右上がりの相関があり、両チャンネルの筋電情報はタイミングが合っていること、などが判ります。

また図2-6-6の「記録(2)」では、多くの区間でシルエットが45度よりも立っていて、上段の筋電レベル(x方向)が下段の筋電レベル(y方向)よりも小さい(劣勢)こと、全体としては右上がりの相関があることが判ります。

そして図2-6-7の「記録(3)」では、シルエットが45度よりも立っていて、上段の筋電レベル(x方向)が下段の筋電レベル(y方向)よりも小さい(劣勢)区間がほとんどであること、全体としては右上がりの相関があるが一部区間で歪んでいること(←センサ電極の接触不良だった区間?)、そして、ほとんど両チャンネルからの信号が消えている(脱力していた)区間がある事が判ります。

このように、筋電記録(30秒間)を8分割することで、時間的にリサジュー図形のパターンが変化することは比較的簡単に求められますが、この8つの正方形ウインドウ内にデータがプロットされている様子(時間的な経過)の情報は失われることになります。

そこで視点を変えて、2チャンネルの筋電情報によるリサジュー図形の特徴パラメータとして、3つのパラメータを演算した解析信号として表示してみたのが、図2-6-8〜図2-6-10のシミュレーション(USBオーディオ治具を使って記録した筋電センサ信号の(1)〜(3)に対応)です。画面に7段並んでいるグラフのいちばん上の2段は、筋電信号(x)と(y)の2チャンネル分です。

(図2-6-8　筋電記録データのリサジュー解析シミュレーション画面例(1))

(図2-6-9　筋電記録データのリサジュー解析シミュレーション画面例(2))

(図2-6-10　筋電記録データのリサジュー解析シミュレーション画面例(3))

3段目のグラフは、2チャンネルの筋電情報を成分とする「ベクトル(x,y)の大きさ」で、リサジュー図形としては、原点(正方形ウインドウの左下端)からプロット点までの距離となり、両チャンネルの信号が大きく異なっていると大きくなる、という意味を持ちます。式としては
	sqrt( x^2 + y^2 )
となる、時間に対応して刻々と得られる計算量です。

4段目のグラフは、この「ベクトル(x,y)の傾き」×「ベクトル(x,y)の大きさ」という計算量です。まず「ベクトル(x,y)の傾き」とは、yに比べてxが優勢(45度の傾きより寝ている)と小さくなり、xに比べてyが優勢(45度の傾きより立っている)と大きくなる、という意味を持ちます。ただしこれだとほとんど両方とも無信号の時に意味なく大きくなる事が判ったので、「ベクトル(x,y)の大きさ」が小さい時には、これと乗算することで小さく抑えるようにしたものです。式としては
	(y / x) * sqrt( x^2 + y^2 )
となる、時間に対応して刻々と得られる計算量です。

5段目のグラフは、ベクトル(x,y)から直線「y=x」までの距離という計算量です。これは、リサジュー図形で傾き45度、つまり直線「y=x」からどれだけ離れているか、というもので、自然に「ベクトル(x,y)の大きさ」の要因も反映されています。式としては、点(x,y)から直線「y=x」に垂線を下ろした足の座標との距離を求めるのですが、演算式を整理すると
	abs( x - y ) / 2
と簡単になりました。そしてこの式を眺めていて、ベクトル(x,y)から直線「y=x」までの距離として絶対値を求めるのはもったいない、正負の符号を付けて中央の「y=x」から上下に振ろう、としたのが6段目のグラフです。

これらのグラフがいちばん上の2段の筋電信号(x)と(y)とともに変化している様子をしばらく眺めていて、6段目のベクトル(x,y)から直線「y=x」までの符号付き距離の振る舞いで不満だったのは、2チャンネルの筋電ともにほぼゼロに近いところで振動していることで、これを抑えるには、4段目のグラフと同じ手法で、「ベクトル(x,y)の大きさ」と乗算すれば、より有効な筋電情報を強調できます。そこで、オフセットを別にすれば
	(y - x) * sqrt( x^2 + y^2 )
となる、時間に対応して刻々と得られる計算量を表示したのが、最下段、7段目のグラフです。

この量は、筋電が両チャンネルともに弱い、「力をあまり入れていない」時には、一種のアーティファクト抑止フィルタが効いた状態で中央値にへばりつき、そして両チャンネルの筋電の「差」、例えば腕を右捻り気味とか左捻り気味とかにアンバランスにした度合いに対応して正負に触れて、さらに筋電(力)の「強さ」にも対応しているので、単なる電極データの「検波+平滑」(パワー)の指標よりも意味のある情報と期待できます。今回は、この「リサジュー解析」のデータを「リサジュー図形によるパターン認識」と規定して、NucleoF401REに実装することにしました。

◆「リサジュー解析」の実装
 
以上の検討から、本システムの2チャンネル筋電情報に「リサジュー図形によるパターン認識」を適用して、実際に動くものとして完成したプログラムがリスト2-6-1です。プログラムを見やすくするために、
	・FIR_LPFの機能をカット
	・レンジ監視のアーティファクト・フィルタリング機能をカット
という変更をしています。これらの機能を再び追加するのは、リスト2-5-5を眺めれば容易です。この一方で、60Hz/50Hzのノッチフィルタは残してあり、実験的に記録した筋電信号のバックグラウンド・ハムノイズに対して、ON/OFFを切り替えることで、有効性を確認できます。

(リスト2-6-1)

また、2チャンネルのA/D変換を同じ割り込みタイミング内で呼ぶことで、時差による差分の誤差を抑えました。そして2チャンネルの筋電信号(「検波+平滑」したもの)を2系統のXBeeワイヤレス情報として15msecで交互に2チャンネル出力していたのを1本にまとめて5msecへ性能を上げました。中核となる処理はたった1行、
	data3 = (float)Lissajous / 10.0 * (data1 - data2) * sqrt( data1*data1 + data2*data2 );
というもので、ここでホストPCからは、入力ゲインとは別に「Lissajous」という定数(0〜127、初期値50)を設けて、この値を10で割ったゲイン(初期値 = 50/10 = 5倍)で、ベクトルの大きさによる重み付けをしています。

このようにNucleoF401REプログラムを実装して、USBオーディオ治具を使って記録した筋電センサ信号の(1)〜(3)に対応した出力を、ホストPCのMaxプログラムを改訂して表示したのが、図2-6-11〜図2-6-13です。それぞれのグラフは、中央値に対して、2チャンネルの筋電信号の片方が強いと上方向/下方向に変動し、その振幅は筋電情報の大きさ(パワー)に対応するとともに、さらに両チャンネルの筋電パワーを乗算することで、低レベルでの信号を抑止しています。ホストに届いたこの生体信号処理情報から、中央より上、下という方向、およびその大きさに対応して意味付けすることで、有効な「筋電ジェスチャー」認識アプリケーションを実現することが出来そうです。

(図2-6-11　筋電記録データに対するNucleoF401REリサジュー解析出力例(1))

(図2-6-12　筋電記録データに対するNucleoF401REリサジュー解析出力例(2))

(図2-6-13　筋電記録データに対するNucleoF401REリサジュー解析出力例3))

●FFTによるスペクトル抽出

生体信号の情報処理の最後のトピックは「FFTによるスペクトル抽出」です。筆者はこれまで、FFTはツール「Max」の便利なオブジェクトを便利な道具として使ってきたので(→Appendix5)、実際にFFTのプログラムをNucleoF401REなとマイコンのファームウェアとして書くのは初めてです。ここでは、簡単なFFTによるスペクトル抽出のアルゴリズムを実装した例を紹介します。

ネットで「FFT」と調べれば、色々な参考資料が豊富に入手できる時代ですが、mbedはオープンソース文化を受けて、以下のようにちょっと調べても14件ほどのサンプルのページ(→★補足資料)がヒットしました。本システムではその中で、
	・対象が筋電情報なのでオーディオ帯域ほどサンプリング周波数が高い必要が無い
	・2チャンネルの筋電信号を処理しているので2チャンネル化する
	・生体情報処理に必須のノッチフィルタと「検波+平滑(移動平均)」と共存させる
	・ホストPCへのXBee転送トラフィックの範囲内で情報をまとめる
というような条件を考慮しつつ、それぞれの内容を検討しました。

==================
★補足資料　ここから↓

Craig Evans / AudioIn
http://developer.mbed.org/users/eencae/code/AudioIn/

Frank Vannieuwkerke / KL25Z_FFT_Demo
http://developer.mbed.org/users/frankvnk/code/KL25Z_FFT_Demo/

不韋 呂 / UIT_FFT_Real
http://developer.mbed.org/users/MikamiUitOpen/code/UIT_FFT_Real/

Tony Abbey / KL25Z_FFT_Demo_tony
http://developer.mbed.org/users/tony1tf/code/KL25Z_FFT_Demo_tony/

Igor Skochinsky / FFT
http://developer.mbed.org/users/igorsk/code/FFT/

Daniel Dobano Fernandez / pruebas
http://developer.mbed.org/users/ddobano/code/pruebas/

Eli Hughes / CMSIS_DSP_401
http://developer.mbed.org/users/emh203/code/CMSIS_DSP_401/

Ale C.- / FFT_Example
http://developer.mbed.org/users/Suky/code/FFT_Example/

James Cobb / audio_FFT
http://developer.mbed.org/users/jcobb/code/audio_FFT/

Yue Kong / WoYaoChengGOng
http://developer.mbed.org/users/KongXiangyue/code/WoYaoChengGOng/

Nathan Lasseter / GraphicEqFFT
http://developer.mbed.org/users/User_4574/code/GraphicEqFFT/

Marcelo Rebonatto / PMED_Tempo
http://developer.mbed.org/users/rebonatto/code/PMED_Tempo/

Elliot Hernandez / SpectrumAnalyzer
http://developer.mbed.org/users/gth646f/code/SpectrumAnalyzer/

CQpub0 Mikami / FFT_Sampling
http://developer.mbed.org/users/CQpub0Mikami/code/FFT_Sampling/

★補足解資料　ここまで↑
==================

そして今回は、「Nathan Lasseter / GraphicEqFFT」(http://developer.mbed.org/users/User_4574/code/GraphicEqFFT/)を参考にして、本システムのNucleoF401REプログラムに組み込んでみました。2010年に公開されているこのコードは、2チャンネルのオーディオ信号をグラフィックイコライザとして周波数帯域ごとのレベル表示するもののようですが、呼び出している「fft」というモジュールがシンプルなので採用しました。

その「fft.cpp」の冒頭には、「この部分は私のオリジナルでなく、"http://www.mit.edu/~emin/source_code/fft/"からいただいたので中身に責任は持てない」と、マサチューセッツ工科大の「source_code/fft」というURLが書かれていましたが、現在は消えているようです。本稿はNucleoF401REのアルゴリズムとしてFFTを実装する目的のために、この中身はブラックボックスとして活用しています。

プログラムとしてはリスト2-6-2のようになりました。ここでは、ホストPCに2チャンネルの筋電FFT情報を送り、さらにプログラムを見やすくするために、
	・FIR_LPFの機能をカット
	・レンジ監視のアーティファクト・フィルタリング機能をカット
	・「リサジュー解析」の部分をカット
という変更をしています。これらの機能を再び追加するのは、リスト2-5-5/2-6-1を眺めれば容易です。この一方で、60Hz/50Hzのノッチフィルタは残してあり、実験的に記録した筋電信号のバックグラウンドハムに対して、ON/OFFを切り替えることで、有効性を確認しています。また平滑処理(移動平均)については、FFT演算をする前の、両チャンネルの筋電情報に対して、0.1msecのサンプリング周期で
	・検波(常にON)
	・ノッチフィルタ
の処理を行った後に移動してあります。

(リスト2-6-2)

FFTのサンプリングとしては1kHzで十分なので1msecとして、ここで16段のバッファに格納して、16段が揃ったところでFFT演算を行って、ホストPCに片チャンネルあたり16msecの間隔でシリアル伝送しています。FFTにつきものの「窓関数」については、もともと筋電信号はインパルスの集合体であり、連続波形信号を切り出す時の窓関数は本質的に不要なので省略しました。

16段のFFTを2チャンネル筋電情報に対して行うということは、毎回のデータは32個の浮動小数点データ(32ビット幅)となります。ここではNucleoF401REのプログラムとして同じ枠組みに入れてリアルタイム処理を実現する、という点にフォーカスしているので、ホストに転送する「7バイトフォーマット」を利用するため、以下のような情報圧縮を行っています。
	・FFTの16チャンネルのうち最初の14段だけを(とりあえず)出力する
	・14段の周波数バンドを2つずつ束ねて(2バンドの平均を出力)→7バンド情報とする
	・7バンドのレベルを3ビット(0〜7)に圧縮する
	・3ビット×7バンド=21ビットの上位3ビットに両チャンネルを区別するビットを持つ
ホストPCにこのようにまとめてXBeeで送るという制約が無ければ、こんなにデータを圧縮せず活用する事も容易ですので、この実装例はあくまでサンプルとして見て下さい。

このようにNucleoF401REプログラムを実装したところで、USBオーディオ治具を使って記録した筋電センサ信号の(1)〜(3)に対応した出力を、ホストPCのMaxプログラムを改訂して表示したのが、図2-6-14〜図2-6-16です。画面内には上下2段、2つの筋電チャンネルに対応したそれぞれ7つのウインドウが並んでいますが、この左端が周波数バンド(1)とバンド(2)の平均レベル、2番目が周波数バンド(2)とバンド(3)の平均レベル・・・と続きます。

(図2-6-14　筋電記録データに対するNucleoF401簡易FFT出力例(1))

(図2-6-15　筋電記録データに対するNucleoF401簡易FFT出力例(2))

(図2-6-16　筋電記録データに対するNucleoF401簡易FFT出力例(3))

かなり波形が激しく振動しているように見えますが、この横幅フルスケールで30秒間の記録筋電センサ信号より長い時間を描画しているので、3種類の記録筋電情報ごとに異なったスペクトル情報をリアルタイムに得ていることが判ります。
そして図2-6-17は、確認のために、このNucleoF401REシステムに対して、USBオーディオ治具から、片方のチャンネルには95Hz、もう片方には29Hzのサイン波を与えてみた確認実験の様子です。それぞれ、異なった周波数バンドのレベルが得られていることを確認しました。

(図2-6-17　NucleoF401簡易FFTの動作検証実験の画面例)

筋電情報についての研究では、筋電信号をスペクトル解析した場合に、周波数特性の変化は「筋肉の疲労」を反映している、という報告があります。上の実験では30秒間の筋電情報なので疲労の影響は検出できていませんが、筋電システムが無意識の筋疲労を検出できるとしたら、面白い応用ができるかもしれません。

●ニューラルネットワークによるパラメータ補間

「生体情報に対するパターン認識」という意味では、本節で紹介した「リサジュー解析」と「スペクトル解析」がそのための複数の情報を抽出するのに役立ちますが、実用を考えた場合にもっとも大きな壁となるのが、「個人差」です。ちょっと筋電電極になった気分で考えてみれば判りますが、例えば「腕」にセンサ電極を取り付けようとして、その人の
	・腕の太さ、手足の長さ
	・筋肉質か脂肪質か
	・皮膚が乾燥系か湿潤系か
	・電極圧迫ベルトがうまく効いているか
などの条件は、まさに千差万別です。その結果、AさんとBさんとで「同じような筋肉ジェスチャー」をしてもらっても、得られる筋電情報はまったく異なります(→Appendix5)。

このように「万人に共通のパターン認識」というのは、キネクトなどの画像処理で「ピース」や「コマネチ」や「ゲッツ!」などの身振りであればほぼ同じように検出できるのに対して、筋電情報では大違いです。そこで必要になるのが「リハーサル」であり、システムとしてはその「個性」「個人らしさ」を学習するアプローチです。

このような場合に有効な手法の一つが、「ニューラルネットワークによるパラメータ補間」です。第1章の1-3で紹介したように、ある種のブラックボックスとして、入力と出力との関係に明確な数理的関係が確立していなくても、「ほどよく」入力に対応した出力を出してくれる(そのために事前に教師情報の学習が必要)、というのがニューラルネットワークです。

リスト2-6-3は、だいぶ昔(22年ほど前)に筆者が作った、実際に走るリアルタイム・ニューラルネットワークのプログラムで、「ニューラルネットによるグラニュラシンセサイザ制御」というものです。グラニュラシンセサイザというのは楽音合成方式の一種なのですが、筋電情報と似ていて、多数のパルス的な要素が時間的にランダムに密集してクラウド的な音響を構成するので、フーリエ合成のような直感的なパラメータを指定できません。

(リスト2-6-3)

そこで筆者は、2次元のジョイスティックをコントローラとして、2次元空間内でジョイスティックをあれこれ動かして、2種類のパラメータに応じたグラニュラシンセサイザのサウンドのいくつかを「お気に入り」として指定して、これをニューラルネットワークの「教師データ」として学習させました。22年前のノートパソコン(PC-9801NT)で2週間ほどかけて、入力層2ノード、中間層16ノード、出力層2ノードの構成の倍精度変数の巨大な行列を「学習」すると、その後は、学習段階で与えたことのない2次元パラメータがやって来ても、そこそこ「妥当な」パラメータをリアルタイムに出力することを確認しました。

そこで、またまたmbedのオープンソース文化を期待して、mbedサイトで「neural network」と検索したのですが、残念ながら、これぞというサンプルはヒットしませんでした。関連するところとしては、巨大な「行列」を扱う必要があるので、「A class to handle Matrices」(http://developer.mbed.org/cookbook/MatrixClass)というページが参考になるかもしれません。

本稿の筋電センサシステムに、ここまで紹介したような「生体情報に対するパターン認識」のアルゴリズムをより実用的に実装した場合には、次に必要となるのは、リアルタイムに動く「柔軟なパラメータ補間」アルゴリズムとなり、この候補としてニューラルネットワークというのは、いまだ現役の有力候補です。いずれ機会があれば、このあたりをより深めた実験に挑戦してみたいと思っています。