Cardiopulmonary coupling quantification method and system based on multimodal coupling analysis

Abstract

The present invention discloses a cardiopulmonary coupling quantification method and system based on multimodal coupling analysis, which relates to the technical field of electrocardiogram signal processing. The method comprises the following steps: collecting electrocardiogram (ECG) signals and respiratory signals of a subject; extracting a heart beat interval (R-R) interval time series from the collected ECG signals; decomposing the R-R interval time series and the respiratory signal using variational mode decomposition to obtain intrinsic mode functions (IMFs) of the two time series; selecting the IMF with the largest power as the dominant IMF among all IMFs of the respiratory signal; selecting the IMF with the frequency matching the dominant IMF of the respiratory signal among the IMFs of the R-R signal; calculating the synchronization index between the dominant IMF of the respiratory signal and the frequency matching IMF of the R-R sequence by using the instantaneous phase difference between the two; and finally, quantifying the instantaneous magnitude of RSA by calculating the power of the frequency matching IMF of the R-R interval sequence during its strong synchronization period.

Description

Technical Field

[0001] The present invention relates to the technical field of electrocardiogram (ECG) signal processing, and in particular to a cardiopulmonary coupling quantification method and system based on multimodal coupling analysis. Background Art

[0002] The autonomic nervous system (ANS) is a critical system in the body that regulates internal environmental stability, responds to external stimuli, and maintains vital signs. It controls the functions of multiple organs, including the cardiovascular, respiratory, and digestive systems, through the interaction of the sympathetic nervous system (SNS) and the parasympathetic nervous system (PNS). An imbalance in the stability of this system, particularly a decrease in the stability and adaptability of the parasympathetic nervous system, can severely impact health. For example, studies have shown that it increases the risk of myocardial infarction and is a prominent feature of patients with congestive heart failure. Diseases associated with autonomic dysfunction and decreased parasympathetic function can severely impact patients' quality of life. High medical costs impose a significant burden on society and have become a prominent medical and public health issue. Accurately quantifying parasympathetic function is fundamental to alleviating this problem.

[0003] Respiratory sinus arrhythmia (RSA) is characterized by an acceleration of the heart rate with inspiration and a decrease with expiration during breathing, reflecting heart rate variability (HRV) and respiratory synchronization. Assessing RSA provides an attractive and noninvasive method for quantifying parasympathetic function. A standard and widely used quantitative RSA metric in clinical research is the high-frequency (HF) power of the heart rate time series, calculated by spectral analysis of heart rate variability within the normal respiratory frequency band (0.15-0.40 Hz).

[0004] However, the use of high-frequency power to quantify RSA has certain limitations. First, to accommodate normal physiological conditions, the predefined HF range restricts the corresponding respiratory rate to 9 to 24 breaths / minute. This premise cannot be met in some clinical settings. For example, the respiratory rate of patients with congestive heart failure sometimes exceeds the conventional HF band, which can lead to deviations in RSA assessment, which is crucial for their health and well-being. Second, even if the respiratory rate is well controlled within the 0.15-0.40 Hz band, the respiratory signal, like other physiological signals, is a time-varying, non-sinusoidal signal due to the inconsistent duration of each breath and the unequal duration of exhalation and inspiration, violating the linear assumption of the above analysis method. Third, in addition to respiratory movement itself, heart rate is also influenced by a complex interaction between central, neural, hormonal, and mechanical feedback mechanisms. Under these complex or nonlinear interactions, intermittent synchronization between the cardiac and respiratory systems can occur in normal subjects, resulting in temporary variations in the amplitude (e.g., envelope) of the RSA oscillations. The above analysis fails to reflect these interactions. The non-stationary and nonlinear nature of these physiological signals creates numerous problems when reconstructing RSA using only a combination of constant-amplitude sinusoidal oscillations and performing spectrum analysis using Fourier transforms. Not only does the finite window effect of the Fourier transform introduce spurious energy around the dominant frequency and distribute power to harmonic frequencies outside the HF band, but other potential interference, such as ectopic heartbeats, data loss, and noise, can also introduce spurious energy into the spectrum, leading to significant deviations in frequency estimation methods.

[0005] Considering (1) the non-stationary and nonlinear nature of biological systems, which adapt to the ever-present internal and external environmental disturbances, and (2) the irrelevant interference that affects RSA evaluation, Lin et al. proposed a new method, Multimodal Coupling Analysis (MMCA), to extract meaningful physiological oscillations (such as RSA) and their dynamic characteristics. This method was applied to RSA quantification, resolving the respiratory frequency constraint problem in traditional frequency-domain RSA quantification. The core algorithm of MMCA is the Hilbert-Huang Transform (HHT), which uses noise-assisted empirical mode decomposition (EMD), also known as ensemble empirical mode decomposition (EEMD), to decompose the RR interval time series and respiratory signal into a finite number of intrinsic mode functions (IMFs). The IMF with the highest power among the respiratory signal IMFs was selected as the dominant IMF to preserve the original respiratory signal while eliminating the effects of non-stationary interference such as body movement, coughing, or swallowing. The minimum average term of the instantaneous frequency difference between the RR interval signal IMFs and the dominant IMF of the respiratory signal was selected as the frequency-matching IMF to reflect the respiratory system's regulation of heart rate fluctuations. The continuous changes in phase synchronization between the two were combined to study cardiopulmonary interactions, thereby eliminating influences that may be unrelated to RSA. MMCA-derived parameters have been shown to perform better in evaluating RSA in terms of frequency power than high-frequency power based on the fast Fourier transform (FFT), phase coherence based on the cross-wavelet transform, and synchronized squeezed wavelet transform-based RSA reconstruction.

[0006] However, the EEMD algorithm used in MMCA only mitigates modal aliasing of EMD to a certain extent by adding noise and repeating the decomposition. Furthermore, the mixing of the decomposed signal with white noise can alter the characteristics of the extreme points, leading to mismatches in the number and frequency of the decomposed intrinsic mode functions. These issues make MMCA less accurate in extracting IMFs from respiratory signals and RR interval time series, potentially leading to biases in subsequent analysis and affecting the time-frequency resolution of RSA assessment results.

[0007] Therefore, how to conduct universal and accurate RSA quantitative evaluation is an urgent problem to be solved. Summary of the Invention

[0008] In view of this, the present invention provides a cardiopulmonary coupling quantification method and system based on multimodal coupling analysis, which can combine variational modal decomposition with multimodal coupling analysis to achieve cardiopulmonary coupling quantification, and has the advantages of universality and accuracy.

[0009] To achieve the above object, the technical solution of the present invention includes the following steps:

[0010] Step 1: Collect the subject's ECG and respiratory signals.

[0011] Step 2: Extract the heart beat interval RR interval time series from the collected ECG signal.

[0012] Step 3: Use variational mode decomposition to decompose the RR interval time series and respiratory signal respectively to obtain the intrinsic mode function (IMF) of the two time series.

[0013] Step 4: Select the IMF with the largest power among all the IMFs of the respiratory signal as the dominant IMF, which is denoted as IMFA.

[0014] Step 5: Select the IMF that matches the dominant IMF frequency of the respiratory signal from all the IMFs corresponding to the RR interval time series, and record it as IMFB.

[0015] Step 6: Calculate the synchronization index between IMFAF and IMFB by using the instantaneous phase difference between them.

[0016] Step 7: Quantify the instantaneous magnitude of RSA by calculating the power of IMFB during its strong synchronization period.

[0017] Furthermore, in step 1, the subject's ECG signal and respiratory signal are collected, and their corresponding sampling frequencies are set to f s1 and f s2 ,The collected ECG signals and respiratory signals are stored in the computer.

[0018] Furthermore, step 2 is specifically as follows:

[0019] The R-peak position is detected in a single-lead ECG signal, and the RR interval time series is calculated based on the difference between adjacent R-peaks. After extracting the normal RR interval time series, a sliding average filter with a fixed data point window is used to remove outliers caused by erroneous heartbeat detection. When the center point in the window lies outside 20% of the mean, it is removed.

[0020] Furthermore, after step 2 and before step 3, the RR interval series and the respiratory signal are uniformly resampled at a frequency of 8 Hz using cubic spline interpolation.

[0021] Furthermore, step 3 is specifically as follows:

[0022] Variational mode decomposition (VMD) is applied to RR interval sequences and respiratory signals respectively to obtain two sets of intrinsic mode functions.

[0023] The intrinsic mode function in VMD is an amplitude-frequency modulated AM-FM signal, defined as u k (t):

[0024] u k (t) = A k (t)cos(φ k (t))(1)

[0025] Among them A k (t) represents the amplitude envelope, φ k (t) represents the phase and is a non-decreasing function.

[0026] At this time, the instantaneous frequency

[0027] The purpose of VMD is to decompose the real-valued input signal f into K discrete sub-signals u k ,k∈{1,2,...,K}, each mode k has a center frequency of ω k The modal function is solved by constructing a constraint equation that minimizes the sum of the bandwidths of each mode.

[0028] Estimate each modal bandwidth, specifically:

[0029] S301: Calculate the modal component u through Hilbert transform k (t) to obtain a one-sided spectrum:

[0030]

[0031] Where δ(t) is the unit impulse function.

[0032] S302: Tuning to the respective estimated center frequencies ω by k Multiply the exponential of to move the spectrum of each modal component to the "baseband", which is recorded as:

[0033]

[0034] S303: Estimate the bandwidth of each modal component by the Gaussian smoothness of the demodulated signal:

[0035]

[0036] The resulting constrained variational problem is as follows:

[0037]

[0038] S304: Introduce a quadratic penalty term α and a Lagrange multiplier λ to transform the problem into an unconstrained equation.

[0039] Enhanced Lagrangian as follows:

[0040]

[0041] S305: Solving the Augmented Lagrangian Equations Using the Alternating Direction Multiplier Method The saddle point of , Equation (6) can be rewritten as an equivalent minimization problem:

[0042]

[0043] Then the modal component u k (t) and their respective center frequencies ω k They are represented as follows:

[0044]

[0045] in Represent the Fourier transform of each component.

[0046] The RR interval series and respiratory signal are decomposed by variational mode decomposition to obtain their IMFs.

[0047] Furthermore, step 5: selecting the frequency matching IMF of the dominant IMF of the respiratory signal in the RR signal IMFs, denoted as IMFB, specifically includes the following specific steps:

[0048] The Hilbert transform is performed on the IMFs decomposed from the RR signal and the dominant IMF of the respiratory signal to obtain the instantaneous amplitude, phase, and frequency information of each IMF. The Hilbert transform of the time domain sequence x(t) is defined as:

[0049]

[0050] Where P is the Cauchy principal value.

[0051] For the signal x(t), its corresponding analytical signal can be constructed from its Hilbert transform and the original signal:

[0052]

[0053] Where A(t) and are the instantaneous amplitude and instantaneous phase of x(t), respectively.

[0054] The IMF of the RR signal with the smallest interpolated average of the instantaneous frequency with the dominant IMF of the respiratory signal is selected as the frequency matching IMF, reflecting the regulation of the respiratory system on heart rate fluctuations in isolation from the influence of other physiological processes.

[0055] Furthermore, step 6: calculate the synchronization index between the two using the instantaneous phase difference between IMFAF and IMFB, specifically:

[0056] The dominant IMF of the respiratory signal and the frequency matching IMF of the RR sequence are divided into multiple signal segments of length T through a sliding window, and the sliding window moves forward 1 second each time.

[0057] According to the definition of instantaneous phase, MMCA quantifies the temporal dependence between two signals as the synchronization strength within a period of time, which is defined as the synchronization index p. i :

[0058]

[0059] in, and They are the instantaneous phases of the RR interval sequence spectrum matching IMF and the respiratory signal dominant IMF, respectively.

[0060] Normalized synchronization index p i Reflects the strength of cardiopulmonary synchronization, p i =0 corresponds to asynchronous, p i =1 corresponds to complete synchronization, that is, the period of strong synchronization of cardiopulmonary coupling Corresponding to a high synchronization index, p i is greater than the threshold Γ.

[0061] Furthermore, step 7: quantify the instantaneous size of RSA by calculating the power of IMFB during its strong synchronization period, specifically:

[0062] When p i >Γ,(13).

[0063] P vagal (t i )=0 when p i ≤Γ.

[0064] Another embodiment of the present invention further provides a cardiopulmonary coupling quantification system based on multimodal coupling analysis, comprising:

[0065] The acquisition module is used to collect the subject's electrocardiogram signal and respiratory signal.

[0066] The extraction module is used to extract the heart beat interval RR interval time series from the collected electrocardiogram signal.

[0067] The decomposition module is used to decompose the RR interval time series and the respiratory signal respectively by using variational mode decomposition to obtain the intrinsic mode function (IMF) of the two time series.

[0068] The dominant IMF selection module is used to select the IMF with the largest power among all the IMFs of the respiratory signal as the dominant IMF, which is denoted as IMFA.

[0069] The frequency matching module is used to select the IMF that matches the frequency of the dominant IMF of the respiratory signal from all the IMFs corresponding to the RR interval time series, which is recorded as IMFB.

[0070] The synchronization index calculation module is used to calculate the synchronization index between the IMFAF and IMFB according to the instantaneous phase difference between the two.

[0071] The quantization module is used to quantify the instantaneous size of RSA by calculating the power of IMFB during its strong synchronization period.

[0072] Furthermore, the extraction module is specifically as follows: detecting the R peak position in the single-lead ECG signal, and calculating the RR interval time series based on the difference between adjacent R peaks; after extracting the normal RR interval time series, using a sliding average filter with a fixed data point window to remove outliers caused by erroneous heartbeat detection, where when the center point in the window is outside 20% of the average value, it will be removed.

[0073] Decomposition modules, specifically:

[0074] Variational mode decomposition (VMD) is applied to RR interval sequences and respiratory signals respectively to obtain two sets of intrinsic mode functions.

[0075] The intrinsic mode function in VMD is an amplitude-frequency modulated AM-FM signal, defined as u k (t):

[0076] u k (t) = A k (t)cos(φ k (t))(1)

[0077] Among them A k (t) represents the amplitude envelope, φ k (t) represents the phase and is a non-decreasing function.

[0078] At this time, the instantaneous frequency

[0079] The purpose of VMD is to decompose the real-valued input signal f into K discrete sub-signals u k ,k∈{1,2,...,K}, each mode k has a center frequency of ω k The modal function is solved by constructing a constraint equation that minimizes the sum of the bandwidths of each mode.

[0080] Estimate each modal bandwidth, specifically:

[0081] S301: Calculate the modal component u through Hilbert transform k (t) to obtain a one-sided spectrum:

[0082]

[0083] Where δ(t) is the unit impulse function.

[0084] S302: Tuning to the respective estimated center frequencies ω by k Multiply the exponential of to move the spectrum of each modal component to the "baseband", which is recorded as:

[0085]

[0086] S303: Estimate the bandwidth of each modal component by the Gaussian smoothness of the demodulated signal:

[0087]

[0088] The resulting constrained variational problem is as follows:

[0089]

[0090] S304: Introduce a quadratic penalty term α and a Lagrange multiplier λ to transform the problem into an unconstrained equation.

[0091] Enhanced Lagrangian as follows:

[0092]

[0093] S305: Solving the Augmented Lagrangian Equations Using the Alternating Direction Multiplier Method The saddle point of , Equation (6) can be rewritten as an equivalent minimization problem:

[0094]

[0095] Then the modal component u k (t) and their respective center frequencies ω k They are represented as follows:

[0096]

[0097] in Represent the Fourier transform of each component.

[0098] The RR interval series and respiratory signal are decomposed by variational mode decomposition to obtain their IMFs.

[0099] The frequency matching module includes the following specific steps:

[0100] The Hilbert transform is performed on the IMFs decomposed from the RR signal and the dominant IMF of the respiratory signal to obtain the instantaneous amplitude, phase, and frequency information of each IMF. The Hilbert transform of the time domain sequence x(t) is defined as:

[0101]

[0102] Where P is the Cauchy principal value.

[0103] For the signal x(t), its corresponding analytical signal can be constructed from its Hilbert transform and the original signal:

[0104]

[0105] Where A(t) and are the instantaneous amplitude and instantaneous phase of x(t), respectively.

[0106] The IMF of the RR signal with the smallest interpolated average of the instantaneous frequency with the dominant IMF of the respiratory signal is selected as the frequency matching IMF, reflecting the regulation of the respiratory system on heart rate fluctuations in isolation from the influence of other physiological processes.

[0107] The synchronization index calculation module specifically performs the following steps:

[0108] The dominant IMF of the respiratory signal and the frequency matching IMF of the RR sequence are divided into multiple signal segments of length T through a sliding window, and the sliding window moves forward 1 second each time.

[0109] According to the definition of instantaneous phase, MMCA quantifies the temporal dependence between two signals as the synchronization strength within a period of time, which is defined as the synchronization index p. i :

[0110]

[0111] in, and They are the instantaneous phases of the RR interval sequence spectrum matching IMF and the respiratory signal dominant IMF, respectively.

[0112] Normalized synchronization index p i Reflects the strength of cardiopulmonary synchronization, p i =0 corresponds to asynchronous, p i =1 corresponds to complete synchronization, that is, the period of strong synchronization of cardiopulmonary coupling Corresponding to a high synchronization index, p i is greater than the threshold Γ.

[0113] The quantization module quantizes the instantaneous size of RSA by calculating the power of IMFB during its strong synchronization period, specifically:

[0114]

[0115] P vagal (t i )=0 when p i ≤Γ.

[0116] Beneficial effects:

[0117] 1: Based on the advantages of variational mode decomposition in high time-frequency resolution and strong noise robustness in non-stationary physiological signals, and the universality and accuracy of multimodal coupling analysis in quantifying respiratory sinus arrhythmia (RSA), the present invention proposes for the first time a new cardiopulmonary coupling quantification algorithm based on variational mode decomposition and multimodal coupling analysis, which includes the following steps: collecting the subject's electrocardiogram (ECG) and respiratory signals; extracting the heart beat interval (RR) interval time series from the collected ECG signals; using variational mode decomposition to decompose the RR interval time series and the respiratory signal respectively to obtain the intrinsic mode functions (IMFs) of the above two time series; selecting the IMF with the largest power as the dominant IMF among all the IMFs of the respiratory signal; selecting the IMF in the RR signal IMFs that matches the frequency of the dominant IMF of the respiratory signal; calculating the synchronization index between the dominant IMF of the respiratory signal and the frequency-matching IMF of the RR sequence through the instantaneous phase difference between the two; and finally, quantifying the instantaneous size of RSA by calculating the power of the frequency-matching IMF of the RR interval sequence during its strong synchronization period. The present invention aims to meet the actual clinical demand for universal and accurate quantitative indicators of parasympathetic nervous function. Based on the proposed new algorithm, it realizes the accurate quantification of respiratory sinus arrhythmia. It breaks through the technical difficulties such as the traditional RSA quantification being limited to a specific respiratory frequency range, the physiological signal violating the linear assumption of the analysis method, and the significant influence of disturbances unrelated to RSA on the quantification results. It is expected to provide a feasible new approach for autonomic nervous function monitoring.

[0118] 2: The present invention proposes a set of cardiopulmonary coupling quantification techniques based on the variational mode decomposition (VMD) algorithm and multimodal coupling analysis. While being suitable for the decomposition of non-stationary and nonlinear physiological signals, VMD adopts a completely non-recursive mode decomposition method. By continuously iteratively searching for the optimal solution for the mode decomposition, it adaptively updates the optimal center frequency and bandwidth of each intrinsic mode function, greatly improving the problems of EEMD such as modal aliasing and the need to introduce white noise to additionally change the signal characteristics. In addition, VMD is an improvement based on the Wiener filter, which is equivalent to adapting the Wiener filter to multiple adaptive frequency bands, and therefore has better robustness to noise.

[0119] 3: This paper optimizes the cardiopulmonary coupling quantification technology based on multimodal coupling analysis based on the variational modal decomposition algorithm, and innovatively proposes a set of cardiopulmonary coupling algorithms based on variational modal decomposition and multimodal coupling analysis. The proposed algorithm is applied to better preserve and quantify respiratory sinus arrhythmia, so as to more accurately quantify parasympathetic nervous function. BRIEF DESCRIPTION OF THE DRAWINGS

[0120] Figure 1 A flow chart of a cardiopulmonary coupling quantification method based on variational mode decomposition and multimodal coupling analysis is proposed for the invention. DETAILED DESCRIPTION

[0121] The present invention is described in detail below with reference to the accompanying drawings and embodiments.

[0122] In order to overcome the technical difficulties of low accuracy caused by the large limitations and energy leakage of traditional respiratory sinus arrhythmia quantification, the present invention proposes a set of cardiopulmonary coupling quantification algorithms based on variational mode decomposition and multimodal coupling analysis to achieve accurate assessment of respiratory sinus arrhythmia and more accurately quantify parasympathetic nerve function. The technical flow chart of the present invention is as follows: Figure 1 The detailed process is as follows:

[0123] Step 1) Collect the subject's ECG signal and respiratory signal, and set the sampling frequency to f s1 and f s2 , stored in the computer.

[0124] Step 2) Extract the beat-to-beat (RR) time series from the ECG signal. Specifically, the R-peak positions are detected in the single-lead ECG signal. After visual inspection to reselect erroneous or missed detections, the RR interval series is calculated based on the difference between adjacent R-peaks. After extracting the normal RR interval time series, a sliding average filter with a fixed (e.g., 41) data point window is used to remove outliers caused by erroneous beat detections. Points within the window whose center lies outside the 20% range of the mean are removed.

[0125] Step 3) The RR interval series and respiratory signal are uniformly resampled at a frequency of 8 Hz using cubic spline interpolation.

[0126] Step 4) Apply variational mode decomposition (VMD) to the RR interval sequence and the respiratory signal to obtain two sets of intrinsic mode functions. The eigenmode functions in VMD are amplitude-frequency modulated (AM-FM) signals and are defined as:

[0127] u k (t) = A k (t)cos(φ k(t)) (1)

[0128] Among them A k (t) represents the amplitude envelope, φ k (t) represents the phase and is a non-decreasing function, so the instantaneous frequency

[0129] The purpose of VMD is to decompose the real-valued input signal f into K discrete sub-signals u k ,k∈{1,2,...,K}, each mode k has a center frequency of ω k The specific bandwidth of the modal function is estimated by constructing a constraint equation that minimizes the sum of the modal bandwidths. Specifically:

[0130] The first step is to calculate the modal component u through Hilbert transform k (t) to obtain a one-sided spectrum:

[0131]

[0132] Where δ(t) is the unit impulse function;

[0133] In the second step, by tuning to the respective estimated center frequencies ω k Multiply the exponential of to move the spectrum of each modal component to the "baseband", which is recorded as:

[0134]

[0135] The third step is to estimate the bandwidth of each modal component by the Gaussian smoothness of the demodulated signal (i.e. the L2 norm of the gradient):

[0136]

[0137] The resulting constrained variational problem is as follows:

[0138]

[0139] The fourth step is to introduce the quadratic penalty term α and the Lagrange multiplier λ to transform the problem into an unconstrained equation. The quadratic penalty term helps the signal to maintain good accuracy in the presence of additive Gaussian noise, while the Lagrange multiplier ensures strict enforcement of the constraints. as follows:

[0140]

[0141] Step 5: Use the Alternating Direction Method of Multipliers (ADMM) to solve the enhanced Lagrangian equation In order to continuously update the modal u k (t), rewrite Equation (6) as an equivalent minimization problem:

[0142]

[0143] Then the modal component u k (t) and their respective center frequencies ω k They are represented as follows:

[0144]

[0145] in Represent the Fourier transform of each component.

[0146] The RR interval series and respiratory signal are decomposed by variational mode decomposition to obtain their IMFs.

[0147] Step 5) Select the dominant IMF among the respiratory signal IMFs. For the respiratory signal, the IMF with the largest power among all IMFs is selected as the dominant IMF to eliminate the effects of non-stationary interference such as body movement, coughing, or swallowing while preserving the original respiratory signal.

[0148] Step 6) Select the frequency matching IMF of the dominant IMF of the respiratory signal in the RR signal IMFs. Specifically, perform Hilbert transform on the IMFs decomposed from the RR signal and the dominant IMF of the respiratory signal to obtain the instantaneous amplitude, phase, and frequency information of each IMF. The Hilbert transform of the time domain sequence x(t) is defined as:

[0149]

[0150] Where P is the Cauchy principal value. For a signal x(t), its corresponding analytical signal can be constructed from its Hilbert transform and the original signal:

[0151]

[0152] Where A(t) and The IMF of the RR signal with the smallest interpolated average of the instantaneous frequency with the dominant IMF of the respiratory signal is selected as the frequency matching IMF, reflecting the regulation of the respiratory system on heart rate fluctuations in isolation from the influence of other physiological processes.

[0153] Step 7) Divide the dominant IMF of the respiratory signal and the frequency matching IMF of the RR sequence into multiple signal segments of length T using a sliding window, with the sliding window moving forward 1 second each time. Based on the instantaneous phase definition, MMCA quantifies the temporal dependence between the two signals as the synchronization strength over a period of time, which is defined as the synchronization index p.i :

[0154]

[0155] in, and are the instantaneous phases of the RR interval sequence spectrum matching IMF and the respiratory signal dominant IMF, respectively. The normalized synchronization index p i Reflects the strength of cardiopulmonary synchronization, p i =0 corresponds to asynchronous, p i =1 corresponds to complete synchronization. That is, the period of strong synchronization of cardiopulmonary coupling Corresponding to a high synchronization index (p i greater than the threshold Γ).

[0156] Step 8) quantify the instantaneous magnitude of RSA by calculating the power of the RR interval sequence frequency matching IMF during its strong synchronization period:

[0157] When p i >Γ,(13)

[0158] P vagal (t i )=0 when p i ≤Γ.

[0159] This invention has been validated on the exposure therapy for specific phobias dataset in the Physionet database. By plotting the extraction results of VMD-MMCA and EEMD-MMCA for the RSA-related IMF components of the RR interval signal during the anxiety phase of the subjects, we can see that the RSA energy obtained by EEMD-MMCA leaks into multiple IMF components, while VMD-MMCA extracts all the RSA energy into IMF2. In addition, the RSA quantified by VMD-MMCA is higher than that of EEMD-MMCA in both the anxiety phase and the resting phase, with a significant statistical difference, proving that VMD-MMCA solves the adjacent IMF energy leakage phenomenon caused by EEMD and achieves accurate RSA assessment. At the same time, the two parameters calculated by VMD-MMCA, the mean phase synchronization index and the proportion of strong coupling time periods, can distinguish between the anxiety phase and the resting phase, providing assistance for emotion recognition. In response to the clinical needs of accurately assessing RSA and quantifying parasympathetic nervous function, there is a wider application space.

[0160] This invention proposes for the first time a new algorithm for extracting cardiopulmonary coupling features based on variational modal decomposition and multimodal coupling analysis. It is suitable for accurately evaluating the function of the parasympathetic nervous system by quantifying respiratory sinus arrhythmia. It is efficient, reliable, and easy to software.

[0161] In summary, the above are only preferred embodiments of the present invention and are not intended to limit the scope of protection of the present invention. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A cardiopulmonary coupling quantification method based on multimodal coupling analysis, characterized in that: include: Step 1: Collect the subject's ECG and respiratory signals; Step 2: Extract the heart beat interval RR interval time series from the collected ECG signal; Step 3: Decompose the RR interval time series and respiratory signal respectively using variational mode decomposition to obtain the intrinsic mode function (IMF) of the RR interval time series and respiratory signal; Step 4: Select the IMF with the largest power among all the IMFs of the respiratory signal as the dominant IMF, which is recorded as IMFA; Step 5: Select the IMF that matches the dominant IMF frequency of the respiratory signal from all the IMFs corresponding to the RR interval time series, and record it as IMFB; Step 6: Calculate the synchronization index between IMFA and IMFB by the instantaneous phase difference between them; Step 7: Quantify the instantaneous magnitude of respiratory sinus arrhythmia RSA by calculating the power of the IMFB during its strong synchronization period.

2. The cardiopulmonary coupling quantification method based on multimodal coupling analysis according to claim 1, characterized in that: In step 1, the subject's ECG signal and respiratory signal are collected, and their corresponding sampling frequencies are set to f s1 and f s2 ,The collected ECG signals and respiratory signals are stored in the computer.

3. The cardiopulmonary coupling quantification method based on multimodal coupling analysis according to claim 1, characterized in that: The step 2 is specifically as follows: The R peak position is detected in the single-lead ECG signal, and the RR interval time series is calculated based on the difference between adjacent R peaks. After extracting the normal RR interval time series, a sliding average filter with a fixed data point window is used to remove outliers caused by erroneous heartbeat detection. When the center point in the window is outside 20% of the mean value, it will be removed.

4. The cardiopulmonary coupling quantification method based on multimodal coupling analysis according to any one of claims 1 to 3, characterized in that: After step 2 and before step 3, the RR interval sequence and the respiratory signal are uniformly resampled at a frequency of 8 Hz using cubic spline interpolation.

5. The cardiopulmonary coupling quantification method based on multimodal coupling analysis according to claim 1, characterized in that: The step 3 is specifically as follows: Variational mode decomposition (VMD) is applied to RR interval series and respiratory signals respectively to obtain two sets of intrinsic mode functions. The intrinsic mode function in VMD is an amplitude-frequency modulated AM-FM signal, defined as u k (t): u k (t)=A k (t)cos(φ k (t))(1) Among them A k (t) represents the amplitude envelope, φ k (t) represents the phase and is a non-decreasing function; At this time, the instantaneous frequency The purpose of VMD is to decompose the real-valued input signal f into K discrete sub-signals u k ,k∈{1,2,...,K}, each mode k has a center frequency of ω k The modal function is continuously updated and solved by constructing a constraint equation that minimizes the sum of the bandwidths of each mode. Estimate each modal bandwidth, specifically: S301: Calculate the modal component u through Hilbert transform k (t) to obtain a one-sided spectrum: Where δ(t) is the unit impulse function; S302: Tuning to the respective estimated center frequencies ω by k Multiply the exponential of , and move the spectrum of each modal component to the "baseband", which is recorded as: S303: Estimate the bandwidth of each modal component by the Gaussian smoothness of the demodulated signal: The resulting constrained variational problem is as follows: S304: Introduce the quadratic penalty term α and the Lagrange multiplier λ to transform the problem into an unconstrained equation; Enhanced Lagrangian as follows: S305: Solving the Augmented Lagrangian Equations Using the Alternating Direction Multiplier Method The saddle point of , Equation (6) can be rewritten as an equivalent minimization problem: Then the modal component u k (t) and their respective center frequencies ω k They are represented as follows: in represent the Fourier transform of their respective components; The RR interval series and respiratory signal are decomposed by variational mode decomposition to obtain their IMFs.

6. The cardiopulmonary coupling quantification method based on multimodal coupling analysis according to claim 1, characterized in that: The step 5: selecting the IMF that matches the dominant IMF frequency of the respiratory signal from all the IMFs corresponding to the RR interval time series, denoted as IMFB, specifically includes the following specific steps: Perform Hilbert transform on the IMF decomposed from the RR signal and the dominant IMF of the respiratory signal to obtain the instantaneous amplitude, phase, and frequency information of each IMF. The Hilbert transform of the time domain sequence x(t) is defined as: Where P is the Cauchy principal value; For the signal x(t), its corresponding analytical signal can be constructed from its Hilbert transform and the original signal: Where A(t) and are the instantaneous amplitude and instantaneous phase of x(t) respectively; The IMF of the RR signal with the smallest average instantaneous frequency difference with the dominant IMF of the respiratory signal is selected as the frequency matching IMF, reflecting the regulation of the respiratory system on heart rate fluctuations in isolation from the influence of other physiological processes.

7. The cardiopulmonary coupling quantification method based on multimodal coupling analysis according to claim 1, characterized in that: Step 6: Calculate the synchronization index between IMFA and IMFB by using the instantaneous phase difference between the two, specifically: The dominant IMF of the respiratory signal and the frequency matching IMF of the RR sequence are divided into multiple signal segments of length T through a sliding window, and the sliding window moves forward 1 second each time; According to the definition of instantaneous phase, multimodal coupling analysis (MMCA) quantifies the temporal dependence between two signals as the synchronization strength within a period of time, which is defined as the synchronization index p. i : in, and are the instantaneous phases of the RR interval sequence spectrum matching IMF and the respiratory signal dominant IMF, respectively; Normalized synchronization index p i Reflects the strength of cardiopulmonary synchronization, p i =0 corresponds to asynchronous, p i =1 corresponds to complete synchronization, that is, the period of strong synchronization of cardiopulmonary coupling Corresponding to a high synchronization index, p i is greater than the threshold Γ.

8. The cardiopulmonary coupling quantification method based on multimodal coupling analysis according to claim 1, characterized in that: Step 7: quantifying the instantaneous size of RSA by calculating the power of IMFB during its strong synchronization period, specifically:

9. Cardiopulmonary coupling quantification system based on multimodal coupling analysis, characterized by: include: An acquisition module, used to acquire the subject's electrocardiogram signal and respiratory signal; An extraction module is used to extract the heart beat interval RR interval time series from the collected electrocardiogram signal; A decomposition module is used to decompose the RR interval time series and the respiratory signal respectively by using variational mode decomposition to obtain the intrinsic mode function IMF of the RR interval time series and the respiratory signal; The dominant IMF selection module is used to select the IMF with the largest power among all the IMFs of the respiratory signal as the dominant IMF, which is denoted as IMFA; The frequency matching module is used to select the IMF that matches the frequency of the dominant IMF of the respiratory signal from all the IMFs corresponding to the RR interval time series, denoted as IMFB; A synchronization index calculation module is used to calculate the synchronization index between IMFA and IMFB based on the instantaneous phase difference between the two; The quantification module is used to quantify the instantaneous magnitude of respiratory sinus arrhythmia RSA by calculating the power of the IMFB during its strong synchronization period.

10. The cardiopulmonary coupling quantification system based on multimodal coupling analysis according to claim 9, characterized in that: The extraction module specifically detects the R peak position in the single-lead ECG signal and calculates the RR interval time series based on the difference between adjacent R peaks; after extracting the normal RR interval time series, uses a sliding average filter with a fixed data point window to remove outliers caused by erroneous heartbeat detection, wherein when the center point in the window is outside 20% of the average value, it will be removed; The decomposition module is specifically: Variational mode decomposition (VMD) is applied to RR interval series and respiratory signals respectively to obtain two sets of intrinsic mode functions. The intrinsic mode function in VMD is an amplitude-frequency modulated AM-FM signal, defined as u k (t): u k (t)=A k (t)cos(φ k (t))(1) Among them A k (t) represents the amplitude envelope, φ k (t) represents the phase and is a non-decreasing function; At this time, the instantaneous frequency The purpose of VMD is to decompose the real-valued input signal f into K discrete sub-signals u k ,k∈{1,2,...,K}, each mode k has a center frequency of ω k The modal function is continuously updated and solved by constructing a constraint equation that minimizes the sum of the bandwidths of each mode. Estimate each modal bandwidth, specifically: S301: Calculate the modal component u through Hilbert transform k (t) to obtain a one-sided spectrum: Where δ(t) is the unit impulse function; S302: Tuning to the respective estimated center frequencies ω by k Multiply the exponential of , and move the spectrum of each modal component to the "baseband", which is recorded as: S303: Estimate the bandwidth of each modal component by the Gaussian smoothness of the demodulated signal: The resulting constrained variational problem is as follows: S304: Introduce the quadratic penalty term α and the Lagrange multiplier λ to transform the problem into an unconstrained equation; Enhanced Lagrangian as follows: S305: Solving the Augmented Lagrangian Equations Using the Alternating Direction Multiplier Method The saddle point of , Equation (6) can be rewritten as an equivalent minimization problem: Then the modal component u k (t) and their respective center frequencies ω k They are represented as follows: in represent the Fourier transform of their respective components; After decomposing the RR interval series and respiratory signal by variational mode decomposition, the IMF of both are obtained; The frequency matching module includes the following specific steps: Perform Hilbert transform on the IMF decomposed from the RR signal and the dominant IMF of the respiratory signal to obtain the instantaneous amplitude, phase, and frequency information of each IMF. The Hilbert transform of the time domain sequence x(t) is defined as: Where P is the Cauchy principal value; For the signal x(t), its corresponding analytical signal can be constructed from its Hilbert transform and the original signal: Where A(t) and are the instantaneous amplitude and instantaneous phase of x(t) respectively; The IMF of the RR signal with the smallest average instantaneous frequency difference with the dominant IMF of the respiratory signal is selected as the frequency matching IMF, reflecting the regulation of the respiratory system on heart rate fluctuations in isolation from the influence of other physiological processes; The synchronization index calculation module specifically performs the following steps: The dominant IMF of the respiratory signal and the frequency matching IMF of the RR sequence are divided into multiple signal segments of length T through a sliding window, and the sliding window moves forward 1 second each time; According to the definition of instantaneous phase, multimodal coupling analysis (MMCA) quantifies the temporal dependence between two signals as the synchronization strength within a period of time, which is defined as the synchronization index p. i : in, and are the instantaneous phases of the RR interval sequence spectrum matching IMF and the respiratory signal dominant IMF, respectively; Normalized synchronization index p i Reflects the strength of cardiopulmonary synchronization, p i =0 corresponds to asynchronous, p i =1 corresponds to complete synchronization, that is, the period of strong synchronization of cardiopulmonary coupling Corresponding to a high synchronization index, p i greater than the threshold Γ; The quantization module calculates the instantaneous size of the power quantization RSA of the IMFB during its strong synchronization period, specifically:

Images