CN103932687A - 一种脉象信号预处理方法和装置 - Google Patents

Abstract

本发明公开了一种脉象信号预处理方法和装置，包括以下步骤：对采集的脉象信号通过小波变换进行信号分解；对分解后的脉象信号进行去除噪声；对去除噪声后的信号消除基线漂移。从而，本发明能够有效提高脉象信号时域特征提取的精度和准确度。

Description

一种脉象信号预处理方法和装置

技术领域

本发明涉及中医脉象测量技术领域，特别是指一种脉象信号预处理方法和装置。

背景技术

在生物医学信号检测技术中，脉象信号是mV级信号，主要的频率范围<40Hz，一般在10Hz以下。由于使用环境的复杂和被测量者活动，采集的脉象信号将会受到多种干扰，具有较强的随机性和背景噪声，而且属于非线性、非平稳的微弱信号。这些干扰主要是交流电引起的工频干扰、肌电干扰、人体的微动与电极接触不良引起的电极接触噪声、运动伪迹(基线变化)和由于呼吸引起的基线漂移，因此在对脉象信号的时域特征进行提取之前消除脉象信号中多种干扰是极为重要，并且能够影响到脉象信号时域特征提取的精度和准确度。

发明内容

有鉴于此，本发明的目的在于提出一种脉象信号预处理方法和装置，能够有效提高脉象信号时域特征提取的精度和准确度。

基于上述目的本发明提供的脉象信号预处理方法，包括以下步骤：

对采集的脉象信号通过小波变换进行信号分解；

对分解后的脉象信号进行去除噪声；

对去除噪声后的信号消除基线漂移。

可选地，所述对采集的脉象信号通过小波变换进行信号分解时，对信号进行低通和高通滤波，将信号分解为位于不同频带和时段内。

进一步地，所述对分解后的脉象信号进行去除噪声时，引入门限来作为甄别受到噪声污染的小波系数，等于和小于门限的小波系数为由噪声产生，置其为零而舍去；对于大于门限的小波系数为是含有用信号成分，给予保留。

进一步地，所述对去除噪声后的信号消除基线漂移时，由于基线漂移的成分为缓变趋势分量，在小波分解中会直接显现较大的尺度，只要在重构过程中将这一尺度下的分量直接去除，即实现基线矫正、恢复去除基线漂移后的脉象信号。

还有，本发明还提供了一种脉象信号预处理装置，包括：

分解单元，用于对采集的脉象信号通过小波变换进行信号分解；

去噪单元，与所述分解单元连接，用于对分解后的脉象信号进行去除噪声；

除漂移单元，与所述去噪单元连接，用于对去除噪声后的信号消除基线漂移。

可选地，所述分解单元对信号进行低通和高通滤波，将信号分解为位于不同频带和时段内。

进一步地，所述去噪单元引入门限来作为甄别受到噪声污染的小波系数，等于和小于门限的小波系数为由噪声产生，置其为零而舍去；对于大于门限的小波系数为是含有用信号成分，给予保留。

进一步地，所述除漂移单元对去除噪声后的信号消除基线漂移时，由于基线漂移的成分为缓变趋势分量，在小波分解中会直接显现较大的尺度，只要在重构过程中将这一尺度下的分量直接去除，即实现基线矫正、恢复去除基线漂移后的脉象信号。

从上面所述可以看出，本发明提供的脉象信号预处理方法和装置，通过对采集的脉象信号通过小波变换进行信号分解，对分解后的脉象信号进行去除噪声，然后对去除噪声后的信号消除基线漂移。从而，本发明能够消除脉象信号中的干扰和噪声，有效提高脉象信号时域特征提取的精度和准确度。

附图说明

图1为本发明实施例一种脉象信号预处理方法的流程示意图；

图2为本发明实施例Daubechies小波函数系中的db8滤波器处理后的示意图；

图3为本发明实施例Mallat算法多分辨率分解树的示意图；

图4为本发明实施例小波包分析的三层分解树的示意图；

图5为本发明实施例一种脉象信号预处理装置的结构示意图；

图6为本发明实施例对三峰平脉信号施加50Hz工频干扰的预处理结果的示意图；

图7为本发明实施例对三峰平脉信号施加线性上升基线漂移和肌电干扰的预处理结果的示意图；

图8为本发明实施例对三峰平脉信号施加正弦基线漂移和白噪声的预处理结果的示意图；

图9为本发明实施例对两峰滑脉信号施加50Hz工频干扰的预处理结果的示意图；

图10为本发明实施例对两峰滑脉信号施加线性下降基线漂移和肌电干扰的预处理结果的示意图；

图11为本发明实施例对两峰滑脉信号施加余弦基线漂移和白噪声的预处理结果的示意图；

图12为本发明实施例对浮相弦脉信号施加50Hz工频干扰的预处理结果的示意图；

图13为本发明实施例对浮相弦脉信号仅施加肌电干扰的预处理结果的示意图；

图14为本发明实施例对浮相弦脉信号仅施加白噪声的预处理结果的示意图；

图15为本发明实施例对平顶涩脉信号施加50Hz工频干扰的预处理结果的示意图；

图16为本发明实施例对平顶涩脉信号仅施加肌电干扰的预处理结果的示意图；

图17为本发明实施例对平顶涩脉信号仅施加白噪声的预处理结果的示意图；

图18为本发明实施例对无胃气脉象信号不添加干扰直接进行预处理的结果示意图；

图19为本发明实施例对无胃气脉象信号仅施加肌电干扰的预处理结果示意图；

图20为本发明实施例对无胃气脉象信号仅施加白噪声的预处理结果示意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。

参阅图1所示，为本发明一种脉象信号预处理方法一个实施例的流程示意图，包括：

步骤101，对采集的脉象信号通过小波变换进行信号分解。具体实施过程如下：

作为本发明的一个实施例，小波分析作为一种时频局部化分析方法，其基本思想是寻求一个满足一定条件的基本小波函数ψ(t)，通过对基本小波函数的伸缩和平移构成小波函数族Ψu,s(t)，然后利用这个函数族来逼近或表示所要研究的信号，并作出相应的分析和处理。

对于函数ψ(t)∈L2(R)(能量有限空间),若满足∫Rψ(t)dt＝0则称之为小波函数。引入参数u、s对ψ(t)进行伸展和平移，有：

${\mathrm{\&psi;}}_{u,s}=\frac{1}{\sqrt{\left|s\right|}}\mathrm{\&psi;}\left(\frac{t-u}{s}\right)$

ψu,s称其为伸展小波。其中s为尺度因子，u为平移因子，于是函数f(t)∈L²(R)关于小波函数ψ(t)的小波变换为：

$W_f=(u,s)=<f,{\mathrm{\&psi;}}_{u,s}>={\&Integral;}_{R}f\left(t\right)\stackrel{\&OverBar;}{{\mathrm{\&Psi;}}_{u,s}\left(t\right)}\mathrm{dt}$

设ψ(t)的Fourier变换为Ψ(ω),当小波函数满足如下关系式：

$C_{\mathrm{\&Psi;}}=2\mathrm{\&pi;}{\&Integral;}_R\frac{|\mathrm{\&Psi;}\left(\mathrm{\&omega;}\right){|}^2}{\left|\mathrm{\&omega;}\right|}\mathrm{d\&omega;}<+\&infin;$

小波变换是可逆的，且有如下逆变换关系式：

$f\left(t\right)=\frac{1}{C_{\mathrm{\&psi;}}}{\&Integral;\&Integral;}_{R^2}\frac{1}{s^2}{W}_f(u,s){\mathrm{\&psi;}}_{u,s}\left(t\right)\mathrm{duds}$

同窗口傅立叶变换一样，小波变换也可以度量频谱成分的时频变化，然而这两种方法在时频平面上的分便率并不相同。小波变换可以写成：

对于由计算机采集获得的离散数字信号f(t)，取其中伸缩步长s₀>1，位移步长u₀≠0,n,m∈Z,由此得到的离散小波：

${\mathrm{\&psi;}}_{m,n}\left(t\right)={s_0}^{-m/2}\mathrm{\&psi;}({s_0}^{-m}t-{\mathrm{nu}}_0)m,n\&Element;Z$

因此得到f(t)的离散小波变换(DWT)：

${\mathrm{WT}}_f(m,n)={\&Integral;}_{R}f\left(t\right)\stackrel{\&OverBar;}{{\mathrm{\&psi;}}_{m,n}\left(t\right)}\mathrm{dt}$

逆变换为：

$f\left(t\right)=C\underset{\&infin;}{\overset{\&infin;}{\mathrm{\&Sigma;}}}\underset{\&infin;}{\overset{\&infin;}{\mathrm{\&Sigma;}}}{\mathrm{WT}}_f(m,n){\mathrm{\&psi;}}_{m,n}\left(t\right)$

C是一个与信号无关的常数。

从频域上看，用不同尺度做小波变换大致相当于用一组带通滤波器对信号进行处理。当s值小时时轴上观察范围小，而在频域上相当于用较高频率做分辨率较高的分析，即用高频小波做精细观察。当s值较大时时轴上考察范围大，而在频域上相当于用低频小波做概貌观察。生理信号被分解后常表现出高频分量持续时间较短，低频分量持续时间较长的特点，这也正和小波分析的性质互相吻合。

在Mallat算法的多分辨率分析中，表明多分辨率分析是按照不同的尺度因子j把Hilbert空间L²(R)分解为所有子空间W_j的正交和的。其中，W_j为小波函数Ψ(t)的闭包(小波子空间)。现在，我们希望进一步对小波子空间W_j按照二进制分式进行频率的细分，以达到提高频率分辨率的目的。

定义子空间是函数u_n(t)的闭包空间，而是函数u_2n(t)的闭包空间，并令u_n(t)满足下面的二尺度方程：

$u_{2n}\left(t\right)=\sqrt{2}\underset{k\&Element;z}{\mathrm{\&Sigma;}}{h}_k{u}_n(2t-k)$

${\mathrm{ui}}_{2n+1}\left(t\right)=\sqrt{2}\underset{k\&Element;z}{\mathrm{\&Sigma;}}{g}_k{u}_n(2t-k)$

由上式构造的序列{u_n(t)}(其中n∈Z₊)称为由基函数u₀(t)＝φ(t)和u₁(t)＝ψ(t)确定的正交小波包。h_k和g_k为分别为低频滤波器组和高频滤波器组，由尺度函数和小波函数构造如下：

$h_k=<\mathrm{\&phi;}(t/2),\mathrm{\&phi;}\left(t\right)>=\frac{1}{\sqrt{2}}\&Integral;\mathrm{\&phi;}(t/2){\mathrm{\&phi;}}^{*}(t-k)\mathrm{dt}$

$g_k=<\mathrm{\&psi;\&phi;}(t/2),\mathrm{\&phi;}\left(t\right)>=\frac{1}{\sqrt{2}}\&Integral;\mathrm{\&psi;\&phi;}(t/2){\mathrm{\&phi;}}^{*}(t-k)\mathrm{dt}$

Daubechies小波函数系构造满足正规性要求的低通滤波器H(z)的方法如下：

令 $H\left(z\right)=\sqrt{2}[\frac{1+z^{-1}}{2}{]}^p{F}_0\left(z\right)$

也就是 $\left(\mathrm{\&omega;}\right)\sqrt{2}\frac{1+z^{-\mathrm{j\&omega;}}}{2}{]}^p{F}_0\left(e^{\mathrm{j\&omega;}}\right)$

规定|F₀(ejω)|²取cosω多项式：

$|F_0\left(e^{\mathrm{j\&omega;}}\right){|}^2=\underset{j=0}{\overset{p-1}{\mathrm{\&Sigma;}}}{(}_j^{p-1+j})(\frac{1-\cos }{2}{)}^j$

式中是从p-1+j中每次取j可能得到的组合数目。

当p给定|F₀(e^jω)|²就可确定，又由于|F₀(e^jω)|²的相应Z变换是F₀(z)F₀(z^-1)，因此从其Z变换的每对互为倒数的零点中取一个便可以得到F₀(z)，从而可以求得低通滤波器H(z)。进一步求得高通滤波器G(z)的方法如下：

二尺度方程的频域形式为：

Φ(2ω)＝H(ω)Φ(ω)

Ψ(2ω)＝G(ω)Φ(ω)

由二尺度方程的频域形式引出的递推形式为：

$\mathrm{\&Phi;}\left(\mathrm{\&omega;}\right)=\underset{i=1}{\overset{\&infin;}{\&prod;}}H\left(2^{-i}\mathrm{\&omega;}\right)$

$\mathrm{\&Psi;}\left(\mathrm{\&omega;}\right)=G\left(\frac{\mathrm{\&omega;}}{2}\right)\underset{i=2}{\overset{\&infin;}{\&prod;}}H\left(2^{-i}\mathrm{\&omega;}\right)$

根据空间二剖分的能量守恒性：

|Φ(ω)|²＝|Φ(2ω)|²+|Ψ(2ω)|²

相应的滤波器表示为：

|H(ω)|²+|G(ω)|²＝1

或：H(z)H(z^-1)+G(z)G(z^-1)＝1

由此可以根据求得的低通滤波器H(z)确定高通滤波器G(z)。对于Daubechies小波函数系中的db8小波构造的h_k和g_k滤波器如图2所示。

应用小波包分解滤波器组对第p层信号进行滤波，相应的离散平滑信号和离散细节信号如下：

$a_k^{p+1}=\underset{n}{\mathrm{\&Sigma;}}{h}_{n-2k}{s}_n^p$

$d_k^{p+1}=\underset{n}{\mathrm{\&Sigma;}}{g}_{n-2k}{s}_n^p$

其中对应Mallat算法多分辨率分析的部分如图3所示。对于小波包的三层分解，其小波包的分解树如图4所示。图中S表示信号，A表示低频，D表示高频，末尾的序号表示小波包分解层数。可见小波包分析Mallat算法多分辨率分析基础上对高频部分的进一步细分。

步骤102，对分解后的脉象信号进行去除噪声，具体的实施过程如下：

在本发明的一个实施例中，数据测量经常受到各种噪声信号的干扰和影响，含干扰的脉象信号可表示为

X(n)＝f(n)+W(n)

其中X(n)是实测的脉象信号，具有较强的随机性和背景噪声，而且属于非线性、非平稳的微弱信号。f(n)是不含噪声的脉象信号，W(n)是各种干扰信号总和，即所谓加性噪声，如50/60Hz工频干扰、肌电干扰等。

通过一个决策算子D来变换含噪声的数据X，可以估计信号：

我们的目标是最小化估计误差，而这个误差可以由损失函数度量：

r(D,f)＝E{||f-DX||²}

基于小波变换分析的多分辨分析即相当于对信号进行低通和高通滤波，可将信号分解为位于不同频带和时段内的各个成分。如果对一个只是含高频噪声的信号进行小波分解，则可以看出高频系数的幅值随着分解层次的增加而很快的衰减，并且高频系数的方差也很快的衰减。

通过小波多分辨分析算法将信号分解后，引入门限来作为甄别受到噪声污染的小波系数。将等于和小于门限的小波系数认为由噪声产生，置其为零而舍去。对于大于门限的小波系数,即认为是含有用信号成分，给予保留。再由小波多分辨分析重构算法根据形成新的信号成分序列来重建信号,从而既获得滤除噪声后的信号，又不致于引起重建结果的明显失真。这就是非线性小波方法用于从噪声中恢复信号的实质。要用小波方法很好地实现信噪分离，关键的问题是如何设计出好的门限。

对于滤波器组分解离散信号时，小波基阀值计算下的信号估计为：

式中ρ_T为小波阀值函数，定义离散正交小波基为ψ_j,m(n)＝ψ_j(n-N2^jm)，信号的支集规范为[0,1]，且有间隔为N^-1的N个样本，因而尺度参数2^j从2^L＝N^-1变化到2^j<1。

小波阀值函数ρ_T的选取原则为

(1)对于高频噪声信号采用硬取阀值

${\mathrm{\&rho;}}_T\left(x\right)=\left\{\begin{array}{cc}x& \left|x\right|>T_j,j<j_h\\ 0& \left|x\right|>T_j,j<j_h\end{array}\right.$

(2)对于中间频率含噪声信号采用软取阀值

${\mathrm{\&rho;}}_T\left(x\right)=\left\{\begin{array}{cc}x-T_j& \left|x\right|>T_j,j_h<j<j_l\\ x+T_j& \left|x\right|>T_j,j_h<j<j_l\\ 0& \left|x\right|>T_j,j_h<j<j_l\end{array}\right.$

(3)对于低频信号部分保持不变

ρ_T(x)＝x   j>j_l

上面式中T_j为对应小波多分辨分析j尺度下的小波阀值，信号高频尺度边界j_h和低频尺度边界j_l，可以根据信号频谱分布的先验知识确定。

要保证较小的阀值估计风险，小波阀值T_j采用如下方法确定：

$T_j={\mathrm{\&sigma;}}_j\sqrt{{2\log }_e{N}_j}$

式中σ_j＝median(abs(d{j}))/0.6745，为小波j尺度下的噪声分解系数标准方差，可以通过小波系数中位近似估计，N_j应小波j尺度下的小波分解系数长度。

为进一步减少阀值估计风险，可以对信号的估计其所有平移且在逆平移后取平均：

上述算法把低分辨率(大尺度)下的小波变换全部保留，高分辨率(小尺度)下的小波变换则只有被确认为边沿附近的各点才予以保留，其余的都加去除。由于噪声的小波变换主要集中在小尺度各层次中，因此经上述处理后，噪声基本去除而边沿信息得以较好保留。

步骤103，对去除噪声后的信号消除基线漂移，具体的实施过程如下：

作为本发明的一个实施例，在信号噪声消除的基础上，可以进一步消除信号的基线干扰。在信号噪声消除过程中，我们利用小波变换多尺度多分辨的特点，即利用小波变换的带通滤波特性和尺度函数的低通滤波特性，将脉象信号进行多尺度小波分解，由于基线漂移的主要成分为缓变趋势分量，在小波分解中会直接显现于某较大的尺度下，只要在重构过程中将这一尺度下的分量直接去除，即可实现基线矫正、恢复去除基线漂移后的脉象信号。在信号的采样频率不变的情况下，由于对应于某一确定的小波变换，其不同尺度下的频窗中心和窗宽是确定的，由此可确定相应去处基线漂移的最大分解尺度。

参阅图5所示，为本发明实施例一种脉象信号预处理装置的结构示意图，根据上面的描述所述的脉象信号预处理装置可以包括：

分解单元501，能够对采集的脉象信号通过小波变换进行信号分解。作为本发明的一个实施例，小波分析作为一种时频局部化分析方法，其基本思想是寻求一个满足一定条件的基本小波函数ψ(t)，通过对基本小波函数的伸缩和平移构成小波函数族Ψ_u,s(t)，然后利用这个函数族来逼近或表示所要研究的信号，并作出相应的分析和处理。

对于函数ψ(t)∈L²(R)(能量有限空间),若满足∫Rψ(t)dt＝0则称之为小波函数。引入参数u、s对ψ(t)进行伸展和平移，有：

${\mathrm{\&psi;}}_{u,s}=\frac{1}{\sqrt{\left|s\right|}}\mathrm{\&psi;}\left(\frac{t-u}{s}\right)$

ψ_u,s称其为伸展小波。其中s为尺度因子，u为平移因子，于是函数f(t)∈L²(R)关于小波函数ψ(t)的小波变换为：

$W_f=(u,s)=<f,{\mathrm{\&psi;}}_{u,s}>={\&Integral;}_{R}f\left(t\right)\stackrel{\&OverBar;}{{\mathrm{\&Psi;}}_{u,s}\left(t\right)}\mathrm{dt}$

设ψ(t)的Fourier变换为Ψ(ω),当小波函数满足如下关系式：

$C_{\mathrm{\&Psi;}}=2\mathrm{\&pi;}{\&Integral;}_R\frac{|\mathrm{\&Psi;}\left(\mathrm{\&omega;}\right){|}^2}{\left|\mathrm{\&omega;}\right|}\mathrm{d\&omega;}<+\&infin;$

小波变换是可逆的，且有如下逆变换关系式：

$f\left(t\right)=\frac{1}{C_{\mathrm{\&psi;}}}{\&Integral;\&Integral;}_{R^2}\frac{1}{s^2}{W}_f(u,s){\mathrm{\&psi;}}_{u,s}\left(t\right)\mathrm{duds}$

同窗口傅立叶变换一样，小波变换也可以度量频谱成分的时频变化，然而这两种方法在时频平面上的分便率并不相同。小波变换可以写成：

对于由计算机采集获得的离散数字信号f(t)，取其中伸缩步长s₀>1，位移步长u₀≠0,n,m∈Z,由此得到的离散小波：

ψ_m,n(t)＝s₀ ^-m/2ψ(s₀ ^-mt-nu₀)   m,n∈Z

因此得到f(t)的离散小波变换(DWT)：

${\mathrm{WT}}_f(m,n)={\&Integral;}_{R}f\left(t\right)\stackrel{\&OverBar;}{{\mathrm{\&psi;}}_{m,n}\left(t\right)}\mathrm{dt}$

逆变换为：

$f\left(t\right)=C\underset{\&infin;}{\overset{\&infin;}{\mathrm{\&Sigma;}}}\underset{\&infin;}{\overset{\&infin;}{\mathrm{\&Sigma;}}}{\mathrm{WT}}_f(m,n){\mathrm{\&psi;}}_{m,n}\left(t\right)$

C是一个与信号无关的常数。

从频域上看，用不同尺度做小波变换大致相当于用一组带通滤波器对信号进行处理。当s值小时时轴上观察范围小，而在频域上相当于用较高频率做分辨率较高的分析，即用高频小波做精细观察。当s值较大时时轴上考察范围大，而在频域上相当于用低频小波做概貌观察。生理信号被分解后常表现出高频分量持续时间较短，低频分量持续时间较长的特点，这也正和小波分析的性质互相吻合。

在Mallat算法的多分辨率分析中，，表明多分辨率分析是按照不同的尺度因子j把Hilbert空间L²(R)分解为所有子空间W_j的正交和的。其中，W_j为小波函数Ψ(t)的闭包(小波子空间)。现在，我们希望进一步对小波子空间W_j按照二进制分式进行频率的细分，以达到提高频率分辨率的目的。

定义子空间是函数u_n(t)的闭包空间，而是函数u_2n(t)的闭包空间，并令u_n(t)满足下面的二尺度方程：

$u_{2n}\left(t\right)=\sqrt{2}\underset{k\&Element;z}{\mathrm{\&Sigma;}}{h}_k{u}_n(2t-k)$

${\mathrm{ui}}_{2n+1}\left(t\right)=\sqrt{2}\underset{k\&Element;z}{\mathrm{\&Sigma;}}{g}_k{u}_n(2t-k)$

由上式构造的序列{u_n(t)}(其中n∈Z₊)称为由基函数u₀(t)＝φ(t)和u₁(t)＝ψ(t)确定的正交小波包。h_k和g_k为分别为低频滤波器组和高频滤波器组，由尺度函数和小波函数构造如下：

$h_k=<\mathrm{\&phi;}(t/2),\mathrm{\&phi;}\left(t\right)>=\frac{1}{\sqrt{2}}\&Integral;\mathrm{\&phi;}(t/2){\mathrm{\&phi;}}^{*}(t-k)\mathrm{dt}$

$g_k=<\mathrm{\&psi;\&phi;}(t/2),\mathrm{\&phi;}\left(t\right)>=\frac{1}{\sqrt{2}}\&Integral;\mathrm{\&psi;\&phi;}(t/2){\mathrm{\&phi;}}^{*}(t-k)\mathrm{dt}$

Daubechies小波函数系构造满足正规性要求的低通滤波器H(z)的方法如下：

令 $H\left(z\right)=\sqrt{2}[\frac{1+z^{-1}}{2}{]}^p{F}_0\left(z\right)$

也就是 $\left(\mathrm{\&omega;}\right)\sqrt{2}\frac{1+z^{-\mathrm{j\&omega;}}}{2}{]}^p{F}_0\left(e^{\mathrm{j\&omega;}}\right)$

规定|F₀(e^jω)|²取cosω多项式：

$|F_0\left(e^{\mathrm{j\&omega;}}\right){|}^2=\underset{j=0}{\overset{p-1}{\mathrm{\&Sigma;}}}{(}_j^{p-1+j})(\frac{1-\cos }{2}{)}^j$

式中是从p-1+j中每次取j可能得到的组合数目。

当p给定|F₀(e^jω)|²就可确定，又由于|F₀(e^jω)|²的相应Z变换是F₀(z)F₀(z^-1)，因此从其Z变换的每对互为倒数的零点中取一个便可以得到F₀(z)，从而可以求得低通滤波器H(z)。进一步求得高通滤波器G(z)的方法如下：

二尺度方程的频域形式为：

Φ(2ω)＝H(ω)Φ(ω)

Ψ(2ω)＝G(ω)Φ(ω)

由二尺度方程的频域形式引出的递推形式为：

$\mathrm{\&Phi;}\left(\mathrm{\&omega;}\right)=\underset{i=1}{\overset{\&infin;}{\&prod;}}H\left(2^{-i}\mathrm{\&omega;}\right)$

$\mathrm{\&Psi;}\left(\mathrm{\&omega;}\right)=G\left(\frac{\mathrm{\&omega;}}{2}\right)\underset{i=2}{\overset{\&infin;}{\&prod;}}H\left(2^{-i}\mathrm{\&omega;}\right)$

根据空间二剖分的能量守恒性：

|Φ(ω)|²＝|Φ(2ω)|²+|Ψ(2ω)|²

相应的滤波器表示为：

|H(ω)|²+|G(ω)|²＝1

或：H(z)H(z^-1)+G(z)G(z^-1)＝1

由此可以根据求得的低通滤波器H(z)确定高通滤波器G(z)。对于Daubechies小波函数系中的db8小波构造的h_k和g_k滤波器如图2所示。

应用小波包分解滤波器组对第p层信号进行滤波，相应的离散平滑信号和离散细节信号如下：

$a_k^{p+1}=\underset{n}{\mathrm{\&Sigma;}}{h}_{n-2k}{s}_n^p$

$d_k^{p+1}=\underset{n}{\mathrm{\&Sigma;}}{g}_{n-2k}{s}_n^p$

其中对应Mallat算法多分辨率分析的部分如图3所示。对于小波包的三层分解，其小波包的分解树如图4所示。图中S表示信号，A表示低频，D表示高频，末尾的序号表示小波包分解层数。可见小波包分析Mallat算法多分辨率分析基础上对高频部分的进一步细分。

去噪单元502，与分解单元501连接，能够对分解后的脉象信号进行去除噪声。

在本发明的一个实施例中，数据测量经常受到各种噪声信号的干扰和影响，含干扰的脉象信号可表示为

X(n)＝f(n)+W(n)

其中X(n)是实测的脉象信号，具有较强的随机性和背景噪声，而且属于非线性、非平稳的微弱信号。f(n)是不含噪声的脉象信号，W(n)是各种干扰信号总和，即所谓加性噪声，如50/60Hz工频干扰、肌电干扰等。

通过一个决策算子D来变换含噪声的数据X，可以估计信号：

我们的目标是最小化估计误差，而这个误差可以由损失函数度量：

r(D,f)＝E{||f-DX||²}

基于小波变换分析的多分辨分析即相当于对信号进行低通和高通滤波，可将信号分解为位于不同频带和时段内的各个成分。如果对一个只是含高频噪声的信号进行小波分解，则可以看出高频系数的幅值随着分解层次的增加而很快的衰减，并且高频系数的方差也很快的衰减。

通过小波多分辨分析算法将信号分解后，就可根据先验知识，引入门限来作为甄别受到噪声污染的小波系数。将等于和小于门限的小波系数认为由噪声产生，置其为零而舍去。对于大于门限的小波系数,即认为是含有用信号成分，给予保留。再由小波多分辨分析重构算法根据形成新的信号成分序列来重建信号,从而既获得滤除噪声后的信号，又不致于引起重建结果的明显失真。这就是非线性小波方法用于从噪声中恢复信号的实质。要用小波方法很好地实现信噪分离，关键的问题是如何设计出好的门限。

对于滤波器组分解离散信号时，小波基阀值计算下的信号估计为：

式中ρ_T为小波阀值函数，定义离散正交小波基为ψ_j,m(n)＝ψ_j(n-N2^jm)，信号的支集规范为[0,1]，且有间隔为N^-1的N个样本，因而尺度参数2^j从2^L＝N^-1变化到2^j<1。

小波阀值函数ρ_T的选取原则为

(1)对于高频噪声信号采用硬取阀值

${\mathrm{\&rho;}}_T\left(x\right)=\left\{\begin{array}{cc}x& \left|x\right|>T_j,j<j_h\\ 0& \left|x\right|>T_j,j<j_h\end{array}\right.$

(2)对于中间频率含噪声信号采用软取阀值

${\mathrm{\&rho;}}_T\left(x\right)=\left\{\begin{array}{cc}x-T_j& \left|x\right|>T_j,j_h<j<j_l\\ x+T_j& \left|x\right|>T_j,j_h<j<j_l\\ 0& \left|x\right|>T_j,j_h<j<j_l\end{array}\right.$

(3)对于低频信号部分保持不变

ρ_T(x)＝x   j>j_l

上面式中T_j为对应小波多分辨分析j尺度下的小波阀值，信号高频尺度边界j_h和低频尺度边界j_l，可以根据信号频谱分布的先验知识确定。

要保证较小的阀值估计风险，小波阀值T_j采用如下方法确定：

$T_j={\mathrm{\&sigma;}}_j\sqrt{{2\log }_e{N}_j}$

式中σj＝median(abs(d{j}))/0.6745，为小波j尺度下的噪声分解系数标准方差，可以通过小波系数中位近似估计，N_j应小波j尺度下的小波分解系数长度。

为进一步减少阀值估计风险，可以对信号的估计其所有平移且在逆平移后取平均：

上述算法把低分辨率(大尺度)下的小波变换全部保留，高分辨率(小尺度)下的小波变换则只有被确认为边沿附近的各点才予以保留，其余的都加去除。由于噪声的小波变换主要集中在小尺度各层次中，因此经上述处理后，噪声基本去除而边沿信息得以较好保留。

除漂移单元503，与去噪单元502连接，能够对去除噪声后的信号消除基线漂移。

作为本发明的一个实施例，在信号噪声消除的基础上，可以进一步消除信号的基线干扰。在信号噪声消除过程中，我们利用小波变换多尺度多分辨的特点，即利用小波变换的带通滤波特性和尺度函数的低通滤波特性，将脉象信号进行多尺度小波分解，由于基线漂移的主要成分为缓变趋势分量，在小波分解中会直接显现于某较大的尺度下，只要在重构过程中将这一尺度下的分量直接去除，即可实现基线矫正、恢复去除基线漂移后的脉象信号。在信号的采样频率不变的情况下，由于对应于某一确定的小波变换，其不同尺度下的频窗中心和窗宽是确定的，由此可确定相应去处基线漂移的最大分解尺度。

另外，本发明对典型三峰平脉、二峰滑脉、浮相弦脉、峰平涩脉，以及脉律和脉力不均匀的“无胃气”脉，进行信号预处理的结果如下：

对三峰平脉的信号预处理结果如图6至图8所示(图中f(t)为原始信号，X(t)为添加干扰的信号，～f(t)为消除干扰后的信号)。图6给出对信号施加50Hz工频干扰的预处理结果，图7给出对信号施加线性上升基线漂移和肌电干扰的预处理结果，图8给出对信号施加正弦基线漂移和白噪声的预处理结果。由上面图中可见，对三峰平脉的信号预处理可以有效消除基线漂移和噪声，预处理后的信号细节大部分得到保留，信号基本不失真。

对两峰滑脉的信号预处理结果如图9至图11所示(图中f(t)为原始信号，X(t)为添加干扰的信号，～f(t)为消除干扰后的信号)。图9给出对信号施加50Hz工频干扰的预处理结果，图10给出对信号施加线性下降基线漂移和肌电干扰的预处理结果，图11给出对信号施加余弦基线漂移和白噪声的预处理结果。由上面图中可见，对两峰滑脉的信号预处理可以有效消除基线漂移和噪声，预处理后的信号细节大部分得到保留，信号基本不失真。

对浮相弦脉的信号预处理结果如图12至图14所示(图中f(t)为原始信号，X(t)为添加干扰的信号，～f(t)为消除干扰后的信号)。图12给出对信号施加50Hz工频干扰的预处理结果，图13给出对信号仅施加肌电干扰的预处理结果，图14给出对信号仅施加白噪声的预处理结果。由上面图中可见，对浮相弦脉的信号预处理可以有效消除各种噪声，预处理后的信号细节大部分得到保留，信号基本不失真。

对平顶涩脉的信号预处理结果如图15至图17所示(图中f(t)为原始信号，X(t)为添加干扰的信号，～f(t)为消除干扰后的信号)。图15给出对信号施加50Hz工频干扰的预处理结果，图16给出对信号仅施加肌电干扰的预处理结果，图17给出对信号仅施加白噪声的预处理结果。由上面图中可见，对浮相弦脉的信号预处理可以有效消除各种噪声，预处理后的信号细节大部分得到保留，信号基本不失真。

对无胃气脉象的信号预处理结果如图18至图20所示(图中f(t)为原始信号，X(t)为添加干扰的信号，～f(t)为消除干扰后的信号)。图18给出对信号不添加干扰直接进行预处理的结果，图19给出对信号仅施加肌电干扰的预处理结果，图20给出对信号仅施加白噪声的预处理结果。由上面图中可见，对无胃气脉象的信号预处理，在不添加任何干扰时，信号预处理对原始信号基本是无损的，对添加干扰的信号，可以基本可以消除各种噪声，预处理后的信号细节大部分得到保留，除主波峰顶会略微尖化外，信号其他部分基本不失真。

从上面的描述可以看出，本发明实现的一种脉象信号时频处理方法和装置，创造性地提出了通过对采集的脉象信号通过小波变换进行信号分解，对分解后的脉象信号进行去除噪声，之后对去除噪声后的信号消除基线漂移；并且，根据脉象时频分布特点，采用软阀值方法进行脉象信号的噪声消除，可以有效消除工频干扰、肌电干扰和白噪声，噪声消除基本不损失信号细节；信号预处理受脉象测量条件、不同人生理特征差异性影响较小；与此同时，在脉象消噪过程中，可以同时有效去除脉象信号的基线漂移。噪声消除和基线漂移消除在小波多分辨分解和综合过程中一次完成，简化了通常首先进行基线漂移消除、然后进行信号消噪的信号预处理过程；最后，整个所述的脉象信号时频处理方法和装置简便、紧凑，易于实现。

所属领域的普通技术人员应当理解：以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种脉象信号预处理方法，其特征在于，包括以下步骤：

对采集的脉象信号通过小波变换进行信号分解；

对分解后的脉象信号进行去除噪声；

对去除噪声后的信号消除基线漂移。

2.根据权利要求1所述的预处理方法，其特征在于，所述对采集的脉象信号通过小波变换进行信号分解时，对信号进行低通和高通滤波，将信号分解为位于不同频带和时段内。

3.根据权利要求2所述的预处理方法，其特征在于，所述对分解后的脉象信号进行去除噪声时，引入门限来作为甄别受到噪声污染的小波系数，等于和小于门限的小波系数为由噪声产生，置其为零而舍去；对于大于门限的小波系数为是含有用信号成分，给予保留。

4.根据权利要求2或3所述的预处理方法，其特征在于，所述对去除噪声后的信号消除基线漂移时，由于基线漂移的成分为缓变趋势分量，在小波分解中会直接显现较大的尺度，只要在重构过程中将这一尺度下的分量直接去除，即实现基线矫正、恢复去除基线漂移后的脉象信号。

5.一种脉象信号预处理装置，其特征在于，包括：

分解单元，用于对采集的脉象信号通过小波变换进行信号分解；

去噪单元，与所述分解单元连接，用于对分解后的脉象信号进行去除噪声；

除漂移单元，与所述去噪单元连接，用于对去除噪声后的信号消除基线漂移。

6.根据权利要求5所述的预处理装置，其特征在于，所述分解单元对信号进行低通和高通滤波，将信号分解为位于不同频带和时段内。

7.根据权利要求6所述的预处理装置，其特征在于，所述去噪单元引入门限来作为甄别受到噪声污染的小波系数，等于和小于门限的小波系数为由噪声产生，置其为零而舍去；对于大于门限的小波系数为是含有用信号成分，给予保留。

8.根据权利要求6或7所述的预处理装置，其特征在于，所述除漂移单元对去除噪声后的信号消除基线漂移时，由于基线漂移的成分为缓变趋势分量，在小波分解中会直接显现较大的尺度，只要在重构过程中将这一尺度下的分量直接去除，即实现基线矫正、恢复去除基线漂移后的脉象信号。