A wavelet denoising method and device based on kernel density estimation and a storage medium

Abstract

The application provides a wavelet denoising method and device based on kernel density estimation and a storage medium, and the method comprises the following steps: firstly, wavelet decomposition is performed on a noisy signal, the kurtosis of each layer coefficient is calculated, and the suspected noise corresponding coefficient is screened; the actual probability density distribution of the noise is fitted through kernel density estimation, the threshold value of the 99.7% confidence interval is solved, and the pulse corresponding coefficient is extracted in combination with the chi-square energy window method; and the denoising signal is obtained through coefficient zeroing and signal reconstruction. The application discards the traditional assumption that the noise conforms to the normal distribution, adapts to various noise scenes, improves the threshold estimation accuracy through mechanisms such as kurtosis interval dynamic adjustment and bandwidth optimization, and avoids the loss of pulse edges and the amplitude attenuation. The method has high pulse waveform fidelity and good noise suppression effect, and is suitable for power equipment monitoring, human physiological signal detection and the like.

Description

Technical Field

[0001] This invention belongs to the field of signal processing technology, specifically relating to a wavelet denoising method, device, and storage medium based on kernel density estimation. Background Technology

[0002] Wavelet thresholding denoising is one of the mainstream techniques for suppressing Gaussian white noise. Its core lies in accurately distinguishing between "noise-corresponding coefficients" and "signal-corresponding coefficients" in wavelet coefficients through reasonable threshold calculation. Existing wavelet denoising techniques (such as Chinese patent CN111723677A, a wavelet denoising method based on adaptive thresholding) mainly design threshold calculation logic based on the statistical characteristics of Gaussian white noise. For example, pure noise coefficients are screened by kurtosis, and then the threshold is determined by combining the 3σ criterion, or the impulse boundary identification is optimized by the chi-square energy window method.

[0003] However, in practical applications, noise often exhibits non-Gaussian and non-stationary characteristics. For example, in power equipment monitoring, in addition to Gaussian white noise, there are periodic noises such as power frequency interference and electromagnetic pulse noise such as spike noise; in human physiological signal detection, non-Gaussian distributed noise such as electromyography interference and time-varying noise such as motion artifacts are easily mixed in. Existing technologies, based on the assumption that "noise follows a normal distribution," use a fixed 3σ criterion to calculate the threshold, which has the following shortcomings:

[0004] 1) Standard deviation estimation bias of non-Gaussian noise: When the noise is non-Gaussian distributed, its amplitude distribution deviates from the normal distribution (e.g., the kurtosis of spike noise is much greater than 3). The threshold calculated directly using the 3σ criterion will over-suppress the signal or leave a large amount of residual noise.

[0005] 2) Insufficient robustness of noise figure screening: Existing technologies only screen pure noise figures using a single kurtosis threshold (e.g., [2, 5]). When the kurtosis of noise overlaps with that of small-amplitude pulse signals (e.g., the kurtosis of small-amplitude pulses is close to 3), the pulse figure is easily misjudged as a noise figure, leading to σ estimation bias.

[0006] 3) Limitations in threshold adaptability: A globally fixed threshold cannot adapt to coefficients of different frequency bands and different local signal-to-noise ratios. For example, the high-frequency layer coefficients have a high noise ratio and the low-frequency layer coefficients have a high signal ratio, resulting in an uneven denoising effect. Summary of the Invention

[0007] The first objective of this invention is to provide a wavelet denoising method based on kernel density estimation to address the aforementioned problems.

[0008] Therefore, the above-mentioned objective of the present invention is achieved through the following technical solution:

[0009] A wavelet denoising method based on kernel density estimation includes the following steps:

[0010] S1. Perform discrete wavelet decomposition on the original sampled pulse signal to obtain the detail coefficients and approximation coefficients of each layer;

[0011] S2. Calculate the kurtosis of each coefficient in the detail coefficients and approximation coefficients of each layer using a sliding window;

[0012] S3. Coefficients whose kurtosis values ​​fall within the preset threshold range [2,5] are identified as suspected noise coefficients;

[0013] S4. Based on the suspected noise correspondence coefficients, kernel density estimation is used to fit the noise probability density distribution. ;

[0014] S5, Based on noise probability density distribution Solve The integral equation yields a threshold thr with a confidence level of 99.7%, which serves as the final denoising threshold for the corresponding layer detail coefficients and approximation coefficients.

[0015] S6. Combine the final denoising threshold thr with the chi-square energy window method to extract the pulse correspondence coefficients;

[0016] S7. Set all coefficients except the pulse corresponding coefficients in the detail coefficients and approximation coefficients of each layer to zero to obtain the processed detail coefficients and approximation coefficients of each layer;

[0017] S8. The denoised signal is reconstructed using the processed detail coefficients and approximation coefficients of each layer.

[0018] While adopting the above technical solutions, the present invention may also adopt or combine the following technical solutions:

[0019] As a preferred technical solution of the present invention: in step S1, the detail coefficients and approximation coefficients of each layer are pre-processed with boundary extension, the boundary extension adopts the mirror extension method, and the extension length is 1 / 2 of the sliding window length N in step S2.

[0020] As a preferred embodiment of the present invention: in step S2, the kurtosis calculation formula is:

[0021]

[0022] In the formula, x i (i=1,2,⋯,N) is a local coefficient sequence centered on the coefficient of the kurtosis to be calculated, N is the sequence length and satisfies N< half-wave length of the pulse signal, μ is the mean of the local coefficient sequence within the sliding window, and σ is the standard deviation of the local coefficient sequence within the sliding window.

[0023] As a preferred technical solution of the present invention: In step S3, a noise coefficient quantity verification mechanism is added: If the number of noise coefficients corresponding to the initial screening n < α × L, where L is the total length of the corresponding layer coefficients, α is the proportional threshold and α ∈ [0.05, 0.1], then the upper limit of the preset threshold range of kurtosis is expanded and the maximum is not more than 8, and step S3 is re-executed until n ≥ α × L.

[0024] As a preferred technical solution of the present invention, step S4 specifically includes the following steps:

[0025] S41. Choose the Gaussian kernel function as the basis function for kernel density estimation. The expression for the Gaussian kernel function is:

[0026]

[0027] In the formula, e is the base of the natural logarithm, and π is the constant of pi.

[0028] S42. Determine the optimal bandwidth h using cross-validation to minimize the mean square error of kernel density estimation;

[0029] S43. Based on the noise corresponding coefficient set {x1, x2, ..., x...} initially selected in step S3. n The kernel density estimation function is calculated using the following formula:

[0030]

[0031] In the formula, n is the total number of suspected noise coefficients obtained from step S3, h is the optimal bandwidth for kernel density estimation, and x i Let K(*) be the magnitude of the coefficient corresponding to the i-th suspected noise, and K(*) be the Gaussian kernel function.

[0032] As a preferred technical solution of the present invention: in step S5, the integral equation is solved by numerical integration method, which is either the trapezoidal integration method or the Simpson integration method, and the integration step size is set to 1 / 10 of the amplitude resolution of the corresponding layer coefficient.

[0033] As a preferred technical solution of the present invention, step S6 specifically includes the following steps:

[0034] S61. Set the width M of the chi-square energy window according to the lower limit of the width of the target pulse signal;

[0035] S62. Compare the absolute value of the coefficients point by point with the final denoising threshold thr, and mark the coefficients x(n) with |x(n)|>thr as candidates for impulse corresponding coefficients;

[0036] S63. Locating the pulse starting point by sliding the energy window forward: Using the candidate pulse corresponding coefficient x(n) as the endpoint, construct a sliding energy window of length M, slide it forward point by point, and calculate the nominal cumulative energy E of the energy window:

[0037]

[0038] Where σ noise is the standard deviation of the noise corresponding coefficients preliminarily screened in step S3, n is the total number of suspected noise corresponding coefficients screened in step S3, and x(i) is the amplitude of the i-th suspected noise corresponding coefficient; when E < energy threshold e_thr, the coefficient corresponding to the end point of the current energy window is the starting point of the pulse; the energy threshold e_thr is determined by querying the chi-square distribution table according to the energy window width M;

[0039] S64. Slide the energy window backward to locate the end point of the pulse: Starting from the candidate x(n) of the pulse corresponding coefficient, construct a sliding energy window with a length of M, slide backward point by point and calculate E according to the formula in step S63; when E < e_thr, the coefficient corresponding to the starting point of the current energy window is the end point of the pulse; all the coefficients between the starting point and the end point of the pulse are determined as the final pulse corresponding coefficients;

[0040] S65. Starting from the next coefficient that is M coefficients away from the end point of the pulse, repeat steps S62 to S64 until all coefficients are traversed.

[0041] The second object of the present invention is to provide an electronic device, which includes a processor, a communication interface, a memory, and a communication bus. The processor, the communication interface, and the memory complete mutual communication through the communication bus.

[0042] Memory, which is used to store computer programs;

[0043] Processor, which is used to execute the computer programs stored on the memory to implement the steps of the wavelet denoising method based on kernel density estimation as described above.

[0044] The third object of the present invention is to provide a non-volatile storage medium, in which an executable program is stored. When the executable program is executed by a processor, the steps of the wavelet denoising method based on kernel density estimation as described above are implemented.

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

[0046] 1) Strong adaptability to non-Gaussian noise: By fitting the actual probability density distribution of noise through kernel density estimation and replacing the traditional normal distribution assumption, the problem of threshold estimation deviation in non-Gaussian noise scenarios is solved;

[0047] 2) High robustness in noise coefficient screening: A kurtosis interval dynamic adjustment and quantity verification mechanism is added to ensure that the sample size of noise coefficients meets the requirements of density estimation and avoid misjudging small-amplitude pulses as noise;

[0048] 3) Wide adaptability to various scenarios: It is applicable to multiple scenarios such as power equipment monitoring and human physiological signal detection. There is no need to adjust the core parameters for different scenarios. Only the wavelet basis and the number of decomposition layers need to be optimized, which reduces the complexity of engineering applications. Attached Figure Description

[0049] Figure 1 The flowchart shows the wavelet denoising method based on kernel density estimation provided by this invention.

[0050] Figure 2 This is a schematic diagram of the original sampled signal for power insulation monitoring that needs to be denoised.

[0051] Figure 3 This is a schematic diagram of the denoising results using the general threshold hard thresholding method.

[0052] Figure 4 This is a schematic diagram of the denoising results using the 3σ sliding window method.

[0053] Figure 5 This is a schematic diagram of the denoising result of the method of the present invention. Detailed Implementation

[0054] The present invention will now be described in further detail with reference to the accompanying drawings and specific embodiments.

[0055] like Figure 1 As shown, a wavelet denoising method based on kernel density estimation specifically includes the following steps:

[0056] S1. Wavelet decomposition of noisy signal: Discrete wavelet decomposition is performed on the original sampled pulse signal to obtain the detail coefficients and approximation coefficients of each level. By optimizing the selection of wavelet basis and the number of decomposition levels, the signal and noise separation of the decomposed coefficients is ensured.

[0057] Step S1, the core of discrete wavelet decomposition of the original sampled pulse signal, is to select a wavelet basis and decomposition level that are suitable for the target signal, ensuring effective separation of signal and noise in coefficients of different frequency bands. The specific method is as follows:

[0058] 1. Wavelet Basis Function Selection: Using "maximizing the cross-correlation coefficient between the mother wavelet and the target signal" as the screening criterion, the optimal basis is selected from the db series, sym series, and coif series. For example, for the "double exponential oscillating decay pulse" caused by insulation damage in power equipment, the db8 wavelet is the preferred basis due to its high waveform similarity and good time-domain locality; for human electrocardiogram signals, the sym5 wavelet is the preferred basis due to its excellent smoothness, which can reduce signal distortion.

[0059] 2. Determining the number of decomposition layers: A method based on "effective cumulative energy in the frequency domain" is adopted to avoid over-decomposition or under-decomposition. Specifically, this includes:

[0060] 1) Perform a Fourier transform on the noisy signal and calculate its frequency domain energy distribution E(f);

[0061] 2) According to frequency from the upper limit f m Calculate the reverse normalized cumulative energy up to 0 ;

[0062] 3) Set the energy ratio threshold TER=5%, and define the frequency corresponding to C(f)=TER as the effective upper limit frequency f. max Energy above this frequency is mainly noise. The frequency corresponding to C(f) = 1 - TER is defined as the effective lower limit frequency f. min Energy below this frequency is negligible;

[0063] 4) According to the formula Calculate the maximum number of decomposition layers J, where F s The signal sampling frequency, This indicates rounding down to the nearest integer.

[0064] In addition, to avoid the distortion of amplitude at both ends of the coefficients due to the boundary effect of wavelet decomposition, boundary extension processing is performed on the detail coefficients and approximation coefficients of each layer in advance. The boundary extension adopts the mirror extension method, and the extension length is 1 / 2 of the sliding window length N in step S2 to ensure that the kurtosis calculation at both ends of the coefficients is accurate.

[0065] S2. Coefficient kurtosis calculation: The kurtosis of each coefficient in each layer is calculated point by point. The kurtosis reflects the characteristic of the amplitude distribution deviating from normality, which provides a basis for subsequent noise figure screening. The kurtosis of each coefficient in the detail coefficient and approximation coefficient of each layer is calculated by using a sliding window.

[0066] In step S2, kurtosis calculation, kurtosis is a statistic reflecting the steepness of the amplitude distribution of a random variable. The kurtosis of a normal distribution is K=3. Because the amplitude of a pulse signal is concentrated in a local area, the kurtosis K>3. The kurtosis of spike noise may be much greater than 3 or less than 3. This step calculates the kurtosis coefficient by coefficient, providing a quantitative basis for distinguishing noise from signals. The kurtosis calculation formula is:

[0067]

[0068] In the formula, x i (i=1,2,⋯,N) is a local coefficient sequence centered on the coefficient of the kurtosis to be calculated, N is the sequence length and satisfies N< half-wave length of the pulse signal, μ is the mean of the local coefficient sequence within the sliding window, and σ is the standard deviation of the local coefficient sequence within the sliding window.

[0069] For binary wavelet decomposition, since the bandwidth of the corresponding frequency bands decreases by a factor of 2 for each layer, the same pulse is more dispersed in the high-frequency layer coefficients and more concentrated in the low-frequency layer. Therefore, the N value of the highest-level wavelet coefficients is equal to that of the highest-level scale coefficients, and the N value of the second-highest-level wavelet coefficients is twice that of the highest-level coefficients. This process is repeated layer by layer to ensure that the N value of each layer matches the distribution length of the pulse in that layer. For wavelet packet decomposition, the bandwidth of the corresponding frequency bands for each layer is consistent, so the N value remains the same for all layers.

[0070] S3. Screening of suspected noise coefficients: Based on the preset threshold range of kurtosis, suspected noise coefficients are screened out. At the same time, a quantity verification mechanism is added to ensure that the sample size of suspected noise coefficients meets the requirements of subsequent density estimation. Coefficients whose kurtosis values ​​fall within the preset threshold range [2,5] are identified as suspected noise coefficients.

[0071] Add a noise figure quantity verification mechanism: Step S3 suspected noise corresponding coefficient screening is based on the statistical characteristics of kurtosis. Set a preset kurtosis threshold range [2,5]. The coefficients whose kurtosis falls within this range are initially judged as "suspected noise corresponding coefficients". Since the kurtosis of Gaussian white noise is K=3, the kurtosis of slight non-Gaussian noise (such as low-intensity power frequency interference) is usually in the range of [2,5], while the kurtosis of pulse signals is mostly greater than 5, and the kurtosis of strong spike noise is mostly less than 2. To avoid insufficient sample size of suspected noise coefficients due to the special noise distribution, which could lead to bias in subsequent kernel density estimation, a quantity verification mechanism is added. Let the total length of the corresponding layer coefficients be L, and α be the proportional threshold and α∈[0.05,0.1]. If the number of noise corresponding coefficients n in the initial screening is less than α×L, then the upper limit of the preset kurtosis threshold range is expanded by 0.5 each time, with a maximum of no more than 8, and step S3 is executed again until n≥α×L. If the condition is still not met after expanding to the upper limit, then it is determined that the layer is dominated by the signal, and the threshold is directly set to 1 / 10 of the maximum value of the coefficients in that layer to avoid excessive suppression of the signal.

[0072] S4. Kernel Density Estimation and Noise Distribution Fitting: Using a Gaussian kernel function and bandwidth optimized by cross-validation, kernel density estimation is performed on the correspondence coefficients of suspected noise. The resulting kernel density estimation function is used as the actual probability density function of the noise. Based on the correspondence coefficients of suspected noise, the kernel density estimation is used to fit the noise probability density distribution. ;

[0073] S41. Kernel Function Selection: The Gaussian kernel function is chosen as the basis function for kernel density estimation. The Gaussian kernel function has good smoothness and differentiability, and its fitting error for non-Gaussian noise distributions is small. This avoids the density function discontinuity problems caused by other kernel functions, such as the rectangular kernel and the Epanechnikov kernel. The expression for the Gaussian kernel function is:

[0074]

[0075] In the formula, e is the base of the natural logarithm, and π is the constant of pi.

[0076] S42. Bandwidth Parameter Optimization: Bandwidth h is a key parameter for kernel density estimation. Too small a bandwidth leads to overfitting of the density function and the appearance of spurious peaks, while too large a bandwidth leads to underfitting of the density function and loss of distribution details. This invention determines the optimal bandwidth h through cross-validation, minimizing the mean square error of kernel density estimation and identifying the coefficient set {x1, x2, ..., x...} corresponding to suspected noise. n The data is randomly divided into 10 parts for cross-validation; each time, 9 parts are used as the training set and 1 part as the test set, for different bandwidths h, such as h∈[0.1σ]. noise ,2σ noise ], σ noise To find the standard deviation of the coefficients corresponding to suspected noise, calculate the kernel density estimation function for the training set; calculate the log-likelihood value of the test set under this density function, and select the bandwidth with the largest log-likelihood value as the optimal bandwidth h. opt ;

[0077] S43. Kernel density estimation function calculation: Based on the noise corresponding coefficient set {x1, x2, ..., x...} initially screened in step S3... n The kernel density estimation function is calculated using the following formula:

[0078]

[0079] In the formula, n is the total number of suspected noise coefficients obtained from step S3, h is the optimal bandwidth for kernel density estimation, and x i Let K(*) be the magnitude of the coefficient corresponding to the i-th suspected noise, and K(*) be the Gaussian kernel function.

[0080] Step S4: Kernel density estimation and fitting of noise distribution. Kernel density estimation is a non-parametric density estimation method that does not require a preset distribution type. It can directly fit the probability density function through sample data and is suitable for modeling the distribution of non-Gaussian noise.

[0081] S5. Confidence Interval Threshold Calculation: Based on Noise Probability Density Distribution Solve The integral equation yields a threshold thr with a confidence level of 99.7%, which replaces the traditional 3σ threshold as the final denoising threshold for the corresponding layer detail coefficients and approximation coefficients.

[0082] 1) Construct the integral equation: That is, the cumulative probability to the left of the threshold thr is 99.7%, and the cumulative probability to the right is 0.3%, which corresponds to a low-probability event being judged as a signal;

[0083] 2) Solve the equation using numerical integration: Select the trapezoidal integration method or Simpson's integration method, and set the integration step size to 1 / 10 of the amplitude resolution of the corresponding layer coefficient, where the amplitude resolution = maximum coefficient amplitude / coefficient quantization bits, to ensure the integration accuracy. For example, if the maximum coefficient amplitude is 10V and the quantization bits are 16 bits, then the amplitude resolution ≈ 0.153mV and the integration step size ≈ 0.0153mV. Take the solution result thr as the final threshold.

[0084] S6. Pulse corresponding coefficient extraction: Combine the final denoising threshold thr with the chi-square energy window method to accurately locate the pulse boundary and extract the pulse corresponding coefficients, avoiding "pulse edge loss" or "noise misjudged as signal" caused by a single threshold;

[0085] S61. Set the energy window width: Set the chi-square energy window width M according to the lower limit of the width of the target pulse signal;

[0086] S62. Identify pulse candidate coefficients: Compare the absolute value of the coefficient with the final denoising threshold thr point by point, and mark the coefficient x(n) with |x(n)| > thr as a pulse corresponding coefficient candidate, initially screening out the area that may contain pulses;

[0087] S63. Identify the pulse start point: Slide the energy window forward to locate the pulse start point: Take the pulse corresponding coefficient candidate x(n) as the end point, construct a sliding energy window with a length of M, slide forward point by point and calculate the nominal cumulative energy E of the energy window:

[0088]

[0089] In the formula, σ noise is the standard deviation of the noise corresponding coefficients preliminarily screened in step S3, n is the total number of suspected noise corresponding coefficients screened in step S3, and x(i) is the amplitude of the i-th suspected noise corresponding coefficient; when E < energy threshold e_thr, the coefficient corresponding to the end point of the current energy window is the pulse start point; the energy threshold e_thr is determined by querying the chi-square distribution table according to the energy window width M;

[0090] S64. Identify the pulse end point: Slide the energy window backward to locate the pulse end point: Take the pulse corresponding coefficient candidate x(n) as the start point, construct a sliding energy window with a length of M, slide backward point by point and calculate E according to the formula in step S63; when E < e_thr, the coefficient corresponding to the start point of the current energy window is the pulse end point; all the coefficients between the pulse start point and the pulse end point are determined as the final pulse corresponding coefficients;

[0091] S65. Traverse and extract pulses: Start from the next coefficient that is M coefficients away from the pulse end point, and repeat steps S62 to S64 until all coefficients are traversed.

[0092] S7. Zeroing Coefficients: Set the non-pulse corresponding coefficients (except for the pulse corresponding coefficients) in each layer of detail coefficients and approximation coefficients to zero to replace the traditional hard threshold / soft threshold function. This avoids signal discontinuity caused by hard thresholding or amplitude attenuation caused by soft thresholding, ensures pulse waveform fidelity, and preserves effective signal components to obtain the processed detail coefficients and approximation coefficients.

[0093] S8. Denoising signal reconstruction: Using the processed detail coefficients and approximation coefficients of each layer, the denoised signal is obtained through wavelet reconstruction.

[0094] Example

[0095] This specific embodiment takes the monitoring of insulation damage in a high-voltage motor as an example. Partial discharge commonly occurs when the insulation of a high-voltage motor is damaged. In industrial settings, high-voltage coupling capacitors or high-frequency current transformers are typically used to measure this on the high-voltage side or neutral point side of the motor. However, due to the harsh electromagnetic environment at the site, characterized by numerous background noise sources and diverse noise types, traditional wavelet threshold denoising methods are not ideal, making it difficult to accurately assess the insulation damage of the high-voltage motor. The partial discharge monitoring signal processed in this embodiment comes from the generator set of a large hydropower station. A partial discharge sensor is installed at the generator outlet to pick up the partial discharge signal. The sampling rate of the monitoring system is 40 MS / s, with each acquisition lasting 20 ms (5 power frequency cycles). This embodiment uses a single acquisition as the processing object, with a data length of 800 KB. The signal contains a "double exponential oscillating decaying pulse" (a characteristic signal of insulation damage, with a pulse width of approximately 50 μs and an amplitude of 0.1-1 V) and mixed noise (Gaussian white noise, power frequency interference, and electromagnetic spike noise).

[0096] like Figure 1 The wavelet denoising method based on kernel density estimation of the present invention, shown below, includes:

[0097] S1, for example Figure 2 The original sampled partial discharge pulse signal is shown to be subjected to discrete wavelet decomposition to obtain detail coefficients and approximation coefficients at each level. The discrete wavelet decomposition is a binary wavelet decomposition, and the mother wavelet used is the db8 wavelet, which has high similarity to the partial discharge pulse waveform. The decomposition has 6 levels, obtaining detail coefficients from level 1 to level 6 and approximation coefficients at level 6; where:

[0098] Mother wavelet determination: The mother wavelet selection is determined by calculating the cross-correlation coefficient between different mother wavelets and the actual partial discharge pulse. The mother wavelet selection range used includes db2~db40, sym2~sym20, coif1~coif5, and gaus1~gaus8. It is found that db8 has the largest cross-correlation coefficient with the partial discharge pulse, which is 0.89. Therefore, in this embodiment, db8 is used as the mother wavelet for wavelet decomposition.

[0099] Determining the Number of Decomposition Levels: When performing wavelet denoising on partial discharge signals, in addition to selecting a suitable mother wavelet, it is also necessary to determine the number of wavelet decomposition levels. If the number of decomposition levels is too large, it will increase the algorithm's time complexity and introduce spurious frequency components, reducing the denoising effect; conversely, if the number of decomposition levels is too small, it will reduce the algorithm's resolution, also reducing the denoising effect. There are many methods for determining the number of wavelet decomposition levels, the most common being the maximum decomposition level method, which decomposes the signal to the maximum number of separable levels. However, most decomposition results contain no effective signals, and excessive division wastes computational power. Therefore, this embodiment adopts a method for determining the number of decomposition levels based on the effective cumulative energy in the frequency domain. First, the reverse-order normalized cumulative energy distribution of the signal in the frequency domain is calculated; then, the energy ratio threshold is used... Identify the effective upper frequency limit of the signal from this distribution. Then recalculate. The normalized frequency domain cumulative energy distribution of the following frequency bands, and using... Identify the effective lower limit frequency of the signal Given The following frequency components contain only a small amount of energy; not further subdividing them does not affect noise reduction. Therefore, they should be eliminated first. The reason for the above frequency band is that this frequency band mainly consists of small-amplitude noise, but its frequency range is relatively wide, which has an impact. The determination of the effective lower limit frequency of the pulse signal. The maximum number of decomposition levels can be calculated using the following formula.

[0100]

[0101] In the formula, The signal sampling frequency, This represents the floor function. According to... feature, Typically set as a percentage of the total energy, in this embodiment it is set to 5%, resulting in a signal decomposition layer number of J=6 (calculated f). min =0.5MHz, The lengths of coefficients for each layer are as follows: the lengths of detail coefficients (cD1-cD6) for layers 1-6 are 400003, 200005, 100006, 50006, 25006, and 12506, respectively; the length of approximation coefficient (cA6) for layer 6 is 12506.

[0102] S2. Extend the boundary of the detail coefficients of layers 1 to 6 and the approximation coefficients of layer 6 to obtain new detail coefficients of layers 1 to 6 and approximation coefficients of layer 6;

[0103] Boundary extension is a common technique in wavelet denoising. In this specific embodiment, mirror extension is used, with the extension length being half the width of the sliding window in step S3. This allows the kurtosis of all wavelet decomposition coefficients to be calculated. The sliding window width needs to be determined based on the signal frequency band size within the detail coefficients or approximation coefficients, as the frequency band size determines the relationship between the length of the coefficients at that level after wavelet decomposition and the length of the original signal. The sliding window width should ideally be the length of half the wave of the pulse to be detected. For binary wavelet decomposition, the frequency band widths corresponding to the detail coefficients and approximation coefficients at each level are different. The frequency band widths of adjacent detail coefficients are twice that of each other, while the frequency band widths of the highest-level detail coefficients and approximation coefficients are the same. Therefore, the N used for the highest-level detail coefficients is the same as the N used for the highest-level approximation coefficients, the N used for the next higher-level detail coefficients is half the N used for the highest-level detail coefficients, and so on, to determine the N used for each level of detail coefficients.

[0104] Taking cD3 (the third layer of detail coefficients, corresponding to the frequency band 2.5-5 MHz) as an example, N=160, the kurtosis is calculated for each coefficient to obtain the kurtosis distribution. The kurtosis of the noise coefficient is concentrated in [2.2, 4.8], and the kurtosis of the impulse coefficient is concentrated in [5.5, 12].

[0105] S3. Screening of suspected noise correspondence coefficients is based on the statistical characteristics of kurtosis. A preset threshold range of kurtosis [2,5] is set. Coefficients whose kurtosis falls within this range are initially determined as "noise correspondence coefficients". Since the kurtosis of Gaussian white noise is K=3, the kurtosis of slight non-Gaussian noise (such as low-intensity power frequency interference) is usually in the range of [2,5], while the kurtosis of pulse signals is mostly greater than 5, and the kurtosis of strong spike noise is mostly less than 2. To avoid insufficient sample size of suspected noise coefficients due to the special noise distribution, which would lead to bias in subsequent kernel density estimation, a quantity verification mechanism is added: Let the total length of the corresponding layer coefficients be L, and the proportion threshold α∈[0.05,0.1], that is, the sample size of noise coefficients should be at least 5%-10% of the total length of coefficients; if the number of noise coefficients n selected in the initial screening is less than α×L, then the upper limit of the preset kurtosis threshold range is gradually expanded, increasing by 0.5 each time, with a maximum of no more than 8, and the screening is repeated until n≥α×L; if it is still not satisfied after expanding to the upper limit, then it is determined that the layer is dominated by signals, and the threshold is directly set to 1 / 10 of the maximum value of the coefficients in that layer to avoid excessive suppression of signals.

[0106] The kurtosis preset threshold range is [2,5]. The number of coefficients corresponding to suspected noise of cD3 is n=9800, the total length of coefficients is L=100006, α=0.05, α×L=5003, and n>5003, which meets the quantity requirements.

[0107] S4. Kernel Density Estimation for Noise Distribution Fitting. Kernel density estimation is a non-parametric density estimation method that does not require a pre-defined distribution type. It can directly fit the probability density function from sample data and is suitable for modeling the distribution of non-Gaussian noise. This step fits the distribution of the coefficients corresponding to suspected noise using kernel density estimation, specifically including:

[0108] 1) Kernel function selection: The Gaussian kernel function is used, and its expression is: The Gaussian kernel function has the characteristics of good smoothness and good differentiability, and has small fitting error for the distribution of non-Gaussian noise, avoiding the density function discontinuity problem caused by other kernel functions, such as the rectangular kernel and the Epanechnikov kernel.

[0109] 2) Bandwidth Parameter Optimization: Bandwidth h is a key parameter for kernel density estimation. Too small a bandwidth leads to overfitting of the density function and the appearance of spurious peaks, while too large a bandwidth leads to underfitting of the density function and loss of distribution details. This invention uses cross-validation to select the optimal bandwidth: the set of coefficients corresponding to suspected noise {x1, x2, ..., x...} is... n The data is randomly divided into K=10 parts for cross-validation; each time, 9 parts are used as the training set and 1 part as the test set, for different bandwidths h, such as h∈[0.1σ]. noise , 2σ noise ], σ noise To find the standard deviation of the coefficients corresponding to suspected noise, calculate the kernel density estimation function for the training set; calculate the log-likelihood value of the test set under this density function, and select the bandwidth with the largest log-likelihood value as the optimal bandwidth h. opt .

[0110] 3) Kernel density estimation function calculation: based on optimal bandwidth h opt The probability density function of the noise is calculated using the formula, corresponding to the set of coefficients of the suspected noise. In the formula, n is the number of coefficients corresponding to suspected noise, and x i Let K(*) be the coefficient corresponding to the i-th suspected noise, and K(*) be the Gaussian kernel function.

[0111] Kernel density estimation was performed on the suspected noise coefficients (n=9800) of cD3, and the distribution of non-Gaussian noise was found to be "bimodal" (corresponding to Gaussian white noise and power frequency interference), which was significantly different from the normal distribution curve, thus verifying the necessity of kernel density estimation.

[0112] S5. In solving the confidence interval threshold, the probability density function based on kernel density estimation is used to solve for the equivalent 99.7% confidence interval threshold thr, which is obtained using Simpson's integral method. The threshold value of cD3 was found to be thr = 0.085V, while the traditional 3σ threshold value was 0.062V (a deviation of 37%).

[0113] S6. In the extraction of the pulse correspondence coefficient, by combining the final threshold thr with the improved chi-square energy window method, the pulse boundary is accurately located to avoid "pulse edge loss" or "noise misjudged as signal" caused by a single threshold:

[0114] S61. Setting the energy window width: Set the chi-square energy window width M according to the lower limit of the width of the target pulse signal; corresponding to the lower limit of the pulse width of 50 μs, 50 μs / 25 ns = 2000, look up the chi-square distribution with a degree of freedom of 2000 and the energy threshold e_thr = 2423 corresponding to a 99.5% confidence level;

[0115] S62. Identifying the pulse candidate coefficient: Compare the absolute value of the coefficient with the final threshold thr point by point, and mark the coefficient x(n) with |x(n)| > thr as "pulse correspondence coefficient candidate", initially screening out the area that may contain pulses.

[0116] S63. Identifying the pulse start point: Slide the energy window forward to locate the pulse start point: Taking the candidate coefficient x(n) as the end point, construct a sliding energy window with a length of M, slide forward point by point and calculate the nominal cumulative energy E of the energy window:

[0117]

[0118] In the formula, σ noise is the standard deviation of the noise correspondence coefficient preliminarily screened in step S3; when E < the energy threshold e_thr, the coefficient corresponding to the end point of the current energy window is the pulse start point;

[0119] S64. Identifying the pulse end point: Slide the energy window backward to locate the pulse end point: Taking the candidate coefficient x(n) as the start point, construct a sliding energy window with a length of M, slide backward point by point and calculate E according to the formula in step S63; when E < e_thr, the coefficient corresponding to the start point of the current energy window is the pulse end point; all the coefficients between the pulse start point and the pulse end point are determined as the final pulse correspondence coefficients;

[0120] S65. Traversing to extract pulses: Starting from the next coefficient that is M coefficients away from the pulse end point, repeat steps S62 to S64 until all the coefficients in the current layer are traversed.

[0121] S7. Coefficient zeroing processing: Set to zero the remaining coefficients in each layer of wavelet coefficients and scale coefficients except the final pulse correspondence coefficients, to replace the traditional hard threshold / soft threshold function, avoid signal discontinuity caused by the hard threshold or amplitude attenuation caused by the soft threshold, and ensure the fidelity of the pulse waveform.

[0122] S8. Signal reconstruction: Perform wavelet inverse transform on the processed detail coefficients and approximation coefficients of each layer to reconstruct the denoised signal, and the result is as Figure 3-5As shown, the detail denoising effect is significantly better than the general threshold hard threshold method and the 3σ energy window method.

[0123] The present invention also provides an electronic device, which includes a processor, a communication interface, a memory, and a communication bus. The processor, the communication interface, and the memory communicate with each other via the communication bus.

[0124] Memory, used to store computer programs;

[0125] The processor executes a computer program stored in memory to implement the wavelet denoising method steps based on kernel density estimation as described above.

[0126] The present invention also provides a non-volatile storage medium storing an executable program, which, when executed by a processor, implements the wavelet denoising method steps based on kernel density estimation as described above.

[0127] The technical solution of the present invention has been described in conjunction with the specific experimental procedures shown in the accompanying drawings. However, the scope of protection of the present invention is not limited to these specific embodiments. Without departing from the principles of the present invention, those skilled in the art can make equivalent changes or substitutions to the relevant technical features, and the technical solutions resulting from such changes or substitutions will all fall within the scope of protection of the present invention.

Claims

1. A wavelet denoising method based on kernel density estimation, characterized in that, It includes the following steps: S1. Perform discrete wavelet decomposition on the original sampling pulse signal to obtain detail coefficients and approximation coefficients of each layer; S2. Use a sliding window to calculate the kurtosis of each coefficient in the detail coefficients and approximation coefficients of each layer respectively; S3. Identify the coefficients whose kurtosis values fall within a preset threshold range as the corresponding coefficients of suspected noise; S4. Based on the corresponding coefficients of suspected noise, use kernel density estimation to fit the noise probability density distribution; S5. Based on the noise probability density distribution, solve the integral equation to obtain the threshold thr as the final denoising threshold for the detail coefficients and approximation coefficients of the corresponding layer; S6. Combine the final denoising threshold thr with the chi-square energy window method to extract the corresponding coefficients of the pulse; S7. Set the remaining coefficients in the detail coefficients and approximation coefficients of each layer except the corresponding coefficients of the pulse to zero to obtain the processed detail coefficients and approximation coefficients of each layer; S8. Use the processed detail coefficients and approximation coefficients of each layer to reconstruct the denoised signal.

2. The method according to claim 1, characterized in that: In step S1, boundary extension processing is pre-executed on the detail coefficients and approximation coefficients of each layer. The boundary extension uses the mirror extension method, and the extension length is 1 / 2 of the sliding window length N in step S2.

3. The method according to claim 1, characterized in that: In step S2, the kurtosis calculation formula is: In the formula, x i (i=1,2,⋯,N) is a local coefficient sequence centered on the coefficient of the kurtosis to be calculated, N is the sequence length and satisfies N< half-wave length of the pulse signal, μ is the mean of the local coefficient sequence within the sliding window, and σ is the standard deviation of the local coefficient sequence within the sliding window.

4. The method according to claim 1, characterized in that: In step S3, a noise coefficient quantity verification mechanism is added: if the number n of initially screened corresponding coefficients of noise < α×L, where L is the total length of the coefficients of the corresponding layer and α is the proportional threshold and α∈[0.05,0.1], then increase the upper limit of the preset threshold range of kurtosis and the maximum does not exceed 8, and re-execute step S3 until n≥α×L.

5. The method according to claim 1, characterized in that: Step S4 specifically includes the following steps: S41. Select the Gaussian kernel function as the basis function for kernel density estimation. The Gaussian kernel function expression is: In the formula, e is the base of the natural logarithm, and π is the constant of the circumference ratio; S42. Determine the optimal bandwidth h through the cross-validation method to minimize the mean square error of kernel density estimation; S43. Based on the noise corresponding coefficient set {x1, x2, ..., x...} initially selected in step S3. n The kernel density estimation function is calculated using the following formula: In the formula, n is the total number of suspected noise coefficients obtained from step S3, h is the optimal bandwidth for kernel density estimation, and x i Let K(*) be the magnitude of the coefficient corresponding to the i-th suspected noise, and K(*) be the Gaussian kernel function.

6. The method according to claim 1, characterized in that: In step S5, the numerical integration method is used to solve the integral equation. The numerical integration method is the self-trapezoidal integration method or the Simpson integration method, and the integration step size is set to 1 / 10 of the amplitude resolution of the coefficients of the corresponding layer.

7. The method according to claim 1, characterized in that: Step S6 specifically includes the following steps: S61. Set the chi-square energy window width M according to the lower limit of the width of the target pulse signal; S62. Compare the absolute value of the coefficient with the final denoising threshold thr point by point, and mark the coefficient x(n) with |x(n)|>thr as a candidate for the corresponding coefficient of the pulse; S63. Slide the energy window forward to locate the pulse start point: taking the candidate x(n) of the corresponding coefficient of the pulse as the end point, construct a sliding energy window with a length of M, slide forward point by point and calculate the nominal cumulative energy E of the energy window: In the formula, σ noise Let be the standard deviation of the noise correspondence coefficients initially screened in step S3, n be the total number of suspected noise correspondence coefficients obtained in step S3, and x(i) be the amplitude of the i-th suspected noise correspondence coefficient. When E < energy threshold e_thr, the coefficient corresponding to the end point of the current energy window is the pulse start point; the energy threshold e_thr is determined by querying the chi-square distribution table according to the energy window width M; S64. Slide the energy window backward to locate the pulse end point: taking the candidate x(n) of the corresponding coefficient of the pulse as the start point, construct a sliding energy window with a length of M, slide backward point by point and calculate E according to the formula in step S63; when E < e_thr, the coefficient corresponding to the start point of the current energy window is the pulse end point; All the coefficients between the pulse start point and the pulse end point are determined as the final corresponding coefficients of the pulse; S65. Starting from the next coefficient after the pulse end interval M coefficients, repeat steps S62 to S64 until all coefficients have been traversed.

8. An electronic device, comprising a processor, a communication interface, a memory, and a communication bus, wherein the processor, the communication interface, and the memory communicate with each other via the communication bus, characterized in that: The memory is used to store computer programs; A processor for executing a computer program stored in a memory to implement the wavelet denoising method steps based on kernel density estimation as described in any one of claims 1-6.

9. A non-volatile storage medium, characterized in that: The non-volatile storage medium stores an executable program, which, when executed by a processor, implements the wavelet denoising method steps based on kernel density estimation as described in any one of claims 1-6.

Images