======================
(野村さんへのメモ)
・★1★のフローは箇条書きよりもフローチャートに図示した方が判りやすいので、できればフローチャートにして下さい
・★2★の部分の10枚の図は、波形のシルエットの部分が判ればいいので、かなり小さな図で並べても結構です
======================

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

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

2-4　パワー検出と移動平均

●シミュレーションによる検波/平滑の検討

心電や脳波の場合には、生体リズムの信号としてある種の「波形」が明確に現われるので、その「周期」「振幅(変動)」「周波数成分」などを分析することになりますが、筋電はちょっと違います。前述のように、筋電信号とは「ほぼ同じレペルの神経繊維インパルスが時間的に多数、出現する」という信号なので、いわばインパルスはAM放送信号の「搬送波」で、「力を入れる」などを筋電のパワーから検出するには、AM放送信号の「変調波」(=包絡線[エンベロープ])を抽出する処理になります。この信号処理のイメージを確認するため、まず本システムのホストPC上の「Max」を使って、整流/検波と平滑の処理についてシミュレーションしました。

なお前章の1-3で紹介したように、平滑処理に使う「n段の移動平均」処理とは、「最新のn個のデータを平均する(累算して個数で除算)」という処理を刻々と行うものです。例えば20段の移動平均で最新入力データをD(n)とすれば、20段の移動平均値M(n)の出力は
	M(n) = { D(n) + D(n-1) + ・・・ + D(n-19) } / 20
となります。この20段のデータはシフトレジスタのように、入力されると逐次シフトして格納されるので、最新データが格納されるとともに、最古である20個前のデータが捨てられます。容れ物として20個分の変数を用意して、最古値を捨てたところにその直後の値を代入して・・・と次々にデータをいちいち1個ずつ古い場所に移動させてもいいのですが、FIFOバッファにデータを到来順に格納して、参照する書き込み/読み出しポインタを移動させる方法もあり、後述するNucleoF401REでの実装では、筆者は後者のFIFO方式を採用しています。

◆低い搬送周波数の場合の検波/平滑

図2-4-1は、画面左上で搬送波(キャリア)として、グラフが見やすいようにかなり低い周波数の[-1.0〜1.0]のサイン波形(a)を生成し、これに画面中央上のさらに低速な「0.2〜0.8の範囲でアップダウンする」振幅値(b)を乗算して、「仮想的な振幅変調波形」(c)を生成し、これを画面下段の3種類の検波プロックに供給して、その出力を比較するようにしたMaxパッチです。

(図2-4-1　Maxによる低いキャリア周波数での検波/平滑のシミュレーション画面例)

図2-4-1の下段左は「半波整流」(d)で、入力信号がゼロ以上の時にだけその値が出力され、負の場合にはゼロ出力となり、ダイオードの整流作用に対応します。 その下に、この出力を20段の移動平均した波形(e)と50段の移動平均した波形(f)が並んでいます。左側の20段の移動平均(e)ではゼロ付近の尖った部分がやや丸まったもののほぼ入力(d)そのままなのに対して、右側の50段の移動平均(f)ではかなり滑らかな変化に「平滑」されているのが判ります。 ただしこの50段の移動平均(平滑)出力(f)を2.0倍に拡大した、画面右端の4つの比較ウインドウの2番目(g)は、その上の「変調信号」(b)と比較すると、まったく似ていません。 この低い周波数で半波整流の平滑処理として移動平均を使うには、50段ではまだまだ不足しているわけです。

図2-4-1の下段中央は「全波整流」(h)であり、入力信号の絶対値が出力され、ダイオードを2個使った全波整流回路に対応します。 その下に、この出力を20段の移動平均した波形(i)と50段の移動平均した波形(j)が並んでいます。左側の20段の移動平均(i)では搬送波の振動がまだ残るものの、右側の50段の移動平均(j)ではほぼ平滑されて、「変調信号」と似た信号になっているのが判ります。 画面右端の4つの比較ウインドウの3番目(k)は、この全波整流の50段の移動平均(平滑)出力(j)を1.8倍に拡大したもので、搬送波周波数が低いためにドリフトしているものの、50段あればまずまず平滑できている事が判ります。

図2-4-1の下段右は「二乗検波」(m)であり、入力データ自身を乗算する処理により入力値の2乗が出力され、2乗特性を持ったトランジスタの二乗検波回路に対応します。変調信号は歪みますが、信号の変化を拡大強調する性質があり、音声のAM放送では歪むとしても、再現性より「変化」に注目したい生体情報処理では、この二乗検波にもメリットがあります。その下には、この出力を20段の移動平均した波形(n)と50段の移動平均した波形(p)が並んでいます。左側の20段の移動平均(n)では搬送波の振動がまだ残るものの、右側の50段の移動平均(p)ではほぼ平滑されて、「変調信号」の変化を強調した信号になっているのが判ります。 画面右端の4つの比較ウインドウのいちばん下(q)は、この全波整流の50段の移動平均(平滑)出力(j)を3.0倍に拡大したもので、ドリフトと2乗演算による線形性の崩れが確認できます。 

◆高い搬送周波数の場合の検波/平滑

図2-4-2は、画面左上の搬送波(キャリア)として、グラフでは周期性が見えないほど高い周波数の[-1.0〜1.0]のサイン波形(a)を生成し、これに画面中央上の低速な「0.2〜0.8の範囲でアップダウンする」振幅値(b)を乗算して、「仮想的な振幅変調波形」(c)を生成したMaxパッチで、図2-4-1との違いは搬送波の周波数だけです。

(図2-4-2　Maxによる高いキャリア周波数での検波/平滑のシミュレーション画面例)

図2-4-2の下段左の「半波整流」(d)と、その下の20段の移動平均した波形(e)と50段の移動平均した波形(f)、さらに50段の平滑出力(f)を2.0倍に拡大した、画面右端の4つの比較ウインドウの2番目(g)を見てみると、この周波数ではまったく包絡線(エンベロープ)が得られていないのが判ります。

図2-4-2の下段中央の「全波整流」(h)と、その下の20段の移動平均した波形(i)と50段の移動平均した波形(j)を見ると、50段の平滑出力(j)を1.8倍に拡大した画面右端の4つの比較ウインドウの3番目(k)では、良好に「変調信号」が再現されています。図2-4-1の例よりもドリフトの影響が低減しているのは、変調波周波数に対して搬送波周波数が十分に高いので平滑(積分)効果が大きくなるからで、筋電情報の信号処理においても、このように神経繊維インパルスの周波数特性を考慮する必要があります。

図2-4-2の下段右の「二乗検波」(m)と、その下の20段の移動平均した波形(n)と50段の移動平均した波形(p)を見ると、50段の平滑出力(j)を3.0倍に拡大した 画面右端の4つの比較ウインドウのいちばん下(q)では、こちらも良好に「変調信号」の変化を強調した信号になっています。 筋電の情報処理においては、筋肉の「力」の特性はもともと非線形なので、飽和しないレベルを考慮しつつ、小レベルでの変化を拡大強調する「二乗検波」も有効であると思われますが、とりあえずの結論としては、
	・全波整流(ゼロレベルよりマイナスでは折り返す)
	・移動平均の段数は50段程度あれば十分
ということになったので、これをNucleoF401REのシステムで実装する事にしました。

◆オフセットが乗った場合の検波/平滑

上述の検波/平滑の処理に関して、全波整流で折り返す「ゼロ点」がずれた場合、つまり筋電信号にオフセットが乗った場合の影響について、引き続き「Max」を使ってシミュレーションしました。図2-4-3と図2-4-4は、上述と同様に違いは「搬送波周波数の違い」だけであり、いろいろな検波/平滑ブロックに供給している「変調信号」もまったく同様ですので、解説は図2-4-3について述べます。

(図2-4-3　低いキャリア周波数でオフセットが乗った場合の検波/平滑のシミュレーション画面例)

(図2-4-4　高いキャリア周波数でオフセットが乗った場合の検波/平滑のシミュレーション画面例)

図2-4-3の画面左上では、低い周波数のサイン波形(a)を生成し、これに低速な「0.2〜0.8までの範囲でアップダウンする」振幅値(d)を乗算して、上述の実験と同様の「仮想的な振幅変調波形」(b)を生成し、さらにここでは8ビット精度の整数値(c)にも変換しました。このグラフの精度では、見た感じは(b)も(c)も同じです。

画面右側に6つ並んだグラフのうち、一番上の(d)はこの「変調信号」、次の(e)は「全波整流+50段の移動平均」を浮動小数点演算として実行した、前項と同じ良好な検波/平滑信号、3番目の(f)は、8ビット整数値としたデータを中央値の127より大きいか小さいかで判定した「全波整流+50段の移動平均」のデータで、グラフは(e)と同じになります。

問題はここからです。画面右側に6つ並んだグラフの4番目の(g)は、入力される8ビット整数値を全波整流として折り返す基準値を「120」としたもので、筋電信号にオフセットとしてわずかにプラスされているため、正しい「絶対値」から誤差が生じた状態を想定しています。その結果、グラフ(g)は理想的な検波/平滑の(f)に比べて、ヘンな波形に歪んでいます。同様に、画面右側に6つ並んだグラフの5番目の(h)は、折り返す基準値を「140」としたもので、筋電信号にオフセットとしてわずかにマイナスされているために、こちらも理想的な検波/平滑の(f)に比べて、余計に上下変動する歪が発生しています。また、画面右側の一番下の(i)は、折り返す基準値を「105」とさらに悪くしたもので、理想的な検波/平滑の(f)に比べて、周波数そのものが2倍と見なせるほど、大きく歪んでいます。

図2-4-4は、搬送波=キャリアとして、グラフでは周期性が見えないほど高い周波数とした以外は、まったく図2-4-3と同じ条件で、オフセットを与えた実験の結果です。搬送波周波数が高くなったので、画面右側に6つ並んだグラフのうち、一番上の「変調信号」(d)に対して、次の浮動小数点演算による「全波整流+50段の移動平均」(e)も、8ビット整数値演算による「全波整流+50段の移動平均」(f)も、ほぼ(d)と同じように「変調信号」を再現できています。

これに対して、画面右側に6つ並んだグラフの4番目の、折り返す基準値を「120」とした(g)は、「変調信号」(d)に細かいランダムノイズが乗っているような信号となり、「50段の移動平均」があまり効かなくなっています。グラフの5番目の折り返す基準値を「140」とした(h)でも同様で、一番下の折り返す基準値を「105」とさらに悪くした(i)では、ほとんどランダムノイズだけになってしまいました。

第3章で紹介するように、筋電センサ信号をアナログの全波整流回路とCR平滑回路を通す場合には、最初に大容量コンデンサを直列に入れる(AC結合、直流成分カットのHPF)ことで折り返しの基準ゼロ電圧に設定します。mbedなど単電源システムのA/Dコンバータに筋電信号を入力して、デジタル信号処理として「検波/平滑」する場合には、入力信号の中点電位を正確に設定することが重要だと確認できました。前節2-3で紹介したように、本システムの筋電センサ回路では、この正確な「+3.3V電源電圧の中点基準電圧」を発生するために、シャントレギュレータLM4051を採用しています。

◆「二乗平均平方根」との違い

ここまで、筋電信号のパワーを検出する「検波/平滑」処理として、
	・半波整流
	・全波整流
	・二乗検波
という検波/整流の処理と、それに続く移動平均の平滑処理を見てきましたが、もう一つ、シミュレーションすべき処理がありました。それが「二乗平均平方根」(root mean square, RMS)です。絶対値というのは、「信号がマイナスである場合に符号を変える」という計算だけでなく、「2乗した後で平方根をとる」という計算でも同じになります。一方、統計学で標準偏差を求めたり、物理学・電気工学で時間的に変化する信号の大きさを評価するための「二乗平均平方根」という手法があります。これは、個々のデータを二乗して総和をとり、最後にその「総和の平均」の平方根をとる、という計算量です。

後述するように、筆者はNucleoF401REで「全波整流」→「移動平均」を実装する際に、例えば50段の移動平均であれば、50個のFIFOバッファのようなデータ領域とそのボインタ、さらに「累算値」という変数だけを定義して、
==========
★1★ここから↓

	・A/Dコンバータから最新の入力データを得る

	・FIFOをボインタでアクセスして読み出し「最古値」とする

	・「累算値」に入力データを加算する

	・その結果から「最古値」を減算して新しい「累算値」とする

	・この「累算値」を50で除算して「検波/平滑」出力データとする

	・FIFOをボインタでアクセスして入力データを書き込む
		(最古値が書かれていた場所に最新値として上書きする)
		(→これが50個の中で最後までアクセスされない「最新値」)

	・FIFOのボインタ値を次に進めておく
		(→2番目に古かったデータが次に「最古値」になる)

★1★ここまで↑
==========
というアルゴリズムにしました。このメリットは、移動平均の「50個のデータの加算」という部分で、いちいち毎回、加算ループを50回も回さずに、加わった最新値の加算と、捨てて行く最古値の減算と、の2回だけで済ませるスマートさにあります。

「二乗平均平方根」の成績が「全波整流→移動平均」よりも明確に良好であるとすれば、このアルゴリズムに、入力データを「2乗する計算」と、移動平均として累算値を50で除算した後に「平方根を求める計算」が加わることになります。そこで上述と同じようなシミュレーションとして、設定条件はこれまでと同じで、「仮想的な振幅変調波形」が0.2〜0.8という振幅値レンジなので、ここに「+0.15」と「-0.15」というオフセットを加えて、中央から上下に「寄せた」波形を作って「二乗平均平方根」を計算して比較してみました。

結果は図2-4-5と図2-4-6のように、画面右側に5つ並んだグラフに、
	・サイン波に低速なアップダウン振幅値を乗算した「仮想的な振幅変調波形」
	・「全波整流」→「50段の移動平均」による検波/平滑信号
	・「二乗平均平方根(50段の移動平均)」による検波/平滑信号
	・オフセットで上に寄せた波形→「二乗平均平方根(50段の移動平均)」による検波/平滑信号
	・オフセットで下に寄せた波形→「二乗平均平方根(50段の移動平均)」による検波/平滑信号
を、低い搬送波周波数と高い搬送波周波数とでグラフ表示してみました。

その結果、とりたてて「二乗平均平方根」の成績が「全波整流→移動平均」よりも良好である事はなく、結論として、シンプルな「全波整流」→「移動平均」でOK、と判断しました。なお、この結果には、筋電信号のエンペロープ(パワー)の特性として「仮想的な振幅変調波形」のようにそれほど急激に変化しない、という部分が効いているので、他の生体情報の場合には、さらに検討する必要があるでしょう。

(図2-4-5　低いキャリア周波数でオフセットが乗った場合の二乗平均平方根のシミュレーション画面例)

(図2-4-6　高いキャリア周波数でオフセットが乗った場合の二乗平均平方根のシミュレーション画面例)

●NucleoF401REシステムでの実現

以上のシミュレーションを受けて、実際にNucleoF401REの筋電センサシステムに「全波整流→移動平均」の処理を実装してみました。リスト2-4-1のようになっていて、これまでの実験で確立してきた、
	・タイマ割り込み
	・XBeeによるシリアル通信
	・FIFOによるシリアル送受信のバッファリング
	・アナログ入力(A0〜A3)
の変数定義部分やメソッド類を、「sub.hpp」として切り分けました。そして「main.cpp」の部分に、「全波整流→移動平均」の処理と、ホストPCからのコマンドでパラメータを変更する「シリアル受信・解釈」処理を置いています。

(リスト2-4-1)

「main.cpp」ではリスト2-4-1のように、50μ秒間隔のタイマ割り込みの100発ごと、つまりサンプリング周期として5ミり秒(200Hz)ごとに、
	・A/Dコンバータからの入力データに「gain」設定値を乗算する
	・「全波整流ON」設定であればマイナスの時に符号を反転して全波整流する
	・移動平均処理メソッドを呼び出して移動平均処理(段数は設定値)
	・結果を整数化してXBeeからホストへの送信FIFOに積む
という処理を実行しています。サンプリング周波数が200Hzとやや低いのは、A/D変換したデータを「7バイトフォーマット」で刻々とそのままXBeeからWiFi送信しているからで、mbed内でより高度な信号処理を行って、その結果をホストPCへまとめて送るのであれば、もっとサンプリング周波数を高くすることは容易で、実際に2-5や2-6ではより高速な処理を実装しています。

また、ホストPCからXBeeシリアル受信したコマンドで設定しているパラメータは、
	・A/Dコンバータからの入力データを持ち上げる「gain」 : 0〜15倍
	・全波整流のON/OFF : {detection」
	・移動平均の段数 : ナシ(0)、5段(1)、10段(2)、20段(3)、50段(4)、100段(5)、200段(6)
というもので、いずれもpublic変数として定義して見やすくしています。

●サイン信号を使った検証実験

図2-4-7は、リスト2-4-1のプログラムを走らせて、前述(2-2節)の図2-2-9のようなセッティングで、USBオーディオ治具を経由してNucleoF401REに「0.2Hzのサイン波でリング変調した880Hzのサイン波」をA/D入力して、XBee経由で戻ってきた「全波整流→移動平均」のデータをグラフ表示している様子です。

(図2-4-7　NucleoF401REで実際に検波→平滑処理を行った結果をXBee伝送した画面例)

ホストPCから設定できるのは上記の3つのパラメータで、この例では「Gain = 3」と3倍増幅していること、「detection ON/OFF」のボタンがONになって全波整流していること、「5 mean」とあるので「5段の移動平均」処理であることが判ります。なお、エンベロープに追従するのに5段もあれば十分のように見えるのは、基準周波数が880Hzの場合なので、筋電の神経パルスに近い2kHzあたりにすると、5段や10段では波形がうまく平滑されませんから注意して下さい。

以下は、共通の入力信号に対して、ホストPCのMaxから設定パラメータを変更させて、NucleoF401REから返ってきた結果データのグラフ(Max画面だけ抽出)で、波形を見やすくするために、点プロットでなく個々のプロット値を高さとするバー表示にしています。
==========
(この10枚の図は、波形のシルエットの部分が判ればいいので、かなり小さな図として並べても結構です)
★2★ここから↓

図2-4-8は、「全波整流OFF」、ゲイン「1倍」、移動平均「ナシ」というもので、「0.2Hzのサイン波でリング変調した880Hzのサイン波」そのものです。

図2-4-9は、「全波整流ON」(これ以降のグラフは全てON)、ゲイン「1倍」、移動平均「ナシ」というもので、マイナス部分が折り返された全波整流波形となっています。

図2-4-10は、「全波整流ON」、ゲイン「1倍」、移動平均「5段」というもので、波形シルエットとして滑らかに平滑されています。

図2-4-11は、「全波整流ON」、ゲイン「2倍」、移動平均「5段」というもので、入力信号が小さい場合には、このように持ち上げる機能によって見やすくなります。

図2-4-12は、「全波整流ON」、ゲイン「3倍」(これ以降のグラフは全て3倍)、移動平均「5段」というもので、このぐらい持ち上げるのが理想的ですが、あまり欲張って持ち上げて「飽和」するのは、歪みとなるので絶対に避ける必要があります。

ここからは移動平均の段数を変化させていて、「全波整流ON」、ゲイン「3倍」は共通で、図2-4-13は移動平均「10段」、図2-4-14は移動平均「20段」、図2-4-15は移動平均「50段」、図2-4-16は移動平均「100段」です。

図2-4-17はおなじ「全波整流ON」、ゲイン「3倍」で、移動平均が「200段」ですが、このように100段とか200段になると、ゼロ点で折り返していたディップ部分まで平滑されて全体のレベルが上がり、変動波形としてのダイナミックレンジが低下してきます。これは明らかに「多過ぎる」という事になります。

★2★ここまで↑
==========

●筋電信号を使った検証実験

上述のサイン波形を使った実験で、NucleoF401REを使った筋電システムの動作が確認できたので、ここで前述(2-2節)の筋電センサ回路の製作時に、USBオーディオ治具を使って「記録」していた、テストの筋電信号(第1章1-2の図1-2-1の波形)を、ホストPCのMaxパッチで「再生」して、NucleoF401REのA/D入力に供給してみました。

図2-4-18は、画面右下に2段に並ぶ筋電信号のうち、下段のチャンネル(以降も全て同じ)を、「全波整流OFF」、ゲイン「5倍」(以降は全て5倍)、移動平均「ナシ」として信号処理(ほとんどA/D入力そのまま)して、XBee経由でホストに送信してきたデータです。筋電信号に特有の、神経パルスが多数、頻発している様子がわかります。

図2-4-19は、下段の筋電信号を、「全波整流ON」として(以降は全てON)、ゲイン「5倍」、移動平均「ナシ」としたものです。マイナス部分がきれいに折り返されています。

図2-4-20は、下段の筋電信号を、「全波整流ON」、ゲイン「5倍」、移動平均「5段」としたものです。上述の880Hzのサイン波と違って、ほとんど変化がない(平滑されていない)ので、移動平均の段数が不足しています。

画面右下の下段の筋電信号を、「全波整流ON」、ゲイン「5倍」は共通として、移動平均の段数を図2-4-21は「10段」、図2-4-22は「20段」、図2-4-23は「50段」として移動平均処理したものです。50段になると、シルエットがかなり滑らかになってきますが、それだけ出力が遅延して「激しい動き」を反映しなくなってくるので、丸く滑らかになればいい、というものでもありません。

そして図2-4-24の「100段」、図2-4-25の「200段」となると、波形シルエットは丸く滑らかになるものの、全体のレベルが上がり、変動波形としてのダイナミックレンジが低下してきます。実際の筋電信号を使ってみての結論としては、「およその力加減」という状態を見るには、「100段の移動平均」が適しているように見える一方で、激しい動きのレスポンスとしては、「50段の移動平均」がかなり良好であると言えます。このサンプリング周波数の条件では、筋電信号の周波数に対して、20段以下の移動平均では「平滑」にならない、という事が判りました。
