心电信号处理昆明理工大学.doc

昆明理工大学信息工程与自动化学院学生实验报告
（ 2016 —20 17 学年第二学期）
课程名称：生物医学信号处理 开课实验室： 设备编号：实验日期：2018年5月30日
年级、专业、班

学号

姓名

成绩

实验项目名称
心电信号处理
指导教师
相艳
教师评语

教师签名：
年月日
一、实验目的
1、对心电信号的记录、处理、心电信号的特点、心电信号的噪声抑制，工频干扰的抑制与基线纠漂有总体了解。
2、能利用MATLAB GUI设计简单的GUI程序。
二、实验原理
1、心电信号的特点：
心电信号具有近场检测的特点，离开人体表微小的距离，就基本上检测不到信号；心电信号通常比较微弱，至多为mV量级，且能量主要在几百赫兹以下；干扰即来自生物体内，如肌电干扰、呼吸干扰等，也来自生物体外，如工频干扰，信号拾取时因不良接地等引入的其他外来干扰等；干扰信号与心电信号本身频带重叠（如工频干扰等）。
2、工频干扰抑制：现在使用较多的方法是使用滤波器对工频干扰进行抑制。
3、基线漂移：基线漂移是因呼吸、肢体活动或运动心电图测试所引起的，故这样使得ECG信号的基准线呈现上下飘逸的情况。
三、实验内容
1、对心电信号处理
主程序：
clear;closeall;clc;
load 100_ECG_0_20 //加载心电信号
%%%Eliminate Baseline Drift //消除基线漂移
s1=ECG_2; //把心电信号ECG-2赋给s1
s2=smooth(s1,150); //利用移动平均法对s1做平滑处理
ecgsmooth=s1-s2; //消除基线漂移
%%%apply Wavelet Transform //应用小波变换
[C,L]=wavedec(ecgsmooth,8,'db4'); //用db4对ecgsmooth进行8层分解，其中返回的近似和细节都存放在C中，L存放是近似和各阶细节系数对应的长度（阶数为4阶）
[d1,d2,d3,d4,d5,d6,d7,d8]=detcoef(C,L,[1,2,3,4,5,6,7,8]); //提取小波的细节系数
%%%Denoise //降噪，消除干扰
[thr,sorh,keepapp]=ddencmp('den','wv',ecgsmooth); //返回小波除噪和压缩后的信号
cleanecg=wdencmp('gbl',C,L,'db4',8,thr,sorh,keepapp);//通过门限阈值处理得到小波系数（执行降噪操作）
%%%thresholding1 //取阈值
max_value=max(cleanecg); //最大值（波峰）
mean_value=mean(cleanecg); //最小值（波谷）
threshold=(max_value-mean_value)/2; //最大值与最小值差的一半作为阈值
%%%R detection algorithm//用R检测算法检测信号
a5=appcoef(C,L,'db4',5);//取分解后的近似部分，也就是第5层低频系数
C1=[a5;d5;d4;d3]; //
L1=[length(a5);length(d5);length(d4);length(d3);length(cleanecg)];
R_detect_signal=waverec(C1,L1,'db4'); //用二维小波分解的结果C1，L1重建信号
R_detect_squared=R_detect_signal.^2; //对R检测信号求平方
%%%Beat_Rate_Extraction_Algorithm //计算心率
for a=1:length(R_detect_squared)
ifR_detect_squared(a)>threshold
R_detect_new(a)=R_detect_squared(a);
Else
R_detect_new(a)=0;
end
end
mean_R_detect=5*mean(R_detect_new);
for q=1:length( R_detect_new)-1
ifR_detect_new(q)<mean_R_detect
R_detect_new(q)=0;
end
end
d=0;
for b=1:length( R_detect_new)-1
if ( R_detect_new(b)==0) & ( R_detect_new(b+1)~=0)
d=d+1;
indext(d)= b+1;
end
end
fs_R_deetect=length(R_detect_new)/20;
time=indext.*1/fs_R_deetect;
ind=0;
for z=1:length(time)-1
ind=ind+1;
time_diff(ind)=time(z+1)-time(z);
end
av_time=mean(time_diff);
Square_Number=av_time/.2;
beat_Rate=300/Square_Number;
high=max(R_detect_new);
subplot(411);plot(s1);title('Orginal Signal'); //绘制原信号波形图
subplot(412);plot(s1-s2);title('Baseline drift Elimination'); //绘制消除基线漂移波形图
subplot(413);plot(cleanecg);title('Main Signal'); //绘制消除噪声之后的波形图
subplot(414);plot(R_detect_new);title('R detected Signal'); //绘制R检测算法检测信号的波形图
text(length(R_detect_new)/2,high,['Beat Rate = ',num2str(fix(beat_Rate))],'EdgeColor','red'); //显示心率

图1 运行结果
在此信号处理过程中用到的处理算法有小波分解以及小波信号重构，用了Daubechies(dbN)小波，Daubechies小波简写为dbN，N是小波的阶数。此程序中db4,表示小波阶数为4。其流程图下：

图2小波算法流程图
2、 利用matlab GUI平台设计相关GUI程序
① GUI界面设置

图3 GUI界面设置图
② GUI程序：
functionvarargout = untitled1(varargin)
gui_Singleton = 1;
gui_State = struct('gui_Name', mfilename, ...
'gui_Singleton', gui_Singleton, ...
'gui_OpeningFcn', @untitled1_OpeningFcn, ...
'gui_OutputFcn', @untitled1_OutputFcn, ...
'gui_LayoutFcn', [] , ...
'gui_Callback', []);
ifnargin&&ischar(varargin{1})
gui_State.gui_Callback = str2func(varargin{1});
end
ifnargout
[varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
gui_mainfcn(gui_State, varargin{:});
end
function untitled1_OpeningFcn(hObject, eventdata, handles, varargin)
handles.output = hObject;
guidata(hObject, handles);

functionvarargout = untitled1_OutputFcn(hObject, eventdata, handles)
varargout{1} = handles.output;
%%设置下拉菜单
function popupmenu1_Callback(hObject, eventdata, handles)load 100_ECG_0_20
a=get(handles.popupmenu1,'value');
switch a
case 1
s1=ECG_1;
s2=smooth(s1,150);
ecgsmooth=s1-s2;
case 2
s1=ECG_2;
s2=smooth(s1,150);
ecgsmooth=s1-s2;
End
%%绘制原信号波形图，并将图线设置为红色
function popupmenu1_CreateFcn(hObject, eventdata, handles)
ifispc&&isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
set(hObject,'BackgroundColor','white');
end
function pushbutton1_Callback(hObject, eventdata, handles)
load100_ECG_0_20
a=get( handles.popupmenu1,'value')
switch a
case 1
s1=ECG_1;
case 2
s1=ECG_2;
end
subplot(411);plot(s1,'--r');title('Orginal Signal');
%%绘制基线漂移的波形图，图线设置为绿色
function pushbutton2_Callback(hObject, eventdata, handles)
load100_ECG_0_20
a=get(handles.popupmenu1,'value');
switch a
case 1
s1=ECG_1;
s2=smooth(s1,150);
case 2
s1=ECG_2;
s2=smooth(s1,150);
end
subplot(412);plot(s1-s2,'--g');title('Baseline drift Elimination');
%%绘制进行小波变换降噪之后的主信号波形图
function pushbutton3_Callback(hObject, eventdata, handles)
load100_ECG_0_20
a=get( handles.popupmenu1,'value')
switch a
case 1
s1=ECG_1;
s2=smooth(s1,150);
ecgsmooth=s1-s2;
case 2
s1=ECG_2;
s2=smooth(s1,150);
ecgsmooth=s1-s2;
end
[C,L]=wavedec(ecgsmooth,8,'db4');
[d1,d2,d3,d4,d5,d6,d7,d8]=detcoef(C,L,[1,2,3,4,5,6,7,8]);
[thr,sorh,keepapp]=ddencmp('den','wv',ecgsmooth);
cleanecg=wdencmp('gbl',C,L,'db4',8,thr,sorh,keepapp);
subplot(413);plot(cleanecg);title('Main Signal');
%%绘制进行R检测之后的波形图
function pushbutton4_Callback(hObject, eventdata, handles)
load100_ECG_0_20
a=get( handles.popupmenu1,'value')
switch a
case 1
s1=ECG_1;
s2=smooth(s1,150);
ecgsmooth=s1-s2;
case 2
s1=ECG_2;
s2=smooth(s1,150);
ecgsmooth=s1-s2;
end
[C,L]=wavedec(ecgsmooth,8,'db4');
[d1,d2,d3,d4,d5,d6,d7,d8]=detcoef(C,L,[1,2,3,4,5,6,7,8]);
[thr,sorh,keepapp]=ddencmp('den','wv',ecgsmooth);
cleanecg=wdencmp('gbl',C,L,'db4',8,thr,sorh,keepapp);
max_value=max(cleanecg);
mean_value=mean(cleanecg);
threshold=(max_value-mean_value)/2;
a5=appcoef(C,L,'db4',5);
C1=[a5;d5;d4;d3];
L1=[length(a5);length(d5);length(d4);length(d3);length(cleanecg)];
R_detect_signal=waverec(C1,L1,'db4');
R_detect_squared=R_detect_signal.^2;

for a=1:length(R_detect_squared)
ifR_detect_squared(a)>threshold
R_detect_new(a)=R_detect_squared(a);
else
R_detect_new(a)=0;
end
end
mean_R_detect=5*mean(R_detect_new);
for q=1:length( R_detect_new)-1
ifR_detect_new(q)<mean_R_detect
R_detect_new(q)=0;
end
end
d=0;
for b=1:length( R_detect_new)-1
if ( R_detect_new(b)==0) & ( R_detect_new(b+1)~=0)
d=d+1;
indext(d)= b+1;
end
end
fs_R_deetect=length(R_detect_new)/20;
time=indext.*1/fs_R_deetect;
ind=0;
for z=1:length(time)-1
ind=ind+1;
time_diff(ind)=time(z+1)-time(z);
end
av_time=mean(time_diff);
Square_Number=av_time/.2;
beat_Rate=300/Square_Number;
high=max(R_detect_new);
subplot(414);plot(R_detect_new);title('R detected Signal');
function pushbutton7_Callback(hObject, eventdata, handles)
globalR_detect_new
d=0;
for b=1:length( R_detect_new)-1
if ( R_detect_new(b)==0) & ( R_detect_new(b+1)~=0)
d=d+1;
indext(d)= b+1;
end
end
fs_R_deetect=length(R_detect_new)/20;
time=indext.*1/fs_R_deetect;
ind=0;
for z=1:length(time)-1
ind=ind+1;
time_diff(ind)=time(z+1)-time(z);
end
av_time=mean(time_diff);
Square_Number=av_time/.2;
beat_Rate=300/Square_Number;
high=max(R_detect_new);
set(handles.edit6,'string',num2str(fix(beat_Rate))); //文本编辑框中显示心率
function edit6_Callback(hObject, eventdata, handles)
function edit6_CreateFcn(hObject, eventdata, handles)
ifispc&&isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
set(hObject,'BackgroundColor','white');
end

图4 GUI界面运行结果图