技术领域
本发明属于医学数字信号处理领域，涉及一种对脉象信号的去噪、基线漂移方法矫正以及利用曲线拟合平滑信号的方法。
背景技术
脉象诊断是中国传统医学中一种十分重要的诊断方法，历史悠久，内容丰富。两千多年来，在临床上得到了广泛的应用，积累了大量的理论的经验。然而，传统中医脉诊需要医生经过长时间的专业训练，诊断结果因人主观判断而异，从而限制了脉诊的发展。近三十年来，脉诊客观化研究取得了许多有益的成果，为脉象信号分析提供了基础。
在脉搏信号的采集过程中，因人体呼吸、肌肉轻微运动、机电噪声等因素客观存在，采集到的脉搏信号中存在基线漂移和大量的噪声。脉搏信号的预处理作为脉诊客观化中重要的一环，就是要对原始脉搏信号进行降噪、矫正、分割、平滑等操作，以有益于进一步的脉象识别、辨识等工作。
小波变换是信号的时间—尺度分析方法，在时、频域具有表征信号局部特征的能力。在脉象信号预处理研究当中，大量的实验表明[1,2]，利用小波变换来对脉象信号去噪、检测信号基线漂移严重程度，矫正脉象信号的基线漂移是十分有效的。
脉搏波的分解分析方法最初由钱伟立等提出[3]，该方法指出一个周期的脉搏波可以分解成三个子波，分别对应一个脉象周期的主峰波、重搏前波、重搏波。脉搏波的分解方法可以看成一个曲线拟合问题，将一个周期的脉搏波分解成一组基函数的线性组合，最常用的基函数是高斯函数。高斯拟合有时也作为一种预处理方法，对一个周期的脉象信号进行平滑处理[4,5]。对于高斯曲线拟合的对应的优化问题，当采用平方误差作为目标函数时，该优化问题是非凸的，在采用基于梯度的方法(如：非线性二乘法。)进行问题求解时，得到的结果就十分依赖于初始解的选择的好坏。在已有的研究当中，都是通过反复实验、对比，来设定迭代的初始点。但是，通过反复实验来设定初始迭代点，不仅费时费力，而且该初始解对于新的脉象信号并不鲁棒。
参考文献：
[1]Xu L,Zhang D,Wang K,et al.《Baseline wander correction in pulse waveforms using wavelet-based cascaded adaptive fiter[J]》.Computers in Biology and Medicine,2007,37(5):716-731；
[2]岳沛平,李训铭.《基于小波分析和BP神经网络识别的中医脉象信号识别系统[J]》.计算机与现代化,2006(12):1-4；
[3]钱伟立,徐兰义.《高斯函数分解法提取脉搏波特征[J]》.中国生物医学工程学报,1994,13(1):1-7；
[4]陈雪峰,《脉搏波特征提取算法及其应用研究》,Ph.D.thesis,大连:大连理工大学,2006；
[5]L Wang,L Xu,S Feng,M Meng et al.,《Multi-gaussian fitting for pulse waveform using weighted least squares and multi-criteria decision making method[J].Computers in Biology and Medicine》,vol.43,No.11,pp:1661-1672,2013。
发明内容
为了克服上述现有技术存在的问题，本发明提出了一种基于小波变换和曲线拟合的脉搏波预处理方法，通过脉搏波降噪、周期提取、信号平滑的步骤，实现脉搏波预处理。
本发明提出一种基于小波变换和曲线拟合的脉搏波预处理方法，该方法包括以下步骤：
第一步，用小波全局阈值法来对采集到的脉象信号进行降噪处理，具体包括：
对收集到的原始信号S进行6层db4小波分解，
计算小波全局阈值：
其中，len(S)表示信号S的长度，
对于所有的小波分解系数，若值小于小波全局阈值TPTR，置零，
执行小波逆变换，得到降噪后的脉象信号；
第二步，矫正脉象信号的基线漂移及周期提取，具体包括：
对降噪后的脉象信号执行六层Meyer小波分解，
计算降噪后的脉象信号的能量比(Energy Ratio,ER)：
其中，norm(A1)和mean(CA1)分别表示小波分解的第一层小波系数所构成向量的二范数和平均值，norm(CA6)和mean(CA6)类似；
若ER＜50dB，执行一次小波低频滤波：置小波分解的第六层近似系数(CA6)为零，执行小波逆变换，得到滤除部分基线漂移的脉象信号，
检测出脉象信号中脉象周期的起始点，
以检测到的周期起始点为输入，执行三次样条插值，
从脉象信号中减去三次样条插值的结果，得到去除基线漂移的脉象信号，
随机地选取一个周期代表该脉象，将幅值归一化到1，长度归一化到150，即一个周期含150个采样点，进行进一步处理；
第三步，基于高斯曲线拟合来平滑脉象信号，即：设一个周期的脉象信号为O(t)，高斯曲线拟合问题就等价于优化问题：
argmin C Σ t = 1 150 1 2 ( O ( t ) - Σ i = 1 n g i ( t ) ) 2 , ]]>
其中，O(t)表示要拟合的一个周期的脉象信号，C表示优化问题的解，即各个高斯函数的参数列表；n表示高斯函数的个数，高斯函数的公式为：
第四步：在采集的信号当中难免会一些由于操作失误而产生的不包含或部分包含脉象信息的脉象信号，对于这一类脉象信号，自动抛出异常，完成异常脉象信号交互处理。
所述检测出脉象信号中脉象周期的起始点的步骤，还包括以下流程：
找出脉象信号中所有的极小值点的位置，按照从小到大的顺序排列，构成有序列表Listmin；
在集合Listmin中找出具有最小脉象幅值，即脉象信号的全局最小值点Lmin的位置，并标记该全局最小值点(Lmin,S(Lmin))为一个周期的起始点；
以Lmin为基准点，在列表中，向其左、右方向顺序搜索，记当前基准点为Lcurrent，读取到的下一个极小值点为Lnext，若70≤|Lnext-Lcurrent|≤120，标记Lnext为某周期的起始点，置Lcurrent＝Lnext为新的基准点，再进行搜索，直至所有的极小值点都被访问了一遍。
针对所述第三步中的优化问题，还提出了一个自动求解对应非凸优化问题的递归算法RGM，具体的步骤如下：
输入当前参数列表C，要拟合的脉象信号O(t),t∈I＝{1,2,...,T}，已使用高斯函数个数的计数k，高斯函数的总个数n；
若k≥n，将参数列表C中的参数，按照bi从小到大的顺序排列，终止；
计算各个高斯函数的值，置
置
按照以下公式计算出ak,bk,ck，三个参数分别表示高斯函数gi(t)的幅值、均值和宽度：
ak＝maxP(t)，其中t∈SS，
bk＝t*，其中P(t*)＝ak，
其中|SS|表示集合SS的势；
置参数列表C＝C∪{ak,bk,ck}；
以参数列表C中的参数为初始迭代点，执行非线性最小二乘，并对参数列表作更新；
递归调用RGM(C,O,k+1,n)。
与现有技术相比，本发明能够有效地提取出脉象信号的周期，所提出的曲线拟合算法，在3.1GHz上的机器上运行，仅仅花费约120s完成对所有脉象数据的拟合，且97.62％的脉象的平均拟合误差小于0.35，对于误差大于0.35的拟合结果，也保持的脉象信号的形态特征，并未丢失诊断信息。
附图说明
图1为本发明的一种基于小波变换和曲线拟合的脉搏波预处理方法的整体流程示意图；
图2为采集到的包含基线漂移和大量噪声的脉象信号波形示意图；
图3为六层小波分解示意图；
图4为经过小波全局阈值去噪后的脉象信号波形示意图；
图5为经过基线漂移矫正后的脉象信号波形示意图；
图6为归一化的一个脉象周期示意图
图7为高斯拟合结果示例图；
图8为异常信号处理交互界面图。
具体实施方式
对于不同的脉象，其初始迭代点应该不同。本发明通过观察脉象数据信号和高斯函数的特性，以及考察一个周期的脉象信号的形成过程，提供了一个贪心递归算法来自动求解高斯拟合问题。在我们收集的脉象数据集上的实验表明，提出的算法能够有效地求解脉搏波的高斯拟合问题。
此外，对于数据集中部分少量或不包含人体脉象信息的信号，本发明方法在提取周期的过程中提供给用户一个简单、快速、有效的交互，以让专业人员做出进一步的判断：丢弃完全不包含脉象信息的数据还是截取采集到脉象信号的一个片段以继续处理。
本算法首先基于小波变换来对信号降噪、矫正基线漂移，再对提取出的一个周期的脉象信号执行高斯曲线拟合平滑信号。
具体技术方案包括以下内容：(如图1所示)
第一步：小波全局阈值降噪：用小波全局阈值法来对采集到的脉象信号进行降噪处理，分为以下步骤：
(1)对收集到的原始信号S进行6层db4小波分解(原始信号波形如图2所示)。
(2)计算小波全局阈值：
T P T R = 2 l o g ( l e n ( S ) ) , ]]>
其中，len(S)表示信号S的长度。
(3)对于所有的小波分解系数，若值小于小波全局阈值TPTR，置零。
(4)执行小波逆变换，得到降噪后的脉象信号(如图3所示)。
第二步：基线漂移矫正及周期提取。对于降噪后的脉象信号，由于人体呼吸等因素，存在着基移，利用一种基于小波变换的方法来矫正脉象信号的基线漂移，分为以下步骤：
(1)对降噪后的脉象信号执行六层Meyer小波分解；
Meyer的小波函数是在频率域中定义的，为：
Ψ ( ω ) = 0 , | ω ∉ [ 2 π 3 , 8 π 3 ] | ( 2 π ) 1 2 e j ω 2 s i n ( π 2 v ( 3 2 π | ω | - 1 ) ) , 2 π 3 ≤ ω ≤ 4 π 3 ( 2 π ) 1 2 e j ω 2 sin ( π 2 v ( 3 2 π | ω | - 1 ) ) , 4 π 3 ≤ ω ≤ 8 π 3 ]]>
其中，v(x)为构造Meyer小波的辅助函数：
v(x)＝x4(35-84x+70x2-20x3),x∈[0,1].
小波分解的层数设置为6，得到第六层的近似信号CA6和每一层的细节信号CD1，...,CD6。
(2)计算信号的能量比(Energy Ratio,ER)：
E R = 20 log 10 ( n o r m ( C A 1 - m e a n ( C A 1 ) ) n o r m ( C A 6 - m e a n ( C A 6 ) ) ) , ]]>
其中，norm(A1)和mean(CA1)分别表示小波分解的第一层小波系数所构成向量的二范数和平均值，norm(CA6)和mean(CA6)类似。
(3)对于此脉象信号，若ER＜50dB，执行一次小波低频滤波：置小波分解的第六层近似系数(CA6)为零，执行小波逆变换，得到滤除部分基线漂移的脉象信号；
(4)检测出脉象信号中脉象周期的起始点：
(4-1)、找出脉象信号中所有的极小值点的位置，按照从小到大的顺序排列，构成有序列表Listmin；
(4-2)、在集合Listmin中找出具有最小脉象幅值，即脉象信号的全局最小值点Lmin的位置，并标记该全局最小值点(Lmin,S(Lmin))为一个周期的起始点；
(4-3)、以Lmin为基准点，在列表中，向其左、右方向顺序搜索。记当前基准点为Lcurrent，读取到的下一个极小值点为Lnext，若标记Lnext为某周期的起始点，置Lcurrent＝Lnext为新的基准点，再进行搜索直至所有的极小值点都被访问了一遍；
(5)以检测到的周期起始点为输入，执行三次样条插值；
(6)从脉象信号中减去三次样条插值的结果，得到去除基线漂移的脉象信号(如图5所示)；
(7)随机地选取出一个周期代表该脉象，将幅值归一化到1，长度归一化到150,进行进一步处理(如图6所示)；
第三步：基于高斯曲线拟合的脉象信号平滑：对于提取出的周期，尽管已经做了降噪处理，但在信号难以避免地还包含着少量噪声点；因此，采用高斯拟合方法来平滑脉象信号，并提出了一个自动求解对应非凸优化问题的递归算法。即：设一个周期的脉象信号为O(t)，高斯曲线拟合问题就等价于优化问题：
argmin C Σ t = 1 150 1 2 ( O ( t ) - Σ i = 1 n g i ( t ) ) 2 , ]]>
其中，O(t)表示要拟合的一个周期的脉象信号。C表示优化问题的解，即各个高斯函数的参数列表；n表示高斯函数的个数。高斯函数的形式如下：
g i ( t ) = a i exp ( - ( t - b i c i ) 2 ) . ]]>
递归算法RGM具体的步骤如下：
(1)输入：当前参数列表C，要拟合的脉象信号O(t),t∈I＝{1,2,...,T}，已使用高斯函数个数的计数k，高斯函数的总个数n；
(2)若k≥n，将参数列表C中的参数，按照bi从小到大的顺序排列，终止；
(3)计算各个高斯函数的值，置
(4)置
(5)按照以下公式计算出ak,bk,ck：
ak＝maxP(t),其中t∈SS,
bk＝t*,其中P(t*)＝ak,
c k = Σ t = 1 150 t - b k l o g a k P ( k ) | S S | , ]]>
其中|SS|表示集合SS的势。
(6)置参数列表C＝C∪{ak,∪{ak，bk,ck}；
(7)以参数列表C中的参数为初始迭代点，执行非线性最小二乘，并对参数列表作更新；
(8)递归调用RGM(C,O,k+1,n)。
本发明中，参数列表为空，调用RGM(C,O,1,3)对一个脉象进行高斯拟合。
第四步：异常脉象交互处理。在采集的信号当中难免会一些由于操作失误而产生的不包含或部分包含脉象信息的脉象信号。对于这一类数据，我们的程序自动抛出异常。
下面结合附图对本发明实施例进一步详细说明。本发明的方法在所采集的共包含671个脉象信号的数据集上运行良好。
本方法由Matlab编写实现，并提供了一个用户友好的图形界面。详细实施步骤见下文描述：
步骤1、导入格式为“.csv”的脉象信号数据文件；
步骤2、从数据文件中读取一个脉象信号数据；
步骤3、对脉象信号进行降噪处理：
3-1、对信号执行六层db4小波分解；
3-2、计算全局阈值，置小于此阈值的小波系数为零；
3-3、执行小波逆变换，得到降噪后的脉象信号。
步骤4、基线漂移矫正：
4-1、执行六层Meyer小波变换；
4-2、计算脉象信号的能量比ER；
4-3、若ER小于50dB，则认为该信号的基线漂移较为严重，先执行一步小波低频滤波；找到各个周期的起始点，并以这些周期起始点为输入进行三次样条插值；
4-4、从脉象信号中减去三次样条插值的结果。
步骤5、周期提取与高斯拟合：
5-1、从脉象信号的所有周期中随机地选取一个周期代表此脉象以作进一步处理；
5-2、执行递归的高斯拟合算法。
步骤6、异常信号处理：
6-1、若在处理过程中，存在异常信号，程序抛出异常，显示最原始的脉象信号，并提示用户做进一步地操作(图8)；
6-2、用户观察脉象信号，选择截取信号的一个片段作处理还是丢弃该不包含、任何人体脉象信息的数据；
6-3、若用户选择丢弃该信号，此脉象信号处理终止。读取下一个脉象信号数据；
6-4、若用户选择截取一段信号作处理，程序根据用户鼠标点击的将要选择所要截取片段的起始位置，截取信号片段，再次尝试周期提取；若仍抛出异常，返回6-1。
否则，预处理成功结束；读取下一个脉象信号数据。
以上所述，仅为本发明的较佳实施例而已，用于帮助理解本发明的方法及核心思想，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，所以本说明书内容不应理解为对本发明的限制。