In-band noise removal method based on difference time-frequency mode decomposition

Abstract

The invention discloses an in-band noise removal method based on differential time-frequency mode decomposition. The method comprises the following steps: acquiring a reference signal xr (t) and a mixed signal xm (t), and calculating normalized time-frequency amplitude spectrums of the reference signal xr (t) and the mixed signal xm (t); and obtaining an optimal difference time-frequency spectrum based on a large-margin convex optimization technology. And taking the normalized time-frequency amplitude spectrum of the reference signal and the normalized time-frequency amplitude spectrum of the mixed signal as samples, and obtaining an optimal difference time-frequency amplitude spectrum by adopting a large-margin convex optimization technology. Identifying time-frequency positions of different modes in the optimal difference time-frequency spectrum by using a threshold adaptive technology; and adaptively decomposing the time-frequency spectrum in the mixed signal into time-frequency spectrums of a sensitive component, a reference component and a noise component according to the time-frequency position. And transforming the corresponding time-frequency spectrum by using a time-frequency inverse transformation technology to obtain a corresponding sensitive component, a reference component and a noise component. Noise reduction analysis is carried out from the perspective of time frequency, the obtained optimal difference weight matrix has physical significance, and characteristics related to sensitive components are more highlighted.

Description

Technical Field

[0001] The present invention belongs to the technical field of signal processing, and relates to a signal noise reduction method, in particular to an in-band noise removal method based on differential time-frequency modal decomposition, which can be used to remove noise in the same frequency band of a signal. Background Art

[0002] Extracting sensitive components from complex multi-component signals has long been a research hotspot in signal processing. Numerous modal decomposition algorithms have been developed. Mainstream modal decomposition methods include empirical mode decomposition (EMD), variational mode decomposition (VMD), empirical wavelet transform (EVD), and singular value decomposition (SVD). Huang et al. (2014) proposed the EMD in "Empirical Mode Decomposition: A Method for Analyzing Neural Data," which decomposes signals into complete and nearly orthogonal eigenmode components. However, Lei et al. (2014) pointed out in "A Review on Empirical Mode Decomposition in Fault Diagnosis of Rotating Machinery" that EMD suffers from endpoint effects and modal aliasing. To mitigate these effects, Wu et al. (2014) proposed ensemble empirical mode decomposition (EEMD) in "ENSEMBLE EMPIRICAL MODE DECOMPOSITION: A NOISE-ASSISTED DATAANALYSIS METHOD." This method adds varying noise to the signal and then homogenizes the scales of the extreme points, resulting in a more accurate decomposition. To overcome the large reconstruction errors and poor decomposition completeness of EEMD, Yeh et al. proposed the complementary ensemble empirical mode decomposition (CEEMD) method in their paper "Complementary ensemble empirical mode decomposition: a novel noise enhanced data analysis method." This method adds pairs of positive and negative auxiliary white noise to the signal to be analyzed, which cancel each other out during ensemble averaging, effectively improving decomposition efficiency. However, the number of IMFs generated during EMD decomposition still varies, making it difficult to align the IMF components during the final ensemble averaging, or causing errors in the ensemble averaging.Torres et al. proposed the adaptive noise complete ensemble empirical mode decomposition (CEEMDAN) in their paper "Fault feature extraction of rotating machinery using a reweighted complete ensemble empirical mode decomposition with adaptive noise and demodulation analysis." During the decomposition process, white noise is added to the individual IMF components obtained by EMD decomposition. This results in a smaller noise residue in the reconstructed signal (compared to the EEMD result), reducing the number of screening operations. Furthermore, after the first-order intrinsic mode components of each signal group are decomposed by CEEMDAN, they are immediately ensemble-averaged. This avoids the problem of poor alignment of the final ensemble average due to differences in the IMF decomposition results of each group in CEEMD. It also prevents a poor decomposition of one order of IMF from passing its influence to the next order, affecting subsequent decompositions. Because these methods are essentially empirical and lack comprehensive theoretical support, Gilles et al. proposed the empirical wavelet transform (EWT) in their paper "Empirical Wavelet Transform." This method decomposes the different modes of a signal by designing appropriate wavelet filter banks, providing theoretical support. Subsequently, Dragomiretskiy et al. proposed variational mode decomposition (VMD) in the article “Variational Mode Decomposition”, which extracts different modes of the signal through a completely non-recursive variational mode decomposition model.

[0003] However, the modal decomposition method mentioned above is essentially a set of bandpass filters, which cannot filter out noise within the same frequency band as sensitive components. Therefore, it is necessary to develop a signal processing method that can remove in-band and out-of-band noise. Summary of the Invention

[0004] The purpose of the present invention is to provide an in-band noise removal method based on differential time-frequency modal decomposition, which addresses the problem that most existing signal processing methods cannot remove in-band noise.

[0005] The technical solution of the present invention:

[0006] A method for removing in-band noise based on differential time-frequency modal decomposition comprises the following steps:

[0007] Step S1: Collect the reference signal x of the rotating machinery r (t) and mixed signal x m (t);

[0008] Step S2: Calculate the normalized time-frequency amplitude spectrum of the reference signal and the mixed signal;

[0009] Taking the short-time Fourier transform (STFT) as an example, other time-frequency transform techniques can also be used;

[0010] Step S3: using the normalized time-frequency amplitude spectrum of the reference signal and the normalized time-frequency amplitude spectrum of the mixed signal as samples, a large margin convex optimization technique is used to obtain the optimal difference time-frequency amplitude spectrum;

[0011] Step S4: Using a threshold adaptive technique to identify the time-frequency positions of different modes in the optimal difference time-frequency amplitude spectrum; adaptively decomposing the time-frequency spectrum in the mixed signal into time-frequency spectra of a sensitive component, a reference component, and a noise component according to the time-frequency positions. The time-frequency spectrum at this time includes amplitude information and phase information.

[0012] Step S5: transform the corresponding time-frequency spectrum using the inverse time-frequency transform technique to obtain the corresponding sensitive component, reference component, and noise component;

[0013] In step S2, the short-time Fourier transform is used to normalize the time-frequency amplitude spectrum as follows:

[0014]

[0015] Where: x(τ) is the signal to be analyzed; ω(τ-t) is the window function; τ,t is time; ω is frequency; x ij is STFT x The original value of (t,ω); x min is STFT x The minimum value in (t,ω); x max is STFT x The maximum value in (t,ω); is the normalized value;

[0016] In step S3, the optimal difference time-frequency spectrum is obtained based on the large margin convex optimization technique; the specific calculation process is as follows:

[0017] Assume the dataset is Where P is the number of healthy data sets, Q is the number of fault data sets; i is the label of the sample, when l i = -1 is a healthy sample, when l i =1 is a fault sample; x i is the sample signal, i represents the i-th sample; based on the time-frequency transformation technology, the data set U is transformed into a normalized time-frequency amplitude spectrum data set F i is x iThe absolute value of the normalized time-frequency matrix with m rows and k columns after STFT transformation is used to define the weighted normalized time-frequency amplitude spectrum WNTFM as follows:

[0018] WNTFM=W⊙F

[0019] Where: W∈R m×k is the weight matrix, F∈R m×k is the sample signal x i The normalized time-frequency amplitude spectrum of , ⊙ is the operator of time-frequency matrix operation; the specific operation is as follows:

[0020] (A⊙B) j，k ∈R m×n =A j，k B j，k , j=1,2,...,m and k=1,2,...,n

[0021] Where: A j，k and B j，k Represent the elements in the jth row and kth column of matrix A and matrix B respectively;

[0022] The optimal weight time-frequency matrix W is obtained by establishing the optimization model shown below: * ;

[0023]

[0024] stW=Y

[0025] Where κ is a hyperparameter; C is a regularization parameter; the trace of the tr() matrix; C is a constant; b is a parameter that offsets the hyperplane; y i is the label of the i-th sample; tr(W T W) represents the margin between the time-frequency matrix and the hyperplane; ||Y|| * It is a measure of the rank of W. Minimizing it can capture or reflect some fault characteristics in the time-frequency domain covered by several large singular values; It is based on the hyperplane separation of normal and fault time-frequency matrices;

[0026] The augmented Lagrange equation of the above formula is:

[0027]

[0028] Where, {·} + is a set of indicator functions of integers; σ i (Y) is the singular value of Y; ρ>0 is a hyperparameter; Ω is the augmented Lagrangian parameter;

[0029] In the above formula, the optimal weighted time-frequency matrix W can be obtained by the multiplication alternation algorithm:* , continuously update the relevant parameters according to the following two formulas;

[0030]

[0031]

[0032] Where W (k) ,b (k) , Y (k) The optimal value of the kth iteration;

[0033] In step S4, the specific process of using the threshold adaptive technology to identify the time-frequency positions of different modes in the optimal difference time-frequency spectrum is as follows:

[0034] Step 1: Find the optimal weight time-frequency matrix W * ∈R m×k The minimum and maximum values ​​of each column are as follows:

[0035]

[0036] Step 2: Find L min Positions greater than 0 are considered to be the locations of sensitive components Γ sc , as shown below:

[0037] Γ sc ={imin_value i >0 for i∈{1,2,...,m}}

[0038] Step 3: Find L max Positions less than 0 are considered to be the positions of the basic components Γ rc , as shown below:

[0039] Γ rc ={imax_value i <0 for i∈{1,2,...,m}}

[0040] Step 4: The remaining time-frequency positions are regarded as the positions of the noise components Γ no , as shown below:

[0041]

[0042] Step 5: Find the corresponding time-frequency positions of the sensitive signal, reference signal and noise from the time-frequency matrix of the mixed signal to be analyzed.

[0043] In step S5, the corresponding time-frequency spectrum is transformed using the inverse time-frequency transform technology to obtain the corresponding sensitive component, reference component and noise component. The specific process is as follows: according to the time-frequency position located in step S4, locate the TFS in the time-frequency spectrum of the mixed signal to be analyzed x sc , TFS rc ,TFS no , and then reconstruct the decomposed modes using the inverse time-frequency transformation technique; the calculation formula is as follows:

[0044] x sc =iSTFT(TFS sc ),x rc =iSTFT(TFS rc ),x no =iSTFT(TFS no )

[0045] Among them, x sc ,x rc ,x no The decomposed sensitive modes, fundamental modes and noise modes respectively.

[0046] Beneficial effects of the present invention:

[0047] (1) The present invention relies on convex optimization techniques and then analyzes each time-frequency position, thereby eliminating in-band noise.

[0048] (2) The present invention performs noise reduction analysis from the perspective of time and frequency, and the obtained optimal difference weight matrix has physical meaning and further highlights the characteristics related to the sensitive components. BRIEF DESCRIPTION OF THE DRAWINGS

[0049] Figure 1 is a flow chart of the diagnostic method of the present invention;

[0050] Figure 2 : is a time domain diagram of the reference signal, sensitive signal and mixed signal of the present invention; wherein (a) is the reference signal without noise in this case, (b) is the sensitive signal without noise in this case, (c) is the reference signal with noise in this case, and (d) is the mixed signal in this case;

[0051] Figure 3 : are frequency domain diagrams of the reference signal, sensitive signal, and mixed signal of the present invention; wherein (a) is the spectrum of the reference signal and sensitive signal without noise, (b) is the spectrum of the reference signal with noise, and (c) is the spectrum of the mixed signal;

[0052] Figure 4 is the optimal weight difference time-frequency spectrum of the present invention;

[0053] Figure 5is the adaptive threshold positioning map of the present invention; wherein, (a) is the adaptive positioning sensitive component position map, and (b) is the adaptive positioning reference component position map;

[0054] Figure 6 : This is a time domain diagram of the reshaped reference signal, sensitive signal and mixed signal of the present invention; wherein (a) is the sensitive component, (b) is the reference component, and (c) is the noise component;

[0055] Figure 7 3 is a frequency domain diagram of the reshaped reference signal and sensitive signal according to the present invention; wherein (a) is the spectrum of the sensitive component, and (b) is the spectrum of the reference component. DETAILED DESCRIPTION

[0056] The specific implementation of the present invention is further described below in conjunction with the accompanying drawings and technical solutions.

[0057] The in-band noise removal method based on differential time-frequency modal decomposition is as follows: Figure 1 As shown, the steps are as follows:

[0058] Step S1: Collect reference signal x r (t) and mixed signal x m (t);

[0059] The specific simulation is as follows: the reference signal, sensitive signal and mixed signal are used as the reference signal x r (t) and mixed signal x m (t):

[0060]

[0061] x mix =x sc +x rc +x no

[0062] Where x sc , x rc , x no and x mix They are sensitive signal, reference signal, noise component and mixed signal respectively. N S and N R represents the number of harmonics, f i S and f i R Represents amplitude and frequency. The sampling frequency and sampling duration are set to 5120Hz and 4s respectively. rc The specific parameter values ​​are shown in Table 1, and the corresponding time domain is as follows Figure 2 As shown in (a), xsc The specific parameter values ​​are shown in Table 2, and the corresponding time domain is as follows Figure 2 (b) is shown. no represents Gaussian white noise with a standard deviation of σ = 2. rc like Figure 2 As shown in (c), the mixed signal x mix like Figure 2 (d) As shown. Furthermore, the frequency domain will give us a more comprehensive understanding of the signal in the simulation case. Figure 3 (a), Figure 3 (b) and Figure 3 (c) are the spectra of the pure reference component and the sensitive component (the red spectrum is the spectrum of SC, and the black spectrum is the spectrum of RC), the spectrum of the reference component after adding noise, and the spectrum of the mixed signal. Our goal is to get the mixed signal x mix Separate x sc , x rc and x no .

[0063] Table 1 Pure reference signal x rc The values ​​of each parameter in

[0064]

[0065] Table 2 Pure sensitive signal x sc The values ​​of each parameter in

[0066]

[0067]

[0068] Step S2: Calculate x r (t) and x m The normalized time-frequency amplitude spectrum of (t). The present invention uses short-time Fourier transform (STFT) as an example, and other time-frequency transform technologies can also be used;

[0069] The normalized time-frequency amplitude spectrum taking STFT as an example is as follows:

[0070]

[0071] Where: x(τ) is the signal to be analyzed; ω(τ-t) is the window function; τ,t is time, and ω is frequency. ij is STFT x The original value of (t,ω), x min is STFT x The minimum value in (t,ω); x max is STFT x The maximum value in (t,ω); is the normalized value.

[0072] Step S3: Obtain the optimal difference time-frequency spectrum based on the large margin convex optimization technique. Using the normalized time-frequency amplitude spectrum of the reference signal and the normalized time-frequency amplitude spectrum of the mixed signal as samples, the large margin convex optimization technique is used to obtain the optimal difference time-frequency amplitude spectrum. The specific steps are as follows:

[0073] Take 50 mixed signals and 50 reference signals after adding noise to form a data set: Where P is the number of reference datasets, Q is the number of mixed datasets, and l i is the label of the sample, when l i = -1 is the reference sample, when l i =1 is a mixed sample. Based on time-frequency transformation technology, is the normalized time-frequency amplitude spectrum dataset. i is the signal x i The absolute value of the normalized time-frequency matrix with m rows and k columns after STFT transformation. The weighted normalized time-frequency spectrum WNTFM is defined as:

[0074] WNTFM=W⊙F

[0075] Where: W∈R m×k is the weight matrix, F∈R m×k is the normalized time-frequency amplitude spectrum of signal x. ⊙ is the operator for time-frequency matrix operations. The specific operations are as follows:

[0076] (A⊙B) j,k ∈R m×n =A j,k B j,k ,j=1,2,...,m and k=1,2,...,n

[0077] Where: A j,k and B j,k Represent the elements in the jth row and kth column of matrix A and matrix B respectively;

[0078] The optimal weight time-frequency matrix W is obtained by establishing the optimization model shown below: * .

[0079]

[0080] stW=Y

[0081] Where κ is a hyperparameter; C is a regularization parameter; the trace of the tr() matrix; C is a constant; b is a parameter that offsets the hyperplane, and y i is the label of the i-th sample; tr(W T W) represents the margin between the time-frequency matrix and the hyperplane. ||Y||* It is a measure of the rank of W, minimizing it so that it can capture or reflect some sensitive features in the time-frequency domain covered by several larger singular values. The purpose is to separate the reference and mixed time-frequency matrices based on a hyperplane.

[0082] The augmented Lagrange equation of the above formula is:

[0083]

[0084] Where, {·} + is a set of indicator functions of integers; σ i (Y) are the singular values ​​of Y; ρ>0 is a hyperparameter, and Ω is the augmented Lagrangian parameter.

[0085] In the above formula, the optimal weight time-frequency matrix W can be obtained by the multiplication alternation algorithm * ,like Figure 4 As shown, the relevant parameters are continuously updated according to the following two formulas.

[0086]

[0087] Step S4: Use the threshold adaptive technology to identify the time-frequency positions of different modes in the optimal difference time-frequency spectrum; according to the time-frequency position, the time-frequency spectrum in the mixed signal is adaptively decomposed into the time-frequency spectrum of the sensitive component, the reference component and the noise component. The time-frequency spectrum at this time includes amplitude information and phase information. Specific calculation steps

[0088] Step 1: Find the optimal difference time-frequency matrix W * ∈R m×k The minimum and maximum values ​​of each column are shown below.

[0089]

[0090] Step 2: Find L min Positions greater than 0 are considered to be the locations of sensitive components Γ sc , as shown below, Figure 5 (a).

[0091] Γ sc ={imin_value i >0 for i∈{1,2,...,m}}

[0092] Step 3: Find L max Positions less than 0 are considered to be the positions of the basic components Γ rc , as shown below, Figure 5 (b).

[0093] Γrc ={i|max_value i <0 for i∈{1,2,...,m}}

[0094] Step 4: The remaining time-frequency positions are regarded as the positions of the noise components Γ no , as shown below.

[0095]

[0096] Step 5: Find the corresponding time-frequency positions of the sensitive signal, reference signal and noise from the time-frequency matrix of the mixed signal to be analyzed.

[0097] Step S5: Use the inverse time-frequency transform technology to transform the corresponding time-frequency spectrum to obtain the corresponding sensitive component, reference component and noise component; according to the time-frequency position located in step S4, locate the TFS in the time-frequency spectrum of the mixed signal to be analyzed x sc , TFS rc ,TFS no , and then reconstruct the decomposed modes using the inverse time-frequency transformation technique. The calculation formula is as follows:

[0098] x sc =iSTFT(TFS sc ),x rc =iSTFT(TFS rc ),x no =iSTFT(TFS no )

[0099] Here x sc ,x rc ,x no The sensitive modes, basic modes and noise modes are decomposed respectively, such as Figure 6 As shown, Figure 6 (a) is the sensitive component, Figure 6 (b) is the reference component, Figure 6 (c) is the noise component. The corresponding spectrum is as follows Figure 7 As shown, Figure 7 (a) is the spectrum of the sensitive component, Figure 7 (b) is the spectrum of the reference component. As can be seen from the above figures, the method of the present invention can successfully remove the in-band noise.

[0100] The above descriptions are merely preferred embodiments of the present invention, and all equivalent changes and modifications made within the scope of the patent application of the present invention should fall within the scope of the present invention.

Claims

1. A method for removing in-band noise based on differential time-frequency modal decomposition, characterized in that: The steps include: Step S1: Collect the reference signal x of the rotating machinery r (t) and the mixed signal x m (t); Step S2: Calculate the normalized time-frequency amplitude spectrum of the reference signal and the mixed signal; Step S3: using the normalized time-frequency amplitude spectrum of the reference signal and the normalized time-frequency amplitude spectrum of the mixed signal as samples, and adopting a large-residue convex optimization technique to obtain an optimal difference time-frequency amplitude spectrum; Step S4: using a threshold adaptive technique to identify the time-frequency positions of different modes in the optimal difference time-frequency amplitude spectrum; according to the time-frequency position, the time-frequency spectrum in the mixed signal is adaptively decomposed into the time-frequency spectrum of the sensitive component, the reference component and the noise component, and the time-frequency spectrum at this time includes amplitude information and phase information; Step S5: transform the corresponding time-frequency spectrum using the inverse time-frequency transform technique to obtain the corresponding sensitive component, reference component and noise component.

2. The method for removing in-band noise based on differential time-frequency modal decomposition according to claim 1, characterized in that: In step S2, the short-time Fourier transform is used to normalize the time-frequency amplitude spectrum as follows: Where: x(τ) is the signal to be analyzed; ω(τ-t) is the window function; τ,t is time; ω is frequency; x ij It is STFT x The original value of (t,ω); x min It is STFT x The minimum value in (t,ω); x max It is STFT x The maximum value in (t,ω); is the normalized time-frequency amplitude spectrum.

3. The in-band noise removal method based on differential time-frequency modal decomposition according to claim 1 is characterized in that: In step S3, the optimal difference time-frequency spectrum is obtained based on the large margin convex optimization technique; the specific calculation process is as follows: Assume the data set is Where P is the number of healthy data sets, Q is the number of fault data sets; i is the label of the sample, when l i = -1 is a healthy sample, when l i =1 is a fault sample; x i is the sample signal, i represents the i-th sample; based on the time-frequency transformation technology, the data set U is transformed into a normalized time-frequency amplitude spectrum data set F i is x i The absolute value of the normalized time-frequency matrix with m rows and k columns after STFT transformation is used to define the weighted normalized time-frequency amplitude spectrum WNTFM as follows: WNTFM=W⊙F Where: W∈R m×k is the weight matrix, F∈R m×k is the sample signal x i The normalized time-frequency amplitude spectrum of ⊙ is the operator of time-frequency matrix operation; the specific operation is as follows: (A⊙B) j，k ∈R m×n =A j，k B j，k ,j=1,2,...,m and k=1,2,...,n,2,...,n Where: A j,k and B j,k Represent the elements in the jth row and kth column of matrix A and matrix B respectively; The optimal weighted time-frequency matrix W is obtained by establishing the following optimization model: * ; stW=Y Where κ is a hyperparameter; C is a regularization parameter; tr() is the trace of the matrix; C is a constant; b is a parameter that offsets the hyperplane; y i is the label of the i-th sample; tr(W T W) represents the margin between the time-frequency matrix and the hyperplane; ||Y|| * It is a measure of the rank of W, which can be minimized to capture or reflect some fault characteristics in the time-frequency domain covered by several larger singular values; It is based on the hyperplane separation of normal and faulty time-frequency matrices; The augmented Lagrangian equation of the above formula is: In the formula, {·} + is a set of indicator functions of integers; σ i (Y) is the singular value of Y; ρ>0 is a hyperparameter; Ω is the augmented Lagrangian parameter; In the above formula, the optimal weighted time-frequency matrix W is obtained by the multiplication alternating algorithm: * , continuously update the relevant parameters according to the following two formulas; Where W (k) ,b (k) , Y (k) The optimal value for the kth iteration.

4. The in-band noise removal method based on differential time-frequency modal decomposition according to claim 1 is characterized in that: In step S4, the specific process of using the threshold adaptive technology to identify the time-frequency positions of different modes in the optimal difference time-frequency spectrum is as follows: Step 1: Find the optimal weight time-frequency matrix W * ∈R m×k The minimum and maximum values ​​of each column are as follows: Step 2: Find L min Positions greater than 0 are considered to be the locations of sensitive components Γ sc , as shown below: Γ sc ={i|min_value i >0 for i∈{1,2,...,m}} Step 3: Find L max Positions less than 0 are considered to be the positions of the basic components Γ rc , as shown below: Γ rc ={i|max_value i <0 for i∈{1,2,...,m}} Step 4: The remaining time-frequency positions are regarded as the positions where the noise components are located Γ no , as shown below: Step 5: Find the corresponding time-frequency positions of the sensitive signal, reference signal and noise from the time-frequency matrix of the mixed signal to be analyzed.

5. The method for removing in-band noise based on differential time-frequency modal decomposition according to claim 1, characterized in that: In step S5, the corresponding time-frequency spectrum is transformed by using the inverse time-frequency transform technology to obtain the corresponding sensitive component, reference component and noise component. The specific process is as follows: according to the time-frequency position located in step S4, the TFS in the time-frequency spectrum of the mixed signal x to be analyzed is located. sc 、TFS rc 、TFS no , and then reconstruct the decomposed modes using the inverse time-frequency change technique; the calculation formula is as follows: x sc =iSTFT(TFS sc ),x rc =iSTFT(TFS rc ),x no =iSTFT(TFS no ) Among them, x sc ,x rc ,x no The decomposed sensitive modes, fundamental modes and noise modes respectively.

Images