滤波器设计及在心电信号滤波中的应用

昆明理工大学信息工程与自动化学院学生实验报告
（ 2017—2018学年 第2学期）
课程名称：数字信号处理实验 开课实验室：信自212 实验日期：
一.设计目的和意义
数字滤波器是指输入,输出均为数字信号，通过数值运算处理改变输入信号所含频率成分的相对比例，或者滤除某些频率成分的数字器件或程序。

因此，数字滤波的概念和模拟滤波相同，只是信号的形式和实现滤波方法不同。

正因为数字滤波通过数值运算实现滤波，所以数字滤波器处理精度高，稳定，体积小，重量轻，灵活，不存在阻抗匹配问题，可以实现模拟滤波器无法实现的特殊滤波功能。

希望学生运用《数字信号处理》课程中所学的理论知识和实验技能，基本掌握数字信号处理的基础理论和处理方法，提高分析和解决信号与信息处理相关问题的能力，为以后的工作和学习打下基础。

二.设计原理
2.1 Butterworth 低通数字滤波器的设计
巴特沃斯低通滤波器的平方幅度响应为
()2
211n
c H jw w w =
⎛⎫
+ ⎪⎝⎭
其中，n 为滤波器的阶数，c w 为低通滤波器的截止频率。

该滤波器具有 一些特殊的性质：
① 对所有的n,都有当0w =时,()2
01
H j =； ②对所有的n,都有当c w w =时,()2
0.5
c H jw =；
③
()
2
H jw 是的单调递减函数,即不会出现幅度响应的起伏；
④当n →+∞时,巴特沃斯滤波器趋向于理想的低通滤波器； ⑤在0w =处平方幅度响应的各级导数均存在且等于0，因此()
2
H jw 在该
点上取得最大值，且具有最大平坦特性。

图1展示了2阶、4阶、8阶巴特沃斯低通滤波器的幅频特性。

可见阶数n 越高，其幅频特性越好，低频检测信号保真度越高，过渡带变窄，即衰减加剧，但半功率点不变。

图1 巴特沃斯低通滤波器的幅频特性
2.2 切比雪夫I 型数字低通滤波器
（1）确定数字低通滤波器的技术指标：通带截止频率ωp 、通带衰减αp 、阻带截止频率ωs 、阻带衰减αs
切比雪夫滤波器的振幅平方特性如图2所示：
图2 切比雪夫滤波器的振幅平方特性
(2)将数字低通滤波器的技术指标转换成模拟低通滤波器的技术指标。

如果采用脉冲响不变法，边界频率的转换关系为:
如果采用双线性变换法，边界频率的转换关系为
(3) 按照模拟低通滤波器的技术指标设计模拟低通滤波器。

(4) 利用双线性变换法将模拟滤波器Ha(s)，从s 平面转换到z 平面，得到数字低通滤波器系统函数H(z)。

（5）数字低通技术指标为：
ω =0.4πrad ，α =1dB ；ω =0.5πrad ，α =40dB （6）模拟低通的技术指标为：
，T=1 ， , 归一化截止角频率wp=2pi*Fs/Ft; ws=2pi*Fs/Ft
(7)利用模拟切比雪夫滤波器设计数字滤波器。

通带截止频率为：wp=0.4*pi; 阻带截止频率为：ws=0.5*pi;通带最大衰减为：Rp=1;阻带最大衰减为：As=15;设定周期为1s ；模拟低通滤波器的生成：[b,a]=cheby1(n,1,Wn,'low','s');满足设计指标的最小阶数和截止频率：
Wn[n,Wn]=cheb1ord(OmegaP,OmegaS,1,40,'s')。

最后实现输入输出、幅频特性、相频特性的图形。

三.详细设计步骤
T
ω
=
Ω)2
1(2ωtg T =
Ω
3.1心电数据的导入
3.2 绘出心电信号的时域图和频谱图
将导入的数据分别用t来替换，通过调用plot函数来画出时域图，然后通过对4000个心电数据的幅值进行FFT运算，再次调用plot函数来绘出频域图，具体设计如下：程序10行
figure(1); %新建图像
subplot(2,1,1); %将2个图画到一个平面
t=a(1:4000,1); %用t替换导入的4000个数据
plot(t);%绘出t的图形
title('原始波形图');xlabel('时间（s）');ylabel('幅值（A）');
y1=fft(a(:,1),4000); %行取全部
f1=100*(0:3999)/4000; %先生成一个0,1,2,...,3999的整数向量，然后对对这个向量的每一项乘以100除以4000.
subplot(2,1,2); %将2个图画到一个平面
plot(f1,abs(y1)); %f1为横坐标，abs(y1)为纵坐标作图
title('原始频谱图');xlabel('频率（Hz）');ylabel('幅度(dB)'); 3.3 加入干扰
加入的干扰是：白噪声高频干扰，50Hz的电源线干扰。

3.3.1 白噪声
通过用s来代表加入白噪声后的信号，并进行数字滤波器的频率响应，对s中4000个频率点调用plot函数画出加入白噪声后的时域图，再对ws/pi,abshs 调用plot函数画出加入白噪声后的频谱图，具体操作如下：程序11行q=50*rand(4000,1); %产生4000行1列的位于（0，1）区间的随机数
s=a(:,1)+q;%产生的随机数与原始信号叠加赋给s
[hs,ws]=freqz(s,1,4096); %求离散系统频响特性的函数freqz()
abshs=abs(hs);%hs取绝对值赋给abshs
figure(2);%新建图像2
subplot(2,1,1);
plot(s(1:4000));%绘出从1取到4000时s的图形
title('加入高频干扰信号后的时域图');xlabel('时间（s）');ylabel('幅值（A）');
subplot(2,1,2);
plot(ws/pi,abshs);%ws/pi为横坐标abshs为纵坐标作图
title('加入高频干扰信号后的频谱图');xlabel('Hz');ylabel('幅度');
通过观察加入白噪声后的时域图和频域图，将它与未加入白噪声进行比较可以发现频谱图在0Hz时的幅度增加的很大，而且又在没有谱线的频率上竟然出现了频谱，这是由于白噪声在所有频率上都有频率造成的。

3.3.2 电源线干扰（50Hz）
用x表示加入电源线干扰后的信号，再对电源线信号的4000个频率点进行FFT运算，在进行相关的运算后，通过调用plot函数直接绘出加入电源线干扰后的时域图和频谱图。

具体步骤如下：程序13行
x2=sin(2*pi*50*t); %x2表示正弦信号
t=0:0.00025:0.00025*(4000-1); %从0开始间隔0.00025取值到
0.00025*(4000-1)
x1=a(:,1);%将矩阵a用x1表示
x=x1+x2; %x1和x2叠加赋给x
y2=fft(x,4000); %对x信号做快速傅里叶变换
f2=100*(0:3999)/4000;%先生成一个0,1,2,...,3999的整数向量，然后对对这个向量的每一项乘以100除以4000.
figure(3); %新建图像3
subplot(2,1,1);
plot(t,x); %t为横坐标，x为纵坐标作图
title('加入电源线干扰后的时域图');xlabel('时间（s）');ylabel('幅值（A）');
subplot(2,1,2);
plot(f2,abs(y2)); %f2为横坐标，abs(y2)为纵坐标作图
title('加入电源线干扰后的频谱图');xlabel('幅');ylabel('Hz'); 3.4 滤波器的设计
3.4.1切比雪夫I型数字低通滤波器
用Wp1,Wp2,Ws1,Ws2表示分别用通带和阻带截止频率的角频率算出频带宽带，计算阶数n1和截止频率WN，再设计切比雪夫I型模拟滤波器，采用双线性法将模拟滤波器系数变为数字滤波器系数，画出切比雪夫I型数字滤波器的频率响应，调用filter实现对白噪声的滤波，再最后调用plot函数画出滤除白噪声后的时域图和频域图。

具体过程如下：程序29行
figure(4); %新建图像4
fs=1000; %采样频率
f11=10;%通带频率
f12=25; %阻带频率
Wp1=(f11/fs)*2*pi; %通带角频率
Wp2=(f12/fs)*2*pi; %阻带角频率
Omegap1=2*fs*tan(Wp1/2); %数字转化为模拟
Omegap2=2*fs*tan(Wp2/2); %数字转化为模拟
BW=Omegap2-Omegap1; %频带宽度
W0=Omegap1*Omegap2;
W00=sqrt(W0);
WP=1;
[n1,WN]=buttord(WP,WS,1,50,'s'); %参数WP和WS分别是通带边界频率和阻带边界频率，通带最大衰减为1阻带最小衰减50，返回的参数n1和WN分别为滤波器的阶数和3dB截止频率
[B,A]=cheby1(n1,1,WN,'s');%计算出阶数为n1、截止频率为WN、通带波纹最大衰减为1的数字低通滤波器，它的返回值a、b分别表示数字低通滤波器的系统函数的分子和分母的多项式系数
[BT,AT]=lp2bp(B,A,W00,BW);
[num,den]=bilinear(BT,AT,0.5); %实现模数的映射
freqz(num,den,64);
y=filter(num,den,s);
figure(5);
subplot(2,1,1);
plot(y);
title('滤除高频信号后的时域图');
xlabel('时间（s）');ylabel('幅值（A）');
s3=fft(y,4000);
f1=100*(0:3999)/4000;
subplot(2,1,2);
plot(f1,abs(s3));
title('滤除高频信号后的频谱图');xlabel('频率(Hz)');ylabel('幅值(dB)');
3.4.2高通滤波器的设计
高通滤波器的设计可通过模拟低通滤波器再经频率变换而实现。

其中模
拟低通滤波器可根据已经存在的典型滤波器，如巴特沃斯滤波器等逼近实现，而由低通到高通转换理论依据在绪论部分已经进行了详细的论证，又
本设计基于MATLAB仿真软件实现，可利用MATLAB信号处理工具箱提供的
各种函数模型实现，可忽略其中的一些复杂的函数变换，从而简化理论设
计和论证，具体如下：程序61行
figure(6)
n=0:0.01:2;
for ii=1:4 %定义循环，产生不同阶数的曲线
switch ii
case 1,N=2;
case 2,N=5;
case 3,N=10;
case 4,N=30;
end
[z,p,k]=buttap(N); %调用Butterworth模拟低通滤波器原型函数
[b,a]=zp2tf(z,p,k); %将零点极点增益形式转换为传递函数形式
[H,w]=freqs(b,a,n); %按n指定的频率点给出频率响应
subplot(2,1,1);
hold on;
plot(w,magH2);
end
xlabel('w/wc');
ylabel('|H(jw)|^2');
title('Butterworth 模拟低通滤波器原型');
text(1.5,0.18,'n=2') %对不同曲线做标记
text(1.3,0.08,'n=5')
text(1.16,0.08,'n=10')
text(0.93,0.98,'n=20')
grid on; %模拟高通滤波器设计
m=0:0.01:2;
for ii=1:4
switch ii
case 1,N=2;
case 2,N=5;
case 3,N=10;
case 4,N=30;
end
[z,p,k]=buttap(N);
[b,a]=zp2tf(z,p,k);
[bt,at]=lp2hp(b,a,0.25*2*pi); %由低通原型滤波器转换为截止频率为
0.25Hz的高通滤波器
[Ht,w]=freqs(bt,at,m);
subplot(2,1,2);
hold on; %hold on 是保存axes内图像用的,如果在新画图像之后不想覆盖原图像就要加上hold on这句话
end
title('模拟高通滤波器');xlabel('w/pi');ylabel('|H(jw)|^2');
text(0.5,0.28,'n=2')
text(1.0,0.12,'n=4')
text(1.5,0.28,'n=10')
text(2.0,0.10,'n=30')
grid on; %模拟高通滤波器性能
pha=angle(Ht); %输出系统的相频特性
figure(7);
subplot(2,1,1);
plot(w,20.*log10(pha));
grid;
title('模拟高通滤波器相频特性'); xlabel('w/wc');ylabel('相位/dB');
mag=abs(Ht); %输出系统的幅频特性
subplot(2,1,2);
plot(w,20*log10(mag));
grid;
title(' 模拟高通滤波器幅频特性');xlabel('w/wc');ylabel('幅度/dB');
3.4.3 带阻滤波器的设计
设计带阻滤波器，再画出滤波器的频率响应图，然后用filter实现对工频干扰信号的滤波，调用plot函数分别画出滤除工频干扰后的时域图和频谱图。

具体操作如下：程序12行
figure(8);
d=fir1(4000,[0.157 0.17],'stop');%
freqz(d,512);
y4=filter(d,1,x2);
figure(9);
subplot(2,1,1);plot(y4);
title('滤除电源线干扰后的时域图');xlabel('时间（s）');ylabel('幅值（A）');
subplot(2,1,2);
y5=fft(y4,4000);
f1=100*(0:3999)/4000;
plot(f1,abs(y5));
title('滤除电源线干扰后的频谱图');xlabel('频率(Hz)');ylabel('幅值(dB)');
四.设计结果及分析
4.1心电数据的导入，心电信号的时域图和频谱图如下：
图3 原始时域和频谱图
通过导入的心电信号数据发现，在其频谱图上的0~20Hz和80~100Hz之间的幅值比较大，而在30~70Hz之间的幅值相对较小。

4.2加入高频干扰信号后的波形图
通过观察加入高频信号后的时域图和频域图，将它与未加入前进行比较，可以发现频谱图在0Hz时的幅度增加的很大，而且又在没有谱线的频率上出现了频谱，这是由于白噪声在所有频率上都有频率造成的。

图4加入高频干扰后的时域和频谱图
图5加入电源线干扰后的时域和频谱图
图6切比雪夫1型数字滤波器的频率响应图
图7滤除高频信号后的时域和频谱图
由于白噪声在整个频段上都存在，对与源信号共存的低频信号用选频滤波器是无法滤除的，故而采用低通滤波器将噪声的高频部分去掉。

图8 Butterworth模拟低通滤波器原型模拟高通滤波器
图9模拟高通滤波器的相频和幅频特性
图10 带阻滤波器的频率响应图
图11滤除电源线干扰后的时域和频域图
带阻滤波器滤除工频干扰在28Hz和74Hz的地方最为显著，带阻滤波器由于工频干扰信号是对称的，而且与源信号的高频部分在相同的频段上，因此加噪
后高频部分的频段上的幅值会更高，故我们选择带阻滤波器来滤除噪声干扰的作用。

五.体会
通过本次课程设计，使我基本懂得了基础理论知识和设计方法，同时也学到了一些数字信号处理的知识，在课设时遇到了不少问题，通过查阅相关书籍和请教老师同学获得了不少的帮助，在此期间也学习了不少知识，提升了知识层面的容量，扩宽了设计的道路。

由于在课设中遇到关于MATLAB的问题，解决之后也学习到了在实践时的不同知识，获益良多。

课程设计中需要具有耐心和解决问题的勇气，没有足够的耐心只会感觉到它的枯燥乏味，而体会不到实践带来的动手能力的提高，这是在以后工作中不可缺少的能力之一，希望以后更有提高和进步。