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
| import numpy as np
from typing import Tuple
class SteeringWheelEMG:
"""
方向盘sEMG疲劳检测系统
核心创新：
1. VsESM (Valid sEMG Selection Machine) — 半监督噪声过滤
2. 二层特征：斜率 + 绝对距离（捕捉"波状"过渡）
3. 在行为指标出现前检测到疲劳
"""
def __init__(self, sampling_rate: int = 1000):
self.sampling_rate = sampling_rate
self.window_sec = 5.0
def preprocess(self, raw_emg: np.ndarray) -> np.ndarray:
"""
预处理：带通滤波+降噪
Args:
raw_emg: (N,) 原始EMG信号
Returns:
clean_emg: (N,) 清理后信号
"""
from scipy.signal import butter, filtfilt
nyquist = self.sampling_rate / 2
low = 20 / nyquist
high = 450 / nyquist
b, a = butter(4, [low, high], btype='band')
clean_emg = filtfilt(b, a, raw_emg)
notch_b, notch_a = butter(4, [48/nyquist, 52/nyquist], btype='bandstop')
clean_emg = filtfilt(notch_b, notch_a, clean_emg)
return clean_emg
def extract_features(self, emg: np.ndarray) -> dict:
"""
特征提取（含二层特征）
二层特征是本文创新：
- Slope of Lempel-Zig Complexity: 疲劳过渡的"斜率"
- Absolute distance of SMA: 信号幅度面积的"绝对距离"
"""
window_size = int(self.window_sec * self.sampling_rate)
n_windows = len(emg) // window_size
features = []
for i in range(n_windows):
w = emg[i * window_size:(i + 1) * window_size]
feat = {
"rms": np.sqrt(np.mean(w ** 2)),
"mae": np.mean(np.abs(w)),
"variance": np.var(w),
"zero_crossings": np.sum(np.diff(np.sign(w)) != 0),
"waveform_length": np.sum(np.abs(np.diff(w))),
"sma": np.sum(np.abs(w)) / len(w),
"mean_freq": self._mean_frequency(w),
"median_freq": self._median_frequency(w),
"spectral_entropy": self._spectral_entropy(w),
}
if i > 0:
prev_w = emg[(i-1) * window_size:i * window_size]
lzc_curr = self._lempel_ziv_complexity(w)
lzc_prev = self._lempel_ziv_complexity(prev_w)
feat["lzc_slope"] = (lzc_curr - lzc_prev) / lzc_prev if lzc_prev > 0 else 0
sma_prev = np.sum(np.abs(prev_w)) / len(prev_w)
feat["sma_abs_distance"] = abs(feat["sma"] - sma_prev)
features.append(feat)
return features
def _lempel_ziv_complexity(self, signal: np.ndarray) -> float:
"""LZ复杂度 — 信号复杂度的信息论度量"""
binary = (signal > np.mean(signal)).astype(int)
n = len(binary)
complexity = 1
i = 0
c = 1
while i < n:
j = 0
while i + j < n:
pattern = binary[i:i+j+1]
if not any(binary[k:k+len(pattern)].tolist() == pattern.tolist()
for k in range(i)):
break
j += 1
c += 1
i += j + 1
return c / n
def _mean_frequency(self, signal: np.ndarray) -> float:
"""平均频率"""
freqs = np.fft.rfftfreq(len(signal), 1/self.sampling_rate)
spectrum = np.abs(np.fft.rfft(signal))
return np.sum(freqs * spectrum) / np.sum(spectrum)
def _median_frequency(self, signal: np.ndarray) -> float:
"""中值频率"""
freqs = np.fft.rfftfreq(len(signal), 1/self.sampling_rate)
spectrum = np.abs(np.fft.rfft(signal))
cumsum = np.cumsum(spectrum)
median_idx = np.searchsorted(cumsum, cumsum[-1] / 2)
return freqs[median_idx]
def _spectral_entropy(self, signal: np.ndarray) -> float:
"""频谱熵"""
spectrum = np.abs(np.fft.rfft(signal))
psd = spectrum ** 2
if np.sum(psd) == 0:
return 0
psd = psd / np.sum(psd)
return -np.sum(psd * np.log2(psd + 1e-10))
if __name__ == "__main__":
system = SteeringWheelEMG(sampling_rate=1000)
np.random.seed(42)
alert_emg = np.random.normal(0, 0.1, 30000)
fatigued_emg = np.random.normal(0, 0.05, 30000)
for i in range(0, 30000, 500):
fatigued_emg[i:i+50] += np.random.normal(0, 0.3, 50)
alert_clean = system.preprocess(alert_emg)
fatigued_clean = system.preprocess(fatigued_emg)
alert_feats = system.extract_features(alert_clean)
fatigued_feats = system.extract_features(fatigued_clean)
print("=== 清醒 vs 疲劳 特征对比 ===")
print(f"RMS: 清醒={alert_feats[0]['rms']:.4f} vs 疲劳={fatigued_feats[0]['rms']:.4f}")
print(f"SMA: 清醒={alert_feats[0]['sma']:.4f} vs 疲劳={fatigued_feats[0]['sma']:.4f}")
print(f"MeanFreq: 清醒={alert_feats[0]['mean_freq']:.1f}Hz vs 疲劳={fatigued_feats[0]['mean_freq']:.1f}Hz")
if len(alert_feats) > 1 and len(fatigued_feats) > 1:
print(f"\n二层特征:")
print(f"LZC斜率: 清醒={alert_feats[1].get('lzc_slope', 0):.4f} vs 疲劳={fatigued_feats[1].get('lzc_slope', 0):.4f}")
print(f"SMA距离: 清醒={alert_feats[1].get('sma_abs_distance', 0):.4f} vs 疲劳={fatigued_feats[1].get('sma_abs_distance', 0):.4f}")
|