mmWave FMCW雷达生命体征：Prony与MUSIC算法首次实现0.81 BPM心跳误差

论文信息

项目	内容
论文标题	Breaking Barriers in Health Monitoring: Multi-Scenario Vital Sign Detection Using Mm-Wave MIMO FMCW Radar
作者	Ehsan Sadegi, Paul Havinga
机构	University of Twente, Pervasive Systems, EEMCS
arXiv	2508.20864
发表	2025年8月

核心创新

首次将Prony和MUSIC算法应用于FMCW雷达生命体征检测，实现心率检测误差0.81 BPM（相比FFT方法7.4%误差大幅提升）。

方法	心率MAE	呼吸率MAE	说明
FFT	5.5 BPM	1.5 BPM	基线
MUSIC	1.8 BPM	1.01 BPM	高分辨率
Prony	0.81 BPM	0.8 BPM	最优

方法详解

1. FMCW雷达信号模型

"""
mmWave FMCW Radar 生命体征检测
论文: arXiv 2508.20864, University of Twente

Prony vs MUSIC vs FFT 对比实现

依赖: numpy, scipy
"""

import numpy as np
from scipy import signal as sig
from scipy.signal import butter, filtfilt, find_peaks

class FMCWRadarVitalSigns:
    """
    FMCW雷达生命体征检测器
    
    信号链:
    1. FMCW chirp发射 → 目标反射 → 混频 → IF信号
    2. range-FFT: 分离不同距离的目标
    3. 相位提取: 从目标距离bin提取相位变化
    4. 信号处理: 带通滤波 → 呼吸/心跳分离
    5. 频率估计: FFT vs MUSIC vs Prony
    """
    
    def __init__(self, fc=77e9, bw=4e9, n_chirps=64, 
                 chirp_duration=0.05, n_samples=256):
        """
        Args:
            fc: 中心频率 77GHz (汽车雷达)
            bw: 带宽 4GHz → 距离分辨率 3.75cm
            n_chirps: 每帧chirp数
            chirp_duration: 单个chirp时长(s)
            n_samples: 每chirp采样数
        """
        self.fc = fc
        self.bw = bw
        self.n_chirps = n_chirps
        self.chirp_duration = chirp_duration
        self.n_samples = n_samples
        self.fs = n_samples / chirp_duration
        
        # 距离分辨率
        self.range_resolution = 3e8 / (2 * bw)  # 3.75cm
        
        # 帧率
        self.frame_rate = 1 / (n_chirps * chirp_duration)  # ~312 Hz
        
    def generate_if_signal(self, target_distance, chest_motion):
        """
        生成FMCW IF信号
        
        IF = f_b * t + 相位变化
        
        f_b = 2 * BW * d / (c * T_c)  (beat frequency)
        """
        t = np.linspace(0, self.chirp_duration, self.n_samples)
        
        # beat频率 (距离信息)
        f_beat = 2 * self.bw * target_distance / (3e8 * self.chirp_duration)
        
        # IF信号 (复数)
        if_signal = np.exp(1j * 2 * np.pi * f_beat * t)
        
        # 添加胸壁运动引起的相位变化
        wavelength = 3e8 / self.fc  # ~3.9mm at 77GHz
        phase_mod = 4 * np.pi * chest_motion / wavelength
        
        if_signal *= np.exp(1j * phase_mod)
        
        return if_signal
    
    def range_fft(self, if_signals):
        """
        Range-FFT: 从IF信号提取距离信息
        
        每个chirp做FFT, 找到目标距离bin
        """
        range_profile = np.fft.fft(if_signals, axis=0)
        range_bins = np.fft.fftfreq(self.n_samples, 1/self.fs) * 3e8 / (2*self.bw)
        
        return range_profile, range_bins
    
    def extract_phase(self, range_profile, target_bin):
        """
        从目标距离bin提取相位变化
        
        跨多个chirp的相位变化 = 胸壁运动
        """
        phase = np.angle(range_profile[target_bin, :])
        
        # 相位解包裹
        phase_unwrapped = np.unwrap(phase)
        
        return phase_unwrapped
    
    def preprocess_phase(self, phase, fs=None):
        """
        相位预处理
        
        1. DC偏移校正
        2. 杂波去除 (低通)
        3. 脉冲噪声抑制 (中值滤波)
        4. 去趋势
        """
        if fs is None:
            fs = self.frame_rate
        
        # 1. DC去除
        phase_dc_removed = phase - np.mean(phase)
        
        # 2. 杂波去除 (去除极低频)
        b_clutter, a_clutter = butter(4, 0.05, btype='high', fs=fs)
        phase_clean = filtfilt(b_clutter, a_clutter, phase_dc_removed)
        
        # 3. 脉冲噪声抑制 (中值滤波)
        phase_clean = sig.medfilt(phase_clean, kernel_size=5)
        
        # 4. 去趋势
        from scipy.signal import detrend
        phase_clean = detrend(phase_clean)
        
        return phase_clean
    
    def extract_breath_heart(self, phase, fs=None):
        """
        分离呼吸和心跳信号
        
        呼吸: 0.05-0.5 Hz (3-30 BPM)
        心跳: 0.8-2.0 Hz (48-120 BPM)
        """
        if fs is None:
            fs = self.frame_rate
        
        # 呼吸带通
        b_br, a_br = butter(4, [0.05, 0.5], btype='band', fs=fs)
        breath = filtfilt(b_br, a_br, phase)
        
        # 心跳带通
        b_hr, a_hr = butter(4, [0.8, 2.0], btype='band', fs=fs)
        heart = filtfilt(b_hr, a_hr, phase)
        
        return breath, heart
    
    def estimate_rate_fft(self, signal_data, fs):
        """FFT频率估计"""
        n = len(signal_data)
        freq = np.fft.rfftfreq(n, 1/fs)
        spectrum = np.abs(np.fft.rfft(signal_data))
        
        peak_idx = np.argmax(spectrum)
        return freq[peak_idx] * 60  # BPM
    
    def estimate_rate_music(self, signal_data, fs, n_harmonics=4):
        """
        MUSIC算法频率估计
        
        高分辨率子空间方法
        将信号自相关矩阵分解为信号子空间+噪声子空间
        在噪声子空间中搜索频率峰值
        """
        n = len(signal_data)
        
        # 构建自相关矩阵
        m = 2 * n_harmonics  # 子空间维度
        
        # Hankel矩阵 → 自相关
        from scipy.linalg import hankel, svd
        
        # 构建Hankel矩阵
        hankel_matrix = hankel(signal_data[:n-m+1], signal_data[n-m:])
        
        # SVD分解
        U, s, Vh = svd(hankel_matrix)
        
        # 噪声子空间 (最后 m-n_harmonics 个奇异向量)
        noise_subspace = Vh[n_harmonics:, :].T @ Vh[n_harmonics:, :]
        
        # MUSIC伪谱
        freq_search = np.linspace(0.05, 2.0, 1000)
        pseudo_spectrum = np.zeros(len(freq_search))
        
        for i, f in enumerate(freq_search):
            steering = np.exp(-1j * 2 * np.pi * f * np.arange(m) / fs)
            pseudo_spectrum[i] = 1.0 / (steering @ noise_subspace @ steering.conj())
        
        # 找峰值
        peak_idx = np.argmax(pseudo_spectrum)
        return freq_search[peak_idx] * 60  # BPM
    
    def estimate_rate_prony(self, signal_data, fs, n_harmonics=4):
        """
        Prony算法频率估计
        
        将信号建模为指数衰减正弦和:
        x[n] = Σ A_i * exp(j*2πf_i*n/fs) * exp(σ_i*n)
        
        通过线性预测方程求解频率
        """
        n = len(signal_data)
        m = 2 * n_harmonics  # 模型阶数
        
        # 构建线性预测方程
        # R * a = -r (Yule-Walker类似)
        from scipy.linalg import toeplitz
        
        # 自相关
        r = np.correlate(signal_data, signal_data, mode='full')[n-1:]
        r = r[:m+1]
        
        # 构建Toeplitz矩阵
        R = toeplitz(r[:m], r[:m])
        
        # 求解AR系数
        try:
            a = np.linalg.solve(R, -r[1:m+1])
        except np.linalg.LinAlgError:
            return self.estimate_rate_fft(signal_data, fs)
        
        # 求解多项式根
        # 1 + a1*z^-1 + a2*z^-2 + ... = 0
        poly_coeffs = np.concatenate([[1], a])
        roots = np.roots(poly_coeffs)
        
        # 提取频率 (单位圆上的根)
        freqs = np.angle(roots) * fs / (2 * np.pi)
        
        # 选择正频率中幅度最大的
        positive_freqs = freqs[(freqs > 0.05) & (freqs < 2.0)]
        
        if len(positive_freqs) == 0:
            return self.estimate_rate_fft(signal_data, fs)
        
        # 按幅值排序 (计算每个频率成分的幅值)
        amplitudes = np.abs(roots[(freqs > 0.05) & (freqs < 2.0)])
        best_freq = positive_freqs[np.argmax(amplitudes)]
        
        return best_freq * 60  # BPM


# ========== 测试 ==========
if __name__ == "__main__":
    np.random.seed(42)
    
    radar = FMCWRadarVitalSigns(
        fc=77e9, bw=4e9, n_chirps=128,
        chirp_duration=0.02, n_samples=256
    )
    
    print("=" * 60)
    print("FMCW雷达生命体征检测: FFT vs MUSIC vs Prony")
    print("=" * 60)
    
    print(f"\n雷达配置:")
    print(f"  中心频率: {radar.fc/1e9:.0f} GHz")
    print(f"  带宽: {radar.bw/1e9:.0f} GHz")
    print(f"  距离分辨率: {radar.range_resolution*100:.2f} cm")
    print(f"  帧率: {radar.frame_rate:.1f} Hz")
    
    # 模拟目标
    target_distance = 0.8  # 80cm (驾驶员距离)
    
    # 胸壁运动
    breath_rate_true = 12   # BPM
    heart_rate_true = 75    # BPM
    breath_amp = 0.004      # 4mm
    heart_amp = 0.0003      # 0.3mm
    
    # 生成10秒数据
    duration = 10
    n_chirps_total = int(duration * radar.frame_rate)
    t = np.arange(n_chirps_total) / radar.frame_rate
    
    chest_motion = (breath_amp * np.sin(2*np.pi*breath_rate_true/60*t) +
                    heart_amp * np.sin(2*np.pi*heart_rate_true/60*t))
    
    # 生成IF信号 (每个chirp)
    if_signals = []
    for i in range(n_chirps_total):
        if_sig = radar.generate_if_signal(target_distance + chest_motion[i], 
                                           chest_motion[i])
        if_signals.append(if_sig)
    
    if_signals = np.array(if_signals).T  # (n_samples, n_chirps)
    
    # Range-FFT
    range_profile, range_bins = radar.range_fft(if_signals)
    target_bin = np.argmin(np.abs(range_bins - target_distance))
    
    # 提取相位
    phase = radar.extract_phase(range_profile, target_bin)
    
    # 预处理
    phase_clean = radar.preprocess_phase(phase, radar.frame_rate)
    
    # 分离呼吸/心跳
    breath, heart = radar.extract_breath_heart(phase_clean, radar.frame_rate)
    
    # 频率估计对比
    print(f"\n真实值: 呼吸={breath_rate_true} BPM, 心率={heart_rate_true} BPM")
    print(f"\n{'方法':<15} {'呼吸(BPM)':<15} {'心跳(BPM)':<15} {'呼吸误差':<12} {'心跳误差'}")
    print("-" * 70)
    
    for method_name, method in [('FFT', radar.estimate_rate_fft),
                                 ('MUSIC', radar.estimate_rate_music),
                                 ('Prony', radar.estimate_rate_prony)]:
        breath_est = method(breath, radar.frame_rate)
        heart_est = method(heart, radar.frame_rate)
        
        breath_err = abs(breath_est - breath_rate_true)
        heart_err = abs(heart_est - heart_rate_true)
        
        print(f"{method_name:<15} {breath_est:<15.2f} {heart_est:<15.2f} "
              f"{breath_err:<12.2f} {heart_err:.2f}")
    
    # 论文报告值对比
    print(f"\n{'='*60}")
    print(f"论文报告性能 (多场景平均):")
    print(f"{'方法':<15} {'心率MAE':<15} {'呼吸MAE':<15}")
    print(f"{'FFT':<15} {'5.5 BPM':<15} {'1.5 BPM':<15}")
    print(f"{'MUSIC':<15} {'1.8 BPM':<15} {'1.01 BPM':<15}")
    print(f"{'Prony':<15} {'0.81 BPM':<15} {'0.8 BPM':<15}")

2. 多场景验证

场景	距离	呼吸率	心率	特殊条件
正常坐姿	50cm	12-20	60-80	基线
远距离	80cm	12-20	60-80	SNR降低
侧向	45°	12-20	60-80	截面积减小
运动后	50cm	25-30	100-140	高频高幅
哮喘模拟	50cm	30-35	80-100	高呼吸率
冥想	50cm	6-8	50-60	低频低幅

3. Prony vs MUSIC vs FFT

特性	FFT	MUSIC	Prony
分辨率	受窗口限制	超分辨率	超分辨率
计算复杂度	O(N log N)	O(N²)	O(N²)
心率MAE	5.5 BPM	1.8 BPM	0.81 BPM
呼吸MAE	1.5 BPM	1.01 BPM	0.8 BPM
谐波抑制	❌ 弱	✅ 强	✅ 最强
实时性	✅ 快	⚠️ 中等	⚠️ 中等
参数敏感性	低	中(子空间维度)	中(模型阶数)
噪声鲁棒	一般	好	最好

IMS开发启示

1. 车载77GHz生命体征方案

2. 硬件部署

组件	型号	参数	说明
雷达	TI AWR2944	77GHz, 4Tx4Rx, 4GHz BW	车规级
备选	Calterah Alps	60GHz, 2T4R	座舱专用
处理器	QCS8255	26 TOPS	Prony/MUSIC需NPU加速
采样率	10Hz (10s窗口)	-	生命体征不需要高频
功耗	<2W	-	持续监测

3. 算法选择建议

使用场景	推荐算法	原因
快速初筛	FFT	计算快，粗略估计
精确检测	Prony	0.81 BPM误差，噪声鲁棒
多目标	MUSIC	可分离多个目标频率
资源受限	FFT → Prony	先FFT粗估，再Prony精修
实时监测	Prony (10s窗口)	计算量可控

4. 与IMS现有模块协同

IMS模块	生命体征补充	融合价值
疲劳检测	心率下降+呼吸变浅	超越PERCLOS，生理直接信号
酒驾检测	心率升高+呼吸不稳	补充BiFuseNet面部检测
压力监测	HRV心率变异性	新功能模块
CPD儿童	呼吸确认活体	补充UWB多站方案
健康预警	异常呼吸/心率	新增值功能

5. Euro NCAP 2026对接

Euro NCAP要求	FMCW生命体征方案	可行性
CPD检测活体	✅ 呼吸检测	已验证
驾驶员状态	✅ 心率+呼吸	已验证
无接触	✅ 雷达	满足
隐私	✅ 无摄像头	满足
全天候	✅ 不受光照影响	满足

参考文献

1. Sadegi, E., Havinga, P. “Breaking Barriers in Health Monitoring: Multi-Scenario Vital Sign Detection Using Mm-Wave MIMO FMCW Radar.” arXiv: 2508.20864 (2025).
2. Prony, R. “Essai expérimental et analytique.” (1795) - 经典信号分解方法
3. Schmidt, R.O. “Multiple Emitter Location and Signal Parameter Estimation.” - MUSIC算法
4. TI AWR2944: https://www.ti.com/product/AWR2944

---

论文的核心贡献是将信号处理领域的高分辨率方法(Prony/MUSIC)首次应用于FMCW雷达生命体征检测，将心率误差从FFT的5.5 BPM降至0.81 BPM。对于IMS，这意味着77GHz车载雷达可以在不增加摄像头的情况下实现医疗级生命体征监测，是”健康座舱”的核心传感器。