An adaptive cascaded ECG signal denoising method for electromyographic interference

Abstract

This invention relates to the field of signal denoising technology, and provides an adaptive cascaded method for denoising electromyographic (EMG) signals resistant to EMG interference. The method includes: inputting a noisy EMG signal; optimizing the signal using a gray wolf optimization algorithm to obtain an optimal parameter combination; using this optimal parameter combination to drive variational mode decomposition (VMD) to perform mode decomposition on the noisy EMG signal, obtaining several intrinsic mode function (IMF) components; calculating the correlation coefficients of the IMF components; selecting noise-dominated IMF components and signal-dominated IMF components based on correlation coefficient thresholds; denoising the selected noise-dominated IMF components using wavelet thresholding; and reconstructing the signal using the denoised noise-dominated IMF components and the retained signal-dominated IMF components to obtain a high-fidelity denoised EMG signal. This invention performs directional wavelet thresholding denoising on the noise-dominated components, preserving the key diagnostic features of the EMG waveform.

Description

Technical Field

[0001] This invention relates to the field of signal denoising technology, and in particular to an adaptive cascaded method for denoising electrocardiogram signals in response to electromyographic interference. Background Technology

[0002] The characteristic morphology of electrocardiogram (ECG) signals, such as the P wave and QRS complex, directly determines the accuracy of cardiovascular disease diagnosis. ECG acquisition is susceptible to contamination from electromyography (EMG) interference frequencies of 5–500 Hz. This wideband characteristic overlaps with the effective ECG frequency range of 0.5–40 Hz, causing spectral aliasing and masking key features such as the QRS complex, posing a significant challenge to both manual interpretation and automated diagnosis. Existing denoising methods have limited effectiveness because they rely on fixed parameters and a single strategy, making it difficult to accurately separate noise from ECG characteristics.

[0003] In traditional denoising methods, wavelet transform, while possessing time-frequency localization capabilities, suffers from performance heavily reliant on the selection of wavelet basis, decomposition level, and threshold function, easily leading to high-frequency information loss or over-smoothing. Empirical Mode Decomposition (EMD) offers adaptive decomposition but suffers from mode aliasing and endpoint effects. Variational Mode Decomposition (VMD), by iteratively solving for the optimal solution of the variational model, yields sparse, quasi-orthogonal intrinsic mode functions (IMFs), effectively suppressing mode aliasing. VMD performance is highly dependent on the decomposition level K and the penalty factor. The settings are: manual parameter tuning is highly subjective and lacks consistency, and fixed parameters are difficult to adapt to different signal characteristics, which can easily lead to over-decomposition or under-decomposition, resulting in mode mixing or loss of effective components.

[0004] To address the aforementioned issues, existing technologies combine intelligent optimization algorithms with VMD, such as using Sparrow Search Algorithm (SSA), Improved Whale Optimization Algorithm (IWOA), and Grasshopper Optimization Algorithm (GOA) to adaptively optimize VMD parameters, thereby improving decomposition accuracy and signal-to-noise ratio. Furthermore, hybrid denoising strategies have become a hot topic, including the fusion of Empirical Mode Decomposition (EMD) and Wavelet Transform, Particle Swarm Optimization (PSO) combined with wavelets, Tuna Swarm Optimization (TSO) optimizing fully integrated Empirical Mode Decomposition (CEEMDAN) and stacked Sparse Autoencoder (SSAE) denoising, and the Star Raven Optimization Algorithm (NOA) optimizing the fusion of Discrete Wavelet Transform and Nonlocal Mean Filtering.

[0005] The overall effectiveness of existing technologies remains limited by three core bottlenecks: optimization algorithms are prone to getting trapped in local optima, inter-module collaboration is weak, and robustness to complex noise is insufficient. Therefore, there is an urgent need for an ECG signal denoising method that can overcome these shortcomings in order to effectively solve the problem of accurately separating noise from ECG features caused by spectral aliasing. Summary of the Invention

[0006] To address the aforementioned technical problems, this invention provides an adaptive cascaded method for denoising electromyographic interference, comprising the following steps: The input noisy ECG signal is optimized using the Grey Wolf optimization algorithm to obtain the optimal parameter combination. This optimal parameter combination drives variational mode decomposition (VMD) to perform mode decomposition on the noisy ECG signal, resulting in a first set of several intrinsic mode function (EMF) components. Local templates of the QRS complex wave, obtained through peak detection and short-time averaging, are pre-extracted from the noisy ECG signal. The time-frequency domain characteristics of these templates, including the main peak position, bandwidth, and energy distribution, are transformed into constraints for variational decomposition. The constrained variational mode decomposition is solved to obtain a second set of EMF components. The quality of the first and second sets of EMF components is evaluated, and the set with the better quality is selected. Calculate the correlation coefficients of the intrinsic mode function components, and filter out the noise-dominated intrinsic mode function components and the signal-dominated intrinsic mode function components based on the correlation coefficient threshold. The noise-dominant intrinsic mode function (IMF) components selected were denoised using the wavelet thresholding method. The denoised noise-dominant IMF components were then combined with the retained signal-dominant IMF components to reconstruct the signal, resulting in a high-fidelity denoised ECG signal.

[0007] Preferably, the process of obtaining the first group of several eigenmode function components includes the following steps: The noisy ECG signal is input into the gray wolf optimization algorithm, and the output signal-to-noise ratio is used as the optimization objective. The position of the gray wolf is iteratively updated within the preset search range of the decomposition layer and the penalty factor. After multiple iterations, the optimal parameter combination consisting of the optimal decomposition layer and the optimal penalty factor is obtained. The optimal decomposition level and the optimal penalty factor in the optimal parameter combination are set as the decomposition level parameter and penalty factor parameter of variational mode decomposition, respectively, to drive variational mode decomposition to perform mode decomposition on the input noisy ECG signal. Variational mode decomposition outputs multiple intrinsic mode function components, each of which corresponds to a component in a noisy electrocardiogram signal with different frequency bands.

[0008] Preferably, the process of obtaining the optimal parameter combination consisting of the optimal number of decomposition levels and the optimal penalty factor after multiple iterations includes the following steps: The intensity of high-frequency disturbance is determined based on the frequency domain amplitude change rate of the noisy electrocardiogram signal. The preset search range of the decomposition layer is divided into three sub-intervals based on the intensity of the high-frequency disturbance. A value is then extracted from each sub-interval. The three values ​​are combined with the three preset fixed values ​​of the penalty factor to form three sets of initial candidate parameter pairs. Each pair corresponds to the initial position coordinates of a gray wolf individual. The decomposition level corresponding to the location coordinates of each individual gray wolf and the penalty factor are input into variational mode decomposition to decompose the noisy ECG signal and obtain all intrinsic mode function components under the optimal parameter combination. The difference between the center frequencies of two adjacent intrinsic mode function components is calculated. If the difference is less than the preset proportional threshold of the lower frequency of the two adjacent components, it is determined that the optimal parameter combination produces mode overlap. The decomposition level K is reduced by one step, with the preset step size being 1. At the same time, the penalty factor α is increased by one step, with the preset step size being 10% of the current α value. The updated location coordinates are obtained by increasing the penalty factor. Repeatedly execute the variational mode decomposition, inputting the decomposition level corresponding to the position coordinates of each individual gray wolf and the penalty factor, until the position coordinates of all gray wolves have not moved in two consecutive iterations; output the numerical combination of the optimal decomposition level and the optimal penalty factor.

[0009] Preferably, the process of calculating the difference in center frequencies between two adjacent intrinsic mode function components includes the following steps: All intrinsic mode function (IMF) components output by variational mode decomposition are input into a zero-crossing detector. The detector detects all zero-crossing points in the waveform of each IMF component that change from negative to positive. The time span between adjacent zero-crossing points is recorded, and the reciprocal of the time span is calculated as an instantaneous frequency value. All instantaneous frequency values ​​corresponding to each component are sorted in ascending order, and the value at the middle position is taken as the center frequency of the IMF component. The center frequencies of all components are then arranged in ascending order to obtain a center frequency sequence. Take two adjacent center frequencies sequentially from the center frequency sequence, divide the value of the next center frequency in the center frequency sequence by the value of the previous center frequency to obtain the frequency ratio of the adjacent intrinsic mode function components; then calculate the absolute value of the difference between the frequency ratio and the value 1, and output the absolute value as the frequency separation index of the adjacent intrinsic mode function components. Collect the frequency separation indices of all adjacent intrinsic mode function components and calculate the arithmetic mean of the frequency separation indices; compare the frequency separation index of each pair of adjacent intrinsic mode function components with the arithmetic mean; if the frequency separation index of a certain pair of adjacent intrinsic mode function components is less than the arithmetic mean, the comparison result is recorded as true, indicating that the adjacent intrinsic mode function components meet the mode overlap determination condition; otherwise, it is recorded as false.

[0010] Preferably, the process of obtaining the second set of intrinsic mode function components includes the following steps: The preset constraints, namely that a certain intrinsic mode function component has the highest correlation coefficient with the local template of the QRS composite wave and the instantaneous frequency of the component is limited to the range of 5-20Hz, are placed into each iteration loop of the variational mode decomposition: before the start of each iteration, an empty set of components is initialized and the correlation coefficient flag and instantaneous frequency flag of all intrinsic mode function components are marked as not meeting the standard. After iteratively updating each intrinsic mode function component and its center frequency, the correlation coefficient between each intrinsic mode function component and the local template of the QRS composite wave is calculated sequentially. The intrinsic mode function component with the largest correlation coefficient is identified. If its correlation coefficient is greater than that of all other intrinsic mode function components, its correlation coefficient flag is set to "met the standard". At the same time, the instantaneous frequency sequence of the intrinsic mode function component is calculated. If every frequency value in the instantaneous frequency sequence is greater than or equal to 5Hz and less than or equal to 20Hz, its instantaneous frequency flag is set to "met the standard". If and only if both the correlation coefficient flag and the instantaneous frequency flag have met the criteria, the result of this iteration is accepted, and all current intrinsic mode function components are output as the second set of intrinsic mode function components; otherwise, the result of this update is discarded, each intrinsic mode function component is restored to the value at the end of the previous iteration, and the update step size is reduced before re-executing the iteration.

[0011] Preferably, the process of reverting each intrinsic mode function component to its value at the end of the previous iteration includes the following steps: After each iteration, it is determined whether both the correlation coefficient flag and the instantaneous frequency flag have met the criteria. If both have met the criteria, the result of this iteration is accepted and output. If either flag is not met, the degree of deviation that does not meet the criteria is calculated: for cases where the correlation coefficient is not met, the ratio of the difference between the largest and second-largest correlation coefficients to the largest correlation coefficient is calculated; for cases where the instantaneous frequency is not met, the ratio of the total number of sampling points with instantaneous frequency values ​​below 5Hz or above 20Hz to the total number of sampling points is calculated. The larger of the two calculated ratios is taken as the total deviation coefficient; the recovery amplitude is determined based on the total deviation coefficient: the number of historical iteration steps to be recovered is equal to the preset base recovery steps multiplied by the total deviation coefficient; at the same time, the decay coefficient of the update step size is equal to the preset base decay coefficient multiplied by the total deviation coefficient; Restore each eigenmode function component to its historical value at time step, subtract the recovery step from the current iteration step; multiply the current update step by the decay coefficient to obtain the new step; use the new step to re-execute the current iteration from the restored state.

[0012] Preferably, the process of determining the recovery amplitude based on the total deviation coefficient includes the following steps: After each iteration, it is determined whether the correlation coefficient flag and the instantaneous frequency flag are both met. If both are met, the result of this iteration is accepted and output. If either flag is not met, the percentage of difference when the correlation coefficient is not met is calculated. The difference between the largest and second largest correlation coefficients is divided by the largest correlation coefficient, and the percentage of sampling points exceeding the limit when the instantaneous frequency is not met is calculated. The total number of sampling points with instantaneous frequencies below 5Hz or above 20Hz is divided by the total number of sampling points, and the maximum value of the two is taken as the total deviation coefficient. The recovery steps and step size decay coefficient are determined based on the total deviation coefficient: if the total deviation coefficient is less than or equal to the preset upper limit threshold, the recovery steps are equal to the preset base recovery steps multiplied by the total deviation coefficient, and the step size decay coefficient is equal to the preset base decay coefficient multiplied by the total deviation coefficient; if the total deviation coefficient is greater than the preset upper limit threshold, the recovery steps are directly set to the preset maximum recovery steps, and the step size decay coefficient is set to the square of the base decay coefficient.

[0013] Preferably, the process of constructing a waveform tracking degree classification mechanism based on the correlation coefficient includes the following steps: The amplitudes of all sampling points on the intrinsic mode function component waveform are compared sequentially from left to right with the amplitudes of the corresponding positions in the noisy ECG signal. If the amplitude of the component is greater than the amplitude of the original signal, a high-order flag is recorded; if it is less, a low-order flag is recorded; if they are equal, an equal-order flag is recorded. All flags are concatenated into a flag string according to the sampling point order. Then, the flag string is traversed, and the number of flags contained in the longest segment of consecutive identical flags is counted. The number of flags is used as the waveform following degree of the intrinsic mode function component. The waveform follow-up values ​​of all intrinsic mode function components are sequentially placed into a variable array; The waveform following degree of each component is extracted sequentially, and then compared with the front threshold and the back threshold.

[0014] Preferably, the process of obtaining a high-fidelity denoised ECG signal includes the following steps: Align the denoised noise dominant component set with the retained signal dominant component set according to the same index position in the original decomposition order, and extract one denoised noise dominant component and one signal dominant component; take the position with the same sampling time in the two components as the alignment point, and read the amplitude field of the signal dominant component and the amplitude field of the noise dominant component at each alignment point to form a dual amplitude input pair; For each pair of two amplitude inputs, compare the absolute values ​​of the two amplitudes: store the candidate reconstructed values ​​of all sampling points into a candidate array in chronological order; Starting from the second sampling point to the penultimate sampling point in the candidate array, for each position, extract the three candidate reconstructed values: the previous, current, and next values. Compare the absolute values ​​of the three values, find the value with the middle absolute value, and then extract the sign bit of the value. Combine the sign bit with the middle absolute value to form a new value, and use the new value to replace the original value at the current position in the candidate array. After all positions have been processed, output the candidate array as a high-fidelity denoised ECG signal.

[0015] Preferably, if the absolute value of the dominant signal component is greater than or equal to the absolute value of the dominant noise component, the amplitude of the dominant signal component is selected as the base value; otherwise, the amplitude of the dominant noise component is selected as the base value. The sign bits of the two amplitudes are compared: if the sign bits are the same, the base value is directly output as the candidate reconstructed value of the sampling point; if the sign bits are different, the sign bit of the base value is inverted and output as the candidate reconstructed value.

[0016] This invention constructs a variational mode decomposition parameter closed-loop adaptive optimization strategy with output signal-to-noise ratio as the optimization objective; and utilizes the Grey Wolf algorithm in... The intelligent search within the parameter space of K overcomes the limitations of traditional variational mode decomposition, which relies on manual experience for parameter tuning, thus improving the accuracy and efficiency of mode decomposition. A new cascaded processing paradigm of mode decomposition-component classification-threshold denoising is proposed. The intrinsic mode function components are accurately classified by using the correlation coefficient threshold, and the noise-dominant components are denoised in a targeted manner by combining wavelet threshold functions, achieving efficient denoising while maintaining the integrity of signal features. A joint experimental verification system covering simulated noise and real noise is established. The performance is evaluated through ablation experiments and comparison with multiple algorithms, which confirms the strong robustness and adaptability of the algorithm.

[0017] Other features and advantages of the invention will be set forth in the following description, and will be apparent in part from the description, or may be learned by practicing the invention. The objects and other advantages of the invention may be realized and obtained by means of the structures particularly pointed out in the written description and the accompanying drawings.

[0018] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. Attached Figure Description

[0019] The accompanying drawings are provided to further illustrate the invention and form part of the specification. They are used in conjunction with embodiments of the invention to explain the invention and do not constitute a limitation thereof. In the drawings: Figure 1 This is a flowchart of the adaptive cascaded ECG signal denoising method for electromyography interference in Embodiment 1 of the present invention; Figure 2 The principle of the adaptive cascaded ECG signal denoising method for electromyography interference in Embodiment 1 of the present invention. Figure 1 ; Figure 3 The principle of the adaptive cascaded ECG signal denoising method for electromyography interference in Embodiment 1 of the present invention. Figure 2 ; Figure 4 This is a process diagram of obtaining the first group of several intrinsic mode function components in Embodiment 2 of the present invention; Figure 5 This is a flowchart illustrating the waveform following degree classification mechanism constructed based on correlation coefficient in Embodiment 11 of the present invention. Figure 6 This is a process diagram of obtaining a high-fidelity denoised electrocardiogram signal in Embodiment 12 of the present invention. Detailed Implementation

[0020] The preferred embodiments of the present invention will be described below with reference to the accompanying drawings. It should be understood that the preferred embodiments described herein are for illustration and explanation only and are not intended to limit the present invention.

[0021] The terminology used in the embodiments of this application is for the purpose of describing particular embodiments only and is not intended to limit the embodiments of this application. The singular forms "a," "say," and "this" used in the embodiments of this application are also intended to include the plural forms unless the context clearly indicates otherwise. It should also be understood that the terms used in these embodiments refer to and / or include any or all possible combinations of one or more associated listed items.

[0022] In the following description, when referring to the accompanying drawings, unless otherwise indicated, the same numbers in different drawings represent the same or similar elements. The embodiments described in the following exemplary embodiments do not represent all embodiments consistent with this application. Rather, they are merely examples of apparatuses and methods consistent with some aspects of this application. In the description of this application, it should be understood that the terms first, second, third, etc., are used only to distinguish similar objects and are not necessarily used to describe a specific order or sequence, nor should they be construed as indicating or implying relative importance. Those skilled in the art can understand the specific meaning of the above terms in this application according to the specific circumstances.

[0023] This embodiment combines the Grey Wolf optimization algorithm, intrinsic mode function correlation coefficient classification, and wavelet thresholding method. The Grey Wolf optimization algorithm uses the output signal-to-noise ratio as the optimization objective and adaptively optimizes within a preset search space of decomposition level and penalty factor. This allows variational mode decomposition to dynamically adjust decomposition parameters according to the real-time intensity and spectral distribution of EMG interference, suppressing mode aliasing at the source and achieving sparse frequency band decomposition. The waveform following degree classification mechanism based on correlation coefficients accurately divides intrinsic mode function components into noise-carrying components and signal-carrying components through sign comparison and dynamic threshold segmentation, avoiding the loss of ECG features caused by fixed thresholds or empirical division in traditional methods. The wavelet thresholding method performs directional processing on the noise-carrying components, utilizing its time-frequency localization capability to preserve abrupt change points, such as steep edges of the R wave, while eliminating residual noise. Finally, the denoised components and signal-carrying components are reconstructed, forming a cascade chain of adaptive decomposition-precise classification-targeted denoising, which significantly improves the denoising fidelity and diagnostic feature identifiability of ECG signals under EMG interference conditions.

[0024] Existing hybrid denoising methods often suffer from local optima, weak inter-module synergy, and insufficient robustness to complex noise. This invention overcomes the subjectivity of manual parameter tuning and the poor adaptability of fixed parameters to dynamic noise by using the Grey Wolf optimization algorithm and closed-loop adaptive parameter tuning based on variational mode decomposition. It replaces the conventional correlation coefficient thresholding method with a sign comparison classification strategy based on waveform following degree, solving the problem of accurate separation of feature components under spectral aliasing. Furthermore, it employs wavelet thresholding for secondary filtering of noise-dominant components instead of direct discarding, preserving high-frequency diagnostic details. The overall method constructs a full-link adaptive cascade framework from parameter optimization, mode decomposition, component classification to targeted denoising. Under real EMG interference conditions, it improves the output signal-to-noise ratio by 3.8–4.52 dB and reduces the root mean square error to approximately 0.03, representing significant improvements in noise suppression, signal fidelity, and adaptability compared to existing technologies.

[0025] Example 1: As Figure 1 As shown, this embodiment of the invention provides an adaptive cascaded ECG signal denoising method for electromyographic interference, comprising the following steps: S100: Input a noisy ECG signal, and obtain the optimal parameter combination by the Grey Wolf Optimization Algorithm. The optimal parameter combination drives variational mode decomposition to perform mode decomposition on the noisy ECG signal, and obtain the first group of several intrinsic mode function components. A local template of the QRS complex, obtained through peak detection and short-time averaging, is pre-extracted from noisy ECG signals. The time-frequency domain characteristics of the local template, including the main peak position, bandwidth, and energy distribution, are transformed into variational decomposition constraints. The constraints include: one of the intrinsic mode function components obtained by forced decomposition has the highest correlation coefficient with the local template of the QRS complex, and the instantaneous frequency of the intrinsic mode function component is constrained to the range of 5-20Hz. The constrained variational mode decomposition is solved to obtain a second set of intrinsic mode function components. The constraints are incorporated into each iteration of the variational mode decomposition: After each iteration, the correlation coefficient between each intrinsic mode function component and the local template of the QRS composite wave is calculated, the component with the highest correlation coefficient is marked, and it is checked whether the instantaneous frequency of the component falls entirely within the 5-20Hz range; only when the component simultaneously satisfies both the highest correlation coefficient and the instantaneous frequency being within the 5-20Hz range, the iteration result of this round is accepted and the second set of intrinsic mode function components is output; otherwise, the update of this round is discarded, each component is restored to the value at the end of the previous iteration, the update step size is reduced, and the iteration is re-executed until the constraints are satisfied. The first set of intrinsic mode function components and the second set of intrinsic mode function components are evaluated for quality. The evaluation indicators include signal-to-noise ratio improvement, QRS peak fidelity, and electromyographic residual energy. The best set of intrinsic mode function components is selected as the input for subsequent steps, either modally or overall. S200: Calculate the correlation coefficient of the intrinsic mode function components, and filter out the noise-dominated intrinsic mode function components and the signal-dominated intrinsic mode function components based on the correlation coefficient threshold. S300: The selected noise-dominant intrinsic mode function (IMF) components are denoised using the wavelet thresholding method. The denoised noise-dominant IMF components are then combined with the retained signal-dominant IMF components to reconstruct the signal, resulting in a high-fidelity denoised ECG signal.

[0026] The wavelet thresholding method for denoising includes: performing multi-scale wavelet decomposition on the selected noise-dominant intrinsic mode function (EMF) components, selecting a predetermined wavelet basis function and decomposition level, and decomposing each noise-dominant component into low-frequency approximation coefficients and multi-level high-frequency detail coefficients; secondly, applying a threshold function to each level of high-frequency detail coefficients, calculating an adaptive threshold based on the noise energy distribution, setting coefficients with amplitudes below the threshold to zero or performing soft thresholding, and retaining or attenuating coefficients above the threshold, thereby suppressing noise components; finally, performing inverse wavelet transform on the thresholded high-frequency detail coefficients and the unprocessed low-frequency approximation coefficients to reconstruct the denoised EMF components.

[0027] In this embodiment, addressing the characteristics of high-frequency, non-stationary, and aliased electromyography (EMG) interference with the ECG spectrum, the Gray Wolf Optimization Algorithm, Variational Mode Decomposition (VMD), Intrinsic Mode Functions (IMFs), and Wavelet Throttling Method are combined to construct a cascaded processing chain of adaptive decomposition, precise classification, and targeted denoising. Their synergistic effects are as follows: The Gray Wolf Optimization Algorithm drives adaptive VMD decomposition. The intensity and spectral distribution of EMG interference vary from person to person and from time to time. Traditional VMD requires manual setting of the number of modes K and the penalty factor α, making it difficult to match dynamic noise. The Gray Wolf Optimization Algorithm uses the output signal-to-noise ratio as the optimization objective and automatically searches for the optimal (K, α), enabling VMD to decompose noisy ECG signals into a series of sparse intrinsic mode functions. This solves the problems of parameter dependence on experience and easy getting trapped in local optima, achieving adaptive mode decomposition for EMG interference levels and laying the foundation for noise separation. The correlation coefficient classification of intrinsic mode functions (IMFs) shows that EMG interference energy is mainly concentrated in high-frequency IMFs, while ECG features such as P waves, QRS waves, and T waves are concentrated in low-frequency IMFs. By calculating the correlation coefficients between each IMF and the original signal, and setting thresholds, noise-dominant IMFs can be automatically distinguished, with low correlation coefficients indicating signal-dominant IMFs and high correlation coefficients indicating signal-dominant IMFs. This accurately locates the modal components where EMG interference occurs, avoiding the loss of ECG details caused by the one-size-fits-all approach of traditional filtering. Wavelet thresholding is used for targeted denoising of noise-dominant IMFs. Even within noise-dominant IMFs, some high-frequency ECG components may still remain, such as steep R-wave edges. Wavelet thresholding utilizes its time-frequency localization capability to perform multi-scale decomposition on these IMFs, setting small-amplitude noise coefficients to zero or shrinking them while retaining large-amplitude abrupt change coefficients, corresponding to ECG spikes. Compared to directly discarding noise-dominant IMFs, wavelet denoising can suppress EMG residues while preserving ECG diagnostic features to the greatest extent possible. Cascaded Reconstruction: The signal-dominant intrinsic mode function (EMF) and the denoised EMF are added together to reconstruct the final ECG signal. This cascaded paradigm first uses variational mode decomposition to coarsely separate the main noise frequency bands, and then uses wavelet filtering to eliminate residual noise within the bands, forming a complementary and enhanced denoising effect. Experiments show that under EMG interference, the output signal-to-noise ratio is significantly improved, increasing by 3.8~4.52dB under real noise, and the root mean square error is reduced to about 0.03, confirming its strong robustness and clinical applicability.

[0028] The working principle and beneficial effects of the above technical solution are as follows: (Refer to Appendix for details.) Figure 2 and attached Figure 3This embodiment adaptively determines the optimal parameter combination for variational mode decomposition using the Grey Wolf optimization algorithm, achieving high-precision mode decomposition of noisy ECG signals. By filtering noise-dominant and signal-dominant components based on correlation coefficient thresholds, electromyographic interference and ECG feature components can be effectively separated. The noise components after wavelet threshold denoising and the retained signal components are reconstructed, suppressing electromyographic noise while maximizing the preservation of ECG waveform details, ultimately obtaining an ECG signal with high signal-to-noise ratio and clinical diagnostic value.

[0029] Compared with traditional intelligent optimization algorithms, the Grey Wolf Optimization (GWO) algorithm in this embodiment has significant advantages such as fewer parameters, faster convergence speed, and stronger global search capability, solving the bottleneck problem of easily getting trapped in local optima and slow convergence in variational mode decomposition parameter optimization. Furthermore, wavelet thresholding is used for deep denoising of noise-dominated intrinsic mode function (IMF) components, mainly based on its accurate time-frequency localization capability and efficient sparse signal separation characteristics. It can achieve noise denoising for specific frequency bands of IMF components, effectively avoiding the problem of loss of effective components or over-smoothing that traditional single denoising methods easily cause when processing spectral aliasing signals. Thus, while suppressing high-frequency noise, it preserves the key diagnostic features of the ECG waveform to the greatest extent. Specifically, the hybrid strategy of deeply combining the Grey Wolf Optimization, variational mode decomposition, and wavelet thresholding has the core advantage of constructing a fully adaptive, functionally complementary cascaded denoising paradigm. Compared to traditional hybrid algorithms, this strategy optimizes variational mode decomposition parameters using the Grey Wolf optimization algorithm, ensuring decomposition accuracy from the source and overcoming the performance bottleneck of optimization algorithms. Furthermore, through the synergistic mechanism of variational mode decomposition coarse screening and wavelet thresholding fine filtering, it achieves targeted noise removal and maximizes waveform feature preservation, effectively solving the problems of weak module synergy and insufficient robustness. This embodiment uses the Grey Wolf optimization algorithm with the output signal-to-noise ratio as the optimization objective, automatically searching for the optimal number of modes and penalty factor in variational mode decomposition to achieve adaptive mode decomposition of noisy ECG signals. This overcomes the subjectivity of manual parameter selection, making the frequency bands of the decomposed intrinsic mode function components sparse, providing a precise foundation for noise separation. The time-frequency domain characteristics of the pre-extracted QRS complex local template are transformed into variational decomposition constraints, forcing a certain intrinsic mode function component to simultaneously satisfy the highest correlation coefficient and instantaneous frequency within the range of five to twenty Hz. This ensures that the component highly matches the core morphology and physiological frequency band of the QRS complex, effectively suppressing high-frequency aliasing of electromyographic interference and low-frequency intrusion of baseline drift. The two sets of decomposition results were evaluated for signal-to-noise ratio improvement, QRS peak fidelity, and residual EMG energy, and the optimal result was selected to further ensure the waveform integrity of the subsequent processed signal. The correlation coefficients between the selected intrinsic mode function components and the original signal were calculated, and the components were accurately classified into noise-dominant and signal-dominant categories based on thresholds, locating the mode of EMG interference and avoiding the loss of ECG details caused by the one-size-fits-all approach of traditional filtering. Wavelet threshold denoising was applied to the noise-dominant components, utilizing the time-frequency localization capability of wavelet transform to suppress residual EMG noise while preserving abrupt changes such as steep ECG edges in the components. Finally, the denoised noise-dominant components and the unprocessed signal-dominant components were reconstructed to output a high-fidelity ECG signal.The entire cascaded processing chain forms a synergistic enhancement mechanism of adaptive decomposition, precise classification, and targeted noise reduction, which significantly improves the output signal-to-noise ratio and reduces the root mean square error under dynamic electromyography interference environment, while preserving diagnostic features such as QRS wave, P wave, and T wave to the maximum extent, achieving strong robustness and clinical applicability.

[0030] This embodiment proposes a denoising algorithm for electrocardiogram (ECG) signals based on GWO-VMD-Wavelet. The Grey Wolf optimization algorithm can more efficiently determine the optimal parameter combination for variational mode decomposition; then, the correlation coefficient method is used to distinguish noise and signal components in the intrinsic mode function (IMF) components. Adaptive wavelet thresholding is applied to the noise-dominated IMF components, and finally, the denoised ECG signal is reconstructed.

[0031] This embodiment constructs a variational mode decomposition parameter closed-loop adaptive optimization strategy with output signal-to-noise ratio as the optimization objective; and uses the Grey Wolf algorithm in... The intelligent search within the parameter space of K overcomes the limitations of traditional variational mode decomposition, which relies on manual experience for parameter tuning, thus improving the accuracy and efficiency of mode decomposition. A new cascaded processing paradigm of mode decomposition-component classification-threshold denoising is proposed. The intrinsic mode function components are accurately classified by using the correlation coefficient threshold, and the noise-dominant components are denoised in a targeted manner by combining wavelet threshold functions, achieving efficient denoising while maintaining the integrity of signal features. A joint experimental verification system covering simulated noise and real noise is established. The performance is evaluated through ablation experiments and comparison with multiple algorithms, which confirms the strong robustness and adaptability of the algorithm.

[0032] This embodiment is based on the Grey Wolf Optimization Algorithm (GWO) to optimize variational mode decomposition and combines it with wavelet thresholding to create a hybrid denoising model, GWO-VMD-Wavelet. First, the Grey Wolf Optimization Algorithm is used to adaptively search for the number of modes in the variational mode decomposition. and penalty factor This approach overcomes the subjectivity of traditional parameter selection, improving both the accuracy and efficiency of mode decomposition. Secondly, it introduces an intrinsic mode function classification method based on correlation coefficient thresholds to accurately distinguish between signal and noise-dominant components, effectively solving the feature extraction problem caused by spectral aliasing. Finally, a cascaded processing strategy of decomposition-classification-denoising is constructed to perform targeted wavelet threshold denoising on the noise-dominant components, preserving the key diagnostic features of the ECG waveform.

[0033] This embodiment conducts experiments on simulated datasets and real noise databases. Under simulated noise environments of varying intensities, the average output signal-to-noise ratio (SNR) reaches 18.15 dB, with a root mean square error (RMSE) of 0.059. Under real noise environments, the output SNR is improved by 3.8 dB and 4.52 dB, respectively, while the RMS is reduced by 0.034 and 0.029, respectively. Comprehensive verification demonstrates that the algorithm proposed in this embodiment possesses strong robustness and adaptability.

[0034] Example 2: As Figure 4 As shown, based on Example 1, the process for obtaining the first group of several intrinsic mode function components provided by this embodiment of the invention specifically includes the following steps: S101: Input the noisy ECG signal into the gray wolf optimization algorithm, take the output signal-to-noise ratio as the optimization target, iteratively update the position of the gray wolf individuals within the preset search range of the decomposition layer and the penalty factor, and after multiple iterations, converge to obtain the optimal parameter combination composed of the optimal decomposition layer and the optimal penalty factor. S102: Set the optimal decomposition level and optimal penalty factor in the optimal parameter combination as the decomposition level parameter and penalty factor parameter of variational mode decomposition, respectively, and drive variational mode decomposition to perform mode decomposition on the input noisy ECG signal. S103: Variational mode decomposition outputs multiple intrinsic mode function components, each of which corresponds to a component in a noisy ECG signal with different frequency bands.

[0035] The working principle and beneficial effects of the above technical solution are as follows: This embodiment achieves effective decomposition of noisy ECG signals through the synergistic cooperation of the Grey Wolf optimization algorithm and variational mode decomposition. First, the Grey Wolf optimization algorithm is used to adaptively optimize the key parameters of variational mode decomposition; by using the output signal-to-noise ratio as the optimization target, the number of decomposition layers and the penalty factor are dynamically adjusted within a preset search space, significantly improving the accuracy and efficiency of parameter selection and avoiding the limitations of traditional manual trial and error methods; the final optimal parameter combination enables variational mode decomposition to achieve the best noise reduction and decomposition performance. Second, variational mode decomposition driven by optimal parameters can accurately separate the frequency bands of ECG signals; during the decomposition process, the optimal penalty factor effectively suppresses mode aliasing, while the optimal number of decomposition layers ensures the complete extraction of signal components, making each intrinsic mode function component strictly correspond to different frequency components. The final output of multiple intrinsic mode function components fully preserves the clinical feature information of the electrocardiogram signal; each component represents a specific frequency band range, where the high-frequency component can reflect noise such as electromyographic interference, while the low-frequency component contains the main features of the electrocardiogram waveform; the frequency band separation characteristic lays an important foundation for noise elimination and feature extraction, and also enhances the pertinence and interpretability of signal processing.

[0036] Example 3: Based on Example 2, the process of obtaining the optimal parameter combination consisting of the optimal number of decomposition layers and the optimal penalty factor after multiple iterations, provided by this embodiment of the invention, specifically includes the following steps: S1011: Determine the high-frequency disturbance intensity based on the frequency domain amplitude change rate of the noisy ECG signal, divide the preset search range of the decomposition layer into three sub-intervals based on the high-frequency disturbance intensity, and extract a value from each sub-interval; combine the three values ​​with the three preset fixed values ​​of the penalty factor to form three sets of initial candidate parameter pairs, each pair corresponding to the initial position coordinates of a gray wolf individual. S1012: Input the decomposition level corresponding to the position coordinates of each individual gray wolf and the penalty factor into variational mode decomposition to decompose the noisy ECG signal and obtain all intrinsic mode function components under the optimal parameter combination; calculate the difference in center frequencies between two adjacent intrinsic mode function components; if the difference is less than the preset proportional threshold of the lower frequency of the two adjacent components, it is determined that the optimal parameter combination produces mode overlap, and the decomposition level K is reduced by one step, with the preset step size being 1; at the same time, the penalty factor α is increased by one step, with the preset step size being 10% of the current α value; and the updated position coordinates are obtained by increasing the penalty factor. S1013: Repeatedly execute the variational mode decomposition, inputting the decomposition level corresponding to the position coordinates of each gray wolf individual and the penalty factor, until the position coordinates of all gray wolves individuals have not moved in two consecutive iterations; output the numerical combination of the optimal decomposition level and the optimal penalty factor.

[0037] The working principle and beneficial effects of the above technical solution are as follows: This embodiment dynamically divides the high-frequency disturbance intensity range by the time-domain amplitude change rate, thereby achieving adaptive adjustment of the parameter search range. The division of the three sub-intervals ensures the balance between global search and local optimization. Combined with the fixed gradient setting of the penalty factor, a three-dimensional covering network of the parameter space is formed, effectively solving the computational redundancy problem caused by parameter coupling in the traditional grid search method. A dynamic threshold determination mechanism based on the center frequency difference is adopted. Through the negative feedback adjustment of the individual coordinates of the gray wolf, the decomposition layer and the penalty factor are adjusted synchronously. A saddle point constraint is established in the signal frequency domain feature space, so that the intrinsic mode components achieve equal-interval approximation in the frequency band distribution, thereby significantly reducing the probability of high-frequency mode overlap. The dual termination conditions of coordinate stability and mode separation are designed to construct a composite convergence criterion, which avoids premature convergence and prevents over-decomposition. In the MIT-BIH noise database test, it shows stable parameter output characteristics and maintains the stability of the scheme for input signal-to-noise ratio changes of 0.5-5mV.

[0038] Example 4: Based on Example 3, the process for calculating the difference in center frequencies of two adjacent intrinsic mode function components provided in this embodiment of the invention specifically includes the following steps: S10121: Input all intrinsic mode function components output by variational mode decomposition into a zero-crossing detector, detect all zero-crossing points in the waveform of each intrinsic mode function component that change from negative to positive, record the time span between adjacent zero-crossing points, and calculate the reciprocal of the time span as an instantaneous frequency value; sort all the instantaneous frequency values ​​corresponding to each component in ascending order, take the value of the middle position, output the value as the center frequency of the intrinsic mode function component, and arrange the center frequencies of all components in ascending order to obtain the center frequency sequence; S10122: Take two adjacent center frequencies sequentially from the center frequency sequence, divide the value of the next center frequency in the center frequency sequence by the value of the previous center frequency to obtain the frequency ratio of the adjacent intrinsic mode function components; then calculate the absolute value of the difference between the frequency ratio and the value 1, and output the absolute value as the frequency separation index of the adjacent intrinsic mode function components. S10123: Collect the frequency separation indices of all adjacent intrinsic mode function components and calculate the arithmetic mean of the frequency separation indices; compare the frequency separation index of each pair of adjacent intrinsic mode function components with the arithmetic mean; if the frequency separation index of a certain pair of adjacent intrinsic mode function components is less than the arithmetic mean, then the comparison result is recorded as true, indicating that the adjacent intrinsic mode function components meet the mode overlap determination condition; otherwise, it is recorded as false.

[0039] The working principle and beneficial effects of the above technical solution are as follows: This embodiment obtains the center frequency of each IMF component through instantaneous frequency calculation, and accurately captures the signal frequency characteristics using a zero-crossing detection method, ensuring the accuracy of frequency parameter extraction. By using adjacent center frequency ratio calculation combined with normalization processing, the dimensional influence caused by absolute frequency value differences can be eliminated, forming a dimensionless modal separation evaluation index. By constructing a mean comparison mechanism for the frequency separation index, an adaptive modal overlap judgment threshold is established, avoiding the problem of insufficient environmental adaptability caused by preset fixed thresholds. The overall method forms a complete modal aliasing detection chain: from single-component frequency feature extraction, quantification of adjacent component frequency relationships, and dynamic discrimination based on statistical distribution, it has good anti-noise interference capability and algorithm robustness. This embodiment can effectively identify modal confusion phenomena in VMD decomposition results, providing a quantitative basis for subsequent selection of the optimal decomposition layer or adjustment of penalty parameters, thereby improving the accuracy of time-frequency analysis.

[0040] Example 5: Based on Example 4, the process of detecting all zero-crossing points in the waveform of each intrinsic mode function component that change from negative to positive values, provided in this embodiment of the invention, specifically includes the following steps: S101211: Represent the intrinsic mode function components as a queue of sampling points arranged in chronological order. Each sampling point contains a time field and an amplitude field. Starting from the head of the queue, take out the current sampling point and the next sampling point in sequence, and determine whether the amplitude of the current sampling point is less than zero and whether the amplitude of the next sampling point is greater than zero. If the conditions are met, pack the time and amplitude of the current sampling point and the time and amplitude of the next sampling point into a quadruple and append it to the end of the quadruple queue. Repeat until the entire queue has been traversed, and output the quadruple queue. S101212: Take each quad tuple sequentially from the head of the quad tuple queue and perform the following operations: Subtract the amplitude of the current sampling point from the amplitude of the next sampling point in the quad tuple, and store the difference in a temporary variable as the amplitude difference; Multiply the absolute value of the amplitude of the current sampling point by the time difference between the current sampling point and the next sampling point, and store the product in a temporary variable as the first product; Divide the first product by the amplitude difference; Store the quotient in a temporary variable as the offset; Add the time of the current sampling point to the offset, and store the sum at the end of the zero-crossing time array; After processing all quad tuples, output the zero-crossing time array. S101213: Starting from the second element of the time span array at the zero point, take out the current element and the previous element in sequence, and subtract the value of the previous element from the value of the current element; append the difference to the end of the time span array; after traversing, output the time span array.

[0041] The working principle and beneficial effects of the above technical solution are as follows: This embodiment accurately captures the zero-crossing feature where the signal amplitude changes from negative to positive, and uses a linear interpolation algorithm to compensate for the quantization error caused by discrete sampling, thereby improving the subsampling accuracy of zero-crossing detection. A four-tuple data structure and a dual-buffering processing mechanism are adopted to maintain the integrity of the original sampled data while establishing a calculation buffer for the zero-crossing neighborhood features, ensuring the numerical stability of the extracted instantaneous parameters. An interpolation model is constructed based on the dual linear operations of amplitude difference and time difference, effectively overcoming the influence of sampling rate limitations on the detection of zero-crossing points of high-frequency signals, making the calculation results closer to the theoretical zero-crossing moment. An automatic generation of adjacent zero-crossing interval sequences provides a basic timescale parameter for instantaneous frequency calculation, accurately reflecting the changing patterns of the signal's local periodic characteristics.

[0042] Example 6: Based on Example 5, the process of outputting a quadruple queue provided in this embodiment of the invention specifically includes the following steps: S1012111: Compare the amplitude of each sampling point in the sampling point queue of the intrinsic mode function component with the value of zero. If the amplitude is less than zero, set the sign bit of the sampling point to 0. If the amplitude is greater than or equal to zero, set it to 1. Store the sign bits into the sign bit array in the original order of the sampling points. S1012112: Starting from the first element of the sign bit array, take the current element and the next element in sequence, and determine whether the current element is equal to 0 and the next element is equal to 1. If the conditions are met, record the index position of the current element in the sign bit array and store the index position in the boundary index array. S1012113: Take each index position in the boundary index array in sequence, extract the sampling point corresponding to the index position from the original sampling point queue as the pre-sampling point, and extract the sampling point corresponding to the index position plus one as the post-sampling point; pack the time and amplitude of the pre-sampling point and the time and amplitude of the post-sampling point into a quadruple, and store all quadruples into the quadruple queue in the order of the boundary index array.

[0043] The working principle and beneficial effects of the above technical solution are as follows: This embodiment uses sign bit conversion technology to discretize continuous amplitude values, significantly reducing the complexity of subsequent calculations through binary state switching detection, while preserving the key phase features of the original waveform. A fast positioning mechanism based on index boundary detection is established, which can accurately capture the key transition interval where the signal amplitude crosses zero, avoiding redundant calculations caused by full data domain scanning. By constructing a feature packing strategy for boundary sampling point pairs, a complete transition state description unit containing the time-amplitude relationship between preceding and following moments is formed, providing structured input for sub-sampling precision interpolation calculations. Automated conversion from the original signal to feature tuples is achieved, ensuring that the transition interval data required for instantaneous frequency calculation has a strict temporal correspondence. This embodiment solves the problems of missed detection and false detection in traditional zero-crossing detection methods in discrete sampling systems, improving the robustness of detection for weak and non-stationary signals.

[0044] Example 7: Based on Example 6, the process of packaging the time and amplitude of the previous sampling point and the time and amplitude of the subsequent sampling point into a quadruple provided in this embodiment of the invention specifically includes the following steps: S10123331: Read the index value of the current position from the boundary index array, use the index value as an offset to retrieve the time field and amplitude field of the previous sampling point from the original sampling point array, and then use the index value plus one as an offset to retrieve the time field and amplitude field of the next sampling point; write the time field, amplitude field of the previous sampling point, time field of the next sampling point, and amplitude field of the next sampling point into four consecutively allocated storage units in chronological order, with each unit storing one field. S10123332: The starting address of four storage units is used as a pointer, and the pointer is bound to a counting flag with a length of four to form a data chain pointing to a contiguous storage area; the data chain is pushed to the tail free position of a pre-built linear storage pool, and the starting address of the tail free position is recorded in the index pointer array. S10123333: Retrieve all recorded addresses in the index pointer array in the order of writing. For each address retrieved, locate the corresponding four consecutive storage units. Concatenate the field values ​​in the four consecutive storage units in order into an indivisible composite data block and append the composite data block to the end of the output structure. The length counter of the output structure is incremented synchronously.

[0045] The working principle and beneficial effects of the above technical solution are as follows: This embodiment solves the problems of missed detection and false detection in the traditional zero-crossing detection method in discrete sampling systems, and improves the robustness of detection of weak signals and non-stationary signals.

[0046] Example 8: Based on Example 1, the process for obtaining the second set of intrinsic mode function components provided in this embodiment of the invention specifically includes the following steps: S104: The preset constraints, namely that a certain intrinsic mode function component has the highest correlation coefficient with the local template of the QRS composite wave and the instantaneous frequency of the component is limited to the range of 5-20Hz, are placed into each iteration loop of the variational mode decomposition: Before the start of each iteration, an empty set of components is initialized, and the correlation coefficient flag and instantaneous frequency flag of all intrinsic mode function components are marked as not meeting the standard. The instantaneous frequency of the QRS complex local template is constrained to the range of 5-20Hz. This is because the main energy of the QRS complex in ECG signals is concentrated in this frequency band, with a peak of about 4-12Hz and an average of about 15Hz. This is a common understanding in the field of ECG signal processing. For example, the Pan-Tompkins algorithm uses 5-15Hz, and clinical equipment is set to 5-20Hz. Its function is to force the intrinsic mode function components extracted by variational mode decomposition in the iteration to strictly satisfy the frequency band constraint, accurately separate the core mode components of the corresponding QRS complex, and avoid the mixing of high-frequency electromyographic interference and low-frequency baseline drift. It significantly improves the targeting and accuracy of mode decomposition, ensures that the output components meet the standards in both waveform matching and physiological frequency band, provides a high-quality signal basis for correlation coefficient classification and targeted denoising, effectively suppresses the erroneous decomposition caused by spectral aliasing, and ultimately enhances the robustness of the denoising algorithm and the fidelity of the ECG signal. S105: After iteratively updating each intrinsic mode function component and its center frequency, calculate the correlation coefficient between each intrinsic mode function component and the local template of the QRS composite wave, and find the intrinsic mode function component with the largest correlation coefficient. If its correlation coefficient is greater than that of all other intrinsic mode function components, set its correlation coefficient flag to "met the standard". At the same time, calculate the instantaneous frequency sequence of the intrinsic mode function component. If each frequency value in the instantaneous frequency sequence is greater than or equal to 5Hz and less than or equal to 20Hz, set its instantaneous frequency flag to "met the standard". The correlation coefficient, specifically the Pearson correlation coefficient, is used to quantify the waveform similarity between each intrinsic mode function component and the pre-extracted local template of the QRS complex wave. By calculating the correlation coefficient between each intrinsic mode function component and the QRS template, the component with the largest correlation coefficient is identified – this component is considered most likely to carry the core morphological information of the QRS group. Setting its correlation coefficient flag to "met the standard" is a prerequisite for forcing the variational mode decomposition iteration to retain this component and satisfy the instantaneous frequency constraint, ensuring that the second set of intrinsic mode function components decomposed achieves optimality under both waveform matching degree and physiological frequency band standards. S106: If and only if both the correlation coefficient flag and the instantaneous frequency flag have met the criteria, accept the results of this iteration and output all current intrinsic mode function components as the second set of intrinsic mode function components; otherwise, discard the results of this update, restore each intrinsic mode function component to the value at the end of the previous iteration, reduce the update step size, and re-execute the iteration.

[0047] The working principle and beneficial effects of the above technical solution are as follows: This embodiment combines correlation coefficient constraints with instantaneous frequency range limitations, and adopts an iterative mechanism with dynamically adjusted step size to achieve the following technical effects: Improved decomposition accuracy: Constraints force the VMD algorithm to prioritize the selection of IMF components with the highest correlation to the QRS complex morphology in each iteration, while ensuring that their instantaneous frequency is within the typical frequency band of ECG characteristic waves (5-20Hz), eliminating interference from high-frequency noise and low-frequency baseline drift. Guaranteed component validity: A dual-flag verification mechanism, with correlation coefficient and instantaneous frequency ensuring that the output IMF components simultaneously meet waveform matching and physiological frequency band requirements, avoiding misselection problems that may be caused by a single indicator criterion. Optimized convergence performance: By discarding iteration results that do not meet the conditions and dynamically reducing the step size, the algorithm gradually approaches the local optimum that meets the constraints while maintaining its stability, preventing oscillations or divergences caused by a fixed step size. Enhanced feature extraction specificity: The second set of output IMF components can more accurately reflect the time-frequency characteristics of the QRS complex, providing a high-quality signal basis for heart rate variability analysis or arrhythmia detection.

[0048] Example 9: Based on Example 8, the process of regressing each intrinsic mode function component to its value at the end of the previous iteration, provided by this embodiment of the invention, specifically includes the following steps: S1061: After each iteration, determine whether the correlation coefficient flag and the instantaneous frequency flag are both met; if both are met, accept the results of this iteration and output them; if either flag is not met, calculate the degree of deviation that does not meet the conditions: for cases where the correlation coefficient is not met, calculate the ratio of the difference between the largest and second largest correlation coefficient to the largest correlation coefficient; for cases where the instantaneous frequency is not met, calculate the ratio of the total number of sampling points with instantaneous frequency values ​​below 5Hz or above 20Hz to the total number of sampling points. S1062: Use the larger of the two calculated ratio values ​​as the total deviation coefficient; determine the recovery amplitude based on the total deviation coefficient: make the number of historical iteration steps to be recovered equal to the preset base recovery steps multiplied by the total deviation coefficient; at the same time, make the decay coefficient of the update step size equal to the preset base decay coefficient multiplied by the total deviation coefficient; S1063: Restore each eigenmode function component to its historical value at the current time step, subtract the recovery step step from the current iteration step step; multiply the current update step size by the decay coefficient to obtain the new step size; use the new step size to re-execute the current iteration from the restored state.

[0049] The working principle and beneficial effects of the above technical solution are as follows: This embodiment uses a mechanism of quantifying the degree of deviation and dynamically adjusting the number of recovery steps and the update step size. Based on the degree of non-compliance of the correlation coefficient and instantaneous frequency, it automatically determines the historical number of recovery steps and the step size decay amplitude, ensuring that the algorithm can back down with an appropriate amplitude when the update is frustrated, avoiding the problem of decreased convergence efficiency or overcorrection caused by a fixed recovery strategy. When the total deviation coefficient is large or the constraint conditions deviate significantly, a larger number of recovery steps and a stronger step size decay are selected to quickly back down to an earlier stable state; when the deviation is small, the recovery amplitude is finely adjusted to retain some effective iteration information and improve convergence efficiency; through real-time decay of the step size and state back-off, the risk of parameter oscillation or divergence caused by unmet constraints is effectively suppressed, ensuring the numerical stability of the iteration process. When re-executing the iteration, the adjusted step size is used to start from a better historical state, prioritizing the correction of constraint terms with significant deviations, and pushing the decomposition results to approach the feasible solution that satisfies the double flag bit more quickly. This improves the robustness and convergence efficiency of the algorithm and ensures that the output components strictly meet the feature extraction requirements of the QRS composite wave.

[0050] Example 10: Based on Example 9, the process for determining the recovery amplitude based on the total deviation coefficient provided in this embodiment of the invention specifically includes the following steps: S10621: After each iteration, determine whether the correlation coefficient flag and the instantaneous frequency flag are both met; if both are met, accept the results of this iteration and output them; if either flag is not met, calculate the percentage of difference when the correlation coefficient is not met, divide the difference between the largest and second largest correlation coefficient by the largest correlation coefficient and the percentage of sampling points exceeding the limit when the instantaneous frequency is not met, divide the total number of sampling points with instantaneous frequencies below 5Hz or above 20Hz by the total number of sampling points, and take the maximum value of the two as the total deviation coefficient; S10622: Determine the number of recovery steps and the step size decay coefficient based on the total deviation coefficient: If the total deviation coefficient is less than or equal to the preset upper limit threshold, the number of recovery steps is equal to the preset basic recovery steps multiplied by the total deviation coefficient, and the step size decay coefficient is equal to the preset basic decay coefficient multiplied by the total deviation coefficient; if the total deviation coefficient is greater than the preset upper limit threshold, the number of recovery steps is directly set to the preset maximum number of recovery steps, and the step size decay coefficient is set to the square of the basic decay coefficient.

[0051] The working principle and beneficial effects of the above technical solution are as follows: This embodiment realizes adaptive convergence control under multi-objective collaborative optimization, which significantly improves the robustness, efficiency and reliability of the iterative process.

[0052] Example 11: As Figure 5 As shown, based on Example 1, the waveform following degree classification mechanism constructed based on correlation coefficient provided in this embodiment of the invention specifically includes the following steps: S201: Compare the amplitudes of all sampling points on the intrinsic mode function component waveform from left to right with the amplitudes of the corresponding positions in the noisy ECG signal: if the amplitude of the component is greater than the amplitude of the original signal, record a high-order flag; if it is less than, record a low-order flag; if they are equal, record an equal-order flag; connect all the flags in the order of the sampling points into a flag string; then traverse the flag string and count the number of flags contained in the longest segment of consecutive identical flags, and use the number of flags as the waveform following degree of the intrinsic mode function component; S202: Place the waveform follower values ​​of all intrinsic mode function components into a variable array. Repeat the following operations: Take one value from the leftmost end of the array and another value from the rightmost end, compare the two values, move the larger value to the left end of the array and occupy its left position, move the smaller value to the right end of the array and occupy its right position, then delete the old positions at both ends. After each comparison and movement, the length of the array decreases by two entries; stop repeating the operation when the array length decreases to the point where it is no longer possible to take a value from each end for pairwise comparison; determine the threshold based on the number of remaining entries in the array. If there is one remaining entry in the array: then the value of the entry is used as the first threshold, and the value is also used as the last threshold; If there are two remaining entries in the array, the smaller value of the two entries will be used as the front threshold, and the larger value will be used as the back threshold. If there are zero remaining entries in the array, meaning the initial length is even and all entries are processed in pairs, it indicates that all waveform follow-through values ​​have been sorted. Tracing back the entire movement process, we have actually obtained an ascending order of all values. We take the value at the 25th percentile or lower quartile in the ascending order as the front threshold and the value at the 75th percentile or upper quartile as the back threshold. If the quartiles cannot be accurately calculated due to the small amount of data, we take the average of the two middle values ​​after sorting as the common value of the front and back thresholds, i.e., the median.

[0053] S203: Extract the waveform tracking degree of each component sequentially, and compare the waveform tracking degree with the front threshold and the back threshold: If the waveform following degree is less than the first threshold, the intrinsic mode function component is included in the noise carrying list; if the waveform following degree is greater than the second threshold, the intrinsic mode function component is included in the signal carrying list.

[0054] If the waveform following degree is between the front threshold and the back threshold (including being equal to either threshold), or if the front threshold is equal to the back threshold (e.g., if there are 0 remaining entries and the median is taken), then calculate the center frequency of this component and the median value of the center frequencies of all intrinsic mode function components. Arrange all center frequencies in ascending order of value, and take the value in the middle of the arrangement. If there are two values ​​in the middle position, take the average of these two values. If the center frequency of the intrinsic mode function component is greater than the median value, then it is added to the noise carrying list; otherwise, it is added to the signal carrying list. Output the noise carrying list and the signal carrying list.

[0055] The working principle and beneficial effects of the above technical solution are as follows: This embodiment establishes a three-state quantization mechanism for waveform tracking degree (high / low / equal position). By statistically analyzing the longest identical sequence of continuous flag strings, the dynamic tracking characteristics of modal components on the original signal are effectively characterized. A dual-end comparison and sorting algorithm is used to construct a dynamic threshold decision boundary, avoiding modal classification bias caused by fixed thresholds and realizing automatic hierarchical identification of noise and effective signals. A median frequency reference system is introduced to process samples in the threshold interval. The relative sorting relationship of the center frequency is used to supplement the energy feature discrimination basis, improving the completeness of modal classification. The output structure forms a binary container for noise components and signal components, providing a directly usable component separation scheme for denoising reconstruction or feature extraction. This embodiment realizes a closed-loop processing from waveform tracking analysis to modal classification, taking into account both local waveform similarity characteristics and global frequency distribution characteristics, and has robust adaptability in non-stationary signal processing.

[0056] Example 12: As Figure 6As shown, based on Example 1, the process for obtaining a high-fidelity denoised ECG signal provided by this embodiment of the invention specifically includes the following steps: S301: Align the denoised noise dominant component set with the retained signal dominant component set according to the same index position in the original decomposition order, and take out one denoised noise dominant component and one signal dominant component; take the position with the same sampling time in the two components as the alignment point, and read the amplitude field of the signal dominant component and the amplitude field of the noise dominant component at each alignment point to form a dual amplitude input pair. S302: For each pair of dual-amplitude inputs, compare the absolute values ​​of the two amplitudes: if the absolute value of the dominant signal component is greater than or equal to the absolute value of the dominant noise component, then select the amplitude of the dominant signal component as the base value; otherwise, select the amplitude of the dominant noise component as the base value; compare the sign bits of the two amplitudes: if the sign bits are the same, directly output the base value as the candidate reconstructed value of the sampling point; if the sign bits are different, invert the sign bit of the base value and output it as the candidate reconstructed value; store the candidate reconstructed values ​​of all sampling points in the candidate array in chronological order; S303: Starting from the second sampling point to the penultimate sampling point on the candidate array, for each position, extract the previous, current, and next candidate reconstructed values, compare the absolute values ​​of the three values, find the value with the middle absolute value, and then extract the sign bit of the value; combine the sign bit with the middle absolute value to form a new value, and use the new value to replace the original value of the current position in the candidate array; after all positions have been processed, output the candidate array as a high-fidelity denoised ECG signal.

[0057] The working principle and beneficial effects of the above technical solution are as follows: This embodiment establishes reconstruction rules based on amplitude comparison and sign bit verification through a dual-component aligned sampling mechanism, eliminating out-of-sign interference components while preserving effective signal strength. A two-stage amplitude selection strategy is adopted to achieve dynamic balance control of noise energy and signal energy in the time domain, avoiding amplitude distortion caused by a single component. Median filtering is introduced to process candidate sequences, and isolated abrupt changes are suppressed through neighborhood amplitude sorting, maintaining the smooth continuity of the signal waveform at the microscopic time scale. A complete denoising-reconstruction processing chain is formed, and the output signal simultaneously possesses high-frequency noise suppression capability and local feature preservation capability, meeting the fidelity requirements of clinical ECG signal analysis. This embodiment takes into account both time-domain amplitude correction and local smoothing optimization, effectively preserving key waveform features such as the QRS complex while eliminating random noise.

[0058] Computer-readable storage medium according to embodiments of the present invention.

[0059] Instructions, such as computer-readable instructions, are stored on a non-transitory computer-readable storage medium. When the computer-readable instructions are executed by a processor, the various methods described above can be performed. The non-transitory computer-readable storage medium includes, but is not limited to, volatile memory and / or non-volatile memory. Volatile memory may include, for example, random access memory (RAM) and / or cache memory. Non-transitory non-volatile memory may include, for example, read-only memory (ROM), hard disk, flash memory, etc. For example, the non-transitory computer-readable storage medium can be connected to a computing device such as a computer, and then, when the computing device executes the computer-readable instructions stored on the non-transitory computer-readable storage medium, the various methods described above can be performed.

[0060] Obviously, those skilled in the art can make various modifications and variations to this invention without departing from its spirit and scope. Therefore, if these modifications and variations fall within the scope of equivalents of this invention, this invention is also intended to include these modifications and variations.

Claims

1. An adaptive cascaded ECG signal denoising method for electromyographic interference, characterized in that, Includes the following steps: The input noisy ECG signal is optimized using the Grey Wolf optimization algorithm to obtain the optimal parameter combination. This optimal parameter combination drives variational mode decomposition (VMD) to perform mode decomposition on the noisy ECG signal, resulting in a first set of several intrinsic mode function (EMF) components. Local templates of the QRS complex wave, obtained through peak detection and short-time averaging, are pre-extracted from the noisy ECG signal. The time-frequency domain characteristics of these templates, including the main peak position, bandwidth, and energy distribution, are transformed into constraints for variational decomposition. The constrained variational mode decomposition is solved to obtain a second set of EMF components. The quality of the first and second sets of EMF components is evaluated, and the set with the better quality is selected. Calculate the correlation coefficients of the intrinsic mode function components, and filter out the noise-dominated intrinsic mode function components and the signal-dominated intrinsic mode function components based on the correlation coefficient threshold. The noise-dominant intrinsic mode function (IMF) components selected were denoised using the wavelet thresholding method. The denoised noise-dominant IMF components were then combined with the retained signal-dominant IMF components to reconstruct the signal, resulting in a high-fidelity denoised ECG signal.

2. The adaptive cascaded ECG signal denoising method for electromyographic interference as described in claim 1, characterized in that, The process of obtaining the first set of several eigenmode function components includes the following steps: The noisy ECG signal is input into the gray wolf optimization algorithm, and the output signal-to-noise ratio is used as the optimization objective. The position of the gray wolf is iteratively updated within the preset search range of the decomposition layer and the penalty factor. After multiple iterations, the optimal parameter combination consisting of the optimal decomposition layer and the optimal penalty factor is obtained. The optimal decomposition level and the optimal penalty factor in the optimal parameter combination are set as the decomposition level parameter and penalty factor parameter of variational mode decomposition, respectively, to drive variational mode decomposition to perform mode decomposition on the input noisy ECG signal. Variational mode decomposition outputs multiple intrinsic mode function components, each of which corresponds to a component in a noisy electrocardiogram signal with different frequency bands.

3. The adaptive cascaded ECG signal denoising method for electromyographic interference as described in claim 2, characterized in that, The process of obtaining the optimal parameter combination consisting of the optimal number of decomposition levels and the optimal penalty factor after multiple iterations includes the following steps: The intensity of high-frequency disturbance is determined based on the frequency domain amplitude change rate of the noisy electrocardiogram signal. The preset search range of the decomposition layer is divided into three sub-intervals based on the intensity of the high-frequency disturbance. A value is then extracted from each sub-interval. The three values ​​are combined with the three preset fixed values ​​of the penalty factor to form three sets of initial candidate parameter pairs. Each pair corresponds to the initial position coordinates of a gray wolf individual. The decomposition level corresponding to the location coordinates of each individual gray wolf and the penalty factor are input into variational mode decomposition to decompose the noisy ECG signal and obtain all intrinsic mode function components under the optimal parameter combination. The difference between the center frequencies of two adjacent intrinsic mode function components is calculated. If the difference is less than the preset proportional threshold of the lower frequency of the two adjacent components, it is determined that the optimal parameter combination produces mode overlap. The decomposition level K is reduced by one step, with the preset step size being 1. At the same time, the penalty factor α is increased by one step, with the preset step size being 10% of the current α value. The updated location coordinates are obtained by increasing the penalty factor. Repeatedly execute the variational mode decomposition, inputting the decomposition level corresponding to the position coordinates of each individual gray wolf and the penalty factor, until the position coordinates of all individual gray wolves have not moved in two consecutive iterations; The output is a numerical combination of the optimal number of decomposition levels and the optimal penalty factor.

4. The adaptive cascaded ECG signal denoising method for electromyographic interference as described in claim 3, characterized in that, The process of calculating the difference in center frequencies between two adjacent intrinsic mode function components includes the following steps: All intrinsic mode function (IMF) components output by variational mode decomposition are input into a zero-crossing detector. The detector detects all zero-crossing points in the waveform of each IMF component that change from negative to positive. The time span between adjacent zero-crossing points is recorded, and the reciprocal of the time span is calculated as an instantaneous frequency value. All instantaneous frequency values ​​corresponding to each component are sorted in ascending order, and the value at the middle position is taken as the center frequency of the IMF component. The center frequencies of all components are then arranged in ascending order to obtain a center frequency sequence. Take two adjacent center frequencies sequentially from the center frequency sequence, and divide the value of the next center frequency in the center frequency sequence by the value of the previous center frequency to obtain the frequency ratio of adjacent intrinsic mode function components. Then calculate the absolute value of the difference between the frequency ratio and the value 1, and output the absolute value as the frequency separation index for adjacent intrinsic mode function components; Collect the frequency separation indices of all adjacent intrinsic mode function components and calculate the arithmetic mean of the frequency separation indices; compare the frequency separation index of each pair of adjacent intrinsic mode function components with the arithmetic mean; if the frequency separation index of a certain pair of adjacent intrinsic mode function components is less than the arithmetic mean, then record the comparison result as true, indicating that the adjacent intrinsic mode function components meet the mode overlap determination condition; Otherwise, record it as false.

5. The adaptive cascaded ECG signal denoising method for electromyographic interference as described in claim 4, characterized in that, The process of obtaining the second set of intrinsic mode function components includes the following steps: The preset constraints, namely that a certain intrinsic mode function component has the highest correlation coefficient with the local template of the QRS composite wave and the instantaneous frequency of the component is limited to the range of 5-20Hz, are placed into each iteration loop of the variational mode decomposition: before the start of each iteration, an empty set of components is initialized and the correlation coefficient flag and instantaneous frequency flag of all intrinsic mode function components are marked as not meeting the standard. After iteratively updating each intrinsic mode function component and its center frequency, the correlation coefficient between each intrinsic mode function component and the local template of the QRS composite wave is calculated sequentially. The intrinsic mode function component with the largest correlation coefficient is identified. If its correlation coefficient is greater than that of all other intrinsic mode function components, its correlation coefficient flag is set to "met the standard". At the same time, the instantaneous frequency sequence of the intrinsic mode function component is calculated. If every frequency value in the instantaneous frequency sequence is greater than or equal to 5Hz and less than or equal to 20Hz, its instantaneous frequency flag is set to "met the standard". If and only if both the correlation coefficient flag and the instantaneous frequency flag have met the criteria, the result of this iteration is accepted, and all current intrinsic mode function components are output as the second set of intrinsic mode function components; otherwise, the result of this update is discarded, each intrinsic mode function component is restored to the value at the end of the previous iteration, and the update step size is reduced before re-executing the iteration.

6. The adaptive cascaded ECG signal denoising method for electromyographic interference as described in claim 5, characterized in that, The process of reverting each intrinsic mode function component to its value at the end of the previous iteration includes the following steps: After each iteration, it is determined whether both the correlation coefficient flag and the instantaneous frequency flag have met the criteria. If both have met the criteria, the result of this iteration is accepted and output. If either flag is not met, the degree of deviation that does not meet the criteria is calculated: for cases where the correlation coefficient is not met, the ratio of the difference between the largest and second-largest correlation coefficients to the largest correlation coefficient is calculated; for cases where the instantaneous frequency is not met, the ratio of the total number of sampling points with instantaneous frequency values ​​below 5Hz or above 20Hz to the total number of sampling points is calculated. The larger of the two calculated ratios is taken as the total deviation coefficient; the recovery amplitude is determined based on the total deviation coefficient: the number of historical iteration steps to be recovered is equal to the preset base recovery steps multiplied by the total deviation coefficient; at the same time, the decay coefficient of the update step size is equal to the preset base decay coefficient multiplied by the total deviation coefficient; Restore each eigenmode function component to its historical value at time step, and subtract the number of restoration steps from the current iteration step; Multiply the current update step size by the decay factor to obtain the new step size; use the new step size to re-execute the current iteration from the recovered state.

7. The adaptive cascaded ECG signal denoising method for electromyographic interference as described in claim 6, characterized in that, The process of determining the recovery magnitude based on the total deviation coefficient includes the following steps: After each iteration, it is determined whether the correlation coefficient flag and the instantaneous frequency flag are both met. If both are met, the result of this iteration is accepted and output. If either flag is not met, the percentage of difference when the correlation coefficient is not met is calculated. The difference between the largest and second largest correlation coefficients is divided by the largest correlation coefficient, and the percentage of sampling points exceeding the limit when the instantaneous frequency is not met is calculated. The total number of sampling points with instantaneous frequencies below 5Hz or above 20Hz is divided by the total number of sampling points, and the maximum value of the two is taken as the total deviation coefficient. The recovery steps and step size decay coefficient are determined based on the total deviation coefficient: if the total deviation coefficient is less than or equal to the preset upper limit threshold, then the recovery steps are equal to the preset base recovery steps multiplied by the total deviation coefficient, and the step size decay coefficient is equal to the preset base decay coefficient multiplied by the total deviation coefficient. If the total deviation coefficient is greater than the preset upper limit threshold, the number of recovery steps is directly set to the preset maximum number of recovery steps, and the step size decay coefficient is set to the square of the basic decay coefficient.

8. The adaptive cascaded ECG signal denoising method for electromyographic interference as described in claim 1, characterized in that, The process of constructing a waveform tracking degree classification mechanism based on correlation coefficients includes the following steps: The amplitudes of all sampling points on the intrinsic mode function component waveform are compared sequentially from left to right with the amplitudes of the corresponding positions in the noisy ECG signal. If the amplitude of the component is greater than the amplitude of the original signal, a high-order flag is recorded; if it is less, a low-order flag is recorded; if they are equal, an equal-order flag is recorded. All flags are concatenated into a flag string according to the sampling point order. Then, the flag string is traversed, and the number of flags contained in the longest segment of consecutive identical flags is counted. The number of flags is used as the waveform following degree of the intrinsic mode function component. The waveform follow-up values ​​of all intrinsic mode function components are sequentially placed into a variable array; The waveform following degree of each component is extracted sequentially, and then compared with the front threshold and the back threshold.

9. The adaptive cascaded ECG signal denoising method for electromyographic interference as described in claim 8, characterized in that, The process of obtaining a high-fidelity, denoised ECG signal includes the following steps: Align the denoised noise dominant component set with the retained signal dominant component set according to the same index position in the original decomposition order, and extract one denoised noise dominant component and one signal dominant component. The positions of the two components with the same sampling time are used as alignment points. At each alignment point, the amplitude field of the dominant signal component and the amplitude field of the dominant noise component are read to form a dual amplitude input pair. For each pair of two amplitude inputs, compare the absolute values ​​of the two amplitudes: store the candidate reconstructed values ​​of all sampling points into a candidate array in chronological order; Starting from the second sampling point to the penultimate sampling point in the candidate array, for each position, extract the three candidate reconstructed values: the previous, current, and next values. Compare the absolute values ​​of the three values, find the value with the middle absolute value, and then extract the sign bit of the value. Combine the sign bit with the middle absolute value to form a new value, and use the new value to replace the original value at the current position in the candidate array. After all positions have been processed, output the candidate array as a high-fidelity denoised ECG signal.

10. The adaptive cascaded ECG signal denoising method for electromyographic interference as described in claim 9, characterized in that, If the absolute value of the dominant signal component is greater than or equal to the absolute value of the dominant noise component, the amplitude of the dominant signal component is selected as the base value; otherwise, the amplitude of the dominant noise component is selected as the base value. The sign bits of the two amplitudes are compared: if the sign bits are the same, the base value is directly output as the candidate reconstructed value of the sampling point; if the sign bits are different, the sign bit of the base value is inverted and output as the candidate reconstructed value.

Images