🐳

レーダーによるバイタルサイン周波数推定 - 静的クラッター除去から呼吸検出まで

に公開
信号処理
センサー
iq
バイタル
呼吸
tech

背景

ミリ波センサーを用いた非接触バイタルサイン検出において、生体信号の抽出は複数の信号処理段階を経て実現されます。レーダーから取得される生のIQデータには、対象となる呼吸や心拍による微小な位相変化の他に、静的な反射物からの強力な信号や各種ノイズが混在しています。これらの不要な成分を除去し、バイタルサインの周波数成分を正確に推定することが、実用的な生体監視システムの構築において重要な課題となります。

従来のバイタルサイン検出手法では、接触式センサーや光学センサーが主流でしたが、60GHzレーダーは非接触での測定が可能であり、プライバシーを保護しながら連続的な監視を実現できます。しかし、レーダー信号には静的クラッターや位相の不連続性といった特有の問題があり、これらを適切に処理することで初めて生体信号の検出が可能になります。

この記事で説明していること

本記事では、ミリ波センサーから取得したIQデータを用いて、実際にバイタルサイン（主に呼吸）の周波数を推定するまでの完全な信号処理フローを解説します。各処理段階における理論的背景と実装コードを組み合わせることで、読者が実際にシステムを構築できるレベルの詳細な説明を提供します。

  • IQデータ取得: 前記事で解説した基本概念の応用
  • 静的クラッター除去: 低域通過フィルタによるDC成分抽出と減算処理
  • 位相アンラップ処理: 位相の不連続性を補正する連続化手法
  • 呼吸領域バンドパスフィルタ: 0.1-1.0Hz帯域での生体信号抽出
  • STFT法による周波数推定: 短時間フーリエ変換を用いた呼吸周波数の特定

これらの処理により、最終的には呼吸回数（BPM: Breaths Per Minute）として実用的な情報を取得できるシステムを構築します。

全体処理フロー

IQデータ取得

IQデータの取得については、前記事「レーダーで始めるバイタルサイン検出 - IQデータの基本理解」で詳細に解説しています。ここでは、そこで説明した内容を前提として、実際の信号処理に進みます。

バイタルサイン検出では、特定の距離ビン（対象物からの反射が最も強い位置）でのIQ値を時系列で取得し、その位相変化を解析します。添付コードでは以下のように実装されています：

# 対象距離ポイントのIQ値を抽出
target_iq = averaged_sweep[self.target_position]

# バッファに追加（時系列データとして蓄積）
self.iq_buffer.append(target_iq)  # 解析用（短期）
self.display_buffer.append(target_iq)  # 表示用（長期）

静的クラッター除去

原理の説明

静的クラッター除去は、バイタルサイン検出において最も重要な前処理の一つです。ミリ波センサーが受信する信号には、壁、家具、床などの静的な反射物からの強力な反射波が含まれています。これらの静的成分は、呼吸による微小な位相変化に比べて数百倍から数千倍強く、そのまま解析すると生体信号が埋もれてしまいます。

静的クラッターの除去は、低域通過フィルタを用いてDC（直流）成分を抽出し、原信号から減算することで実現されます。呼吸による位相変化は0.1-1.0Hz程度の比較的低周波成分ですが、静的反射による位相は時間的にほぼ一定であるため、0.1Hz以下の成分として現れます。

実際の除去コード

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')

def _perform_vital_analysis(self):
    """バイタルサイン解析を実行"""
    # バッファを配列に変換
    iq_array = np.array(list(self.iq_buffer))
    
    # 0.1Hz低域通過フィルタでDC成分を抽出
    lpfout = signal.lfilter(self.lpf_b, self.lpf_a, iq_array)
    
    # 静的成分除去
    sig_dcoff = iq_array - lpfout

この処理により、静的な反射成分が除去され、時間的に変化する成分（主に呼吸による位相変化）が強調されます。

位相アンラップ処理

原理の説明

位相アンラップ処理は、複素数信号から抽出した位相の不連続性を補正する重要な処理です。np.angle()関数で計算される位相は-π〜π（-180°〜180°）の範囲で周期的に変化するため、連続的な位相変化が不連続な値として観測される問題があります。

例えば、位相がπ-0.1からπ+0.1に変化する場合、実際の変化量は0.2ラジアンですが、np.angle()の出力では3.04→-3.04のように約6.28ラジアンの急激な変化として記録されます。この不連続性により、後段のフィルタリングやFFT解析で偽の高周波成分が発生してしまいます。

位相アンラップ処理では、隣接するサンプル間の位相差が±πを超える場合に2πの整数倍を加減することで、位相の連続性を回復します。

実際のコード

def _perform_vital_analysis(self):
    """バイタルサイン解析を実行"""
    # 静的成分除去後の信号から位相を抽出
    sig_angle_raw = np.angle(sig_dcoff)
    
    # 位相アンラップ処理（連続化）
    sig_angle = np.unwrap(sig_angle_raw)

NumPyのunwrap()関数は、隣接点間の位相差の絶対値がπを超える場合に自動的に2πの調整を行い、連続的な位相変化を復元します。この処理により、呼吸による真の位相変化パターンが明確に観測できるようになります。

呼吸領域のバンドパスフィルタ

原理の説明

位相アンラップ後の信号には、目的とする呼吸成分（0.1-1.0Hz）の他に、センサーノイズや電源ノイズ、体動による低周波成分など、様々な周波数成分が含まれています。呼吸領域バンドパスフィルタは、これらの不要な成分を除去し、呼吸による位相変化のみを抽出するために使用されます。

人間の呼吸周波数は安静時で約0.2-0.5Hz（12-30回/分）の範囲にあり、運動時でも1.0Hzを超えることは稀です。一方、心拍による微小な位相変化は1.0-2.0Hz程度で現れるため、0.1-1.0Hzのバンドパスフィルタにより呼吸成分を選択的に抽出できます。

実際のコード

def _initialize_filters(self):
    """フィルタ係数の初期化"""
    # 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')

def _perform_vital_analysis(self):
    """バイタルサイン解析を実行"""
    # バンドパスフィルタリング（呼吸帯域）
    bpfout = signal.lfilter(self.bpf_b, self.bpf_a, sig_angle)

4次Butterworth バンドパスフィルタを使用することで、通過帯域での平坦な周波数応答と阻止帯域での急峻な減衰特性を実現しています。

STFT処理による周波数推定

原理の説明

STFT（短時間フーリエ変換）法は、時間的に変化する信号の周波数成分を解析する手法です。従来のFFTが信号全体の周波数成分を一括で解析するのに対し、STFTは短い時間窓を設定してその範囲での周波数解析を行うため、時間的に変化する周波数成分の推定が可能になります。

バイタルサイン検出では、呼吸パターンが時間と共に変化する可能性があるため、固定長のFFTよりもSTFTの方が適している場合があります。ただし、添付コードでは処理の簡素化とリアルタイム性を重視し、固定窓長でのFFT解析を採用しています。

実際のコード

def _perform_phase_fft_analysis(self, bpf_signal, time_axis):
    """位相信号のFFT解析を実行（呼吸周波数解析）"""
    try:
        if len(bpf_signal) < 32:  # 最小長チェック
            return None
        
        # FFT実行
        fft_result = np.fft.rfft(bpf_signal)
        fft_magnitude = np.abs(fft_result)
        
        # 周波数軸作成
        dt = time_axis[1] - time_axis[0] if len(time_axis) > 1 else 1.0 / self.frame_freq
        frequencies = np.fft.rfftfreq(len(bpf_signal), d=dt)
        
        # パワースペクトラム（dB）
        power_linear = fft_magnitude ** 2
        power_db = 10 * np.log10(power_linear + 1e-12)
        
        # 呼吸帯域（0.1-1.0Hz）でピーク検出
        valid_range = (frequencies >= 0.1) & (frequencies <= 1.0)
        if np.any(valid_range):
            valid_frequencies = frequencies[valid_range]
            valid_power = power_db[valid_range]
            
            # ピーク周波数を検出
            peak_idx = np.argmax(valid_power)
            peak_frequency = valid_frequencies[peak_idx]
            peak_power = valid_power[peak_idx]
            
            # BPM計算（呼吸数）
            bpm = peak_frequency * 60.0
        else:
            peak_frequency = 0.0
            peak_power = 0.0
            bpm = 0.0
        
        return {
            'frequencies': frequencies,
            'power_db': power_db,
            'peak_frequency': peak_frequency,
            'peak_power': peak_power,
            'bpm': bpm
        }
        
    except Exception as e:
        print(f"位相FFT解析エラー: {e}")
        return None

この実装では、実時間FFT（np.fft.rfft）を使用してパワースペクトラムを計算し、呼吸帯域内での最大パワーを持つ周波数を呼吸周波数として推定しています。最終的に、この周波数に60を掛けることでBPM（Breaths Per Minute）を算出します。

まとめ

本記事では、ミリ波センサーを用いたバイタルサイン検出における信号処理の全工程を解説しました。静的クラッター除去から始まり、位相アンラップ、バンドパスフィルタリング、そして最終的な周波数推定まで、各段階が連携して動作することで、ノイズに埋もれた微小な生体信号から実用的な呼吸情報を抽出できます。
この処理フローにより呼吸周波数の推定が可能になります。

Discussion