【原】脑电信号滤波-代码实现

脑机接口社区 2021-07-01

展开全文

更多技术，第一时间送达

可以通过对脑电信号进行分析得到较为准确反映用户的行为以及思想方式，在进行分析前需要提取脑电信号中的Delta波，Theta波，Alpha波，Beta波以及Gamma波。

Delta波：0.5-4Hz;

Theta波：4-8Hz;

Alpha波：8-13Hz;

Beta波：13-32Hz;

Gamma波：>32Hz；

由于在脑电信号的采集过程中，会受到周围噪声的干扰以及工频噪声的影响，因此，需要首先对原始信号进行噪声的滤除，通过利用带通滤波器，将0.5Hz - 50Hz频率之间的信号取出。使用巴特沃斯带通滤波器进行信号的滤波。使用巴特沃斯滤波器，需要设置如下参数：

采样频率fc;

采样点数N；

设置频率序列f；

设置通带频率Wn；

滤波案例

1.Delta波

[header,data] = edfread('Affaf Ikram 20121020 1839.L1.edf'); fp_data=data(1,:); fp_data=resample(fp_data,100,500); %降采样率为100hz fc = 600; %采样频率600Hz N = 30000; %采样点数 n = 0:N-1; f = n*fc/N; %频率序列 Wn = [0.5*2 4*2]/fc;%设置通带为0.5-4Hz Delta波 [k,l] = butter(2,Wn);%4阶IIR滤波器 result = filtfilt(k,l,fp_data); tmp = fft(result); fp1_Delta= ifft(tmp); figure, subplot(211),plot(f(1:N/2),abs(tmp(1:N/2)*2/N)); title('滤波后频谱结果，频率:0.5-4Hz'); axis([0,100,-inf,inf]) subplot(212),plot(fp1_Delta); title('脑波滤波后Delta波结果图像');

2.Theta波

[header,data] = edfread('Affaf Ikram 20121020 1839.L1.edf'); fp_data=data(1,:); fp_data=resample(fp_data,100,500); %降采样率为100hz fc = 600; %采样频率600Hz N = 30000; %采样点数 n = 0:N-1; f = n*fc/N; %频率序列 Wn = [4*2 8*2]/fc;%设置通带为4-8Hz Theta波 [k,l] = butter(2,Wn);%4阶IIR滤波器 result = filtfilt(k,l,fp_data); tmp = fft(result); fp1_Delta= ifft(tmp); figure, subplot(211),plot(f(1:N/2),abs(tmp(1:N/2)*2/N)); title('滤波后频谱结果，频率:4-8Hz'); axis([0,100,-inf,inf]) subplot(212),plot(fp1_Delta); title('脑波滤波后Theta波结果图像');

3.Alpha波

[header,data] = edfread('Affaf Ikram 20121020 1839.L1.edf'); fp_data=data(1,:); fp_data=resample(fp_data,100,500); %降采样率为100hz fc = 600; %采样频率600Hz N = 30000; %采样点数 n = 0:N-1; f = n*fc/N; %频率序列 Wn = [8*2 13*2]/fc;%设置通带为8-13Hz Alpha波 [k,l] = butter(2,Wn);%4阶IIR滤波器 result = filtfilt(k,l,fp_data); tmp = fft(result); fp1_Delta= ifft(tmp); figure, subplot(211),plot(f(1:N/2),abs(tmp(1:N/2)*2/N)); title('滤波后频谱结果，频率:8-13Hz'); axis([0,100,-inf,inf]) subplot(212),plot(fp1_Delta); title('脑波滤波后Alpha波结果图像');

4.Beta波

[header,data] = edfread('Affaf Ikram 20121020 1839.L1.edf'); fp_data=data(1,:); fp_data=resample(fp_data,100,500); %降采样率为100hz fc = 600; %采样频率600Hz N = 30000; %采样点数 n = 0:N-1; f = n*fc/N; %频率序列 Wn = [13*2 32*2]/fc;%设置通带为13-32Hz Beta波 [k,l] = butter(2,Wn);%4阶IIR滤波器 result = filtfilt(k,l,fp_data); tmp = fft(result); fp1_Delta= ifft(tmp); figure, subplot(211),plot(f(1:N/2),abs(tmp(1:N/2)*2/N)); title('滤波后频谱结果，频率:13-32Hz'); axis([0,100,-inf,inf]) subplot(212),plot(fp1_Delta); title('脑波滤波后 Beta波结果图像');

5.Gamma波

[header,data] = edfread('Affaf Ikram 20121020 1839.L1.edf'); fp_data=data(1,:); fp_data=resample(fp_data,100,500); %降采样率为100hz fc = 600; %采样频率600Hz N = 30000; %采样点数 n = 0:N-1; f = n*fc/N; %频率序列 Wn = [32*2 40*2]/fc;%设置通带为32-50Hz Gamma波 [k,l] = butter(2,Wn);%4阶IIR滤波器 result = filtfilt(k,l,fp_data); tmp = fft(result); fp1_Delta= ifft(tmp); figure, subplot(211),plot(f(1:N/2),abs(tmp(1:N/2)*2/N)); title('滤波后频谱结果，频率:32-50Hz'); axis([0,100,-inf,inf]) subplot(212),plot(fp1_Delta); title('脑波滤波后Gamma波结果图像');