技术领域technical field
本发明涉及中医脉象测量技术领域，特别是指一种脉象信号的提取方法和装置。The invention relates to the technical field of pulse condition measurement in traditional Chinese medicine, in particular to a method and device for extracting pulse condition signals.
背景技术Background technique
在生物医学信号检测技术中，脉象波检测的误差主要不仅仅是来源于肌电干扰、高频噪声、基线漂移等干扰，而且还来源于实际中脉象波形不够典型。脉象信号是时变、非线性的低频弱信号，脉象信号检测时受压力变化、测量装置干扰等多种因素影响，同一点的脉象幅度和脉形走势每次测量往往都会发生变化，很难有一个确定的基准进行特征点定位。In biomedical signal detection technology, the error of pulse wave detection is not only caused by interference such as myoelectric interference, high-frequency noise, and baseline drift, but also from the fact that the actual pulse wave is not typical. The pulse signal is a time-varying, non-linear low-frequency weak signal. The detection of the pulse signal is affected by various factors such as pressure change and interference of the measuring device. A definite datum is used to locate the feature points.
传统的脉象特征方法提取主要基于时域或频域提取脉象特征，对于一些时域或频域特征不明显的脉象，脉象特征提取具有很多不确定性，很难准确地进行脉象特征的提取和定量化，因此现有的脉象特征提取方法不能够准确处理一些特征不明显的时变弱信号。Traditional pulse feature extraction methods are mainly based on time domain or frequency domain extraction of pulse features. For some pulse features that are not obvious in time domain or frequency domain, pulse feature extraction has many uncertainties, and it is difficult to accurately extract and quantify pulse features. Therefore, the existing pulse condition feature extraction methods cannot accurately deal with some time-varying weak signals with inconspicuous features.
发明内容Contents of the invention
有鉴于此，本发明的目的在于提出一种脉象信号的提取方法和装置，能够实现一种简单、快速和准确的提取方法。In view of this, the object of the present invention is to propose a pulse condition signal extraction method and device, which can realize a simple, fast and accurate extraction method.
基于上述目的本发明提供的脉象信号的提取方法，包括以下步骤：The extraction method of the pulse condition signal that the present invention provides based on above-mentioned purpose, comprises the following steps:
对采集的脉象信号进行小波模极大变换；Carry out wavelet modulus maximal transformation to the collected pulse signal;
对小波模极大变换后的脉象信号进行周期分割；Periodically segment the pulse signal after wavelet modulus maximal transformation;
对变换后的脉象信号的特征点进行定位。The characteristic points of the transformed pulse signal are located.
可选地，所述对采集的脉象信号进行小波模极大变换，包括步骤：Optionally, said carrying out wavelet modulus maximal transformation to the pulse signal of collection, comprises the steps:
采集脉象信号，并得到脉象离散数据；Collect pulse signal and get discrete pulse data;
对脉象离散数据进行高斯一阶导数小波模极大变换；Carry out Gaussian first-order derivative wavelet modulus maximal transformation on discrete pulse condition data;
对脉象离散数据进行高斯二阶导数小波模极大变换。Gaussian second-order derivative wavelet modulo-maximum transformation is performed on discrete pulse condition data.
进一步地，所述对脉象离散数据进行高斯一阶导数小波模极大变换时，为得到脉象离散数据对应的脉象边沿的高斯一阶导数小波变换模极大曲线，需要对噪声产生的模极大曲线进行过滤。Further, when the Gaussian first-order derivative wavelet modulus maximal transformation is carried out to the pulse image discrete data, in order to obtain the Gaussian first-order derivative wavelet transform modulus maximal curve of the pulse image edge corresponding to the pulse image discrete data, the modulus of the noise generation needs to be maximized. Curves are filtered.
进一步地，所述对小波模极大变换后的脉象信号进行周期分割时，首先设置模阀值上限、模阀值下限，以及长度阀值上限和长度阀值下限；然后，保留模大于模阀值上限以及模长度大于长度阀值上限的模极大曲线，即可完成脉象周期分割。Further, when the pulse signal after the wavelet modulo maximal transformation is periodically divided, at first the upper limit of the modulus threshold, the lower limit of the modulus threshold, and the upper limit of the length threshold and the lower limit of the length threshold are set; value upper limit and the modulus maximal curve whose modulus length is greater than the length threshold upper limit, the pulse cycle segmentation can be completed.
进一步地，所述对变换后的脉象信号的特征点进行定位时，采用连续小波变换进行由粗到精的跟踪，通过将由粗到精的模极大连接成模极大曲线，来完成脉象信号的特征点定位。Further, when the feature points of the transformed pulse signal are positioned, continuous wavelet transform is used to track from coarse to fine, and the pulse signal is completed by connecting the modulus maxima from coarse to fine to a modulus maxima curve. feature point location.
还有，本发明还提供了一种脉象信号的提取装置，包括：In addition, the present invention also provides a device for extracting pulse signal, comprising:
变换单元，用于对采集的脉象信号进行小波模极大变换；A transformation unit is used to carry out wavelet modulus maximum transformation to the collected pulse signal;
分割单元，与所述变换单元相连，用于对小波模极大变换后的脉象信号进行周期分割；The segmentation unit is connected with the transformation unit, and is used to carry out periodic segmentation to the pulse signal after the wavelet modulus maximum transformation;
定位单元，与所述分割单元相连，用于对变换后的脉象信号的特征点进行定位。The locating unit is connected with the dividing unit and is used for locating the characteristic points of the transformed pulse signal.
可选地，所述变换单元在对采集的脉象信号进行小波模极大变换时，包括步骤：Optionally, when said transformation unit carries out wavelet modulus maximal transformation to the pulse condition signal of collection, comprise the step:
采集脉象信号，并得到脉象离散数据；Collect pulse signal and get discrete pulse data;
对脉象离散数据进行高斯一阶导数小波模极大变换；Carry out Gaussian first-order derivative wavelet modulus maximal transformation on discrete pulse condition data;
对脉象离散数据进行高斯二阶导数小波模极大变换。Gaussian second-order derivative wavelet modulo-maximum transformation is performed on discrete pulse condition data.
进一步地，所述对脉象离散数据进行高斯一阶导数小波模极大变换时，为得到脉象离散数据对应的脉象边沿的高斯一阶导数小波变换模极大曲线，需要对噪声产生的模极大曲线进行过滤。Further, when the Gaussian first-order derivative wavelet modulus maximal transformation is carried out to the pulse image discrete data, in order to obtain the Gaussian first-order derivative wavelet transform modulus maximal curve of the pulse image edge corresponding to the pulse image discrete data, the modulus of the noise generation needs to be maximized. Curves are filtered.
进一步地，所述分割单元进行对小波模极大变换后的脉象信号进行周期分割时，首先设置模阀值上限、模阀值下限，以及长度阀值上限和长度阀值下限；然后，保留模大于模阀值上限以及模长度大于长度阀值上限的模极大曲线，即可完成脉象周期分割。Further, when the segmentation unit carries out periodical segmentation to the pulse signal after the wavelet modulus maximal transformation, at first the upper limit of the modulus threshold value, the lower limit of the modulus threshold value, and the upper limit of the length threshold value and the lower limit of the length threshold value are set; If the modulus maximum curve is greater than the upper limit of the modulus threshold value and the modulus length is greater than the upper limit of the length threshold value, the pulse period segmentation can be completed.
进一步地，所述定位单元在对变换后的脉象信号的特征点进行定位时，采用连续小波变换进行由粗到精的跟踪，通过将由粗到精的模极大连接成模极大曲线，来完成脉象信号的特征点定位。Further, when the positioning unit locates the characteristic points of the transformed pulse signal, it adopts continuous wavelet transform to carry out the tracking from coarse to fine, and connects the modulus maxima from coarse to fine into a modulus maxima curve to obtain Complete the feature point positioning of the pulse signal.
从上面所述可以看出，本发明提供的脉象信号的提取方法和装置，通过对采集的脉象信号进行小波模极大变换；对小波模极大变换后的脉象信号进行周期分割；对变换后的脉象信号的特征点进行定位。从而，所述的脉象信号的提取方法和装置能够准确处理特征不明显的时变弱信号，并且实现了一种简单、快速和准确的提取方法。As can be seen from the above, the extraction method and device of the pulse signal provided by the present invention, by carrying out wavelet modulus maximum transformation to the pulse signal of collection; Carry out periodic segmentation to the pulse signal after wavelet modulus maximum transformation; The characteristic points of the pulse signal are located. Therefore, the pulse condition signal extraction method and device can accurately process time-varying weak signals with inconspicuous features, and realize a simple, fast and accurate extraction method.
附图说明Description of drawings
图1为本发明一种脉象信号的提取方法的流程示意图；Fig. 1 is the schematic flow sheet of the extracting method of a kind of pulse condition signal of the present invention;
图2为本发明实施例小波模极大变换的方法流程示意图；Fig. 2 is a schematic flow chart of a method for wavelet modulus maximal transformation in an embodiment of the present invention;
图3为本发明实施例过滤后的高斯一阶导数小波变换的模极大曲线的示意图；Fig. 3 is the schematic diagram of the maximal modulus curve of Gaussian first-order derivative wavelet transform after filtering in the embodiment of the present invention;
图4为本发明实施例采用高斯一阶导数小波变换的分形谱的示意图；Fig. 4 is the schematic diagram of the fractal spectrum that adopts Gaussian first-order derivative wavelet transform for the embodiment of the present invention;
图5为本发明实施例过滤后的高斯二阶导数小波变换的模极大曲线的示意图；Fig. 5 is the schematic diagram of the maximal modulus curve of Gaussian second-order derivative wavelet transform filtered by the embodiment of the present invention;
图6为本发明实施例采用高斯二阶导数小波变换的分形谱的示意图；6 is a schematic diagram of a fractal spectrum using Gaussian second-order derivative wavelet transform in an embodiment of the present invention;
图7为本发明实施例脉图的脉象特征参数的示意图；Fig. 7 is the schematic diagram of the characteristic parameter of pulse condition of pulse map of the embodiment of the present invention;
图8为本发明实施例对两峰滑脉采用高斯一阶导数进行的小波变换的示意图；Fig. 8 is a schematic diagram of the wavelet transform performed on the two-peak Huamai using Gaussian first-order derivatives according to an embodiment of the present invention;
图9为本发明实施例对两峰滑脉采用高斯二阶导数进行的小波变换的示意图；Fig. 9 is a schematic diagram of the wavelet transformation performed on the two-peak Huamai using Gaussian second-order derivatives according to an embodiment of the present invention;
图10为本发明一种脉象信号的提取装置的结构示意图。Fig. 10 is a schematic structural diagram of a device for extracting a pulse condition signal according to the present invention.
具体实施方式Detailed ways
为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。In order to make the object, technical solution and advantages of the present invention clearer, the present invention will be described in further detail below in conjunction with specific embodiments and with reference to the accompanying drawings.
参阅图1所示，为本发明一种脉象信号的提取方法的流程示意图，包括：Referring to shown in Fig. 1, be the schematic flow sheet of the extracting method of a kind of pulse condition signal of the present invention, comprise:
步骤101，对采集的脉象信号进行小波模极大变换。具体实施过程如下(如图2所示)：Step 101, performing wavelet modulo-maximum transformation on the collected pulse signal. The specific implementation process is as follows (as shown in Figure 2):
步骤201：采集脉象信号，并得到脉象离散数据。Step 201: Collect pulse signal and obtain discrete pulse data.
作为本发明的一个实施例，可以使用仪器脉搏图形记录仪采集脉象信号，然后通过对脉象信号的标准脉图进行数据化得到脉象数据。较佳地，将脉象信号转化为脉象数据是通过数字笔将脉象图采集到计算机后，以x轴象素间隔dx(mm)和y轴象素间隔dy(mm)为单位进行数据化，离散化后的脉象数据形成脉象标准数据库。优选地，形成脉象标准数据库可以将离散化后的脉象数据通过转化为以实际时间单位表示图形，只需将x轴和y轴分别以实际采样率(25mm/dx)/sec和(5mm/dy)/mv转化即可。As an embodiment of the present invention, the pulse pattern recorder can be used to collect the pulse condition signal, and then the pulse condition data can be obtained by digitizing the standard pulse pattern of the pulse condition signal. Preferably, the pulse signal is converted into pulse data after the pulse image is collected into the computer by a digital pen, and then digitized in units of x-axis pixel interval dx (mm) and y-axis pixel interval dy (mm), discrete The transformed pulse data forms a pulse standard database. Preferably, forming the pulse condition standard database can convert the discretized pulse condition data into a graph expressed in actual time units, and only need to convert the x-axis and y-axis to the actual sampling rate (25mm/dx)/sec and (5mm/dy )/mv conversion.
步骤202，对脉象离散数据进行高斯一阶导数小波模极大变换。其具体的实施过程如下：Step 202, performing Gaussian first-order derivative wavelet modulo-maximum transformation on the pulse condition discrete data. Its specific implementation process is as follows:
首先，为得到脉象离散数据对应的脉象边沿的高斯一阶导数小波变换模极大曲线，需要对噪声产生的模极大曲线进行过滤，方法如下：First of all, in order to obtain the Gaussian first-order derivative wavelet transform modulus maximum curve corresponding to the pulse image edge corresponding to the pulse image discrete data, it is necessary to filter the modulus maximum curve generated by the noise, the method is as follows:
1)计算所有模极大曲线起点的模均值vm＝E(moj)(j＝1,…,n)和方差dm2＝D(moj)(j＝1,…,n),模阀值上限thmh＝vm+dm,模阀值下限thml＝vm-dm。1) Calculate the modulus mean vm =E(moj )(j=1,...,n) and variance dm2 =D(moj )(j=1,...,n) of the starting points of all modulus maximal curves, The upper limit of the modulus threshold thmh =vm +dm , the lower limit of the modulus threshold thml =vm -dm .
2)计算所有模极大曲线长度的均值vl＝E(lmj)(j＝1,…,n)和方差dl2＝D(lmj)(j＝1,…,n)，长度阀值上限thlh＝vl+dl,长度阀值下限thll＝vl-dl。2) Calculate the mean vl =E(lmj )(j=1,...,n) and variance dl2 =D(lmj )(j=1,...,n) of the lengths of all modulus maximal curves, the length The upper limit of the threshold value thlh =vl +dl , the lower limit of the length threshold thll =vl -dl .
3)计算脉象主波1/3宽度区间Xw＝[twb,twe],即距离主波顶部1/3主波高度的宽度区间。其中主波波峰位置为主波极大值位置，令w＝(tp-tb)/3，则主波宽度区间Xw近似为主波顶点为中心的w对称宽度区间。3) Calculate the 1/3 width interval of the main wave of the pulse condition Xw =[twb , twe ], that is, the width interval of 1/3 the height of the main wave from the top of the main wave. Wherein the peak position of the main wave is the position of the maximum value of the main wave, if w=(tp -tb )/3, then the width interval Xw of the main wave is approximately the w symmetric width interval centered on the apex of the main wave.
4)对主波1/3宽度附近内的模极大曲线，为研究峰值附近的小隆起不平，保留模大于模阀值下限thml和模长度大于长度阀值下限thll的模极大曲线。4) For the maximum modulus curve in the vicinity of the 1/3 width of the main wave, in order to study the small bumps near the peak, the maximum modulus curve with a mode greater than the lower limit of the mode threshold value thml and a mode length greater than the lower limit of the length threshold value thll is reserved .
5)对主波1/3宽度区间外的模极大曲线，保留模大小大于模阀值上限thmh和模长度大于长度阀值上限thlh的模极大曲线。5) For the maximal modulus curves outside the 1/3 width interval of the main wave, retain the maximal modulus curves whose modulus size is greater than the upper limit of the modulus threshold thmh and whose modulus length is greater than the upper limit of the length threshold thlh .
6)进一步滤除Lipschitz指数α显著小(α<＝－100)和平滑因子σ显著大(σ>＝0.25)的奇异性模极大曲线。6) Further filter out curves with a significantly smaller Lipschitz exponent α (α<=-100) and a significantly larger smoothing factor σ (σ>=0.25).
在本发明的实施例中，可以用小波描述信号的局部奇异性，奇异点小波模极大衰减性可以通过Lipschitz指数α描述，奇异点磨光平滑性可以通过平滑因子σ描述，奇异点分形性可以通过分形维数D(α)描述。In the embodiment of the present invention, the local singularity of the signal can be described by the wavelet, the extreme attenuation of the wavelet mode of the singular point can be described by the Lipschitz exponent α, the smoothness of the singular point can be described by the smoothing factor σ, and the fractal property of the singular point It can be described by fractal dimension D(α).
其中，小波检测信号奇异点衰减特性用Lipschitz指数来描述：Among them, the attenuation characteristics of the singular point of the wavelet detection signal are described by Lipschitz exponent:
设n是一正整数，n＜α≤n+1，f(t)在t0点的Lipschitz指数α充要条件是存在两个常数A和h0(A＞0，h0＞0)，且有一n次多项式Pn(h)使得对h≤h0有：Suppose n is a positive integer, n<α≤n+1, the Lipschitz exponent α of f(t) at point t0 is a necessary and sufficient condition that there are two constants A and h0 (A>0, h0 >0), And there is an n-degree polynomial Pn (h) such that for h≤h0 :
|f(t0+h)-Pn(t0+h)|≤A|h|α，则称f(t)在点t0处为Lipschitzα。实际上多项式Pn(h)是f(t)在t0点的taylor展开式的前n项：|f(t0 +h)-Pn (t0 +h)|≤A|h|α , then f(t) is called Lipschitzα at point t0 . In fact, the polynomial Pn (h) is the first n terms of the Taylor expansion of f(t) at point t0 :
f(t)＝f(t0)+a1h+a2h2+--+anhn+O(hn+1)＝Pn(t)+O(hn+1)。f(t)=f(t0 )+a1 h+a2 h2 +--+an hn +O(hn+1 )=Pn (t)+O(hn+1 ).
由此可见：如果f(t)n次可微，但n阶导数不连续，因此n+1次不可微，则n<α≤n+1，相对于taylor公式中对信号分解误差的n+1次描述，Lipschitzα正则性通过使用非整数的指数改进了这个误差上界。如果上式对所有t0和t0+h在区间[t1,t2]内均成立，则称f(t)在此区间为一致Lipschitzα。It can be seen that if f(t) is differentiable for n times, but the nth order derivative is discontinuous, so n+1 times are not differentiable, then n<α≤n+1, compared to the n+ of the signal decomposition error in Taylor's formula 1 description, Lipschitz α regularity improves this error upper bound by using non-integer exponents. If the above formula is true for all t0 and t0 +h in the interval [t1 ,t2 ], then f(t) is said to be consistent Lipschitzα in this interval.
若一个函数在一点连续且一次可微，但导数不连续，那么它在这点就是Lipschitz指数为1。若f(t)在点t0处不连续但有界，即间断的情形,则Lipschitz指数为0。Lipschitz指数越大，f(t)越光滑。例如，斜坡函数的Lipschitz指数为1，阶跃函数的Lipschitz指数为0，δ函数的Lipschitz指数为-1，白噪的Lipschitz指数为-0.5-ε(ε>0)。If a function is continuous and once differentiable at a point, but its derivative is discontinuous, then it has a Lipschitz exponent of 1 at that point. If f(t) is discontinuous but bounded at point t0 , that is, discontinuous, the Lipschitz exponent is 0. The larger the Lipschitz exponent, the smoother f(t) is. For example, the Lipschitz exponent of the slope function is 1, the Lipschitz exponent of the step function is 0, the Lipschitz exponent of the delta function is -1, and the Lipschitz exponent of the white noise is -0.5-ε (ε>0).
当f(t)完全是正则函数时，同奇异点一样其小波变换也可以有一系列取模极大的点趋于横坐标点v，因此仅沿着尺度搜索小波模极大对于奇异性检测是不够的，需要从模极大幅值的衰减性计算Lipschitz正则性。When f(t) is completely a regular function, its wavelet transform, like the singular point, can also have a series of points with the maximum modulus tending to the abscissa point v, so only searching for the maximum modulus of the wavelet along the scale is useful for singularity detection Not enough, Lipschitz regularity needs to be calculated from the attenuation of the modulo extremely large value.
当v点是孤立奇异点时，对s<s0，收敛于v的所有模极大包含在如下的锥中：|u-v|≤CsWhen point v is an isolated singular point, for s<s0 , all modulus maxima converging on v are contained in the following cone: |uv|≤Cs
式中C为ψ的紧支集[-C,C]，在尺度－时间平面上点v的影响锥是所有点(u,v)的集合。这意味着f在v的邻域内没有快速振荡，|Wf(u,s)|在v的邻域里的衰减性由含于锥中的模极大衰减性控制。f在v的邻域内是一致Lipschitzα的，并且仅当存在K>0，使得锥中的模极大点(u,v)满足：In the formula, C is the compact support set [-C, C] of ψ, and the influence cone of point v on the scale-time plane is the set of all points (u, v). This means that there is no fast oscillation of f in the neighborhood of v, and the decay of |Wf (u,s)| in the neighborhood of v is controlled by the decay of the modulus maxima contained in the cone. f is uniform Lipschitzα in the neighborhood of v, and only if there exists K>0 such that the modulus maxima (u,v) in the cone satisfy:
上式是Lipschitzα对未磨光的孤立奇异点的小波模极大变化影响，在v点的Lipschitz正则性就是log2|Wf(u,s)|作为log2s的函数沿着收敛于v点的极大曲线的最大斜率。The above formula is the influence of Lipschitzα on the wavelet mode maximum change of the unpolished isolated singular point. The Lipschitz regularity at point v is log2 |Wf (u,s)| as a function of log2 s and converges to v The maximum slope of the maximal curve for a point.
另外，即使无限次连续可微的信号，也可以有很大的变化，如图像中阴影处的边界上，图像的亮度有很快的变化，但由于衍射作用，它却是不连续的。将这些跳跃处光滑化就是用高斯核将之模拟成一个扩散过程，其根方差就是表征边沿平滑特性的平滑因子σ，可以通过小波模极大的衰减性来测量。在峰变处v的邻域内，设f(t)＝f0*gσ(t)，gσ是方差为σ2的高斯函数：In addition, even infinitely continuous differentiable signals can have great changes, such as the image brightness changes rapidly on the boundary of shadows in the image, but it is discontinuous due to the effect of diffraction. To smooth these jumps is to use Gaussian kernel to simulate it as a diffusion process, and its root variance is the smoothing factor σ that characterizes the smoothness of the edge, which can be measured by the extreme attenuation of the wavelet mode. In the neighborhood of v where the peak changes, let f(t)=f0 *gσ (t), gσ is a Gaussian function with variance σ2 :
如果f0在v点的Lipschitzα，则它在v的邻域内是一致Lipschitzα的，当小波函数取为高斯函数的导数时ψ＝(-1)nθ(n)，其中θ(t)＝λexp(-t2/(2β2))，则f在v的邻域内是一致Lipschitzα的，并且仅当存在K>0，使得锥中的模极大点(u,v)满足：If f0 is Lipschitzα at point v, it is consistent Lipschitzα in the neighborhood of v, when the wavelet function is taken as the derivative of the Gaussian function, ψ=(-1)n θ(n) , where θ(t)=λexp (-t2 /(2β2 )), then f is consistent Lipschitzα in the neighborhood of v, and only if there is K>0, so that the modulus maximum point (u,v) in the cone satisfies:
上式对孤立奇异点考虑了磨光平滑性，最后一项表达了对小波模极大变化的影响。方差β2依赖于小波的选取，是预先给定的。对于较大的s>>σ/β，小波变化模按sα+1/2衰减，高斯平均对它的影响还感觉不到。当s≤σ/β，由于高斯平均的影响，f在v点的变化不非常强的关联于s。在这些细的尺度下，小波变换模以sn+1/2级衰减。The above formula considers the smoothness of the smoothness of the isolated singular point, and the last term expresses the influence on the extreme change of the wavelet mode. Variance β2 depends on the selection of wavelet and is given in advance. For a larger s>>σ/β, the wavelet change mode decays according to sα+1/2 , and the influence of Gaussian averaging on it cannot be felt. When s≤σ/β, the change of f at v is not very strongly related to s due to the Gaussian average. At these fine scales, the wavelet transform modulus decays by stepsn+1/2 .
另外，当信号具有非孤立奇异性时，信号离散化后的每一个采样点来自于一个小的时间区间，信号在该区间中含有无穷多个互异的奇异性，由于数值分辨率的有限性，Lipschitz指数的点态度量并不可行，因此这种奇异性分布必须利用多分形的自相似整体度量，以给出了不同Lipschitz正则性的奇异特性的整体分配，这种度量即奇异谱分形维数D(α)。In addition, when the signal has non-isolated singularities, each sampling point after the discretization of the signal comes from a small time interval, and the signal contains infinitely many different singularities in this interval, due to the limited numerical resolution , the point-attitude measure of the Lipschitz exponent is not feasible, so this singularity distribution must use a multifractal self-similar overall measure to give the overall distribution of the singular properties of different Lipschitz regularities, this measure is the singular spectrum fractal dimension Number D(α).
设Sα是使得信号f在其上的Lipschitz正则性为α的所有点的集合，f的奇异谱D(α)定义为Sα的分形维数，即Nα(s)～s-D(α)。式中Nα(s)表示Sα中点的个数，奇异谱给出了在任意尺度s下Lipschitzα奇异性的比例值。当自相似信号的谱D(α)是凸的，其计算过程如下：对每一尺度s计算出Wf(u,s)和模极大，随s的变化链接小波模极大。计算分解函数：然后，利用log2Z(s,q)作为log2s的函数的线性衰减性计算τ(q)，即log2Z(q,s)≈τ(q)log2s+C(q)。最后计算谱D(α)：对D(α)可采用四元组[αmin,αmax,αDmax,Dmax]描述，表示D(α)的支撑域为[αmin,αmax]，当α＝αDmax时D(α)取最大值Dmax(aDmax)。Let Sα be the set of all points on which the Lipschitz regularity of the signal f is α, and the singular spectrum D(α) of f is defined as the fractal dimension of Sα , that is, Nα (s)～s-D( a) . In the formula, Nα (s) represents the number of points in Sα , and the singularity spectrum gives the proportional value of the Lipschitz α singularity at any scale s. When the spectrum D(α) of the self-similar signal is convex, the calculation process is as follows: Wf (u, s) and modulus are calculated for each scale s, and the wavelet modulus is maximized with the change of s. Compute the factorization function: Then, τ(q) is calculated using the linear decay of log2 Z(s,q) as a function of log2 s, ie log2 Z(q,s)≈τ(q)log2 s+C(q). Finally calculate the spectrum D(α): D(α) can be described by the quadruple [αmin ,αmax ,αDmax ,Dmax ], indicating that the support domain of D(α) is [αmin ,αmax ], when α=αDmax D( α) takes the maximum value Dmax (aDmax ).
如图3给出过滤后的高斯一阶导数小波变换的模极大曲线，以及图4给出采用高斯一阶导数小波变换的分形谱，支撑区间为[0.10000.2040]，分形谱位于坐标0.1560的最大值为0.3490。Figure 3 shows the modulus maximum curve of the filtered Gaussian first-order derivative wavelet transform, and Figure 4 shows the fractal spectrum using the Gaussian first-order derivative wavelet transform, the support interval is [0.10000.2040], and the fractal spectrum is located at the coordinate 0.1560 The maximum value of is 0.3490.
步骤203，对脉象离散数据进行高斯二阶导数小波模极大变换。其具体的实施过程如下：Step 203, performing Gaussian second-order derivative wavelet modulus maximal transformation on the pulse condition discrete data. Its specific implementation process is as follows:
为得到脉象离散数据对应的脉象局部极大极小点的高斯二阶导数小波变换模极大曲线，需要对噪声产生的模极大曲线进行过滤，方法如下：In order to obtain the Gaussian second-order derivative wavelet transform modulus max curve of the pulse image local maxima and minima points corresponding to the pulse image discrete data, it is necessary to filter the modulus max curve generated by the noise, the method is as follows:
1)计算所有模极大曲线起点的模均值vm＝E(moj)(j＝1,…,n)和方差dm2＝D(moj)(j＝1,…,n),模阀值上限thmh＝vm+dm,模阀值下限thml＝vm-dm。1) Calculate the modulus mean vm =E(moj )(j=1,...,n) and variance dm2 =D(moj )(j=1,...,n) of the starting points of all modulus maximal curves, The upper limit of the modulus threshold thmh =vm +dm , the lower limit of the modulus threshold thml =vm -dm .
2)计算所有模极大曲线长度的均值vl＝E(lmj)(j＝1,…,n)和方差dl2＝D(lmj)(j＝1,…,n),长度阀值上限thlh＝vl+dl,长度阀值下限thll＝vl-dl。2) Calculate the mean vl =E(lmj )(j=1,...,n) and variance dl2 =D(lmj )(j=1,...,n) of the lengths of all modulus maximal curves, the length The upper limit of the threshold value thlh =vl +dl , the lower limit of the length threshold thll =vl -dl .
3)计算脉象主波1/3宽度区间Xw＝[twb,twe],即距离主波顶部1/3主波高度的宽度区间。其中主波波峰位置为主波波峰极大值位置，令w＝(tp-tb)/3则主波宽度区间Xw近似为主波顶点为中心的w对称宽度区间。3) Calculate the 1/3 width interval of the main wave of the pulse condition Xw =[twb , twe ], that is, the width interval of 1/3 the height of the main wave from the top of the main wave. Wherein the peak position of the main wave is the position of the maximum value of the main wave peak, if w=(tp -tb )/3, the width interval Xw of the main wave is approximate to the w-symmetrical width interval centered on the apex of the main wave.
4)对主波宽度区间内的模极大曲线，保留模大小大于模阀值下限thml和模长度，大于长度阀值下限thll的模极大曲线。4) For the maximal modulus curves in the main wave width interval, retain the maximal modulus curves whose mode size is greater than the lower limit of the mode threshold thml and the mode length, and greater than the lower limit of the length threshold thll .
5)对主波宽度区间外的模极大曲线，保留模大小大于模阀值上限thmh和模长度大于长度阀值上限thlh的模极大曲线。5) For the maximal modulus curves outside the main wave width interval, retain the maximal modulus curves whose mode size is greater than the upper limit of the mode threshold value thmh and whose mode length is greater than the upper limit of the length threshold value thlh .
作为本发明的一个实施例，为了对噪声产生的模极大曲线过滤效果更好，进一步可以滤除Lipschitz指数α显著小(α<＝－100)和平滑因子σ显著大(σ>＝0.25)的奇异性模极大曲线。As an embodiment of the present invention, in order to have a better filtering effect on the maximum modulus curve generated by noise, it is further possible to filter out significantly small Lipschitz exponent α (α<=-100) and significantly large smoothing factor σ (σ>=0.25) The singularity modulus maximal curve of .
如图5给出过滤后的高斯二阶导数小波变换的模极大曲线，图6给出采用高斯二阶导数小波变换的分形谱，支撑区间为[0.1000，0.5640]，分形谱位于坐标0.3160的最大值为0.4004。Figure 5 shows the modulus maximal curve of the filtered Gaussian second-order derivative wavelet transform, and Figure 6 shows the fractal spectrum using the Gaussian second-order derivative wavelet transform. The maximum value is 0.4004.
在本发明的具体实施例中，采用高斯一阶导数小波变换计算的模极大曲线数据如表4-1所示，采用高斯二阶导数小波变换计算的模极大曲线数据如表4-2所示，表中pos_s表示模极大曲线的起始位置，pos_e表示模极大曲线的终止位置，scale表示极大曲线的终止位置对应的尺度，mo表示模极大曲线的起始位置对应的模极大值，sign表示模极大曲线的起始位置对应的模极大值符号，len表示模极大曲线的起止长度，alfapd表示不考虑平滑特征的Lipschitz指数，alfa表示考虑平滑特征的Lipschitz指数，delta表示考虑平滑特征的平滑因子，residual表示采用多尺度近似公式计算多尺度参数的残留误差，result表示对模极大曲线过滤后的结果，其中为1表示为保留的曲线，为0表示被过滤的曲线。In a specific embodiment of the present invention, the modulus maximal curve data calculated by Gaussian first-order derivative wavelet transform is shown in Table 4-1, and the modulus maximal curve data calculated by Gaussian second-order derivative wavelet transform is shown in Table 4-2 As shown in the table, pos_s indicates the starting position of the maximum modulus curve, pos_e indicates the end position of the maximum modulus curve, scale indicates the scale corresponding to the end position of the maximum modulus curve, and mo indicates the corresponding value of the starting position of the maximum modulus curve Modulus maximum value, sign indicates the modulus maximum value sign corresponding to the starting position of the modulus maximum curve, len indicates the start and end length of the modulus maximum curve, alfapd indicates the Lipschitz exponent that does not consider the smooth feature, alfa indicates the Lipschitz that considers the smooth feature Index, delta means the smoothing factor considering smooth features, residual means the multi-scale approximation formula is used to calculate the residual error of multi-scale parameters, result means the result after filtering the modulus maximal curve, where 1 means the reserved curve, and 0 means The filtered curve.
表4-1采用高斯一阶导数小波变换的模极大曲线数据Table 4-1 Modular maximum curve data using Gaussian first derivative wavelet transform
表4-2采用高斯二阶导数小波变换的模极大曲线数据Table 4-2 Modular maximum curve data using Gaussian second derivative wavelet transform
如果将表4-1和表4-2中保留的模极大曲线按照起始点pos_s顺序排列，用“1”表示表4-1中的保留的模极大曲线，用“2”表示表4-2中的保留的模极大曲线，则表4-1和表4-2中的保留的模极大曲线顺序排列如下：If the modulus maximum curves reserved in Table 4-1 and Table 4-2 are arranged in the order of the starting point pos_s, use "1" to represent the reserved modulus maximum curves in Table 4-1, and use "2" to represent Table 4 For the maximum modulus curves reserved in -2, the order of the maximum modulus curves reserved in Table 4-1 and Table 4-2 is as follows:
[1212122221112212112121121111111212]在上面对保留的模极大曲线顺序排列基础上：用“+”表示类型为“1”的模极大曲线中的上升沿(模极大为正)，用“－”表示类型为“1”的模极大曲线中的下降沿(模极大为负)。用“|”表示类型为“2”的模极大曲线中的模极大大于模值上限(mo>thmh)和模长度大于长度上限(l>thlh)的波峰(模符号为正)，用“*”表示类型为“2”的模极大曲线中的模极大小于模值上限大于模值下限(thmh>mo>thml)和模长度小于长度上限但大于长度下限(thlh>l>thll)的波峰。用“_”表示类型为“2”的模极大曲线中的模极大大于模值上(mo>thmh)和模长度大于长度上限(l>thlh)的波谷(模符号为负)，用“.”表示类型为“2”的模极大曲线中的模极大小于模值上限大于模值下限(thmh>mo>thml)和模长度小于长度上限但大于长度下限(thlh>l>thll)的波谷。用“^”表示主波轮廓峰顶(当主波峰有多个局部接近极大值时，)。则脉象边沿和波峰谷特征点可进一步顺序表示如下：[1212122221112212112121121111111212]Based on the sequence of reserved modulus max curves above: use "+" to indicate the rising edge of the modulus max curve whose type is "1" (modulo max is positive), use "-" to indicate Falling edge in a modulo-maximum curve of type "1" (modulomaximum is negative). Use "|" to indicate the peak (the modulus symbol is positive) in the modulus maximal curve of type "2" where the modulus is greater than the upper limit of the modulus (mo>thmh) and the modulus length is greater than the upper limit of the length (l>thlh), and "*" indicates that the modulus maximum in the modulus max curve of type "2" is less than the modulus upper limit than the modulus lower limit (thmh>mo>thml) and the modulus length is less than the length upper limit but greater than the length lower limit (thlh>l>thll ) peak. Use "_" to indicate that the modulus maximum in the modulus maximal curve of type "2" is greater than the modulus value (mo>thmh) and the modulus length is greater than the trough of the upper limit (l>thlh) (the modulus sign is negative), use "." indicates that in the modulus max curve of type "2", the modulus max is smaller than the modulus upper limit and greater than the modulus lower limit (thmh>mo>thml) and the modulus length is less than the length upper limit but greater than the length lower limit (thlh>l>thll ) trough. Use "^" to indicate the peak of the main wave profile (when there are multiple local peaks close to the maximum value of the main wave). Then the edge of the pulse condition and the peak-to-valley feature points can be further expressed in sequence as follows:
-_-._+|***^|-.*-.*-_.._+*|--._.+*|-*-_.-.-----*-._+|-_-._+|***^|-.*-.*-_.._+*|--._.+*|-*-_.-.-----*-._ +|
由此我们可以识别出脉象各波的特征点，方法如下：周期脉象主波的起点对应第一个“|”，脉象主波顶点对应“^”，脉象主波的终点对应最后一个“|”。对照脉象小波包分解的基本轮廓近似位置，首先排除脉象尾波或抖动影响，并进一步根据小波模极大特征点，依次确定脉象其它波。通过确定脉象波波谷起点(大“_”或小“.”)、脉象波上升沿“+”、脉象波顶点(大“|”或小“*”)、脉象波下降沿“-”、脉象波波谷终点(大“_”或小“.”)，可以依次确定脉象主波、重搏前波(潮波)和重搏波位置。在主波附近1/3宽度外，从第一个“_”到“^”对应主波的上升沿上，顺序“+”间出现的“.”表示上升沿上有小不平，如果顺序连续“+”间出现“_”则表示上升沿有大不平。在主波附近1/3宽度内，顺序“+_+”间出现的“.”表示主波峰上有小不平，如果顺序顺序“+_+”间出现“_”则表示主波峰上有大不平，如果出现“+-”或“-+”则表示主波峰上有切迹。如果主波最大值“^”附近顺序出现“|^|”，如果“|”距离较远高度相近，则表示主波波峰为平顶。如果主波最大值“^”附近，相近高度范围内没有“|”或者“*”，则表明主波峰顶为尖相。从“^”到最后一个“*”为主波的下降沿，在主波的下降沿上顺序出现“_+|-_”表示存在明显重搏前波或重搏波，顺序出现“.-|-_”表示存在紧张(由”|”表示)的重搏前波，顺序出现“.-*-_”表示存在较小紧张(由”*”表示)的重搏前波，可进一步依据“*”或“|”点的平滑性和重搏前波相对宽度，可以确定波形是否为宽大弦型。确定是否有重搏波的方法是在主波附近1/3宽度范围外，从右边寻找到的第一个的脉波，则第第一个“.”或“_”对应重搏波起点(即降中峡)，顺序出现“-”前的“*”或“|”对应重搏波顶点，顺序对应“-”后的的第一个“.”或“_”为应重搏波终点。如果在主波附近1/3宽度范围外，不能从右边寻找脉波，则表示没有重搏波。确定是否有重搏前波的方法是在存在重搏波前提下，从主波下降沿右边寻找第二个脉波，则第第一个“.”或“_”对应重搏前波起点，顺序出现“-”前的“*”或“|”对应重搏前波顶点，顺序对应“-”后的的第一个“.”或“_”为应重搏前波终点。如果在主波附近1/3宽度范围外，不能从右边寻找脉波，则表示没有重搏前波。通过对照前面小波包轮廓分解可以进行波形确认，如果脉象主波顶点“^”和小波包分解轮廓顶点”～”基本重合，则表明脉象主波宽大涩性，否则表明脉象主波窄小顺应性好，从而可以区分重搏前波和主波顶峰接近重合时的涩性和滑性。From this we can identify the characteristic points of each wave of the pulse condition, the method is as follows: the starting point of the main wave of the periodic pulse condition corresponds to the first "|", the apex of the main wave of the pulse condition corresponds to "^", and the end point of the main wave of the pulse condition corresponds to the last "|" . Comparing with the approximate position of the basic outline of the wavelet packet decomposition of the pulse condition, the influence of the coda wave or jitter of the pulse condition is excluded first, and further according to the maximum characteristic point of the wavelet modulus, other waves of the pulse condition are determined in turn. By determining the starting point of the pulse wave trough (big "_" or small "."), the rising edge of the pulse wave "+", the apex of the pulse wave (big "|" or small "*"), the falling edge of the pulse wave "-", the pulse wave The end point of the wave trough (big "_" or small "."), can sequentially determine the position of the main pulse wave, pre-dicrotic wave (tidal wave) and dicrotic wave. Outside the 1/3 width of the main wave, from the first "_" to "^" corresponding to the rising edge of the main wave, the "." appearing between the sequence "+" indicates that there is a small unevenness on the rising edge. If the sequence is continuous A "_" between "+" means that there is a big unevenness on the rising edge. Within 1/3 of the width of the main wave, the "." appearing between the sequence "+_+" indicates that there is a small unevenness on the main wave peak. If the "_" appears between the sequence "+_+", it indicates that there is a large unevenness on the main wave peak. Uneven, if "+-" or "-+" appears, it means that there is a notch on the main peak. If "|^|" appears sequentially near the maximum value "^" of the main wave, and if "|" is far away and close in height, it means that the peak of the main wave is flat-topped. If there is no "|" or "*" in the vicinity of the main wave maximum value "^", it indicates that the main wave peak is a sharp phase. From "^" to the last "*" on the falling edge of the main wave, "_+|-_" appears sequentially on the falling edge of the main wave, indicating that there is an obvious dicrotic front wave or dicrotic wave, and the sequential appearance of ".- |-_" indicates that there is a tense (indicated by "|") dicrotic front wave, and the sequence of ".-*-_" indicates that there is a less tense (indicated by "*") dicrotic front wave, which can be further based on The smoothness of the "*" or "|" point and the relative width of the dicrotic front wave can determine whether the waveform is a broad chord type. The method to determine whether there is a dicrotic wave is to find the first pulse wave from the right outside the 1/3 width range near the main wave, then the first "." or "_" corresponds to the starting point of the dicrotic wave ( That is, Jiangzhongxia), the "*" or "|" before the "-" in sequence corresponds to the peak of the dicrotic wave, and the first "." or "_" after the sequence corresponding to the "-" is the end point of the dicrotic wave . If the pulse wave cannot be found from the right outside the 1/3 width range near the main wave, it means that there is no dicrotic wave. The method to determine whether there is a dicrotic wave is to look for the second pulse wave from the right side of the falling edge of the main wave under the premise of the dicrotic wave, then the first "." or "_" corresponds to the starting point of the dicrotic wave, The "*" or "|" before the "-" in sequence corresponds to the apex of the dicrotic anterior wave, and the first "." or "_" after the sequence corresponding to the "-" is the end point of the dicrotic anterior wave. If the pulse wave cannot be found from the right outside the 1/3 width range near the main wave, it means that there is no dicrotic front wave. Waveform confirmation can be carried out by comparing the wavelet packet contour decomposition above. If the pulse condition main wave apex "^" basically coincides with the wavelet packet decomposition contour apex "~", it indicates that the pulse condition main wave is broad and astringent; otherwise, it indicates that the pulse condition main wave is narrow and small. OK, so we can distinguish the astringency and slipperiness when the dicrotic prewave and main wave peaks are close to coincident.
根据上述原则脉象特征点符号分解结果如下所示，潮波和重搏前波出现“+”表示波形明显，尾部由于有部分上升沿因此也出现“+”。According to the above principles, the symbol decomposition results of pulse characteristic points are shown below. The appearance of "+" in tidal wave and dicrotic front wave indicates that the waveform is obvious, and "+" also appears in the tail because there is a part of the rising edge.
-._+|***^|-.*-.*-_.._ _+*|--._ .+*|-*-_ -*-._+|-._+|***^|-.*-.*-_.._ _+*|--._ .+*|-*-_ -*-._+|
主 潮 重 尾main tide tide heavy tail
大+ 大+ 大+Big+ Big+ Big+ Big+
进一步可以依次得到脉象主波起点(0点)、主波顶点(1点)、重搏前波起点(2点)、重搏前波顶点(3点)、降中峡点(4点)、重搏波顶点(5点)、重搏波终点(6点)、主波终点(7点)：Further, the starting point of main pulse wave (0 point), main wave peak (1 point), dicrotic front wave starting point (2 points), dicrotic front wave peak (3 points), Jiangzhongxia point (4 points), Dicrotic wave apex (5 points), dicrotic wave end point (6 points), main wave end point (7 points):
[0.1313 0.2469 0.5031 0.5781 0.6844 0.7531 0.8719 1.2156][0.1313 0.2469 0.5031 0.5781 0.6844 0.7531 0.8719 1.2156]
在确定上面脉象关键特征点后，我们可以很方便的计算下面脉象特征参数(如图7所示)。为准确计算脉象特征参数，可以根据周期脉象主波的起止点，进一步进行准确的基线偏移调整，脉象特征参数可以在此基准上准确计算得到。脉图的主要特征参数包括：After determining the key feature points of the pulse condition above, we can easily calculate the characteristic parameters of the pulse condition below (as shown in Figure 7). In order to accurately calculate the characteristic parameters of the pulse condition, further accurate baseline offset adjustment can be performed according to the start and end points of the main wave of the periodic pulse condition, and the characteristic parameters of the pulse condition can be accurately calculated on this basis. The main characteristic parameters of the pulse map include:
t1：为脉搏波图起始点主波峰顶之间的时值，t1对应于左心室的快速射血期；t1 : is the time value between the peaks of the main wave at the starting point of the pulse wave diagram, and t1 corresponds to the rapid ejection period of the left ventricle;
t2：为脉搏波图起始点到重搏前波起点之间的时值；t2 : the time value between the starting point of the pulse wave pattern and the starting point of the dicrotic pre-wave;
t3：为脉搏波图起始点到重搏前波峰顶之间的时值；t3 : the time value between the starting point of the pulse wave diagram and the peak of the pre-dicrotic wave;
t4：为脉搏波图起始点到降中峡之间的时值，t4对应于左心室的收缩期，此后时间对应于左心室的舒张期；t4 : the time value between the starting point of the pulse wave diagram and the descending gorge, t4 corresponds to the systolic period of the left ventricle, and the time thereafter corresponds to the diastolic period of the left ventricle;
t5：为脉搏波图起始点到重搏波峰顶之间的时值；t5 : the time value between the starting point of the pulse wave pattern and the peak of the dicrotic wave;
t6：为脉搏波图起始点到重搏波峰结束点之间的时值；t6 : the time value between the start point of the pulse wave pattern and the end point of the dicrotic peak;
T：为脉搏波图的起始点到终止点的时值。T对应于左心室的一个心动周期，对应于脉搏，亦即一个脉动周期；T: the time value from the start point to the end point of the pulse wave diagram. T corresponds to a cardiac cycle of the left ventricle, corresponding to the pulse, that is, a pulsation cycle;
h1：主波幅度，为主波峰顶到脉搏波图基线的高度(基线与时间轴平行)。主要反映左心室的射血功能和大动脉的射血功能和大动脉的顺应性，即左心室的射血功能强，大动脉顺应性好的状态下则h1高大；h1: main wave amplitude, the height from the main wave peak to the baseline of the pulse wave diagram (the baseline is parallel to the time axis). It mainly reflects the ejection function of the left ventricle, the ejection function of the aorta and the compliance of the aorta, that is, the ejection function of the left ventricle is strong, and the h1 is high when the compliance of the aorta is good;
h1e：主波峡幅度，是主波与重搏前波之间的一个低谷的幅度。其主要的生理意义与h2e一致，脉图分析时往往可以略去；h1e: the amplitude of the main wave gorge, which is the amplitude of a trough between the main wave and the front wave of the dicrosis. Its main physiological meaning is consistent with that of h2e, and it can often be omitted in the pulse diagram analysis;
h2：重搏前波幅度，为重搏前波峰顶到脉搏波图基线的高度。主要反映动脉血管收缩或者硬化导致张力升高，或外周阻力增高时，均可导致h2的幅度增高。重搏前波高度的抬高一般伴随其时相的提高，反映了动脉血管高张力、高阻力状态时，脉搏反射波传导速度的增快；h2: Pre-dicrotic wave amplitude, which is the height from the peak of the pre-dicrotic wave to the baseline of the pulse wave diagram. It mainly reflects that arterial vasoconstriction or hardening lead to increased tension, or increased peripheral resistance, which can lead to an increase in the amplitude of h2. The elevation of dicrotic anterior wave height is generally accompanied by the increase of its phase, which reflects the acceleration of the pulse reflection wave transmission velocity in the state of high tension and high resistance of arteries;
h2e：降中峡幅度，为降中峡谷底到脉搏波图基线的高度。降中峡高度主要反映动脉血管外周阻力大小，外周阻力增加时，表现为h2e增加；h2e: The amplitude of the descending canyon, which is the height from the bottom of the descending canyon to the baseline of the pulse wave map. The height of the descending isthmus mainly reflects the peripheral resistance of the arteries, and when the peripheral resistance increases, it is manifested as an increase in h2e;
h3：重搏波幅度，为重搏波峰顶到降中峡谷底所做的基线平行线之间的高度。重搏波幅度主要反映大动脉的弹性(顺应性)情况。当大动脉顺应性降低时，h3减少，或者为0(重搏波峰顶与降中峡谷底同有一水平)，甚至为负(重搏波峰顶重搏低于降中峡谷底水平)；h3: Dicrotic wave amplitude, the height between the baseline parallel lines from the peak of the dicrotic wave to the bottom of the descending canyon. Dicrotic wave amplitude mainly reflects the elasticity (compliance) of the aorta. When the compliance of the great artery decreases, h3 decreases, or it is 0 (the peak of the dicrotic wave is at the same level as the bottom of the descending canyon), or even negative (the peak of the dicrotic wave is lower than the level of the bottom of the descending canyon);
hw1：相对主波上1/3处宽度的主波高度；hw1: the height of the main wave relative to the width at 1/3 of the main wave;
hw2：相对重搏前波始终点的重搏前波高度；hw2: relative to the height of the dicrotic front wave at the point of the dicrotic front wave;
hw3：相对重搏波始终点的重搏波高度；hw3: the height of the dicrotic wave relative to the point of the dicrotic wave;
w1：主波上距顶部1/3高度处的宽度，相当于动脉内高压力水平维持的时间；w1: the width of the main wave at the height of 1/3 from the top, which is equivalent to the time when the high pressure level in the artery is maintained;
w2：为重搏前波的宽度；w2: the width of the dicrotic front wave;
w3：为重搏波的宽度；w3: the width of the dicrotic wave;
步骤102，对小波模极大变换后的脉象信号进行周期分割，其具体的实施过程如下：Step 102, the pulse condition signal after the wavelet modulus maximum transformation is carried out periodical segmentation, and its specific implementation process is as follows:
作为本发明的一个实施例，脉象主波起点附近具有最大上升斜率，因此高斯一阶导数小波变换地模极大最大，需要进行脉象周期分割可以将阀值设置为thmh＝0.8momax,thml＝0.8lmax。保留模大于模阀值上限thmh和模长度大于长度阀值上限thlh的模极大曲线，即可完成脉象周期分割。As an embodiment of the present invention, there is a maximum rising slope near the starting point of the main wave of the pulse condition, so the Gaussian first-order derivative wavelet transform has a maximum modulus, and it is necessary to divide the pulse period and the threshold value can be set to thmh =0.8momax ,thml = 0.8lmax . Keeping the maximum modulus curve whose modulus is greater than the upper limit of the modulus threshold value thmh and whose modulus length is greater than the upper limit of the length threshold value thlh can complete the division of the pulse cycle.
作为本发明的另一个实施例，为了进一步准确确定脉象主波起点，可在脉象边沿进行脉象周期分割。由于脉象主波起点具有最大谷底模极大，因此高斯二阶导数小波变换的模极大最大，需要进行脉象周期分割。优选地，在进行脉象周期分割时可以将阀值设置为thmh＝0.8momax,thml＝0.8lmax，然后保留模大于模阀值上限thmh并且模长度大于长度阀值上限thlh的模极大曲线。As another embodiment of the present invention, in order to further accurately determine the starting point of the main wave of the pulse condition, the pulse period can be divided at the edge of the pulse condition. Since the starting point of the main pulse wave has the maximum modulus at the bottom, the Gaussian second-order derivative wavelet transform has the largest modulus, and pulse cycle segmentation is required. Preferably, the threshold can be set to thmh = 0.8momax , thml = 0.8lmax when performing pulse cycle segmentation, and then retain the models whose modulus is greater than the upper limit of the modulus threshold thmh and whose modulus length is greater than the upper limit of the length threshold thlh modulus curve.
步骤103，对变换后的脉象信号的特征点进行定位，具体的实施过程如下：Step 103, locate the characteristic points of the transformed pulse signal, the specific implementation process is as follows:
为了对脉形不拘的脉象信号进行特征点定位，我们采用连续小波变换进行由粗到精的跟踪，通过将由粗到精的模极大连接成模极大曲线，来完成脉象信号的特征点定位。在实施例中，在进行奇异点定位时，采用二进小波变换进行由粗到精的策略跟踪各尺度下的小波变换极大值，方法是先从j最粗尺度的一级开始，找到这一尺度上属于信号的小波变换模极大值，然后逐步减小j尺度值，每次以高一级已找到的极值点位置为先验知识，寻找其在本级的对应极值点，直到逐级搜索到最细尺度为止。In order to locate the feature points of the pulse signal with an unconstrained pulse shape, we use the continuous wavelet transform to track from coarse to fine, and complete the feature point location of the pulse signal by connecting the modulus maxima from coarse to fine into a modulus maximal curve . In the embodiment, when performing singular point location, the binary wavelet transform is used to track the maximum value of the wavelet transform at each scale from coarse to fine. The maximum value of the wavelet transform modulus belonging to the signal on the first scale, and then gradually reduce the j-scale value, and each time use the position of the extreme point found at the higher level as the prior knowledge to find its corresponding extreme point at the current level, Until the smallest scale is searched step by step.
作为本发明的另一个实施例，由于在上面特征点定位的过程中，由粗到精的跟踪受到各种干扰因素影响，特别是对于脉象信号而言由于具有一定分形性，此时跟踪过程往往需要经验确定。因此，本发明首先通过对脉象周期的确定，可以获得连续脉象信号的频率和节律特征。即当小波取为高斯函数的一阶导数时，小波模极大表示脉象信号各组成波的上升沿和下降沿的急剧变化点；当小波取为高斯函数的二阶导数时，小波模极大表示脉象信号各组成波的谷底和波峰的急剧变化点。As another embodiment of the present invention, in the process of locating the above feature points, the tracking from coarse to fine is affected by various interference factors, especially for the pulse signal because of its fractal nature, the tracking process is often Needs to be determined empirically. Therefore, the present invention can first obtain the frequency and rhythm characteristics of the continuous pulse signal by determining the pulse cycle. That is, when the wavelet is taken as the first-order derivative of the Gaussian function, the maximum wavelet modulus indicates the sharp change point of the rising edge and the falling edge of each component wave of the pulse signal; when the wavelet is taken as the second-order derivative of the Gaussian function, the wavelet modulus is extremely large Indicates the sharp change point of the valley and peak of each constituent wave of the pulse signal.
然后，在小波包分解确定脉象信号基本轮廓基础上，通过排列脉象波上升沿与下降沿特征点、脉象波的谷底和波峰特征点，可以确定不明显的重搏前波、重搏波，以及在主波上升沿和下降沿上出现的不规则小波，从而获得对脉象波形的整体认识，进一步准确获得脉象信号的多尺度特征和各种时域特征。Then, on the basis of wavelet packet decomposition to determine the basic outline of the pulse signal, by arranging the rising and falling edge feature points of the pulse wave, the valley and peak feature points of the pulse wave, it is possible to determine the non-obvious dicrotic pre-wave, dicrotic wave, and Irregular wavelets appear on the rising and falling edges of the main wave, so as to obtain the overall understanding of the pulse waveform, and further accurately obtain the multi-scale characteristics and various time-domain characteristics of the pulse signal.
其中，多尺度特征补充了一些时域难于表达的特征。例如通过衰减度参数可以区别噪声性奇点和信号特征点；通过脉波峰顶特征点的模极大值可以研究脉象的紧张程度；通过平滑参数可以了解弧大弦长程度；通过分形参数可以获知脉象信号分形形态，如当脉象由滑至涩时分形程度增大。从小波模极大沿趋于v点的极大曲线s上取不同点，由模极大参数公式可以计算出参数α及σ，以及分形维数D(α)，这些参数从不同角度描述了脉形的多尺度特征，建立了不同尺度下特征点的模极大联系。Among them, multi-scale features supplement some features that are difficult to express in the time domain. For example, the noise singularity and the signal characteristic point can be distinguished through the attenuation parameter; the tension of the pulse can be studied through the modulus maximum value of the characteristic point of the pulse wave peak; the degree of arc chord length can be understood through the smoothing parameter; The fractal shape of the pulse signal, for example, the degree of fractal increases when the pulse changes from slippery to astringent. Different points are taken from the maximal curve s that tends to point v from the wavelet modulus maxima, and the parameters α and σ, as well as the fractal dimension D(α), can be calculated from the modulus maxima parameter formula. These parameters describe the The multi-scale feature of the vein shape establishes the modulus-maximum connection of the feature points at different scales.
下面对一两峰脉象提取脉形参数，此两峰脉象重搏前波基本消失。图8消扰脉象信号采用高斯一阶导数进行的小波变换。图9采用高斯二阶导数进行的小波变换。由此可得:脉象主波、潮波和重搏波的起点坐标、顶点坐标和终点时间点坐标。脉象主波、重搏前波(潮波)和重搏波的绝对高度和相对高度(主波相对W31)，以及绝对宽度和相对宽度(主波指W31)。脉象主波、重搏前波(潮波)和重搏波顶点的模极大、衰减参数、平滑参数和分维。Next, the pulse shape parameters are extracted for a two-peak pulse condition, and the dicrotic anterior wave of the two-peak pulse condition basically disappears. Fig. 8 is the wavelet transform performed by the Gaussian first-order derivative for the pulse signal elimination. Fig. 9 Wavelet transform using Gaussian second derivative. Thus it can be obtained: the starting point coordinates, apex coordinates and end time point coordinates of pulse condition main wave, tidal wave and dicrotic wave. Absolute height and relative height (main wave relative to W31) and absolute width and relative width (main wave refers to W31) of main wave of pulse condition, dicrotic front wave (tidal wave) and dicrotic wave. The modulus of main wave of pulse condition, dicrotic front wave (tidal wave) and dicrotic wave apex, attenuation parameter, smoothing parameter and fractal dimension.
参阅图10所示，为本发明一种脉象信号的提取装置的结构示意图，所述的脉象信号的提取装置根据上面的描述可以包括：Referring to shown in Fig. 10, it is the structural representation of the extracting device of a kind of pulse condition signal of the present invention, the extracting device of described pulse condition signal can comprise according to above description:
变换单元1001，能够对采集的脉象信号进行小波模极大变换，其具体的功能包括：Transformation unit 1001 can carry out wavelet modulus maximum transformation to the pulse signal of collection, and its concrete function comprises:
首先，采集脉象信号，并得到脉象离散数据。Firstly, collect the pulse signal and obtain discrete pulse data.
作为本发明的一个实施例，可以使用仪器脉搏图形记录仪采集脉象信号，然后通过对脉象信号的标准脉图进行数据化得到脉象数据。较佳地，将脉象信号转化为脉象数据是通过数字笔将脉象图采集到计算机后，以x轴象素间隔dx(mm)和y轴象素间隔dy(mm)为单位进行数据化，离散化后的脉象数据形成脉象标准数据库。优选地，形成脉象标准数据库可以将离散化后的脉象数据通过转化为以实际时间单位表示图形，只需将x轴和y轴分别以实际采样率(25mm/dx)/sec和(5mm/dy)/mv转化即可。As an embodiment of the present invention, the pulse pattern recorder can be used to collect the pulse condition signal, and then the pulse condition data can be obtained by digitizing the standard pulse pattern of the pulse condition signal. Preferably, the pulse signal is converted into pulse data after the pulse image is collected into the computer by a digital pen, and then digitized in units of x-axis pixel interval dx (mm) and y-axis pixel interval dy (mm), discrete The transformed pulse data forms a pulse standard database. Preferably, forming the pulse condition standard database can convert the discretized pulse condition data into a graph expressed in actual time units, and only need to convert the x-axis and y-axis to the actual sampling rate (25mm/dx)/sec and (5mm/dy )/mv conversion.
其次，对脉象离散数据进行高斯一阶导数小波模极大变换。Secondly, Gaussian first-order derivative wavelet modulo-maximum transformation is performed on the pulse condition discrete data.
为得到脉象离散数据对应的脉象边沿的高斯一阶导数小波变换模极大曲线，需要对噪声产生的模极大曲线进行过滤，方法如下：In order to obtain the Gaussian first-order derivative wavelet transform modulus max curve of the pulse image edge corresponding to the pulse image discrete data, it is necessary to filter the modulus max curve generated by the noise, the method is as follows:
1)计算所有模极大曲线起点的模均值vm＝E(moj)(j＝1,…,n)和方差dm2＝D(moj)(j＝1,…,n),模阀值上限thmh＝vm+dm,模阀值下限thml＝vm-dm。1) Calculate the modulus mean vm =E(moj )(j=1,...,n) and variance dm2 =D(moj )(j=1,...,n) of the starting points of all modulus maximal curves, The upper limit of the modulus threshold thmh =vm +dm , the lower limit of the modulus threshold thml =vm -dm .
2)计算所有模极大曲线长度的均值vl＝E(lmj)(j＝1,…,n)和方差dl2＝D(lmj)(j＝1,…,n)，长度阀值上限thlh＝vl+dl,长度阀值下限thll＝vl-dl。2) Calculate the mean vl =E(lmj )(j=1,...,n) and variance dl2 =D(lmj )(j=1,...,n) of the lengths of all modulus maximal curves, the length The upper limit of the threshold value thlh =vl +dl , the lower limit of the length threshold thll =vl -dl .
3)计算脉象主波1/3宽度区间Xw＝[twb,twe],即距离主波顶部1/3主波高度的宽度区间。其中主波波峰位置为主波极大值位置，令w＝(tp-tb)/3，则主波宽度区间Xw近似为主波顶点为中心的w对称宽度区间。3) Calculate the 1/3 width interval of the main wave of the pulse condition Xw =[twb , twe ], that is, the width interval of 1/3 the height of the main wave from the top of the main wave. Wherein the peak position of the main wave is the position of the maximum value of the main wave, if w=(tp -tb )/3, then the width interval Xw of the main wave is approximately the w symmetric width interval centered on the apex of the main wave.
4)对主波1/3宽度附近内的模极大曲线，为研究峰值附近的小隆起不平，保留模大于模阀值下限thml和模长度大于长度阀值下限thll的模极大曲线。4) For the maximum modulus curve in the vicinity of the 1/3 width of the main wave, in order to study the small bumps near the peak, the maximum modulus curve with a mode greater than the lower limit of the mode threshold value thml and a mode length greater than the lower limit of the length threshold value thll is reserved .
5)对主波1/3宽度区间外的模极大曲线，保留模大小大于模阀值上限thmh和模长度大于长度阀值上限thlh的模极大曲线。5) For the maximal modulus curves outside the 1/3 width interval of the main wave, retain the maximal modulus curves whose modulus size is greater than the upper limit of the modulus threshold thmh and whose modulus length is greater than the upper limit of the length threshold thlh .
6)进一步滤除Lipschitz指数α显著小(α<＝－100)和平滑因子σ显著大(σ>＝0.25)的奇异性模极大曲线。6) Further filter out curves with a significantly smaller Lipschitz exponent α (α<=-100) and a significantly larger smoothing factor σ (σ>=0.25).
在本发明的实施例中，可以用小波描述信号的局部奇异性，奇异点小波模极大衰减性可以通过Lipschitz指数α描述，奇异点磨光平滑性可以通过平滑因子σ描述，奇异点分形性可以通过分形维数D(α)描述。In the embodiment of the present invention, the local singularity of the signal can be described by the wavelet, the extreme attenuation of the wavelet mode of the singular point can be described by the Lipschitz exponent α, the smoothness of the singular point can be described by the smoothing factor σ, and the fractal property of the singular point It can be described by fractal dimension D(α).
其中，小波检测信号奇异点衰减特性用Lipschitz指数来描述：Among them, the attenuation characteristics of the singular point of the wavelet detection signal are described by Lipschitz exponent:
设n是一正整数，n＜α≤n+1，f(t)在t0点的Lipschitz指数α充要条件是存在两个常数A和h0(A＞0，h0＞0)，且有一n次多项式Pn(h)使得对h≤h0有：Suppose n is a positive integer, n<α≤n+1, the Lipschitz exponent α of f(t) at point t0 is a necessary and sufficient condition that there are two constants A and h0 (A>0, h0 >0), And there is an n-degree polynomial Pn (h) such that for h≤h0 :
|f(t0+h)-Pn(t0+h)|≤A|h|α，则称f(t)在点t0处为Lipschitzα。实际上多项式Pn(h)是f(t)在t0点的taylor展开式的前n项：|f(t0 +h)-Pn (t0 +h)|≤A|h|α , then f(t) is called Lipschitzα at point t0 . In fact, the polynomial Pn (h) is the first n terms of the Taylor expansion of f(t) at point t0 :
f(t)＝f(t0)+a1h+a2h2+--+anhn+O(hn+1)＝Pn(t)+O(hn+1)。f(t)=f(t0 )+a1 h+a2 h2 +--+an hn +O(hn+1 )=Pn (t)+O(hn+1 ).
由此可见：如果f(t)n次可微，但n阶导数不连续，因此n+1次不可微，则n<α≤n+1，相对于taylor公式中对信号分解误差的n+1次描述，Lipschitzα正则性通过使用非整数的指数改进了这个误差上界。如果上式对所有t0和t0+h在区间[t1,t2]内均成立，则称f(t)在此区间为一致Lipschitzα。It can be seen that if f(t) is differentiable for n times, but the nth order derivative is discontinuous, so n+1 times are not differentiable, then n<α≤n+1, compared to the n+ of the signal decomposition error in Taylor's formula 1 description, Lipschitz α regularity improves this error upper bound by using non-integer exponents. If the above formula is true for all t0 and t0 +h in the interval [t1 ,t2 ], then f(t) is said to be consistent Lipschitzα in this interval.
若一个函数在一点连续且一次可微，但导数不连续，那么它在这点就是Lipschitz指数为1。若f(t)在点t0处不连续但有界，即间断的情形,则Lipschitz指数为0。Lipschitz指数越大，f(t)越光滑。例如，斜坡函数的Lipschitz指数为1，阶跃函数的Lipschitz指数为0，δ函数的Lipschitz指数为-1，白噪的Lipschitz指数为-0.5-ε(ε>0)。If a function is continuous and once differentiable at a point, but its derivative is discontinuous, then it has a Lipschitz exponent of 1 at that point. If f(t) is discontinuous but bounded at point t0 , that is, discontinuous, the Lipschitz exponent is 0. The larger the Lipschitz exponent, the smoother f(t) is. For example, the Lipschitz exponent of the slope function is 1, the Lipschitz exponent of the step function is 0, the Lipschitz exponent of the delta function is -1, and the Lipschitz exponent of the white noise is -0.5-ε (ε>0).
当f(t)完全是正则函数时，同奇异点一样其小波变换也可以有一系列取模极大的点趋于横坐标点v，因此仅沿着尺度搜索小波模极大对于奇异性检测是不够的，需要从模极大幅值的衰减性计算Lipschitz正则性。When f(t) is completely a regular function, its wavelet transform, like the singular point, can also have a series of points with the maximum modulus tending to the abscissa point v, so only searching for the maximum modulus of the wavelet along the scale is useful for singularity detection Not enough, Lipschitz regularity needs to be calculated from the attenuation of the modulo extremely large value.
当v点是孤立奇异点时，对s<s0，收敛于v的所有模极大包含在如下的锥中：|u-v|≤CsWhen point v is an isolated singular point, for s<s0 , all modulus maxima converging on v are contained in the following cone: |uv|≤Cs
式中C为ψ的紧支集[-C,C]，在尺度－时间平面上点v的影响锥是所有点(u,v)的集合。这意味着f在v的邻域内没有快速振荡，|Wf(u,s)|在v的邻域里的衰减性由含于锥中的模极大衰减性控制。f在v的邻域内是一致Lipschitzα的，并且仅当存在K>0，使得锥中的模极大点(u,v)满足：In the formula, C is the compact support set [-C, C] of ψ, and the influence cone of point v on the scale-time plane is the set of all points (u, v). This means that there is no fast oscillation of f in the neighborhood of v, and the decay of |Wf (u,s)| in the neighborhood of v is controlled by the decay of the modulus maxima contained in the cone. f is uniform Lipschitzα in the neighborhood of v, and only if there exists K>0 such that the modulus maxima (u,v) in the cone satisfy:
上式是Lipschitzα对未磨光的孤立奇异点的小波模极大变化影响，在v点的Lipschitz正则性就是log2|Wf(u,s)|作为log2s的函数沿着收敛于v点的极大曲线的最大斜率。The above formula is the influence of Lipschitzα on the wavelet mode maximum change of the unpolished isolated singular point. The Lipschitz regularity at point v is log2 |Wf (u,s)| as a function of log2 s and converges to v The maximum slope of the maximal curve for a point.
另外，即使无限次连续可微的信号，也可以有很大的变化，如图像中阴影处的边界上，图像的亮度有很快的变化，但由于衍射作用，它却是不连续的。将这些跳跃处光滑化就是用高斯核将之模拟成一个扩散过程，其根方差就是表征边沿平滑特性的平滑因子σ，可以通过小波模极大的衰减性来测量。在峰变处v的邻域内，设f(t)＝f0*gσ(t)，gσ是方差为σ2的高斯函数：In addition, even infinitely continuous differentiable signals can have great changes, such as the image brightness changes rapidly on the boundary of shadows in the image, but it is discontinuous due to the effect of diffraction. To smooth these jumps is to use Gaussian kernel to simulate it as a diffusion process, and its root variance is the smoothing factor σ that characterizes the smoothness of the edge, which can be measured by the extreme attenuation of the wavelet mode. In the neighborhood of v where the peak changes, let f(t)=f0 *gσ (t), gσ is a Gaussian function with variance σ2 :
如果f0在v点的Lipschitzα，则它在v的邻域内是一致Lipschitzα的，当小波函数取为高斯函数的导数时ψ＝(-1)nθ(n)，其中θ(t)＝λexp(-t2/(2β2))，则f在v的邻域内是一致Lipschitzα的，并且仅当存在K>0，使得锥中的模极大点(u,v)满足：If f0 is Lipschitzα at point v, it is consistent Lipschitzα in the neighborhood of v, when the wavelet function is taken as the derivative of the Gaussian function, ψ=(-1)n θ(n) , where θ(t)=λexp (-t2 /(2β2 )), then f is consistent Lipschitzα in the neighborhood of v, and only if there is K>0, so that the modulus maximum point (u,v) in the cone satisfies:
上式对孤立奇异点考虑了磨光平滑性，最后一项表达了对小波模极大变化的影响。方差β2依赖于小波的选取，是预先给定的。对于较大的s>>σ/β，小波变化模按sα+1/2衰减，高斯平均对它的影响还感觉不到。当s≤σ/β，由于高斯平均的影响，f在v点的变化不非常强的关联于s。在这些细的尺度下，小波变换模以sn+1/2级衰减。The above formula considers the smoothness of the smoothness of the isolated singular point, and the last term expresses the influence on the extreme change of the wavelet mode. Variance β2 depends on the selection of wavelet and is given in advance. For a larger s>>σ/β, the wavelet change mode decays according to sα+1/2 , and the influence of Gaussian averaging on it cannot be felt. When s≤σ/β, the change of f at v is not very strongly related to s due to the Gaussian average. At these fine scales, the wavelet transform modulus decays by stepsn+1/2 .
另外，当信号具有非孤立奇异性时，信号离散化后的每一个采样点来自于一个小的时间区间，信号在该区间中含有无穷多个互异的奇异性，由于数值分辨率的有限性，Lipschitz指数的点态度量并不可行，因此这种奇异性分布必须利用多分形的自相似整体度量，以给出了不同Lipschitz正则性的奇异特性的整体分配，这种度量即奇异谱分形维数D(α)。In addition, when the signal has non-isolated singularities, each sampling point after the discretization of the signal comes from a small time interval, and the signal contains infinitely many different singularities in this interval, due to the limited numerical resolution , the point-attitude measure of the Lipschitz exponent is not feasible, so this singularity distribution must use a multifractal self-similar overall measure to give the overall distribution of the singular properties of different Lipschitz regularities, this measure is the singular spectrum fractal dimension Number D(α).
设Sα是使得信号f在其上的Lipschitz正则性为α的所有点的集合，f的奇异谱D(α)定义为Sα的分形维数，即Nα(s)～s-D(α)。式中Nα(s)表示Sα中点的个数，奇异谱给出了在任意尺度s下Lipschitzα奇异性的比例值。当自相似信号的谱D(α)是凸的，其计算过程如下：对每一尺度s计算出Wf(u,s)和模极大，随s的变化链接小波模极大。计算分解函数：然后，利用log2Z(s,q)作为log2s的函数的线性衰减性计算τ(q)，即log2Z(q,s)≈τ(q)log2s+C(q)。最后计算谱D(α)：对D(α)可采用四元组[αmin,αmax,αDmax,Dmax]描述，表示D(α)的支撑域为[αmin,αmax]，当α＝αDmax时D(α)取最大值Dmax(aDmax)。Let Sα be the set of all points on which the Lipschitz regularity of the signal f is α, and the singular spectrum D(α) of f is defined as the fractal dimension of Sα , that is, Nα (s)～s-D( a) . In the formula, Nα (s) represents the number of points in Sα , and the singularity spectrum gives the proportional value of the Lipschitz α singularity at any scale s. When the spectrum D(α) of the self-similar signal is convex, the calculation process is as follows: Wf (u, s) and modulus are calculated for each scale s, and the wavelet modulus is maximized with the change of s. Compute the factorization function: Then, τ(q) is calculated using the linear decay of log2 Z(s,q) as a function of log2 s, ie log2 Z(q,s)≈τ(q)log2 s+C(q). Finally calculate the spectrum D(α): D(α) can be described by the quadruple [αmin ,αmax ,αDmax ,Dmax ], indicating that the support domain of D(α) is [αmin ,αmax ], when α=αDmax D( α) takes the maximum value Dmax (aDmax ).
然后，对脉象离散数据进行高斯二阶导数小波模极大变换。Then, Gaussian second-order derivative wavelet modulo-maximum transformation is performed on the pulse condition discrete data.
为得到脉象离散数据对应的脉象局部极大极小点的高斯二阶导数小波变换模极大曲线，需要对噪声产生的模极大曲线进行过滤，方法如下：In order to obtain the Gaussian second-order derivative wavelet transform modulus max curve of the pulse image local maxima and minima points corresponding to the pulse image discrete data, it is necessary to filter the modulus max curve generated by the noise, the method is as follows:
1)计算所有模极大曲线起点的模均值vm＝E(moj)(j＝1,…,n)和方差dm2＝D(moj)(j＝1,…,n),模阀值上限thmh＝vm+dm,模阀值下限thml＝vm-dm。1) Calculate the modulus mean vm =E(moj )(j=1,...,n) and variance dm2 =D(moj )(j=1,...,n) of the starting points of all modulus maximal curves, The upper limit of the modulus threshold thmh =vm +dm , the lower limit of the modulus threshold thml =vm -dm .
2)计算所有模极大曲线长度的均值vl＝E(lmj)(j＝1,…,n)和方差dl2＝D(lmj)(j＝1,…,n),长度阀值上限thlh＝vl+dl,长度阀值下限thll＝vl-dl。2) Calculate the mean vl =E(lmj )(j=1,...,n) and variance dl2 =D(lmj )(j=1,...,n) of the lengths of all modulus maximal curves, the length The upper limit of the threshold value thlh =vl +dl , the lower limit of the length threshold thll =vl -dl .
3)计算脉象主波1/3宽度区间Xw＝[twb,twe],即距离主波顶部1/3主波高度的宽度区间。其中主波波峰位置为主波波峰极大值位置，令w＝(tp-tb)/3则主波宽度区间Xw近似为主波顶点为中心的w对称宽度区间。3) Calculate the 1/3 width interval of the main wave of the pulse condition Xw =[twb , twe ], that is, the width interval of 1/3 the height of the main wave from the top of the main wave. Wherein the peak position of the main wave is the position of the maximum value of the main wave peak, if w=(tp -tb )/3, the width interval Xw of the main wave is approximate to the w-symmetrical width interval centered on the apex of the main wave.
4)对主波宽度区间内的模极大曲线，保留模大小大于模阀值下限thml和模长度，大于长度阀值下限thll的模极大曲线。4) For the maximal modulus curves in the main wave width interval, retain the maximal modulus curves whose mode size is greater than the lower limit of the mode threshold thml and the mode length, and greater than the lower limit of the length threshold thll .
5)对主波宽度区间外的模极大曲线，保留模大小大于模阀值上限thmh和模长度大于长度阀值上限thlh的模极大曲线。5) For the maximal modulus curves outside the main wave width interval, retain the maximal modulus curves whose mode size is greater than the upper limit of the mode threshold value thmh and whose mode length is greater than the upper limit of the length threshold value thlh .
作为本发明的一个实施例，为了对噪声产生的模极大曲线过滤效果更好，进一步可以滤除Lipschitz指数α显著小(α<＝－100)和平滑因子σ显著大(σ>＝0.25)的奇异性模极大曲线。As an embodiment of the present invention, in order to have a better filtering effect on the maximum modulus curve generated by noise, it is further possible to filter out significantly small Lipschitz exponent α (α<=-100) and significantly large smoothing factor σ (σ>=0.25) The singularity modulus maximal curve of .
分割单元1002，与变换单元1001连接，能够对小波模极大变换后的脉象信号进行周期分割。The division unit 1002 is connected with the transformation unit 1001, and can perform periodic division on the pulse condition signal after wavelet modulo-maximum transformation.
作为本发明的一个实施例，脉象主波起点附近具有最大上升斜率，因此高斯一阶导数小波变换地模极大最大，需要进行脉象周期分割可以将阀值设置为thmh＝0.8momax,thml＝0.8lmax。保留模大于模阀值上限thmh和模长度大于长度阀值上限thlh的模极大曲线，即可完成脉象周期分割。As an embodiment of the present invention, there is a maximum rising slope near the starting point of the main wave of the pulse condition, so the Gaussian first-order derivative wavelet transform has a maximum modulus, and it is necessary to divide the pulse period and the threshold value can be set to thmh =0.8momax ,thml = 0.8lmax . Keeping the maximum modulus curve whose modulus is greater than the upper limit of the modulus threshold value thmh and whose modulus length is greater than the upper limit of the length threshold value thlh can complete the division of the pulse cycle.
作为本发明的另一个实施例，为了进一步准确确定脉象主波起点，可在脉象边沿进行脉象周期分割。由于脉象主波起点具有最大谷底模极大，因此高斯二阶导数小波变换的模极大最大，需要进行脉象周期分割。优选地，在进行脉象周期分割时可以将阀值设置为thmh＝0.8momax,thml＝0.8lmax，然后保留模大于模阀值上限thmh并且模长度大于长度阀值上限thlh的模极大曲线。As another embodiment of the present invention, in order to further accurately determine the starting point of the main wave of the pulse condition, the pulse period can be divided at the edge of the pulse condition. Since the starting point of the main pulse wave has the maximum modulus at the bottom, the Gaussian second-order derivative wavelet transform has the largest modulus, and pulse cycle segmentation is required. Preferably, the threshold can be set to thmh = 0.8momax , thml = 0.8lmax when performing pulse cycle segmentation, and then retain the models whose modulus is greater than the upper limit of the modulus threshold thmh and whose modulus length is greater than the upper limit of the length threshold thlh modulus curve.
定位单元1003，与分割单元1002连接，能够对变换后的脉象信号的特征点进行定位。The positioning unit 1003 is connected with the segmentation unit 1002 and can locate the characteristic points of the transformed pulse signal.
为了对脉形不拘的脉象信号进行特征点定位，我们采用连续小波变换进行由粗到精的跟踪，通过将由粗到精的模极大连接成模极大曲线，来完成脉象信号的特征点定位。在实施例中，在进行奇异点定位时，采用二进小波变换进行由粗到精的策略跟踪各尺度下的小波变换极大值，方法是先从j最粗尺度的一级开始，找到这一尺度上属于信号的小波变换模极大值，然后逐步减小j尺度值，每次以高一级已找到的极值点位置为先验知识，寻找其在本级的对应极值点，直到逐级搜索到最细尺度为止。In order to locate the feature points of the pulse signal with an unconstrained pulse shape, we use the continuous wavelet transform to track from coarse to fine, and complete the feature point location of the pulse signal by connecting the modulus maxima from coarse to fine into a modulus maximal curve . In the embodiment, when performing singular point location, the binary wavelet transform is used to track the maximum value of the wavelet transform at each scale from coarse to fine. The maximum value of the wavelet transform modulus belonging to the signal on the first scale, and then gradually reduce the j-scale value, and each time use the position of the extreme point found at the higher level as the prior knowledge to find its corresponding extreme point at the current level, Until the smallest scale is searched step by step.
作为本发明的另一个实施例，由于在上面特征点定位的过程中，由粗到精的跟踪受到各种干扰因素影响，特别是对于脉象信号而言由于具有一定分形性，此时跟踪过程往往需要经验确定。因此，本发明首先通过对脉象周期的确定，可以获得连续脉象信号的频率和节律特征。即当小波取为高斯函数的一阶导数时，小波模极大表示脉象信号各组成波的上升沿和下降沿的急剧变化点；当小波取为高斯函数的二阶导数时，小波模极大表示脉象信号各组成波的谷底和波峰的急剧变化点。As another embodiment of the present invention, in the process of locating the above feature points, the tracking from coarse to fine is affected by various interference factors, especially for the pulse signal because of its fractal nature, the tracking process is often Needs to be determined empirically. Therefore, the present invention can first obtain the frequency and rhythm characteristics of the continuous pulse signal by determining the pulse cycle. That is, when the wavelet is taken as the first-order derivative of the Gaussian function, the maximum wavelet modulus indicates the sharp change point of the rising edge and the falling edge of each component wave of the pulse signal; when the wavelet is taken as the second-order derivative of the Gaussian function, the wavelet modulus is extremely large Indicates the sharp change point of the valley and peak of each constituent wave of the pulse signal.
然后，在小波包分解确定脉象信号基本轮廓基础上，通过排列脉象波上升沿与下降沿特征点、脉象波的谷底和波峰特征点，可以确定不明显的重搏前波、重搏波，以及在主波上升沿和下降沿上出现的不规则小波，从而获得对脉象波形的整体认识，进一步准确获得脉象信号的多尺度特征和各种时域特征。Then, on the basis of wavelet packet decomposition to determine the basic outline of the pulse signal, by arranging the rising and falling edge feature points of the pulse wave, the valley and peak feature points of the pulse wave, it is possible to determine the non-obvious dicrotic pre-wave, dicrotic wave, and Irregular wavelets appear on the rising and falling edges of the main wave, so as to obtain the overall understanding of the pulse waveform, and further accurately obtain the multi-scale characteristics and various time-domain characteristics of the pulse signal.
其中，多尺度特征补充了一些时域难于表达的特征。例如通过衰减度参数可以区别噪声性奇点和信号特征点；通过脉波峰顶特征点的模极大值可以研究脉象的紧张程度；通过平滑参数可以了解弧大弦长程度；通过分形参数可以获知脉象信号分形形态，如当脉象由滑至涩时分形程度增大。从小波模极大沿趋于v点的极大曲线s上取不同点，由模极大参数公式可以计算出参数α及σ，以及分形维数D(α)，这些参数从不同角度描述了脉形的多尺度特征，建立了不同尺度下特征点的模极大联系。Among them, multi-scale features supplement some features that are difficult to express in the time domain. For example, the noise singularity and the signal characteristic point can be distinguished through the attenuation parameter; the tension of the pulse can be studied through the modulus maximum value of the characteristic point of the pulse wave peak; the degree of arc chord length can be understood through the smoothing parameter; The fractal shape of the pulse signal, for example, the degree of fractal increases when the pulse changes from slippery to astringent. Different points are taken from the maximal curve s that tends to point v from the wavelet modulus maxima, and the parameters α and σ, as well as the fractal dimension D(α), can be calculated from the modulus maxima parameter formula. These parameters describe the The multi-scale feature of the vein shape establishes the modulus-maximum connection of the feature points at different scales.
从上面的描述可以看出，本发明实现的一种脉象信号的提取方法和装置，创造性地提出了对采集的脉象信号进行小波模极大变换、周期分割、定位来完成准确处理特征不明显的时变弱信号；同时，当小波取为高斯函数的一阶导数时，小波模极大表示脉象信号各组成波的上升沿和下降沿的急剧变化点；当小波取为高斯函数的二阶导数时，小波模极大表示脉象信号各组成波的谷底和波峰的急剧变化点。通过对脉象周期的确定，可以获得连续脉象信号的频率和节律特征；而且，在小波包分解确定脉象信号基本轮廓基础上，通过排列脉象波上升沿与下降沿特征点、脉象波的谷底和波峰特征点，可以确定不明显的重搏前波、重搏波，以及在主波上升沿和下降沿上出现的不规则小波，从而获得对脉象波形的整体认识，进一步准确获得脉象信号的多尺度特征和各种时域特征；与此同时，多尺度特征补充了一些时域难于表达的特征。通过衰减度参数可以区别噪声性奇点和信号特征点；通过脉波峰顶特征点的模极大值可以研究脉象的紧张程度；通过平滑参数可以了解弧大弦长程度；通过分形参数可以获知脉象信号分形形态，如当脉象由滑至涩时分形程度增大；最后，整个所述的脉象信号的提取方法和装置简便、紧凑，易于实现。As can be seen from the above description, a method and device for extracting a pulse signal realized by the present invention creatively proposes to carry out wavelet modulus maximum transformation, periodic segmentation, and positioning to the collected pulse signal to complete accurate processing. At the same time, when the wavelet is taken as the first-order derivative of the Gaussian function, the wavelet modulus is extremely large, indicating the sharp change point of the rising edge and falling edge of each component wave of the pulse signal; when the wavelet is taken as the second-order derivative of the Gaussian function When the wavelet modulus is extremely large, it indicates the sharp change point of the bottom and peak of each component wave of the pulse signal. By determining the pulse period, the frequency and rhythm characteristics of the continuous pulse signal can be obtained; moreover, on the basis of determining the basic outline of the pulse signal by wavelet packet decomposition, by arranging the rising and falling edge feature points of the pulse wave, the valley and the peak of the pulse wave Feature points can determine the inconspicuous dicrotic pre-wave, dicrotic wave, and irregular wavelets appearing on the rising and falling edges of the main wave, so as to obtain an overall understanding of the pulse waveform and further accurately obtain the multi-scale pulse signal features and various temporal features; at the same time, multi-scale features supplement some features that are difficult to express in the temporal domain. The attenuation parameter can be used to distinguish the noise singularity and the signal feature point; the pulse condition can be studied through the modulus maximum value of the pulse wave peak feature point; the degree of arc chord length can be understood through the smoothing parameter; the pulse condition can be known through the fractal parameter The fractal form of the signal, for example, the degree of fractal increases when the pulse condition changes from slippery to astringent; finally, the entire pulse signal extraction method and device are simple, compact and easy to implement.
所属领域的普通技术人员应当理解：以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。Those of ordinary skill in the art should understand that: the above descriptions are only specific embodiments of the present invention, and are not intended to limit the present invention. Any modifications, equivalent replacements, and improvements made within the spirit and principles of the present invention etc., should be included within the protection scope of the present invention.
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201410163098.2ACN103932686A (en) | 2014-04-22 | 2014-04-22 | Method and device for extracting pulse condition signal |
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201410163098.2ACN103932686A (en) | 2014-04-22 | 2014-04-22 | Method and device for extracting pulse condition signal |
| Publication Number | Publication Date |
|---|---|
| CN103932686Atrue CN103932686A (en) | 2014-07-23 |
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201410163098.2APendingCN103932686A (en) | 2014-04-22 | 2014-04-22 | Method and device for extracting pulse condition signal |
| Country | Link |
|---|---|
| CN (1) | CN103932686A (en) |
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107736880A (en)* | 2017-10-24 | 2018-02-27 | 新绎健康科技有限公司 | A kind of pulse analysis method and system |
| CN107811619A (en)* | 2017-12-08 | 2018-03-20 | 西安科技大学 | Portable pulse-taking instrument and its analysis method |
| CN110595956A (en)* | 2019-08-09 | 2019-12-20 | 浙江工业大学 | A wear state mutation detection method based on fractal characteristics of abrasive particles |
| CN114145722A (en)* | 2021-12-07 | 2022-03-08 | 西安邮电大学 | A pulse pathological feature mining method for patients with pancreatitis |
| CN116109818A (en)* | 2023-04-11 | 2023-05-12 | 成都中医药大学 | Traditional Chinese medicine pulse condition distinguishing system, method and device based on facial video |
| CN116889388A (en)* | 2023-09-11 | 2023-10-17 | 长春理工大学 | An intelligent detection system and method based on rPPG technology |
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN1792319A (en)* | 2005-11-03 | 2006-06-28 | 浙江大学 | Automatic testing method for traditional Chinese medical pulse manifestation characteristics parameter |
| CN101408912A (en)* | 2008-11-21 | 2009-04-15 | 天津师范大学 | Method for automatically extracting characteristic function of traditional Chinese medicine pulse manifestation |
| US20140016840A1 (en)* | 2008-06-30 | 2014-01-16 | Nellcor Puritan Bennett Ireland | Systems and methods for ridge selection in scalograms of signals |
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN1792319A (en)* | 2005-11-03 | 2006-06-28 | 浙江大学 | Automatic testing method for traditional Chinese medical pulse manifestation characteristics parameter |
| US20140016840A1 (en)* | 2008-06-30 | 2014-01-16 | Nellcor Puritan Bennett Ireland | Systems and methods for ridge selection in scalograms of signals |
| CN101408912A (en)* | 2008-11-21 | 2009-04-15 | 天津师范大学 | Method for automatically extracting characteristic function of traditional Chinese medicine pulse manifestation |
| Title |
|---|
| 杨力华 等译: "《信号处理的小波导引》", 30 September 2002* |
| 王燕 等: "基于小波模极大原理的脉象特征提取研究", 《航天医学与医学工程》* |
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107736880A (en)* | 2017-10-24 | 2018-02-27 | 新绎健康科技有限公司 | A kind of pulse analysis method and system |
| CN107736880B (en)* | 2017-10-24 | 2024-01-30 | 新绎健康科技有限公司 | Pulse analysis method and system |
| CN107811619A (en)* | 2017-12-08 | 2018-03-20 | 西安科技大学 | Portable pulse-taking instrument and its analysis method |
| CN110595956A (en)* | 2019-08-09 | 2019-12-20 | 浙江工业大学 | A wear state mutation detection method based on fractal characteristics of abrasive particles |
| CN114145722A (en)* | 2021-12-07 | 2022-03-08 | 西安邮电大学 | A pulse pathological feature mining method for patients with pancreatitis |
| CN114145722B (en)* | 2021-12-07 | 2023-08-01 | 西安邮电大学 | Pulse pathological feature mining method for pancreatitis patients |
| CN116109818A (en)* | 2023-04-11 | 2023-05-12 | 成都中医药大学 | Traditional Chinese medicine pulse condition distinguishing system, method and device based on facial video |
| CN116109818B (en)* | 2023-04-11 | 2023-07-28 | 成都中医药大学 | Traditional Chinese medicine pulse condition distinguishing system, method and device based on facial video |
| CN116889388A (en)* | 2023-09-11 | 2023-10-17 | 长春理工大学 | An intelligent detection system and method based on rPPG technology |
| CN116889388B (en)* | 2023-09-11 | 2023-11-17 | 长春理工大学 | An intelligent detection system and method based on rPPG technology |
| Publication | Publication Date | Title |
|---|---|---|
| CN103932686A (en) | Method and device for extracting pulse condition signal | |
| CN108670209A (en) | Method and system for automatically identifying traditional Chinese medicine pulse condition | |
| CN105286815A (en) | Pulse wave signal feature point detection method based on waveform time domain features | |
| CN102293639B (en) | Pulse condition signal time domain feature extraction method | |
| CN106037694A (en) | Continuous blood pressure measuring device based on pulse waves | |
| CN109124610B (en) | Anti-interference method and device for non-invasive blood pressure measurement | |
| CN103584854B (en) | Extraction method of electrocardiosignal R waves | |
| CN101732033A (en) | Method and device for extracting characteristic parameter in human body waveform | |
| CN110141196A (en) | Peripheral arterial vessel elasticity evaluation method and system based on double triangle blood flow model | |
| CN104680186A (en) | Automatic classification method for ST-segment evaluation patterns in electrocardiograph signals | |
| CN108403094A (en) | Method for identifying pulse wave crest | |
| CN103829944B (en) | Based on the thoracic impedance signal processing method of pattern recognition | |
| CN104622440B (en) | The method and device of punctuate during a kind of extraction pulse wave | |
| CN104173030A (en) | Pulse wave starting point real-time detection method resisting waveform change interference and application thereof | |
| CN107898443B (en) | Method, device and computer storage medium for detecting dichotomous wave | |
| WO2020114448A1 (en) | Photoplethysmography signal feature point detecting method and device | |
| CN103284703B (en) | Aortic pulse wave transfer time measuring method based on upper extremity artery information | |
| CN107616782A (en) | A kind of electrocardiosignal quality determining method and device | |
| CN114652276A (en) | Pulse wave velocity detection method and device based on video image | |
| Lyle et al. | Beyond HRV: Analysis of ECG signals using attractor reconstruction | |
| CN100493445C (en) | Automatic Detection Method of TCM Pulse Condition Characteristic Parameters | |
| CN114795153A (en) | Fusion extraction method and device for pulse wave and electrocardiosignal high-quality waveform | |
| CN104622446A (en) | Heart rate variability signal optimization method based on KHM clustering algorithm | |
| CN102217931A (en) | Method and device for acquiring heart rate variation characteristic parameter | |
| CN102755154B (en) | Calculation method for extracting conduction time from pulse wave |
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| RJ01 | Rejection of invention patent application after publication | ||
| RJ01 | Rejection of invention patent application after publication | Application publication date:20140723 |