使用巴特沃斯滤波器从EEG信号中分出五个频段的波

巴特沃斯滤波器的特点是通频带内的频率响应曲线最大限度平坦，没有起伏，而在阻频带则逐渐下降为零。 

一、数据集

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)