質問:
スライディング方式の離散周期変換について教えてください。リアルタイム・アプリケーションにおいて、そのアルゴリズムは心拍数や信号品質の変動をどのように処理するのでしょうか？
回答:
離散周期変換（DPT：Discrete Period Transform）のアルゴリズムは、リアルタイム・アプリケーションにおいて心拍数や信号品質の変動に適応できるように設計されています。このアルゴリズムでは、スライディング・ウィンドウという手法を使用します。つまり、新しいサンプルに対して新しいウィンドウを適用してデータを分析します。DPTのアルゴリズムでは、このスライディング・ウィンドウ方式によって時間の経過に伴う信号の変化を追跡します。また、この方式を活用することにより、分析方法を調整することが可能になります。更に、DPTアルゴリズムでは、自己相関とアンサンブル平均化の手法を組み合わせて使用します。それにより、アーチファクトの影響を低減し、S/N比を改善することができます。結果として、厳しい条件下でも精度と信頼性に優れる結果を得ることが可能になります。
はじめに
生理学的な起源から生じる生体信号は、ノイズやモーション・アーチファクトによって変動する可能性があり、多くの場合、信号自体と同じ通過帯域を共有します1。また、生体信号はほぼ定常的（準定常）なものですが、時間の経過と共に周期と振幅が変化する可能性があります2。そうした信号（データ）に対し、単純なフィルタ処理を適用すべきではありません。それ以外の方法によって情報を抽出しなければ、求められる結果は得られないでしょう。有効な方法の1つは、データに対して時間的にリンクされ、アンサンブル平均の時間枠として機能する別の信号を使用することです。アンサンブル平均は、心電図（ECG）のソースからの外部トリガを使用することで、オキシメトリ（血中の酸素濃度の測定）の信号に効果的に適用されます3。ただ、ECGのソースを利用できないケースも少なくありません。筆者らは、ECGのソースからのトリガを使うことなく信号を処理し、トリガを使用する場合と同等の結果を得ることに成功しました。
当初、筆者らは一種の自己相関処理とアンサンブル平均の処理を実行するアルゴリズムを開発していました4。しかし、時間領域でのアンサンブル平均は不要であることがすぐに判明しました。なぜなら、関連するすべての情報は周期領域のデータ自体に含まれていることがわかったからです。心拍数と血中酸素飽和度（SpO2）は、スライディング・ウィンドウ方式のDPT（以下、スライディングDPT）によって生成された結果から直接計算することができました。
筆者らの研究は、離散フーリエ変換（DFT：Discrete Fourier Transform）を見直すことから始まりました。DFTを使用すれば、信号の周波数スペクトルを生成できるだけでなく、それを利用して周期を決めることもできます5、6。この研究のもう1つの目標は、非常に高い分解能でサンプリングを行えるようにすることでした。DFTで高い分解能を得るには、数多くのデータ・サンプルを収集する必要があります。生体信号はほぼ定常的なものだと言えますが、完全に定常的なわけではありません。そのため、多数のサンプルを使用してDFTを実行すると、スペクトルが不鮮明なものになることがよくあります7。筆者らが必要としていたのは、分解能が高く、DFTを使用する場合よりも少ないサンプルで良好な結果が得られるアルゴリズムです。そのアルゴリズムによって、不定長かつリアルタイムのデータを処理することが目標でした。この目標を達成するために、スライディングDFTに似たスライディングDPTを採用することにしたのです。
開発した手法
ここからは、筆者らが開発したスライディングDPTについて詳しく説明していきます。
アルゴリズムの要件
当初の目的は、データが確率的かつ非定常的なものだとしても、その根底にある基本周期を決定することが可能なアルゴリズムを見いだすことでした。そのアルゴリズムの要件は以下のようなものでした。
- ECGやSpO2など、あらゆる生体信号の基本周期を決定できること
- 応答時間は、心拍数の周期や振幅の変化をリアルタイムに追跡できるくらい短いこと（十分に高速であること）
- 信号の中断や、過剰なノイズ／モーション・アーチファクトが生じた後に迅速に回復できること
- サンプリング・レートを決定する際の制限要因にならないよう、十分に高い演算速度が得られること
- 低消費電力であることが求められる携帯型機器のように、記憶容量が低～中程度であっても実行できること
アルゴリズムの開発
筆者らの目的は、周期を見いだせるようにすることでした。その起点として、まずはDFTの式に注目し、周波数の項を周期に置き換えることにしました。そして、DFTのように周波数をインクリメントする代わりに周期をインクリメントすることにしました。DFTでは、f0を基本波とすると、周波数を1f0、2f0、3f0、…といった具合に直線的に増加させます。それに対し、DPTでは、周期をサンプリング周期T0の倍数として直線的に増加させます。2つのアルゴリズムの式には類似点がありますが、DFTによってDPTと同じ結果を得ることはできません。これらは基本的に異なるアルゴリズムだからです。DFTとDPTは、それぞれの実装の基になる式を分析することで比較できます。サンプリング周波数をfSとすると、NポイントのDFTの周波数ビンkは、周波数fK = k×fS/N〔Hz〕に対応します。以下に示す式（1）は、サンプル・シーケンスXi... Xi + n - 1のk番目の周波数ビンのスペクトルを表します（k = 0、1、2、... N - 1）。
DFTのi番目のサンプルは、以下の式のように表されます。
DFTの基底関数は、図1に示す各高調波の前のN番目の縦座標値と同じ縦座標値を持ちます。DFTで得られたすべての高調波は、高次の高調波が低次の高調波の正確な倍数になるという性質（以下、通約性）を備えているからです。
DPTのNの項は、周期ごとに修正する必要があります。周期はそれぞれの単純な倍数ではなく、サンプル周期ごとに変化するからです（図2）。
スライディングDFTにもスライディングDPTにも、一定の数の最新サンプルを保持するための循環バッファ（または反復バッファ）が必要です。実装にあたってはこの要件を満たさなければなりません。入力されるデータが実数の場合、1つのバッファが使用されます。入力されるデータが複素数の場合には2つのバッファが必要です。DPTのi番目のサンプルは、以下の式で表されます。
ここで、RBSは循環バッファのサイズ、TLは最長周期の長さ、Tnはその時点で処理を行っている基底要素の周期です。このようにすれば、各基底周期の開始時と終了時の縦座標値が等しくなります。周期sは、サンプリングされたデータの周期をカバーするために最小周期から最大周期までの範囲内で選択します。この実装では、図2に示した複素正弦波の増分位相角である一連の基底関数を使用します。
DPTでは、基底関数が複素関数の集合で構成されます。そのほとんどには通約性がなく、サンプル周期ごとに異なるものになります。そのため、DPTの実装の難易度はやや高くなります。効率的なDPTでは、図3に示す基底位相角を使用します。本稿で紹介する実装でも、この方法を採用しました。
フェーザは、以下の式を使用することで簡単に導出できます。
ここでsは、サンプル周期刻みで選択した最小周期から最大周期までの周期の集合を表します。
アルゴリズムの実装
スライディングDPTは、IIR（Infinite Impulse Response：無限インパルス応答）フィルタの形で実装しました。信号の流れを表す図で言えば、スライディングDFTの実装と同様に、コム・フィルタの後段に共振回路が配置される形になります。コム・フィルタによってN番目のサンプルに生じる遅延により、過渡応答はN- 1番目のサンプルの長さになります。なお、心拍数に対して調整されたコム・フィルタは他の例でも一定の成果を収めています8。DPTの複素基底関数または位相角の成分のすべてが調和関係にあるわけではありません。したがって、これらの関数の端点はDFTを使用する場合のようにサンプル空間で連続関数になるとは限りません。しかし、スライディングDPTを実装すれば、基底関数がラップされて基底関数が連続したものになります。データと基底関数をスライドして相関を算出することにより、基底関数の連続性が維持されます。
スライディング・ウィンドウのアルゴリズムでは、長さがNのウィンドウを不定長のデータ配列上でスライドします。DPTの場合、実数と虚数の両方の入力データを処理できます。それに向けて、2つの循環バッファを使用します。ただ、入力に実数成分しかなく、1つの循環バッファだけが使われるケースも少なくありません。とはいえ、入力と基底関数の位相関係によっては、結果が複素数になることがあり得ます。その結果は、最大周期の長さに対応する2つのアンサンブル・バッファに格納されます。
MATLABで実現した概念実証モデル
筆者らは、MATLAB®のスクリプトを使用して式（4）を実装しました。図4に、振幅が±1、周期が45ミリ秒／79ミリ秒／175ミリ秒の正弦関数を入力として使用した場合のプロットを示します。MATLABのスクリプトでは、400ミリ秒（200周期/分）から2秒（40周期/分）までに周期を制限しています。この例では、計5000個の特定のデータ・サンプルを対象として処理を実行しました。入力したのは振幅が1の正弦波のデータなので、各周期の振幅も1になっています。図4から、この変換は高い分解能で実行されていることがわかります。
図5に示したのは、1分あたり73周期で振幅が4.5の余弦波を使用した場合の結果です。この例では、1500データ・ポイントの長さの循環バッファを使用しました。振幅の誤差は0.366%、周期の誤差は0.234%です。誤差をここまで小さく抑えられている点は注目に値します。一般的なバイオメディカル・アプリケーションでは、この程度の誤差は許容範囲内です。例えば、末梢毛細血管のSpO2の測定において、これらの誤差が重大な問題につながることはないでしょう。SpO2は、赤色分光信号と赤外線分光信号の「比の比」（RoR：Ratio of Ratios）を用いてレシオメトリックに計算されるからです9、10。
得られた結果
ここからは、開発したアルゴリズムによって得られる結果を示していきます。
スライディングDPTをパルス・オキシメトリに適用する
パルス・オキシメトリでスライディングDPTを使用する方法についてもう少し詳しく説明します。この用途で適切な結果を得るためには2つの循環配列が必要になります。1つは赤色光の履歴用に使用し、もう1つは赤外光の履歴用に使用します。スライディングDPTのアルゴリズムを実行する際、処理される周期のビンの長さを備える循環バッファの更新内容は、その周期に対応する基底関数によって順に変化していきます。このバッファの長さによって全体の分解能が決まります。また、これらのバッファを満たせるだけの十分なデータを処理すると、変換結果は安定限界に達します。その結果、入力データの変化に応じて振幅または周期が変化するだけの状態になります。本稿で紹介するデータ処理では、循環バッファによって10秒間分の最新データを保持することにしました。
未処理のデータは、アナログ・デバイセズの研究者が収集しました。データ処理には、MATLABのスクリプトによって実装したスライディングDPTのソフトウェアを使用しました。図6に示したのが結果の一例です。この図には、ある被験者から取得した未処理のデータ、1Hz～4Hzのバンドパス・フィルタで処理したデータ、合計幅が200ミリ秒の平滑化用移動平均フィルタを使用したデータをプロットしています。一方、図7に示すスペクトルは、循環バッファが満たされ、スペクトルが安定した振幅に達した後に取得したものです。DPTでは、新たなデータがサンプリングされたら、未処理のデータの変化を追跡し続けることになります。また、それに応じてスペクトルが更新されます。
SpO2の値の推定には、一般的に比の比と呼ばれている式を使用しました。AC成分の入力には、図7に示したスペクトルのプロットのピーク値を使用しています。また、DC成分の入力には図6に示したフィルタ処理を適用していない信号の平均値を使用しました。
筆者らは、スライディングDPTの性能を評価するために1つの比較を行いました。比較の対象としたのは、Masimoのオキシメータ（酸素濃度計）とSET（Signal extraction Technology）のアルゴリズムを使用して収集したSpO2と心拍数のデータです。スライディングDPTによるデータの取得には、アナログ・デバイセズの「MAX30101」を使用しました。これは、PPG（Photoplethysmogram）に対応するパルス・オキシメータ用のセンサーです。これを使用して、SpO2と心拍数のデータを同時に取得しました。図8と図9は、被験者からのデータを無作為に選択してプロットしたものです。
医療分野では、同じパラメータの値を測定する2種類の計測器を使用し、それぞれによって取得された値を比較／評価するということが行われます。一方の計測器によって正しい結果が得られると仮定し、それを基準（標準）として使用します。
稿末の参考文献11、12を執筆したBland氏とAltman氏は、2つの定量的な測定値の間の一致性に基づく分析方法を考案しました。2人は平均差について研究し、その成果として一致性の許容範囲を定める方法を生み出しました。このBland-Altman分析は、平均差の間の偏りを評価し、一致性がある範囲を推定する簡単な方法だと言えます。例として、2つの医療用計測器のテストを実施し、一方を基準として扱うケースを考えます。その場合、臨床の用途でもう一方の計測器を使用したとして、基準になる計測器と同等のレベルの結果が得られるのか否かを判断するにはどうすればよいのでしょうか。同等のレベルにあると見なすには、得られた結果が標準偏差（以下、SD）の±2倍つまり95%の範囲内に収まっている必要があります。
Bland-Altman分析を使用する場合、2つの変数の間の差を扱う統計的手法が適用されることになります。この点は、2つの変数の間の関係を特定する相関分析とは対照的です。
DPTのアルゴリズムの精度と再現性の評価も実施しました。そのために、まずMAX30101を使用して26人の健康な成人被験者のデータを収集しました。また、最新のSignal extraction Technology®13を搭載したMasimoのオキシメータでも同じ条件でデータを取得しました。比較の対象（コホート）となったのは、20歳から40歳までの15人の男性と11人の女性です。ここで、SpO2の値は性別によって多少異なることがわかっています。ある研究では、健康な成人（比較的年齢は低い）の場合、男性のSpO2は平均で97.1±1.2%、女性のSpO2は平均で98.6±1.0%という結果が得られています14。ただ、筆者らの評価では、各被験者を対象とし、2つのオキシメータによる測定結果を単純に比較しています。男女間の違いを比較することは目的としていません。
Bland-Altman分析の基準を使用した結果を図10と図11に示しました。各円は、個々の被験者の分析結果を表しています。図10を見ると、SpO2のすべての値はBland-Altman分析の基準を満たしていることがわかります。
図11で矢印によって指し示した心拍数の測定結果は、標準偏差の2倍の範囲外に位置しています。図12に、時間に対するこの被験者の心拍数のプロットを示しました。Masimoのオキシメータで測定した場合、標準偏差は1.7892です。一方、DPTアルゴリズムを実装したMAX30101で測定した場合の標準偏差は0.8935でした。どちらの測定器が正確なのかを特定するのは困難ですが、標準偏差の値は1つの目安になるでしょう。
DPTのアルゴリズムを実装したオキシメータのプロトタイプ
筆者らは、DPTのアルゴリズムを実装したオキシメータのプロトタイプも設計しました。SpO2と心拍数を測定するセンサーとしては「MAX30102」を採用しました。また、「Raspberry Pi Zero」をコンピュータ・プラットフォームとして選択し、ベアメタルOSで動作するArm®プロセッサを使用することにしました。OSとスライディングDPTのアルゴリズムは、標準的なC言語によって記述しました。図13に示したのがプロトタイプの外観です。オキシメータ全体に対しては、USB 3.0を使用して給電します。ソフトウェアを使用した制御により、2つのD/Aコンバータからリボン・ケーブルを介してTektronixのオシロスコープ「DPO4032」にデータを送信し、プロットを生成するようにしています。その後、プロットのデータはネットワーク接続を介してデスクトップ型PCに送信されます。図15に示したのが、ある被験者を対象として取得した結果です。循環バッファを満たすために10秒間経過した後、約9秒間にわたってデータを取得しました。
赤色光と赤外光のDC信号は、1次のローパス・フィルタ（IIRフィルタ）を使用して未処理の信号から抽出します。一方、AC信号は1次のハイパス・フィルタ（IIRフィルタ）を使用して抽出します（図14）。両フィルタの時定数は約1秒に設定しました。サンプル・データは、MAX30102からの割り込みをタイミング信号として100SPSのレートで取得します。赤色光と赤外光に対応する信号が、12ビット、固定小数点のデジタル・データとして取得されます。
この例では、赤色光と赤外光のAC信号をフィルタによって抽出し、それ以上の前処理を行うことなくDPTによって処理を実施しました。その結果である図16を見ると、分光信号の基本成分に対応したピークが発生しています。心拍数は横軸上のピークの位置によって決まります。比の比の式によってSpO2の値を直接計算するために、赤色光と赤外光のピークの振幅を使用しました。
考察
オキシメータによって生成された未処理の光信号には、安定した大きなDC成分と、DC成分の約1%で小さく振動するAC成分が含まれています。この振動成分は、毛細血管内の脈動活動を表しています。モーション・アーチファクトやその他のアーチファクトが存在すると、そうした信号が覆い隠されてしまう可能性があります。つまり、信号を正確に取得できなくなるということです。そのため、アーチファクトと本来の信号を分離する方法を確立することを目指して膨大な時間が費やされてきました。ただ、それらの方法の多くは非常に複雑で、実装が困難であることが明らかになっています16、17。
本稿で紹介した研究は、上記の問題を解消することを目的としたものです。DPTのアルゴリズムは少数のサンプルしか必要としません。それでも高精度の測定を実現でき、多くの課題が解消されます。周期の領域で測定を行い、各ビンをサンプル周期の分だけ離すことにより、必要な分解能が得られます。DPTによって得た周期と振幅の情報を使用することにより、時間領域に戻ることなく、心拍数とSpO2の値を直接計算することができます。
まとめ
本稿では、インクリメンタル方式のDPTのアルゴリズムを利用した周期領域の分析方法について解説しました。この方法は、周期的な生体信号を処理してスペクトル成分を求めるための効率的かつ効果的な手段だと言えます。また、周波数領域の分析機能を利用できることに加え、実装上の長所も備えています。アナログ・デバイセズのセンサーICであるMAX30101を使用すれば、DPTのアルゴリズムを実行することができます。それによって、本稿で示したとおり、医療現場で使われているMasimoのオキシメータを置き換えられるレベルの十分な精度が得られます。
謝辞
建設的な案を提示していただくと共に、本稿の校正にも協力してくださったウィスコンシン大学医用生体工学部のAmit Nimunkar博士に感謝します。また、貴重な支援と激励をいただいたEverett Smith博士、本稿のプロジェクト向けに未処理の被験者データを提供してくれたアナログ・デバイセズに感謝します。
References
- 1 Amal Jubran「Pulse Oximetry（パルス・オキシメトリ）」Critical Care（救命救急医療）、Vol. 3、No. 2、1999年2月.
- 2 Han-Wook Lee、Ju-Won Lee、Won-Geun Jun、Gun-Ki Lee「The Periodic Moving Average Filter for Removing Motion Artifacts from PPG Signals（PPG信号からモーション・アーチファクトを除去するための周期的移動平均フィルタ）」International Journal of Control, Automation, and Systems（制御、オートメーション、システムに関する国際ジャーナル）、Vol. 5、No. 6、2007年12月
- 3 Brendan Conlon、James A. Devine、James A. Dittma「 ECG Synchronized Pulse Oximeter（ECGに同期するパルス・オキシメータ）」U.S. Patent（米国特許）第4,960,126号、1990年10月
- 4 James Reuss、Dennis Bahr「Period Domain Analysis in Fetal Pulse Oximetry（胎児用のパルス・オキシメトリにおける周期領域の解析）」Proceedings of the Second Joint 24th Annual Conference and the Annual Fall Meeting of the Biomedical Engineering Society、2002年
- 5 Eric Jacobsen、Richard Lyons「The Sliding DFT（スライディングDFT）」IEEE Signal Processing Magazine（IEEE信号処理誌）、Vol. 20、No. 2、2003年3月
- 6 Eric Jacobsen、Richard Lyons.「An Update to the Sliding DFT（スライディングDFTに関する最新情報）」IEEE Signal Processing Magazine（IEEE信号処理誌）、Vol. 21、No. 1、2004年1月
- 7 Lawrence R. Rabiner、Bernard Gold「Theory and Application of Digital Signal Processing（デジタル信号処理の理論と応用）」Prentice-Hall、1975年1月
- 8 Ludvik Alkhoury、Ji-Won Choi、Chizhong Wang、Arjun Rajasekar、Sayandeep Acharya、Sean Mahoney、Barry S. Shender、Leonid Hrebien、Mose Kam「Heart-Rate Tuned Comb Filters For Processing Photoplethysmogram (PPG)Signals in Pulse Oximetry（パルス・オキシメトリで光電式容積脈波（PPG）信号を処理するために心拍数に対して調整されたコム・フィルタ）」Journal of Clinical Monitoring and Computing（臨床モニタリングおよびコンピューティング・ジャーナル）、Vol. 35、No. 4、2021年8月
- 9 「Recommended Configurations and Operating Profiles for MAX30101/MAX30102 EV Kits（MAX30101/MAX30102の評価用キット、推奨される構成と動作プロファイル）」Maxim Integrated、2018年3月
- 10 Sang-Soo Oak、Praveen Aroul「How to Design Peripheral Oxygen Saturation (SpO2) and Optical Heart Rate Monitoring (OHRM) Systems Using the AFE4403（AFE4403を使用し、末梢酸素飽和度（SpO2）と光学式心拍数を監視（OHRM）するためのシステムを設計）」Texas Instruments、2015年3月
- 11 Douglas Altman、J. Martin Bland「Measurement in Medicine: The Analysis of Method Comparison Studies（医学向けの測定：方法比較研究の分析）」Journal of the Royal Statistical Society: Series D (The Statistician)（王立統計学会誌：シリーズD（統計学者））、Vol. 32、No. 3、1983年9月
- 12 J. Martin Bland、Douglas G. Altman「Statistical Methods for Assessing Agreement Between Two Methods of Clinical Measurement（2つの臨床的測定法の間の一致度を評価するための統計的手法）」The Lancet、1986年2月
- 13 Julian M. Goldman、Michael T. Petterson、Robert J.Kopotic、Steven J. Barker「Masimo Signal Extraction Pulse Oximetry（Masimoの信号抽出パルス・オキシメトリ）」Journal of Clinical Monitoring and Computing（臨床モニタリングおよびコンピューティング・ジャーナル）、Vol. 16、No.7、2000年
- 14 Sagi Levental、Elie Picard、Francis Mimouni、Leon Joseph、Tal Y Samuel、Reuben Bromiker、Dror Mandel、Nissim Arish、Shmuel Goldberg「Sex-Linked Differencein Blood Oxygen Saturation（血中酸素飽和度の性差）」The Clinical Respiratory Journal（臨床呼吸器ジャーナル）、Vol.12、No. 5、2018年5月
- 15 Sam Koblenski「Everyday DSP for Programmers: DC and Impulsive Noise Removal（プログラマのための日常的なDSP：直流ノイズとインパルス・ノイズの除去）」2015年11月
- 16 Surekha Palreddy「Signal Processing Algorithms（信号処理のアルゴリズム）」Design of Pulse Oximeters（パルス酸素濃度計の設計）、first edition、1997年10月
- 17 Terry L. Rusch、Ravi Sankar、John E. Scharf「Signal Processing Methods for Pulse Oximetry（パルス・オキシメトリのための信号処理方法）」Computers in Biology and Medicine（生物学と医学におけるコンピュータ）、Vol. 26、No.2、1996年3月
