基于MATLAB的脑电信号带通滤波器的仿真与比较

张井想++梁雨++许朋

摘 要 脑电信号的频率范围一般在0.5～35Hz。针对脑电信号的频率特性，用MATLAB的信号处理工具箱的函数设计巴特沃斯、切比雪夫和椭圆函数带通滤波器。通过仿真及结果的分析得出更加适合脑电信号软件处理的带通滤波器。

【关键词】脑电信号MATLAB 带通滤波器 仿真

1 引言

脑部疾病长期以来一直威胁着人类的健康，因此对其预防和及时发现在减少脑部疾病危害中极为重要。脑电信号受到工频50Hz信号等噪声及生理信号等干扰，故能够设计一种可以分离出脑电信号的带通滤波器就显得尤为重要。本文从巴特沃斯带通滤波器、切比雪夫带通滤波器和椭圆带通滤波器的MATLAB仿真及结果的分析得出更加适合脑电信号软件编程处理的带通滤波器，从而获取更加纯净的脑电信号。

2 MATLAB简介

MATLAB语言是一种面向工程与科学的计算语言。MATLAB信号处理工具箱提供了设计巴特沃斯、切比雪夫和椭圆函数滤波器等函数，本文利用这些函数，进行了巴特沃斯、切比雪夫和椭圆函数滤波器的程序设计，并通过仿真结果并分析这三种滤波器的优缺点及适用场合。

2.1 巴特沃斯带通滤波器的设计与仿真

2.1.1 butter函数

butter函数是用于设计巴特沃斯的滤波器

[k， l] = butter（n， Wn）；

当Wn = [W1W2]时，它可以设计2n 阶的巴特沃斯带通滤波器，其通带为 W1 < W < W2。

2.1.2 巴特沃斯带通滤波器设计

设计一个巴特沃斯带通滤波器，绘制原始信号、滤波后信号FFT、原始信号FFT和归一化的信噪比图。20Hz和50Hz 正弦波组成原始信号，将噪声 50Hz 的正弦波滤掉，通过函数 butter 设计一组带通滤波器系数，其阶数是2，0.5Hz 到 35Hz的通带频率，采样率 1Kbps。程序代码如下：

fc=1000； %设置1k的采样频率

N=1024； %采样点数

n=0：N-1；

t=0：1/fc：1-1/fc； %时间序列

f=n*fc/N； %频率序列

X1=sin（2*pi*50*t）； %50Hz的噪声

X2=sin（2*pi*20*t）； %20Hz的信号

X=X1+X2； %信号混合

subplot（221）；

plot（t，X）； %绘制原始信号

xlabel（时间）；

ylabel（幅值）；

title（原始信号）；

grid on；

subplot（222）；

Y=fft（X，N）； %绘制原始信号的幅频响应

plot（f，abs（Y））；

xlabel（频率/Hz）；

ylabel（振幅）；

title（原始信号 FFT）；

grid on；

subplot（223）；

Wn=[0.5*2 35*2]/fs； %设置通带 0.5Hz 到 35Hz

[k，l]=butter（1，Wn）； %注意第一个参数虽然是 1，但生成的却是 2 阶 IIR 滤波器系数

Y2=filtfilt（k，l，X）； %计算滤波后的波形 y2

Y3=fft（Y2，N）； %滤波后波形的幅频响应

plot（f，abs（Y3））；

xlabel（频率/Hz）；

ylabel（振幅）；

title（滤波后信号 FFT）；

grid on；

[H，F]=freqz（k，l，512）；

subplot（224）；

plot（F/pi，abs（H））；

xlabel（归一化频率）； %绘制绝对幅频响应

ylabel（幅度）；

P1=sum（X2.^2）； %脑电信号的总功率

P2=sum（（Y2-X2）.^2）； %剩余噪声的功率

SNR=10*log10（P1/P2）； %脑电信号的总功率和剩余噪声的功率比值

title（[Order=，int2str（2）， SNR=，num2str（SNR）]）；

grid on；

Matlab 的仿真结果知，SNR=4.5311。

2.2 切比雪夫带通滤波器设计与仿真

2.2.1 cheby1 函数

cheby1 函数是用来设计切比雪夫I 型的滤波器。

[k，l] = cheby1（n， Rp， Wn）；cheby1 函数可以设计带通切比雪夫I 型数字滤波器，其阻带内为单调，通带内为等纹波。

[k，l] = cheby1（n， Rp， Wn）；Rp 用来确定通带内的纹波，该滤波器的截止频率是Wn，可以设计 n 阶低通切比雪夫I 型数字滤波器。当 Wn=[W1， W2]时，cheby1 函数可以设计出其通带为 W1

2.2.2 切比雪夫带通滤波器设计

设计一个巴特沃斯带通滤波器，绘制原始信号、滤波后信号FFT、原始信号FFT和归一化的信噪比图。20Hz和50Hz 正弦波组成原始信号，将噪声 50Hz 的正弦波滤掉，通过函数 cheby1 设计一组带通滤波器系数，其阶数是2，0.5Hz 到 35Hz的通带频率，采样率 1Kbps，通带纹波 1db。 Matlab程序代码同上，不同部分如下：

Wn=[0.5*2 35*2]/fs； %设置通带 0.5Hz 到 35Hz

[k，l]=cheby1（1，1，Wn）； %注意第一个参数虽然是 1，但生成的却是 2 阶 IIR 滤波器系数

Matlab 的仿真结果知，SNR=8.4301。

2.3 椭圆带通滤波器设计与仿真

2.3.1 ellip函数

ellip函数是用来设计椭圆型的滤波器

[k，l] = ellip（n， Rp， RS， Wn）；

2.3.2 椭圆带通滤波器设计

设计一个巴特沃斯带通滤波器，绘制原始信号、滤波后信号FFT、原始信号FFT和归一化的信噪比图。20Hz和50Hz 正弦波组成原始信号，将噪声 50Hz 的正弦波滤掉，通过函数 ellip 设计一组带通滤波器系数，其阶数是2，0.5Hz 到 35Hz的通带频率，采样率 1Kbps，阻带40db，通带纹波 1db。Matlab程序代码同上，不同部分如下：

Wn=[0.5*2 35*2]/fs； %设置通带 0.5Hz 到 35Hz

[k，l] = ellip（2， 1， 40， Wn）；

Matlab 的仿真结果知，SNR=8.4434。

3 结束语

通过三个带通滤波器的仿真图的信噪比比较可知，巴特沃斯函数设计的数字带通滤波器具有最大的平坦幅度，但其截止频率处的下降斜度会受到损失，使幅度响应衰减较慢。与切比雪夫和巴特沃斯滤波器相比，椭圆函数只需较低的阶数可以设计出衰减更快、下降斜度更大的滤波器，但通带和阻带内均为等纹波。综上，高阶巴特沃斯带通滤波器能满足脑电信号对通带内的幅度响应平坦及截止频率处的下降斜率的要求。

参考文献

[1]李钟慎.基于MATLAB设计巴特沃斯低通滤波器[J].信息技术，2003（03）：49-50+52.

[2]刘凌云，赵鹏宇，弓美桃.基于MATLAB的低通巴特沃斯滤波器仿真[J].数字技术与应用，2013（02）：124.

作者简介

张井想（1991-），男，工程硕士学位。主要研究方向为医疗电子。

作者单位

江苏师范大学 江苏省徐州市 221116