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
| """
MPVMD: 多变量生理变分模态分解
论文核心：将多台雷达的多通道位移信号联合分解
- 公共频率：呼吸/心率在所有通道共享
- 谐波约束：呼吸谐波 + 心率谐波
- 间隙分量：填充呼吸和心率之间的噪声频段
"""
import numpy as np
from typing import Tuple, List
class MPVMD:
"""
多变量生理变分模态分解
输入：M个雷达通道的位移信号
输出：呼吸频率、心率频率、各通道的模态分量
"""
def __init__(
self,
num_radars: int = 4,
fs: float = 100,
resp_range: Tuple[float, float] = (0.1, 0.5),
hr_range: Tuple[float, float] = (0.8, 3.0),
max_harmonics: int = 3,
max_iter: int = 500,
alpha: float = 2000,
tol: float = 1e-6
):
self.M = num_radars
self.fs = fs
self.resp_range = resp_range
self.hr_range = hr_range
self.K = 2 * max_harmonics + 1
self.max_harmonics = max_harmonics
self.max_iter = max_iter
self.alpha = alpha
self.tol = tol
def decompose(self, signals: np.ndarray) -> dict:
"""
多通道信号联合分解
Args:
signals: (M, N) M个雷达通道，N个采样点
Returns:
result: 包含呼吸频率、心率、各模态
"""
M, N = signals.shape
assert M == self.M, f"期望{self.M}通道，得到{M}"
S = np.fft.fft(signs, axis=1)
omega = np.zeros(self.K)
omega[0] = 0.2
omega[1] = 0.4
omega[2] = 0.6
omega[3] = 1.5
omega[4] = 3.0
if self.K > 5:
omega[5] = 4.5
gap_freq = 0.7
u = np.zeros((self.M, self.K, N), dtype=complex)
for iter in range(self.max_iter):
for k in range(self.K):
freq_k = omega[k]
for m in range(self.M):
other_u = u[m].sum(axis=0) - u[m, k]
residual = S[m] - np.fft.fft(other_u)
freqs = np.fft.fftfreq(N, 1/self.fs)
weight = self.alpha / (alpha + (freqs - freq_k)**2)
u[m, k] = np.fft.ifft(weight * residual)
for k in range(self.K):
power = np.zeros(N)
for m in range(self.M):
power += np.abs(np.fft.fft(u[m, k]))**2
omega[k] = freqs[np.argmax(power)]
if iter > 0 and np.max(np.abs(omega - omega_old)) < self.tol:
break
omega_old = omega.copy()
resp_freq = omega[0]
hr_freq = omega[3] if self.K > 3 else omega[-1]
return {
'resp_freq': resp_freq,
'resp_rate_bpm': resp_freq * 60,
'hr_freq': hr_freq,
'heart_rate_bpm': hr_freq * 60,
'modes': u,
'converged': iter < self.max_iter - 1,
'iterations': iter
}
if __name__ == "__main__":
np.random.seed(42)
fs = 100
duration = 30
t = np.arange(fs * duration) / fs
N = len(t)
resp_freq_true = 0.25
hr_freq_true = 1.2
snr_channels = [1.0, 0.7, 0.5, 0.3]
signals = np.zeros((4, N))
for m in range(4):
snr = snr_channels[m]
resp = snr * np.sin(2 * np.pi * resp_freq_true * t)
hr = snr * 0.3 * np.sin(2 * np.pi * hr_freq_true * t)
noise = np.random.randn(N) * 0.1
signals[m] = resp + hr + noise
mpvmd = MPVMD(num_radars=4, fs=fs)
result = mpvmd.decompose(signals)
print("=" * 60)
print("MPVMD 分布式雷达生命体征检测结果")
print("=" * 60)
print(f"真实呼吸频率: {resp_freq_true * 60:.1f} BPM")
print(f"估计呼吸频率: {result['resp_rate_bpm']:.1f} BPM")
print(f"误差: {abs(result['resp_rate_bpm'] - resp_freq_true * 60):.2f} BPM")
print()
print(f"真实心率: {hr_freq_true * 60:.1f} BPM")
print(f"估计心率: {result['heart_rate_bpm']:.1f} BPM")
print(f"误差: {abs(result['heart_rate_bpm'] - hr_freq_true * 60):.2f} BPM")
print(f"\n收敛: {result['converged']} (迭代{result['iterations']}次)")
print("\n--- 单雷达对比 ---")
for m in range(4):
single_result = MPVMD(num_radars=1, fs=fs).decompose(signals[m:m+1])
print(f" 雷达{m+1} (SNR={snr_channels[m]}): "
f"呼吸={single_result['resp_rate_bpm']:.1f}, "
f"心率={single_result['heart_rate_bpm']:.1f}")
|
