巴特沃斯滤波器的特点是通频带内的频率响应曲线最大限度平坦,没有起伏,而在阻频带则逐渐下降为零。
一、数据集
UCDDB:下载链接如下St. Vincent's University Hospital / University College Dublin Sleep Apnea Database v1.0.0 (physionet.org)
也可以在终端通过以下命令下载
wget -r -N -c -np https://physionet.org/files/ucddb/1.0.0/
选择EEG通道:C3A2
二、数据处理
1.去除伪影和噪声
使用Butterworth滤波器将信号频率限制在0.25Hz-40Hz的范围内
2.分成五个频段
EEG信号分为五个频段
Delta | 0.25Hz~4Hz |
Theta | 4Hz~8Hz |
Alpha | 8Hz~12Hz |
Sigma | 12Hz~16Hz |
Beta | 16Hz~40Hz |
通过Butterworth滤波器将它们分离
三、代码实现
1.滤波器
def filters(data,sample_rate):
# 去除伪影和噪声
# 由于采样频率是128Hz 带通划分 0.25Hz到40Hz 所以w1=2*0.25/128 w2=2*40/128
w1=2*0.25/sample_rate
w2=2*40/sample_rate
b, a = scipy.signal.butter(2, [w1,w2], 'bandpass') #配置滤波器 2 表示滤波器的阶数
filtedData = scipy.signal.filtfilt(b, a, data) #data为要过滤的信号
# 得到不同频段的波
#delta波 0.25-4
w1=2*0.25/sample_rate
w2=2*4/sample_rate
b, a = scipy.signal.butter(2, [w1,w2], 'bandpass') #配置滤波器 2 表示滤波器的阶数
delta_wave = scipy.signal.filtfilt(b, a, filtedData)
#theta波 4-8
w1=2*4/sample_rate
w2=2*8/sample_rate
b, a = scipy.signal.butter(2, [w1,w2], 'bandpass') #配置滤波器 2 表示滤波器的阶数
theta_wave = scipy.signal.filtfilt(b, a, filtedData)
# alpha 8-12
w1=2*8/sample_rate
w2=2*12/sample_rate
b, a = scipy.signal.butter(2, [w1,w2], 'bandpass') #配置滤波器 2 表示滤波器的阶数
alpha_wave = scipy.signal.filtfilt(b, a, filtedData)
# sigma 12-16
w1=2*12/sample_rate
w2=2*16/sample_rate
b, a = scipy.signal.butter(2, [w1,w2], 'bandpass') #配置滤波器 2 表示滤波器的阶数
sigma_wave = scipy.signal.filtfilt(b, a, filtedData)
#beta 16-40
w1=2*16/sample_rate
w2=2*40/sample_rate
b, a = scipy.signal.butter(2, [w1,w2], 'bandpass') #配置滤波器 2 表示滤波器的阶数
beta_wave = scipy.signal.filtfilt(b, a, filtedData)
return delta_wave,theta_wave,alpha_wave,sigma_wave,beta_wave
2.数据处理
def pre_process(path,file,select_ch = "C3A2"):
psg_fname = path + file + ".rec"
apnea_name = path + file + "_respevt.txt"
psg_f = pyedflib.EdfReader(psg_fname)
# Extract signal from the selected channel
start_datetime = psg_f.getStartdatetime()
print(start_datetime)
ch_names = psg_f.getSignalLabels()
ch_samples = psg_f.getNSamples()
epoch_duration = psg_f.datarecord_duration
# 得到选择的通道的index
select_ch_idx = -1
for s in range(psg_f.signals_in_file):
if ch_names[s] == select_ch:
select_ch_idx = s
break
if select_ch_idx == -1:
raise Exception("Channel not found.")
# 采样频率(一秒钟的数据点数量)
sampling_rate = psg_f.getSampleFrequency(select_ch_idx)
n_epoch_samples = int(epoch_duration * sampling_rate)
delta,theta,alpha,sigma,beta=filters(psg_f.readSignal(select_ch_idx),sampling_rate)
#将几种波reshape成(记录时长,采样频率)
delta=delta.reshape(-1, n_epoch_samples)
theta=theta.reshape(-1, n_epoch_samples)
alpha=alpha.reshape(-1, n_epoch_samples)
sigma=sigma.reshape(-1, n_epoch_samples)
beta=beta.reshape(-1, n_epoch_samples)
n_epochs = psg_f.datarecords_in_file
assert len(delta) == n_epochs, f"signal: {delta.shape} != {n_epochs}"
delta,theta,alpha,sigma,beta=np.asarray(delta),np.asarray(theta),np.asarray(alpha),np.asarray(sigma),np.asarray(beta)
# 画图
plt.subplot(5,1,1)
plt.plot(delta.flatten()[0:1000])
plt.subplot(5,1,2)
plt.plot(theta.flatten()[0:1000])
plt.subplot(5,1,3)
plt.plot(alpha.flatten()[0:1000])
plt.subplot(5,1,4)
plt.plot(sigma.flatten()[0:1000])
plt.subplot(5,1,5)
plt.plot(beta.flatten()[0:1000])
return delta,theta,alpha,sigma,beta
四、完整代码
import pyedflib
import numpy as np
import scipy
import matplotlib.pyplot as plt
import os
def filters(data,sample_rate):
# 去除伪影和噪声
# 由于采样频率是128Hz 带通划分 0.25Hz到40Hz 所以w1=2*0.25/128 w2=2*40/128
w1=2*0.25/sample_rate
w2=2*40/sample_rate
b, a = scipy.signal.butter(2, [w1,w2], 'bandpass') #配置滤波器 2 表示滤波器的阶数
filtedData = scipy.signal.filtfilt(b, a, data) #data为要过滤的信号
# 得到不同频段的波
#delta波 0.25-4
w1=2*0.25/sample_rate
w2=2*4/sample_rate
b, a = scipy.signal.butter(2, [w1,w2], 'bandpass') #配置滤波器 2 表示滤波器的阶数
delta_wave = scipy.signal.filtfilt(b, a, filtedData)
#theta波 4-8
w1=2*4/sample_rate
w2=2*8/sample_rate
b, a = scipy.signal.butter(2, [w1,w2], 'bandpass') #配置滤波器 2 表示滤波器的阶数
theta_wave = scipy.signal.filtfilt(b, a, filtedData)
# alpha 8-12
w1=2*8/sample_rate
w2=2*12/sample_rate
b, a = scipy.signal.butter(2, [w1,w2], 'bandpass') #配置滤波器 2 表示滤波器的阶数
alpha_wave = scipy.signal.filtfilt(b, a, filtedData)
# sigma 12-16
w1=2*12/sample_rate
w2=2*16/sample_rate
b, a = scipy.signal.butter(2, [w1,w2], 'bandpass') #配置滤波器 2 表示滤波器的阶数
sigma_wave = scipy.signal.filtfilt(b, a, filtedData)
#beta 16-40
w1=2*16/sample_rate
w2=2*40/sample_rate
b, a = scipy.signal.butter(2, [w1,w2], 'bandpass') #配置滤波器 2 表示滤波器的阶数
beta_wave = scipy.signal.filtfilt(b, a, filtedData)
return delta_wave,theta_wave,alpha_wave,sigma_wave,beta_wave
def pre_process(path,file,select_ch = "C3A2"):
psg_fname = path + file + ".rec"
apnea_name = path + file + "_respevt.txt"
psg_f = pyedflib.EdfReader(psg_fname)
# Extract signal from the selected channel
start_datetime = psg_f.getStartdatetime()
print(start_datetime)
ch_names = psg_f.getSignalLabels()
ch_samples = psg_f.getNSamples()
epoch_duration = psg_f.datarecord_duration
# 得到选择的通道的index
select_ch_idx = -1
for s in range(psg_f.signals_in_file):
if ch_names[s] == select_ch:
select_ch_idx = s
break
if select_ch_idx == -1:
raise Exception("Channel not found.")
# 采样频率(一秒钟的数据点数量)
sampling_rate = psg_f.getSampleFrequency(select_ch_idx)
n_epoch_samples = int(epoch_duration * sampling_rate)
delta,theta,alpha,sigma,beta=filters(psg_f.readSignal(select_ch_idx),sampling_rate)
#将几种波reshape成(记录时长,采样频率)
delta=delta.reshape(-1, n_epoch_samples)
theta=theta.reshape(-1, n_epoch_samples)
alpha=alpha.reshape(-1, n_epoch_samples)
sigma=sigma.reshape(-1, n_epoch_samples)
beta=beta.reshape(-1, n_epoch_samples)
n_epochs = psg_f.datarecords_in_file
assert len(delta) == n_epochs, f"signal: {delta.shape} != {n_epochs}"
delta,theta,alpha,sigma,beta=np.asarray(delta),np.asarray(theta),np.asarray(alpha),np.asarray(sigma),np.asarray(beta)
# 画图
plt.subplot(5,1,1)
plt.plot(delta.flatten()[0:1000])
plt.subplot(5,1,2)
plt.plot(theta.flatten()[0:1000])
plt.subplot(5,1,3)
plt.plot(alpha.flatten()[0:1000])
plt.subplot(5,1,4)
plt.plot(sigma.flatten()[0:1000])
plt.subplot(5,1,5)
plt.plot(beta.flatten()[0:1000])
return delta,theta,alpha,sigma,beta
if __name__ == '__main__':
path = "UCDDB/" #更换成自己的路径
files_=[]
for root, dirs, files in os.walk(path):
for file in files:
if file[-4:] == '.rec':
files_.append(file[:-4])
files_.sort()
for file in files_:
d,t,a,s,b=pre_process(path,file)