数字信号处理选做实验--心电信号数字滤波系统设计

1、实验题目心电信号数字滤波系统设计实验时间2013年12 月16日实验地点A1118实验成绩实验性质验证性 设计性 综合性教师评语： 出勤率好 原理正确 方案合理 实验结果正确 回答问题正确 报告规范一、实验目的本实验项目主要研究数字心电信号的初步分析及滤波器的应用。通过完成本实验的设计，主要达到以下几个目的：1、了解MATLAB软件的特点和使用方法，熟悉基于Simulink的动态建模和仿真的步骤和过程。2、了解人体心电信号的时域特征和频谱特征。3、进一步了解数字信号的分析方法。 4、通过应用具体的滤波器进一步加深对滤波器理解。 5、通过本课题的设计，培养学生运用所学知识分析和解决实际问题的能力

2、。二、实验原理本文是利用MATLAB软件对美国麻省理工学院提供的MIT-BIH数据库的122号心电信号病例进行分析，利用MATLAB软件及simulink平台对122号心电信号的病例进行读取、插值、高通滤波、低通滤波等的处理。将心电信号中的高频和低频的杂波进行滤除后对插值前后滤波前后的时域波形及频谱进行分析。同时也将滤波器的系统函数进行读取，分析，画出滤波的信号流程图，也分析各个系统及级联后的系统的幅频响应、相位响应来判断系统的稳定性，并用MATLAB软件将图形画出，以便于以后的对系统进行分析。三、使用仪器、器材Matlab软件四、实验步骤1.对原始心电信号做线性插值 运算公式如下: 其中 是

3、第i个数据时间点，Ai是与之对应的数据，N是两数据之间需要的插值数， 是需要插值的两点数据差， ，2.根据心电信号的频域特征，设计相应的低通和带通滤波器 一般正常人的心电信号频率在0.7100HZ范围内，人体心电信号微弱，信噪比小，因此，在采集心电信号时，易受到仪器、人体活动等因素的影响，而且所采集的心电信号常伴有干扰。采集心电数据时，由于人的说话呼吸，常常会混有约为0.1Hz到0.25Hz频段的干扰，对于这些低频干扰，可以让信号通过一个高频滤波器，低截止频率设置为0.25,来滤除低频信号，对于高频信号干扰，可以让信号再通过一个低频滤波器，其中截止频率设置为99H3.对处理前后的心电信号分别做

4、频谱分析 利用Matlab软件对处理前后的心电信号编程显示其频谱，分析比较滤波前后的频谱，得出结论。如果分析频谱，滤波效果不明显，则需变动滤波器参数指标，重新设计滤波器。通过频谱分析，多次试验确定最合适的滤波器。4 .Simulink仿真. 根据前面的设计，进行基于Simulink的动态仿真设计。实现心电信号的分析和处理。给出系统的基于Simulink的动态建模和仿真的系统方框图，同时记录系统的各个输出点的波形和频谱图。 五、实验过程原始记录（数据、图标、计算等）1. 提取txt格式心电信号：fid=fopen(122.txt);C=textscan(fid,%8c %f %*f,header

5、lines,2);f close(fid);a=C1;y=C2; k=length(a)for i=1:k c(i)=strread(a(i,:),%*s %f,delimiter,:);endx=c;plot(x,y)2. 对原始心电信号进行线性插值t=0.006:0.001:2.5;F=interp1(x,y,t);F=F;t=t;plot(t,F) 3. 把数据读到txt中fid = fopen(t.txt,wt);fprintf(fid,%gn,t); fclose(fid);fid = fopen(F.txt,wt);fprintf(fid,%gn,F); fclose(fid);4

6、. 插值前后波形比较subplot(2,2,1)plot(x,y)title(初始信号时域波形)axis(0 2.5 -2 1)subplot(2,2,2)fs=1000;N=length(y)n=1:N;f1=n*fs/N;Y1=fft(y);plot(f1,abs(Y1)title(初始信号频谱)axis(0 1000 0 200)subplot(2,2,3)plot(t,F)title(差值后信号时域波形)axis(0 2.5 -2 1)M=length(F);m=1:M;f 2=m*fs/M;Y2=fft(F);subplot(2,2,4)plot(f2,abs(Y2)title(插值

7、后信号频谱)axis(0 1000 0 200) 5 模拟低通滤波器：wp=60*2*pi;ws=99*2*pi;Rp=1;As=40;N,wc=buttord(wp,ws,Rp,As,s)B,A=butter(N,wc,s)k=0:511;fk=0:1000/512:1000;wk=2*pi*fk;Hk=freqs(B,A,wk);plot(fk,20*log10(abs(Hk); grid on Untitled5N = 11wc = 409.2596B = 1.0e+28 * Columns 1 through 10 0 0 0 0 0 0 0 0 0 0 Columns 11 thro

8、ugh 12 0 5.3949A = 1.0e+28 * Columns 1 through 10 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0008 Columns 11 through 12 0.0926 5.39496模拟高通滤波器：wp=0.7*2*pi;ws=0.25*2*pi;Rp=0.1;As=40;N,wc=buttord(wp,ws,Rp,As,s)B0,A0=butter(N,wc,s);wph=2*pi*0.25;hk=freqs(B0,A0,wph);BH,AH=lp2hp(B,A,

9、wph);h,w=freqs(BH,AH);plot(w,20*log10(abs(h);axis(0,1,-80,5);grid on N = 7wc = 3.0327B0 = 1.0e+003 * 0 0 0 0 0 0 0 2.3595A0 = 1.0e+003 * 0.0010 0.0136 0.0929 0.4070 1.2343 2.5905 3.4964 2.35957滤波前后图形对比：t=0.006:0.001:59.997;F=interp1(x,y,t);F=F;t=t;figure(1)subplot(3,1,1);plot(1000*t,F);wp=0.7*2*pi;w

10、s=0.25*2*pi;Rp=0.1;As=40;T=1;N,wc=buttord(wp,ws,Rp,As,s)B,A=butter(N,wc,s);b,a=imp_invr(B,A,T)db,mag,pha,w=freqz_m(b,a);y1=filter(b,a,F);subplot(3,1,2);plot(y1);title(高通滤波后)wp1=2*pi*60;ws1=2*pi*99;Rp1=0.1;As1=40;T1=1000;OmegaP1=wp1/T1;OmegaS1=ws1/T1; cs1,ds1=afd_butt(OmegaP1,OmegaS1,Rp1,As1);b1,a1=i

11、mp_invr(cs1,ds1,T)db1,mag1,pha1,w1=freqz_m(b1,a1);y2=filter(b1,a1,y1);subplot(3,1,3);plot(y2);title(低通滤波后)M=length(F);m=1:M;fs=1000;f2=m*fs/M;F1=fft(F);Y1=fft(y1);Y2=fft(y2)figure(2)subplot(3,1,1)plot(f2,abs(F1)axis(0,1000,0,200)title(原始信号频谱_59.997)subplot(3,1,2)plot(f2,abs(Y1)axis(0,1000,0,200)titl

12、e(高通滤波后信号频谱_59.997)subplot(3,1,3)plot(f2,abs(Y2)axis(0,1000,0,200)title(低通滤波后信号频谱_59.997)8.截取到2.5s对截取的部分进行滤波及频谱分析t1=0.006:0.001:2.5;F0=interp1(x,y,t1);F0=F0;t1=t1;figure(8)subplot(3,1,1);plot(1000*t1,F0);wp=0.7*2*pi;ws=0.25*2*pi;Rp=0.1;As=40;T=1;N,wc=buttord(wp,ws,Rp,As,s)B,A=butter(N,wc,s);b,a=imp_

13、invr(B,A,T)db,mag,pha,w=freqz_m(b,a);y11=filter(b,a,F0);subplot(3,1,2);plot(y11);title(高通滤波后_2.5)wp1=2*pi*60;ws1=2*pi*99;Rp1=0.1;As1=40;T1=1000;OmegaP1=wp1/T1;OmegaS1=ws1/T1;cs1,ds1=afd_butt(OmegaP1,OmegaS1,Rp1,As1);b1,a1=imp_invr(cs1,ds1,T)db1,mag1,pha1,w1=freqz_m(b1,a1);y21=filter(b1,a1,y11);subplot(3,1,3);plot(y21);title(低通滤波后_2.5)M=length(F0);m=1:M;fs=1000;f2=m*fs/M;F01=fft(F0);Y11=fft(y11);Y21=fft(y21)figure(9)subplot(3,1,1)plot(f2