基于MATLAB 的脉搏信号处理软件系统

基于MATLAB 的脉搏信号处理软件系统
摘要: 本文根据在实验室里测得的脉搏数据，基于MATLBA设计一个脉搏信号的GUI处理界面，并利用MATLAB强大数字信号处理功能还原脉搏波形，并对波形的特征信息进行提取及存储。

原始信号进行了去除基线漂移、通过巴特沃斯带通滤波器以及二阶切比雪夫滤波器去除50HZ工频干扰，并且能计算实时的脉率并更新，显示脉率变化趋势曲线，进行频谱分析和输出文档。

此软件有两个GUI界面，第一个为密码登陆界面，第二个为脉搏信号处理系统GUI界面。

第二个GUI界面主要分为五大模块:1.打开与退出模块包括打开数据和退出系统;2.信号回放模块包括对原信号和滤波信号的回放、暂停回放、继续回放、关闭窗口;3.信号放大与缩小模块包括对信号的X轴和Y轴的放大、缩小处理;信号快进退模块包括对信号的快进、慢进、快退、慢退处理;4.脉率实时处理模块包括输出脉率曲线、暂停回放、输出脉搏信息、脉搏频谱分析、清除波形、输出文档;5.脉率信号输出模块包括输出实时的脉率更新、以及脉搏数据的信息，诸如脉搏采样频率、采样时间、最大脉率值、最小脉率等。

关键词：脉搏;脉率;Matlab ;GUI ;
1 引言
人体内部各个生理系统之间(如循环系统、呼吸系统等)是相互耦合的。

反映人身体健康状态相对最重要、最全面的是心脏血液循环系统,因此通过采集脉搏波进而分析心脏循环系统功能,能从一个方面较全面反映人体的健康情况。

从脉搏波中提取人体的生理病理信息作为临床诊断和治疗的依据，历来都受到中外医学界的重视。

几乎世界上所有的民族都用过“摸脉”作为诊断疾病的手段。

脉搏波所呈现出的形态(波形)、强度(波幅)、速率(波速)和节律(周期)等方面的综合信息，在很大程度上反映出人体心血管系统中许多生理病理的血流特征,因此对脉搏波采集和处理具有很高的医学价值和应用前景。

目前脉搏信息的研究已经应用于以下几个方面:(1)中医脉象信息的检测与识别;(2)血压的临床检测;(3)心率稳定性的一种简便估计方法;(4)心输出量的一种测量方法;(5)血管功能的一种早期、无创检测方法。

MATLAB(Matrix Laboratory,矩阵实验室)是由美国MathWorks公司开发的一种功能强、效率高、简单易学的可视化软件,覆盖面包括控制、通讯、金融、图像处理、建筑、生物学等几乎所有的行业与科学领域。

除了经典的一些算法外,MATLAB
还提供了丰富的数据分析和处理功能模块,如神经网络、小波分析、信号处理、图
像处理、自动控制、模糊控制、系统仿真等,因此MATLAB是一种高效的编程软件。

本文介绍利用MATLAB软件作为技术平台,实现对统脉搏波快速、准确实时显示,而
且实现方法简单有效,有一定的实用性。

2 软件总体界面与总体设计思路
2.1密码登陆系统
图2.1 未运行前密码登陆界面
图2.2 运行后输入错误的密码登陆界面
图2.3 运行后输入正确的密码登陆界面
2.2 脉搏信号处理系统界面
图2.4 未运行前的脉搏信号处理系统界面
图2.4运行后的脉搏原信号与滤波信号回放
图2.6运行后的脉率曲线回放与判断脉率正常、输出脉搏信息
图2.7运行后的频谱分析
2.3 软件总体设计思路
3 信号回放模块设计
3.1 脉搏波信号滤波
脉搏信号是强干扰下的微弱信号，脉搏信号幅度很小，大约是微伏到毫伏的数量级范围，其主要频带范围在0～20Hz之间。

一般情况下为1Hz左右，脉搏信号可看成一个准直流信号，也可看成是一个甚低频交变信号。

根据脉搏功率谱能量分析，健康人脉搏能量绝大多数分布于0.5～5Hz，而病人脉搏在1Hz以下和较高频段(如5Hz以上或10Hz以上)仍有相当一部分的能量分布。

脉搏信号极容易引入干扰，一般情况下，由体表检测到的脉搏波主要有三种干扰：基线漂移，高频随机干扰和运动伪差噪声，还有一些干扰，比如工频干扰，等，这些干扰较之前三种干扰比较小，而且很容易去掉]。

1).基线漂移：这种噪声是由于被测对象的呼吸运动和身体移位而产生的，呼吸的频率分量通常在0.8Hz以下，每个人的呼吸频率也不同；身体移位也可以用低频分量来表征，因此基线漂移基本上都是趋势分量和低频分量。

2).高频随机干扰：脉搏信号的采集过程中有好多随机噪声和环境干扰的影响，这些噪声都表现为高频噪声，在5Hz以上。

3).运动伪差噪声：由于对病人进行血氧检测时，病人手指(或脚趾)经常会发生运动，使手指(或脚趾)与传感器之间的距离发生位移，从而得到导致测量的病人脉搏波形很不稳定的一种噪声，这种噪声和信号在同一个频带范围内，是脉搏信号去噪的一个重点和难点。

所以，为了对信号做准确的分析，在分析处理之前必须做一些必要的预处理。

针对信号中存在噪声的特点，基线漂移和呼吸等低频干扰在1Hz 以下，而脉搏信号主要在低频范围，所以可以设计让信号先通过一个巴特沃斯带通滤波器，借以滤除基线漂移、呼吸引起的干扰（考虑到不丢失太多的其他信息，通带截止频率设置为Wp=[0.9,50]，阻带截止频率设置为Ws=[0.3,140],通带波纹系数Rp=3,阻带波纹系数Rs=10），如图3.2，然后再通过一个切比雪夫二型滤波器滤除固定的工频干扰（通带截止频率设置为Wp=[30,65],阻带截止频率Ws=[45,55],通带波纹系数Rp=3,阻带波纹系数Rs=60），如图3.3。

由于Butterworth滤波器通带内有最大的平滑特性，信号经过后衰减小，因此我们选用Butterworth带通滤波器滤除基线漂移和呼吸等引起的干扰，但由于IIR 滤波器本身固有的缺点，信号通过Butterworth带通滤波器后相位会失真，故我们可设计零相位Butterworth带通滤波器去噪。

50Hz固定工频在频域中是一个点，因而要求设计的带阻滤波器有好的截止特性，而切比雪夫II型滤波器有较好的截止
特性，并且在其通带内单调，故而设计零相位切比雪夫II型滤波器滤除工频干扰。

流程图如下：
图3.1 脉搏信号滤波流程图
图3.2 巴特沃斯带通滤波器图3.3 带阻切比雪夫II型滤波器3.2 脉搏信号回放按钮设计思路
信号回放模块设计如图3.4，考虑到用户使用的方便性，信号回放过程中才可对信号进行放大与缩小操作，或者快进退操作，在回放过程中，退出系统按钮处于不激活状态，关闭窗口和继续回放也处于不激活状态，只有当信号暂停时，才
使前面三个按钮激活，或者进行下一步，输出脉率曲线功能。

图 3.4 信号回放模块设计
1).信号回放：主要采用了Matlab中的while、if、plot、axis、global、set、pause 命令。

思路先是将信号和滤波信号用plot命令分别在axes1、和axes2上画出图形，接着用axis移动坐标轴，用while进行死循环，让信号不停的回放。

2).信号暂停：主要利用if-break语句来实现，在信号回放的控件Callback函数中，定义一个全局变量flg1,初始flg1=1，接着在while每次循环中，检测flg1的值，为1继续循环回放，为0的话跳出循环，停止回放。

所以，要信号暂停回放，只需要在信号暂停按钮的Callback函数中，定义一个全局变量flg1，并命flg1=0即可。

3).继续回放：也是利用全局变量来实现，当信号暂停回放时，记录下些时的X 轴和Y轴的长度大小，点击继续回放就还是用plot命令分别在axes1、和axes2上画出图形，接着用axis移动坐标轴，用while进行死循环，让信号不停的回放。

只不过这里要注意的是，此时坐标轴的移动不是从0开始，而是从暂停时的记录下来的全局变量X轴的长度大小xmin、xmax。

4).关闭窗口：直接用delete(allchild(handles.axes1)); delete(allchild(handles.axes2)); 语句
图 3.5原信号和滤波后的脉搏信号回放
4信号放大与缩小模块设计
4.1 设计界面
考虑到信号回放过程中每个人的脉搏的曲线都会有所不一样，两个波峰之间的时间，波峰的数值大小，为了能方便看到每个点的数据，我在信号放大与缩小模块中加了四个按钮，可沿X轴缩小放大和沿Y轴缩小放大，如图4.1。

另外，考虑到操作的方便性，设计了在信号回放或者继续回放的过程中，这四个按钮才被
激活，信号暂停时，这四个按钮都处于不激活的状态。

图4.1 信号放大与缩小模块
4.2 信号放大与缩小按钮设计思路
1).X轴缩小: 即将信号沿X轴压缩，如图4.3，X轴的长度变长，但是Y轴保持不变，这里，主要也是利用axis命令，实时的改变X轴长度，并且用while进行死
循环，直到遇到flg1=0，才跳出循环。

程序流程图如下：
图4.2 信号放大流程图
图4.3 X轴缩小
2).X轴放大: 即将信号沿X轴拉伸，如图4.4，X轴的长度变短，但是Y轴保持不变，这里，主要也是利用axis命令，实时的改变X轴长度，并且用while进行死循环，直到遇到flg1=0，才跳出循环。

图4.4 X轴放大
3).Y轴缩小: 即将信号沿Y轴压缩，如图4.5，Y轴的长度变长，但是Y轴保持不变，这里，主要也是利用axis命令，实时的改变X轴长度，并且用while进行死循环，直到遇到flg1=0，才跳出循环。

图4.5Y轴缩小
4).Y轴放大: 即将信号沿Y轴拉伸，如图4.6，Y轴的长度变短，但是Y轴保持
不变，这里，主要也是利用axis命令，实时的改变X轴长度，并且用while进行死循环，直到遇到flg1=0，才跳出循环。

图4.6 Y轴放大
5 信号快进退模块设计
5.1 设计界面
考虑到有时信号的采样时间很长，为了能快速浏览信号的波形或者看之前的波形，我特别又设计此模块。

图 5.1 信号快进退模块界面
5.2信号快进退按钮设计思路
1).快进：设计全局变量pauset=pauset*0.5,即让每次移动坐标轴的时间间隔减小，再令Xmin=Xmin+0.1,Xmax=Xmax+0.1，移动X轴，这样回放的速度就会变快，按钮每按一次，回放速度加快一倍。

2).慢进：设计全局变量pauset=pauset*2,即让每次移动坐标轴的时间间隔加大，再令Xmin=Xmin+0.1,Xmax=Xmax+0.1，移动X轴，这样回放的速度就会变慢，按钮每按一次，回放速度减少一倍。

3).快退：设计全局变量pauset=pauset*0.5,即让每次移动坐标轴的时间间隔减小，再令Xmin=Xmin-0.1,Xmax=Xmax-0.1，移动X轴，这样回放的速度就会变快，按钮
每按一次，回放速度加快一倍。

4).慢退：设计全局变量pauset=pauset*2，即让每次移动坐标轴的时间间隔减小，再令Xmin=Xmin-0.1,Xmax=Xmax-0.1，移动X 轴，这样回放的速度就会变快，按钮每按一次，回放速度加快一倍。

6 脉率实时处理模块设计
6.1 脉搏信号峰值提取
脉搏信号的特征提取是脉搏信号分析中的关键，能不能准确的提取出来脉搏信号的特征关系到能不能准确的分析脉搏信号，其时域中提取的主要特征有脉搏的峰值、周期等。

而在这其中最重要的是峰值的提取，准确的提取出来峰值，可为进一步分析其周期，脉率等其他 参数打下基础，如图6.1 显示的脉搏波，h 为脉搏波形的高度。

图6.1 脉搏波形的高度
提取信号的峰值基本思想是基于阈值的方法，统计分析采集到得脉搏信号发现，虽然每个脉搏波的波峰值大小并不相同，但他们总在一个范围内波动，波动的范围基本上不超过最大波形高度的0.3倍，为了更为可靠的检测波峰，以0.4倍为参考。

因此可以认为波峰点是在每个脉搏周期中波形的最大值附近，大于其邻域内所有点的点。

由于所用数据采样频率为800HZ ，通过对大量脉搏波测量数据分析，把邻域长度定义为400个数据点可以准确识别波峰点。

假设采集到的脉搏波为数字信号序列][j P , ],0[N j ∈实际波峰][j P 的判别条件为:
⎪⎩⎪⎨⎧+-∈≥∆≤-]200,200[],[][]}[][max{p p p p j j j j P j P j P j P
公式中4.0]})[min{]}[(max {⨯-=∆j P j P ， p
j 为波峰点的位置。

提取出波峰点的位置后，其他时域的特征，如周期、脉率等可相应而出。

峰值提取程序流程如下:
6.2 脉率实时处理界面设计
界面设计如图6.2，脉率信号输出窗口可输出平均脉率、采样频率、采样时长、最大脉率、最小脉率等。

图6.2 脉率实时处理模块
6.3 脉率实时处理按钮设计思路
1).输出脉率曲线：实时动态输出脉率变化趋势曲线与脉率峰值变化曲线，如图6.3同时在文本框1中实时平均脉率，这也是实时更新输出。

图6.3 软件实时输出脉率趋势曲线
2).暂停回放：停止脉率曲线的更新和平均脉率的更新，同时激活输出脉搏信息按钮
3).输出脉搏信息：输出所有脉率点的平均脉率，同时在列表中输出采样频率、采样时长、最大脉率、最小脉率，并跳出判断脉率是否正常的信息框（如图6.4）。

之后激活频谱分析按钮。

图6.4 输出脉搏信号
图6.5 输出脉搏信息到列表和文本框
4).频谱分析：输出滤波后信号的幅频图和信号功率谱图
图6.7 频谱分析
5).清除波形：清除窗口的绘图
6).输出文档：输出TXT文档，以电脑当前的日期和时间命名，包含有所有计算出的脉率点，以及采样频率、采样时长、最大脉率、最小脉率。

7 参考文献
[1] 罗志昌,张松,杨益民,脉搏波的工程分析与临床应用[M].北京：科学出版社，2006：1-20
[2] 郑阿奇,MATLAB实用教程[M].北京：电子工业出版社，2004．
[3] 王炳和,罗建,相敬林等. 人体脉搏功率谱分析与中医脉诊机理研究[J]. 西北大学学报,2001,31( 1): 21-25.
8 致谢
这里，首先要感谢候文生老师提供让我设计此软件的机会。

软件完成后，感触颇多，觉得又学到了点了知识，很有成就感。

软件前前后后花了总共一个礼拜的时间，这星期里起床后吃饭到睡觉时都在想着怎么来编写代码，怎么来完善软件功能，一整天开机就Matlab，脑子就直想着怎么写程序。

设计每个控件按钮时，总会碰到很多问题,因为调试时各种Bug,光是界面就重画了10来次。

但我不放弃不后退，每次通过网上搜集资料、反复查看代码、和同学之间进行交流总能把问题解决好，而且又想到新的功能设计。

虽说软件还不是很强大，设计也还不是很好，但还是算是自己设计的第一个软件。

最后，衷心感谢候老师安排的此次实验，希
望以后还能有更多这样的实践机会。

附录:
（程序太多了，只附上主要部分的程序）
1.打开数据
function pushbutton5_Callback(hObject, eventdata, handles)
kong='';
set(handles.listbox1,'string',kong);
set(handles.edit1,'string',kong);
axes(handles.axes1);
grid on;
axis([0,1,-2,4])
axes(handles.axes2);
grid on;
axis([0,1,-2,4])
global str;
[filename,pathname]=...
uigetfile({'*.txt'},'choose');
str=[pathname filename];
set(handles.pushbutton1,'enable','on');
2.退出系统
function pushbutton6_Callback(hObject, eventdata, handles)
button1=questdlg('你确定退出吗？','退出程序','确定','取消','确定'); if strcmp(button1,'确定')
close(gcf);
end
3.原信号与滤波信号回放
function pushbutton1_Callback(hObject, eventdata, handles)
set(handles.pushbutton2,'enable','on');
set(handles.pushbutton5,'enable','off');
set(handles.pushbutton6,'enable','off');
set(handles.pushbutton10,'enable','off');
set(handles.pushbutton11,'enable','off');
set(handles.pushbutton15,'enable','on');
set(handles.pushbutton16,'enable','on');
set(handles.pushbutton17,'enable','on');
set(handles.pushbutton18,'enable','on');
set(handles.pushbutton19,'enable','on');
set(handles.pushbutton20,'enable','on');
set(handles.pushbutton21,'enable','on');
set(handles.pushbutton22,'enable','on');
global str
global t
global maibo
global ymin
global ymax
global xmin
global xmax
global pauset
pauset=0.5;
data=load(str);
t=data(:,1);maibo=data(:,2);
ymin=min(data(:,2))-0.5;
ymax=max(data(:,2))+0.5;
%%
global s21
global dt
yy=data;
x=yy(:,2)';
dt=yy(:,1)';
fs=800;
%%去除基线漂移
k = .7; % cut-off value
fc=0.3/fs;
alpha = (1-k*cos(2*pi*fc)-sqrt(2*k*(1-cos(2*pi*fc))-k^2*sin(2*pi*fc)^2))/(1-k); y = zeros(size(x));
for i = 1:size(x,1)
y(1,:) = filtfilt(1-alpha,[1 -alpha],x(1,:));
end
x1=x-y;
%%butter带通滤波
Wp=[0.9 50]/400;Ws=[0.3 140]/400; [n,Wn] = buttord(Wp,Ws,3,10); [b,a] = butter(n,Wn);
%%cheby去除50hz工频干扰
Wp1=[30 65]/400;Ws1=[45 55]/400; rp=3;rs=60;
[n1,Wn1] = cheb2ord(Wp1,Ws1,rp,rs); [b1,a1] = cheby2(n1,rs,Wn1,'stop'); %%信号之间的对比
x2=filtfilt(b,a,x1);
s21=filtfilt(b1,a1,x2);
axes(handles.axes1);
plot(t,maibo);
title('原始脉搏信号');
xlabel('时间/s');
ylabel('幅值/mv');
grid
set(gca,'position',[10 27.5 120 20]); xmin=0;
xmax=max(t)/10;
axis([xmin,xmax,ymin,ymax])
axes(handles.axes2);
plot(dt,s21);grid;
title('滤波去噪后脉搏信号');
xlabel('时间/s');
ylabel('幅值/mv');
set(gca,'position',[10 3 120 20]); axis([xmin,xmax,ymin,ymax])
%%%
global flg1
基于MATLAB的脉搏信号处理软体系统
flg1=1;
while 1
if (flg1==0)
break;
end
xmin=xmin+0.1;
xmax=xmax+0.1;
axes(handles.axes1);
plot(t,maibo);
title('原始脉搏信号');
xlabel('时间/s');
ylabel('幅值/mv');
grid
set(gca,'position',[10 27.5 120 20]);
axis([xmin,xmax,ymin,ymax]); %移动坐标系
axes(handles.axes2);
plot(dt,s21);grid;
title('滤波去噪后脉搏信号');
xlabel('时间/s');
ylabel('幅值/mv');
set(gca,'position',[10 3 120 20]);
axis([xmin,xmax,ymin,ymax])
pause(pauset);
end
guidata(hObject,handles);
4.x轴缩小
function pushbutton17_Callback(hObject, eventdata, handles) global xmin
global xmax
global s21
global dt
global t
global maibo
global ymin
global ymax
global pauset
xmin=1/2*xmin;
global flg1
flg1=1;
while 1
if (flg1==0)
break;
end
xmin=xmin+0.1;
xmax=xmax+0.1;
axes(handles.axes1);
plot(t,maibo);
title('原始脉搏信号');
xlabel('时间/s');
ylabel('幅值/mv');
grid
set(gca,'position',[10 27.5 120 20]);
axis([xmin,xmax,ymin,ymax])
axes(handles.axes2);
plot(dt,s21);grid;
title('滤波去噪后脉搏信号');
xlabel('时间/s');
ylabel('幅值/mv');
set(gca,'position',[10 3 120 20]);
axis([xmin,xmax,ymin,ymax])
pause(pauset);
end
guidata(hObject,handles);
5.快进
function pushbutton15_Callback(hObject, eventdata, handles) global pauset
pauset=pauset*0.5;
6.输出脉率曲线
function pushbutton7_Callback(hObject, eventdata, handles) set(handles.pushbutton1,'enable','off');
set(handles.pushbutton5,'enable','off');
set(handles.pushbutton6,'enable','off');
set(handles.pushbutton8,'enable','off');
set(handles.pushbutton9,'enable','off');
set(handles.pushbutton2,'enable','off');
set(handles.pushbutton10,'enable','off');
set(handles.pushbutton11,'enable','off');
set(handles.pushbutton17,'enable','off');
set(handles.pushbutton18,'enable','off');
set(handles.pushbutton19,'enable','off');
set(handles.pushbutton20,'enable','off');
set(handles.pushbutton15,'enable','off');
set(handles.pushbutton16,'enable','off');
set(handles.pushbutton21,'enable','off');
set(handles.pushbutton22,'enable','off');
set(handles.pushbutton23,'enable','off');
set(handles.pushbutton14,'enable','on');
global str
global pauset
global t
global TT
yy=load(str);
%HH=1./yy(3:1)
x=yy(:,2)';
fs=1000;
%%去除基线漂移
k = .7; % cut-off value
fc=0.3/fs;
alpha = (1-k*cos(2*pi*fc)-sqrt(2*k*(1-cos(2*pi*fc))-k^2*sin(2*pi*fc)^2))/(1-k); y = zeros(size(x));
for i = 1:size(x,1)
y(1,:) = filtfilt(1-alpha,[1 -alpha],x(1,:));
end
x1=x-y;
%%butter带通滤波
Wp=[0.9 50]/500;Ws=[0.3 140]/500;
[n,Wn] = buttord(Wp,Ws,3,10);
[b,a] = butter(n,Wn);
%%cheby去除50hz工频干扰
Wp1=[30 65]/500;Ws1=[45 55]/500;
rp=3;rs=60;
[n1,Wn1] = cheb2ord(Wp1,Ws1,rp,rs);
[b1,a1] = cheby2(n1,rs,Wn1,'stop');
%%信号之间的对比
x2=filtfilt(b,a,x1);
s21=filtfilt(b1,a1,x2);
%%频谱图
F=fft(s21);
%%R波的提取
A=s21;
PM=max(A);
MM=min(A);
G=(PM-MM)*0.3; %峰值大小波动范围不超过最大波形高度的0.3倍
if max(A(1:200))==max(A(1:550))
[P1(1),t1(1)]=max(A(1:200));
cnt1=1;
elseif max(A(1:200))~=max(A(1:550))
cnt1=0;
end
for i=201:length(A)-200 %因为周期大于600，波峰大于其左右各200个点内的所有值
if PM-A(i)<G&&A(i)==max(A(i-200:i+200))
P1(cnt1+1)=A(i);
t1(cnt1+1)=i*1.25; %%%时间！！！
cnt1=cnt1+1;
end
end
if max(A(length(A)-199:length(A)))==max(A(length(A)-550:length(A)))
[P1(cnt1+1),t1(cnt1+1)]=max(A(length(A)-199:length(A)));
cnt1=cnt1+1;
end
% %求周期T
j=1;
for m=2:length(t1)
T(j)=t1(m)-t1(m-1);
j=j+1;
MMT(j-1)=(1./mean(T))*1000*60;
end
TT=(1./T)*1000*60;
TT=(1./T)*1000*60;
yymin=1;
yymax=2;
mmax=max(P1)+2;
axes(handles.axes1);
plot(P1);
axis([yymin yymax 0 mmax]);
grid on;
xlabel('峰值数');ylabel('峰值');
title('脉搏峰值变化曲线');
set(gca,'position',[10 27.5 120 20]);
axes(handles.axes2);
plot(TT);
axis([yymin yymax 0 100]);
grid on;
xlabel('脉搏数');ylabel('脉率（次/min）'); title('脉率变化曲线');
set(gca,'position',[10 3 120 20]);
mT=mean(TT);
global maibo
maibo=mT;
global flg3
flg3=1;
global xxmin
global xxmax
i=1;
while 1
if (flg3==0)
xxmin=yymin;
xxmax=yymax;
break;
end
if (yymax==cnt1);
yymin=0;
yymax=5;
end
mai=round(MMT(i));
mai=num2str(mai);
i=i+1;
if (i==cnt1)
i=1;
end
set(handles.edit1,'string',mai,'fontsize',18);
yymin=yymin+1;
yymax=yymax+1;
axes(handles.axes1);
plot(P1);
axis([yymin yymax 0 mmax]);
grid on;
xlabel('峰值数');ylabel('峰值');
title('脉搏峰值变化曲线');
set(gca,'position',[10 27.5 120 20]);
axes(handles.axes2);
plot(TT);
axis([yymin yymax 0 100]);
grid on;
xlabel('脉搏数');ylabel('脉率（次/min）');
title('脉率变化曲线');
set(gca,'position',[10 3 120 20]);
pause(pauset);
end
7.输出脉搏信息
function pushbutton8_Callback(hObject, eventdata, handles)
global maibo
global TT
global t
str1='脉搏采样频率：';
HH1=1./t(2);
str2=num2str(HH1);
str3='Hz';
strr=[str1 str2 str3];
ktr1='脉搏采样时长：';
HH2=max(t);
ktr2=num2str(HH2);
ktr3='S';
ktrr=[ktr1 ktr2 ktr3];
htr1='脉率最大值:';
htr2=num2str(round(max(TT)));
htr3='次/min';
htrr=[htr1 htr2 htr3];
gtr1='脉率最小值:';
gtr2=num2str(round(min(TT)));
gtr3='次/min';
gtrr=[gtr1 gtr2 gtr3];
trrr=strvcat(strr,ktrr,htrr,gtrr);
set(handles.listbox1,'string',trrr);
maibo1=num2str(round(maibo));
if maibo<60
msgbox('健康成人在安静状态下脉率为60～100次/min,您的脉率过慢，建议到
医院检查','脉率正常与否判断');
end
if maibo>100
msgbox('健康成人在安静状态下脉率为60～100次/min,您的脉率过快，建议到医院检查','脉率正常与否判断');
else
msgbox('健康成人在安静状态下脉率为60～100次/min,您的脉率属于正常范围，请继续保持','脉率正常与否判断');
end
set(handles.pushbutton24,'enable','on');
set(handles.pushbutton8,'enable','off');
set(handles.edit1,'string',maibo1,'fontsize',18);
7.频谱分析
function pushbutton24_Callback(hObject, eventdata, handles)
set(handles.pushbutton9,'enable','on');
set(handles.pushbutton23,'enable','on');
global str
yy=load(str);
x=yy(:,2)';
fs=1000;
%%去除基线漂移
k = .7; % cut-off value
fc=0.3/fs;
alpha = (1-k*cos(2*pi*fc)-sqrt(2*k*(1-cos(2*pi*fc))-k^2*sin(2*pi*fc)^2))/(1-k);
y = zeros(size(x));
for i = 1:size(x,1)
y(1,:) = filtfilt(1-alpha,[1 -alpha],x(1,:));
end
x1=x-y;
%%butter带通滤波
Wp=[0.9 50]/500;Ws=[0.3 140]/500;
[n,Wn] = buttord(Wp,Ws,3,10);
[b,a] = butter(n,Wn);
%%cheby去除50hz工频干扰
Wp1=[30 65]/500;Ws1=[45 55]/500;
rp=3;rs=60;
[n1,Wn1] = cheb2ord(Wp1,Ws1,rp,rs);
[b1,a1] = cheby2(n1,rs,Wn1,'stop');
%信号对比
x2=filtfilt(b,a,x1);
s21=filtfilt(b1,a1,x2);
%%R波的提取
A=s21;
PM=max(A);
MM=min(A);
G=(PM-MM)*0.3; %峰值大小波动范围不超过最大波形高度的0.3倍
if max(A(1:200))==max(A(1:550))
[P1(1),t1(1)]=max(A(1:200));
cnt1=1;
elseif max(A(1:200))~=max(A(1:550))
cnt1=0;
end
for i=201:length(A)-200 %因为周期大于600，波峰大于其左右各200个点内的所有值
if PM-A(i)<G&&A(i)==max(A(i-200:i+200))
P1(cnt1+1)=A(i);
t1(cnt1+1)=i*1.25;
cnt1=cnt1+1;
end
end
if max(A(length(A)-199:length(A)))==max(A(length(A)-550:length(A)))
[P1(cnt1+1),t1(cnt1+1)]=max(A(length(A)-199:length(A)));
cnt1=cnt1+1;
end
j=1;
for m=2:length(t1)
T(j)=t1(m)-t1(m-1);
j=j+1;
end
%%频谱图
F=fft(s21);
axes(handles.axes1);
fumax=max(abs(F)*2/length(s21))+0.3;
plot([0:length(s21)-1]/(length(s21)*(1/fs)),abs(F)*2/length(s21));grid;
xlabel('频率/HZ');ylabel('幅值');
title('幅频图');
axis([0 100 0 fumax]);
set(gca,'position',[10 27.5 120 20]);
%信号的功率谱估计
axes(handles.axes2);
w=hanning(128);
nfft=1024;
noverlap=64;
[p,f]=spectrum(x,nfft,noverlap,w,fs);
plot(f,p);grid;
xlabel('频率(Hz)');
ylabel('功率谱(dB)');
title('信号的功率谱图');
axis([0 50 0 7]);
set(gca,'position',[10 3 120 20]);
8.输出文档
function pushbutton9_Callback(hObject, eventdata, handles) global TT
str=datestr(now);
str=strrep(str,':','_');
str=strrep(str,' ','_');
str=strrep(str,'-','_');
filename=strcat(str,'.txt');
fid=fopen(filename,'wt');%写入文件路径
TTT=TT';
[m,n]=size(TTT);
for i=1:1:m
for j=1:1:n
if j==n
fprintf(fid,'%g\n',TTT(i,j));
else
fprintf(fid,'%g\t',TTT(i,j));
end
end
end
fclose(fid);。