73 lines
1.9 KiB
Python
73 lines
1.9 KiB
Python
# -*- coding: utf-8 -*-
|
||
"""
|
||
Created on Mon Sep 29 16:14:17 2025
|
||
|
||
@author: 23749
|
||
"""
|
||
|
||
|
||
import numpy as np
|
||
from scipy.signal import butter, filtfilt
|
||
|
||
## 1.Bandpass Filter
|
||
def butter_bandpass(lowcut, highcut, fs, order=4):
|
||
# 滤波器
|
||
nyq = 0.5 * fs #ny:Nyquist频率，即能表示的最大有效频率
|
||
low = lowcut / nyq
|
||
high = highcut / nyq
|
||
b, a = butter(order, [low, high], btype='band') #巴特沃斯滤波器，order=4阶
|
||
return b, a
|
||
|
||
def bandpass_filter(data, lowcut, highcut, fs, order=4):
|
||
b, a = butter_bandpass(lowcut, highcut, fs, order)
|
||
return filtfilt(b, a, data)
|
||
|
||
## 2.Eye Blink Dectection
|
||
def blink_detection(F, fs, Dmin, Dmax, Emin, Emax):
|
||
"""
|
||
波形检测
|
||
输入: 差分特征向量 F, 采样率 fs
|
||
输出: b (0/1), 以及计算出的 d, e
|
||
"""
|
||
|
||
if F is None or len(F) < 3:
|
||
return 0, None, None
|
||
|
||
# 找最大时间(peak) & 最小时间(valley)
|
||
t_peak = np.argmax(F)
|
||
t_valley = np.argmin(F)
|
||
|
||
# 要求 peak 在 valley 之前（符合 blink 形态），否则交换
|
||
if t_valley < t_peak:
|
||
t_peak, t_valley = t_valley, t_peak
|
||
|
||
# 计算持续时间 d (ms)
|
||
d = (t_valley - t_peak) * 1000.0 / fs
|
||
|
||
# 计算能量 e (差分平方和)
|
||
e = np.sum(F[t_peak:t_valley + 1] ** 2)
|
||
|
||
# 阈值判定
|
||
if Dmin <= d <= Dmax and Emin <= e <= Emax:
|
||
b = 1 # 检测到眨眼
|
||
else:
|
||
b = 0 # 否则 no blink
|
||
return b, d, e
|
||
|
||
if __name__ == '__main__':
|
||
import matplotlib.pyplot as plt
|
||
|
||
fs = 250 # 采样率
|
||
t = np.arange(0, 5, 1/fs)
|
||
eog = 0.01 * np.random.randn(len(t)) # 基线+噪声
|
||
|
||
# 模拟眨眼（在 2.0s 注入脉冲）
|
||
center = int(2.0 * fs)
|
||
eog[center:center+5] += 0.5
|
||
eog[center+5:center+15] -= 0.4
|
||
|
||
# 测试 blink_detection
|
||
F = np.diff(eog)
|
||
b, d, e = blink_detection(F, fs, 70, 500, 0.1, 10)
|
||
print(f"Detected: {b}, Duration: {d}ms, Energy: {e}")
|