CN111528825A - Photoelectric volume pulse wave signal optimization method - Google Patents

Abstract

The invention discloses a method for optimizing a photoplethysmography signal, which comprises the following steps: s01, synchronously acquiring photoplethysmographic signals and motion state signals of the tested person; s02, performing a first round of optimization on the pulse wave signals, namely judging whether the tested person is in a static state or a slow walking state or a violent movement state according to the movement state signals, and if the tested person is in the violent movement state, rejecting the pulse wave signals in a corresponding time period; s03, denoising the pulse wave signals after the first round of optimization; s04, extracting characteristic indexes of the de-noised pulse wave signals; and S05, comparing the extracted characteristic indexes with a threshold range, and when at least one characteristic index falls outside the threshold range, marking the pulse wave period corresponding to the characteristic index as an abnormal period and removing the abnormal period to finally obtain the optimized pulse wave signal. This scheme is applicable to wearing formula medical equipment.

Description

A Photoplethysmographic Signal Optimization Method

technical field

The invention relates to the field of human blood pressure parameter processing, in particular to a photoelectric volume pulse wave signal optimization method.

Background technique

Photoplethysmography (PPG) is obtained by using photoelectric means to detect changes in blood volume based on the Beer-Lamber Law, and is referred to as PPG for short. It contains a large amount of useful information related to the physiological system, and has been widely used in clinical vascular assessment and monitoring of physiological parameters such as cardiac output, blood pressure, respiratory rate, heart rate and blood oxygen saturation.

After the fixed wavelength light is irradiated vertically to the skin surface through the LED, the light intensity obtained from the photoelectric receiver will be weakened due to the absorption of the light by the skin, muscles, blood and bones. By measuring the transmitted or reflected light to the receiver Changes in light intensity detect pressure-induced changes in blood volume, which are then plotted as a curve known as a photoplethysmography (PPG).

The PPG signal collected by the front-end hardware usually contains a lot of noise, generally including high-frequency noise caused by ambient light, dark current, power frequency interference, electromagnetic interference, electromyographic interference, etc., low-frequency baseline drift noise caused by physiological activities such as human respiration, and Movement interference.

SUMMARY OF THE INVENTION

The invention mainly solves the technical problem that the photoplethysmography (PPG) signal in the prior art contains more noise, and provides a photoplethysmographic optimization method for eliminating and reducing interference and noise.

The present invention mainly solves the above-mentioned technical problems through the following technical solutions: a method for optimizing a photoplethysmographic wave signal, comprising the following steps:

S01, synchronously acquiring the photoplethysmographic signal and the motion state signal of the subject;

S02, performing the first round of optimization on the pulse wave signal, that is, judging whether the subject is in a static state, a slow-walking state or a vigorous exercise state according to the motion state signal, and if it is in a vigorous exercise state, the pulse wave signal of the corresponding time period is eliminated;

S03, denoising the pulse wave signal after the first round of optimization;

S04, extracting characteristic indexes of the denoised pulse wave signal;

S05. Compare the extracted feature index with the threshold range, when at least one feature index falls outside the threshold range, mark the pulse wave cycle corresponding to this feature index as an abnormal cycle and remove it, and finally obtain an optimized pulse wave wave signal.

Preferably, in step S01, the motion state signal is obtained by a triaxial acceleration sensor fixed on the non-limb of the subject.

Preferably, in step S02, determining whether the subject is in a vigorous exercise state according to the exercise state signal is specifically:

Compute the resultant acceleration in three directions:

where a _t represents the resultant acceleration, and a _x , a _y and a _z represent the acceleration values on the X-axis, Y-axis and Z-axis measured by the three-axis acceleration sensor, respectively;

When the resultant acceleration is less than the upper limit of the first threshold and greater than the lower limit of the first threshold, it is determined that the subject is in a stationary state, the upper limit of the first threshold is 1.1g, and the lower limit of the first threshold is 0.9g, and the signal is a normal signal at this time;

When the combined acceleration is less than the lower limit of the second threshold or greater than the upper limit of the second threshold, it is determined that the subject is in a state of vigorous exercise. The lower limit of the second threshold is 0.8g, and the upper limit of the second threshold is 1.5g. At this time, the signal is an abnormal signal and is discarded. ;

When the resultant acceleration is less than or equal to the lower limit of the first threshold and greater than or equal to the lower limit of the second threshold, or the resultant acceleration is less than or equal to the upper limit of the second threshold and greater than the upper limit of the first threshold, it is determined that the subject is in a slow-walking state, and the signal at this time is a sub-normal signal .

Since the coordinate system of the three-axis accelerometer will change with the change of the subject's posture, only using the single-axis change to judge the exercise intensity will produce a large error, and cannot reflect the exercise state well, and the acceleration of the three axes is used at the same time. It will increase the complexity of the algorithm. Therefore, in order to better reflect the exercise intensity and at the same time make the calculation simpler, this scheme chooses to use a combination of accelerations in three directions to judge the exercise intensity. Two upper threshold limits and two lower threshold limits are used to complete the judgment of normal signals, sub-normal signals and abnormal signals at one time, which is convenient for subsequent processing such as blood pressure estimation model construction.

Preferably, in step S03, performing denoising processing on the pulse wave signal after the first round of optimization is specifically:

S301, remove high-frequency noise by a double-density wavelet threshold method, specifically:

The pulse wave signal after the first round of optimization is passed through the filtering system, and decomposed to obtain a low-frequency coefficient and two high-frequency coefficients; each layer of filtering system includes a low-pass filter h ₀ (n) and two high-pass filters h ₁ ( n) and h ₂ (n);

The low-frequency coefficients are passed through the filtering system again, and the process is repeated three times to complete the three-layer decomposition of the double-density wavelet;

The wavelet coefficients (that is, the signal decomposed by the low-frequency and high-frequency functions) are processed through the threshold function;

Inverse transform the processed wavelet coefficients to reconstruct the denoised signal;

S302, remove baseline drift noise by a cubic spline interpolation method.

Preferably, the threshold function is:

x>0, sgn(x)=1

x=0, sgn(x)=0

x<0, sgn(x)=-1

In the formula, s is the processed wavelet coefficient, x is the unprocessed wavelet coefficient, and T is the denoising threshold.

Preferably, the denoising threshold T is determined by the following formula:

where σ is the estimated noise variance, N is the signal length, and ω _b is obtained by the following process:

Square the wavelet coefficients of a certain layer, and then sort them from small to large to obtain the sequence W;

W=[ω ₁ ,ω ₂ ,...ω _N-1 ]

Calculate the risk value of each element in W in turn, where the risk value of the kth element ω _k is:

k=0,1,...N-1;

The minimum risk calculated by taking the wavelet coefficients of this layer is ω _b ;

The parameters θ and μ are obtained by the following formulas:

θ=(W-n)/n

Preferably, in the step S302, the removal of baseline drift noise by cubic spline interpolation is specifically:

The trough points of the pulse wave signal are identified by the findpeaks function, and then fitted by the cubic spline interpolation method to obtain the base drift of the signal. Finally, the pulse wave signal is subtracted from the base drift to obtain the signal after the baseline drift has been removed.

Preferably, in the step S04, the feature index extraction of the denoised pulse wave signal is specifically:

S401. Obtain the local peak value of the pulse wave signal through the findpeaks function;

S402. If the time difference between the two local peaks is less than 0.5 seconds, delete the point with a smaller amplitude; the remaining points are the main wave peak points of each cycle;

S403, multiplying the main wave crest by -1 to obtain the starting point of the cycle where the main wave crest is located, and subtracting 1 from the starting point to obtain the end point of the previous cycle;

S404. Calculate the characteristic index in units of cycles, and the calculation formula is as follows:

H1=H _{main wave peak point-} H _{starting point}

H2=H _{main wave peak point} - H _{end point}

ST=T _{main wave peak point -} T _{starting point}

DT=T _{end point-} T _{main wave peak point}

表示信号均值，std表示信号标准差。H表示位置点，T表示时间点。Among them, H1 is the rising range of the cycle, H2 is the falling range of the cycle, ST is the systolic time of the cycle, DT is the diastolic time of the cycle, K is the kurtosis of the cycle, and _xi represents the current pulse wave signal value,

Represents the signal mean, and std represents the signal standard deviation. H represents a location point, and T represents a time point.

Preferably, the threshold range in step S05 is determined in the following manner:

S501. For the subject, after being processed by a denoising algorithm, select template data of the subject in a stationary state and a slow-walking state with relatively small noise;

S502. Calculate the feature index sequence of the template data, and obtain the median line of each feature index sequence by using median filtering with a window length of 50;

S503. For each feature index sequence, the upper threshold is the sum of the offset of the median line and the upper threshold, the lower threshold is the sum of the offset of the median line and the lower threshold, and the range between the upper threshold and the lower threshold is the threshold range; the upper threshold offset is μ _seq +3σ _seq , the lower threshold offset is μ _seq -3σ _seq , μ _seq is the normal mean of the feature index sequence, and 3σ _seq is its normal standard deviation.

The substantial effect brought by the present invention is to eliminate abnormal signals seriously disturbed by motion and ensure the quality of PP signals; use double-density wavelet transform to remove high-frequency noise, use cubic spline interpolation to remove baseline drift noise, effectively remove noise, and compare Well preserved waveform characteristics.

Description of drawings

Fig. 1 is a kind of flow chart of the present invention;

Fig. 2 is a schematic diagram of a double density wavelet transform of the present invention.

Detailed ways

The technical solutions of the present invention will be further described in detail below through embodiments and in conjunction with the accompanying drawings.

Embodiment: A photoplethysmographic signal optimization method of this embodiment, as shown in FIG. 1 , includes the following steps:

S01, synchronously acquiring the photoplethysmographic signal and the motion state signal of the subject;

S02, performing the first round of optimization on the pulse wave signal, that is, judging whether the subject is in a static state, a slow-walking state or a vigorous exercise state according to the motion state signal, and if it is in a vigorous exercise state, the pulse wave signal of the corresponding time period is eliminated;

S03, denoising the pulse wave signal after the first round of optimization;

S04, extracting characteristic indexes of the denoised pulse wave signal;

S05. Compare the extracted feature index with the threshold range, when at least one feature index falls outside the threshold range, mark the pulse wave cycle corresponding to this feature index as an abnormal cycle and remove it, and finally obtain an optimized pulse wave wave signal.

In step S01, the motion state signal is obtained by a three-axis acceleration sensor fixed on the non-limb of the subject.

In step S02, judging whether the subject is in a vigorous exercise state according to the motion state signal is specifically: calculating the resultant acceleration in three directions:

where a _t represents the resultant acceleration, and a _x , a _y and a _z represent the acceleration values on the X-axis, Y-axis and Z-axis measured by the three-axis acceleration sensor, respectively;

When the resultant acceleration is less than the upper limit of the first threshold and greater than the lower limit of the first threshold, it is determined that the subject is in a stationary state, the upper limit of the first threshold is 1.1g, and the lower limit of the first threshold is 0.9g, and the signal is a normal signal at this time;

When the combined acceleration is less than the lower limit of the second threshold or greater than the upper limit of the second threshold, it is determined that the subject is in a state of vigorous exercise. The lower limit of the second threshold is 0.8g, and the upper limit of the second threshold is 1.5g. At this time, the signal is an abnormal signal and is discarded. ;

When the resultant acceleration is less than or equal to the lower limit of the first threshold and greater than or equal to the lower limit of the second threshold, or the resultant acceleration is less than or equal to the upper limit of the second threshold and greater than the upper limit of the first threshold, it is determined that the subject is in a slow-walking state, and the signal at this time is a sub-normal signal .

Since the coordinate system of the three-axis accelerometer will change with the change of the subject's posture, only using the single-axis change to judge the exercise intensity will produce a large error, and cannot reflect the exercise state well, and the acceleration of the three axes is used at the same time. It will increase the complexity of the algorithm. Therefore, in order to better reflect the exercise intensity and at the same time make the calculation simpler, this scheme chooses to use a combination of accelerations in three directions to judge the exercise intensity. Two upper threshold limits and two lower threshold limits are used to complete the judgment of normal signals, sub-normal signals and abnormal signals at one time, which is convenient for subsequent processing such as blood pressure estimation model construction.

In step S03, denoising the pulse wave signal after the first round of optimization is specifically:

S301, remove high-frequency noise by a double-density wavelet threshold method, specifically:

The pulse wave signal after the first round of optimization is passed through the filtering system, and decomposed to obtain a low-frequency coefficient and two high-frequency coefficients; each layer of filtering system includes a low-pass filter h ₀ (n) and two high-pass filters h ₁ ( n) and h ₂ (n);

The low-frequency coefficients are passed through the filtering system again, and the process is repeated three times to complete the three-layer decomposition of the double-density wavelet; Figure 2 shows a schematic diagram of this process;

The wavelet coefficients (that is, the signal decomposed by the low-frequency and high-frequency functions) are processed through the threshold function;

Inverse transform the processed wavelet coefficients to reconstruct the denoised signal;

S302, remove baseline drift noise by a cubic spline interpolation method.

The threshold function is:

x>0, sgn(x)=1

x=0, sgn(x)=0

x<0, sgn(x)=-1

In the formula, s is the processed wavelet coefficient, x is the unprocessed wavelet coefficient, and T is the denoising threshold.

The denoising threshold T is determined by the following formula:

where σ is the estimated noise variance, N is the signal length, and ω _b is obtained by the following process:

Square the wavelet coefficients of a certain layer, and then sort them from small to large to obtain the sequence W;

W=[ω ₁ ,ω ₂ ,...ω _N-1 ]

Calculate the risk value of each element in W in turn, where the risk value of the kth element ω _k is:

k=0,1,...N-1;

The minimum risk calculated by taking the wavelet coefficients of this layer is ω _b ;

The parameters θ and μ are obtained by the following formulas:

θ=(W-n)/n

In step S302, the removal of baseline drift noise by cubic spline interpolation is specifically:

The trough points of the pulse wave signal are identified by the findpeaks function, and then fitted by the cubic spline interpolation method to obtain the base drift of the signal. Finally, the pulse wave signal is subtracted from the base drift to obtain the signal after the baseline drift has been removed.

In the step S04, the characteristic index extraction of the denoised pulse wave signal is as follows:

S401. Obtain the local peak value of the pulse wave signal through the findpeaks function;

S402. If the time difference between the two local peaks is less than 0.5 seconds, delete the point with a smaller amplitude; the remaining points are the main wave peak points of each cycle;

S403, multiplying the main wave crest by -1 to obtain the starting point of the cycle where the main wave crest is located, and subtracting 1 from the starting point to obtain the end point of the previous cycle;

S404. Calculate the characteristic index in units of cycles, and the calculation formula is as follows:

H1=H _{main wave peak point-} H _{starting point}

H2=H _{main wave peak point} - H _{end point}

ST=T _{main wave peak point -} T _{starting point}

DT=T _{end point-} T _{main wave peak point}

表示信号均值，std表示信号标准差。H表示位置点，T表示时间点。Among them, H1 is the rising range of the cycle, H2 is the falling range of the cycle, ST is the systolic time of the cycle, DT is the diastolic time of the cycle, K is the kurtosis of the cycle, and _xi represents the current pulse wave signal value,

Represents the signal mean, and std represents the signal standard deviation. H represents a location point, and T represents a time point.

The threshold range in the step S05 is determined in the following manner:

S501. For the subject, after being processed by a denoising algorithm, select template data of the subject in a stationary state and a slow-walking state with relatively small noise;

S502. Calculate the feature index sequence of the template data, and obtain the median line of each feature index sequence by using median filtering with a window length of 50;

S503. For each feature index sequence, the upper threshold is the sum of the offset of the median line and the upper threshold, the lower threshold is the sum of the offset of the median line and the lower threshold, and the range between the upper threshold and the lower threshold is the threshold range; the upper threshold offset is μ _seq +3σ _seq , the lower threshold offset is μ _seq -3σ _seq , μ _seq is the normal mean of the feature index sequence, and 3σ _seq is its normal standard deviation.

The final optimized PPG signal can be used for blood pressure estimation and other purposes.

The specific embodiments described herein are merely illustrative of the spirit of the invention. Those skilled in the art to which the present invention pertains can make various modifications or additions to the described specific embodiments or substitute in similar manners, but will not deviate from the spirit of the present invention or go beyond the definitions of the appended claims range.

Although this paper uses more terms such as double density wavelet threshold method and feature index extraction, it does not exclude the possibility of using other terms. These terms are used only to more conveniently describe and explain the essence of the present invention; it is contrary to the spirit of the present invention to interpret them as any kind of additional limitation.

Claims (9)

1. a photoplethysmographic signal optimization method, is characterized in that, comprises the following steps: S01, synchronously acquiring the photoplethysmographic signal and the motion state signal of the subject; S02, performing the first round of optimization on the pulse wave signal, that is, judging whether the subject is in a static state, a slow-walking state or a vigorous exercise state according to the motion state signal, and if it is in a vigorous exercise state, the pulse wave signal of the corresponding time period is eliminated; S03, denoising the pulse wave signal after the first round of optimization; S04, extracting characteristic indexes of the denoised pulse wave signal; S05. Compare the extracted feature index with the threshold range, when at least one feature index falls outside the threshold range, mark the pulse wave cycle corresponding to this feature index as an abnormal cycle and remove it, and finally obtain an optimized pulse wave wave signal. 2 . The method for optimizing a photoplethysmography signal according to claim 1 , wherein, in step S01 , the motion state signal is obtained by a triaxial acceleration sensor fixed on the non-limb of the subject. 3 . 3. a kind of photoelectric volume pulse wave signal optimization method according to claim 2, is characterized in that, in step S02, according to motion state signal, judge whether the tested person is in vigorous motion state is specifically: calculate the resultant acceleration of three directions :

where a _t represents the resultant acceleration, and a _x , a _y and a _z represent the acceleration values on the X-axis, Y-axis and Z-axis measured by the three-axis acceleration sensor, respectively; When the resultant acceleration is less than the upper limit of the first threshold and greater than the lower limit of the first threshold, it is determined that the subject is in a stationary state, the upper limit of the first threshold is 1.1g, and the lower limit of the first threshold is 0.9g; When the combined acceleration is less than the lower limit of the second threshold or greater than the upper limit of the second threshold, it is determined that the subject is in a state of vigorous exercise, the lower limit of the second threshold is 0.8g, and the upper limit of the second threshold is 1.5g; When the combined acceleration is less than or equal to the first lower threshold and greater than or equal to the second lower threshold, or the combined acceleration is less than or equal to the second upper threshold and greater than the first upper threshold, it is determined that the subject is in a slow walking state. 4. a kind of photoplethysmographic pulse wave signal optimization method according to claim 1, is characterized in that, in step S03, the pulse wave signal after the first round of optimization is carried out denoising processing specifically: S301, remove high-frequency noise by a double-density wavelet threshold method, specifically: The pulse wave signal after the first round of optimization is passed through the filtering system, and decomposed to obtain a low-frequency coefficient and two high-frequency coefficients; each layer of filtering system includes a low-pass filter h ₀ (n) and two high-pass filters h ₁ ( n) and h ₂ (n); The low-frequency coefficients are passed through the filtering system again, and the process is repeated three times to complete the three-layer decomposition of the double-density wavelet; Process the wavelet coefficients through a threshold function; Inverse transform the processed wavelet coefficients to reconstruct the denoised signal; S302, remove baseline drift noise by a cubic spline interpolation method. 5. a kind of photoplethysmographic signal optimization method according to claim 4, is characterized in that, described threshold value function is:

x>0, sgn(x)=1 x=0, sgn(x)=0 x<0, sgn(x)=-1 In the formula, s is the processed wavelet coefficient, x is the unprocessed wavelet coefficient, and T is the denoising threshold. 6. a kind of photoplethysmography signal optimization method according to claim 5, is characterized in that, described denoising threshold T is determined by following formula:

where σ is the estimated noise variance, N is the signal length, and ω _b is obtained by the following process: Square the wavelet coefficients of a certain layer, and then sort them from small to large to obtain the sequence W; W=[ω ₁ ,ω ₂ ,...ω _N-1 ] Calculate the risk value of each element in W in turn, where the risk value of the kth element ω _k is:

The minimum risk calculated by taking the wavelet coefficients of this layer is ω _b ; The parameters θ and μ are obtained by the following formulas: θ=(W-n)/n

7. a kind of photoplethysmography signal optimization method according to claim 4 or 5 or 6, is characterized in that, described step S302 removes baseline drift noise by cubic spline interpolation and is specifically: The trough points of the pulse wave signal are identified by the findpeaks function, and then fitted by the cubic spline interpolation method to obtain the base drift of the signal. Finally, the pulse wave signal is subtracted from the base drift to obtain the signal after the baseline drift has been removed. 8. The method for optimizing a photoplethysmographic pulse wave signal according to claim 7, wherein, in the step S04, the characteristic index extraction of the denoised pulse wave signal is specifically: S401. Obtain the local peak value of the pulse wave signal through the findpeaks function; S402. If the time difference between the two local peaks is less than 0.5 seconds, delete the point with a smaller amplitude; the remaining points are the main wave peak points of each cycle; S403, multiplying the main wave crest by -1 to obtain the starting point of the cycle where the main wave crest is located, and subtracting 1 from the starting point to obtain the end point of the previous cycle; S404. Calculate the characteristic index in units of cycles, and the calculation formula is as follows: H1=H _{main wave peak point-} H _{starting point} H2=H _{main wave peak point} - H _{end point} ST=T _{main wave peak point -} T _{starting point} DT=T _{end point-} T _{main wave peak point}

Among them, H1 is the rising range of the cycle, H2 is the falling range of the cycle, ST is the systolic time of the cycle, DT is the diastolic time of the cycle, K is the kurtosis of the cycle, and _xi represents the current pulse wave signal value,

Represents the signal mean, and std represents the signal standard deviation. 9. The photoplethysmographic signal optimization method according to claim 8, wherein the threshold range is determined in the following manner: S501. For the subject, after processing by a denoising algorithm, select template data of the subject in a stationary state and a slow-walking state with relatively small noise; S502. Calculate the feature index sequence of the template data, and obtain the median line of each feature index sequence by using median filtering with a window length of 50; S503. For each feature index sequence, the upper threshold is the sum of the offset of the median line and the upper threshold, the lower threshold is the sum of the offset of the median line and the lower threshold, and the range between the upper threshold and the lower threshold is the threshold range; the upper threshold offset is μ _seq +3σ _seq , the lower threshold offset is μ _seq -3σ _seq , μ _seq is the normal mean of the feature index sequence, and 3σ _seq is its normal standard deviation.

Images