======================
(野村さんへのメモ)
・図2-5-11はネットから落としたので、そのまま使わずトレースして下さい
======================

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

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

2-5　ノイズ除去フィルタ

本節では、前節2-4で「検波+平滑」までの処理を行った筋電センサシステムのmbedマイコンNucleoF401REのプログラムをより拡張していく事例として、いくつかのノイズ除去フィルタのアルゴリズムを紹介します。

第1章(1-3)でも紹介したように、CQ出版から多くの信号処理の書籍を出版している三上直樹氏が、ちょうどインターフェース誌でNucleoF401REをオーディオ信号処理に活用した事例を連載し、そのNucleoF401REプログラムをmbedのサイトで公開しています(筆者も本稿のプログラムを同様に公開中)。

そこで本節では、まずFIRフィルタによってローパスフィルタを組み込んだ事例として、mbedのサイトにある三上氏のNucleoF401REソースをこちらのシステムで動くように埋め込んでみましたので、オープンソースのコードを自分なりにカスタマイズする参考としてみて下さい。
この他に本節では、「電源ノイズ(ハム)除去ノッチフィルタ」と、筋電など生体信号に特有の「アーティファクトのフィルタ」のアルゴリズムを実装していきます。

●高域ノイズ抑止と積分(平滑) - LPF

リスト2-5-1では、前節2-4での「検波+平滑」までの処理加えてに、上述の三上氏(CQpub0 Mikami)がmbedのサイトで公開しているコード(FIR_LPF_Direct)から、ディジタルフィルタ演算の係数などをそのまま活用しました。三上氏のサンプルコードはオーディオ帯域用で、本システムの筋電信号帯域ではほぼ「素通し」ですので、係数を変更したい場合には、CQ出版サイトで三上氏が公開しているという設計ツールを使って下さい。ここでは、FIR LPFのディジタル演算を並行して呼び出してもNucleoF401REはまだまだ余裕で仕事できている、という事を確認しています。

(リスト2-5-1)

このNucleoF401REソフトの改訂点としては、タイマ割り込みの時間間隔を実戦的な5μsecに設定し、メインループ中のカウント値を三上氏のプログラムに合わせて、同じ「0.1msec」のサンプリング周期でディジタルフィルタ演算を行いました。結果をホストPCにXBee経由でシリアル伝送する間隔は10msecとして、前節で200段の移動平均は不要と判ったので、100段までを選択するようにしました。なお、FIRフィルタの部分を「サブ FIR_LPF.hpp」として切り出して、「サブ sub.hpp」の部分はリスト2-4-1から変更ありません。

移動平均については、データをFIFOバッファに積んで、読み出しポインタをうまく操作することで、「多段のデータをいちいち隣りに移動」というループ処理を省略したのですが、こちらのFIRフィルタの処理は教科書通りに、

    float acc = 0;	累算値を初期化

    xn[0] = data;	入力A/D値をデータ[0]へ

    for (i=0; i<=order; i++) acc = acc + fc[i] * xn[i];	係数と乗算して累算

    for (i=order; i>0; i--) xn[i] = xn[i-1];	データを1つずつ隣りに移動

    return(acc);	累算値を返す

というように、見やすさを尊重してそのまま愚直に計算しています。

これを、ホストPCのMaxパッチで動作確認したのが、図2-5-1から図2-5-4です。波形のグラフは左端から入力されて右にスクロールしていくので、右ほど古い、というように眺めて下さい。

図2-5-1は、右半分は「検波オフ」(移動平均もオフ)の状態で、入力信号はバイポーラになっています。時間的に後半の左半分は「検波オン」(移動平均はオフ)にしたところで、整流された波形が確認できますが、このFIRフィルタのカットオフ周波数はずっと高いので、平滑・積分の気配はまったくありません。これ以降は全て「検波オン」としています。

(図2-5-1　プログラム[リスト2-5-1]の動作画面例(1))

図2-5-2は、右半分は「5段の移動平均」の状態で、時間的に後半の左半分は「10段の移動平均」の状態です。図2-5-3は、右端から1/3は「10段の移動平均」の状態、中央付近の1/3は「20段の移動平均」の状態、時間的に後半の左端1/3は「50段の移動平均」の状態です。ディジタルフィルタによるローパスフィルタは実質的には素通しなので、移動平均による平滑化が効いてくるのが判ります。

(図2-5-2　プログラム[リスト2-5-1]の動作画面例(2))

(図2-5-3　プログラム[リスト2-5-1]の動作画面例(3))

図2-5-4は、右端から1/3は「20段の移動平均」の状態、中央付近の1/3は「60段の移動平均」の状態、時間的に後半の左端1/3は「100段の移動平均」の状態です。細かい振動成分をカットして丸めるには、100段の移動平均は良好ですが、反応が鈍くなるというトレードオフになります。

(図2-5-4　プログラム[リスト2-5-1]の動作画面例(4))

●直流成分カット - HPF(5Hz)

ノイズ除去フィルタのアルゴリズムでハイパスフィルタを使う目的としては、
	・信号の直流成分をカットする
	・電極を貼った身体が動くことによる変動成分(アーティファクト)をカットする
という2つがあります。

前者の直流成分カットHPFについては、カットオフ周波数は5Hzもあればよく、これはJIS規格の「筋電」で規定されているカットオフです。ただし本システムでは、筋電電極からのアナログ筋電センサで、既にコンデンサによる直流成分カット(F=4.8Hz)のHPFが入っているので、ここでさらにHPFを入れる意味がなく、本稿では省略します。

筋電システムにおいて。デジタルフィルタで直流成分カットHPFを入れるとすれば、筋電センサのフロントエンドに敢えてHPFを入れずにそのままA/D変換入力して、数Hzオーダの超低周波帯域の変動成分を計測・監視したい、というような特異なケースに限られそうです。ただしこれは、アーティファクトの直撃を受けるので、信号の飽和対策などが必要になります。

◆1チャンネルから2チャンネルへの改訂

ここで、本システム(NucleoF401RE側と、ホストPC側の両方)を実験用の1チャンネル(USB治具からのアナログ(3)入力)から、次のステップのために2チャンネル(USB治具からのアナログ(3)(4)入力)に改訂しました。さらに「検波のON/OFF」実験が終ったので「常に検波ON」とするとともに、上述のFIRフィルタによるローパスフィルタをON/OFFできるようにしました。また、移動平均については刻む必然性もないので、設定バッファ最大までの間で「0段〜100段」までの数値を任意にホストPCから与えるようにしました。

NucleoF401REからホストPCに対しては、2チャンネルのA/D入力のそれぞれを、約15msecごとに返すようにしました。リスト2-5-2がそのNucleoF401REのソースコードで、ここでは以降のフィルタ処理を見やすくするために、「sub.hpp」に周辺ルーチンを移動しました。これ以前のバージョンと同じ名前ですので区別するよう注意して下さい。

図2-5-5は、このように改訂した2チャンネル版のソフトを走らせて、2チャンネル分の筋電(記録)信号をUSBオーディオ治具から供給し、リアルタイムにNucleoF401REが「検波+平滑(50段の移動平均)」、FIR_LPF処理を行っている様子です。近くに配置した筋電センサ電極によって、両方のチャンネルがほぼ同期して信号を出していること、それでも波形の様子は微妙に異なっていることから、この2チャンネルの情報をうまくリアルタイム処理して、何らかの「ジェスチャー」・「ニュアンス」・「コントロール」・「コミュニケーション」のような情報を抽出できる可能性を示しています。

(リスト2-5-2)

(図2-5-5　2チャンネル化したプログラム[リスト2-5-2]の動作画面例)

●ハム抑止 - ノッチフィルタ

上述のように本システムを2チャンネルA/D化したので、USBオーディオ治具からNucleoF401REのA/D入力に与える筋電信号として、敢えて片チャンネルの共通の筋電信号を用いて、両チャンネルに60Hz/50Hzのサイン波(電源ハムノイズのイメージ)を加えて入力し、片チャンネルだけノッチフィルタ「有り」(もう一方は「素通し」)として効果を比較検証できることになりました。

ここでは、第1章の1-3で紹介した、図1-3-5のノッチフィルタのアルゴリズムと係数をNucleoF401REに実装しました。まずは文献に対応して、図1-3-5のアルゴリズムに臨時の変数を定義した図2-5-6に対して数式と定数を確認すると、サンプリング周期ごとに、
	・入力x[n]は、図の左側の加算出力(★1)と加算され → y[n]となる
	・y[n]はサンプリングごとに2段にわたって保存され → y[n-1]と → y[n-2]となる
	・y[n-1]が定数b1と乗算され、さらにy[n-2]が定数a2と乗算された結果と加算されて → ★1になる
	・y[n-1]が定数b1と乗算された結果と、y[n-2]とが加算され
		 → さらにy[n]と加算され、その結果が定数C0と乗算された結果が出力となる
というような流れとなります。

(図2-5-6　ノッチフィルタのアルゴリズムを実装するため変数を定義)

定数については、レゾナンス特性は良好な「鋭い」方を採用するとして、
	サンプリング周期 T = 0.1msec = 0.0001sec
	サンプリング周波数 fs = 1/T = 10kHz = 10000Hz
	レゾナンス帯域幅 B0 = fs/100 = 100
	ノッチ周波数 F0 = 60Hz/50Hz　※それぞれ計算して切り替える
としました。

そこで、リスト2-5-3(Unixターミナルの画面より)のように、図2-5-6そのままに係数を計算するプログラムを作って、ノッチ周波数 F0 = 60Hz/50Hzのそれぞれに対して、4つの定数を計算しました。この両方をNucleoF401REのプログラムに定数テーブルとして持ち、ホストPCからのコマンドで切り替えることにしました。

(リスト2-5-3)

そして、NucleoF401REのプログラムとして2チャンネルの両方にノッチフィルタを入れたものが、リスト2-5-4です。ここではリスト2-5-2に対して、「sub.hpp」と「FIR_LPF」の2つは変更なく共通となっています。ホストPCのMaxプログラムからXBee経由のシリアルで設定できるパラメータとしては、
	・FIR LPFのON/OFF
	・A/D入力を持ち上げるゲイン(1倍〜15倍)
	・平滑の移動平均の段数(0段〜100段)
	・各チャンネルのノッチフィルタのON/OFF
	・ノッチフィルタ周波数として50Hz/60Hz
というように定義しています。

(リスト2-5-4)

この動作を、ホストPC側のMaxパッチとUSBオーディオ治具から受けとった筋電記録信号のスクリーンショットで確認してみましょう。図2-5-7は、両チャンネルに共通の筋電信号を与えて、まだ60Hzのノイズのレベルをゼロに絞っているので、同じ筋電信号だけです。画面がスクロールしていく半分あたりで、上チャンネルだけノッチフィルタをONにすると、画面左側の無信号レベルが下がりました。これは、筋電信号レコーディングの際に乗っていた環境ハムノイズが、ノッチフィルタで低減されていることの確認です。

(図2-5-7　ノッチフィルタを適用させた[リスト2-5-4]の動作画面例(1))

図2-5-8の画面右側1/3のあたりで、両チャンネルに加える60Hzの正弦波(ノイズ)のレベルを上げてみました。画面左側2/3あたりを見ると、ノッチフィルタを入れていない下チャンネルの信号レベルがハムノイズに押し上げられているのに対して、上チャンネルの上昇は抑えられていて、ハムノイズ抑止に効いていることがわかります。

(図2-5-8　ノッチフィルタを適用させた[リスト2-5-4]の動作画面例(2))

図2-5-9ではさらに画面中央付近から、両チャンネルに加える60Hzの正弦波のレベルを上げて、筋電信号レベルを越えるほどまでにすると、ここでもノッチフィルタを入れた上チャンネルの上昇は抑えられて、ハムノイズ抑止に効いていることがわかります。

(図2-5-9　ノッチフィルタを適用させた[リスト2-5-4]の動作画面例(3))

図2-5-10では、画面右1/3あたりで付加ノイズのレベルをちょっと下げて、今度はノッチフィルタの周波数を「60Hz」から「50Hz」に変更してみました。ノッチ周波数がずれているのですが、それでも上チャンネルの方がバックグラウンドノイズのレベルを抑えていることが判りました。

(図2-5-10　ノッチフィルタを適用させた[リスト2-5-4]の動作画面例(4))

●アーティファクトのフィルタリング

さて、本節「ノイズ除去フィルタ」の最後のトピックは、筋電をはじめとする生体信号に特有の「アーティファクト」ノイズを抑止するためのフィルタリングの信号処理的なアルゴリズムです。音声信号処理や映像信号処理におけるノイズ除去フィルタというのは、ノイズの電気的な性質が確定しているので、ピンポイントにその電気的・物理的・数学的な性質に対応したアルゴリズムを設計できます。例えば、周囲の騒音を逆相でキャンセルするノイズキャンセラーなどはその典型で、環境音を入力して、その信号をマイナスにして本来のサウンドに加算する、とルールが明確です。

これに対して、生体信号のアーティファクトというのは、偶発的な接触不良とか、身体の動きに伴って電極の接触状態が変化したりするので、ルールとしてアルゴリズム化するのはとても困難です。ここでは2つの例を紹介しておきます。

◆変動成分カット - HPF(50Hz)

アーティファクトのうち、人間の身体の振動、つまり「震え」のような成分は、およそ10数Hzあたりの周波数成分になるために、カットオフ50Hzあたりのハイパスフィルタをかけることで抑止する、というアプローチがあります。あと少しだけカットオフを上げて70Hzあたりにすることで、電源ハム(50Hz/60Hz)までカットしよう、という欲張った手法もありますが、これだと肝心の筋電情報のうち低域成分についても消してしまうので、悩ましいところです。

ハイパスフィルタ(HPF)を構成するには、例えば図2-5-11のような構成の2次のIIRフィルタである「Biquadフィルタ」を使うことができます。係数の計算方法は、有名な「Cookbook formulae for audio EQ biquad filter coefficients」(http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt)という文献があり、関連した日本語サイトや計算ソフトもたくさんあります。実装の手法は上記のノッチフィルタとほぼ同様ですので、本稿では省略しています。

(図2-5-11 Biquadフィルタのアルゴリズム例)

◆レンジ監視

アーティファクトのフィルタリングのアルゴリズム例として、ヒューリスティクス(経験則)に基づく例として、ここでは「信号レンジ監視」というアイデアを紹介します。これは筋電センサの実験でよく遭遇してきた、「電極が外れる」・「電極ケーブルが接触不良/断線する」・「他の信号線と接触して異常電圧がかかる」というようなトラブルを、リアルタイムで動く筋電システムが回避するためのものです。

その原理はとても簡単で、筋電信号のレベルが適切に設定されている状況で、入力信号レベルが「ゼロまたは最大値」に振り切れてしまっている状態を監視して、
	・設定した下限スレショルド以下、あるいは上限スレショルド以上のであるか監視・検出する
	・検出したら、入力データを無視するだけでなく、その時には前回のデータを入力データとして用いる
	・この異常状態が続くかをカウントして、カウントが設定値以上ならシステムに報告する
というようなアルゴリズムとなります。一般的な「デジタル信号処理」とは様子が違いますが、生体信号の情報処理アルゴリズムとしては、これはとても重要です。つまり、生体とともに稼働するのでストップする事は許されず、「前回のデータを入力データとして用いる」という対応をで稼働を続け、さらに「異常状態が続くいたらシステムに報告する」という動作を行うものです。

実際にNucleoF401REのプログラムに実装してみたのがリスト2-5-5です。ここでは、上述のノッチフィルタとともに、汎用(共用)部分については「sub.hpp」に移動させて、メインの「main.cpp」を今後の改訂に備えて整理しましたので、前のソースプログラムと同じモジュール名でも互換でなくなっている(「FIR_LPF.hpp」だけは同じ)ことに注意して下さい。ただし、USBオーディオ治具を経由すると、最大(飽和)出力でもフルスケールに届かず、出力をOFFにしても中点電位に戻るという回路特性上、このアルゴリズムにひっかかる異常信号が発生できないので、実用に備えてホストPCからコマンドで「ON/OFF」設定できるようにしました。

(リスト2-5-5)
