ミリ波センサの呼吸信号解析：BPF・EMAスムージング・カルマンフィルタによる3段階ノイズ除去
はじめに
ミリ波センサを用いた非接触呼吸モニタリングでは、取得した信号に様々なノイズが含まれるため、安定した呼吸数の抽出が課題となります。本記事では、ノイズの多い生体信号から安定した呼吸数を抽出するための3段階の信号処理アプローチについて解説します。
BPF（バンドパスフィルタ）、EMA（指数移動平均）スムージング、カルマンフィルタを組み合わせることで、ロバストな呼吸信号解析を実現する手法を紹介します。
本記事で説明する内容
本記事では、ミリ波センサから取得した呼吸信号の段階的なノイズ除去手法について説明します。BPFからEMAスムージング、そしてカルマンフィルタへと続く3段階の処理により、それぞれ異なる特性のノイズに対処します。各処理の役割と効果を可視化しながら、実装例とパラメータ調整のポイントについても解説します。
各信号処理手法の原理と実装
BPF（バンドパスフィルタ）処理
目的
呼吸周波数帯域（0.1-1.0Hz）以外のノイズを除去します。
原理
人間の呼吸は通常6-60回/分（0.1-1.0Hz）の範囲にあります。この帯域外の信号はノイズとして除去することで、呼吸信号を抽出します。バターワースフィルタを使用することで、通過帯域では平坦な特性を保ちながら、遮断帯域では急峻な減衰特性を実現できます。
def _initialize_filters(self):
"""フィルタ係数の初期化"""
nyquist_freq = self.frame_freq / 2
# 0.1Hz低域通過フィルタ（DC成分抽出用）
lpf_cutoff = 0.1 / nyquist_freq
if lpf_cutoff >= 1.0:
lpf_cutoff = 0.95
self.lpf_b, self.lpf_a = butter(2, lpf_cutoff, btype='low')
# 0.1-1.0Hzバンドパスフィルタ（呼吸帯域）
max_freq = min(1.0, nyquist_freq * 0.9)
low_freq = 0.1 / nyquist_freq
high_freq = max_freq / nyquist_freq
if high_freq >= 1.0:
high_freq = 0.95
if low_freq >= high_freq:
low_freq = high_freq * 0.1
self.bpf_b, self.bpf_a = butter(4, [low_freq, high_freq], btype='band')
実装のポイント
カットオフ周波数の設定では、ナイキスト周波数を考慮する必要があります。サンプリング周波数の半分がナイキスト周波数となり、これを超える周波数は正しく表現できません。フィルタ次数は4次としており、これはリップルと遮断特性のバランスを考慮した選択です。高次にするほど急峻な遮断特性が得られますが、計算負荷と位相歪みが増加するため、実用的な範囲で設定しています。
EMA（指数移動平均）スムージング処理
目的
高周波ノイズの平滑化と信号の連続性向上を図ります。
原理
EMAは過去の値に指数的に減衰する重みを付けて平均化する手法です。最新の値により大きな重みを与えながら、過去の値も考慮することで、急激な変化を抑制しつつ信号のトレンドを保持できます。時定数τによってスムージングの強度を調整でき、τが大きいほど平滑化効果が強くなります。
def apply_ema_smoothing(self, signal_value):
"""
指数移動平均（EMA）スムージングを適用（リアルタイム版）
Parameters:
signal_value: float - 現在の信号値
Returns:
float - スムージング後の値
"""
if not self.ema_initialized:
# 初回は入力値をそのまま使用
self.ema_previous_value = signal_value
self.ema_initialized = True
return signal_value
# EMA計算: EMA[n] = α × X[n] + (1-α) × EMA[n-1]
ema_value = self.ema_alpha * signal_value + (1.0 - self.ema_alpha) * self.ema_previous_value
self.ema_previous_value = ema_value
return ema_value
実装のポイント
リアルタイム処理では、逐次的に新しいデータが入力されるたびに計算を行います。スムージング係数αは、α = 1 - exp(-dt/τ)で計算され、サンプリング間隔dtと時定数τの関係を適切に反映します。本実装では時定数を1.0秒に設定しており、これは呼吸信号の変化に追従しながら、細かなノイズを効果的に除去できる値として選定しています。
カルマンフィルタ処理
目的
呼吸の周期的特性を考慮した最適推定を行います。
原理
カルマンフィルタは、システムの状態空間モデルに基づいて、ノイズを含む観測値から最適な状態推定を行う手法です。呼吸信号を正弦波としてモデル化し、位相、角周波数、振幅を状態変数として扱います。予測ステップと更新ステップを繰り返すことで、観測ノイズとプロセスノイズの両方を考慮した最適推定が可能になります。
class KalmanBreathingFilter:
"""呼吸信号用カルマンフィルタ"""
def __init__(self, dt=0.05, initial_frequency=0.3):
self.dt = dt
self.initialized = False
# 状態ベクトル: [位相, 角周波数, 振幅]
self.state_dim = 3
self.x = np.zeros(self.state_dim)
self.P = np.eye(self.state_dim)
# 初期値設定
self.x[0] = 0.0 # 初期位相
self.x[1] = 2.0 * np.pi * initial_frequency # 初期角周波数
self.x[2] = 1.0 # 初期振幅
# プロセスノイズ共分散行列 Q
self.Q = np.diag([0.001, 0.0001, 0.001])
# 観測ノイズ分散 R
self.R = 0.05
状態空間モデル
状態遷移モデルでは、位相が角周波数に従って時間とともに進行し、角周波数と振幅は基本的に一定と仮定しますが、わずかなプロセスノイズを許容します。観測モデルは h(x) = A × sin(φ) となり、非線形な関係を表現しています。
実装のポイント
非線形観測モデルに対応するため、拡張カルマンフィルタ（EKF）として実装しています。観測ヤコビアンを計算し、線形化を行うことで標準的なカルマンフィルタの枠組みで処理できます。プロセスノイズQと観測ノイズRのバランスが重要で、Qが大きいと状態の変化に素早く追従しますが推定が不安定になり、Rが大きいと観測値への信頼度が下がり推定が保守的になります。
処理結果と効果
各処理段階での効果
白点線がBPF処理後の結果です。緑線がEMAスムージング、青線がカルマンフィルタの処理結果です。
緑線と青線が遅れ気味です。呼吸信号というベースが正弦波の信号であり、今回の信号処理と相性がいいと思います。
BPF処理により、呼吸帯域外のノイズが大幅に削減されます。DC成分や体動による低周波ノイズ、測定系由来の高周波ノイズが除去され、0.1-1.0Hzの呼吸信号のみが抽出されます。この段階で呼吸パターンの大まかな形状が明確になります。
EMAスムージングは、BPF処理後に残る細かなノイズを平滑化します。時定数1秒の設定により、呼吸の基本的なパターンを保持しながら、測定ノイズによる急激な変化を抑制します。信号の連続性が向上し、後段のカルマンフィルタ処理がより安定して動作する基盤を作ります。
カルマンフィルタ処理では、呼吸の周期性を積極的に活用した推定を行います。単純な平滑化では失われがちな呼吸の位相情報を保持しながら、ノイズに埋もれた呼吸信号を抽出します。位相、周波数、振幅を同時に推定することで、呼吸数だけでなく呼吸パターンの詳細な特徴も把握できます。
まとめ
本記事では、ミリ波センサから取得した呼吸信号に対して、BPF、EMAスムージング、カルマンフィルタを組み合わせた3段階のノイズ除去手法を紹介しました。各手法はそれぞれ異なる特性のノイズに対処し、組み合わせることで相乗効果を発揮します。
BPFが周波数領域でのノイズ除去を担当し、EMAスムージングが時間領域での平滑化を行い、カルマンフィルタが信号のモデルベース推定を実現するという役割分担により、ロバストな呼吸信号解析が可能になります。
Discussion