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
| import numpy as np
from scipy.signal import find_peaks
from typing import Tuple
class ECGFeatureExtractor:
"""
ECG/HRV 特征提取器
论文方法：2秒短窗内的 ECG 统计特征 + HRV
"""
def __init__(self, fs: int = 250):
self.fs = fs
def extract(self, ecg_signal: np.ndarray) -> dict:
"""
提取 ECG 特征
Args:
ecg_signal: 2秒 ECG 信号
Returns:
features: ECG 特征字典
"""
r_peaks = find_peaks(ecg_signal, distance=self.fs*0.3,
height=0.5*np.max(ecg_signal))[0]
rr_intervals = np.diff(r_peaks) / self.fs
features = {}
if len(rr_intervals) >= 2:
features['rr_mean'] = np.mean(rr_intervals)
features['rr_std'] = np.std(rr_intervals)
features['rr_min'] = np.min(rr_intervals)
features['rr_max'] = np.max(rr_intervals)
features['hr_mean'] = 60 / np.mean(rr_intervals)
features['rmssd'] = np.sqrt(np.mean(np.diff(rr_intervals)**2))
features['pnn50'] = np.sum(np.abs(np.diff(rr_intervals)) > 0.05) / len(rr_intervals)
else:
features['ecg_mean'] = np.mean(ecg_signal)
features['ecg_std'] = np.std(ecg_signal)
features['ecg_skew'] = float(np.sign(np.mean((ecg_signal - np.mean(ecg_signal))**3)))
features['ecg_kurtosis'] = np.mean((ecg_signal - np.mean(ecg_signal))**4) / (np.var(ecg_signal)**2 + 1e-8)
features['ecg_energy'] = np.sum(ecg_signal**2)
return features
class EEGFeatureExtractor:
"""EEG 频域特征提取"""
def __init__(self, fs: int = 250):
self.fs = fs
self.bands = {
'delta': (0.5, 4),
'theta': (4, 8),
'alpha': (8, 13),
'beta': (13, 30)
}
def extract(self, eeg_signal: np.ndarray) -> dict:
"""提取多通道 EEG 频域特征"""
from scipy.signal import welch
features = {}
n_channels = eeg_signal.shape[0]
for ch in range(n_channels):
freqs, psd = welch(eeg_signal[ch], fs=self.fs, nperseg=self.fs*2)
for band_name, (f_low, f_high) in self.bands.items():
mask = (freqs >= f_low) & (freqs <= f_high)
power = np.sum(psd[mask])
features[f'ch{ch}_{band_name}_power'] = power
features[f'ch{ch}_{band_name}_rel'] = power / (np.sum(psd) + 1e-10)
theta_power = features[f'ch{ch}_theta_power']
alpha_power = features[f'ch{ch}_alpha_power']
beta_power = features[f'ch{ch}_beta_power']
features[f'ch{ch}_theta_alpha'] = theta_power / (alpha_power + 1e-10)
features[f'ch{ch}_alpha_beta'] = alpha_power / (beta_power + 1e-10)
return features
class FatigueDetectionPipeline:
"""
飞行员疲劳检测完整管线
EEG + ECG → 特征提取 → ANOVA-SVM选择 → XGBoost分类
"""
def __init__(self):
self.eeg_extractor = EEGFeatureExtractor(fs=250)
self.ecg_extractor = ECGFeatureExtractor(fs=250)
self.selected_features = None
def extract_features(
self,
eeg: np.ndarray,
ecg: np.ndarray
) -> dict:
"""提取全部特征"""
eeg_feat = self.eeg_extractor.extract(eeg)
ecg_feat = self.ecg_extractor.extract(ecg)
return {**eeg_feat, **ecg_feat}
def anova_select(self, features_list: list, labels: np.ndarray, top_k: int = 50) -> list:
"""ANOVA 特征筛选"""
from sklearn.feature_selection import f_classif
feature_names = list(features_list[0].keys())
X = np.array([[f[k] for k in feature_names] for f in features_list])
F, p = f_classif(X, labels)
ranked = np.argsort(F)[::-1][:top_k]
selected = [feature_names[i] for i in ranked]
return selected
if __name__ == "__main__":
pipeline = FatigueDetectionPipeline()
eeg = np.random.randn(14, 500)
ecg = np.random.randn(500)
features = pipeline.extract_features(eeg, ecg)
print(f"提取特征数: {len(features)}")
batch = [pipeline.extract_features(np.random.randn(14, 500), np.random.randn(500))
for _ in range(32)]
labels = np.random.randint(0, 2, 32)
selected = pipeline.anova_select(batch, labels, top_k=20)
print(f"ANOVA选择特征: {selected[:5]}...")
print(f"\n=== 论文报告结果 ===")
print(f"被试数: 32")
print(f"时间窗: 2秒")
print(f"跨被试准确率: 88.42%")
print(f"跨片段准确率: 98.36%")
print(f"训练时间: 39.3秒")
print(f"模型: XGBoost + ANOVA-SVM")
|
