CN116763275B - Time-frequency analysis method and device for processing photoplethysmography signals - Google Patents
Time-frequency analysis method and device for processing photoplethysmography signals Download PDFInfo
- Publication number
- CN116763275B CN116763275B CN202311054338.0A CN202311054338A CN116763275B CN 116763275 B CN116763275 B CN 116763275B CN 202311054338 A CN202311054338 A CN 202311054338A CN 116763275 B CN116763275 B CN 116763275B
- Authority
- CN
- China
- Prior art keywords
- frequency
- order
- signal
- representing
- time
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000004458 analytical method Methods 0.000 title claims abstract description 28
- 238000012545 processing Methods 0.000 title claims abstract description 28
- 238000013186 photoplethysmography Methods 0.000 title claims abstract description 19
- 238000001228 spectrum Methods 0.000 claims abstract description 46
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 36
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 20
- 230000009466 transformation Effects 0.000 claims abstract description 15
- 230000001052 transient effect Effects 0.000 claims abstract description 13
- 238000004364 calculation method Methods 0.000 claims abstract description 9
- 238000001514 detection method Methods 0.000 claims abstract description 8
- 238000000605 extraction Methods 0.000 claims abstract description 8
- 230000006870 function Effects 0.000 claims description 35
- 238000000034 method Methods 0.000 claims description 23
- 230000001419 dependent effect Effects 0.000 claims description 21
- 238000004590 computer program Methods 0.000 claims description 6
- 230000001678 irradiating effect Effects 0.000 claims description 3
- 238000013519 translation Methods 0.000 claims description 3
- 230000000149 penetrating effect Effects 0.000 claims description 2
- 230000002776 aggregation Effects 0.000 abstract description 14
- 238000004220 aggregation Methods 0.000 abstract description 14
- 238000009826 distribution Methods 0.000 abstract description 14
- 230000009467 reduction Effects 0.000 abstract description 3
- 230000000694 effects Effects 0.000 description 7
- 230000008859 change Effects 0.000 description 5
- 238000005516 engineering process Methods 0.000 description 5
- 230000008901 benefit Effects 0.000 description 3
- 239000008280 blood Substances 0.000 description 2
- 210000004369 blood Anatomy 0.000 description 2
- 210000004204 blood vessel Anatomy 0.000 description 2
- 238000011156 evaluation Methods 0.000 description 2
- 238000002474 experimental method Methods 0.000 description 2
- 230000008569 process Effects 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 238000005070 sampling Methods 0.000 description 2
- 206010011409 Cross infection Diseases 0.000 description 1
- 206010029803 Nosocomial infection Diseases 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000036772 blood pressure Effects 0.000 description 1
- 238000007405 data analysis Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000000691 measurement method Methods 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 239000000758 substrate Substances 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
Classifications
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/02—Detecting, measuring or recording pulse, heart rate, blood pressure or blood flow; Combined pulse/heart-rate/blood pressure determination; Evaluating a cardiovascular condition not otherwise provided for, e.g. using combinations of techniques provided for in this group with electrocardiography or electroauscultation; Heart catheters for measuring blood pressure
- A61B5/021—Measuring pressure in heart or blood vessels
- A61B5/02108—Measuring pressure in heart or blood vessels from analysis of pulse wave characteristics
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
- A61B5/7203—Signal processing specially adapted for physiological signals or for diagnostic purposes for noise prevention, reduction or removal
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
- A61B5/7235—Details of waveform analysis
- A61B5/7253—Details of waveform analysis characterised by using transforms
- A61B5/7257—Details of waveform analysis characterised by using transforms using Fourier transforms
Landscapes
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Molecular Biology (AREA)
- Surgery (AREA)
- Biophysics (AREA)
- Pathology (AREA)
- Signal Processing (AREA)
- Biomedical Technology (AREA)
- Heart & Thoracic Surgery (AREA)
- Medical Informatics (AREA)
- Veterinary Medicine (AREA)
- Physiology (AREA)
- Animal Behavior & Ethology (AREA)
- General Health & Medical Sciences (AREA)
- Public Health (AREA)
- Cardiology (AREA)
- Artificial Intelligence (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Psychiatry (AREA)
- Mathematical Physics (AREA)
- Vascular Medicine (AREA)
- Measuring Pulse, Heart Rate, Blood Pressure Or Blood Flow (AREA)
Abstract
The invention discloses a time-frequency analysis method for processing a photoplethysmography signal, which takes the photoplethysmography signal as an original time domain signal; carrying out fractional Fourier transformation on the signal under different orders, carrying out kurtosis coefficient calculation on the obtained fractional spectrum under each order, carrying out peak detection, and determining the number of peaks as the number of modal decomposition in variable-fraction modal decomposition; constructing a variation modal decomposition algorithm to obtain eigenmode components, calculating a correlation coefficient between each eigenmode component and an original time domain signal, and obtaining a reconstructed signal according to the relation between the eigenmode components and the original time domain signal to realize noise reduction of the signal; and performing generalized wavelet transformation on the reconstructed signal, and performing second-order transient extraction to obtain a two-dimensional time-frequency characteristic representation. The invention can well remove noise interference of the photoelectric volume pulse signals, improve the aggregation degree of time-frequency distribution of the signals, and is well applied to photoelectric volume pulse signal analysis under the background of baseline drift, power frequency noise, electromagnetic interference and the like.
Description
Technical Field
The present invention relates to the field of photoelectric signal processing, and in particular, to a time-frequency analysis method and apparatus for processing a photoplethysmography signal.
Background
The photoelectric volume pulse method is a non-invasive physiological signal measurement technique, and is mainly used for measuring the volume change and pulse waveform of blood vessels. The technology is based on the mechanism of pulse wave, places a photodiode and a light source on the skin of a human body, and then measures the change of the blood volume in a blood vessel through a sensor. The signals generated by the change of the blood volume are converted into electric signals, and related physiological parameters such as blood pressure, heart rate, heartbeat and the like can be obtained through data processing and analysis. The photoelectric volume pulse method has the advantages of no wound, convenient detection, simple operation, stable performance, good repeatability, safety, no cross infection and the like, so that the photoelectric volume pulse method can be used for clinical detection, monitoring and first-aid physical ability test in hospitals, and can also be applied to community and household medical care.
Because the human body is a complex system and the photoelectric signal is a relatively weak signal, the small change in the human body or the weak interference from the outside can influence the waveform of the pulse wave finally acquired, and the influence of baseline drift, power frequency noise, electromagnetic interference and the like is easy to occur in the process of acquiring the volume pulse wave.
Disclosure of Invention
In order to solve the above problems, the present invention provides a time-frequency analysis method and apparatus for processing photoplethysmography signals, wherein the method includes the following steps:
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, combining the reconstructed signal with generalized wavelet transformation and second-order transient extraction transformation, namely second-order transient extraction generalized wavelet transformation (STEGWT) to obtain two-dimensional time-frequency characteristic representation.
Further, the light source in step S1 is 650nm red light or 940nm near infrared light.
Further, in step S2, the fractional fourier transform performed on the original time domain signal under different orders is expressed as:
（1）
wherein ,representing a p-order fractional Fourier transform spectrum, p being a rotation order, n being any integer, alpha being the rotation angle of the time-frequency plane, the relation to the rotation order p being +.>2，/>The amplitude of the fractional Fourier transform result, u is the fractional Fourier transform domain, which is simply called the fractional domain,>j is an imaginary unit for the original time domain signal;
based on the above formula, the signals are calculated at different ordersThe fractional Fourier transform under the condition to obtain a fractional spectrum，/>Representing the square of the absolute value.
Further, in step S3, the calculation formula of kurtosis coefficient for the fractional spectrum under each order is:
（2）
wherein ,is the kurtosis coefficient at order p, E represents the desired value operator, < >>Representing a p-order fractional Fourier transform spectrum, < >>Is->Mean value of->Representing the absolute value.
Further, 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 ,for the j-th eigenmode component and the original time domain signal +.>Phase of (2)The off-coefficient, N is the number of sampling points of the signal, < ->Representation->Mean value of->Representing the j-th eigenmode component, +.>Representation->N is equal to the number of correlation coefficients, +.>Representation->Is a mean value of (c).
Further, the step S6 specifically includes:
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 standardized real window with non-negative symmetry is defined, j is the imaginary unit, ++>The calculation formula is as follows:
（6）
wherein ,representing the reconstructed signal, < >>Time variable representing the reconstructed signal, +.>Represents an intermediate variable,/-> and />Respectively a frequency translation operator and a frequency rotation operator, j is an imaginary unit, m represents the total number of sine functions or cosine functions,> and />For 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 ,the impact component representing the strong pulse component is represented by the dirac delta function,/i>Representing the dirac delta function, a being the impact component amplitude,/->Is the moment at which the impact occurs;
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, i is an imaginary unit,/->For the sign of the derivative of ω, re () represents the real part;
substituting formula (8) into formula (9) includes:
（10）
establishing a second-order frequency-dependent model:
（11）
wherein ,representing the second orderFrequency-variant model->Representing the magnitude of the second order frequency-dependent model,for the phase of the second order frequency-variant model ω represents the frequency in the time-frequency domain of the reconstructed signal,/->Representing the phase of the original time domain signal in the frequency domain, i being the imaginary unit, ζ representing the frequency domain of the 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 a generalized wavelet transform in the reconstructed signal frequency domain under a second order frequency-variant model,representing a second order frequency-dependent model +.>Representing a window function in the frequency domain, i being the 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 ,symbols that are derived from t;
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 for second-order frequency-variant modelLate estimation, ->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 the dirac delta function, i being the 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.
The invention also proposes a time-frequency analysis device for processing photoplethysmography signals, said 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.
The technical scheme provided by the invention has the beneficial effects that:
the method comprises the steps of denoising and time-frequency analysis of photoelectric volume pulse signals, wherein in the implementation process of FrVMD (FrVMD algorithm is formed by combining FrFT technology and VMD), peak detection is carried out by selecting the kurtosis coefficient of FrFT (fractional Fourier transform) under different rotation orders, the obtained peak number is taken as the modal decomposition number of VMD (variable modal decomposition) to carry out modal decomposition and reconstruction on the photoelectric volume pulse signals, and the photoelectric volume pulse signals can be decomposed and reconstructed to remove false and nonsensical eigen mode components, so that interference noise is filtered to a certain extent, and generalized wavelet transform (GWT: generalized Warblet transform) is selected to obtain time-frequency joint distribution of the signals; the Second order transient extraction transformation (STET) is applied to the optimization of the frequency spectrum at the GWT so as to improve the aggregation degree of the time-frequency distribution and greatly reduce the influence of noise interference on the time-frequency representation of the signal. The invention can well remove noise interference of the photoplethysmogram data, improve the aggregation degree of the time-frequency distribution of the photoplethysmogram data, and is well applied to the photoplethysmogram data analysis under the background of baseline drift, power frequency noise, electromagnetic interference and the like.
Drawings
FIG. 1 is a block flow diagram of a time-frequency analysis method for processing photoplethysmography signals according to an embodiment of the present invention;
FIG. 2 is a FrVMD decomposition of the photoplethysmogram 1 according to an embodiment of the present invention;
fig. 3 is a post-FrVMD decomposition reconstruction result of the photoplethysmogram 1 according to an embodiment of the present invention, and fig. 3 (a) is a post-reconstruction time domain signal 1; fig. 3 (b) is an enlarged version at the rectangular box of the reconstructed time domain signal 1; fig. 3 (c) is an original time domain signal 1; fig. 3 (d) is an enlarged version at the rectangular box of the original time domain signal 1;
fig. 4 is a FrVMD decomposed reconstruction result of the photoplethysmogram 1 according to an embodiment of the present invention, and (a) in fig. 4 is a spectrum of the reconstructed signal 1; fig. 4 (b) is an enlarged version at the spectrum rectangular box of the reconstructed signal 1; fig. 4 (c) is the spectrum of the original signal 1; fig. 4 (d) is an enlarged version at the spectrum rectangular box of the original signal 1;
FIG. 5 is a FrVMD decomposition of the photoplethysmogram 2 according to an embodiment of the present invention;
fig. 6 is a post-FrVMD decomposition reconstruction result of the photoplethysmogram 2 according to an embodiment of the present invention, where (a) in fig. 6 is the post-reconstruction time domain signal 2; fig. 6 (b) is an enlarged version at the rectangular box of the reconstructed time domain signal 2; fig. 6 (c) is the original time domain signal 2; fig. 6 (d) is an enlarged version at the rectangular box of the original time domain signal 2;
fig. 7 is a FrVMD decomposition reconstruction result of the photoplethysmogram 2 according to an embodiment of the present invention, and (a) in fig. 7 is a spectrum of the reconstructed signal 2; fig. 7 (b) is an enlarged version at the spectrum rectangular box of the reconstructed signal 2; fig. 7 (c) is the spectrum of the original signal 2; fig. 7 (d) is an enlarged version at the spectrum rectangular box of the original signal 2.
Detailed Description
For the purpose of making the objects, technical solutions and advantages of the present invention more apparent, embodiments of the present invention will be further described with reference to the accompanying drawings.
The invention effectively combines the fractional Fourier transform (FrFT) technology and the Variation Modal Decomposition (VMD) to form a FrVMD algorithm, realizes the denoising treatment of the photo-capacitive product pulse signal, effectively combines the generalized wavelet transform (GWT: generalized Warblet transform) technology and the Second order transient-extraction transform (STET: second-order transient-extracting transform) to form a STEGWT algorithm, and realizes the time-frequency analysis of the photo-capacitive product pulse signal. The signals are decomposed through the FrVMD, whether the modal components are left or not is determined according to whether the correlation between the modal components and the original signals exceeds a threshold value, and finally the left modal components are overlapped to obtain a reconstructed photoplethysmogram signal so as to achieve the purpose of denoising. The method comprises the steps of processing signals by a GWT method to obtain generalized wavelet transformation time-frequency distribution of the signals, and then carrying out second-order transient extraction by STET technology based on the processing result of the GWT to obtain time-frequency distribution with high aggregation.
A flow chart of a time-frequency analysis method for processing photoplethysmography signals according to an embodiment of the present invention is shown in FIG. 1, and includes the following steps:
s1, vertically irradiating the finger of the human body by using a light source, and receiving the transmitted light by using a photoelectric detector to obtain a photoelectric volume pulse wave signal serving as an original time domain signal.
The light source in this embodiment adopts 650nm red light or 940nm near infrared light.
S2, carrying out fractional Fourier transform (FrFT) on the original time domain signal under different orders to obtain fractional spectrums under different orders. Determining the search step size of the FrFT, at [0,2]Obtaining according to the step lengthObtain N rotation orderspThe method comprises the steps of carrying out a first treatment on the surface of the Calculating the signal at different orderspThe fractional Fourier transform under the condition to obtain a fractional spectrum。
The fractional Fourier transform of the original time domain signal under different orders is expressed 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, ">J is an imaginary unit for the original time domain signal.
Based on the above formula, the signals are calculated at different ordersThe fractional Fourier transform under the condition to obtain a fractional spectrum，/>Representing the square of the absolute value.
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 the input parameters in the variable modal decomposition algorithm.
The calculation formula of kurtosis coefficient for the fractional spectrum under each order is as follows:
（2）
wherein ,is thatpThe kurtosis coefficient at the level of the order,Erepresenting the expected value operator, ++>Representation ofpOrder fractional order Fourier transform spectrum, ++>Is->Mean value of->Representing the absolute value.
Carrying out kurtosis coefficient calculation on the fractional spectrum under each order to obtain N kurtosis coefficients; and carrying out peak detection on the kurtosis coefficient, determining the number of peaks, and taking the number of peaks as the K value of the modal decomposition number of the input parameters in the VMD 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 the 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. The reconstruction method can determine which eigenmode components are real components of the signal and which are spurious, nonsensical eigenmode components. Heavy weightThe constructed signal is recorded asx(τ)。
The correlation coefficient between each eigenmode component and the original time domain signal and the set threshold TH are formulated as:
（3）
（4）
wherein ,is the firstjThe individual eigenmode components are +.>Is used for the correlation coefficient of (a),Nfor the number of sampling points of the signal, ">Representation->Mean value of->Represent the firstjThe eigenmode components>Representation->Is used for the average value of (a),nequal to the number of correlation coefficients>Representation->Is a mean value of (c).
S6, combining the reconstructed signal with generalized wavelet transformation and second-order transient extraction transformation to obtain two-dimensional time-frequency characteristic representation.
The reconstructed signal is expressed as follows by generalized wavelet transform:
（5）
wherein ,a generalized wavelet transform representing the reconstructed signal,ωrepresenting the frequency in the time-frequency domain of the reconstructed signal, < >>Is a sliding window function, t represents the time of window center when the time window slides, +.>A normalized real window, typically a gaussian window,jis imaginary unit, ++>The calculation formula is as follows:
（6）
wherein ,representing the reconstructed signal, < >>Time variable representing the reconstructed signal, +.>Represents an intermediate variable,/-> and />A frequency translation operator and a frequency rotation operator respectively,mrepresenting a sineThe total number of functions or cosine functions, and />For Fourier coefficients +.>For the corresponding harmonic component frequencies;
will beIs abbreviated as->，
The impact component of the strong pulse component is expressed as a dirac delta function:
（7）
wherein s (t) represents dirac for the impact component of the strong pulse componentδThe representation of the function, i.e. the impact function of the pulse,representing diracδThe function of the function is that,Ａfor the impact component amplitude, +.>Is the moment at which the impact occurs;
substituting the impact property represented by the formula (7) into the formula (5)The impact characteristic representing the signal is then integrated to give:
（8）
wherein at this timeDelta represents the time support range of the window function, < +.>Is a window function;
ideally, the GWT transformation results of the Dirac delta function should be focused at the moment of impact. Estimating the diffuse energy of a generalized wavelet transform that redistributes reconstructed signals using two-dimensional (2D) Group Delay (GD), written as:
（9）
wherein ,representing a first order two-dimensional group delay,iis imaginary unit, ++>For the sign of the derivative of ω, re () represents the real part;
substituting formula (8) into formula (9) includes:
（10）
equation (10) shows the two-dimensional GD estimate of the Dirac delta function calculated using equation (7) and the time instant of the occurrence of the pulseAnd consistent.
Establishing a second-order frequency-dependent model:
（11）
wherein ,second order generalized wavelet transform representing reconstructed signal, i.e. second order frequency-variant model, +.>Representing the magnitude of the second order frequency-dependent model, +.>For the phase of the second order frequency-dependent model, +.>Representing the phase of the original time domain signal in the frequency domain,iis an imaginary unit, ζ represents the frequency domain of the 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 a generalized wavelet transform in the reconstructed signal frequency domain under a second order frequency-variant model,representing a second order frequency-dependent model +.>Representing a window function in the frequency domain,iis imaginary unit, ++>，，/>Represents the standard deviation of the window function +.>Is->Is obtained by substituting the formula (11) into the formula (12):
（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）
equation (15) provides an indication of how to calculate a new 2D GD estimate in the second order signal analysis. Thus, a new two-dimensional estimate is defined：
（16）
wherein ,to get a pair oftSymbols of partial derivatives;
the new two-dimensional group delay calculation method defined according to the formula (16) takes the correlation constraint into consideration, and the second-order two-dimensional group delay of the second-order frequency-variant model is expressed as follows:
（17）
wherein ,namely, second-order two-dimensional group delay estimation of second-order frequency-variant model,/->A first-order two-dimensional group delay estimation representing a second-order frequency-dependent model;
the second order transient transformation of equation (11) is:
（18）
thus, the first and second substrates are bonded together,the result of the generalized wavelet transform is extracted for the FrVMD-second order transient state, and a high-resolution time-frequency representation can be generated for the second-order frequency-varying signal.
The algorithm can well remove the noise of the signal, obtain the time-frequency joint distribution condition of the signal, help to further understand the time-frequency distribution change of the photoelectric volume pulse data under the background of baseline drift, power frequency noise, electromagnetic interference and the like, and provide basis and help for the recognition and research of the photoelectric volume pulse signal. In order to verify the effectiveness of the algorithm in photoelectric volume pulse data processing research under the background of baseline drift, power frequency noise, electromagnetic interference and the like, two groups of measured data 1 containing noise and measured data 2 containing baseline drift are selected for experimental verification.
1. FrVMD denoising experiment
(1) Noisy measured data 1
First, the noise reduction processing is performed on the actual measurement data 1 containing noise by using the FrVMD, and fig. 2 shows the decomposition result of the FrVMD, and table 1 shows the correlation between each modal component and the original signal and the threshold value.
TABLE 1
As can be seen from fig. 2 and table 1, the photoplethysmogram signal is decomposed into 3 modal components, and only modal component 1 exceeds the threshold, i.e. has a high correlation with the original signal, modal component 1 is considered to be an effective component, and modal component 2 and modal component 3 are considered to be ineffective components to be removed. And then summing the left modal components to obtain a final reconstruction result.
Fig. 3 is a post-FrVMD decomposition reconstruction result of the photoplethysmogram 1 according to an embodiment of the present invention, and fig. 3 (a) is a post-reconstruction time domain signal 1; fig. 3 (b) is an enlarged version at the rectangular box of the reconstructed time domain signal 1; fig. 3 (c) is an original time domain signal 1; fig. 3 (d) is an enlarged version of the original time domain signal 1 at the rectangular frame, fig. 4 is a FrVMD decomposed reconstruction result of the photoplethysmogram signal 1 according to the embodiment of the present invention, and fig. 4 (a) is a spectrum of the reconstructed signal 1; fig. 4 (b) is an enlarged version at the spectrum rectangular box of the reconstructed signal 1; fig. 4 (c) is the spectrum of the original signal 1; fig. 4 (d) is an enlarged version at the spectrum rectangular box of the original signal 1.
As can be seen from fig. 3, the original time domain signal contains more noise and is more glitch in the appearance of the time domain curve. The FrVMD method can effectively remove noise of signals, so that the reconstructed time domain signal curve is smoother and free of burrs. As can be seen from fig. 4, the spectrum of the reconstructed signal does not contain the unwanted noise spectrum, whereas the original signal obviously contains the unwanted noise spectrum in addition to the spectrum of the real component at the center. The FrVMD method can effectively remove noise of a signal.
In order to further compare the performance of the FrVMD reconstruction method, the average absolute percentage error (MAPE), root Mean Square Error (RMSE), and determination coefficient (R2) are selected to determine the decomposition reconstruction denoising effect, and the definition is shown in the formulas (19), (20), and (21). The larger the Mean Absolute Error (MAE) and Root Mean Square Error (RMSE), the less desirable the denoising effect of the signal. The closer the determination coefficient (R2) is to 1, the better the denoising effect.
（19）
（20）
（21）
The different error values of the signals are shown in table 2, and as can be seen from table 2, the error of the FrVMD algorithm in decomposing and reconstructing the photoplethysmogram signal is smaller, the reconstructed signal has high correlation with the original signal, and the FrVMD algorithm has a relatively obvious effect of removing noise on the basis of retaining the effective and real components of the original signal.
TABLE 2
(2) Measured data 2 containing baseline wander
The measured data 2 including baseline drift was subjected to noise reduction processing by using the FrVMD, and fig. 5 shows the decomposition result of the FrVMD, and fig. 3 shows the correlation between each modal component and the original signal and the threshold value. As can be seen from fig. 5 and table 3, the photoplethysmogram signal is decomposed into 10 modal components, and only modal component 1 and modal component 2 exceed the threshold, i.e. have high correlation with the original signal, the modal component 1 and modal component 2 are considered to be effective components, and the remaining modal components are ineffective components to be removed. And then summing the left modal components to obtain a final reconstruction result.
TABLE 3 Table 3
Fig. 6 is a post-FrVMD decomposition reconstruction result of the photoplethysmogram 2 according to an embodiment of the present invention, where (a) in fig. 6 is the post-reconstruction time domain signal 2; fig. 6 (b) is an enlarged version at the rectangular box of the reconstructed time domain signal 2; fig. 6 (c) is the original time domain signal 2; fig. 6 (d) is an enlarged version at the rectangular box of the original time domain signal 2.
Fig. 7 is a FrVMD decomposition reconstruction result of the photoplethysmogram 2 according to an embodiment of the present invention, and (a) in fig. 7 is a spectrum of the reconstructed signal 2; fig. 7 (b) is an enlarged version at the spectrum rectangular box of the reconstructed signal 2; fig. 7 (c) is the spectrum of the original signal 2; fig. 7 (d) is an enlarged version at the spectrum rectangular box of the original signal 2.
As can be seen from fig. 6, the original time domain signal contains baseline wander, which is shown in fig. 6 (d) in the time domain curve. The FrVMD method can effectively remove baseline drift of the signal, so that the reconstructed time domain signal curve is more regular, as shown in (b) in fig. 6. As can be seen from fig. 7, the spectrum of the reconstructed signal does not contain the unwanted noise spectrum, whereas the original signal obviously contains the unwanted noise spectrum in addition to the spectrum of the real component at the center. The FrVMD method can effectively remove noise of a signal.
Table 4 shows that the error of the FrVMD algorithm in decomposing and reconstructing the photoplethysmogram signal is smaller, and the reconstructed signal has higher correlation with the original signal, which indicates that the FrVMD algorithm has obvious effect of removing the interference such as baseline drift on the basis of retaining the effective and real components of the original signal.
TABLE 4 Table 4
2. STEGWT time-frequency analysis experiment
(1) Noisy measured data 1
Firstly, STEGWT is adopted to carry out time-frequency analysis on the noise-containing measured data 1, the time-frequency curve obtained by the FrVMD-STEGWT method is more concentrated, and the time-frequency aggregation degree is higher.
The higher the time-frequency aggregation degree is, the more accurate the time-frequency distribution is for the local positioning and resolution of the signal frequency components, so that the better the processing effect of the time-frequency analysis method is. Herein, rayleigh Li Shang (renyi entopy) is taken as a performance evaluation index for measuring the aggregation degree of time-frequency distribution, and the performance evaluation index is calculatedThe formula is:, wherein ,/>Representing the number of rayleigh entropies, taken here +.>；/>Representing the time-frequency distribution of the signal. Rayleigh entropy->The smaller the value of (c) is, the higher the aggregation degree of the time-frequency distribution is. />
Table 5 shows that the Rayleigh entropy value of the FrVMD-STEGWT of the signal 1 is 9.6842, which is far smaller than that of the 3 time-frequency analysis methods of SST, GWT, WVD, namely the FrVMD-STEGWT has higher time-frequency aggregation degree, and the improved method can effectively improve the time-frequency aggregation degree when processing photoelectric volume pulse signals containing noise interference.
TABLE 5
(2) Measured data 2 containing baseline wander
The FrVMD-STEGWT method has the advantages of concentrated time-frequency curve and high time-frequency aggregation degree.
Table 6 shows the rayleigh entropy of different time-frequency distributions of the signal 2, and the higher the time-frequency aggregation degree, the more accurate the local positioning and resolution of the time-frequency distribution on the frequency components of the signal, thereby indicating that the better the processing effect of the time-frequency analysis method. The Rayleigh entropy value of the FrVMD-STEGWT is 11.3464, which is far smaller than SST, GWT, WVD, namely the time-frequency aggregation degree of the FrVMD-STEGWT is higher, and the improvement method can effectively improve the time-frequency aggregation degree of the FrVMD-STEGWT when the photoelectric volume pulse signal with baseline drift is processed.
TABLE 6
The invention also proposes a time-frequency analysis device for processing photoplethysmography signals, 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.
The previous description of the disclosed embodiments is provided to enable any person skilled in the art to make or use the present invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the generic principles defined herein may be applied to other embodiments without departing from the spirit or scope of the invention. Thus, the present invention is not intended to be limited to the embodiments shown herein but is to be accorded 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 value of->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),nequal 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, and />A frequency translation operator and a frequency rotation operator respectively,jin units of imaginary numbers,mrepresenting the total number of sine functions or cosine functions, +.> and />For 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, +.>Is the moment at which the impact occurs;
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, +.>Being a second-order frequency-dependent modelThe phase of the light emitted from the light source is,ω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 as a second order frequency-dependent modelTwo-dimensional group delay estimation,/->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.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311054338.0A CN116763275B (en) | 2023-08-22 | 2023-08-22 | Time-frequency analysis method and device for processing photoplethysmography signals |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311054338.0A CN116763275B (en) | 2023-08-22 | 2023-08-22 | Time-frequency analysis method and device for processing photoplethysmography signals |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116763275A CN116763275A (en) | 2023-09-19 |
| CN116763275B true CN116763275B (en) | 2023-10-31 |
Family
ID=88011972
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202311054338.0A Active CN116763275B (en) | 2023-08-22 | 2023-08-22 | Time-frequency analysis method and device for processing photoplethysmography signals |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116763275B (en) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117473233B (en) * | 2023-12-27 | 2024-03-01 | 华东交通大学 | Sample component analysis method, system, storage medium and computer |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110309817A (en) * | 2019-07-19 | 2019-10-08 | 北京理工大学 | A kind of pulse wave motion artifacts minimizing technology of parameter adaptive optimization VMD |
| CN113627375A (en) * | 2021-08-16 | 2021-11-09 | 北京信息科技大学 | Planetary gear fault diagnosis method and system, storage medium and computing device |
| CN116449328A (en) * | 2023-04-24 | 2023-07-18 | 中国人民解放军陆军工程大学 | Radar working mode identification method and system based on time-frequency analysis feature fusion |
| CN116522073A (en) * | 2023-04-10 | 2023-08-01 | 哈尔滨工业大学 | Adaptive signal decomposition method of fractional Fourier transform domain, storage medium and computer |
-
2023
- 2023-08-22 CN CN202311054338.0A patent/CN116763275B/en active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110309817A (en) * | 2019-07-19 | 2019-10-08 | 北京理工大学 | A kind of pulse wave motion artifacts minimizing technology of parameter adaptive optimization VMD |
| CN113627375A (en) * | 2021-08-16 | 2021-11-09 | 北京信息科技大学 | Planetary gear fault diagnosis method and system, storage medium and computing device |
| CN116522073A (en) * | 2023-04-10 | 2023-08-01 | 哈尔滨工业大学 | Adaptive signal decomposition method of fractional Fourier transform domain, storage medium and computer |
| CN116449328A (en) * | 2023-04-24 | 2023-07-18 | 中国人民解放军陆军工程大学 | Radar working mode identification method and system based on time-frequency analysis feature fusion |
Non-Patent Citations (2)
| Title |
|---|
| A novel solution for improved performance of Time-frequency concentration;Juan Guo等;Mechanical Systems and Signal Processing;第1-23页 * |
| SE-KSD优化的FRFT-VMD轴承故障诊断方法;韦明辉;机械科学与技术;第1-9页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116763275A (en) | 2023-09-19 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Wang et al. | Adaptive Fourier decomposition based ECG denoising | |
| Sathyapriya et al. | Analysis and detection R-peak detection using Modified Pan-Tompkins algorithm | |
| Bahoura et al. | DSP implementation of wavelet transform for real time ECG wave forms detection and heart rate analysis | |
| Nagendra et al. | Application of wavelet techniques in ECG signal processing: an overview | |
| Petrenas et al. | An echo state neural network for QRST cancellation during atrial fibrillation | |
| Castillo et al. | Noise suppression in ECG signals through efficient one-step wavelet processing techniques | |
| CN116763275B (en) | Time-frequency analysis method and device for processing photoplethysmography signals | |
| TW201717845A (en) | Devices, systems, and methods for determining heart rate of a subject from noisy electrocardiogram data | |
| Awodeyi et al. | Median filter approach for removal of baseline wander in photoplethysmography signals | |
| Luengo et al. | Blind analysis of atrial fibrillation electrograms: a sparsity-aware formulation | |
| CN114027804A (en) | Pulse condition diagnosis method, device and readable storage medium | |
| Bhogeshwar et al. | To verify and compare denoising of ECG signal using various denoising algorithms of IIR and FIR filters | |
| Sahoo et al. | Autocorrelation and Hilbert transform-based QRS complex detection in ECG signal | |
| Barmase et al. | Wavelet transform-based analysis of QRS complex in ECG signals | |
| Tajane et al. | Comparative analysis of mother wavelet functions with the ecg signals | |
| Wu et al. | An improved method for ECG signal feature point detection based on wavelet transform | |
| Sikkandar et al. | Analysis of cardiac abnormalities using Hilbert-Huang transform | |
| Bognár et al. | Geometric interpretation of QRS complexes in ECG signals by rational functions | |
| Zhang et al. | Design of a real-time ECG filter for resource constraint computer | |
| Hung et al. | A method for suppressing respiratory noise in impedance cardiography and comprehensive assessment of noise reduction performance | |
| Wang et al. | A novel frequency domain method for estimating blood pressure from photoplethysmogram | |
| Zhou et al. | Fetal Electrocardiogram Extraction and Performance Analysis. | |
| Yanık et al. | Detection of ECG characteristic points using multiresolution analysis | |
| Arvinti et al. | Adaptive thresholding algorithm for noisy electrocardiograms using reverse biorthogonal mother wavelets | |
| Ouali et al. | Electrocardiogram Signal Denoising by Hilbert Transform and Synchronous Detection |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |