CN106645947A - Time-frequency analysis method based on nonlinear mode decomposition and adaptive optimal kernel - Google Patents

Abstract

The invention discloses a time-frequency analysis method based on nonlinear mode decomposition and adaptive optimal kernel; the time-frequency analysis method combines the advantages in nonlinear mode decomposition analysis method and those in adaptive optimal kernel and comprises: decomposing a multi-component nonstable signal into a group of physically significant nonlinear mode components by using nonlinear mode composition algorithm, wherein the requirement on the anti-noise performance of subsequent analysis method is lowered since the algorithm has high noise robustness; then, enabling a kernel function to change adaptively with signal changes by using the time-frequency analysis method of adaptive optimal kernel, so that cross terms are effectively inhibited and time-frequency concentrating capacity is improved. The novel analysis method inherits the advantages in both nonlinear mode composition method and adaptive optimal kernel analysis method, and has excellent performances.

Description

A kind of Time-Frequency Analysis Method decomposed based on nonlinear model with adaptive optimal kernel

Technical field

The present invention relates to the time frequency analysis field of non-stationary signal, especially a kind of to be decomposed based on nonlinear model and adaptive Answer the Time-Frequency Analysis Method of optimum core.

Background technology

Effect of one good Time-Frequency Analysis Method in the estimation analysis of non-stationary signal is self-evident.Signal itself has There are many attributes, and for Signal estimation, frequency domain characteristic is critically important attribute.

There are many signal transactings and the method estimated, most classical certainly Fourier transform at present, it can reflect letter Number frequency characteristic, the Fourier transform frequency domain resolution capability good for stationary signal possesses can obtain signal after conversion Frequency spectrum, but Fourier transform lost temporal information, i.e., the time that each frequency spectrum occurs is not aware that, so Fourier becomes Change and be not suitable for processing non-stationary signal.For the non-stationary signal in actual environment, conventional linear time-frequency representation method has short When Fourier transformation, wavelet transformation, S-transformation etc., but their time-frequency precision has much room for improvement；And Wigner-Ville distribution, Although the bilinearity time-frequency distributions such as Cohen class time-frequency distributions were inevitably being calculated with higher time-frequency precision Cross term is introduced in journey.

Sophisticated signal can be decomposed into one and be by empirical mode decomposition (Empirical Mode Decomposition, EMD) Row intrinsic mode functions (Intrinsic Mode Function, IMF).Although EMD has in the cross term for suppressing multicomponent data processing Original performance, but its shortcoming it is also obvious that noiseproof feature is poor, and meeting is because the presence of discontinuous signal or noise And bring modal overlap phenomenon.In order to overcome the shortcoming of EMD, a kind of analysis method aided in based on noise is occurred in that, be called collection Close empirical mode decomposition (Ensemble Empirical Mode Decomposition, EEMD).Although EEMD is in noiseproof feature On improve than EMD, but still reach to less than desired value.

Time-frequency representation (Time-Frequency Representation, TFR) is that a kind of highly effective instantaneous frequency is estimated Meter method, especially for multicomponent data processing.General TFR methods only have the window function or kernel function of fixation, therefore they are only It is relatively good to the treatment effect of a certain class signal.Adaptive optimal kernel time-frequency representation (AOK TFR) is a kind of wink based on signal When frequency estimating methods, using the radially Gaussian kernel function that can change with signal intensity, so, when AOK has preferable Frequency focusing performance and suppressing crossterms ability.

EMD+AOK can reduce the impact of cross term when multicomponent data processing is analyzed, but excessively sensitive to noise. EEMD+AOK can to a certain extent improve noiseproof feature, but also far from enough.

Nonlinear model decomposes (Nonlinear Mode Decomposition, NMD) What is arranged has the nonlinear model component of physical significance, with very strong noise robustness.It is that time frequency analysis, data are substituted into inspection Test and a kind of New Algorithm that the fusion of various methods gets up such as harmonic wave differentiates.However, decomposing the nonlinear model for obtaining by NMD Formula component is not the simple component as IMF, need to be changed with the change of signal with the kernel function energy self adaptation of AOK methods Feature is being optimized.

Based on the new time-frequency representation method of NMD and AOK, NMD is not only make use of to the effective decomposability of multicomponent data processing Can, while also inherit the outstanding time-frequency focusing performance of AOK and the ability of effective suppressing crossterms, solve EMD+AOK and The noiseproof feature difference of EEMD+AOK and the still problem containing a certain amount of cross term.

The content of the invention

Goal of the invention：To solve above-mentioned technical problem, the present invention provides a kind of decomposition based on nonlinear model and self adaptation The Time-Frequency Analysis Method of optimum core, (Nonlinear Mode Decomposition, NMD) is decomposed treating using nonlinear model Process signal is decomposed into a series of nonlinear model components, suppresses the interference of noise；Then each component of signal for obtaining is passed through Adaptive optimal kernel (AOK) analysis method processes suppressing crossterms, because the analysis method (AOK) of adaptive optimal kernel can be automatically The change of trace analysis signal, kernel function can the adaptive change according to the change of signal, the core letter that different signal is produced Several always optimums, so its time frequency analysis Performance comparision is good.

Technical scheme：To realize above-mentioned technique effect, the technical scheme that the present invention is provided is：

A kind of Time-Frequency Analysis Method based on nonlinear model decomposition and adaptive optimal kernel includes step：

(1) it is s (t) to define pending signal, and the sample frequency of s (t) is f_s, data length be N；Using nonlinear model Pending signal s (t) is decomposed into one group of nonlinear model component by decomposition method, i.e.,：

Wherein, c_iT () is i-th nonlinear model component of s (t), n (t) represents noise；

(2) each nonlinear model component that step (1) is decomposited is analyzed using adaptive optimal nuclear analytical method, Including step：

(2-1) radially Gaussian kernel function is built

In formula,Represent for controlling radial direction Gaussian function in radial angleThe spread function in direction；

(2-2) to obtain as the optimum core of signal adaptive change is as target problem, building optimization problem model is：

Arrange optimization problem model constraints be：

Wherein,For i-th nonlinear model component c_iT the ambiguity function of () under polar coordinates, β is optimum core Volume；Expression formula in rectangular coordinate system is：

Wherein, A_i(t, θ, τ) isMap amount in rectangular coordinate system, s^*(t),w^*U () is s (t), c_i (t), the conjugate complex number form of w (u)；Window function w (u) be with the center of t, width for 2T symmetrical rhombus window function, and | u | ＞ During T, w (u)=0, variable τ and θ are the parameters of general fuzzy field { τ, u }, | τ | ＜ 2T；

(2-3) according to A_i(t, θ, τ) solving-optimizing problem model, obtains nonlinear model component c_iThe optimum kernel function of (t) Φ_(i)opt(t, θ, τ)；

(2-4) nonlinear model component c is calculated_iT the time-frequency representation of the adaptive optimal kernel of () is：

In formula,For nonlinear model component c_iThe energy value of (t) in t；

(3) time-frequency representation of the adaptive optimal kernel of all nonlinear model components obtained according to step (2), during calculating The result of frequency analysis is：

Further, pending signal s (t) is decomposed into one group using nonlinear model decomposition method in the step (1) The method of nonlinear model component includes step：

(1-1) wavelet transformation expression formula W of signal s (t) is calculated_s(ω, t)：

Wherein,It is the Fourier transform of s (t), i.e.,s⁺T () is the positive frequency of s (t) signals Part,It is wavelet function, withFourier transform pair each other, and meet conditionψ^*(t),It is respectively ψ (t),Conjugate complex number；ω_ψSmall echo crest frequency is represented, ω_ψ=1, f₀For resolution parameter, for weighing in conversion process Middle time and the resolution ratio of frequency；

(1-2) judge that whether the calculated wavelet transformation of step (1-1) is the optimal time-frequency representation of s (t)；If judging knot Fruit is no, then calculate adding window Fourier transformation expression formula G of signal s (t)_s(ω, t)：

Wherein, g (t) is the window function of adding window Fourier transformation,For the Fourier transformation of g (t), following bar is met Part：

And

(1-3) all of ridge curve of the optimal time-frequency representation of signal s (t) is found out, and harmonic component is reconstructed with ridge method, its In, h order harmonic components are：

x^(h)(t)=A^(h)(t)cosφ^(h)(t), h ∈ [1,2 ..., N]

In formula, A^(h)(t)、φ^(h)T () is respectively the amplitude and phase place of h order harmonic components；v^(h)T () is h subharmonic The frequency of component, y^(h)(t)≡φ′^(h)(t)；φ′^(h)T () is φ^(h)The derivative of (t) to time t；N represents harmonic component most High reps, i.e. data length；

(1-4) true and false for substituting harmonic component of the method for inspection to extracting using noise immunity differentiates, filters out institute The real harmonic component having,

(1-5) all of real harmonic component is added, obtains a nonlinear model component c₁(t)；

(1-6) nonlinear model is cut from original signal s (t) and decomposes the c for obtaining₁(t), and residual components are repeated to walk Rapid 1-1 obtains all of each nonlinear model component c to step 1-6_i(t)；

Finally, former echo signal s (t) can be expressed as：

In formula, n (t) represents noise.

Further, the method that all of ridge curve of the optimal time-frequency representation of signal s (t) is found out in the step (1-3) For：

In moment t_i, h maximum point is found out, and the h maximum point found out is connected, form moment t_iRidge curve For：

I=1 in above formula, 2 ..., N, N are data length；H_s(ω, t) be signal s (t) optimal time-frequency representation, as W_s (ω, t) or G_s(ω, t)；H_s(ω, t) in the ridge point line at all moment found out, that is, constitute N vallate curves.

Further, h order harmonic components x are reconstructed with ridge method in the step (1-3)^(h)(t)=A^(h)(t)cosφ^(h) T () comprises the following steps：

If the optimal time-frequency representation of pending signal s (t) is wavelet transformation,

If the optimal time-frequency representation of pending signal s (t) is windowed FFT,

In formula,WithIt is by improved discretization factor of influence obtained from parabola interpolation.

Further, the method for the real harmonic component of screening includes step in the step (1-4)：

(1-4-1) identification statistic D (α is built_A, α_ν), for weigh each harmonic component amplitude and frequency it is orderly Degree；D(α_A, α_v) expression formula be：

In formula,WithA is represented respectively^(h)(t) and v^(h)The spectrum entropy of (t)；α_A, α_vFor weight coefficient；

(1-4-2) N is created for original signal s (t)_sIndividual Fourier transform alternate data；Wherein, j-th Fourier transform is replaced The expression formula of codes or data is：

(1-4-3) the optimal time-frequency representation formula of each Fourier transform alternate data is calculated, and respectively from wherein carrying Each harmonic component is taken out, calculating the corresponding identification statistic of each Fourier transform alternate data is：

Define significance indexIn formulaTo meet D_s＞ D₀ Alternate data number；D₀For the degree of order of the h subharmonic of original signal；

It is p to arrange level of signifiance index, and when there is x alternate data D is met_s＞ D₀And x >=p × N_sWhen, judge corresponding Harmonic component is not noise；Otherwise, it is determined that corresponding harmonic component is noise；

(1-4-4) the comprehensive measurement value of the degree of correlation between harmonic wave is calculated

Wherein,

In formula, w_A, w_φ, w_vRepresentWeights, ρ^(h)It is that amplitude and phase equalization distribute equal weights ρ^(h)≡ρ^(h)(1,1,0)；

(1-4-5) threshold value of comprehensive measurement value is defined

(1-4-6) when a harmonic component meets ρ^(h)≥ρ_minAnd during significance index p >=95%, judge this harmonic wave Component, by inspection, is real harmonic component.

Further, judge whether an alternate data meets D in the step (1-4-3)_s＞ D₀Method be：

Choose the parameter (α of any three groups of different values_A, α_v), three groups of corresponding identifications of parameter are calculated respectively counts value, if There is any one identification statistics value and be more than D₀, then judge that the alternate data meets D_s＞ D₀。

Beneficial effect：Compared with prior art, the present invention has the advantage that：

NMD methods have good inhibition not only for white Gaussian noise, while also can effectively suppress other classes The noise signal of type, with very strong noise robustness and extensive adaptability；However, by NMD decompose obtain it is non-linear Mode component is not the simple component as IMF, need to be changed with the change of signal with the kernel function energy self adaptation of AOK methods The characteristics of being optimized, AOK methods can effectively, adaptively track the change of non-stationary signal.The present invention is based on NMD With the new time-frequency representation method of AOK, NMD is not only make use of to the effective decomposability of multicomponent data processing, while also inheriting The outstanding time-frequency focusing performance and the ability of effective suppressing crossterms of AOK, solves the noiseproof feature of EMD+AOK and EEMD+AOK Difference and the still problem containing a certain amount of cross term.

Description of the drawings

Fig. 1 is the flow chart of the present invention.

Specific embodiment

Below by taking multi -components non-stationary emulation signal as an example, embodiments of the present invention and superior are described in detail Property.

Assume that s (t) is a multicomponent data processing containing white Gaussian noise：

S (t)=s_d(t)+n(t)

s_d(t)=cos (20 π t)+sin (200 π t)+sin (400 π t)+sin (100 π (t-0.5)²)

In formula, n (t) represents white Gaussian noise；s_dT () is preferable multicomponent data processing.Setting sample frequency is 1kHz, is adopted Sample time 1s, data length is 1000.Window length 2T=128 is set, and kernel function volumetric constraint is β=5.

The flow process of the present invention is as shown in figure 1, comprise the following steps：

Step A：Prepare pending signal s (t), its sample frequency is f_s, data length be N；

Step B：NMD analyses are carried out to signal s (t)；

Step B-1：Calculate wavelet transformation (Wavelet Transform, the WT) W of signal s (t)_s(ω, t),

Wherein,It is the Fourier transform of s (t), i.e.,s⁺T () is the positive frequency of s (t) signals Part,ψ (t) is wavelet function, withFourier transform pair each other, and meet conditionψ^*(t),It is respectively ψ (t),Conjugate complex number.Represent small echo crest frequency.

ω_ψ=1

Wherein, f₀It is a resolution parameter, for weighing resolution ratio (the usual feelings of time and frequency in conversion process It is defaulted as 1) under condition.The high then frequency resolution of temporal resolution is reduced, and vice versa.

Step B-2：The wavelet transformation W that checking step B-1 is obtained_s(ω, t) be whether signal s (t) optimal frequency schedule Show, if it is not, then adopting adding window Fourier transformation (Windowed Fourier Transform, WFT) G_s(ω, t), WFT is fixed Justice is：

Wherein, g (t) is the window function of WFT,Fourier for g (t) is converted, and meets condition：Select Gaussian window as the window function of WFT, its expression formula is：

Step B-3 finds out all of ridge curve of the optimal time-frequency representation of signal s (t)HereIt is the h time The ridge curve of harmonic wave.So-called ridge curve is exactly the curve that some Local modulus maximas are formed by connecting on time-frequency figure.Each pole Big value point is just ridge point.

In moment t_i, using following formula algorithm, h maximum point can be found out：

I=1 in above formula, 2 ..., N, N are data length.H_s(ω, t) be signal s (t) optimal time-frequency representation, as W_s (ω, t) or G_s(ω, t).

By H_s(ω, t) in the ridge point line at all moment found out, that is, constitute H vallate curves.

Step B-4：It is x with ridge method reconstruct h order harmonic components^(h)(t)=A^(h)(t)cosφ^(h)(t), wherein A^(h) (t)、φ^(h)T () is respectively the amplitude and phase place of h order harmonic components, ν^(h)T () is the frequency of h order harmonic components, ν^(h)(t) ≡φ′^(h)(t), φ '^(h)T () is φ^(h)The derivative of (t) to time t.

If the optimal time-frequency representation of pending signal s (t) is W_s(ω, t), then can be calculated by below equation To h order harmonic components x of signal s (t)^(h)(t)。

If the optimal time-frequency representation of pending signal s (t) is G_s(ω, t), then calculate the harmonic wave correlation of h time of s (t) Component x^(h)T () is：

In formula,WithIt is by improved discretization factor of influence obtained from parabola interpolation.

Step B-5：The true and false for substituting harmonic component of the method for inspection to extracting using noise immunity differentiates, filters out All of real harmonic component, and stop decomposable process when continuous three harmonic components are judged as fictitious time.Including following Step：

(1) harmonic component is extracted from TFR and corresponding identification statistic D (α is calculated_A, α_v)；

Amplitude A of each harmonic component for extracting^(h)(t) and frequency ν^(h)T the degree of order of () can compose entropy with itWithQuantitatively to weigh, identification statistic D and spectrum entropy Q are defined as follows：

Wherein, α_A, α_vIt is to calculate D (α_A, α_v) weight coefficient.

(2) N is created for original signal s (t)_sIndividual Fourier transform alternate data；

OrderAmplitude keep constant, phase place be changed into being evenly distributed on [0,2 π) on N_sIndividual random phaseIt is this The corresponding Fourier inversion of each phase place of random distribution is a Fourier transform alternate data s ' of s (t)_j(t)。

Here j=1,2 ..., N_s。

(3) TFR related to each alternate data is calculated, and respectively from wherein extracting each harmonic component, so as to Calculate the corresponding identification statistic of each alternate data

Define significance indexIn formulaIt is D_s＞ D₀Alternate data Number；D₀For the degree of order of the h subharmonic of original signal.

Assume to create N_sIndividual alternate data and level of signifiance setup measures are p, i.e. at least p × N_sIndividual alternate data expires Sufficient D_s＞ D₀Can just think that the component is not noise, so as to continue decomposable process.

Judge whether an alternate data meets D_s＞ D₀Method be：Choose the parameter (α of any three groups of different values_A, α_v), three groups of corresponding identifications of parameter are calculated respectively and counts value, if there is any one identification statistics value is more than D₀, then judge The alternate data meets D_s＞ D₀。

(4) the comprehensive measurement value of the degree of correlation between harmonic wave is calculated

Wherein,

In formula, w_A, w_φ, w_vRepresentWeights, here, acquiescence use ρ^(h)≡ρ^(h)(1,1, it is 0) width Degree and phase equalization distribute equal weights, and do not distribute weights to frequency invariance.

(5) in order to reduce the false judgment to true harmonic components, the threshold value of comprehensive measurement value is defined

(6) when a harmonic component meets ρ^(h)≥ρ_minAnd during significance index p >=95%, then it is assumed that this harmonic wave point Amount, by inspection, is real harmonic component.

Step B-6：All of real harmonic component is added so as to obtain a nonlinear model component c₁(t)。

Step B-7：Nonlinear model is cut from original signal s (t) and decomposes the c for obtaining₁(t), and to residual components weight Multiple step B-1 obtains all of each nonlinear model component c to step B-6_i(t)。

Finally, former echo signal s (t) can be expressed as：

In formula, n (t) represents noise.

Step C：Each nonlinear model component c that NMD methods are decomposited_iT () carries out AOK analyses, its process is as follows：

Step C-1：Radially Gaussian kernel function is selected first, and it is defined as follows：

In formula,Represent for controlling radial direction Gaussian function in radial angleThe spread function in direction.

Step C-2：For the optimum core for obtaining changing with signal adaptive, following constrained optimization should be solved and asked Topic.

Constraints is：

Wherein,It is the ambiguity function under polar coordinates,It is corresponding for each nonlinear model componentβ is the volume of optimum core.

Definition A in rectangular coordinate system_i(t；θ, τ) be：

Wherein, s^*(t),w^*U () is s (t), c_i(t), the conjugate complex number form of w (u)；Window function w (u) is with t's Center, width is the symmetrical rhombus window function of 2T, and during | u | ＞ T, w (u)=0；Variable τ and θ are general fuzzy field { τ, u } Parameter, and | τ | ＜ 2T.

Step C-3：(4) in step C-2 are substituted into, (5) formula can be by solving this constrained optimization Problem obtains an optimum kernel functionIt can change over time as short-time ambiguity function.

Step C-4：Calculate current point in time (t) certain nonlinear model component c_iThe adaptive optimal kernel of (t) Time-frequency representation AOK TFR_i：

In formula,For a certain component c of signal_iThe energy value of (t) in some t.

Step C-5：Repeat step C-1 obtains the self adaptation of all nonlinear model components for decompositing most to step C-4 Excellent core time-frequency representation.

Step D：By the AOK TFR of all component of signals for obtaining_i(i.e.) summation, obtain final time-frequency Result P of analysis_NMD-AOK(t, ω)：

The result of NMD+AOK time frequency analysis can be finally given by above formula.

By technical scheme provided by the present invention compared with prior art, analysis below result can be obtained：

Only a certain amount of cross term can be produced with AOK methods, and noiseproof feature is very poor；EMD+AOK and EEMD+AOK side Method can suppressing crossterms to a certain extent, but it is still excessively sensitive to noise；And NMD+AOK methods of the present invention Whether in terms of noise or suppressing crossterms are removed, hence it is evident that better than other several methods.

The above is only the preferred embodiment of the present invention, it should be pointed out that：For the ordinary skill people of the art For member, under the premise without departing from the principles of the invention, some improvements and modifications can also be made, these improvements and modifications also should It is considered as protection scope of the present invention.

Claims (6)

1. it is a kind of based on nonlinear model decompose and adaptive optimal kernel Time-Frequency Analysis Method, it is characterised in that including step：

(1) it is s (t) to define pending signal, and the sample frequency of s (t) is f_s, data length be N；Decomposed using nonlinear model Pending signal s (t) is decomposed into one group of nonlinear model component by method, i.e.,：

$s\left(t\right)=\underset{i}{\Σ}{c}_i\left(t\right)+n\left(t\right)$

Wherein, c_iT () is i-th nonlinear model component of s (t), n (t) represents noise；

(2) each nonlinear model component that step (1) is decomposited is analyzed using adaptive optimal nuclear analytical method, including Step：

(2-1) radially Gaussian kernel function is built

In formula,Represent for controlling radial direction Gaussian function in radial angleThe spread function in direction；

(2-2) to obtain as the optimum core of signal adaptive change is as target problem, building optimization problem model is：

Arrange optimization problem model constraints be：

Wherein,For i-th nonlinear model component c_iT the ambiguity function of () under polar coordinates, β is the body of optimum core Product；Expression formula in rectangular coordinate system is：

$A_i(t,\mathrm{\θ},\mathrm{\τ})={\∫}_{-\mathrm{\∞}}^{+\mathrm{\∞}}{c}_i^{*}(u-\frac{\mathrm{\τ}}{2})w^{*}(u-t-\frac{\mathrm{\τ}}{2})\·c_i(u-\frac{\mathrm{\τ}}{2})w(u-t-\frac{\mathrm{\τ}}{2})e^{j\mathrm{\θ}u}du$

Wherein, A_i(t, θ, τ) isMap amount in rectangular coordinate system, s^*(t),w^*U () is s (t), c_i(t), w The conjugate complex number form of (u)；Window function w (u) is that, with the center of t, width is the symmetrical rhombus window function of 2T, and during | u | ＞ T, w U ()=0, variable τ and θ are the parameters of general fuzzy field { τ, u }, | τ | ＜ 2T；

(2-3) according to A_i(t, θ, τ) solving-optimizing problem model, obtains nonlinear model component c_iThe optimum kernel function of (t) Φ_(i)opt(t, θ, τ)；

(2-4) nonlinear model component c is calculated_iT the time-frequency representation of the adaptive optimal kernel of () is：

$P_{\left(i\right)AOK}(t,\mathrm{\ω})=\frac{1}{4\mathrm{\π}}{\∫}_{-\mathrm{\∞}}^{+\mathrm{\∞}}{\∫}_{-\mathrm{\∞}}^{+\mathrm{\∞}}{A}_i(t,\mathrm{\θ},\mathrm{\τ}){\mathrm{\Φ}}_{\left(i\right)opt}(t,\mathrm{\θ},\mathrm{\τ})e^{-j\mathrm{\θ}t-j\mathrm{\τ}\mathrm{\ω}}d\mathrm{\θ}d\mathrm{\τ}$

In formula, P_(i)AOK(t, ω) is nonlinear model component c_iThe energy value of (t) in t；

(3) time-frequency representation of the adaptive optimal kernel of all nonlinear model components obtained according to step (2), frequency division during calculating The result of analysis is：

$P_{NMD-AOK}(t,\mathrm{\ω})=\underset{i}{\Σ}{P}_{{\left(i\right)}_{AOK}}(t,\mathrm{\ω}).$

2. it is according to claim 1 it is a kind of based on nonlinear model decompose and adaptive optimal kernel Time-Frequency Analysis Method, Characterized in that, in the step (1) using nonlinear model decomposition method by pending signal s (t) be decomposed into one group it is non-linear The method of mode component includes step：

(1-1) wavelet transformation expression formula W of signal s (t) is calculated_s(ω, t)：

$W_s(\mathrm{\ω},t)\≡{\∫}_{-\mathrm{\∞}}^{+\mathrm{\∞}}{s}^{+}\left(u\right){\mathrm{\ψ}}^{*}\[\frac{\mathrm{\ω}(u-t)}{{\mathrm{\ω}}_{\mathrm{\ψ}}}\]\frac{\mathrm{\ω}du}{{\mathrm{\ω}}_{\mathrm{\ψ}}}=\frac{1}{2\mathrm{\π}}{\∫}_0^{\mathrm{\∞}}{e}^{i\mathrm{\ξ}t}\hat{s}\left(\mathrm{\ξ}\right){\widehat{\mathrm{\ψ}}}^{*}\left(\frac{{\mathrm{\ω}}_{\mathrm{\ψ}}\mathrm{\ξ}}{\mathrm{\ω}}\right)d\mathrm{\ξ}$

Wherein,It is the Fourier transform of s (t), i.e.,s⁺T () is the positive frequency portion of s (t) signals Point,ψ (t) is wavelet function, withFourier transform pair each other, and meet condition It is respectively ψ (t),Conjugate complex number；ω_ψSmall echo crest frequency is represented, ω_ψ=1, f₀For resolution parameter, for weighing the resolution of time and frequency in conversion process Rate；

(1-2) judge that whether the calculated wavelet transformation of step (1-1) is the optimal time-frequency representation of s (t)；If judged result is It is no, then calculate adding window Fourier transformation expression formula G of signal s (t)_s(ω, t)：

$G_s(\mathrm{\ω},t)\≡\underset{-\mathrm{\∞}}{\overset{\mathrm{\∞}}{\∫}}{s}^{+}\left(u\right)g(u-t)e^{-i\mathrm{\ω}(u-t)}dt=\frac{1}{2\mathrm{\π}}\underset{0}{\overset{\mathrm{\∞}}{\∫}}{e}^{i\mathrm{\ξ}t}\hat{s}\left(\mathrm{\ξ}\right)\hat{g}(\mathrm{\ω}-\mathrm{\ξ})d\mathrm{\ξ}$

Wherein, g (t) is the window function of adding window Fourier transformation,For the Fourier transformation of g (t), following condition is met：

And

(1-3) all of ridge curve of the optimal time-frequency representation of signal s (t) is found out, and harmonic component is reconstructed with ridge method, wherein, H order harmonic components are：

x^(h)(t)=A^(h)(t)cosφ^(h)(t), h ∈ [1,2 ..., N]

In formula, A^(h)(t)、φ^(h)T () is respectively the amplitude and phase place of h order harmonic components；ν^(h)T () is h order harmonic components Frequency, ν^(h)(t)≡φ^′(h)(t)；φ^′(h)T () is φ^(h)The derivative of (t) to time t；N represents the most high order of harmonic component Number, i.e. data length；

(1-4) true and false for substituting harmonic component of the method for inspection to extracting using noise immunity differentiates, filters out all of Real harmonic component,

(1-5) all of real harmonic component is added, obtains a nonlinear model component c₁(t)；

(1-6) nonlinear model is cut from original signal s (t) and decomposes the c for obtaining₁(t), and to residual components repeat step 1-1 To step 1-6, all of each nonlinear model component c is obtained_i(t)；

Finally, former echo signal s (t) can be expressed as：

$s\left(t\right)=\underset{i}{\Σ}{c}_i\left(t\right)+n\left(t\right)$

In formula, n (t) represents noise.

3. it is according to claim 2 it is a kind of based on nonlinear model decompose and adaptive optimal kernel Time-Frequency Analysis Method, Characterized in that, finding out all of ridge curve of the optimal time-frequency representation of signal s (t) in the step (1-3)Method For：

In moment t_i, h maximum point is found out, and the h maximum point found out is connected, form moment t_iRidge curve be：

${\mathrm{\ω}}_p\left(t\right)=\underset{\mathrm{\ω}\∈\[{\mathrm{\ω}}_{-}\left(t_i\right),{\mathrm{\ω}}_{+}\left(t_i\right)\]}{\arg max}\left|H_s(\mathrm{\ω},t)\right|$

I=1 in above formula, 2 ..., N, N are data length；H_s(ω, t) be signal s (t) optimal time-frequency representation, as W_s(ω, Or G t)_s(ω, t)；H_s(ω, t) in the ridge point line at all moment found out, that is, constitute N vallate curves.

4. it is according to claim 3 it is a kind of based on nonlinear model decompose and adaptive optimal kernel Time-Frequency Analysis Method, Characterized in that, reconstructing h order harmonic components x with ridge method in the step (1-3)^(h)(t)=A^(h)(t)cosφ^(h)T () includes Following steps：

If the optimal time-frequency representation of pending signal s (t) is wavelet transformation,

If the optimal time-frequency representation of pending signal s (t) is windowed FFT,

$\left\{\begin{array}{c}{v}^{\left(h\right)}\left(t\right)={\mathrm{\ω}}_p^{\left(h\right)}\left(t\right)+\mathrm{\δ}{v}_d^{\left(h\right)}\left(t\right)\\ A^{\left(h\right)}{e}^{{\mathrm{i\φ}}^{\left(h\right)}\left(t\right)}=\frac{2G_s\left({\mathrm{\ω}}_p^{\left(h\right)}\right(t),t)}{\hat{g}\[{\mathrm{\ω}}_p^{\left(h\right)}\left(t\right)-v^{\left(h\right)}\left(t\right)\]}\end{array}\right.$

In formula,KindIt is by improved discretization factor of influence obtained from parabola interpolation.

5. it is according to claim 4 it is a kind of based on nonlinear model decompose and adaptive optimal kernel Time-Frequency Analysis Method, Characterized in that, the method for the real harmonic component of screening includes step in the step (1-4)：

(1-4-1) identification statistic D (α is built_A, α_ν), for weighing the amplitude of each harmonic component and the degree of order of frequency；D (α_A, α_v) expression formula be：

$D({\mathrm{\α}}_A,{\mathrm{\α}}_v)\≡{\mathrm{\α}}_{A}Q\[{\hat{A}}^{\left(h\right)}\left(\mathrm{\ξ}\right)\]+{\mathrm{\α}}_{v}Q\[{\hat{v}}^{\left(h\right)}\left(\mathrm{\ξ}\right)\],$

$Q\[f\left(x\right)\]\≡-\∫\frac{|f\left(x\right){|}^2}{\∫|f\left(x\right){|}^{2}dx}\ln \frac{|f\left(x\right){|}^2}{\∫|f\left(x\right){|}^{2}dx}dx$

In formula,WithA is represented respectively^(h)(t) and v^(h)The spectrum entropy of (t)；α_A, α_vFor weight coefficient；

(1-4-2) N is created for original signal s (t)_sIndividual Fourier transform alternate data；Wherein, j-th Fourier transform alternate data Expression formula be：

$s_j^{\′}\left(t\right)=\frac{1}{2\mathrm{\π}}\∫e^{i\mathrm{\ξ}t}\left|\hat{s}\left(\mathrm{\ξ}\right)\right|e^{{\mathrm{i\φ}}_{{\mathrm{\ξ}}_j}}d\mathrm{\ξ},j=1,2,...,N_s;$

(1-4-3) the optimal time-frequency representation formula of each Fourier transform alternate data is calculated, and respectively from wherein extracting Each harmonic component, calculating the corresponding identification statistic of each Fourier transform alternate data is：

$D_j({\mathrm{\α}}_A,{\mathrm{\α}}_v)\≡{\mathrm{\α}}_{A}Q\[{\hat{A}}^{\left(j\right)}\left(\mathrm{\ξ}\right)\]+{\mathrm{\α}}_{v}Q\[{\hat{v}}^{\left(j\right)}\left(\mathrm{\ξ}\right)\],j=1,2,...,N_s$

Define significance indexIn formulaTo meet D_s＞ D₀Replacement The number of data；D₀For the degree of order of the h subharmonic of original signal；

It is p to arrange level of signifiance index, and when there is x alternate data D is met_s＞ D₀And x >=p × N_sWhen, judge corresponding harmonic wave Component is not noise；Otherwise, it is determined that corresponding harmonic component is noise；

(1-4-4) the comprehensive measurement value of the degree of correlation between harmonic wave is calculated

${\mathrm{\ρ}}^{\left(h\right)}(w_A,w_{\mathrm{\φ}},w_v)={\left(q_A^{\left(h\right)}\right)}^{w_A}{\left(q_{\mathrm{\φ}}^{\left(h\right)}\right)}^{w_{\mathrm{\φ}}}{\left(q_v^{\left(h\right)}\right)}^{w_v}$

Wherein,

$q_A^{\left(h\right)}\&equiv;\exp \{-\frac{\sqrt{<{\&lsqb;A^{\left(h\right)}\left(t\right)<A^{\left(1\right)}\left(t\right)>-A^{\left(1\right)}\left(t\right)<A^{\left(h\right)}\left(t\right)>\&rsqb;}^2>}}{<A^{\left(1\right)}\left(t\right)A^{\left(h\right)}\left(t\right)>}\},$

$q_{\mathrm{\&phi;}}^{\left(h\right)}\&equiv;a|<\exp \{i\&lsqb;{\mathrm{\&phi;}}^{\left(h\right)}\left(t\right)-{\mathrm{h\&phi;}}^{\left(1\right)}\left(t\right)\&rsqb;\}>|,$

$q_v^{\left(h\right)}\&equiv;\exp \{-\frac{\sqrt{<{\&lsqb;v^{\left(h\right)}\left(t\right)-{\mathrm{hv}}^{\left(1\right)}\left(t\right)\&rsqb;}^2>}}{<v^{\left(h\right)}\left(t\right)>}\}.$

In formula, w_A, w_φ, w_vRepresentWeights, ρ^(h)It is that amplitude and phase equalization distribute equal weights ρ^(h)≡ ρ^(h)(1,1,0)；

(1-4-5) threshold value of comprehensive measurement value is defined

(1-4-6) when a harmonic component meets ρ^(h)≥ρ_minAnd significance index p 95% when, judge this harmonic component It is real harmonic component by inspection.

6. it is according to claim 5 it is a kind of based on nonlinear model decompose and adaptive optimal kernel Time-Frequency Analysis Method, Characterized in that, judging whether an alternate data meets D in the step (1-4-3)_s＞ D₀Method be：

Choose the parameter (α of any three groups of different values_A, α_v), three groups of corresponding identifications of parameter are calculated respectively and counts value, if existing Any one identification statistics value is more than D₀, then judge that the alternate data meets D_s＞ D₀。