1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
| """
mmWave雷达心率估计：克服呼吸谐波干扰
论文: Xu et al., Sensors 2026, 26(17), 5587
核心挑战:
- 心脏位移: ~0.1-0.5mm (60-100 BPM)
- 呼吸位移: ~1-12mm (12-20 BPM)
- 呼吸谐波: 2x/3x/4x呼吸频率 → 24-80 BPM → 与心率重叠!
依赖: numpy, scipy
"""
import numpy as np
from scipy.signal import find_peaks, butter, filtfilt
from scipy.fft import fft, fftfreq
class RadarHeartRateEstimator:
"""
雷达心率估计器
方法: 呼吸谐波抑制 + 心率频谱提取
架构:
1. 雷达相位解调 → 胸壁位移信号
2. 自适应带通滤波 → 分离呼吸和心搏
3. 谐波干扰消除 → 去除呼吸谐波
4. 心率频谱估计 → FFT峰值检测
5. 追踪与平滑 → 滑动窗口心率
"""
def __init__(self,
sample_rate: float = 100,
window_sec: float = 10,
hr_range: tuple = (50, 120),
resp_range: tuple = (10, 30)):
self.fs = sample_rate
self.window = int(window_sec * sample_rate)
self.hr_min, self.hr_max = hr_range
self.resp_min, self resp_max = resp_range
def extract_phase(self, radar_iq: np.ndarray) -> np.ndarray:
"""
从雷达IQ数据提取相位（正比于位移）
Args:
radar_iq: (N,) complex array, I+jQ
Returns:
phase: (N,) 胸壁位移信号
"""
phase = np.arctan2(radar_iq.imag, radar_iq.real)
phase_unwrapped = np.unwrap(phase)
phase_unwrapped -= np.mean(phase_unwrapped)
return phase_unwrapped
def adaptive_bandpass(self, signal: np.ndarray,
low_hz: float, high_hz: float) -> np.ndarray:
"""自适应带通滤波"""
nyq = self.fs / 2
b, a = butter(4, [low_hz/nyq, high_hz/nyq], btype='band')
return filtfilt(b, a, signal)
def suppress_respiratory_harmonics(self, signal: np.ndarray,
resp_freq_hz: float) -> np.ndarray:
"""
抑制呼吸谐波干扰
呼吸频率的2/3/4倍频可能与心率频率重叠
需要自适应陷波滤波
"""
nyq = self.fs / 2
for harmonic in [2, 3, 4, 5]:
notch_freq = resp_freq_hz * harmonic
if notch_freq < nyq - 0.5:
q_factor = 30
b, a = self._notch_filter(notch_freq, q_factor)
signal = filtfilt(b, a, signal)
return signal
def _notch_filter(self, freq: float, q: float = 30):
"""设计IIR陷波滤波器"""
w0 = freq / (self.fs / 2)
b, a = butter(2, [w0 - 0.5/self.fs, w0 + 0.5/self.fs], btype='bandstop')
return b, a
def estimate_respiratory_rate(self, signal: np.ndarray) -> float:
"""估计呼吸频率"""
resp_sig = self.adaptive_bandpass(signal,
self.resp_min/60,
self.resp_max/60)
freqs = fftfreq(len(resp_sig), 1/self.fs)
spectrum = np.abs(fft(resp_sig))
mask = (freqs > self.resp_min/60) & (freqs < self.resp_max/60)
peak_idx = np.argmax(spectrum[mask])
return freqs[mask][peak_idx]
def estimate_heart_rate(self, radar_iq: np.ndarray) -> dict:
"""
完整心率估计流程
Args:
radar_iq: (N,) complex IQ数据
Returns:
dict with heart_rate, respiratory_rate, confidence
"""
phase = self.extract_phase(radar_iq)
resp_freq = self.estimate_respiratory_rate(phase)
resp_rate_bpm = resp_freq * 60
hr_low = self.hr_min / 60
hr_high = self.hr_max / 60
cardiac_sig = self.adaptive_bandpass(phase, hr_low, hr_high)
cardiac_clean = self.suppress_respiratory_harmonics(
cardiac_sig, resp_freq
)
freqs = fftfreq(len(cardiac_clean), 1/self.fs)
spectrum = np.abs(fft(cardiac_clean))
mask = (freqs > hr_low) & (freqs < hr_high)
spectrum_hr = spectrum[mask]
freqs_hr = freqs[mask]
peaks, properties = find_peaks(spectrum_hr, height=np.max(spectrum_hr)*0.3)
if len(peaks) > 0:
best_peak = peaks[np.argmax(spectrum_hr[peaks])]
hr_freq = freqs_hr[best_peak]
hr_bpm = hr_freq * 60
confidence = spectrum_hr[best_peak] / (np.sum(spectrum_hr) + 1e-6)
else:
hr_bpm = 0
confidence = 0
return {
'heart_rate_bpm': round(hr_bpm, 1),
'respiratory_rate_bpm': round(resp_rate_bpm, 1),
'confidence': round(confidence, 3),
'snr_db': round(10 * np.log10(
np.max(spectrum_hr) / (np.mean(spectrum_hr) + 1e-6) + 1e-6
), 1)
}
if __name__ == "__main__":
np.random.seed(42)
fs = 100
T = 30
t = np.arange(0, T, 1/fs)
resp_freq = 0.3
resp_signal = 4.0 * np.sin(2 * np.pi * resp_freq * t)
hr_freq = 1.2
cardiac_signal = 0.3 * np.sin(2 * np.pi * hr_freq * t)
for h in [2, 3, 4]:
resp_signal += (4.0 / h) * np.sin(2 * np.pi * resp_freq * h * t)
phase = resp_signal + cardiac_signal
radar_iq = np.exp(1j * phase) + np.random.randn(len(t)) * 0.01
estimator = RadarHeartRateEstimator(sample_rate=fs, window_sec=30)
result = estimator.estimate_heart_rate(radar_iq)
print(f"真实心率: 72 BPM")
print(f"真实呼吸: 18 BPM")
print(f"---")
print(f"估计心率: {result['heart_rate_bpm']} BPM")
print(f"估计呼吸: {result['respiratory_rate_bpm']} BPM")
print(f"置信度: {result['confidence']}")
print(f"SNR: {result['snr_db']} dB")
|
