A signal denoising processing method and system based on wavelet transform

Abstract

The application relates to the technical field of signal processing, in particular to a signal denoising processing method and system based on wavelet transform, which comprises the following steps: performing multi-scale wavelet transform on an original noisy signal to obtain a wavelet coefficient sequence of each scale; for each scale, defining a time inversion sequence of the wavelet coefficient sequence; calculating the normalized similarity of the wavelet coefficient sequence and the time inversion sequence; calculating a causal asymmetry index according to the normalized similarity; constructing a shrinkage factor according to the causal asymmetry index, and performing weighted shrinkage on the original wavelet coefficient sequence by using the shrinkage factor to obtain a denoised wavelet coefficient sequence; and performing inverse wavelet transform on the denoised wavelet coefficient sequence to reconstruct a denoised original signal. The application solves the technical problems in the prior art that the distinction between signals and noises completely depends on the amplitude or energy characteristics of wavelet coefficients, and that weak fault impact pulses cannot be effectively extracted under a low signal-to-noise ratio condition.

Description

Technical Field

[0001] This invention relates to the field of signal processing technology, and more specifically, to a signal denoising method and system based on wavelet transform. Background Technology

[0002] Rotating machinery is core equipment in industrial production. The operating status of key transmission components such as gearboxes, rolling bearings, and wind turbine drive chains directly affects the safety and reliability of the entire equipment. When these components suffer localized damage, such as pitting on gear teeth or peeling of the inner and outer rings of bearings, periodic impact pulses are generated when meshing or rolling over the damage point. The arrival time, amplitude variation, and attenuation pattern of these impact pulses contain rich fault information, which can be used to determine the fault type, locate the damaged area, and assess the severity of the damage. However, in actual industrial field signal acquisition, effective vibration signals are often submerged by strong background noise. In the early stages of a fault, the energy of the impact pulses generated at the damage point is relatively weak, and their amplitude may even be far lower than the root mean square value of the noise, and their impact characteristics are almost completely submerged in the background noise.

[0003] In existing technologies, wavelet transform-based denoising methods are mainly used. Wavelet transform has excellent time-frequency localization capabilities, enabling the decomposition of signals into different frequency sub-bands for analysis, making it the mainstream tool in the field of mechanical vibration signal processing. There are many improved wavelet-based denoising schemes. The first type is wavelet threshold denoising. This method thresholds the wavelet coefficients. The problem with this method is that it relies solely on the amplitude of the wavelet coefficients to distinguish between signal and noise. When the noise intensity is comparable to a weak impact pulse, the amplitude of the wavelet coefficients corresponding to the impact pulse may be lower than that of the wavelet coefficients corresponding to the noise. Threshold denoising leads to the impact pulse being treated as noise and suppressed, while high-intensity noise spikes are treated as signals and preserved, making it unsuitable for low signal-to-noise ratio scenarios. The second type is denoising methods based on inter-scale correlation. This type of method utilizes the strong correlation between signals and adjacent wavelet scales, while the correlation of noise rapidly decreases with increasing scale. It distinguishes between signals and noise by calculating the product or correlation coefficient of wavelet coefficients at adjacent scales. However, when the energy of the impact pulse is weak, its cross-scale correlation also weakens, making it difficult to distinguish from noise. The third type of method is a denoising method based on statistical models. This typically involves using methods such as Hidden Markov Models (HMMs) and Bayesian estimation frameworks to model wavelet coefficients and achieve signal-to-noise separation by estimating model parameters. However, this type of method requires prior assumptions about the probability distributions of signal and noise, making computation complex, parameter estimation difficult, and the model assumptions hard to meet under complex industrial field conditions.

[0004] In summary, existing wavelet transform-based denoising methods rely entirely on the amplitude or energy characteristics of wavelet coefficients to distinguish between signals and noise. When the amplitude of the impulse pulse is lower than the noise level, weak impulse pulses are incorrectly suppressed while high-intensity noise spikes are misjudged as signals, resulting in weak fault impulse pulses being unable to be effectively extracted under low signal-to-noise ratio conditions. Summary of the Invention

[0005] The purpose of this application is to provide a signal denoising processing method and system based on wavelet transform, which solves the technical problems in the prior art, such as the fact that the distinction between signal and noise depends entirely on the amplitude or energy characteristics of wavelet coefficients, resulting in the inability to effectively extract weak fault impulse pulses under low signal-to-noise ratio conditions.

[0006] To solve the above-mentioned technical problems, the solution adopted in this application is as follows:

[0007] A signal denoising method based on wavelet transform includes:

[0008] S1: Perform multi-scale wavelet transform on the original noisy signal to obtain... A sequence of wavelet coefficients at various scales;

[0009] S2: For each scale, define the time-reversal sequence of the wavelet coefficient sequence;

[0010] S3: Calculate the normalized similarity between the wavelet coefficient sequence and the time-inverted sequence;

[0011] S4: Calculate the causal asymmetry index based on standardized similarity;

[0012] S5: Construct a shrinkage factor based on the causal asymmetry index, and use the shrinkage factor to perform weighted shrinkage on the original wavelet coefficient sequence to obtain the denoised wavelet coefficient sequence.

[0013] S6: Perform inverse wavelet transform on the denoised wavelet coefficient sequence to reconstruct the original denoised signal.

[0014] Preferably, the specific implementation method of S1 includes the following steps:

[0015] S11: Acquire the raw, noisy vibration signal sequence of rotating machinery;

[0016] S12: Determine the maximum number of wavelet transform decomposition levels based on the sampling frequency and lowest frequency of the original noisy vibration signal sequence. ;

[0017] S13: Utilizing dual-tree complex wavelet transform to process the original noisy vibration signal sequence. Layer-by-layer recursive wavelet decomposition yields the scaling coefficients for each layer.

[0018] S14: Extract wavelet coefficient sequences from the scaling coefficients of each layer;

[0019] S15: Extract the original phase information from the scale coefficients of each layer.

[0020] Preferably, the specific implementation method of S2 includes the following steps:

[0021] S21: Obtain the wavelet coefficient sequences and their lengths at each scale output from step S1;

[0022] S22: For each scale of the wavelet coefficient sequence, with the last point of the wavelet coefficient sequence as the center of symmetry, the entire sequence is time-flipped to generate a time-inverted sequence.

[0023] Preferably, the specific implementation method of S3 includes the following steps:

[0024] S31: Set a local window;

[0025] S32: At each scale, starting from the beginning position of the wavelet coefficient sequence at that scale, slide the local window point by point along the direction of the wavelet coefficient sequence at that scale, and at each local window position, extract a segment of the wavelet coefficient sequence within the local window range;

[0026] S33: At each scale, starting from the beginning position of the time inversion sequence at that scale, slide a local window point by point along the direction of the time inversion sequence at that scale, and at each local window position, extract a time inversion sequence segment within the local window range;

[0027] S34: Calculate the normalized similarity between the wavelet coefficient sequence fragment and the time-inverted sequence fragment corresponding to each local window position.

[0028] Preferably, the specific implementation method of S4 includes the following steps:

[0029] S41: Obtain the standardized similarity corresponding to each local window position at each scale output in step S3;

[0030] S42: For each local window location at each scale, subtract the standardized similarity corresponding to that local window location from 1, and use the difference as the causal asymmetry index of that local window location;

[0031] The causal asymmetry indices at different local window locations at this scale constitute a causal asymmetry index sequence;

[0032] S43: For the positions at both ends of a sequence whose standardized similarity could not be calculated due to the inability of a local window to fully cover them at each scale, their causal asymmetry index is assigned a value of 0.

[0033] Preferably, the specific implementation method of S5 includes the following steps:

[0034] S51: Obtain the causal asymmetric exponential sequence and wavelet coefficient sequence at each scale;

[0035] S52: For the index position of each wavelet coefficient at each scale, input the causal asymmetric exponential sequence of the index position into the preset shrinkage factor mapping function to obtain the shrinkage factor corresponding to the index position;

[0036] S53: For each index position at each scale, multiply the shrinkage factor corresponding to the index position by the wavelet coefficient at that position to obtain the denoised wavelet coefficient.

[0037] Preferably, the specific implementation method of S52 includes the following steps:

[0038] S521: Determine the type of basis function for the contraction factor mapping function, i.e., use the hyperbolic tangent function as the basis function for the mapping function;

[0039] S522: Introduce a gain constant, and use the product of the gain constant and the causal asymmetry index as the independent variable of the hyperbolic tangent function to obtain the complete form of the contraction factor mapping function;

[0040] S523: Determine the value of the gain constant to obtain the final shrinkage factor mapping function.

[0041] Preferably, the shrinkage factor mapping function is specifically:

[0042] ;

[0043] in, The preset gain constant, and ; For the first Scale No. The causal asymmetry index of each index position; This is the contraction factor corresponding to the index position.

[0044] Preferably, the specific implementation method of S6 includes the following steps:

[0045] S61: Obtain the denoised wavelet coefficient sequence at each scale and the original phase information at each scale;

[0046] S62: Recombine the denoised wavelet coefficient sequence with the original phase information to restore it to the denoised complex wavelet coefficients;

[0047] S63: From the Starting from the first layer, the signal is reconstructed layer by layer using multi-scale wavelet inverse transform. That is, in each layer of reconstruction, the denoised complex wavelet coefficients of the current layer and the reconstruction approximation coefficients of the previous layer are used to recover the approximation coefficients of the current layer through a synthesis filter bank until the approximation coefficients of the 0th layer are recovered.

[0048] S64: Average the approximation coefficients of layer 0 to obtain the reconstructed denoised original signal.

[0049] A signal denoising system based on wavelet transform, used to implement the aforementioned signal denoising method based on wavelet transform, includes:

[0050] The wavelet transform module is used to perform multi-scale wavelet transform on the original noisy signal to obtain... A sequence of wavelet coefficients at various scales;

[0051] The time inversion module is used to define the time inversion sequence of the wavelet coefficient sequence for each scale;

[0052] The similarity calculation module is used to calculate the standardized similarity between the wavelet coefficient sequence and the time-inverted sequence;

[0053] The causal index calculation module is used to calculate the causal asymmetry index based on standardized similarity.

[0054] The shrinkage and denoising module is used to construct a shrinkage factor based on the causal asymmetry index, and to use the shrinkage factor to perform weighted shrinkage on the original wavelet coefficient sequence to obtain the denoised wavelet coefficient sequence.

[0055] The inverse transform reconstruction module is used to perform inverse wavelet transform on the denoised wavelet coefficient sequence to reconstruct the denoised original signal.

[0056] The technical solution of this application has at least the following advantages and beneficial effects:

[0057] This invention determines the presence of a causal signal structure at a given location by comparing the differences in the variation trends of wavelet coefficient sequence segments and their time-reversed sequence segments. Since impulse pulses have a definite physical occurrence time, their rising edges undergo fundamental changes in morphology under time reversal; while noise, as a random process, retains its statistical characteristics unchanged under time reversal. Therefore, by calculating the causal asymmetry index, the consistency of the variation trends between wavelet coefficient sequence segments and their time-reversed sequence segments can be determined. A weak impulse pulse with a very small amplitude, as long as its rising edge exhibits a unidirectional increasing trend from low to high within the local window, will have a standardized similarity significantly deviating from 1, thus being assigned a higher causal asymmetry index and being fully preserved in subsequent shrinkage processing. Conversely, a noise spike with a large amplitude, because its variation trend within the local window is non-directional, has a standardized similarity close to 1, and a causal asymmetry index close to 0, thus being effectively suppressed in subsequent shrinkage processing. Furthermore, this invention maps the causal asymmetry index to a contraction factor using a hyperbolic tangent function. The causal asymmetry index itself is adaptively calculated from the local morphological features of the signal, without requiring prior knowledge of noise, thus solving the problems in the background technology. Attached Figure Description

[0058] Figure 1 This is a flowchart of the method of the present invention. Detailed Implementation

[0059] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0060] Please refer to Figure 1 This invention provides a signal denoising processing method based on wavelet transform, comprising:

[0061] S1: Perform multi-scale wavelet transform on the original noisy signal to obtain... A sequence of wavelet coefficients at various scales;

[0062] S2: For each scale, define the time-reversal sequence of the wavelet coefficient sequence;

[0063] S3: Calculate the normalized similarity between the wavelet coefficient sequence and the time-inverted sequence;

[0064] S4: Calculate the causal asymmetry index based on standardized similarity;

[0065] S5: Construct a shrinkage factor based on the causal asymmetry index, and use the shrinkage factor to perform weighted shrinkage on the original wavelet coefficient sequence to obtain the denoised wavelet coefficient sequence.

[0066] S6: Perform inverse wavelet transform on the denoised wavelet coefficient sequence to reconstruct the original denoised signal.

[0067] Rotating machinery, such as gearboxes and rolling bearings, generates periodic, weak impact pulses when critical components suffer localized damage during operation. The arrival time, amplitude, and attenuation pattern of these impact pulses are typically used to determine the type and severity of the fault. However, when acquiring these impact pulse signals, the effective signal is often drowned out by background noise. In the early stages of component damage, the amplitude of the impact pulses can even be far lower than the root mean square value of the noise. Therefore, traditional amplitude thresholding methods for noise reduction are no longer applicable.

[0068] In this embodiment, in S1, a multi-scale wavelet transform is performed on the original noisy signal to obtain... The specific implementation method of the wavelet coefficient sequence at multiple scales includes the following steps:

[0069] S11: Acquire the raw, noisy vibration signal sequence of rotating machinery;

[0070] S12: Determine the maximum number of wavelet transform decomposition levels based on the sampling frequency and lowest frequency of the original noisy vibration signal sequence. ;

[0071] S13: Utilizing dual-tree complex wavelet transform to process the original noisy vibration signal sequence. Layer-by-layer recursive wavelet decomposition yields the scaling coefficients for each layer.

[0072] S14: Extract wavelet coefficient sequences from the scaling coefficients of each layer;

[0073] S15: Extract the original phase information from the scale coefficients of each layer;

[0074] Specifically, in S11, the original noisy vibration signal sequence of the rotating machinery. A vibration acceleration sensor mounted on the equipment housing or bearing housing is used to obtain a fixed sampling frequency. The acquired signal contains periodic weak impulse signals to be retained and background noise to be suppressed.

[0075] In S12, the impulse pulse signal has a steep rising edge and an exponentially decaying low-frequency tail in the time domain, and its energy is distributed over a relatively wide frequency band. After layer wavelet decomposition, the first The bandwidth corresponding to the layer approximation coefficient is approximately If the lowest frequency of the original noisy vibration signal sequence is ,but This allows us to calculate the maximum number of decomposition levels. ;

[0076] In S13, during the analysis of the original noisy vibration signal sequence of rotating machinery, the rotational speed fluctuations of the machinery cause a time shift in the impact pulse signal. If a conventional discrete wavelet transform is used, the time shift of the impact pulse signal will cause a large jump in the amplitude of the wavelet coefficients;

[0077] Therefore, this invention employs dual-tree complex wavelet transform, using two parallel real wavelet trees to form the real and imaginary parts of the complex coefficients, respectively. The filter bank of the two real wavelet trees is designed as an approximate Hilbert transform pair, thus possessing approximate translation invariance. When the impact pulse experiences a slight time shift due to rotational speed fluctuations, the amplitude of the wavelet coefficients at each scale remains stable, with only the phase changing.

[0078] The specific process of wavelet decomposition is as follows: starting from level 0, for the original noisy vibration signal sequence... The process involves filtering and downsampling layer by layer, with each layer outputting its scaling coefficients. These scaling coefficients are complex numbers and include both detail coefficients and approximation coefficients for that layer. The approximation coefficients are then passed to the next layer for further wavelet decomposition, yielding the wavelet coefficients and approximation coefficients for that layer, and so on, until... All layers have been decomposed.

[0079] In S14, no. Wavelet coefficients output by layer decomposition Actually, the department is The imaginary part is In this invention, the extracted amplitude sequence is: That is, wavelet coefficients The model;

[0080] In S15, the first Wavelet coefficients output by layer decomposition phase .

[0081] In this step, the original noisy vibration signal sequence is decomposed into time-domain waveforms of different frequency bands by performing multi-scale wavelet decomposition and amplitude extraction. The high-frequency and low-frequency signals in the original noisy vibration signal sequence are presented at different scales. Based on this, the vibration energy of the rotating machinery is transformed into waveform shape by extracting the amplitude sequence of wavelet coefficients.

[0082] In this embodiment, in S2, for each scale, a time-reversal sequence of wavelet coefficients is defined, and its specific implementation method includes the following steps:

[0083] S21: Obtain the wavelet coefficient sequences and their lengths at each scale output from step S1;

[0084] S22: For the wavelet coefficient sequence at each scale, with the last point of the wavelet coefficient sequence as the center of symmetry, the entire sequence is time-reversed to generate a time-reversed sequence;

[0085] Specifically, in S22, the first coefficient of the wavelet coefficient sequence is transformed into the last coefficient of the time inversion sequence, the second coefficient is transformed into the second-to-last coefficient of the time inversion sequence, and so on until each coefficient in the wavelet coefficient sequence has completed the time inversion process. The wavelet coefficient sequence and the time inversion sequence form a pair of positive and negative coefficients.

[0086] More specifically, in S21, the first The amplitude sequence of the scale output is ,in , The sampling point number ranges from 1 to... , This is the length of the wavelet coefficient sequence at this scale, since downsampling is performed during each wavelet decomposition. Divide the length of the original noisy vibration signal sequence by .

[0087] In S22, for the first The scale, and its time-reversal sequence are:

[0088] ;

[0089] in, The sampling point number ranges from 1 to... That is, the first coefficient of the wavelet coefficient sequence. After inversion, it becomes the last coefficient of the time-inversion sequence. The last coefficient of the wavelet coefficient sequence After inversion, it becomes the first coefficient of the time-inversion sequence. The rising edge in the wavelet coefficient sequence becomes the falling edge in the time inversion sequence.

[0090] It should be noted that the impact pulse signal in the vibration signal of rotating machinery has a clear temporal direction. Its waveform rises rapidly from a low amplitude at the baseline to a peak, and then slowly decays back down. Reversing this waveform along the time axis results in a waveform that rises slowly and then falls rapidly, completely different from the original waveform. Background noise, on the other hand, is randomly fluctuating and has no fixed time direction. After reversing it, the statistical characteristics of its local waveform are not significantly different from those before the reversal. This difference can be used to separate the impact pulse signal from the background noise.

[0091] In this embodiment, in S3, for each scale, the normalized similarity between the wavelet coefficient sequence and the time-reversed sequence is calculated. The specific implementation method includes the following steps:

[0092] S31: Set a local window;

[0093] S32: At each scale, starting from the beginning position of the wavelet coefficient sequence at that scale, slide the local window point by point along the direction of the wavelet coefficient sequence at that scale, and at each local window position, extract a segment of the wavelet coefficient sequence within the local window range;

[0094] S33: At each scale, starting from the beginning position of the time inversion sequence at that scale, slide a local window point by point along the direction of the time inversion sequence at that scale, and at each local window position, extract a time inversion sequence segment within the local window range;

[0095] S34: Calculate the normalized similarity between the wavelet coefficient sequence segment and the time-inverted sequence segment corresponding to each local window position;

[0096] In S31, the length of the local window is determined based on the duration of a single impact pulse signal in the rotating machinery from the start of oscillation to the decay to the base level, so that the local window can completely cover the process of an impact pulse signal from the start of oscillation to the decay. In S32, the wavelet coefficient sequence segment reflects the changing trend of the original vibration energy distribution within the local window, and the time inversion sequence segment reflects the changing trend of the original vibration energy within the local window after time reversal.

[0097] Furthermore, in S34, the normalized similarity between the wavelet coefficient sequence segment and the time-reversed sequence segment corresponding to each local window position is calculated. The specific calculation method is as follows:

[0098] Dividing the inner product of the wavelet coefficient sequence segment and the time-reversed sequence segment by the product of their respective moduli, we obtain a range of values. The standardized similarity value;

[0099] The closer the normalized similarity is to 1, the more consistent the changing trends of the wavelet coefficient sequence segment and its time-inverted sequence segment are within the local window; the closer the normalized similarity is to 0, the greater the difference in the changing trends of the wavelet coefficient sequence segment and its time-inverted sequence segment within the local window.

[0100] Specifically, the length of the local window This determines the time resolution and reliability of the causal asymmetry test. In the scenario of rotating machinery fault diagnosis, the duration of the impact pulse can be estimated by the resonant frequency and damping characteristics of the structure under test. For example, when the sampling frequency is 25600Hz and the structural resonant frequency is 3000Hz, the damped oscillation excited by a single impact pulse lasts approximately 5 to 10 cycles, corresponding to 40 to 80 sampling points. Therefore, the window length... The local window is set to 1.5 to 2 times the duration of the pulse, such as 64 or 128 points, to ensure that the local window can completely cover the rising edge, peak region, and decay tail of the impact pulse. This avoids the loss of causal information due to truncating the impact pulse shape because the local window is too short, and also avoids a decrease in the accuracy of causality judgment due to covering adjacent impact pulses or introducing too much noise because the local window is too long. The length of the local window... Once set, it remains fixed during the execution of step S3.

[0101] The length is The local window slides simultaneously from the starting positions of the wavelet coefficient sequence and the time inversion sequence, moving one sampling point at a time. The starting index of the local window is... Take in sequence In the first The window stays at a position, and the index range covered by the window is... to The wavelet coefficient sequence segment extracted within this index range is:

[0102] ;

[0103] The time-reversal sequence segments extracted from the same index range of the time-reversal sequence are:

[0104] ;

[0105] The standardized similarity is denoted as . Then the first Standardized similarity of the dwell positions of individual local windows The specific calculation formula is as follows:

[0106] ;

[0107] in, It is a very small positive number used to prevent division by zero errors when all time-reversal sequence segments and wavelet coefficient sequence segments within a local window are zero.

[0108] In the vibration signal of rotating machinery, for the rising edge position of the impact pulse, the wavelet coefficient sequence segment... Within this local window, it shows a rapid increasing trend, while the time-reversal sequence fragments... Within the same time window, due to time reversal, they exhibit a rapid decreasing trend, and their trends are opposite. The inner product calculation result is small, therefore the standardized similarity is... Less than 1. For pure noise locations, wavelet coefficient sequence segments and time-reversal sequence fragments All are random sequences with similar trends, resulting in a large inner product calculation; therefore, the standardized similarity is low. The similarity is close to 1. It should be noted that for a weak impulse pulse with a very small amplitude, as long as its rising edge shows a unidirectional increasing trend from low to high within the local window, its standardized similarity will deviate significantly from 1; while for a noise spike with a large amplitude, since the trend of change within the local window is non-directional, its standardized similarity is still close to 1. Compared with the limitation of the traditional amplitude threshold method failing at low signal-to-noise ratios, the processing effect of this method is better.

[0109] In this embodiment, in S4, the causal asymmetry index is calculated based on the standardized similarity. The specific implementation method includes:

[0110] S41: Obtain the standardized similarity corresponding to each local window position at each scale output in step S3;

[0111] S42: For each local window location at each scale, subtract the standardized similarity corresponding to that local window location from 1, and use the difference as the causal asymmetry index of that local window location;

[0112] The causal asymmetry indices at different local window locations at this scale constitute a causal asymmetry index sequence;

[0113] S43: For the positions at both ends of a sequence whose standardized similarity could not be calculated due to the inability of a local window to fully cover them at each scale, their causal asymmetry index is assigned a value of 0.

[0114] Specifically, the causal asymmetry index is denoted as... The calculation formula is as follows:

[0115] ;

[0116] The closer the causal asymmetry index is to 1, the greater the morphological difference between the wavelet coefficient sequence segment within the local window and its time-inverted sequence segment, meaning the higher the degree of morphological mismatch of the wavelet coefficient sequence segment within the local window under time inversion. Conversely, the closer the causal asymmetry index is to 0, the more similar the wavelet coefficient sequence segment within the local window is to its time-inverted sequence segment, meaning the wavelet coefficient sequence segment within the local window is basically the same before and after time inversion, indicating that the local location may be background noise.

[0117] It is important to emphasize that the calculation of the causal asymmetry index is based on the difference in the variation trend between a wavelet coefficient sequence segment and its time-reversed sequence segment, and is independent of the amplitude of the wavelet coefficient sequence segment. Therefore, even a weak impulse pulse with a very small amplitude will have a causal asymmetry index significantly greater than 0 if its rising edge exhibits a unidirectional increasing trend from low to high within the local window. Conversely, a noise spike with a large amplitude, because its variation trend within the local window is non-directional, will have a causal asymmetry index very close to 0, and thus will not be misidentified as an impulse pulse. This characteristic allows the causal asymmetry index to effectively identify weak impulse pulses with amplitudes far below the noise level from a strong noise background.

[0118] In this embodiment, in step S5, a shrinkage factor is constructed based on the causal asymmetry index, and the original wavelet coefficient sequence is weighted and shrunk using the shrinkage factor to obtain a denoised wavelet coefficient sequence. The specific implementation method includes:

[0119] S51: Obtain the causal asymmetric exponential sequence and wavelet coefficient sequence at each scale;

[0120] S52: For the index position of each wavelet coefficient at each scale, input the causal asymmetric exponential sequence of the index position into the preset shrinkage factor mapping function to obtain the shrinkage factor corresponding to the index position;

[0121] S53: For each index position at each scale, multiply the shrinkage factor corresponding to the index position by the wavelet coefficient at that position to obtain the denoised wavelet coefficients;

[0122] After processing all the wavelet coefficients at each index position at each scale, we obtain the denoised wavelet coefficient sequence at that scale.

[0123] Furthermore, in S52, before inputting the causal asymmetry exponential sequence into the preset contraction factor mapping function, it is necessary to construct the contraction factor mapping function. The specific construction process is as follows:

[0124] S521: Determine the type of basis function for the contraction factor mapping function, i.e., use the hyperbolic tangent function as the basis function for the mapping function;

[0125] S522: Introduce a gain constant, and use the product of the gain constant and the causal asymmetry index as the independent variable of the hyperbolic tangent function to obtain the complete form of the contraction factor mapping function;

[0126] S523: Determine the value of the gain constant to obtain the final shrinkage factor mapping function.

[0127] Specifically, the causal asymmetry index sequence is denoted as... ,in This is the scale index for wavelet transform. This refers to the index position, i.e., the sampling point number. This represents the length of the wavelet coefficient sequence at this scale. (Causal asymmetric exponential sequence) With wavelet coefficient sequence At the same scale and index position One-to-one correspondence, that is, each index position This also corresponds to a causal asymmetry index value. and wavelet coefficient values ;

[0128] No. Scale No. Shrinkage factor at each index position The causal asymmetry index at this location The shrinkage factor mapping function is generated by a contraction factor mapping function. This function is a monotonically increasing function with a domain and range of [0,1], and satisfies the condition that the output is 0 when the input is 0, and the output approaches 1 when the input approaches 1. This invention uses a combination of a hyperbolic tangent function and a gain constant as the contraction factor mapping function, and its calculation formula is as follows:

[0129] ;

[0130] in, The preset gain constant, and ; For the first Scale No. The causal asymmetry index of each index position; This is the contraction factor corresponding to the index position;

[0131] Gain constant Its function is to adjust the contraction factor mapping function by scaling the causal asymmetry exponent and then inputting the hyperbolic tangent function. The steepness of the transition within the value range, The larger the value, the steeper the transition of the contraction factor mapping function in the middle section of the causal asymmetry exponent; The smaller the value, the smoother the transition of the contraction factor mapping function in the middle section of the causal asymmetry index;

[0132] More specifically, the gain constant The value of needs to be set according to the noise level and impulse characteristics of the actual signal. If the gain constant... If the value is too large, pseudo-Gibberish oscillations can easily be introduced at the boundary between the impulse pulse and noise; if the gain constant is too large... If the value is too small, the wavelet coefficients in the noisy region cannot be sufficiently suppressed, leading to incomplete denoising. In this invention, the gain constant... A gain constant ranging from 3 to 8 can achieve good noise reduction results. As an example, the gain constant... The value is 5. When the gain constant... At that time, if =0.5, then When the causal asymmetry index reaches 0.5, the contraction factor is close to 1; if =0.1, then The contraction factor is in an intermediate state. This value ensures that wavelet coefficients are preserved only when the causal asymmetry exponent is large, and are attenuated when the causal asymmetry exponent is low. When the causal asymmetry exponent is in a gentle transition region, the wavelet coefficients change smoothly.

[0133] In this embodiment, in step S6, the denoised wavelet coefficient sequence is subjected to inverse wavelet transform to reconstruct the denoised original signal. The specific implementation method includes:

[0134] S61: Obtain the denoised wavelet coefficient sequence at each scale and the original phase information at each scale;

[0135] S62: Recombine the denoised wavelet coefficient sequence with the original phase information to restore it to the denoised complex wavelet coefficients;

[0136] S63: From the Starting from the first layer, the signal is reconstructed layer by layer using multi-scale wavelet inverse transform. That is, in each layer of reconstruction, the denoised complex wavelet coefficients of the current layer and the reconstruction approximation coefficients of the previous layer are used to recover the approximation coefficients of the current layer through a synthesis filter bank until the approximation coefficients of the 0th layer are recovered.

[0137] S64: Average the approximation coefficients of layer 0 to obtain the reconstructed denoised original signal.

[0138] Specifically, in S5, the wavelet coefficient sequences at each scale were subjected to shrinkage processing based on the causal asymmetry exponent to obtain the denoised amplitude sequences. ,in, This is the index position. In the dual-tree complex wavelet transform, the wavelet coefficients are complex numbers, containing both amplitude and phase information. In S5, only the amplitude of the wavelet coefficients is shrunk; their phase information remains unchanged.

[0139] Therefore, in this step, the denoised amplitude sequence With phase Recombining them yields the denoised complex wavelet coefficients. :

[0140] ;

[0141] Actually, the department is The imaginary part is Thus, the denoised complex wavelet coefficients retain both the noise amplitude suppression effect of S5 and the phase structure of the impulse pulse in the original signal, ensuring that the impulse pulse remains in the correct time position after reconstruction.

[0142] More specifically, the multi-scale wavelet inverse transform employs the dual-tree complex wavelet inverse transform corresponding to the forward decomposition in S1. During the forward decomposition, trees A and B respectively use their respective decomposition filter banks for filtering and downsampling; during the inverse transform, trees A and B respectively use their respective synthesis filter banks for upsampling and filtering. Let the synthesis low-pass filter of tree A be denoted as... The synthesized high-pass filter is The synthesized low-pass filter for tree B is: The synthesized high-pass filter is The synthesis filter and the decomposition filter used in S1 satisfy the complete reconstruction condition;

[0143] Multiscale wavelet inverse transform from the first The layer begins execution, the first The approximation coefficient of the layer is and Taking tree A as an example, the first The reconstruction formula for the layer approximation coefficients is:

[0144] ;

[0145] in, For the first The approximate coefficients obtained from the layer reconstruction For the noise reduction The real part of the layer detail factor;

[0146] This formula means: [The formula is incomplete and requires further context.] The approximation coefficients of the first layer are upsampled and filtered using a synthesis low-pass filter, while the denoised detail coefficients of the same layer are upsampled and filtered using a synthesis filter. The two are then added together to obtain the first layer. Approximation coefficients of the layer.

[0147] The reconstruction method of tree B is symmetrical to that of tree A, specifically as follows:

[0148] ;

[0149] from Begin by performing the above refactoring operations layer by layer upwards, until... To recover the approximate coefficients of the 0th layer and ;

[0150] Since the dual-tree complex wavelet transform uses two parallel real wavelet trees for processing, to eliminate the phase deviation between the two trees and ensure the accuracy of the reconstructed signal, the final denoised original signal is the average of the approximation coefficients of tree A and tree B reconstructed to the 0th level, that is:

[0151] ;

[0152] in, This is the final denoised vibration signal. , The length of the original noisy vibration signal sequence.

[0153] This step completes the signal reconstruction from wavelet coefficients in the transform domain to the time-domain vibration signal. By preserving the original phase information and only shrinking the amplitude, it is ensured that the impact pulse retains its original energy attenuation characteristics and is accurately located at the correct occurrence time after reconstruction.

[0154] Another aspect of the present invention discloses a signal denoising system based on wavelet transform, applicable to the aforementioned signal denoising method based on wavelet transform, comprising:

[0155] The wavelet transform module is used to perform multi-scale wavelet transform on the original noisy signal to obtain... A sequence of wavelet coefficients at various scales;

[0156] The time inversion module is used to define the time inversion sequence of the wavelet coefficient sequence for each scale;

[0157] The similarity calculation module is used to calculate the standardized similarity between the wavelet coefficient sequence and the time-inverted sequence;

[0158] The causal index calculation module is used to calculate the causal asymmetry index based on standardized similarity.

[0159] The shrinkage and denoising module is used to construct a shrinkage factor based on the causal asymmetry index, and to use the shrinkage factor to perform weighted shrinkage on the original wavelet coefficient sequence to obtain the denoised wavelet coefficient sequence.

[0160] The inverse transform reconstruction module is used to perform inverse wavelet transform on the denoised wavelet coefficient sequence to reconstruct the denoised original signal.

[0161] The various embodiments of the present invention have now been described in detail. To avoid obscuring the concept of the invention, some details known in the art have not been described. Those skilled in the art will fully understand how to implement the technical solutions of this invention based on the above description, and the scope of the invention is defined by the appended claims.

Claims

1. A signal denoising method based on wavelet transform, characterized in that, Includes the following steps: S1: Perform multi-scale wavelet transform on the original noisy signal to obtain... A sequence of wavelet coefficients at various scales; S2: For each scale, define the time-reversal sequence of the wavelet coefficient sequence; S3: Calculate the normalized similarity between the wavelet coefficient sequence and the time-inverted sequence; S4: Calculate the causal asymmetry index based on standardized similarity; S5: Construct a shrinkage factor based on the causal asymmetry index, and use the shrinkage factor to perform weighted shrinkage on the original wavelet coefficient sequence to obtain the denoised wavelet coefficient sequence. S6: Perform inverse wavelet transform on the denoised wavelet coefficient sequence to reconstruct the original denoised signal.

2. The signal denoising method based on wavelet transform according to claim 1, characterized in that, The specific implementation method of S1 includes the following steps: S11: Acquire the raw, noisy vibration signal sequence of rotating machinery; S12: Determine the maximum number of wavelet transform decomposition levels based on the sampling frequency and lowest frequency of the original noisy vibration signal sequence. ; S13: Utilizing dual-tree complex wavelet transform to process the original noisy vibration signal sequence. Layer-by-layer recursive wavelet decomposition yields the scaling coefficients for each layer. S14: Extract wavelet coefficient sequences from the scaling coefficients of each layer; S15: Extract the original phase information from the scale coefficients of each layer.

3. The signal denoising method based on wavelet transform according to claim 2, characterized in that, The specific implementation method of S2 includes the following steps: S21: Obtain the wavelet coefficient sequences and their lengths at each scale output from step S1; S22: For each scale of the wavelet coefficient sequence, with the last point of the wavelet coefficient sequence as the center of symmetry, the entire sequence is time-flipped to generate a time-inverted sequence.

4. The signal denoising method based on wavelet transform according to claim 3, characterized in that, The specific implementation method of S3 includes the following steps: S31: Set a local window; S32: At each scale, starting from the beginning position of the wavelet coefficient sequence at that scale, slide the local window point by point along the direction of the wavelet coefficient sequence at that scale, and at each local window position, extract a segment of the wavelet coefficient sequence within the local window range; S33: At each scale, starting from the beginning position of the time inversion sequence at that scale, slide a local window point by point along the direction of the time inversion sequence at that scale, and at each local window position, extract a time inversion sequence segment within the local window range; S34: Calculate the normalized similarity between the wavelet coefficient sequence fragment and the time-inverted sequence fragment corresponding to each local window position.

5. The signal denoising method based on wavelet transform according to claim 4, characterized in that, The specific implementation method of S4 includes the following steps: S41: Obtain the standardized similarity corresponding to each local window position at each scale output in step S3; S42: For each local window location at each scale, subtract the standardized similarity corresponding to that local window location from 1, and use the difference as the causal asymmetry index of that local window location; The causal asymmetry indices at different local window locations at this scale constitute a causal asymmetry index sequence; S43: For the positions at both ends of a sequence whose standardized similarity could not be calculated due to the inability of a local window to fully cover them at each scale, their causal asymmetry index is assigned a value of 0.

6. The signal denoising method based on wavelet transform according to claim 5, characterized in that, The specific implementation method of S5 includes the following steps: S51: Obtain the causal asymmetric exponential sequence and wavelet coefficient sequence at each scale; S52: For the index position of each wavelet coefficient at each scale, input the causal asymmetric exponential sequence of the index position into the preset shrinkage factor mapping function to obtain the shrinkage factor corresponding to the index position; S53: For each index position at each scale, multiply the shrinkage factor corresponding to the index position by the wavelet coefficient at that position to obtain the denoised wavelet coefficient.

7. The signal denoising method based on wavelet transform according to claim 6, characterized in that, The specific implementation method of S52 includes the following steps: S521: Determine the type of basis function for the contraction factor mapping function, i.e., use the hyperbolic tangent function as the basis function for the mapping function; S522: Introduce a gain constant, and use the product of the gain constant and the causal asymmetry index as the independent variable of the hyperbolic tangent function to obtain the complete form of the contraction factor mapping function; S523: Determine the value of the gain constant to obtain the final shrinkage factor mapping function.

8. The signal denoising method based on wavelet transform according to claim 7, characterized in that, The specific shrinkage factor mapping function is as follows: ； in, The preset gain constant, and ; For the first Scale No. The causal asymmetry index of each index position; This is the contraction factor corresponding to the index position.

9. A signal denoising method based on wavelet transform according to claim 8, characterized in that, The specific implementation method of S6 includes the following steps: S61: Obtain the denoised wavelet coefficient sequence at each scale and the original phase information at each scale; S62: Recombine the denoised wavelet coefficient sequence with the original phase information to restore it to the denoised complex wavelet coefficients; S63: From the Starting from the first layer, the signal is reconstructed layer by layer using multi-scale wavelet inverse transform. That is, in each layer of reconstruction, the denoised complex wavelet coefficients of the current layer and the reconstruction approximation coefficients of the previous layer are used to recover the approximation coefficients of the current layer through a synthesis filter bank until the approximation coefficients of the 0th layer are recovered. S64: Average the approximation coefficients of layer 0 to obtain the reconstructed denoised original signal.

10. A signal denoising system based on wavelet transform, used to implement the signal denoising method based on wavelet transform as described in any one of claims 1-9, characterized in that, include: The wavelet transform module is used to perform multi-scale wavelet transform on the original noisy signal to obtain... A sequence of wavelet coefficients at various scales; The time inversion module is used to define the time inversion sequence of the wavelet coefficient sequence for each scale; The similarity calculation module is used to calculate the standardized similarity between the wavelet coefficient sequence and the time-inverted sequence; The causal index calculation module is used to calculate the causal asymmetry index based on standardized similarity. The shrinkage and denoising module is used to construct a shrinkage factor based on the causal asymmetry index, and to use the shrinkage factor to perform weighted shrinkage on the original wavelet coefficient sequence to obtain the denoised wavelet coefficient sequence. The inverse transform reconstruction module is used to perform inverse wavelet transform on the denoised wavelet coefficient sequence to reconstruct the denoised original signal.

Images