基于MATLAB的心电信号分析
摘要：
本课题设计了一个简单的心电信号分析系统。直接采用Matlab语言编程的静态仿真方式、采用Simulink进行动态建模和仿真的方式，对输入的原始心电信号，进行线性插值处理，并通过matlab语言编程设计对其进行时域和频域的波形频谱分析，根据具体设计要求完成系统的程序编写、调试及功能测试。得出一定的结论。
关键字：matlab、心电信号提取、线性插值、滤波、simulink仿真。
一、课题目的及意义
心电信号是人类最早研究并应用于医学临床的生物信号之一，它比其它生物电信号更易于检测，并且具有较直观的规律性，因而心电图分析技术促进了医学的发展。
然而，心电图自动诊断还未广泛应用于临床，从国内外的心电图机检测分析来看，自动分析精度还达不到可以替代医生的水平，仅可以为临床医生提供辅助信息。其主要原因是心电波形的识别不准，并且心电图诊断标准不统一。因此，探索新的方法以提高波形识别的准确率，寻找适合计算机实现又具诊断价值的诊断标准，是改进心电图自动诊断效果，扩大其应用范围的根本途径。如何把心电信号的特征更加精确的提取出来进行自动分析，判断出其异常的类型成了亟待解决的焦点问题。本课题通过matlab语言编程，对原始心电信号进行一定的分析处理。
二、课题任务及要求
1、必做部分
（1）利用Matlab对MIT-BIH数据库提供的数字心电信号进行读取，并还原实际波形。
（2）对原始心电信号做线性插值
（3）对处理前后的心电信号分别做频谱分析
利用Matlab软件对处理前后的心电信号编程显示其频谱，分析比较滤 频谱，得出结论。
（4） Simulink仿真
根据前面的设计，进行基于Simulink的动态仿真设计。实现心电信号 处理。
2、选作部分
（1）只截取大约2.5s，三个周期左右，大约800个采样数据进行分析。（2）60Hz工频陷波器设计
波前后的的分析和
三、设计技术指标
四、设计方案论证
1、必做部分
2、选作部分
五、设计内容及结果分析
1、基于matlab编写的程序如下：
%读取心电信号并转化成数组形式
function [t,Xn]=duquexinhao1(w)
fid=fopen(w);
C=textscan(fid,'%8c %f %*f','headerlines',2);%去除前两行
fclose(fid);
a=C{2};
b=C{1};
k=length(b);
for i=1:k
c(i)=strread(b(i,:),'%*s %f','delimiter',':');
end
c=c';
d=[c,a];
t=d(:,1); %时间
Xn=d(:,2); %幅度
%线性插值
function [t3,Xn3]=xianxingchazhi(t,Xn)
m=max(t);
t3=0:0.001:m;
t3=t3';
Xn3=interp1(t,Xn,t3);
%保存插值前的信号
function baocun1(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);
%保存插值后的信号
function baocun2(t1,Xn1)
fid = fopen('t1.txt','wt');
fprintf(fid,'%g\\n',t1);
fclose(fid);
fid = fopen('Xn1.txt','wt');
fprintf(fid,'%g\\n',Xn1);
fclose(fid);
%画初始信号和即插值后信号频谱
function keshehuatu(t,Xn,t1,Xn1)
f=1000;
T=1/f;
m=1:length(Xn);
k1=length(Xn1);
m1=1:k1;
q=f*m/length(Xn);
q1=f*m1/k1;
subplot(2,2,1)
plot(t,Xn)
title('初始信号时域波形')
subplot(2,2,2)
Y=fft(Xn);
plot(q,abs(Y))
title('初始信号频谱')
subplot(2,2,3)
axis([0,1000,0,1000])
plot(t1,Xn1)
title('插值信号时域波形')
Y1=fft(Xn1);
subplot(2,2,4)
axis([0,1000,0,5000])
plot(q1,abs(Y1))
title('插值信号频谱')
%低通滤波器
function [H,f]=kesheditonglvboqi(wp,ws,Rp,As,Xn1)
T=0.001;
f=1/T;
[N,Wc]=buttord(wp,ws,Rp,As,'s');
[b,a]=butter(N,Wc,'s');
f=(0:length(Xn1)-1)*f/length(Xn1);w=f*2*pi;
H=freqs(b,a,w);
%高通滤波器
function [H,f]=keshegaotonglvboqi(wp,ws,Rp,As,Xn1)
T=0.001;
fs=1/T;
[N,Wc]=buttord(wp,ws,Rp,As,'s');
[b,a]=butter(N,Wc,'high','s');
f=(0:length(Xn1)-1)*fs/length(Xn1);w=f*2*pi;
H=freqs(b,a,w);
%带阻滤波器
function [H,f]=keshedaizulvboqi(wp,ws,p,s,Xn1)
T=0.001;
f=1/T;
[N,Wc]=buttord(wp,ws,p,s,'s');
[b,a]=butter(N,Wc,'stop','s');
f=(0:length(Xn1)-1)*f/length(Xn1);w=f*2*pi;
H=freqs(b,a,w);
主函数如下
（1）、将信号通过低通、高通、带阻滤波器程序
[t,Xn]=duquexinhao1('117.txt');
baocun1(t,Xn) %保存读取信号
[t1,Xn1]=xianxingchazhi(t,Xn);
baocun2(t1,Xn1)%保存插值后信号
xy=[t1,Xn1]; %仿真输入二维数组
figure(1)
keshehuatu(t,Xn,t1,Xn1) %画原始信号和插值后信号波形和频谱
wp=90*2*pi; %低通滤波器滤波
ws=99*2*pi;
p=1;
s=35;
[H1,f]=kesheditonglvboqi(wp,ws,p,s,Xn1);
wp=4*2*pi; %高通滤波器滤波
ws=0.25*2*pi;
p=1;
s=35;
[H2,f]=keshegaotonglvboqi(wp,ws,p,s,Xn1);
wp=[58,62]*2*pi; %带阻滤波器
ws=[59.9,60.1]*2*pi;
[H3,f]=keshedaizulvboqi(wp,ws,p,s,Xn1);
H=abs(H1).*abs(H2).*abs(H3); %低通和高通和带阻组合的滤波器
Y=H'.*abs(fft(Xn1)); %经过滤波后心电信号频谱
y=ifft(Y); %滤波后心电信号时域波形
figure(2)
subplot(2,2,1)
plot(f,abs(H1))
axis([0,150,0,1.5])
title('低通滤波器')
subplot(2,2,2)
plot(f,abs(H2))
axis([0,50,0,1.5])
title('高通滤波器')
subplot(2,2,3)
plot(f,abs(H3))
axis([0,150,0,1.5])
title('带阻滤波器')
subplot(2,2,4)
plot(f,abs(H))
axis([0,100,0,1.5])
title('组合后滤波器')
figure(3)
plot(f,abs(Y))
axis([0,100,0,80])
title('滤波后心电信号频谱')
figure(4)
subplot(2,1,1)
plot(t1,Xn1)
title('滤波前信号')
subplot(2,1,2)
plot(t1,y)
title('滤波后信号')
所出图形如下
结果分析：
（2）、直接通过带通滤波器程序
[t,Xn]=duquexinhao1('117.txt');
baocun1(t,Xn) %保存读取信号
[t1,Xn1]=xianxingchazhi(t,Xn);
baocun2(t1,Xn1)%保存插值后信号
figure(1)
keshehuatu(t,Xn,t1,Xn1) %画原始信号和插值后信号波形和频谱
wp=[2,80]*2*pi;
ws=[0.25,99]*2*pi;
p=1;
s=35;
[H1,f]=kesheditonglvboqi(wp,ws,p,s,Xn1);
H=abs(H1) ; %带通
Y=H'.*abs(fft(Xn1));%经过滤波后心电信号频谱
y=ifft(Y); %滤波后心电信号时域波形
figure(2)
subplot(1,2,1)
plot(f,abs(H1))
axis([0,200,0,1.5])
title('带通滤波器')
subplot(1,2,2)
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,-1.5,1.5])
title('滤波前截取一部分信号')
subplot(2,2,4)
plot(t1,y)
axis([0,1.5,-1.5,1.5])
title('滤波后截取一部分信号')
所出图形如下
结果分析：
（3）、将信号通过低通、高通组合成的带通滤波器程序
[t,Xn]=duquexinhao1('117.txt');
baocun1(t,Xn) %保存读取信号
[t1,Xn1]=xianxingchazhi(t,Xn);
baocun2(t1,Xn1)%保存插值后信号
figure(1)
keshehuatu(t,Xn,t1,Xn1) %画原始信号和插值后信号波形和频谱
xy=[t1,Xn1];
wp=0.52*2*pi; %低通滤波器滤波
ws=0.62*2*pi;
p=1;
s=35;
[H1,f]=kesheditonglvboqi(wp,ws,p,s,Xn1);
wp=0.10*2*pi; %高通滤波器滤波
ws=0.25*2*pi;
p=1;
s=35;
[H2,f]=keshegaotonglvboqi(wp,ws,p,s,Xn1);
H=abs(H1).*abs(H2); %低通和高通组合的带通
Y=H'.*abs(fft(Xn1)); %经过滤波后心电信号频谱
y=ifft(Y); %滤波后心电信号时域波形
figure(2)
subplot(2,2,1)
plot(f,abs(H1))
axis([0,1,0,1.5])
title('低通滤波器')
subplot(2,2,2)
plot(f,abs(H2))
axis([0,1,0,1.5])
title('高通滤波器')
subplot(2,2,3)
plot(f,abs(H))
axis([0,1,0,1.5])
title('组合 带通滤波器')
subplot(2,2,4)
plot(f,abs(Y))
axis([0,1,0,260])
title('滤波后心电信号频谱')
figure(3)
subplot(2,1,1)
plot(t1,Xn1)
title('滤波前信号')
subplot(2,1,2)
plot(t1,y)
title('滤波后信号')
所出图形如下
结果分析：
三种方案比较分析：
（4）系统零极点分析（在此只以高通滤波器为例）
%求高通滤波器的阶数及分子分母系数
wp=0.1*pi;ws=0.25*pi;Rp=1;As=35;T=1;%数字指标
OmegaP=(2/T)*tan(wp/2);%通带模拟频率
OmegaS=(2/T)*tan(ws/2);%阻带模拟频率
[cs,ds]=afd_butt(OmegaP,OmegaS,Rp,As);%归一化巴特沃斯滤波器原型系统函数
N=ceil((log10((10^(Rp/10)-1)/(10^(As/10)-1)))/(2*log10(wp/ws)))
OmegaC=wp/((10^(Rp/10)-1)^(1/(2*N))); %求对应于N的3db截止频率；
[b,a]=u_buttap(N,OmegaC);%去归一化巴特沃斯滤波器原型系统函数
[db,mag,pha,w]=freqz_m(b,a);
subplot(2,1,1);plot(w/pi,mag);
title('digital filter Magnitude Response'); axis([0,1,0,0.01])
subplot(2,1,2);plot(w/pi,db);
title('digital filter Magnitude in DB'); axis([0,1,-30,5]);
%结果：N = 5
% Butterworth Filter Order= 5
%OmegaC = 0.3626
%b = 0.0063
%a = 1.0000 1.1734 0.6884 0.2496 0.0559 0.0063
%N=6
（5）求上述高通滤波器的系统函数及其频谱
b=0.0063;
a=[ 1.0000 1.1734 0.6884 0.2496 0.0559 0.0063];
h=impz(b,a); %系统的单位取样响应
figure(1);plot(h) %画出单位取样响应
title('h(n)')
figure(2)
fs=1000;
[H,f]=freqz(b,a,256,fs); %求出系统的频率响应
mag=abs(H); %幅度响应
ph=angle(H); %相位响应
ph=ph*180/pi;
subplot(2,1,1),plot(f,mag);grid %画出幅度响应
xlabel('frequency(Hz)');ylabel('magnitude');
title('|H(jw)|');
subplot(2,1,2);plot(f,ph);grid %画出相位响应
xlabel('frequency(Hz)');
ylabel('phase');
title('相位');
figure(3)
zr=roots(b) %求出系统的零点
pk=roots(a) %求出系统的极点
zplane(b,a); %zplane函数画出零极点图
%结果：zr = Empty matrix: 0-by-1 %pk = -0.1120 + 0.3438i
% -0.1120 - 0.3438i %-0.3674
%-0.2910 + 0.2156i %-0.2910 - 0.2156i
图形如下
系统函数及级联图：
结果分析：
（7）Simulink仿真：（在此只取第一种方案）
图形如下：
技术指标及结果分析：
选作部分：
1、 只截取大约2.5s，三个周期左右，大约800个采样数据进行分析
程序如下：
function [t,Xn]=duqu2(w)
fid=fopen(w);
C=textscan(fid,'%8c %f %*f','headerlines',2);
fclose(fid);
a=C{2};
b=C{1};
k=length(b);
for i=1:k
c(i)=strread(b(i,:),'%*s %f','delimiter',':');
end
c=c';
d=[c,a];
%截取2.5s的心电信号
for i=1:k
if c(i)<=2.5
e(i,:)=d(i,:);
else break;
end
end
t=e(:,1); %时间
Xn=e(:,2); %幅度
调用程序：
图形如下：
结果分析：
2、60Hz工频陷波器设计：
程序如下：
%60Hz工频陷波器设计
wp=[58,62]*2*pi;
ws=[59.9,60.1]*2*pi;
[H3,f]=keshedaizulvboqi(wp,ws,p,s,Xn1);
plot(f,abs(H3))
axis([0,150,0,1.5])
title('60Hz工频陷波器设计')
图形如下：
分析：
课题总结
附录
一、参考文献
[1] 北京迪阳正泰科技发展公司.综合通信实验系统——信号与系统指导书二版）. 2006,6
[2] 丁玉美.数字信号处理（第二版）.西安电子科技大学出版社，2001
[3] 吴大正. 信号与线性系统分析（第四版）. 高等教育出版社，2005,8
[4] 谢嘉奎. 电子线路--线性部分(第四版). 高等教育出版社，2003,2
[5] 陈后金. 信号分析与处理实验. 高等教育出版社，2006，8
二、 附录——设计原理
1.心电信号的读取
txt格式的数据文件内容及格式如图1-1所示（以100.txt为例）。
（第
图1-1 txt格式的心电数据文件
其中文件的第一列为采样时间，第二列是在以MLII这种导联方式所得到的采样数据，第三列是以V5这种导联方式所得到的采样数据，全文件记录了约为10s的心电数据，3600个采样数据，每一行数据之间用Tab符分隔。
由于数据文件中后两列数据是对同一种心电信号进行不同的导联方式所得到的采样数据，所以可以只采用其中的一种采样数据，摒弃另外一种，即可完成对此心电信号的分析。全部的心电文件记录时间约为10s，共计12个左右周期的心电信号。
实际设计心电信号数据文件时应注意：
2.心电信号的线性插值处理
根据上文中提到的插值公式，以此为原理，设计Matlab程序，对心电信号数据做线性插值处理。插值完以后的数据应该是时间均匀的、以0.001秒为间隔的。
此步骤的实现主要是基于Matlab中的数组操作函数来实现，建议一定首先熟悉并掌握Matlab中的所有数组操作函数的作用和操作方法。
其中一种插值方法的思路是：将第一步中读取的心电信号数据的时间数据和幅值数据分别存放在一个一维数组中。然后利用for循环结构把所有数据依次读取进来。判断时间数据数组中前后两个相邻的数据间隔是否为0.001s，如果是则判断下一对相邻两个数据；如果间隔大于0.001s则进行一维插值处理。
注意对时间数据做插值的同时一定不要忘记对幅值数据同样做插值处理，时间数据和幅值数据一定是相互对应的。
3. 滤波器设计
3.1 模拟滤波器设计原理
（1）模拟巴特沃思滤波器原理
巴特沃斯滤波器具有单调下降的幅频特性：在小于截止频率
的范围内，具有最平幅度的响应，而在
后，幅频响应迅速下降。
巴特沃思低通滤波器幅度平方函数为：
(3-1)
式中N为滤波器阶数，
为3dB截止角频率。将幅度平方函数写成s的函数形式：
(3-2)
该幅度平方函数有2N个等间隔分布在半径为
的圆上的极点
，
为了形成稳定的滤波器，取左半平面的N个极点构成
，即：
(3-3)
为使设计统一，将频率归一化，得到归一化极点
，相应的归一化系统函数为：
(3-4)
多项式形式为：
(3-5)
（2）模拟切比雪夫滤波器原理
切比雪夫滤波器的幅频特性具有等波纹特性，有两种形式，在通带内等波纹、阻带单调的是
型滤波器，在通带内单调、在阻带内等波纹的是Ⅱ滤波器。以
型滤波器为例。
切比雪夫滤波器的幅度平方函数为：
(3-6)
ε为小于1的正数，表示通带内幅度波动的程度。Ωp称为通带截止频率。令λ=Ω/Ωp，称为对Ωp的归一化频率。CN(x)为N阶切比雪夫多项式。幅度平方函数的极点是分布在bΩp为长半轴，aΩp为短半轴的椭圆上的点。同样取s平面左半平面的极点构成
：
(3-7)
进行归一化，得到：
其中
，
3.2 模拟滤波器数字化原理
将模拟滤波器转化为数字滤波器在工程上常用的有脉冲响应不变法和双线性变换法。
脉冲响应不变法是一种时域上的转换方法，它是数字滤波器的单位取样响应在采样点上等于模拟滤波器的单位冲激响应，即：
(3-8)
设模拟滤波器只有单阶极点，其系统函数为：
(3-9)
对
进行拉氏反变换得到
，对
进行等间隔采样，得到，对
进行Z变换，得到数字滤波器系统函数：
(3-10)
这种方法s和z的关系是：
。该方法的优点是频率坐标变换时线性的切数字滤波器的单位脉冲响应完全模仿模拟滤波器的单位冲激响应，时域特性逼近好；缺点是会产生频谱混叠现象，适合低通、带通滤波器的设计，不适合高通、带阻滤波器的设计。
双线性变换法为了克服频谱混叠现象，采用非线性频率压缩方法，将整个频率轴上的频率范围压缩到
之间，再用
转换到Z平面上。
这种方法s和z的关系是：
。该方法克服了频谱混叠现象，但带来了频率坐标变换的非线性：
，由模拟滤波器系统函数转换为数字滤波器系统函数公式为：
(3-11)
3.3数字高通、带通、带阻滤波器的设计
上述滤波器可以借助于模拟滤波器的频率变换设计一个所需类型的模拟滤波器，通过双线性变换法将其转换成所需类型的数字滤波器。
首先确定所需类型数字滤波器的技术指标；然后将数字滤波器技术指标按照公式
再
转换成所需类型滤波器的模拟域技术指标；将所需类型滤波器的模拟域技术指标转换成低通滤波器技术指标；设计归一化模拟低通滤波器；去归一化得到模拟低通滤波器的系统函数；将模拟低通滤波器转换为所需类型的模拟滤波器；最后通过双线性变换法转换成所需类型的数字滤波器。
课程设计心得体会
因篇幅问题不能全部显示，请点此查看更多更全内容