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
| import numpy as np
from scipy.signal import butter, filtfilt, find_peaks
class RadarHeartRateEstimator:
"""
雷达心率估计器（基于论文方法概念实现）
解决呼吸谐波对心率信号的干扰
"""
def __init__(self, fs: float = 20.0):
"""
Args:
fs: 采样率 (Hz)，雷达帧率
"""
self.fs = fs
self.resp_range = (0.15, 0.5)
self.hr_range = (1.0, 2.0)
def estimate_heart_rate(self,
chest_displacement: np.ndarray) -> dict:
"""
从胸壁位移信号估计心率
Args:
chest_displacement: 雷达提取的胸壁位移, shape=(N,)
Returns:
{
'heart_rate_bpm': float,
'respiratory_rate_bpm': float,
'confidence': float,
}
"""
freqs, spectrum = self._compute_spectrum(chest_displacement)
resp_freq = self._find_respiratory_freq(freqs, spectrum)
harmonics = self._find_harmonics(freqs, spectrum, resp_freq)
clean_spectrum = self._remove_harmonics(
spectrum, freqs, resp_freq, harmonics
)
hr_freq = self._find_heart_rate(freqs, clean_spectrum)
return {
'heart_rate_bpm': hr_freq * 60,
'respiratory_rate_bpm': resp_freq * 60,
'confidence': self._compute_confidence(
clean_spectrum, freqs, hr_freq
),
}
def _compute_spectrum(self, signal):
"""计算频谱"""
n = len(signal)
freqs = np.fft.rfftfreq(n, 1/self.fs)
spectrum = np.abs(np.fft.rfft(signal))
return freqs, spectrum
def _find_respiratory_freq(self, freqs, spectrum):
"""识别呼吸基频"""
mask = (freqs >= self.resp_range[0]) & (freqs <= self.resp_range[1])
resp_band = spectrum.copy()
resp_band[~mask] = 0
peak_idx = np.argmax(resp_band)
return freqs[peak_idx]
def _find_harmonics(self, freqs, spectrum, fundamental):
"""识别呼吸谐波"""
harmonics = []
for n in range(2, 6):
harmonic_freq = fundamental * n
if harmonic_freq > self.fs / 2:
break
mask = (freqs >= harmonic_freq * 0.9) & \
(freqs <= harmonic_freq * 1.1)
if np.any(mask):
peak_idx = np.argmax(spectrum[mask])
harmonics.append({
'order': n,
'freq': freqs[mask][peak_idx],
'amplitude': spectrum[mask][peak_idx],
})
return harmonics
def _remove_harmonics(self, spectrum, freqs,
fundamental, harmonics):
"""消除呼吸谐波"""
clean_spectrum = spectrum.copy()
for freq_target in [fundamental] + [h['freq'] for h in harmonics]:
mask = (freqs >= freq_target * 0.95) & \
(freqs <= freq_target * 1.05)
clean_spectrum[mask] *= 0.1
return clean_spectrum
def _find_heart_rate(self, freqs, clean_spectrum):
"""从清理后的频谱提取心率"""
mask = (freqs >= self.hr_range[0]) & (freqs <= self.hr_range[1])
hr_band = clean_spectrum.copy()
hr_band[~mask] = 0
peak_idx = np.argmax(hr_band)
return freqs[peak_idx]
def _compute_confidence(self, spectrum, freqs, hr_freq):
"""计算置信度"""
mask = (freqs >= hr_freq * 0.95) & (freqs <= hr_freq * 1.05)
hr_power = np.max(spectrum[mask])
total_power = np.sum(spectrum[(freqs >= self.hr_range[0]) &
(freqs <= self.hr_range[1])])
return min(1.0, hr_power / (total_power + 1e-6))
if __name__ == "__main__":
np.random.seed(42)
estimator = RadarHeartRateEstimator(fs=20.0)
duration = 30
t = np.arange(duration * 20) / 20.0
resp_signal = 4.0 * np.sin(2 * np.pi * 0.3 * t)
resp_h2 = 1.0 * np.sin(2 * np.pi * 0.6 * t)
resp_h3 = 0.5 * np.sin(2 * np.pi * 0.9 * t)
heart_signal = 0.3 * np.sin(2 * np.pi * 1.2 * t)
noise = np.random.randn(len(t)) * 0.1
chest_displacement = resp_signal + resp_h2 + resp_h3 + heart_signal + noise
result = estimator.estimate_heart_rate(chest_displacement)
print("=== 雷达心率估计 ===")
print(f"实际心率: 72.0 次/分")
print(f"估计心率: {result['heart_rate_bpm']:.1f} 次/分")
print(f"实际呼吸: 18.0 次/分")
print(f"估计呼吸: {result['respiratory_rate_bpm']:.1f} 次/分")
print(f"置信度: {result['confidence']:.2f}")
error = abs(result['heart_rate_bpm'] - 72.0)
print(f"\n心率误差: {error:.1f} 次/分")
print(f"准确率: {(1 - error/72)*100:.1f}%")
|
