張井想++梁雨++許朋

摘 要 腦電信號的頻率范圍一般在0.5～35Hz。針對腦電信號的頻率特性，用MATLAB的信號處理工具箱的函數(shù)設(shè)計巴特沃斯、切比雪夫和橢圓函數(shù)帶通濾波器。通過仿真及結(jié)果的分析得出更加適合腦電信號軟件處理的帶通濾波器。

【關(guān)鍵詞】腦電信號MATLAB 帶通濾波器 仿真

1 引言

腦部疾病長期以來一直威脅著人類的健康，因此對其預防和及時發(fā)現(xiàn)在減少腦部疾病危害中極為重要。腦電信號受到工頻50Hz信號等噪聲及生理信號等干擾，故能夠設(shè)計一種可以分離出腦電信號的帶通濾波器就顯得尤為重要。本文從巴特沃斯帶通濾波器、切比雪夫帶通濾波器和橢圓帶通濾波器的MATLAB仿真及結(jié)果的分析得出更加適合腦電信號軟件編程處理的帶通濾波器，從而獲取更加純凈的腦電信號。

2 MATLAB簡介

MATLAB語言是一種面向工程與科學的計算語言。MATLAB信號處理工具箱提供了設(shè)計巴特沃斯、切比雪夫和橢圓函數(shù)濾波器等函數(shù)，本文利用這些函數(shù)，進行了巴特沃斯、切比雪夫和橢圓函數(shù)濾波器的程序設(shè)計，并通過仿真結(jié)果并分析這三種濾波器的優(yōu)缺點及適用場合。

2.1 巴特沃斯帶通濾波器的設(shè)計與仿真

2.1.1 butter函數(shù)

butter函數(shù)是用于設(shè)計巴特沃斯的濾波器

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

當Wn = [W1W2]時，它可以設(shè)計2n 階的巴特沃斯帶通濾波器，其通帶為 W1 < W < W2。

2.1.2 巴特沃斯帶通濾波器設(shè)計

設(shè)計一個巴特沃斯帶通濾波器，繪制原始信號、濾波后信號FFT、原始信號FFT和歸一化的信噪比圖。20Hz和50Hz 正弦波組成原始信號，將噪聲 50Hz 的正弦波濾掉，通過函數(shù) butter 設(shè)計一組帶通濾波器系數(shù)，其階數(shù)是2，0.5Hz 到 35Hz的通帶頻率，采樣率 1Kbps。程序代碼如下：

fc=1000； %設(shè)置1k的采樣頻率

N=1024； %采樣點數(shù)

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； %設(shè)置通帶 0.5Hz 到 35Hz

[k，l]=butter（1，Wn）； %注意第一個參數(shù)雖然是 1，但生成的卻是 2 階 IIR 濾波器系數(shù)

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

Y3=fft（Y2，N）； %濾波后波形的幅頻響應

plot（f，abs（Y3））；

xlabel（頻率/Hz）；

ylabel（振幅）；

title（濾波后信號 FFT）；

grid on；

[H，F(xiàn)]=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 的仿真結(jié)果知，SNR=4.5311。

2.2 切比雪夫帶通濾波器設(shè)計與仿真

2.2.1 cheby1 函數(shù)

cheby1 函數(shù)是用來設(shè)計切比雪夫I 型的濾波器。

[k，l] = cheby1（n， Rp， Wn）；cheby1 函數(shù)可以設(shè)計帶通切比雪夫I 型數(shù)字濾波器，其阻帶內(nèi)為單調(diào)，通帶內(nèi)為等紋波。

[k，l] = cheby1（n， Rp， Wn）；Rp 用來確定通帶內(nèi)的紋波，該濾波器的截止頻率是Wn，可以設(shè)計 n 階低通切比雪夫I 型數(shù)字濾波器。當 Wn=[W1， W2]時，cheby1 函數(shù)可以設(shè)計出其通帶為 W1

2.2.2 切比雪夫帶通濾波器設(shè)計

設(shè)計一個巴特沃斯帶通濾波器，繪制原始信號、濾波后信號FFT、原始信號FFT和歸一化的信噪比圖。20Hz和50Hz 正弦波組成原始信號，將噪聲 50Hz 的正弦波濾掉，通過函數(shù) cheby1 設(shè)計一組帶通濾波器系數(shù)，其階數(shù)是2，0.5Hz 到 35Hz的通帶頻率，采樣率 1Kbps，通帶紋波 1db。 Matlab程序代碼同上，不同部分如下：

Wn=[0.5*2 35*2]/fs； %設(shè)置通帶 0.5Hz 到 35Hz

[k，l]=cheby1（1，1，Wn）； %注意第一個參數(shù)雖然是 1，但生成的卻是 2 階 IIR 濾波器系數(shù)

Matlab 的仿真結(jié)果知，SNR=8.4301。

2.3 橢圓帶通濾波器設(shè)計與仿真

2.3.1 ellip函數(shù)

ellip函數(shù)是用來設(shè)計橢圓型的濾波器

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

2.3.2 橢圓帶通濾波器設(shè)計

設(shè)計一個巴特沃斯帶通濾波器，繪制原始信號、濾波后信號FFT、原始信號FFT和歸一化的信噪比圖。20Hz和50Hz 正弦波組成原始信號，將噪聲 50Hz 的正弦波濾掉，通過函數(shù) ellip 設(shè)計一組帶通濾波器系數(shù)，其階數(shù)是2，0.5Hz 到 35Hz的通帶頻率，采樣率 1Kbps，阻帶40db，通帶紋波 1db。Matlab程序代碼同上，不同部分如下：

Wn=[0.5*2 35*2]/fs； %設(shè)置通帶 0.5Hz 到 35Hz

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

Matlab 的仿真結(jié)果知，SNR=8.4434。

3 結(jié)束語

通過三個帶通濾波器的仿真圖的信噪比比較可知，巴特沃斯函數(shù)設(shè)計的數(shù)字帶通濾波器具有最大的平坦幅度，但其截止頻率處的下降斜度會受到損失，使幅度響應衰減較慢。與切比雪夫和巴特沃斯濾波器相比，橢圓函數(shù)只需較低的階數(shù)可以設(shè)計出衰減更快、下降斜度更大的濾波器，但通帶和阻帶內(nèi)均為等紋波。綜上，高階巴特沃斯帶通濾波器能滿足腦電信號對通帶內(nèi)的幅度響應平坦及截止頻率處的下降斜率的要求。

參考文獻

[1]李鐘慎.基于MATLAB設(shè)計巴特沃斯低通濾波器[J].信息技術(shù)，2003（03）：49-50+52.

[2]劉凌云，趙鵬宇，弓美桃.基于MATLAB的低通巴特沃斯濾波器仿真[J].數(shù)字技術(shù)與應用，2013（02）：124.

作者簡介

張井想（1991-），男，工程碩士學位。主要研究方向為醫(yī)療電子。

作者單位

江蘇師范大學 江蘇省徐州市 221116