PPG signal decomposition method and device based on downsampling short-time variational modal decomposition

Abstract

The present invention discloses a PPG signal decomposition method based on down-sampling short-time variational modal decomposition, comprising: step 1, selecting a function window, multiplying the function window with the acquired signal at each sampling moment, and obtaining the signal frame of the signal at each sampling moment; step 2, setting the initial value of the instantaneous frequency; step 3, according to the corresponding initial value of the instantaneous frequency, using the variational modal decomposition method to decompose the signal frame, and obtain each signal component fragment corresponding to the signal frame; step 4, overlapping and adding each signal component fragment to obtain each complete signal component and the corresponding instantaneous frequency; step 5, according to the instantaneous frequency corresponding to the signal component at each sampling moment, reconstructing the instantaneous frequency-time function corresponding to each signal component in the PPG signal. The present invention also provides a PPG signal decomposition device. The method of the present invention does not need to rely on the ability of RD, and can process PPG signals with severe spectrum aliasing to obtain pure PPG signals.

Description

Technical Field

[0001] The present invention relates to the technical field of electrical signal decomposition and reconstruction, and in particular to a PPG signal decomposition method and device based on down-sampling short-time variational modal decomposition. Background Art

[0002] Photoplethysmography (PPG) is a pulse signal collected by photoplethysmography. Specifically, this method uses a sensor to emit low-intensity infrared light to the skin. The light passes through the skin tissue, arteries, and veins, and is eventually absorbed and reflected into the sensor. Since the absorption rate of light by tissues such as skin, fat, and muscle remains basically unchanged during blood circulation, and the blood volume of subcutaneous blood vessels is very sensitive to the absorption rate of light, the intensity of the reflected light received by the sensor can reflect the changes in blood flow in the subcutaneous blood vessels of the human body over time, thereby providing information about the physiological state of the human body.

[0003] At present, PPG signals are widely used in sports or health detection functions of wrist-worn wearable devices. It can provide important physiological indicators to measure health status, such as heart rate, blood oxygen, etc. However, due to human movement factors, the PPG signal collected by wrist-worn wearable devices is easily interfered by motion artifacts. Motion artifacts and PPG signals are often in similar frequency bands, and it is almost impossible to distinguish the two signal components in the frequency domain. Therefore, it is difficult to remove motion artifacts and restore clean PPG signals through frequency domain signal processing methods. On the other hand, both PPG signals and motion artifact signals are nonlinear, non-stationary time-varying signals, and their dynamic characteristics over time can be expressed in the time-frequency domain. As a result, the two can be more clearly distinguished in the time-frequency plane.

[0004] Patent document CN106691425A discloses a wrist heart rate monitoring method for a sports bracelet. The method first performs threshold processing on the original optical pulse wave signal to eliminate signals unrelated to the heart rate. Then, the synchronous acceleration signal is processed by a singular spectrum analysis method to obtain a motion artifact reference signal. The motion artifact signal and the pre-processed optical signal are adaptively filtered to eliminate the motion artifact components in the optical signal to obtain a pure pulse wave signal. This method processes the optical signal by generating a reference signal, but there may be a problem of signal loss when generating the reference signal.

[0005] Patent document CN113367676A discloses a PPG motion noise filtering method and device for wearable devices, including: performing fast Fourier transform on the PPG signal of photoplethysmography to obtain the PPG signal spectrum, and correcting the PPG signal spectrum; performing inverse Fourier transform on the corrected PPG signal spectrum to obtain a noise reference signal; and using an adaptive filtering algorithm to adaptively filter the noise reference signal to obtain a filtered PPG signal. This method relies on the ability of RD to complete the signal separation work. Summary of the invention

[0006] In order to solve the above problems, the present invention provides a PPG signal decomposition method based on down-sampling short-time variational modal decomposition. The method has strong noise robustness, does not need to rely on the ability of RD, can process PPG signals with severe spectral aliasing, and obtain pure PPG signals.

[0007] A PPG signal decomposition method based on down-sampling short-time variational mode decomposition, comprising:

[0008] Step 1, selecting a function window, multiplying the function window with the acquired signal at each sampling time, and obtaining a signal frame of the signal at each sampling time;

[0009] Step 2: Initialize the signal frame obtained in step 1 according to the preset mode number to obtain the corresponding instantaneous frequency initial value;

[0010] Step 3: Decompose the signal frame using the variational mode decomposition method according to the initial value of the instantaneous frequency obtained in step 2, and set the initial estimated value of the instantaneous frequency for the signal component in the next signal frame, repeat the operation until the window function reaches the end of the signal, and obtain the signal component fragments corresponding to the signal frame;

[0011] Step 4: Based on the time axis order, overlap and add the signal component fragments obtained in step 3 to obtain each complete signal component and the corresponding instantaneous frequency;

[0012] Step 5: Reconstruct the instantaneous frequency-time function corresponding to each signal component in the PPG signal according to the instantaneous frequency corresponding to the signal component at each sampling moment in step 4.

[0013] The present invention is based on the framework of short-time Fourier transform (STVMD), performs instantaneous frequency extraction and pattern decomposition in a time window manner, and uses the frequency spectrum of each signal component and the estimated value of the instantaneous frequency to guide iteration to obtain the final true value, and then reconstructs the instantaneous frequency-time function corresponding to each signal component by obtaining the short-time inverse time-frequency transform.

[0014] Specifically, the signal in step 1 is a PPG signal, which is acquired by photoplethysmography and is used to evaluate physiological indicators of human movement health status.

[0015] Preferably, the expression of the signal frame at each sampling moment in step 1 is as follows:

[0016] x t [n]=x[n]g[nt]

[0017] Among them, x represents the signal, g represents the window function, x t Represents the signal frame obtained by multiplying the signal by the window function at sampling time t.

[0018] Specifically, the initial value of the instantaneous frequency in step 2 is selected to be 0.

[0019] Specifically, the specific process of step 3 is as follows:

[0020] Step 3-1, performing discrete Fourier transform on the signal frame to obtain the corresponding spectrum;

[0021] Step 3-2, according to the spectrum obtained in step 3-1, setting the spectrum estimation value of the signal component segment, and updating the spectrum estimation value of the signal component segment through iterative calculation;

[0022] Step 3-3, updating the instantaneous frequency estimation value corresponding to the signal component segment through iterative calculation according to the frequency spectrum obtained in step 3-1;

[0023] Step 3-4, repeat steps 3-2 to 3-3 until the spectrum estimation value of the signal component segment meets the termination condition, and obtain the spectrum and instantaneous frequency of the signal component segment.

[0024] Preferably, the step 3-2 uses a spectrum iteration calculation formula to update the spectrum estimation value, and its expression is as follows:

[0025]

[0026] in, Represents signal frame x t The spectrum obtained by discrete Fourier transform of the pth signal component fragment, λ t represents the Lagrange multiplier, α represents the penalty term coefficient, f s represents the sampling frequency of the signal, represents the instantaneous frequency of the pth signal component fragment at the sampling time t, and N represents the length of the window function g.

[0027] Preferably, the step 3-3 updates the instantaneous frequency estimation value by using an iterative calculation formula according to the center frequency, and the expression is as follows:

[0028]

[0029] in, represents the instantaneous frequency of the pth signal component fragment at sampling time t, f s represents the sampling frequency of the signal, N represents the length of the window function g, Represents signal frame x t The spectrum obtained by discrete Fourier transform of the pth signal component segment.

[0030] Preferably, the expression of the initial estimated value of the instantaneous frequency in step 3 is as follows:

[0031]

[0032] in, represents the instantaneous frequency of the pth signal component at sampling time t, represents the initial estimate of the instantaneous frequency, It represents the parameter that characterizes the instantaneous frequency change rate.

[0033] Specifically, the calculation process of the instantaneous frequency initial estimate is as follows:

[0034] Step 3-5, calculate the estimated value of the frequency modulation parameter of the signal component in the signal frame at the current moment:

[0035]

[0036] in, represents the estimated value of the frequency modulation parameter of the pth signal component at sampling time t, represents the time-frequency of the short-time Fourier transform, g′, g″, represents the window function coefficient;

[0037] Step 3-6: Calculate the bandwidth estimate of the signal component in the signal frame at the current moment:

[0038]

[0039] in, represents the bandwidth estimate of the pth signal component at sampling time t, represents the instantaneous frequency of the pth signal component at sampling time t, f s Indicates the signal sampling frequency;

[0040] Step 3-7: Calculate the parameter that characterizes the instantaneous frequency change rate:

[0041]

[0042] in, represents the parameter characterizing the instantaneous frequency change rate of the pth signal component at the sampling time t;

[0043] Step 3-8: Set the initial estimated value of the instantaneous frequency according to the parameters obtained from step 3-8 to step 3-7.

[0044] Specifically, the complete signal component in step 4 is expressed as follows:

[0045]

[0046] Among them, x p represents the pth signal component, Represents the pth signal component fragment of the signal frame.

[0047] The present invention also provides a PPG signal decomposition device, comprising a computer memory, a computer processor, and a computer program stored in the computer memory and executable on the computer processor, wherein the computer memory executes the above-mentioned PPG signal decomposition method based on down-sampling short-time variational modal decomposition; when the computer processor executes the computer program, the following steps are implemented: inputting an initial PPG signal acquired, analyzing and reconstructing it according to the PPG signal decomposition method, and outputting a signal amplitude diagram corresponding to each signal component in the PPG signal.

[0048] Compared with the prior art, the present invention has the following beneficial effects:

[0049] (1) Using the short-time Fourier transform framework, the PPG signal is separated and the cross components are restored in a time-window manner. The separation of PPG signals and motion artifacts can be achieved without relying on the ability of RD.

[0050] (2) When there is noise in the signal, the method of the present invention can still completely separate the different signal components and obtain a pure PPG signal amplitude diagram. BRIEF DESCRIPTION OF THE DRAWINGS

[0051] Figure 1 The original PPG signal diagram provided for this embodiment;

[0052] Figure 2 A real motion artifact signal diagram provided for this embodiment;

[0053] Figure 3 A pure PPG signal diagram provided for this embodiment;

[0054] Figure 4 The PPG signal decomposition method based on down-sampling short-time variational modal decomposition provided by the present invention;

[0055] Figure 5A signal component diagram obtained by decomposing the PPG signal provided in this embodiment;

[0056] Figure 6 A time-frequency representation diagram obtained by decomposing the PPG signal provided in this embodiment;

[0057] Figure 7 This is a schematic diagram of the output results of the PPG signal decomposition device provided in this embodiment. DETAILED DESCRIPTION

[0058] Taking the problem of removing motion artifacts from PPG signals as an example, the implementation plan and effect of the cross-signal decomposition method based on down-sampling variational mode decomposition are described below. During human motion, the PPG signal collected by the wearable device will be contaminated by motion artifacts, and its instantaneous frequency will obviously cross with the frequency of the PPG signal. This example is carried out in the MATLAB simulation environment:

[0059] like Figure 1 As shown, the original PPG signal image containing motion artifacts is acquired by the wearable device;

[0060] like Figure 2 As shown, it is the real motion artifact signal diagram existing in the original PPG signal diagram;

[0061] like Figure 3 As shown, this is a pure PPG signal diagram after removing motion artifacts using traditional methods.

[0062] like Figure 4 As shown, a PPG signal decomposition method based on down-sampling short-time variational mode decomposition includes:

[0063] Step 1, select a function window g[nt], multiply the function window with the acquired PPG signal at each sampling time, and obtain the signal frame of the PPG signal at each sampling time;

[0064] x t [n]=x[n]g[nt]

[0065] Among them, x represents the signal, g represents the window function, x t Represents the signal frame obtained by multiplying the signal by the window function at sampling time t;

[0066] Step 2: Initialize the signal frame obtained in step 1 according to the preset modal number to obtain the corresponding instantaneous frequency initial value. The instantaneous frequency initial value in this embodiment is selected as 0:

[0067] Step 3: According to the initial value of instantaneous frequency set in step 2, the signal frame is decomposed using the variational mode decomposition method:

[0068] Step 3-1, performing discrete Fourier transform on the signal frame to obtain the corresponding spectrum;

[0069] Step 3-2: according to the spectrum obtained in step 3-1, set the spectrum estimation value of the signal component segment, and update the spectrum estimation value of the signal component segment using the spectrum iteration calculation formula:

[0070]

[0071] in, Represents signal frame x t The spectrum obtained by discrete Fourier transform of the pth signal component fragment, λ t represents the Lagrange multiplier, α represents the penalty term coefficient, f s represents the sampling frequency of the signal, represents the instantaneous frequency of the pth signal component fragment at sampling time t, and N represents the length of the window function g;

[0072] Step 3-3: Based on the spectrum obtained in step 3-1, the center frequency iteration calculation formula is used to update the instantaneous frequency estimate corresponding to the signal component segment:

[0073]

[0074] in, represents the instantaneous frequency of the pth signal component fragment at sampling time t, f s represents the sampling frequency of the signal, N represents the length of the window function g, Represents signal frame x t The spectrum obtained by discrete Fourier transform of the pth signal component segment;

[0075] Step 3-4, repeating steps 3-2 to 3-3 until the spectrum estimation value of the signal component segment meets the termination condition, and obtaining the spectrum and instantaneous frequency of the signal component segment;

[0076] And set the initial estimated value of instantaneous frequency for the signal component in the next signal frame, repeat the operation until the window function reaches the end of the signal, and obtain the signal component fragments corresponding to the signal frame;

[0077] The calculation process of the initial estimate of the instantaneous frequency is as follows:

[0078] Step 3-5, calculate the estimated value of the frequency modulation parameter of the signal component in the signal frame at the current moment:

[0079]

[0080] in, represents the estimated value of the frequency modulation parameter of the pth signal component at sampling time t, represents the time-frequency of the short-time Fourier transform, g′, g″, represents the window function coefficient;

[0081] Step 3-6: Calculate the bandwidth estimate of the signal component in the signal frame at the current moment:

[0082]

[0083] in, represents the bandwidth estimate of the pth signal component at sampling time t, represents the instantaneous frequency of the pth signal component at sampling time t, f s Indicates the signal sampling frequency;

[0084] Step 3-7: Calculate the parameter that characterizes the instantaneous frequency change rate:

[0085]

[0086] in, represents the parameter characterizing the instantaneous frequency change rate of the pth signal component at the sampling time t;

[0087] Step 3-8: According to the parameters obtained from steps 3-5 to 3-7, set the initial estimate of the instantaneous frequency:

[0088]

[0089] in, represents the instantaneous frequency of the pth signal component at sampling time t, represents the initial estimate of the instantaneous frequency, It represents the parameter that characterizes the instantaneous frequency change rate;

[0090] Step 4: Based on the time axis order, overlap and add the signal component fragments obtained in step 3:

[0091]

[0092] Obtain each complete signal component and the corresponding instantaneous frequency;

[0093] like Figure 5 As shown in the figure, the signal components of the initial PPG signal are shown. The waveform of u2 in the figure is similar to Figure 2 The real motion artifact signal graph in has the highest similarity.

[0094] like Figure 6 As shown, it is the instantaneous frequency diagram corresponding to each signal component of the initial PPG signal.

[0095] Step 5: Reconstruct the instantaneous frequency-time function corresponding to each signal component in the PPG signal according to the instantaneous frequency corresponding to the signal component at each sampling moment in step 4.

[0096] The present invention also provides a PPG signal decomposition device, comprising a computer memory, a computer processor, and a computer program stored in the computer memory and executable on the computer processor, wherein the computer memory executes the above-mentioned PPG signal decomposition method based on downsampled short-time variational modal decomposition.

[0097] like Figure 7 As shown, to remove Figure 5 The PPG signal diagram after u2 corresponds to the signal component Figure 1 The comparison chart of the pure PPG signal graph in the figure shows that according to calculation, the mean square error of the two is only 0.224, which means that the output result of the device is highly consistent with the theoretical result, which can show that the purification effect of the device is good.

Claims

1. A PPG signal decomposition method based on down-sampling short-time variational mode decomposition, characterized in that: include: Step 1: Select a function window, multiply the function window with the acquired signal at each sampling time, and obtain the signal frame of the signal at each sampling time. The expression of the signal frame at each sampling time is as follows: x t [n]=x[n]g[n-t] Among them, x represents the signal, g represents the window function, x t Represents the signal frame obtained by multiplying the signal by the window function at sampling time t; Step 2: Initialize the signal frame obtained in step 1 according to the preset mode number to obtain the corresponding instantaneous frequency initial value; Step 3: According to the initial value of instantaneous frequency obtained in step 2, the signal frame is decomposed by using the variational mode decomposition method, and the initial estimated value of instantaneous frequency is set for the signal component in the next signal frame. The operation is repeated until the window function reaches the end of the signal to obtain the signal component fragments corresponding to the signal frame. The expression of the initial estimated value of instantaneous frequency is as follows: in, represents the instantaneous frequency of the pth signal component at sampling time t, represents the initial estimate of the instantaneous frequency, It represents the parameter that characterizes the instantaneous frequency change rate; The specific process of step 3 is as follows: Step 3-1, performing discrete Fourier transform on the signal frame to obtain the corresponding spectrum; Step 3-2: according to the spectrum obtained in step 3-1, set the spectrum estimation value of the signal component segment, update the spectrum estimation value of the signal component segment by iterative calculation, and update the spectrum estimation value by using the spectrum iterative calculation formula, which is expressed as follows: in, Represents signal frame x t The spectrum obtained by discrete Fourier transform of the pth signal component fragment, λ t represents the Lagrange multiplier, α represents the penalty term coefficient, f s represents the sampling frequency of the signal, represents the instantaneous frequency of the pth signal component fragment at sampling time t, and N represents the length of the window function g; Step 3-3: According to the spectrum obtained in step 3-1, the instantaneous frequency estimation value corresponding to the signal component segment is updated by iterative calculation, and the instantaneous frequency estimation value is updated by using the center frequency iterative calculation formula, which is expressed as follows: in, represents the instantaneous frequency of the pth signal component fragment at sampling time t, f s represents the sampling frequency of the signal, N represents the length of the window function g, Represents signal frame x t The spectrum obtained by discrete Fourier transform of the pth signal component segment; Step 3-4, repeating steps 3-2 to 3-3 until the spectrum estimation value of the signal component segment meets the termination condition, and obtaining the spectrum and instantaneous frequency of the signal component segment; Step 4: Based on the time axis order, overlap and add the signal component fragments obtained in step 3 to obtain each complete signal component and the corresponding instantaneous frequency; Step 5: Reconstruct the instantaneous frequency-time function corresponding to each signal component in the PPG signal according to the instantaneous frequency corresponding to the signal component at each sampling moment in step 4.

2. The PPG signal decomposition method based on down-sampling short-time variational modal decomposition according to claim 1, characterized in that: The signal in step 1 is a PPG signal, which is acquired by photoplethysmography and is used to evaluate the physiological indicators of human movement health status.

3. The PPG signal decomposition method based on down-sampling short-time variational modal decomposition according to claim 1, characterized in that: The complete signal component in step 4 is expressed as follows: Among them, x p represents the pth signal component, Represents the pth signal component fragment of the signal frame.

4. A PPG signal decomposition device, comprising a computer memory, a computer processor, and a computer program stored in the computer memory and executable on the computer processor, characterized in that: The computer memory executes the PPG signal decomposition method based on down-sampled short-time variational modal decomposition as described in any one of claims 1 to 3; the computer processor implements the following steps when executing the computer program: inputting the acquired initial PPG signal, analyzing and reconstructing it according to the PPG signal decomposition method, and outputting a signal amplitude diagram corresponding to each signal component in the PPG signal.

Images