そのためには,デジタル信号処理の基本を理解しておく必要がある。前回の特集記事『Raspberry Pi3を使ったリアル波形描画(6)Qtでマルチスレッド処理をやってみる で』4msぴったりのサンプリングにこだわっていたのは,例えば,心電図などのデジタル信号処理としてブレのない(ジッターのない)信号サンプリングは必須の条件だからだ。
Raspberry Pi3 ではアナログの信号を取り込む際,多くの場合 A/D変換器は外付けで用意することになる。
実際の工作については『はじめてのAD変換(RaspberryPi3で試すアナログ・デジタル変換)』の記事を参考にしていただきたい。
問題は,A/D変換器に対して,定期的なA/D変換のタイミングを作るのは誰か?ということだ。例えば,部屋の温度が何度なんかを知りたいのならば,A/D変換のタイミングはかなり遅くて適当(例えば5分おきとか)でもいいだろう。もともとの信号の変化がそんなに早くないからだ。
一方で,例えば生体信号である心電図などを計測する際には,早いサンプリングが必要となる。なぜなら,心電図の変化,例えば,心臓がドクンとしたときの急峻な立ち上がりの部分(左図のQRS波形)はおよそ幅が10msくらいなので,それよりも早いサンプリングでないとこのQRS波形をA/D変換器で取り込んで再現することができなくなる。
QRS波形の幅が10msなのに10msでサンプリングしていたら,タイミングによってはQRS波形をすっぽり取りこぼしてしまうかもしれない。
だから心電図をA/D変換するならば最低でも4msのサンプリング間隔(250Hz)は必要と言われる。
そして,そのサンプリング間隔にブレ(ジッター)があると,本当の心電図が正確に再現できないということは容易に想像できるだろう。
したがって,A/D変換器に対して,例えば心電図の取り込みならば,正確に4ms間隔の割り込みを発生させる必要がある。4ms, 3ms, 5ms, 4ms というような感じではだめだ。ブレは 4ms ± 10μs (0.25%)ぐらいには押さえたい。
このタイムインターバルのトリガーを Raspberry Pi3 の外側で発生させるのであれば,水晶発振子で定期的な割り込みを発生させれば,相当正確なサンプリング間隔を作れる。(数十キロHz~)
ただ『Raspberry Pi3を使ったリアル波形描画(6) Qtでマルチスレッド処理をやってみる』の記事に書いたように,Raspberry Pi3自身で(ある程度)正確なサンプリング間隔を作れるのであれば,わざわざ外部でサンプリングの間隔を作る回路は必要なくなるので都合がよい。Raspberry Pi3 で4ms間隔のタイミングを作ることができれば,それをトリガーにして,外部のA/D変換器から GPIOを経由して,SPI通信やI2C通信を使って A/D変換した値を取り込む例はいろいろな資料で公開されている。
外付けのA/D変換器とアナログアンプと電極があれば,心電図を取り込むこと自体はそう難しいことではない。
ただし,デジタル信号処理の基本を理解していないと,正しい計測はできない。
デジタル信号処理の基本
そこで,この本に書いてあることの一部を紹介していこうと思う。こちらの記事『アナログ信号のサンプリング処理』も参照されたし。
まず,アナログとデジタルの違いについて。
【アナログとデジタル】
- アナログとは連続的に変化する信号を無段階に表すことを指し、デジタルとはある値を離散的な数値によって表すことを指します。組込み機器では一般に、アナログデータをセンサによって電流や電圧などに変換し、その変化をADC(アナログ-デジタル変換器)などでデジタルデータにして処理を行います。変換されたデジタルデータはデジタルのまま利用されたり、DAC(デジタル-アナログ変換器)などで再度アナログデータに変換したりします。
- アナログシステムでは使用する部品の性能に誤差があるため、温度などの周辺の環境が変わっても同じような動作をさせるためには誤差を見込んだ設計や、機器生産時の調整や校正を行う必要があります。
- これに対してデジタルシステムでは、安定した動作、高再現性、変更のしやすさなどのメリットがあり、最近の組み込みシステムではセンサからのアナログ出力をできるだけ早い段階でデジタルに変換して処理や制御に用いるケースが多くなっています。
- しかしながら、デジタルシステムでは量子化時の誤差やプロセッサの処理速度の限界などがあるため、システムのリアルタイム性や連続性、コストメリットなどを考慮すると要求仕様をアナログシステムで実現した方がよい場合もあります。
【シャノンのサンプリング定理】
- 関数 f(t) がW Hz 以上の周波数を含まないとき、その関数は(1/2W)秒間隔で抽出した標本値により完全に表現できる。
たぶん,これだけ読んでも直感的に理解できないと思うので,『ディジタル信号処理の基礎』の図(p40 図3.8)を使って説明したい。
シャノンのサンプリング定理を簡単に言えば「計測したい信号の周波数の倍以上の周波数でサンプリングしないと元の信号を正しく再現できない」ということになる。
これを理解するために左図の右上の波形を見て欲しい。
右上の実線のsin波に対して,□のポイントでサンプリングしている。(見た目,サンプリングが遅い!)元波形(実線のsin波)を再現しようとすると,点線のsin波もまったく同じ□の点を通っていることが分かる。
これがいわゆる「エリアシング」という現象で,A/D変換する前に元波形からこのような高調波を取り除いてあげなければ,元にはなかったはずの波形が「ゴースト」のように現れてしまうのだ。(音声ではビート音,画像ではモアレ)
まとめると
まとめると
- アナログ信号を正確に表現するためには、最低サンプリング周波数は原信号の最高周波数の2倍に等しいか、それ以上でなければならない。
- サンプリング周波数の1/2の周波数をナイキスト周波数と呼ぶ。
- ナイキスト周波数を超える周波数成分は標本化した際に折り返し(エリアシングとも言う)という現象を生じ、再生時に元の信号として忠実には再現されない。
ということになる。実際には,上図の左下のように元波形の周波数のちょうど倍くらいのサンプリングにすることはなく,上図の右下のようにターゲットとなる波形の早い周波数成分の数倍の周波数でサンプリングをすることが多い。
なお,それでもシャノンのサンプリング定理は生きているので,必ずサンプリング周波数の1/2以上の周波数(ナイキスト周波数)を超える周波数成分はA/D変換する前に限りなくゼロにしておく必要がある。
これはアナログフィルタで取り除くしかない。このフィルタを「アンチエリアシングフィルタ」という。
デジタル信号処理を行うにあたっては,シャノンのサンプリング定理は必ず知っていないといけないし,A/D変換器の前にアナログのアンチエリアシングフィルタは必要だ。
それを知らないでいると 元波形にはなかった「ゴースト波形」に悩まされることになる。ソフトウェアエンジニアだからといってデジタル信号処理やアナログ回路のことは知らなくてもいい訳ではなく,何をやりたいかによっては,広い範囲の知識や原理原則を知る必要がある。(手段そのものの理解ではなく,目的と目的を実現する原理と,目的を実現するための方法を理解しないと,原理から外れた変更をしてしまい失敗したりする)
シャノンのサンプリング定理やアンチエリアシングフィルタのことは,ハード屋さんが考えることだと決めつけていると,ソフトウェア技術者がよかれと思ってサンプリング周波数を早くしたり,遅くしたりしたときに,アンチエリアシングフィルタとの相性が合わない事態になり,そのことに気が付かず,「ゴースト波形」が出て「なぜ?」と悩むことが起こる。
デジタル信号を扱うのならばデジタル信号の基礎は,ハード屋さんも,ソフト屋さんも知っている必要がある。
なお,本当に人間の心電図を計測するならば,電気安全対策は必ず理解して実施しておかないと危ないのでこちらの記事【トランジスタ技術2006年01月号】 投稿記事『心電計の製作』をよく読んでおくべし。
デジタル信号処理の流れ
『ディジタル信号処理の基礎』の図(p70 図3.55) を引用してデジタル信号処理の流れを説明する。
【デジタル信号処理の流れ】
- まず,低域通過フィルタ(「アンチエリアシングフィルタ」)にて,サンプリング周波数の1/2以上の周波数をカットする。
- A/D変換器で A/D変換する間アナログ信号をサンプル・ホールドする。(最近のA/D変換器はサンプル・ホールドも込みでやってくれる場合が多い)
- A/D変換したデータをCPU(上図ではDSP)に取り込む。(Raspberry Pi3 では,GPIOを経由して,SPI通信やI2C通信で取り込むのが簡単)
- 取り込んだデジタルデータに対して,デジタルフィルタなどの信号処理を行う。
- 加工が済んだら,デジタルデータをD/A変換する。(加工データを外部にアナログ出力する場合)
- D/A変換したアナログ信号に対して,低域通過フィルタをかける。
4の,デジタルフィルタや信号処理のところで,何をするのかは,何の信号を取扱い,その信号で何をしたいのかによって,いろいろ変わる。
例えば,心電図を取り扱うのであれば,次のようなことをソフトウェアで実施することになる。(【トランジスタ技術2006年01月号】 投稿記事『心電計の製作』を参照のこと)
- 目的の周波数成分(エリアシングフィルタよりも低い周波数成分)だけを取り出すために,デジタルフィルタでハイカット(低域通過)フィルタをかける。(ノイズ除去フィルタ)
- 電源ノイズ(関東:50Hz,関西:60Hz)だけを取り除く,フィルタ(電源ノイズ除去フィルタ)をかける。
- 基線のゆっくりした動揺を取り除くための ローカット(高域通過)フィルタをかける。
- 心拍数を計測するためには,心電図波形のQRS波形をカウントする必要がある。QRS波形をカウントするには,QRSを強調するようなデジタルフィルタをかけて,その上でスレッショルドレベルを設定して,スレッショルドを超えたらQRS波形としてカウントする。
なお,上記のハイカットフィルタやローカットフィルタ(ドリフトフリーフィルタ)をかけたことにより,もとの心電図波形が大きくひずんでしまっては,元も子もないので,目的に合った「有限インパルス応答」(FIR)フィルタを適用することが重要となる。
デジタルフィルタとリングバッファ
A/D変換した信号に継続的にデジタルフィルタをかける場合,データを格納する領域は有限なので,データを格納するためにリングバッファを用意する。『ディジタル信号処理の基礎』の図(p103 図4.33)
図はリング上になっているが,メモリ領域はリング上にはなっていないので,データの格納ポイントをインクリメントしていってバッファの範囲を超えたらポインタを元に戻すという処理が必要になる。
リングバッファ処理やデジタルフィルタはサンプリング毎に実施しなければいけないため,いくら今どきのCPUの処理スピードが早くなったとはいえ,4ms以下で毎回やらなければいけない処理はできるだけ軽くしておいた方がよい。
蛇足だが,整数と浮動小数点の違いを意識せずに,処理時間があまりない中で浮動小数点の演算をバンバン気にせずにやってしまうソフトウェアの新人技術者をよく見かける。CやC++で書いたソースコードが機械語に翻訳されたコードの量を見たことがないのだろう。浮動小数点の演算を専用のプロセッサを使わずに行うとものすごく機械語のコードが増えるので,それを知っていると,できるだけ整数演算で済ませられないかといつも考えるようになる。
さて,リングバッファを用意する際には,バッファサイズは2のn乗であると都合がよい。
リングバッファのポインタは,バッファの上限を超えたら,バッファ数ぶん引かないとポイントが誤った位置を指し示してしまう。バッファサイズは2のn乗であれば,簡単な計算で知りたいポインタの位置を知ることができる。
そのことを説明しよう。
例えばバッファサイズが8で、ポインタが現在 5 で、+4 の位置を知りたいなら
i = (( i + 4 ) & 0x07 ) を計算すればよい。 答え (5 + 4)- 8 = 1
ポインタが現在 5 で、-6 の位置を知りたいなら
i = (( i - 6 ) & 0x07 )を計算すればよい。 答え (5 - 6) + 8 = 7
この計算ならば処理が高速で,意味さえ理解できれいればコーディングミスが少ない。
欲しいビットフィールドを1にしてANDを取る意味を考えてみよう。
( Index + i ) & 0x07 で 5 + 4 → 1 になるしくみ
10進数の5は2進数なら 0101
10進数の4は2進数なら 0100
足すと 9 1001
10進数の7(8-1)は2進数なら 0111
(1001 を 0111=7でマスクする。すなわちANDを取る)
0001 = 10進数の1 になる。
( Index + i ) & 0x07 で 5 - 6 → 7 になるしくみ
10進数の5は2進数なら 0101
10進数の6は2進数なら 0110
引くと -1 1111
10進数の7(8-1)は2進数なら 0111
(1111 を 0111=7でマスクする。すなわちANDを取る)
0111 = 10進数の7 になる。
このように,バッファサイズが 2のn乗ならば,現在のポインタに対して 欲しい相対位置を足したり引いたりして,バッファサイズ-1 と AND を取ると 正しい絶対位置を計算できる。
これを,バッファサイズが2のn乗でないときと,2のn乗のときのC言語のコードの違いは次のようになる。
signed short TempIndex1:
unsigned short TempIndex2:
InputBuffer[index] = InputData;
// バッファサイズが2のn乗でない=13などのとき
TempIndex1 = (signed short)(index - 12);
if ( TempIndex1 >= BUFFER_SIZE ) (
TempIndex1 = TempIndex1 - BUFFER_SIZE;
)
else if ( TempIndex1 < 0 ) (
TempIndex1 = TempIndex1 + BUFFER_SIZE;
)
// バッファサイズが2のn乗=8のときの例
TempIndex2 = ( (signed short)Index + 4 ) & 0x07;
TempIndex2 = ( (signed short)Index - 6 ) & 0x07;
ちょっとした違いなのだが,下の方がずっといいと思うのは自分だけだろうか?
この計算をマクロに登録しておくと,すごくスッキリしたコードになる。
この計算をマクロに登録しておくと,すごくスッキリしたコードになる。
デジタルフィルタの効果を視覚的に確認する
今度は,Interface誌 2003年1月号に投稿した記事『オブジェクト思考を使ったリアルタイム信号計測システムの開発』の記事から引用する。
左図は 1Hz のsin波に10Hzのノイズに見立てた sin波を重畳させた波形だ。
デジタルフィルタを作成して 10Hz のノイズだけを取り除いてみる。
このとき Filter A と Filter B を試してみて,その結果を画面上で確認してみたのが次の図になる。
デジタルフィルタがいいのは,こういった「試しにやってみる」の効果を確認するのが簡単でできることだ。
左図の左側がノイズを重畳させた元波形とカットオフ周波数 8.4Hz のハイカットフィルタを通した結果の波形(赤色)で,右波形が元波形とカットオフ周波数4.3Hzのハイカットフィルタを通した波形(グリーン)だ。
ちなみに,カットオフ周波数(遮断周波数)とは,その周波数を越えると(あるいは下回ると)回路の利得が通常値の 3 dB 低下する値のことである。(Wikipediaの解説)
ようするに欲しい周波数帯域におけるsin波の振幅を1としたときに,ハイカットフィルタをかけると,周波数高域になるにしたがって振幅が小さくなってくる。(フィルタの特性によっては単調に減少するとは限らない)
そのときに,振幅1に対して-3dB(およそ 0.73)になる周波数をカットオフ周波数という。
紫の波形が Filter A(カットオフ周波数 8.4Hz)をかけた結果で,10Hzのノイズは低減しているものの取り切れていない。
一方,黄色の波形は Filter B(カットオフ周波数 4.3Hz)をかけた結果で,10Hzのノイズは奇麗に取れて 1Hzの sin波だけが残っている。
左図は,5Hzの sin波を Filter A と Filter B に通したときに,どれくらい時間遅れが発生しているかを見たものだ。
元波形に対して 振幅が減衰し,波形のピーク点が僅か左にずれていることが分かる。
位相ひずみがないフィルタ(FIRフィルタ)をかけると,時間遅れが発生する。ドリフトフリーフィルタ(高域通過フィルタ)などでは,この時間遅れが数百msec にもなる場合があるので,どれくらい遅れているのかが分かるようにしておくもとも重要だ。
次回は,Interface誌 2003年1月号に投稿した記事『オブジェクト思考を使ったリアルタイム信号計測システムの開発』をベースに,波形の入力を A/D変換器や サンプルデータやファイルデータに 簡単に切り替えるオブジェクト指向設計を使った事例を紹介する。
デジタル信号処理のポイントは次のようなことになる。
デジタル信号処理のポイントは次のようなことになる。
- 入力したい信号の最高周波数を調べる。
- シャノンのサンプリング定理を理解し,サンプリング周波数を定め,アンチエリアシングフィルタにて,サンプリング周波数の1/2以上の周波数成分をアナログフィルタで取り除く。
- 正確なサンプリング間隔を作り,デジタルデータを取り込む。
- デジタルデータをバッファリングして,デジタルフィルタをかける。
- その他,二次処理を行う。
そして,デジタル信号処理のシステムの保守性や再利用性を高め,新しい機能を安全に追加できるようにするためには,ソフトウェアのアーキテクチャを工夫する必要がある。
C言語でゴリゴリ書いてしまい「動くソフトウェア」を作ることはできる。しかし,それは,『ソフトウェアアーキテクチャとは何かを考える』の記事で書いたように,丸太組工法で高級マンションを建てようとしているように思える。
どのようにスマートに設計するか,できるのかを次回説明しようと思う。