CN116763275A - Time-frequency analysis method and device for processing photoplethysmography signals - Google Patents

Abstract

The invention discloses a time-frequency analysis method for processing photoplethysm pulse signal. The photoplethysm wave signal is used as the original time domain signal; the signal is subjected to fractional Fourier transform at different orders, and each obtained Calculate the kurtosis coefficient of the fractional spectrum under the order, and perform peak detection to determine the number of peaks as the number of mode decompositions in the variational mode decomposition; construct a variational mode decomposition algorithm, obtain the intrinsic mode components, and calculate each The correlation coefficient between the eigenmode components and the original time domain signal is used to obtain the reconstructed signal based on the relationship between the two to achieve signal denoising; the reconstructed signal is subjected to generalized Warblet transformation and second-order transient extraction is performed to obtain Two-dimensional time-frequency feature representation. The invention can better remove the noise interference of the photoplethysmographic pulse signal, improve the aggregation degree of the signal's time-frequency distribution, and is well applied in the analysis of the photoplethysmographic pulse signal under the background of baseline drift, power frequency noise, electromagnetic interference and other backgrounds.

Description

A time-frequency analysis method and device for processing photoplethysmography signals

Technical Field

The present invention relates to the field of photoelectric signal processing, and in particular to a time-frequency analysis method and device for processing photoplethysmography signals.

Background Art

Photoplethysmography is a non-invasive physiological signal measurement technology, mainly used to measure the volume changes and pulse waveforms of blood vessels. This technology is based on the mechanism of pulse waves. Photodiodes and light sources are placed on human skin, and then the blood volume changes in blood vessels are measured through sensors. The signals generated by the changes in blood volume are converted into electrical signals, which can then be used to obtain relevant physiological parameters such as blood pressure, heart rate, heartbeat and other information through data processing and analysis. The photoplethysmography method has many advantages such as non-invasiveness, convenient detection, simple operation, stable performance, good repeatability, safety and no cross-infection, etc., which makes it not only suitable for clinical testing, monitoring, emergency physical fitness testing in hospitals, but also can be applied to community and family health care.

Since the human body is a complex system and the photoelectric signal is a relatively weak signal, slight changes in the human body or weak interference from the outside world will affect the final collected pulse wave waveform. In the process of collecting volume pulse waves, it is easy to be affected by baseline drift, power frequency noise and electromagnetic interference.

Summary of the invention

In order to solve the above problems, the present invention provides a time-frequency analysis method and device for processing photoplethysmography signals, wherein the method comprises the following steps:

S1. Use a light source to vertically illuminate a human finger, use a photodetector to receive the transmitted light, and obtain a photoelectric volume pulse wave signal as the original time domain signal;

S2, performing fractional Fourier transform of different orders on the original time domain signal to obtain fractional spectra of different orders;

S3, calculating the kurtosis coefficient of the fractional spectrum at each order, and performing peak detection on the kurtosis coefficient, and determining the number of peaks as the input parameter modal decomposition number in the variational modal decomposition algorithm;

S4, constructing a variational mode decomposition algorithm, taking the number of peaks as the mode decomposition number K in the variational mode decomposition algorithm, and obtaining K eigenmode components;

S5, calculating the correlation coefficient between each eigenmode component and the original time domain signal, setting a threshold TH, when the correlation coefficient between a certain eigenmode component and the original signal is greater than the threshold TH, the eigenmode component is retained, otherwise, the eigenmode component is removed, and the retained eigenmode components are summed to obtain the reconstructed signal;

S6. The reconstructed signal is combined with the generalized Warblet transform and the second-order transient extraction transform, namely the second-order transient extraction generalized Warblet transform (STEGWT), to obtain a two-dimensional time-frequency feature representation.

Furthermore, the light source in step S1 is 650nm red light or 940nm near-infrared light.

Furthermore, in step S2, fractional Fourier transform of different orders is performed on the original time domain signal and is expressed by the formula:

(1)

in, represents the p-order fractional Fourier transform spectrum, p is the rotation order, n is an arbitrary integer, α is the rotation angle in the time-frequency plane, and its relationship with the rotation order p is 2, is the amplitude of the fractional Fourier transform result, u is the fractional Fourier transform domain, referred to as the fractional domain, is the original time domain signal, j is the imaginary unit;

Based on the above formula, calculate the signal at different orders The fractional Fourier transform under the condition is used to obtain the fractional spectrum. , Represents the square of the absolute value.

Furthermore, in step S3, the kurtosis coefficient of the fractional spectrum at each order is calculated as follows:

(2)

in, is the kurtosis coefficient at order p, E represents the expected value operator, represents the p-order fractional Fourier transform spectrum, for The mean of Indicates absolute value.

Further, in step S5, the correlation coefficient between each eigenmode component and the original time domain signal is calculated, and the threshold TH is set to be expressed by the formula:

(3)

(4)

in, is the jth eigenmode component and the original time domain signal The correlation coefficient of N is the number of sampling points of the signal. express The mean of represents the jth eigenmode component, express The mean of , n is equal to the number of correlation coefficients, express The mean of .

Further, step S6 is specifically as follows:

The generalized Warblet transform expression for the reconstructed signal is:

(5)

in, represents the generalized Warblet transform of the reconstructed signal, ω represents the frequency of the reconstructed signal in the time-frequency domain, is a sliding window function, t represents the time when the window center is sliding. A non-negative symmetric normalized real window is defined, j is the imaginary unit, The calculation formula is as follows:

(6)

in, represents the reconstructed signal, represents the time variable of the reconstructed signal, represents an intermediate variable, and are the frequency translation operator and the frequency rotation operator respectively, j is the imaginary unit, m represents the total number of sine functions or cosine functions, and are the Fourier coefficients, is the corresponding harmonic component frequency;

The impact component of the strong pulse component is expressed by the Dirac delta function:

(7)

in, The impulsive component of the strong pulse component is represented by the Dirac delta function. represents the Dirac delta function, A is the amplitude of the impulse component, The moment when the impact occurs;

Substituting the impact expressed by equation (7) into equation (5) yields:

(8)

in, is the window function;

The spread energy of the generalized Warblet transform of the reconstructed signal is redistributed using the two-dimensional group delay estimate and is written as:

(9)

in, represents the two-dimensional group delay estimate, i is the imaginary unit, is the symbol for partial derivative with respect to ω, Re() means taking the real part;

Substituting formula (8) into formula (9), we have:

(10)

Build a second-order frequency-dependent model:

(11)

in, represents the second-order frequency-dependent model, represents the amplitude of the second-order frequency-dependent model, is the phase of the second-order frequency-varying model, ω represents the frequency of the reconstructed signal in the time-frequency domain, represents the phase of the original time domain signal in the frequency domain, i is an imaginary unit, ξ represents the frequency domain of the second-order frequency-varying model, It means to find the first-order derivative, It means to find the second-order derivative;

The generalized Warblet transform expression in the frequency domain of the reconstructed signal under the second-order frequency-varying model is:

(12)

in, represents the generalized Warblet transform in the frequency domain of the reconstructed signal under the second-order frequency-varying model, represents the second-order frequency-dependent model, represents the window function in the frequency domain, i is the imaginary unit, for The generalized Warblet transform of , , represents the standard deviation of the window function. Substituting formula (11) into formula (12) yields:

(13)

The first-order two-dimensional group delay estimate for the second-order frequency-dependent model is:

(14)

According to formula (14), we get:

(15)

Define a new two-dimensional estimate :

(16)

in, is the sign of the partial derivative with respect to t;

Therefore, the second-order two-dimensional group delay of the second-order frequency-dependent model is expressed as follows:

(17)

in, is the second-order two-dimensional group delay estimate of the second-order frequency-dependent model, Represents the first-order two-dimensional group delay estimate of the second-order frequency-dependent model;

Therefore, the second-order transient transformation of formula (11) is:

(18)

in, represents the Dirac delta function, i is the imaginary unit, The results of the generalized Warblet transform for second-order transient extraction are obtained by effectively combining the fractional Fourier transform technique with variational mode decomposition.

The present invention also provides a time-frequency analysis device for processing a photoplethysmogram signal, the device comprising:

processor;

a memory having stored thereon a computer program executable on the processor;

Wherein, when the computer program is executed by the processor, a time-frequency analysis method for processing photoplethysmography signals is implemented.

The beneficial effects brought by the technical solution provided by the present invention are:

The scheme of the present invention is to perform denoising processing and time-frequency analysis on a photoplethysmogram signal. In the process of realizing FrVMD (combining FrFT technology and VMD to form an FrVMD algorithm), the kurtosis coefficient of FrFT (fractional Fourier transform) under different rotation orders is selected for peak detection. The obtained number of peaks is used as the number of modal decomposition of VMD (variational modal decomposition) as an input condition to perform modal decomposition and reconstruction on the photoplethysmogram signal. The photoplethysmogram signal can be decomposed and reconstructed to remove false and meaningless intrinsic mode components, thereby filtering out interference noise to a certain extent, and selecting a generalized Warblet transform (GWT: generalized Warblet transform) to obtain the time-frequency joint distribution of the signal; the second-order transient-extracting transform (STET: Second-order transient-extracting transform) is applied to the optimization of the GWT time-frequency spectrum to improve the aggregation of the time-frequency distribution and greatly reduce the influence of noise interference on the time-frequency representation of the signal. The present invention can effectively remove noise interference of photoplethysmography signal data, improve the aggregation of time-frequency distribution of photoplethysmography signal data, and is well applied in photoplethysmography signal data analysis under the background of baseline drift, power frequency noise and electromagnetic interference.

BRIEF DESCRIPTION OF THE DRAWINGS

1 is a flowchart of a time-frequency analysis method for processing a photoplethysmography signal according to an embodiment of the present invention;

FIG2 is a FrVMD decomposition result of a photoplethysmography signal 1 according to an embodiment of the present invention;

FIG3 is a reconstructed result of the photoplethysmography signal 1 after FrVMD decomposition according to an embodiment of the present invention, wherein (a) in FIG3 is the reconstructed time domain signal 1; (b) in FIG3 is an enlarged version of the reconstructed time domain signal 1 at the rectangular frame; (c) in FIG3 is the original time domain signal 1; and (d) in FIG3 is an enlarged version of the original time domain signal 1 at the rectangular frame;

FIG4 is a reconstructed result of the FrVMD decomposition of the photoplethysmography signal 1 according to an embodiment of the present invention, wherein (a) in FIG4 is the spectrum of the reconstructed signal 1; (b) in FIG4 is an enlarged version of the spectrum of the reconstructed signal 1 at the rectangular frame; (c) in FIG4 is the spectrum of the original signal 1; and (d) in FIG4 is an enlarged version of the spectrum of the original signal 1 at the rectangular frame;

FIG5 is a FrVMD decomposition result of a photoplethysmography signal 2 according to an embodiment of the present invention;

FIG6 is a reconstructed result after FrVMD decomposition of the photoplethysmography signal 2 according to an embodiment of the present invention, wherein (a) in FIG6 is the reconstructed time domain signal 2; (b) in FIG6 is an enlarged version of the reconstructed time domain signal 2 at the rectangular frame; (c) in FIG6 is the original time domain signal 2; and (d) in FIG6 is an enlarged version of the original time domain signal 2 at the rectangular frame;

Figure 7 is the reconstruction result of the photoplethysmography signal 2 after FrVMD decomposition according to an embodiment of the present invention, where (a) in Figure 7 is the spectrum of the reconstructed signal 2; (b) in Figure 7 is an enlarged version of the spectrum of the reconstructed signal 2 at the rectangular frame; (c) in Figure 7 is the spectrum of the original signal 2; and (d) in Figure 7 is an enlarged version of the spectrum of the original signal 2 at the rectangular frame.

DETAILED DESCRIPTION

In order to make the objectives, technical solutions and advantages of the present invention more clear, the embodiments of the present invention will be further described below in conjunction with the accompanying drawings.

The present invention effectively combines the fractional Fourier transform (FrFT) technology and the variational mode decomposition (VMD) to form the FrVMD algorithm, realizes the denoising of the photoplethysmography signal, and effectively combines the generalized Warblet transform (GWT: generalized Warblet transform) technology and the second-order transient-extracting transform (STET: Second-order transient-extracting transform) to form the STEGWT algorithm, realizes the time-frequency analysis of the photoplethysmography signal. The signal is decomposed by FrVMD, and whether the modal component is retained is determined according to whether the correlation between the modal component and the original signal exceeds the threshold. Finally, the modal components retained are superimposed to obtain the reconstructed photoplethysmography signal to achieve the purpose of denoising. The signal is processed by the GWT method to obtain the generalized Warblet transform time-frequency distribution of the signal, and then the STET technology performs second-order transient extraction based on the processing result of the GWT to obtain a time-frequency distribution with high concentration.

A flow chart of a time-frequency analysis method for processing a photoplethysmography signal according to an embodiment of the present invention is shown in FIG1 , and includes the following steps:

S1. Use a light source to vertically illuminate the human finger, use a photodetector to receive the transmitted light, and obtain a photoelectric volume pulse wave signal as the original time domain signal.

The light source in this embodiment uses 650nm red light or 940nm near-infrared light.

S2. Perform fractional Fourier transform (FrFT) of different orders on the original time domain signal to obtain fractional spectra of different orders. Determine the search step of FrFT, and obtain N rotation orders p in [0,2] according to the step; calculate the fractional Fourier transform of the signal at different orders p to obtain the fractional spectrum .

The fractional Fourier transform of the original time domain signal at different orders is expressed by the formula:

(1)

in, represents the p -order fractional Fourier transform spectrum, p is the rotation order, n is an arbitrary integer, α is the rotation angle in the time-frequency plane, and its relationship with the rotation order p is 2, is the amplitude of the fractional Fourier transform result, u is the fractional Fourier transform domain, referred to as the fractional domain, is the original time domain signal, and j is the imaginary unit.

Based on the above formula, calculate the signal at different orders The fractional Fourier transform under the condition is used to obtain the fractional spectrum. , Represents the square of the absolute value.

S3. Calculate the kurtosis coefficient of the fractional spectrum at each order, perform peak detection on the kurtosis coefficient, and determine the number of peaks as the input parameter modal decomposition number in the variational mode decomposition algorithm.

The calculation formula for the kurtosis coefficient of the fractional spectrum at each order is:

(2)

in, is the kurtosis coefficient at order p , E represents the expected value operator, represents the p- order fractional Fourier transform spectrum, for The mean of Indicates absolute value.

The kurtosis coefficient of the fractional spectrum under each order is calculated to obtain N kurtosis coefficients; the peak value of the kurtosis coefficient is detected to determine the number of peaks, and the number of peaks is used as the input parameter mode decomposition number K value in the VMD algorithm.

S4. Construct a variational mode decomposition algorithm, use the number of peaks as the mode decomposition number K in the variational mode decomposition algorithm, and obtain K eigenmode components.

S5. Calculate the correlation coefficient between each eigenmode component and the original time domain signal, set the threshold TH, and when the correlation coefficient between a certain eigenmode component and the original signal is greater than the threshold TH, retain the eigenmode component, otherwise, remove the eigenmode component, and sum the retained eigenmode components to obtain the reconstructed signal. This reconstruction method can determine which eigenmode components are the real components of the signal and which are false and meaningless eigenmode components. The reconstructed signal is recorded as x(τ) .

The correlation coefficient between each eigenmode component and the original time domain signal and the set threshold TH are expressed as follows:

(3)

(4)

in, is the jth eigenmode component and the original time domain signal The correlation coefficient of N is the number of sampling points of the signal. express The mean of represents the jth eigenmode component, express The mean of , n is equal to the number of correlation coefficients, express The mean of .

S6. The reconstructed signal is subjected to a two-dimensional time-frequency feature representation by combining the generalized Warblet transform with the second-order transient extraction transform.

The generalized Warblet transform expression for the reconstructed signal is:

(5)

in, represents the generalized Warblet transform of the reconstructed signal, ω represents the frequency of the reconstructed signal in the time-frequency domain, is a sliding window function, t represents the time when the window center is sliding. Defines a non-negative symmetric normalized real window, usually a Gaussian window, j is an imaginary unit, The calculation formula is as follows:

(6)

in, represents the reconstructed signal, represents the time variable of the reconstructed signal, represents an intermediate variable, and are the frequency translation operator and the frequency rotation operator respectively, m represents the total number of sine functions or cosine functions, and are the Fourier coefficients, is the corresponding harmonic component frequency;

Will Abbreviated as ,

The impact component of the strong pulse component is expressed by the Dirac delta function:

(7)

Among them, s(t) represents the impulsive component of the strong pulse component expressed by the Dirac delta function, that is, the impulse function of the pulse, represents the Dirac delta function, A is the amplitude of the impulse component, The moment when the impact occurs;

Substituting the impact expressed by formula (7) into formula (5) Represents the impulse characteristics of the signal and then integrates to obtain:

(8)

Among them, at this time , Δ represents the time support range of the window function, is the window function;

Ideally, the GWT transformation of the Dirac delta function should be concentrated at the time of the impact. The diffusion energy of the generalized Warblet transform that redistributes the reconstructed signal is estimated using the two-dimensional (2D) group delay (GD), written as:

(9)

in, represents the first-order two-dimensional group delay, i is the imaginary unit, is the symbol for partial derivative with respect to ω, Re() means taking the real part;

Substituting formula (8) into formula (9), we have:

(10)

Equation (10) shows that the 2D GD estimate of the Dirac delta function calculated using equation (7) is related to the pulse occurrence time instant Consistent.

Build a second-order frequency-dependent model:

(11)

in, represents the second-order generalized Warblet transform of the reconstructed signal, that is, the second-order frequency-varying model, represents the amplitude of the second-order frequency-dependent model, is the phase of the second-order frequency-variable model, represents the phase of the original time domain signal in the frequency domain, i is an imaginary unit, ξ represents the frequency domain of the second-order frequency-varying model, It means to find the first-order derivative, It means to find the second-order derivative.

The generalized Warblet transform expression in the frequency domain of the reconstructed signal under the second-order frequency-varying model is:

(12)

in, represents the generalized Warblet transform in the frequency domain of the reconstructed signal under the second-order frequency-varying model, represents the second-order frequency-dependent model, represents the window function in the frequency domain, i is the imaginary unit, , , represents the standard deviation of the window function, for The generalized Warblet transform of , substituting formula (11) into formula (12) yields:

(13)

The first-order two-dimensional group delay estimate for the second-order frequency-dependent model is:

(14)

According to formula (14), we get:

(15)

Equation (15) provides an indication of how to compute the new 2D GD estimate in second-order signal analysis. Therefore, a new two-dimensional estimate is defined as :

(16)

in, To find the sign of the partial derivative with respect to t ;

According to the new two-dimensional group delay calculation method defined in equation (16), and considering the relevant constraints, the second-order two-dimensional group delay of the second-order frequency-variant model is expressed as follows:

(17)

in, This is the second-order two-dimensional group delay estimate of the second-order frequency-varying model, Represents the first-order two-dimensional group delay estimate of the second-order frequency-dependent model;

The second-order transient transformation of formula (11) is:

(18)

therefore, That is, FrVMD is the result of the generalized Warblet transform for second-order transient extraction, which can generate a high-resolution time-frequency representation for the second-order frequency-varying signal.

The algorithm of the present invention can effectively remove the noise of the signal and obtain the joint distribution of the time-frequency of the signal, which is helpful to further understand the changes in the time-frequency distribution of photoplethysmography data under the background of baseline drift, power frequency noise and electromagnetic interference, and provide a basis and help for the research on photoplethysmography signal recognition. In order to verify the effectiveness of the algorithm in the research on photoplethysmography data processing under the background of baseline drift, power frequency noise and electromagnetic interference, this paper selects two groups of measured data 1 with noise and 2 with baseline drift for experimental verification.

1. FrVMD denoising experiment

(1) Measured data with noise 1

Firstly, FrVMD is used to perform noise reduction on the noisy measured data 1. Figure 2 is the decomposition result of FrVMD. Table 1 shows the correlation and threshold of each modal component with the original signal.

Table 1

As shown in Figure 2 and Table 1, the photoplethysmography signal is decomposed into three modal components. Only modal component 1 exceeds the threshold, that is, it has a high correlation with the original signal. Therefore, modal component 1 is considered to be a valid component, and modal components 2 and 3 are invalid components that should be removed. Then the remaining modal components are summed to obtain the final reconstruction result.

Figure 3 is the reconstruction result after FrVMD decomposition of the photoplethysmography signal 1 according to an embodiment of the present invention, (a) in Figure 3 is the reconstructed time domain signal 1; (b) in Figure 3 is an enlarged version of the reconstructed time domain signal 1 at the rectangular frame; (c) in Figure 3 is the original time domain signal 1; (d) in Figure 3 is an enlarged version of the original time domain signal 1 at the rectangular frame, and Figure 4 is the reconstruction result after FrVMD decomposition of the photoplethysmography signal 1 according to an embodiment of the present invention, (a) in Figure 4 is the spectrum of the reconstructed signal 1; (b) in Figure 4 is an enlarged version of the spectrum of the reconstructed signal 1 at the rectangular frame; (c) in Figure 4 is the spectrum of the original signal 1; (d) in Figure 4 is an enlarged version of the spectrum of the original signal 1 at the rectangular frame.

As shown in Figure 3, the original time domain signal contains more noise, and the time domain curve shows more burrs. The FrVMD method can effectively remove the noise of the signal, making the reconstructed time domain signal curve smoother and free of burrs. As shown in Figure 4, the spectrum of the reconstructed signal does not contain redundant noise spectrum, while the original signal, except for the spectrum of the true component at the center, obviously contains the spectrum of redundant noise. Therefore, the FrVMD method can effectively remove the noise of the signal.

In order to further compare the performance of the FrVMD reconstruction method, the mean absolute percentage error (MAPE), root mean square error (RMSE), and determination coefficient (R2) are selected to judge the decomposition and reconstruction denoising effect, and their definitions are shown in equations (19), (20), and (21). The larger the mean absolute error (MAE) and root mean square error (RMSE), the less ideal the denoising effect of the signal. The closer the determination coefficient (R2) is to 1, the better the denoising effect.

(19)

(20)

(twenty one)

Table 2 shows the different error values of the signal. It can be seen from Table 2 that the error of the FrVMD algorithm in decomposing and reconstructing the photoplethysmography signal is small, and the reconstructed signal has a high correlation with the original signal, indicating that the FrVMD algorithm has a more obvious noise removal effect on the basis of retaining the effective real components of the original signal.

Table 2

(2) Measured data including baseline drift 2

FrVMD is used to perform noise reduction on the measured data 2 containing baseline drift. Figure 5 shows the decomposition result of FrVMD, and Table 3 shows the correlation and threshold of each modal component with the original signal. As shown in Figure 5 and Table 3, the photoelectric volumetric pulse signal is decomposed into 10 modal components. Only modal component 1 and modal component 2 exceed the threshold, that is, they are highly correlated with the original signal. Then modal component 1 and modal component 2 are considered to be valid components, and the remaining modal components are invalid components and should be removed. Then the remaining modal components are summed to obtain the final reconstruction result.

Table 3

Figure 6 is the reconstruction result after FrVMD decomposition of the photoplethysmography signal 2 according to an embodiment of the present invention, where (a) in Figure 6 is the reconstructed time domain signal 2; (b) in Figure 6 is an enlarged version of the reconstructed time domain signal 2 at the rectangular frame; (c) in Figure 6 is the original time domain signal 2; and (d) in Figure 6 is an enlarged version of the original time domain signal 2 at the rectangular frame.

Figure 7 is the reconstruction result of the photoplethysmography signal 2 after FrVMD decomposition according to an embodiment of the present invention, where (a) in Figure 7 is the spectrum of the reconstructed signal 2; (b) in Figure 7 is an enlarged version of the spectrum of the reconstructed signal 2 at the rectangular frame; (c) in Figure 7 is the spectrum of the original signal 2; and (d) in Figure 7 is an enlarged version of the spectrum of the original signal 2 at the rectangular frame.

As shown in Figure 6, the original time domain signal contains baseline drift, as shown in Figure 6 (d) in the time domain curve. The FrVMD method can effectively remove the baseline drift of the signal, making the reconstructed time domain signal curve more regular, as shown in Figure 6 (b). As shown in Figure 7, the spectrum of the reconstructed signal does not contain redundant noise spectrum, while the original signal, except for the spectrum of the true component at the center, obviously contains the spectrum of redundant noise. Therefore, the FrVMD method can effectively remove the noise of the signal.

Table 4 shows the different error values of signal 2. It can be seen from Table 4 that the error of the FrVMD algorithm in decomposing and reconstructing the photoplethysmography signal is small, and the reconstructed signal has a high correlation with the original signal, indicating that the FrVMD algorithm has a more obvious effect in removing interference such as baseline drift while retaining the effective real components of the original signal.

Table 4

2. STEGWT time-frequency analysis experiment

(1) Measured data with noise 1

Firstly, STEGWT is used to perform time-frequency analysis on the noisy measured data 1. The time-frequency curves obtained by the FrVMD-STEGWT method are relatively concentrated and the time-frequency aggregation is high.

The higher the time-frequency aggregation, the more accurate the local positioning and resolution of the signal frequency components by the time-frequency distribution, which means the better the processing effect of the time-frequency analysis method. This paper uses Rayleigh entropy as a performance evaluation index to measure the aggregation of time-frequency distribution, and its calculation formula is: ,in, represents the number of Rayleigh entropy, and in this paper, ; Represents the time-frequency distribution of the signal. Rayleigh entropy The smaller the value of , the higher the aggregation of the time-frequency distribution.

Table 5 shows the Rayleigh entropy of signal 1 with different time-frequency distributions. The Rayleigh entropy value of FrVMD-STEGWT is 9.6842, which is much smaller than that of the three time-frequency analysis methods, SST, GWT, and WVD. That is, the time-frequency aggregation of FrVMD-STEGWT is relatively high, which proves that this improved method can effectively improve the time-frequency aggregation when processing photoelectric volume pulse signals with noise interference.

Table 5

(2) Measured data including baseline drift 2

The time-frequency curves obtained by the FrVMD-STEGWT method are relatively concentrated, and the time-frequency aggregation is high.

Table 6 shows the Rayleigh entropy of different time-frequency distributions of signal 2. The higher the time-frequency aggregation, the more accurate the local positioning and resolution of the signal frequency components by the time-frequency distribution, which indicates that the processing effect of the time-frequency analysis method is better. The Rayleigh entropy value of FrVMD-STEGWT is 11.3464, which is much smaller than that of the three time-frequency analysis methods of SST, GWT, and WVD, that is, the time-frequency aggregation of FrVMD-STEGWT is higher, which proves that this improved method can effectively improve its time-frequency aggregation when processing photoelectric volume pulse signals with baseline drift.

Table 6

The present invention also provides a time-frequency analysis device for processing a photoplethysmogram signal, comprising:

processor;

a memory having stored thereon a computer program executable on the processor;

When the computer program is executed by the processor, a time-frequency analysis method for processing photoplethysmography signals is implemented.

The above description of the disclosed embodiments enables those skilled in the art to implement or use the present invention. Various modifications to these embodiments will be apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the present invention. Therefore, the present invention will not be limited to the embodiments shown herein, but rather to the widest scope consistent with the principles and novel features disclosed herein.

Claims (7)

1. A time-frequency analysis method for processing photoplethysmography signals, comprising the steps of:

s1, vertically irradiating a finger penetrating through a human body by using a light source, and receiving the penetrated light by using a photoelectric detector to obtain a photoelectric volume pulse wave signal serving as an original time domain signal;

s2, carrying out fractional Fourier transform on the original time domain signal under different orders to obtain fractional spectrums under different orders;

s3, calculating kurtosis coefficients of the fractional spectrums under each order, carrying out peak detection on the kurtosis coefficients, and determining the number of peaks as the number of modal decomposition of input parameters in a variable modal decomposition algorithm;

s4, constructing a variation modal decomposition algorithm, and taking the number of peaks as the number K of modal decomposition in the variation modal decomposition algorithm to obtain K eigenmode components;

s5, calculating a correlation coefficient between each eigenmode component and an original time domain signal, setting a threshold value TH, reserving the eigenmode component when the correlation coefficient between one eigenmode component and the original signal is larger than the threshold value TH, otherwise, eliminating the eigenmode component, and summing the reserved eigenmode components to obtain a reconstructed signal;

s6, converting the reconstructed signal by combining generalized wavelet transformation and second-order transient extraction transformation to obtain the time-frequency characteristic representation of the reconstructed signal.

2. A time-frequency analysis method for processing photoplethysmography signals according to claim 1, wherein the light source in step S1 is 650nm red light or 940nm near infrared light.

3. A time-frequency analysis method for processing photoplethysmography signals according to claim 1, wherein the fractional fourier transform of the original time-domain signal at different orders in step S2 is formulated as:

（1）

wherein ,representation ofpThe order-fractional order fourier transform spectrum,pin order to be able to rotate the order,nis an arbitrary integer number of the whole,αis the rotation angle and rotation order of the time-frequency planepThe relation of->2，For the magnitude of the fractional fourier transform result,uis a fractional Fourier transform domain, which is simply called a fractional domain, ">As the original time-domain signal is to be processed,jis an imaginary unit;

based on the above formula, the signals are calculated at different ordersThe lower fractional Fourier transform to obtain a fractional spectrum +.>，Representing the square of the absolute value.

4. The method according to claim 1, wherein in step S3, the kurtosis coefficient calculation formula for the fractional spectrum at each order is:

（2）

wherein ,is thatpThe kurtosis coefficient at the level of the order,Erepresenting the expected value operator, ++>Representation ofpThe order-fractional order fourier transform spectrum,is->Mean of (2)，Representing the absolute value.

5. A time-frequency analysis method for processing photoplethysmography signals according to claim 1, characterized in that in step S5, a correlation coefficient between each eigenmode component and the original time-domain signal is calculated, and a set threshold TH is formulated as:

（3）

（4）

wherein ,is the firstjThe individual eigenmode components are +.>Is used for the correlation coefficient of (a),Nthe number of samples to be taken for a signal,representation->Mean value of->Represent the firstjThe eigenmode components>Representation->Is used for the average value of (a),nis equal to the number of correlation coefficients,representation->Is a mean value of (c).

6. A time-frequency analysis method for processing photoplethysmography signals according to claim 1, characterized in that step S6 is specifically:

the reconstructed signal is expressed as follows by generalized wavelet transform:

（5）

wherein ,generalized wavelet transform representing a reconstructed signal, ω representing the frequency in the time-frequency domain of the reconstructed signal,/v>Is a sliding window function, t represents the time of window center when the time window slides, +.>A normalized real window of non-negative symmetry is defined,jis imaginary unit, ++>The calculation formula is as follows:

（6）

wherein ,representing the reconstructed signal, < >>Time variable representing the reconstructed signal, +.>Which represents an intermediate variable that is used to control the operation of the device, andA frequency translation operator and a frequency rotation operator respectively,jin units of imaginary numbers,mrepresenting the total number of sine functions or cosine functions, +.> andFor Fourier coefficients +.>For the corresponding harmonic component frequencies;

the impact component of the strong pulse component is expressed as a dirac delta function:

（7）

wherein ,dirac for impact component representing strong pulse componentδRepresentation of a function->Representing diracδThe function of the function is that,Ａfor the impact component amplitude, +.>For the moment of occurrence of the impact；

Substituting the impact represented by formula (7) into formula (5) to obtain:

（8）

wherein ,is a window function;

redistributing the diffuse energy of the reconstructed signal generalized wavelet transform using two-dimensional group delay estimation, written as:

（9）

wherein ,representing a two-dimensional group delay estimate,iis imaginary unit, ++>To pair(s)ωSymbols of the bias derivatives are calculated, and Re () represents a real part;

substituting formula (8) into formula (9) includes:

（10）

establishing a second-order frequency-dependent model:

（11）

wherein ,representing a second order frequency-dependent model +.>Representing the magnitude of the second order frequency-dependent model, +.>For the phase of the second order frequency-dependent model,ωrepresenting the frequency in the time-frequency domain of the reconstructed signal, < >>Representing the phase of the original time domain signal in the frequency domain,iis an imaginary unit of number and is,ξfrequency domain representing a second order frequency-variant model, +.>Representing the first derivative,Representing the second derivative;

the generalized wavelet transform expression in the reconstructed signal frequency domain under the second-order frequency-variant model is:

（12）

wherein ,representing generalized wavelet transform in reconstructed signal frequency domain under second order frequency-variant model, +.>Representing a second order frequency-dependent model +.>Representing a window function in the frequency domain,iis imaginary unit, ++>Is->Is a generalized wavelet transform of +.>，，Expressing the standard deviation of the window function, substituting equation (11) into equation (12) yields:

（13）

the first-order two-dimensional group delay estimation of the second-order frequency-dependent model is as follows:

（14）

obtained according to formula (14):

（15）

defining a new two-dimensional estimate：

（16）

wherein ,to pair(s)tCalculating the symbols of the bias guide;

thus, the second-order two-dimensional group delay of the second-order frequency-dependent model is expressed as follows:

（17）

wherein ,second order two-dimensional group delay estimation for second order frequency-dependent model,/->A first-order two-dimensional group delay estimation representing a second-order frequency-dependent model;

thus, the second order transient transformation of equation (11) is:

（18）

wherein ,representing diracδThe function of the function is that,iis imaginary unit, ++>The generalized wavelet transform results are extracted for the second order transient under the effective combination of the fractional Fourier transform technique and the variant modal decomposition.

7. A time-frequency analysis device for processing photoplethysmography signals, the device comprising:

a processor;

a memory having stored thereon a computer program executable on the processor;

wherein the computer program, when executed by the processor, implements a time-frequency analysis method for processing photoplethysmography signals according to any one of claims 1 to 6.