数字滤波器设计及在心电信号滤波中的应用.doc

1、 目 录1 课程设计的目的12 课程设计的原理12.1 Butterworth低通数字滤波器的设计12.2切比雪夫I型数字低通滤波器2 2.3 IIR数字滤波器的性质.23 课程设计设计步骤及结果分析43.1心电数据的导入43.2绘出心电信号的时域图和频谱图43.3加入噪声干扰53.3.1白噪声.53.3.2工频干扰（50Hz）.73.3.3带限chirp噪声.83.4滤波器的设计93.4.1 Butterwort型低通数字滤波器.93.4.2切比雪夫I型数字低通滤波器.123.4.3带阻滤波器的设计.154 心得体会171 课程设计的目的数字滤波器是指输入,输出均为数字信号，通过数值运算处理

2、改变输入信号所含频率成分的相对比例，或者滤除某些频率成分的数字器件或程序。因此，数字滤波的概念和模拟滤波相同，只是信号的形式和实现滤波方法不同。正因为数字滤波通过数值运算实现滤波，所以数字滤波器处理精度高，稳定，体积小，重量轻，灵活，不存在阻抗匹配问题，可以实现模拟滤波器无法实现的特殊滤波功能。希望学生运用数字信号处理课程中所学的理论知识和实验技能，基本掌握数字信号处理的基础理论和处理方法，提高分析和解决信号与信息处理相关问题的能力，为以后的工作和学习打下基础。2 课程设计的原理2.1 Butterworth低通数字滤波器的设计巴特沃斯低通滤波器的平方幅度响应为其中，n为滤波器的阶数，为低通滤

3、波器的截止频率。 该滤波器具有 一些特殊的性质： 对所有的n,都有当时, ； 对所有的n,都有当时, ；是的单调递减函数,即不会出现幅度响应的起伏； 当时,巴特沃斯滤波器趋向于理想的低通滤波器； 在处平方幅度响应的各级导数均存在且等于0，因此在该点上取得最大值，且具有最大平坦特性。 图1展示了2阶、4阶、8阶巴特沃斯低通滤波器的幅频特性。可见阶数n越高，其幅频特性越好，低频检测信号保真度越高，过渡带变窄，即衰减加剧，但半功率点不变。图1 巴特沃斯低通滤波器的幅频特性2.2 切比雪夫I型数字低通滤波器（1）确定数字低通滤波器的技术指标：通带截止频率p、通带衰减ap、阻带截止频率s、阻带衰减as切

4、比雪夫滤波器的振幅平方特性如图2所示：图2 切比雪夫滤波器的振幅平方特性(2)将数字低通滤波器的技术指标转换成模拟低通滤波器的技术指标。如果采用脉冲响不变法，边界频率的转换关系为:如果采用双线性变换法，边界频率的转换关系为(3) 按照模拟低通滤波器的技术指标设计模拟低通滤波器。(4) 利用双线性变换法将模拟滤波器Ha(s)，从s平面转换到z平面，得到数字低通滤波器系统函数H(z)。（5）数字低通技术指标为：p=0.4rad, a p=1dB; s=0.5rad, a s=40Db（6）模拟低通的技术指标为：归一化截止角频率wp=2pi*Fs/Ft; ws=2pi*Fs/Ft(7)利用模拟切比雪

5、夫滤波器设计数字滤波器。通带截止频率为：wp=0.4*pi; 阻带截止频率为：ws=0.5*pi;通带最大衰减为：Rp=1;阻带最大衰减为：As=15;设定周期为1s；模拟低通滤波器的生成：b,a=cheby1(n,1,Wn,low,s);满足设计指标的最小阶数和截止频率：Wnn,Wn=cheb1ord(OmegaP,OmegaS,1,40,s)。最后实现输入输出、幅频特性、相频特性的图形。2.3 IIR数字滤波器的性质无限长冲激响应( IIR) 数字滤波器是数字滤器的一种，数字滤波器根据其冲激响应函数的时域特性,还包括有限长冲激响应( FIR)数字滤波器。IIR数字滤波器的特征是：具有无限持

6、续时间冲激响应；需要用递归模型来实现，这可以从其差分方程得出，也可以从其系数函数为：得出。 数字巴特沃思滤波器属于IIR滤波器，该类滤波器具有特定的性质和设计方法。目前比较成熟的IIR数字滤波器设计方法有两种：1）直接法目前所用的方法主要是：零极点累试法、频域幅度平方误差最小法和时域单位脉冲响应逼近法。直接法的最大优点在于可以设计任意幅频特征的滤波器。2)间接法，目前所用的方法主要是：冲激响应不变法、阶跃响应不变法和双线性法。它们都是借助于 已经成熟的现有低通滤波器原型进行设计，即对数字低通数字滤波器，先将数字低通滤波器的技术指标按希望的设计方法转换为模拟低通滤波器的技术指标，再按指定的模拟低

7、通滤波器的类型设计模拟滤波器H(s)，然后，将模拟滤波器的系统函数H(s)从s平面转换到z平面，得到数字低通滤波器的系统函数H(z)；如所设计的数字滤波器为高通、带通或带阻滤波器，则可借助模拟滤波器的频带变换转换为低通模拟滤波器。 由于直接法设计巴特沃思滤波器相对复杂，在不需要任意幅频特征的情况下，一般采用问接法，同时由于冲激响应不变法和脉冲响应不变法，从s平面转换到z平面的映射为多值映射，容易造成频谱混叠，故而本文采用不会产生频谱混叠的双线性变换法。 3 课程设计步骤及结果分析3.1心电数据的导入将老师给的心电信号原始数据存于桌面，然后我选用桌面上的第一组数据a=load(C:Documen

8、ts and Settingshnie桌面心电信号数据a01.txt);3.2 绘出心电信号的时域图和频谱图将导入的两行数据分别用t,b来替换，然后通过调用plot函数来画出时域图，然后通过对1000个心电数据的幅值进行FFT运算，再次调用plot函数来绘出频域图，具体设计如下：figure(1);subplot(2,1,1);t=a(1:1000,1);b=a(1:1000,2);plot(t,b);title(原始波形图);xlabel(时间（s）);ylabel(幅值（A）);y1=fft(a(:,2),1000);f1=100*(0:999)/1000;subplot(2,1,2);p

9、lot(f1,abs(y1);title(原始频谱图);xlabel(频率（Hz）);ylabel(幅度(dB);图图1.1 原始时域和频谱图通过导入的心电信号数据发现在其频谱图上的020Hz和80100Hz之间的幅值比较大，而在3070Hz之间的幅值相对较小。3.3 加入噪声干扰这里我加入的噪声是：白噪声，50Hz的工频噪声，带限chirp噪声。3.3.1 白噪声通过用s来代表加入白噪声后的信号，并进行数字滤波器的频率响应，对s中的1000个频率点调用plot函数画出加入白噪声后的时域图，再对ws/pi,abshs调用plot函数画出加入白噪声后的频谱图，具体操作如下：q=0.8*rand(

10、1000,1);s=a(:,2)+q;hs,ws=freqz(s,1,1024);abshs=abs(hs);figure(2);subplot(2,1,1);plot(s(1:1000);title(加入白噪声后的时域图);xlabel(时间（s）);ylabel(幅值（A）);subplot(2,1,2);plot(ws/pi,abshs);title(加入白噪声后的频谱图);xlabel(幅度);ylabel(Hz);图1.2 加入白噪声后的时域图和频谱图通过观察加入白噪声后的时域图和频域图，将它与未加入白噪声进行比较，可以发现频谱图在0Hz时的幅度增加的很大，而且又在没有谱线的频率上竟

11、然出现了频谱，这是由于白噪声在所有频率上都有频率造成的。3.3.2 工频干扰（50Hz）用x表示加入工频干扰后的信号，再对工频信号的1000个频率点进行FFT运算，在进行相关的运算后，通过调用plot函数直接绘出加入工频干扰后的时域图和频谱图。具体步骤如下：x2=sin(2*pi*50*t); t=0:0.001:0.001*(1000-1);x1=a(:,2);x=x1+x2;y2=fft(x2,1000);f2=100*(0:999)/1000;figure(3);subplot(2,1,1);plot(t,x);title(加入工频干扰后的时域图);xlabel(时间（s）);ylabe

12、l(幅值（A）);subplot(2,1,2);plot(f2,abs(y2);title(加入工频干扰后的频谱图);xlabel(幅度);ylabel(Hz);图1.3 加入工频干扰后的时域图和频谱图比较图1.1和图1.3，发现加噪声后的幅值有比较微小的变化，而在频谱图在38Hz和63Hz附近增值很大，这是由于50Hz的工频噪声造成的，而在其他频率范围内也有比较明显的变化。3.3.3 带限chirp噪声 同样用k表示加入带限chirp噪声后的信号，再用freqz进行滤波器的频率响应，调用plot对k画出加入带限噪声后的时域和频域图，具体描述如下：p=0.5*chirp(a(:,1),0,a(

13、1000,1),200);k=x1+p;hc,wc=freqz(k,1,1024);abshc=abs(hc);figure(4);subplot(2,1,1);plot(k(1:1000);title(加入带限噪声后的时域图);xlabel(时间（s）);ylabel(幅值（A）);subplot(2,1,2);plot(wc/pi,abshc);title(加入带限噪声后的频谱图);xlabel(幅度);ylabel(Hz);图1.4 加入带限噪声后的时域图和频谱图加入带限chirp噪声干扰后，发现加噪之前时域图是均匀分布的，而加噪后，则变为前面稀疏后面密集的情况了，对于频谱图而言，加噪前

14、是中间凹两边凸，加噪后是前凸后平，在个别点上幅值增加很多，这是由于加入chirp噪声的结果。3.4 滤波器的设计3.4.1 Butterworth型低通数字滤波器用wp和ws表示分别将通带，阻带截止频率的角频率表示，在分别计算阶数n1和截止频率Wn,再设计低通Butterworth型模拟滤波器，然后采用双线性法将模拟滤波器系数变为数字滤波器系数，画出滤波器频谱图，调用filter实现对工频干扰的滤波，用plot函数画出滤除工频干扰后的时域图和滤除白噪声后的频谱图。具体操作如下：figure(5);fs=100;f1=5;f2=10;wp=(f1/fs)*2*pi;ws=(f2/fs)*2*pi

15、;Omegap=2*fs*tan(wp/2);Omegas=2*fs*tan(ws/2);n1,Wn=buttord(Omegap,Omegas,1,50,s);b,a=butter(n1,Wn,s);bz,az=bilinear(b,a,fs);freqz(bz,az,512,fs);y=filter(bz,az,x2);figure(6);subplot(2,1,1)plot(y);title(滤除工频干扰后的时域图);xlabel(时间（s）);ylabel(幅值（A）);y3=fft(y,1000);f1=100*(0:999)/1000;subplot(2,1,2);plot(f1,

16、abs(y3);title(滤除白噪声后的频谱图);xlabel(频率(Hz);ylabel(幅值(dB);图1.5 Butterworth型低通数字滤波器频率响应图图1.6 滤除单频正弦波后的时域图和滤除噪声后的频谱图由于信号处于频段的低频部分，而工频信号的频谱在整个上是呈对称分布的，与源信号的频段分布是相似的，采用低通滤波器将噪声信号的高频部分滤掉，由图1.3和1.6会发现加入工频干扰后的源信号，滤波后时域图的幅值比原来变小了，滤波后的频谱图频率在010Hz和90100Hz之间的幅值变化较大，而在1090Hz的区间，频谱图基本趋于直线。3.4.2切比雪夫I型数字低通滤波器用Wp1,Wp2,

17、Ws1,Ws2表示分别将通带，阻带截止频率的角频率表示,算出频带宽带，计算阶数n1和截止频率WN，再设计切比雪夫I型模拟滤波器，采用双线性法将模拟滤波器系数变为数字滤波器系数，画出切比雪夫I型数字滤波器的频率响应，调用filter实现对白噪声的滤波，再最后调用plot函数画出滤除白噪声后的时域图和频域图。具体过程如下：figure(7);fs=100;f11=10;f12=25;f21=5;f22=30;Wp1=(f11/fs)*2*pi;Ws1=(f21/fs)*2*pi;Wp2=(f12/fs)*2*pi;Ws2=(f22/fs)*2*pi;Omegap1=2*fs*tan(Wp1/2);

18、Omegap2=2*fs*tan(Wp2/2);Omegas1=2*fs*tan(Ws1/2);Omegas2=2*fs*tan(Wp2/2);BW=Omegap2-Omegap1;W0=Omegap1*Omegap2;W00=sqrt(W0);WP=1;WS=WP*(W02-Ws12)/(Ws1*BW);n1,WN=buttord(WP,WS,1,50,s);B,A=cheby1(n1,1,WN,s);BT,AT=lp2bp(B,A,W00,BW);num,den=bilinear(BT,AT,0.5);freqz(num,den,64);y=filter(num,den,s);figure

19、(8);subplot(2,1,1);plot(y);title(滤除白噪声后的时域图);xlabel(时间（s）);ylabel(幅值（A）);s3=fft(y,1000);f1=100*(0:999)/1000;subplot(2,1,2);plot(f1,abs(s3);title(滤除白噪声后的频谱图);xlabel(频率(Hz);ylabel(幅值(dB);图1.7 切比雪夫1型数字滤波器的频率响应图图1.8 滤除白噪声后的时域图和频域图由于白噪声在整个频段上都存在，对与源信号共存的低频信号用选频滤波器是无法滤除的，故而采用低通滤波器将噪声的高频部分去掉。比较图1.2和图1.7，它们

20、在时域图上差不多是一致的，但在频谱图中的图1.2在0Hz上的幅值是很搞的，在滤波后的就明显变小了。3.4.3 带阻滤波器的设计 设计带阻滤波器，再画出滤波器的频率响应图，然后用filter实现对工频干扰信号的滤波，调用plot函数分别画出滤除工频干扰后的时域图和频谱图。具体操作如下：figure(9);d=fir1(1000,0.157 0.17,stop);freqz(d,512);y4=filter(d,1,x2);figure(10);subplot(2,1,1);plot(y4);title(滤除工频干扰后的时域图);xlabel(时间（s）);ylabel(幅值（A）);subplo

21、t(2,1,2);y5=fft(y4,1000);f1=100*(0:999)/1000;plot(f1,abs(y5);title(滤除工频干扰后的频谱图);xlabel(频率(Hz);ylabel(幅值(dB);图1.9带阻滤波器的频率响应图图1.10 滤除工频干扰后的时域图和频谱图带阻滤波器滤除工频干扰在28Hz和74Hz的地方最为显著，带阻滤波器由于工频干扰信号是对称的，而且与源信号的高频部分在相同的频段上，因此加噪后高频部分的频段上的幅值会更高，故我们选择带阻滤波器来滤除噪声干扰的作用。4 心得体会通过本次课程设计，使我基本懂得了数字信号处理的基础理论知识和设计方法，同时在课设时遇到

22、了不少的问题，通过查阅相关书籍或请教老师同学获得了不少的帮助，在此期间也学到了不少的知识，提升了知识层面的容量，扩宽了设计的道路。由于在课设中也遇到了关于MATLAB的相关问题，但也比较轻松的解决了，同时也学习到了在实践时的不同知识，真的让我获益良多。课程设计中需要自我的耐心和解决问题的勇气，正好让我体验了一把，没有足够的耐心只会感觉到它的枯燥乏味，而没有体会到实践带来的动手能力的提高，这是在以后工作中不可缺少的能力之一，在此感谢辛勤的老师和帮助过我的同学，希望我们一起提高和进步。评分表课题名称： 项 目评 价设计方案的合理性与创造性设计与调试结果设计说明书的质量答辩陈述与回答问题情况课程设计周表现情况综合成绩 教师签名： 日 期： 19