見出し画像

NIRSの信号処理について考えてみる〜②HPFとLPF〜

N.

①でdeoxy-Hbも見ろやと言っておきながら，ここからはoxy-Hbだけ登場します（だって作図が面倒なんだもの）。まあ信号処理の方法に注目しているわけだから問題ないんだけどね。

（今回のHPFとLPFの設計・適用には小野先生のテキスト[1]を参考にしています）

サンプルデータ

細かい説明は抜きにして，早速本題へ。

心理学実験だとブロックデザインという一定の周期で刺激を呈示するデザインを採用することが多いから，その刺激に対する脳の応答も周期的なものになる。

だから計測時間と刺激呈示の時間間隔から脳応答が大体何Hzの成分になるのか予想できるわけなんだよね。

ここでは，実際に「Rest (20秒) → Task (20秒)」という流れを9回繰り返したブロックデザインで得られたデータに，ノイズ除去の結果をわかりやすくするためにあえて線形トレンドと0.1 Hz (振幅0.02)の正弦波を加えたものをサンプルとして扱ってみる（図1）。

ここからはMATLABのコードも載せていくので一緒にやってみてもいいかもしれない。それぞれ手元にあるデータは違うと思うけど，とりあえず「いじるデータ（oxy-Hb）」，「時間情報」，「サンプリング周波数」の3つを用意しておけば同じようにできるはず。

data = readmatrix('sampleData.csv'); %サンプルデータ

time = data(:, 1); %時間情報がデータの1列目にあるので抽出
oxy = data(:, 2); %サンプルデータが2列目にあるので抽出
fs = 1/time(2, :); %サンプリング周波数(このデータでは55.556 Hz)
図1　サンプルデータ
%図1の作り方

plot(time, oxy, 'r');
axis tight
ylim([-0.2 0.4])
xlabel('Time (s)'); 
ylabel('\DeltaHb concentration (mM mm)');
title('fNIRS data');


%（なお，図1の縦波線は刺激呈示のタイミングを示すが，自作関数を作成して表示しているため省略）


360秒の中で9回の周期的な反応が得られるから，脳の応答成分は0.025 Hzだと予想される。

だから0.025 Hzよりも低い周波数成分と高い周波数成分を除去してあげればいいわけ。

このデータに含まれる周波数成分を調べてみると，付与したノイズと脳の応答成分がきちんと含まれていることがわかる（図2）。

（ここではフーリエ変換をしているんだけど，これについては次の記事で少しだけ詳しく説明するつもり）

図2　サンプルデータの周波数成分
%図2の作り方

x = oxy';

L = length(x);
win = hann(L)'; %ハン窓を使用

xWin = x.*win; %窓関数の畳み込み

acf = 1/(sum(win)/L); %Amplitude Correction Factorの算出

Y = fft(xWin); %フーリエ変換
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = fs*(0:(L/2))/L;

plot(f, P1*acf, 'r', 'LineWidth', 2)
title('Single-Sided Amplitude Spectrum')
xlabel('f (Hz)')
ylabel('|P1(f)|')
xlim([0 0.15])
ylim([0 0.07])


ハイパスフィルタで低周波成分を除去

ではまず，HighをPassするハイパスフィルタ (HPF)を使って応答成分よりも下の周波数成分を除去してみよう。

MATLABだと色々な方法があるけど，ここでは"fir1"という関数を使用してフィルタを作成する。次数3200，遮断周波数0.01 Hzでやってみる（何故フィルタの次数がこんなに高いのかについては後述）。

HPF = fir1(3200, 0.01/(fs/2), 'high');
oxy_H = filtfilt(HPF, 1, oxy);

freqz(HPF, 1) %フィルタの周波数応答を確認（図5を参照）

作成したフィルタを関数"filtfilt"で適用すると，ドリフト補正ができていることがわかる（図3）。ただし，元のデータに含まれていた約0.005 Hz未満の成分の振幅が大きかったため，これを減衰したものの完全に取り除けたわけではない。

これを取り除く方法もないわけではないけど，今回はとりあえずこのまま行く。

図3　HPF適用後

そもそもフィルタを通すと遮断周波数以降の成分が完全に消えるわけではなくて，あくまで「減衰」するということには注意が必要。

また，遮断周波数の周辺にもフィルタの影響が出ることにも注意（図4）。

図4　フィルタの理想と現実
（https://synapse.kyoto/glossary/cutoff_frequency/page001.html [2]を参考に作成）

で，次数についてだけど，今回のデータで十分なドリフト補正ができるものを探していたら3200くらいだったんだよね。

"fir1"で設計したフィルタは遮断周波数での正規化ゲインが–6dB（つまり遮断周波数の部分の振幅を0.5倍にする）となるらしいんだけど [3]，設定した遮断周波数が小さいからか–6dBの減衰が得られたのが3200次だったというわけで（図5）。

図5　設計したHPFの遮断周波数での応答


ローパスフィルタで高周波成分を除去

では次に，LowをPassするローパスフィルタ (LPF)を使って応答成分よりも上の周波数成分を除去してみよう。

やり方はHPFとほとんど同じ。関数"fir1"でフィルタを作成して，HPFを適用したデータに適用する。次数500，遮断周波数0.07 Hzでやってみる。

LPF = fir1(500, 0.07/(fs/2));
oxy_HL = filtfilt(LPF, 1, oxy_H);

作成したフィルタを関数"filtfilt"で適用すると，高周波成分がある程度除去ができていることがわかる（図6）。

図6　HPF + LPF適用後

実際にはこの後でブロックごとのデータを切り出してベースライン補正，加算平均，標準化という流れになると思うけど，そのあたりは小野先生のテキスト[1]などをご参照くださいませませ（←ミスじゃないよ）。


線形トレンドの除去について

さてさて，上のLPF適用後の波形を見てもわかるとおり，線形トレンドが除去しきれていない問題が残る。

そもそも応答成分の周波数が0.025 Hzと小さく，その上トレンド成分がさらに小さい上に応答成分に近いということもあって，HPFでは限界がある気がする。

MATLABには線形トレンドを除去する"detrend"という関数があって，これを使うともっと簡単に補正ができる（ただし，低周波成分を減衰させるのではなく「線形」トレンドを除去するという点には注意）。

線形トレンド以外の低周波成分もあるということを考えるとdetrendとHPFを併用してもいいと思うけど，とりあえずここではHPFの代わりにdetrendを使ってみる。

ここではdetrendの第3引数（"on"）でブレークポイントを指定しているけど[4]，これまでの図の縦波線同様に刺激が呈示されたタイミングがこれに該当する。

oxy_dt = detrend(oxy, 'linear', on) %第3引数はなくてもOK

はい。これだけ。結果を見ると，線形トレンドが消えていることがわかる（図7）。これHPFを使った時よりも理想的な感じでは？

図7　detrend適用後

これに先ほどのLPFを適用するとこんな感じ（図8）。やっぱりこっちの方がいい気がする。

図8　detrend + LPF適用後


今回はこんな感じかな。お疲れさまでした。

引用・参考文献

[1] 小野弓絵 (2018) MATLABで学ぶ生体信号処理，コロナ社，pp112-113
[2] しなぷすのハード制作記　「カットオフ周波数の解説」https://synapse.kyoto/glossary/cutoff_frequency/page001.html (2023年3月13日)
[3] MathWorksヘルプセンター fir1 https://jp.mathworks.com/help/signal/ref/fir1.html (2023年3月13日)
[4] MathWorksヘルプセンター detrend https://jp.mathworks.com/help/matlab/ref/detrend.html (2023年3月13日)

この記事が気に入ったらサポートをしてみませんか？