bandpass filterの帯域指定値

R2019aでbandpass関数を使用してフィルターをかけましたが、フィルターのかかっていないデータが返ってきました。

256HZでサンプリングされたdataに対して0.1~30Hzのバンドパスをかけたいため、

A=bandpass(data,[0.1 30],256);

とコードしましたが、フィルタリングされていないデータが返ってきました。

B=bandpass(data,[1 30],256);

と1～30Hzでフィルタリングしたデータは問題ありません。

下図は同じデータをフィルタリングしたもののパワースペクトルでAが赤線、Bが青線となっています。

サポートページにも通過周波数を指定するfpassのベクトルに関して整数値とは指定されていなかったのですが、

整数値でないといけないのか、コードのエラーなのか教えていただければ幸いです。

[Akira Agata]

bandpass関数に入力されている配列 data のサイズ（＝サンプル数）は、いくつでしょうか？

おそらく、通過帯域の低周波側 (0.1 Hz) に対してbandpass関数内部でIIRフィルタを設計するにあたり、入力信号のサンプル数が足りていないことが原因ではないかと推測します。

[Naoya]

バンドパスフィルタはローパスフィルタとハイパスフィルタの組み合わせになると思いますが、低周波側の通過帯域端が 0.1 Hz ということで、低周波側の阻止帯域幅がとても狭く特性的に厳しいフィルタ仕様となっているかもしれません。

また、bandpass関数のフィルタリングは filtfilt関数によるゼロ位相フィルタリングを行うことになりますが、このようなフィルタ特性が厳しい場合、特に信号の両端部分において大きくゆがむような時刻歴出力が得られる傾向にあります。

例えば、

>> [~,d] = bandpass(data, [0.1,30],256);

で フィルタオブジェクト d を得ることができます。

bandpass 関数内部のfiltfilt処理ではなく、 得られたフィルタオブジェクト d を使って、filter処理でより理想に近い時刻歴出力が得られるかもしれません。

>> A = filter(d, data);

[Shoumei]

この問題そのものの解ではありませんが、こういった意図しない動作をした場合の解決方法です。

該当するMATLAB関数をエディターで開く。（ソースコードが見れるのもMATLABの良いところ）

>> edit bandpass

確認したい箇所にブレークポイントを設定して、その関数を実行する。

例えば

162行目：フィルタ設計で与えられる各パラメータ

253行目：設計されたフィルタの仕様や特性を確認

125行目：設計されたフィルタでフィルタリングを行った結果を確認

といったところを確認すると、内部で処理が意図した通りに行われているか確認できます。

入力信号、出力信号が意図した通りに定義されているか時系列プロットして確認してみてはいかがでしょうか。

MathWorksさんは推奨されない使い方かもしれませんが・・・

[kazuhiko takabatake]

>>Akira Agata

dataは脳波を256Hzでサンプリングした、サンプルサイズ46080になるデータです。

>>Naoya

アドバイスいただいたフィルタ関数を得て、filter処理を行う事で改善はしましたが、まだ30Hz以上の周波数帯の成分が含まれる結果でした。（下図赤線）

青線はbandpass関数で1~30Hzに設定したもので、黄色はScott Makeig（SCCN/INC/UCSD [e-mail]）によって配布されているeegfiltというコードを比較のため使用したものになります。

ご指摘のように0.1Hzのローカット設定ではフィルタの歪みが含まれてしまうため、bandpass関数をそのまま使用することは困難なようでした。ほかのコードを使用するかフィルターを設計してかける方法も考慮してみます。

>>shoumei

アドバイスありがとうございます。あまりやったことがないのですが、editでコードの内部も一度確認してみます。