NIRSの信号処理について考えてみる〜⑤多重解像度解析〜
多重解像度解析
前回の記事では，ウェーブレット変換によって手元にある波の時間情報と周波数情報をいい感じに得ることができるということを書いた。スカログラムを出力することで「どんな波がいつ」発生したのかを知ることができるが，個々の波の形までは知ることができなかった。
この「個々の波の形」，すなわち，ある周波数帯域に含まれる波形を得る方法の1つに「多重解像度解析（MRA: Multi-resolution Analysis）」というものがある。多重解像度解析では，ある波が様々な種類の波で構成されているとして，その波を個々の波に分解するが，分解した波を足し合わせればまた元の波に戻るようにする。イメージとしては，元の波に対して遮断周波数を変えながらフィルタをかけることを繰り返して様々な周波数帯域の波を得るというようなものだろうか。
厳密には「ウェーブレット関数」や「スケーリング関数」など，考慮すべき要素が多々あるのだが，本記事ではそれらに触れず，多重解像度解析によってどのようなことができるのかを実演によって説明する。
NIRS データに対する適用
応答成分を約 0.025 Hz としたサンプルデータ（図1）に対して，MATLAB で MRA を行う。
MRA ではウェーブレット変換を使用する。マザーウェーブレットの選択に加えて，分解するレベルを設定する必要がある。今回は4次の symlet （「マザーウェーブレットの比較」を参照）でレベル12まで分解する。図2は分解した結果である。
なお，本記事の方法では，各波形が含まれる周波数帯域は以下のように算出される（図3）。
MRA の結果から，D11 が応答成分を含む波であると考えられる。確認のために，D11 の波を抽出し，周波数成分を調べた。抽出された波が図3に示した範囲の周波数の波を含んでいることがわかる（図4）。
元の波には応答成分以外の高周波・低周波ノイズも含まれているが，応答成分の周波数が特定または推測可能である場合はそれが含まれる周波数帯域の階層の波を抽出することで，結果として前の記事で紹介したフィルタリングと同等の処理がなされたことになる。
なお，MRA の適用と結果の出力，および特定の階層の波の抽出は以下のようなコードで実現できる。
x = data;
wv = 'sym4';
level = 12;
mra = modwtmra(modwt(x, wv, level), wv);
helperMRAPlot(x, mra, time, 'wavelet', 'Wavelet MRA')
%D11の抽出
ext = mra(11, :)';結果の出力について補足
MAR の結果の出力には，WaveletToolBox に含まれる "helperMRAPlot.m" を使用する。デフォルトでは波の色がデフォルトの第1カラーであり，また階層数が多くなると図が1ページに収まらない（スクロールが必要になり，特に図の保存が手間になると思われる）。
波の色を変えたい場合は，helperMRAPlot.m に記述された"plot(ax(kk),time,mra(kk,:));" の部分（おそらく31および41行目）に色のオプションを追加すればよい。また，"% Resize each subplot to almost touch the one above." 以下の "Height = 0.07"（おそらく92行目）を小さくすることで出力される図を1ページで表示できるようになる。
マザーウェーブレットの違いの影響
前の記事で述べたとおり，マザーウェーブレットの選択おける明確な基準はないが，マザーウェーブレットの違いで結果が大きく変わることはない。上の MRA について，マザーウェーブレットを変えて実行した結果について，D11 を例にして比較すると，結果に大きな差はないことがわかる（図5）。
引用文献
[1] MathWorks ヘルプセンター Practical Introduction to Multiresolution Analysis https://jp.mathworks.com/help/wavelet/ug/practical-introduction-to-multiresolution-analysis.html (2023年3月22日)