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

1、西南科技大学课 程 设 计 报 告课程名称： 通信系统课程设计 设计名称：滤波器设计及在心电信号滤波中的应用姓 名： 学 号: 班 级： 指导教师： 起止日期： 2013.6.24-2013.7.6 西南科技大学信息工程学院制课 程 设 计 任 务 书学生班级： 学生姓名： 学号： 设计名称： 滤波器设计及在心电信号滤波中的应用 起止日期： 2013.6.24-2013.7.6 指导教师： 设计要求：一般正常人的心电信号频率在0.7100HZ范围内，幅度为10v（胎儿）5mv（成人）。人体心电信号微弱，信噪比小，因此，在采集心电信号时，易受到仪器、人体活动等因素的影响，而且所采集的心电信号常伴

2、有干扰。 采集心电数据时，由于人的说话呼吸，常常会混有约为0.1Hz到0.25Hz频段的干扰，对于这些低频干扰，可以让信号通过一个高频滤波器，低截止频率设置为0.25，来滤波低频信号； 对于高频信号干扰，可以让信号再通过一个低频滤波器，其中截止频率设置为99Hz；由于电子设备采集到的信号经常会混有电源线干扰，电源线干扰是以50 Hz为中心的窄带噪声，带宽小于1Hz。对于电源线干扰，可以让信号通过一个带阻滤波器。 某一个实际心电图信号的采样序列存于文件ecg.txt中，共4000个采样点。要求根据以上指标进行如下设计： 1、设计一个低通滤波器滤除心电信号中的高频信号干扰，对滤波前后的心电信号的频

3、谱进行分析比较。 2、设计一个高通滤波器滤除心电信号中的低频干扰，对滤波前后的心电信号的频谱进行分析比较。 3、设计一个带阻滤波器（50Hz陷波器）滤除心电信号中的电源线干扰，对滤波前后的心电信号的频谱进行分析比较。说明：1、使用MATLAB或其它编程软件编程实现。2、如果滤波效果不明显，则需变动滤波器参数指标，重新设计滤波器，通过频谱分析，多次试验确定最合适的滤波器。课 程 设 计 学 生 日 志时间设计内容2013.6.24-2013.6.30查阅资料，确定初步方案2013.6.31-2013.7.2编写代码，调试改错2013.7.3-2013.7.5完善代码并撰写课程设计报告 2013.

4、7.6答辩课 程 设 计 评 语 表指导教师评语： 成绩： 指导教师： 年 月 日滤波器设计及在心电信号滤波中的应用一 设计目的和意义数字滤波器是指输入,输出均为数字信号，通过数值运算处理改变输入信号所含频率成分的相对比例，或者滤除某些频率成分的数字器件或程序。因此，数字滤波的概念和模拟滤波相同，只是信号的形式和实现滤波方法不同。正因为数字滤波通过数值运算实现滤波，所以数字滤波器处理精度高，稳定，体积小，重量轻，灵活，不存在阻抗匹配问题，可以实现模拟滤波器无法实现的特殊滤波功能。希望学生运用数字信号处理课程中所学的理论知识和实验技能，基本掌握数字信号处理的基础理论和处理方法，提高分析和解决信号

5、与信息处理相关问题的能力，为以后的工作和学习打下基础。二 设计原理 2.1 Butterworth低通数字滤波器的设计巴特沃斯低通滤波器的平方幅度响应为其中，n为滤波器的阶数，为低通滤波器的截止频率。 该滤波器具有 一些特殊的性质： 对所有的n,都有当时, ； 对所有的n,都有当时, ；是的单调递减函数,即不会出现幅度响应的起伏； 当时,巴特沃斯滤波器趋向于理想的低通滤波器； 在处平方幅度响应的各级导数均存在且等于0，因此在该点上取得最大值，且具有最大平坦特性。 图1展示了2阶、4阶、8阶巴特沃斯低通滤波器的幅频特性。可见阶数n越高，其幅频特性越好，低频检测信号保真度越高，过渡带变窄，即衰减加

6、剧，但半功率点不变。图1 巴特沃斯低通滤波器的幅频特性2.2 切比雪夫I型数字低通滤波器（1）确定数字低通滤波器的技术指标：通带截止频率p、通带衰减ap、阻带截止频率s、阻带衰减as切比雪夫滤波器的振幅平方特性如图2所示：图2 切比雪夫滤波器的振幅平方特性(2)将数字低通滤波器的技术指标转换成模拟低通滤波器的技术指标。如果采用脉冲响不变法，边界频率的转换关系为:如果采用双线性变换法，边界频率的转换关系为(3) 按照模拟低通滤波器的技术指标设计模拟低通滤波器。(4) 利用双线性变换法将模拟滤波器Ha(s)，从s平面转换到z平面，得到数字低通滤波器系统函数H(z)。（5）数字低通技术指标为： p=

7、0.4rad，ap=1dB； s=0.5rad，as=40dB（6）模拟低通的技术指标为： p=2Ttan12p，T=1 p=2tan0.1=0.65rad/s，p=1dB s=2tan0.15=1.019/s, s=15d归一化截止角频率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');满足设计

8、指标的最小阶数和截止频率：Wnn,Wn=cheb1ord(OmegaP,OmegaS,1,40,'s')。最后实现输入输出、幅频特性、相频特性的图形。三 详细设计步骤3.1心电数据的导入将老师给的心电信号原始数据存于桌面文件夹，导入数据代码为：a=load('C:UsersAdministratorDesktop题目七数据文件题目七数据文件ecg.txt');3.2 绘出心电信号的时域图和频谱图将导入的数据分别用t来替换，通过调用plot函数来画出时域图，然后通过对4000个心电数据的幅值进行FFT运算，再次调用plot函数来绘出频域图，具体设计如下：figur

9、e(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为横

10、坐标，abs(y1)为纵坐标作图title('原始频谱图');xlabel('频率（Hz）');ylabel('幅度(dB)');3.3 加入干扰这里我加入的干扰是：白噪声高频干扰，50Hz的电源线干扰。3.3.1 白噪声通过用s来代表加入白噪声后的信号，并进行数字滤波器的频率响应，对s中4000个频率点调用plot函数画出加入白噪声后的时域图，再对ws/pi,abshs调用plot函数画出加入白噪声后的频谱图，具体操作如下：q=50*rand(4000,1); %产生4000行1列的位于（0，1）区间的随机数s=a(:,1)+q; %产生的随机

11、数与原始信号叠加赋给shs,ws=freqz(s,1,4096); %求离散系统频响特性的函数freqz()abshs=abs(hs); %hs取绝对值赋给abshsfigure(2); %新建图像2subplot(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('加入高频干扰

12、信号后的频谱图');xlabel('Hz');ylabel('幅度');通过观察加入白噪声后的时域图和频域图，将它与未加入白噪声进行比较可以发现频谱图在0Hz时的幅度增加的很大，而且又在没有谱线的频率上竟然出现了频谱，这是由于白噪声在所有频率上都有频率造成的。3.3.2 电源线干扰（50Hz）用x表示加入电源线干扰后的信号，再对电源线信号的4000个频率点进行FFT运算，在进行相关的运算后，通过调用plot函数直接绘出加入电源线干扰后的时域图和频谱图。具体步骤如下：x2=sin(2*pi*50*t); %x2表示正弦信号t=0:0.00025:0.000

13、25*(4000-1); %从0开始间隔0.00025取值到 0.00025*(4000-1)x1=a(:,1); %将矩阵a用x1表示x=x1+x2; %x1和x2叠加赋给xy2=fft(x,4000); %对x信号做快速傅里叶变换f2=100*(0:3999)/4000; %先生成一个0,1,2,.,3999的整数向量，然后对对这个向量的每一项乘以100除以4000.figure(3); %新建图像3subplot(2,1,1);plot(t,x); %t为横坐标，x为纵坐标作图title('加入电源线干扰后的时域图');xlabel('时间（s）');yl

14、abel('幅值（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实现对白噪

15、声的滤波，再最后调用plot函数画出滤除白噪声后的时域图和频域图。具体过程如下：figure(4); %新建图像4fs=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

16、,'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)

17、;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 高通滤波器的设计 高通滤波器的设计可通过模拟低通滤波器再经频率变换而实现。其中模拟低通滤波器

18、可根据已经存在的典型滤波器，如巴特沃斯滤波器等逼近实现，而由低通到高通转换理论依据在绪论部分已经进行了详细的论证，又本设计基于MATLAB仿真软件实现，可利用MATLAB信号处理工具箱提供的各种函数模型实现，可忽略其中的一些复杂的函数变换，从而简化理论设计和论证，具体如下：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; endz,p,k=buttap(N); %调用Butterworth模拟低通滤波器原型函数b,a=zp2tf(

19、z,p,k); %将零点极点增益形式转换为传递函数形式H,w=freqs(b,a,n); %按n指定的频率点给出频率响应magH2=(abs(H).2; subplot(2,1,1);hold on;plot(w,magH2); endxlabel('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

20、9;)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; endz,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内图像用的,如果在新

21、画图像之后不想覆盖原图像就要加上hold on这句话plot(w,abs(Ht);endtitle('模拟高通滤波器');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

22、,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实现对工频

23、干扰信号的滤波，调用plot函数分别画出滤除工频干扰后的时域图和频谱图。具体操作如下：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(

24、y5);title('滤除电源线干扰后的频谱图');xlabel('频率(Hz)');ylabel('幅值(dB)');四 设计结果及分析4.1 心电数据的导入，心电信号的时域图和频谱图如下：图3 原始时域和频谱图通过导入的心电信号数据发现，在其频谱图上的020Hz和80100Hz之间的幅值比较大，而在3070Hz之间的幅值相对较小。4.2 加入高频干扰信号后的波形图通过观察加入高频信号后的时域图和频域图，将它与未加入前进行比较，可以发现频谱图在0Hz时的幅度增加的很大，而且又在没有谱线的频率上出现了频谱，这是由于白噪声在所有频率上都有频率造成的。图4 加入高频干扰后的时域和