数字信号处理课程设计二心电信号分析系统设计-matlab

数字信号处理课程设计报告课题名称：心电信号分析系统的设计–matlab班级: 通信103姓名：滕冲学号： 102189成绩：指导教师：王宝珠日期： 2011年12月28日目录摘要 (3)一、课程设计的目的、意义 (3)二、设计任务及技术指标 (3)三、设计方案论证 (3)三、相关函数及程 (3)（1）读取txt形式的心电信号M文件 (3)(2) 数据保存：分别用来保存插值前和插值后的数据 (4)(3)线性插值： (4)(4)对以上步骤进行画图 (4)(5) 两个滤波器 (5)（6）主程序 (5)五、实验结果与分析 (10)六、结论（只对设计课题内容进行总结并给出结论） (14)七、课程设计心得与体会 (15)八、参考文献 (16)摘要：设计一个心电信号分析系统。

读取原始心电信号txt文件并将其前两行内容用程序的方法去除，对其做一次线性插值并保存在新的txt文件里面。

由心电信号特点得知需要设计带通滤波器还有一个带阻滤波器，设置初步滤波技术指标画出初步滤波图形后看出的结果好不好，出图是否与设计指标相符，若计数指标不好应该怎样调节，对处理完的时域和频域波形进行分析，得出结论，做过程记录。

该系统做完后，启用simulink功能，建立模型，进行连线，设置参数，最后启动仿真。

观察示波器模块的输入和输出波形是否真确。

关键词：数字心电信号，线性插值，滤波器，功能仿真一、课程设计的目的、意义通过本课题的设计学会用matlab软件实现对心电信号进行读取，设计滤波器并进行滤波，分析心电信号的特点并反馈给滤波器的设计指标，同时强化对matlab语句应用的训练更好地掌握该软件；提高对实际问题的分析解决能力，锻炼综合能力；通过实际的动手过程发现问题解决问题，加深对所学知识的理解和掌握，重在过程。

二、设计任务及技术指标(1)对原始数字心电信号进行读取，由数字信号数据绘制出其时域波形。

(2)对数字信号数据做一次线性插值，使其成为均匀数字信号，以便后面的信号分析。

(3)根据心电信号的频域特征（自己查阅相关资料），设计相应的低通和高通滤波器。

(4)编程绘制实现信号处理前后的频谱，做频谱分析，得出相关结论。

(5)对系统进行综合测试，整理数据，撰写设计报告。

三、设计方案论证（1）采用线性插值函数interp1完成插值处理（2）设计一个带通滤波器，经画图分析，需加上一个带阻滤波器，滤除60Hz处的干扰信号（3）初始技术指标：带通：wp=[1,90]*2*pi（滤波器通带截止频率）;ws=[0.25,99]*2*pi（阻带截止频率）; Rp=1db（通带衰减）； As=25db（阻带衰减）；带阻：wp1=[58,62]*2*pi;（滤波器通带截止频率）;ws1=[59.9,60.1]*2*pi;（阻带截止频率）; Rp=1db（通带衰减）； As=25db（阻带衰减）；四、相关函数及程序（1）读取txt形式的心电信号M文件%读取心电信号并转化成数组形式function [t,Xn]=duquxinhao(w)fid=fopen(w);A=textscan(fid,'%8c %f %*f','headerlines',2); %将其前两行内容去除fclose(fid);a=A{1};b=A{2};k=length(a)for i=1:kc(i)=strread(a(i,:),'%*s %f','delimiter',':');endc=c'd=[c,b];%为了减少分析数据的工作量截取2.5s(三个周期左右)的心电信号for i=1:kif c(i)<=2.5e(i,:)=d(i,:);else break;endendt=e(:,1); %时间Xn=e(:,2); %幅度（2）数据保存：分别用来保存插值前和插值后的数据function shujubaocun(t,Xn)fid = fopen('t.txt','wt');fprintf(fid,'%g\n',t);fclose(fid);fid = fopen('Xn.txt','wt');fprintf(fid,'%g\n',Xn);fclose(fid);（3）线性插值：function [t3,Xn3]=xianxingchazhi(t,Xn)m=max(t);t3=0:0.001:m;t3=t3';Xn3=interp1(t,Xn,t3);（4）对以上步骤进行画图：%画初始信号和即插值后信号频谱function huizhiboxing(t,Xn,t1,Xn1)fs=1000;T=1/fs;n1=length(Xn);n=1:n1;m1=length(Xn1);m=1:m1;q=fs*n/n1;q1=fs*m/m1;subplot(2,2,1)plot(t,Xn)title('初始信号时域波形')subplot(2,2,2)Y=fft(Xn);plot(q,abs(Y))axis([0 100 0 200]);title('初始信号频谱')subplot(2,2,3)plot(t1,Xn1)title('插值后信号时域波形')Y1=fft(Xn1);subplot(2,2,4)plot(q1,abs(Y1))title('均匀信号频谱')axis([0 100 0 200]);（5）两个滤波器：%带通滤波器function [H,f]=daitonglvboqi(wp,ws,Rp,As,Xn2) T=0.001;f=1/T;[N,Wc]=buttord(wp,ws,Rp,As,'s');[b,a]=butter(N,Wc,'s');%[b0,B,A]=dir2cas(b,a)f=(0:length(Xn2)-1)*f/length(Xn2);w=f*2*pi; H=freqs(b,a,w);%plot(w/(2*pi),20*log10(abs(H)));%xlabel('频率Hz')%ylabel('db')%title(['巴特沃斯N=',num2str(N)])%grid on%带阻滤波器function [H,f]=daizulvboqi(wp,ws,p,s,Xn2)T=0.001;fs=1/T;[N,Wc]=buttord(wp,ws,p,s,'s');[b,a]=butter(N,Wc,'stop','s');f=(0:length(Xn2)-1)*fs/length(Xn2);w=f*2*pi; H=freqs(b,a,w);（6）主程序：[t,Xn]=duquxinhao('113.txt')shujubaocun(t,Xn) %保存读取信号[t1,Xn1]=xianxingchazhi(t,Xn);xin=[t1,Xn1]; %用二维数组读取信号源shujubaocun(t1,Xn1) %保存插值后信号figure(1)huizhiboxing(t,Xn,t1,Xn1) %画原始信号和插值后信号波形和频谱wp=[5,80]*2*pi;ws=[0.25,99]*2*pi;wp1=[58,62]*2*pi;ws1=[59.9,60.1]*2*pi;Rp=1;As=20;[H1,f]=daitonglvboqi(wp,ws,Rp,As,Xn1);[H2,f]=daizulvboqi(wp1,ws1,Rp,As,Xn1);H=abs(H1) ;%带通HH=abs(H2)%带阻Y=(H.*HH)'.*abs(fft(Xn1));%经过滤波后心电信号频谱y=ifft(Y); %滤波后心电信号时域波形figure(2)subplot(3,1,1)plot(f,abs(H1))axis([-20,140,0,1.5])title('带通滤波器')subplot(3,1,2)plot(f,abs(HH))axis([0,200,0,1.5])title('带阻滤波器')subplot(3,1,3)plot(f,abs(Y))axis([0,100,0,80])title('滤波后心电信号频谱')figure(3)subplot(2,2,1)plot(t1,Xn1)title('滤波前信号')subplot(2,2,2)plot(t1,y)title('滤波后信号')subplot(2,2,3)plot(t1,Xn1)axis([0,1.5,0,1.5])title('滤波前截取一部分信号')subplot(2,2,4)plot(t1,y)axis([0,1.5,0,1.5])title('滤波后截取一部分信号')% 分析带通滤波器figure(4)wp=[5,80]*pi*2 %模拟通带角频率ws=[0.25,99]*pi*2 %模拟阻带角频率Rp=1,As=5%通带波动、阻带衰减%巴特沃斯型模拟带通滤波器设计[N,wc]=buttord(wp,ws,Rp,As,'s')%计算滤波器阶数N和3db截止频率wc[b,a]=butter(N,wc,'bandpass','s') %计算模拟带通滤波器函数系数%巴特沃斯型模拟带通滤波器频率响应[H,w]=freqs(b,a)mag=abs(H);[bn,an]=impinvar(b,a,1000) %用冲击响应不变法将模拟滤波器转换成数字滤波器；数字带通滤波器系数%巴特沃斯型数字带通滤波器频率响应ht=impz(bn,an) %数字滤波器冲击响应曲线subplot(3,1,1)plot(w/(2*pi),mag)%模拟带通滤波器的幅频响应axis([-50,150,0,1.5])title('模拟带通滤波器幅频曲线')grid %模拟幅频曲线subplot(3,1,2);plot(ht);title('数字带通滤波器冲击响应曲线')axis([0,100,-0.5,0.5]);gridsubplot(3,1,3);zplane(bn,an);axis([-2,2,-1.5,1.5])title('数字带通滤波器的零极点图')%zplane函数画出零极点图grid% 分析带阻滤波器figure(5)wp1=[58,62]*2*pi;%模拟通带角频率ws1=[59.9,60.1]*2*pi; %模拟阻带角频率Rp=1,As=20%通带波动、阻带衰减%巴特沃斯型模拟带通滤波器设计[N,wc]=buttord(wp1,ws1,Rp,As,'s')%计算滤波器阶数N和3db截止频率wc[b,a]=butter(N,wc,'stop','s') %计算模拟带通滤波器函数系数%巴特沃斯型模拟带通滤波器频率响应[H,w]=freqs(b,a)mag=abs(H);[bn,an]=impinvar(b,a,1000) %用冲击响应不变法将模拟滤波器转换成数字滤波器；数字带阻滤波器系数ht=impz(bn,an) %数字滤波器冲击响应曲线subplot(3,1,1)plot(w/(2*pi),mag)%模拟带阻滤波器的幅频响应axis([0,150,0,1.5])title('模拟带阻滤波器幅频曲线')grid %模拟幅频曲线subplot(3,1,2);plot(ht);title('数字带阻滤波器冲击响应曲线')axis([0,600,-0.05,0.05]);gridsubplot(3,1,3);zplane(bn,an);title('数字带阻滤波器的零极点图')%zplane函数画出零极点图grid五、实验结果与分析分析：由图可以看出，经过线性插值后的时域信号和均匀信号频谱图，在均匀信号频谱中可以看出在以60Hz为中心的窄带处存电源线干扰，所以在后面添加限制60Hz频率通过的带阻滤波器。