CN106546818A - A kind of harmonic signal detection method based on DNL Mode Decomposition - Google Patents

Abstract

The invention discloses a kind of harmonic signal detection method based on DNL Mode Decomposition, difference algorithm is first applied to primary signal by this new harmonic detecting method, then obtains a series of significant nonlinear model components using nonlinear model decomposition algorithm；To nonlinear model decomposition being carried out again after each nonlinear model component integration, and extract first nonlinear model component as the nonlinear model component of primary signal；Finally, the various harmonic signals of signal are extracted with Hilbert marginal spectrum.The noise robustness characteristic of nonlinear model decomposition method is inherited based on the harmonic signal detection method of DNL Mode Decomposition, and the interference of chaos and noise signal can be suppressed, it is advantageous in terms of higher harmonic components by a small margin are extracted.

Description

A kind of harmonic signal detection method based on DNL Mode Decomposition

Technical field

The present invention relates to Power estimation and signal processing applications field, more particularly to it is a kind of based on DNL Mode Decomposition Harmonic signal detection method.

Background technology

Adaptive Time Frequency Analysis method is widely used in numerous areas, such as speech signal analysis, Signal processing of sonar With mechanical fault diagnosis etc..There are many observable chaotic signals in natural phenomena, such as the bad border noise in ocean, by electromagnetic pulse Clutter produced in ocean surface etc..Signal is detected and analyzed with certain difficulty under Chaotic Background.

Fourier transformation has simple and calculates efficient advantage, is to analyze the most frequently used algorithm of harmonic wave.But, it is present Three major defects：Aliasing, fence effect and spectrum leakage, are not suitable for analyzing non-stationary signal.Based on Radial Basis Function neural The method of network can also be used for the harmonic amplitude for detecting measured signal.However, the major defect of neural network algorithm is to train Journey is complicated, convergence rate is slow and is easily trapped into local minimum.

Hilbert-Huang transform (Hilbert-Huang Transform, HHT) is a kind of new Non-stationary Signal Analysis Sophisticated signal is decomposed into a series of by method using empirical mode decomposition (Empirical Mode Decomposition, EMD) Intrinsic mode functions (Intrinsic Mode Function, IMF), further make Hilbert conversion, are closed to each IMF In the temporal frequency Joint Distribution of primary signal.EMD is widely used in non-linear and Non-stationary Signal Analysis.But, in reality Using in, EMD occurs mode mixing phenomenon, and which shows two aspects：One is that to contain yardstick in same IMF components poor Different larger component of signal；Two is that the component of signal of same yardstick has been occurred in different IMF.Ask to solve this Topic, adds white Gaussian noise to signal and applies EMD, then, gathers empirical mode decomposition (Ensemble Empirical Mode Decomposition, EEMD) arise at the historic moment.Difference empirical mode decomposition (Differential Empirical Mode Decomposition, DEMD) then EMD is applied to primary signal calculus of differences first, EMD can be extracted and EEMD cannot Detached higher harmonic components by a small margin.However, all there is the defect to noise-sensitive in three of the above method.

Signal decomposition can be significant nonlinear model point on series of physical by nonlinear model decomposition method (NMD) Amount, while eliminating noise.Relative to past method, NMD makes an uproar because parameter can be adaptive selected with fabulous Sound robust property.But, performances of the NMD when harmonic signal by a small margin is detected not is fine, is particularly deposited in Chaotic Background In case, its performance is general.

The content of the invention

It is an object of the invention to provide a kind of harmonic signal detection method based on DNL Mode Decomposition, the party Method combines the advantage of existing method, overcomes empirical mode decomposition method, set empirical mode decomposition method and difference experience Mode decomposition method is to chaos and noise signal sensitivity and nonlinear model decomposition method performance in small-signal harmonic detecting Bad shortcoming.

For achieving the above object, the technical solution used in the present invention is：It is a kind of based on the humorous of DNL Mode Decomposition Ripple signal detecting method, comprises the following steps：

Step A：Prepare primary signal s (t) of process to be detected, its sample rate is f_s, data length is N；

Step B：Calculus of differences is carried out to primary signal s (t) and obtains new signal s ' (t)；

Step C：Nonlinear model is carried out to signal s ' (t) and decomposes acquisition nonlinear model component c_i′(t)；

Step C-1：The wavelet transformation W of signal calculated s ' (t)_s′(ω, t), wavelet transformation is defined as：

Wherein,It is the Fourier transform of s ' (t)；s′⁺T () is the positive frequency part of signal s ' (t), its expression formula For：Wavelet functions of the ψ (t) for wavelet transformation,For the Fourier transform of ψ (t), condition is metSubscript * represents conjugate operation；It is the crest frequency of small echo；Wavelet function is just adopting logarithm State is distributed small echo, itsAnd ω_ψIt is expressed as follows：

Wherein, f₀It is the resolution parameter for weighing time and frequency resolution in conversion process, generally gives tacit consent to f₀=1；

Step C-2：Check whether wavelet transformation is optimal time-frequency representation, if it is not, then adopting adding window Fourier transformation G_s′ (ω, t), adding window Fourier transform definition is as follows：

Wherein, g (t) is the window function of adding window Fourier transformation,For the Fourier transform of g (t), condition is met：Select Gaussian window as the window function of adding window Fourier transformation, its expression formula is：

Step C-3：Find out all of ridge curve in the time-frequency representation of signal s ' (t)Here defineIt is h The ridge curve of subharmonic；

In sometime t_n, using following formula algorithm, h maximum point can be found out；

N=1,2 in above formula ..., N；N is data length；H_s′(ω is t) in above-mentioned Fourier transformation or adding window Fu Time-frequency representation after leaf transformation, as W_s′(ω, t) or G_s′(ω,t)；

By H_s′(ω, t) in the ridge point line at all moment found out, may make up h vallate curves

Step C-4：Using ridge curve Reconstruction h order harmonic componentsWherein A^(h)(t)、WithIt is its amplitude, phase place and frequency respectively；Calculating process is as follows：

If the time-frequency representation of signal s ' (t) adopts Fourier transformation, h order harmonic components x^(h)T () is obtained by formula (1) Arrive：

If the time-frequency representation of signal s ' (t) adopts adding window Fourier transformation, h order harmonic components x^(h)T () is by formula (2) obtain：

Wherein,WithRespectively Fourier transformation and adding window Fourier transformation are produced by parabola interpolation The amendment that raw discretization affects；

Step C-5：The method of inspection is substituted with noise immunity and determines effective harmonic component；The method is differentiated using alternate data The true and false of the harmonic component for extracting, filters out all of real harmonic component, and when continuous three harmonic components are judged to Break and stop catabolic process for fictitious time；Comprise the following steps that：

(1) calculate identification statistic D of a certain harmonic component for extracting₀(α_A,α_ν)；

Amplitude A of each harmonic component for extracting^(h)(t) and frequency ν^(h)T the degree of order of () can compose entropy with whichWithQuantitatively to weigh, wherein,WithRespectively A^(h)(t) and ν^(h)The Fourier of (t) Leaf transformation, identification statistic D are defined as follows：

Wherein, α_AAnd α_vRespectivelyWithWeight coefficient；

(3) N is created for signal s ' (t)_sIndividual Fourier transform alternate data, its production method is：

Wherein, φ_ξObey [0,2 π) on be uniformly distributed, each φ_ξOne Fourier transform alternate data of correspondence；

(3) time-frequency representation corresponding with each alternate data is calculated, and therefrom extracts each harmonic component respectively, counted Calculate the identification statistic of each alternate dataHere defining significance level index is：

In formula,To meet D_s>D₀Alternate data number；Assume to create N_sIndividual alternate data and by significance Horizontal setup measures are p, i.e. at least N_s× p alternate data meets D_s>D₀Just think that the component is not noise, so as to continue Catabolic process；Parameter (α of this detection using three groups of different values_A,α_ν), i.e., calculate respectively D (1,1), D (0,1) and D (1,0) Value, as long as wherein at least one value does not meet null hypothesises, then it is assumed that meet D_s>D₀；

(4) calculate the comprehensive measurement value of the degree of association between harmonic wave

Wherein,

In formula, w_A,w_φ,w_vRepresentWeights；Acquiescence uses ρ^(h)≡ρ^(h)(1,1,0) it is amplitude and phase place Concordance distributes equal weights, and does not distribute weights to frequency invariance；

(5) in order to reduce the false judgment to true harmonic components, the threshold value for defining comprehensive measurement value is：

(6) when the comprehensive measurement value index of a harmonic component meets ρ^(h)≥ρ_minAnd significance level index During significance-level >=p=95%, then it is assumed that this harmonic component, by inspection, is real harmonic component；If Can not be substituted by noise immunity and be checked, then be stopped nonlinear model and decompose；

Step C-6：All of real harmonic component is added and constitutes a nonlinear model component c₁′(t)；

Step C-7：The nonlinear model component is deducted from signal s ' (t), repeat step C-1 obtains institute to step C-6 Some nonlinear model component c_i′(t)；

Step D：To the nonlinear model component c for obtaining_i' (t) integrations obtain b_i(t)；

Step E：To b_iT () carries out nonlinear model decomposition with nonlinear model decomposition method again；

Step E-1：By the process of step C-1 to step C-6 to each b_iT () carries out nonlinear model decomposition；

Step E-2：Extract each b_iT first nonlinear model component of () is used as the non-linear of primary signal s (t) Mode component c_i(t)；Finally obtain the nonlinear model decomposition result of original signal：

Step F：S (t) to obtaining carries out analysis of spectrum, extracts harmonic signal；

Step F-1：To c_iT () does Hilbert transform and produces quadrature component

Step F-2：Construction complex signal z_iT (), its expression formula is：

Step F-3：By complex signal z_iT () is converted into polar form, try to achieve instantaneous envelope a_i(t) and instantaneous frequency ω_i (t)；z_iT the polar form of () is：

Wherein, instantaneous envelope a_i(t) and instantaneous phase φ_iT () is expressed as follows：

Its instantaneous frequency ω_iT () is instantaneous phase φ_iT the derivative of () to the time, is expressed as：

Step F-4：Seek all nonlinear model component c_i(t) hilbert spectrum H (ω, t)；

Hilbert spectrum is defined as：

Step F-5：Hilbert marginal spectrum is sought, the spectral peak of signal harmonic component is extracted；

The expression formula of Hilbert marginal spectrum is：

Wherein, T represents the sampling time, and its computing formula is：T=N/f_s。

Compared with prior art, its remarkable advantage is the present invention：

1. DNL Mode Decomposition method can more effectively suppress noise and the interference produced by chaotic signal；

2., for the signal with multiple harmonic components, DNL Mode Decomposition method can extract wherein amplitude Less harmonic component.

Description of the drawings

Fig. 1 is the time domain waveform of chaotic signal d (t) and primary signal s (t)；

The Hilbert marginal spectrum that Fig. 2 is obtained using DNL Mode Decomposition method when being SNR=25dB；

The Hilbert marginal spectrum that Fig. 3 is obtained using nonlinear model decomposition method when being SNR=25dB；

The Hilbert marginal spectrum that Fig. 4 is obtained using empirical mode decomposition method when being SNR=25dB；

The Hilbert marginal spectrum obtained using set empirical mode decomposition method when Fig. 5 is SNR=25dB；

The Hilbert marginal spectrum that Fig. 6 is obtained using difference empirical mode decomposition method when being SNR=25dB.

Specific embodiment

Below by taking the primary signal added with white Gaussian noise and Chaotic Background as an example, with reference to accompanying drawing, the present invention is described in detail Embodiment.

The present invention carries out calculus of differences to primary signal s (t) first and obtains new signal, then carries out nonlinear model to which Formula decomposition obtains one group of nonlinear model component c_i′(t)；Nonlinear model component to obtaining is integrated and obtains b_i(t)；Again Nonlinear model decomposition is carried out to which, b is extracted_iT first nonlinear model component of () is used as the non-linear of primary signal Mode component, is designated as c_i(t)；Finally, using spectral analysis technology field of signal processing advantage, it is non-linear to what is finally given Mode component carries out Hilbert limit analysis of spectrum, realizes the detection of harmonic signal.

Assume that primary signal s (t) is：

S (t)=cos (2 π ft)+0.3cos (2 π 3ft)+0.07cos (2 π 5ft)+n (t)+d (t)

Wherein, n (t) is white Gaussian noise, chaotic signal d (t)=0.05x (t).Here, if fundamental frequency f=30Hz, Sample frequency f_s=600Hz, data length are N=650.

X (t) is produced by typical Lorenz chaos systems, and its system model is described as follows：

Wherein, σ=10, r=28, b=8/3, initial value x₀=y₀=z₀=0.1.D (t) and s (t) time domain waveforms such as Fig. 1 It is shown.

The present invention is that (Differential Nonlinear Mode Decomposition, difference are non-based on DNMD for one kind Linear model decompose) harmonic signal detection method, comprise the steps：

Step A：Prepare primary signal s (t) of process to be detected, its sample rate is f_s=600Hz, data length are N= 650；

Step B：Calculus of differences is carried out to primary signal s (t) and obtains new signal s ' (t)；

Step C：Nonlinear model is carried out to signal s ' (t) and decomposes acquisition nonlinear model component c_i′(t)；

Step C-1：The wavelet transformation W of signal calculated s ' (t)_s′(ω, t), wavelet transformation is defined as：

Wherein,It is the Fourier transform of s ' (t)；s′⁺T () is the positive frequency part of signal s ' (t), its expression formula For：Wavelet functions of the ψ (t) for wavelet transformation,For the Fourier transform of ψ (t), condition is metSubscript * represents conjugate operation；It is the crest frequency of small echo；Wavelet function is just adopting logarithm State is distributed small echo, itsAnd ω_ψIt is expressed as follows：

Wherein, f₀It is the resolution parameter for weighing time and frequency resolution in conversion process, generally gives tacit consent to f₀=1；

Step C-2：Check whether wavelet transformation is optimal time-frequency representation, if it is not, then adopting adding window Fourier transformation G_s′ (ω, t), adding window Fourier transform definition is as follows：

Wherein, g (t) is the window function of adding window Fourier transformation,For the Fourier transform of g (t), condition is met：Select Gaussian window as the window function of adding window Fourier transformation, its expression formula is：

Step C-3：Find out all of ridge curve in the time-frequency representation of signal s ' (t)When so-called ridge curve is exactly Frequency schemes the curve are formed by connecting by upper some Local modulus maximas, and each maximum point is just ridge point；Here defineIt is h The ridge curve of subharmonic；

In sometime t_n, using following formula algorithm, h maximum point can be found out；

N=1,2 in above formula ..., N；N is data length；H_s′(ω is t) in above-mentioned Fourier transformation or adding window Fu Time-frequency representation after leaf transformation, as W_s′(ω, t) or G_s′(ω,t)；

By H_s′(ω, t) in the ridge point line at all moment found out, may be constructed h vallate curves

Step C-4：Using ridge curve Reconstruction h order harmonic componentsWherein A^(h)(t)、WithIt is its amplitude, phase place and frequency respectively；Calculating process is as follows：

If the time-frequency representation of signal s ' (t) adopts Fourier transformation, h order harmonic components x^(h)T () is obtained by formula (1) Arrive：

If the time-frequency representation of signal s ' (t) adopts adding window Fourier transformation, h order harmonic components x^(h)T () is by formula (2) obtain：

Wherein,WithRespectively Fourier transformation and adding window Fourier transformation are produced by parabola interpolation The amendment that raw discretization affects；

Step C-5：The method of inspection is substituted with noise immunity and determines effective harmonic component；The method is differentiated using alternate data The true and false of the harmonic component for extracting, filters out all of real harmonic component, and when continuous three harmonic components are judged to Break and stop catabolic process for fictitious time；Comprise the following steps that：

(1) calculate identification statistic D of a certain harmonic component for extracting₀(α_A,α_ν)；

Amplitude A of each harmonic component for extracting^(h)(t) and frequency ν^(h)T the degree of order of () can compose entropy with whichWithQuantitatively to weigh, wherein,WithRespectively A^(h)(t) and ν^(h)The Fourier of (t) Leaf transformation, identification statistic D are defined as follows：

Wherein, α_AAnd α_vRespectivelyWithWeight coefficient；

(4) N is created for signal s ' (t)_sIndividual Fourier transform alternate data, its production method is：

Wherein, φ_ξObey [0,2 π) on be uniformly distributed, each φ_ξOne Fourier transform alternate data of correspondence；

(3) time-frequency representation corresponding with each alternate data is calculated, and therefrom extracts each harmonic component respectively, counted Calculate the identification statistic of each alternate dataHere defining significance level index is：

In formula,To meet D_s>D₀Alternate data number；Assume to create N_sIndividual alternate data and by significance Horizontal setup measures are p, i.e. at least N_s× p alternate data meets D_s>D₀Just think that the component is not noise, so as to continue Catabolic process；Parameter (α of this detection using three groups of different values_A,α_ν), i.e., calculate respectively D (1,1), D (0,1) and D (1,0) Value, as long as wherein at least one value does not meet null hypothesises, then it is assumed that meet D_s>D₀；

(4) calculate the comprehensive measurement value of the degree of association between harmonic wave

Wherein,

In formula, w_A,w_φ,w_vRepresentWeights；Here, acquiescence uses ρ^(h)≡ρ^(h)(1,1,0) it is 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 for defining comprehensive measurement value is：

(6) when the comprehensive measurement value index of a harmonic component meets ρ^(h)≥ρ_minAnd significance level index During significance-level >=p=95%, then it is assumed that this harmonic component, by inspection, is real harmonic component；If Can not be substituted by noise immunity and be checked, then be stopped nonlinear model and decompose；

Step C-6：All of real harmonic component is added and constitutes a nonlinear model component c₁′(t)；

Step C-7：The nonlinear model component is deducted from signal s ' (t), repeat step C-1 obtains institute to step C-6 Some nonlinear model component c_i′(t)；

Step D：To the nonlinear model component c for obtaining_i' (t) integrations obtain b_i(t)；

Step E：To b_iT () carries out nonlinear model decomposition with nonlinear model decomposition method again；

Step E-1：By the process of step C-1 to step C-6 to each b_iT () carries out nonlinear model decomposition；

Step E-2：Extract each b_iT first nonlinear model component of () is used as the non-linear of primary signal s (t) Mode component c_i(t)；Finally obtain the nonlinear model decomposition result of original signal：

Step F：S (t) to obtaining carries out analysis of spectrum, extracts harmonic signal；

Step F-1：To c_iT () does Hilbert transform and produces quadrature component

Step F-2：Construction complex signal z_iT (), its expression formula is：

Step F-3：By complex signal z_iT () is converted into polar form, try to achieve instantaneous envelope a_i(t) and instantaneous frequency ω_i (t)；z_iT the polar form of () is：

Wherein, instantaneous envelope a_i(t) and instantaneous phase φ_iT () is expressed as follows：

Its instantaneous frequency ω_iT () is instantaneous phase φ_iT the derivative of () to the time, is expressed as：

Step F-4：Seek all nonlinear model component c_i(t) hilbert spectrum H (ω, t)；

Hilbert spectrum is defined as：

Step F-5：Hilbert marginal spectrum is sought, the spectral peak of signal harmonic component is extracted；

The expression formula of Hilbert marginal spectrum is：

Wherein, T represents the sampling time, and its computing formula is：T=N/f_s。

Fig. 2 to Fig. 6 is followed successively by as signal to noise ratio snr=25dB, and DNL Mode Decomposition method, non-thread is respectively adopted Sexual norm decomposition method, empirical mode decomposition method, set empirical mode decomposition method and difference empirical mode decomposition method are obtained The Hilbert marginal spectrum for arriving.From figure 2 it can be seen that DNL Mode Decomposition method clearly can extract it is original The fundamental frequency (30Hz) of signal, and three times and quintuple harmonics frequency (90Hz, 150Hz), and the spectrum of each frequency component Peak is very sharp；And nonlinear model decomposition method (shown in Fig. 3) has only extracted fundametal compoment (30Hz) and triple-frequency harmonics point Amount (90Hz), fails to extract quintuple harmonics (150Hz) component；Find out from Fig. 4-Fig. 6, empirical mode decomposition method, set Empirical mode decomposition method and difference empirical mode decomposition method are serious by the interference effect of noise and chaos, are only capable of extracting The fundametal compoment (30Hz) of signal, it is impossible to accurately extract other harmonic components of signal, only at third harmonic frequencies (90Hz) Nearby there is a certain degree of projection.

Quantitatively to compare the performance of difference nonlinear model decomposition method and other methods, parameter R_k, it is defined asWherein, number of times of the k for harmonic wave, A_kIt is the amplitude of correspondence harmonic frequency,It is to be decomposed point The meansigma methodss of the spectrum of amount.R of every kind of method under different signal to noise ratios_kValue is as shown in table 1.

R of the 1 every kind of method of table under different signal to noise ratios_kValue

In table 1, DNMD is DNL Mode Decomposition method, and NMD is nonlinear model decomposition method, and EMD is Empirical Mode State decomposition method, EEMD set empirical mode decomposition methods and DEMD are difference empirical mode decomposition method, and SNR is signal to noise ratio.

Can be seen that from 1 column data of table：As SNR=20dB, DNL Mode Decomposition method can be extracted The fundamental frequency (30Hz) of primary signal, three times and quintuple harmonics frequency (90Hz, 150Hz), nonlinear model decomposition method is only Fundametal compoment (30Hz) and third-harmonic component (90Hz) are extracted, and empirical mode decomposition method, set empirical modal have divided Solution method and difference empirical mode decomposition method have only extracted fundametal compoment (30Hz).Although in SNR=15dB, difference is non- Linear model decomposition method is lost the ability for extracting quintuple harmonics (150Hz), but its R₁=52.74dB and R₃=28.30dB It is all higher than the R of nonlinear model decomposition method₁=51.90dB and R₃=26.28dB.

In sum, above example is not intended to limit the guarantor of the present invention only to illustrate technical scheme Shield scope.All any modification, equivalent substitution and improvements within the spirit and principles in the present invention, done etc., which all should be covered In the middle of scope of the presently claimed invention.

Claims (1)

1. a kind of harmonic signal detection method based on DNL Mode Decomposition, it is characterised in that：Comprise the following steps：

Step A：Prepare primary signal s (t) of process to be detected, its sample rate is f_s, data length is N；

Step B：Calculus of differences is carried out to primary signal s (t) and obtains new signal s ' (t)；

Step C：Nonlinear model is carried out to signal s ' (t) and decomposes acquisition nonlinear model component c_i′(t)；

Step C-1：The wavelet transformation W of signal calculated s ' (t)_s′(ω, t), wavelet transformation is defined as：

$\begin{array}{c}{W}_{s^{\′}}(\mathrm{\ω},t)=\underset{-\mathrm{\∞}}{\overset{\mathrm{\∞}}{\∫}}{s}^{\′+}\left(u\right){\mathrm{\ψ}}^{*}\[\frac{\mathrm{\ω}(u-t)}{{\mathrm{\ω}}_{\mathrm{\ψ}}}\]\frac{\mathrm{\ω}du}{{\mathrm{\ω}}_{\mathrm{\ψ}}}\\ =\frac{1}{2\mathrm{\π}}\underset{0}{\overset{\mathrm{\∞}}{\∫}}{e}^{i\mathrm{\ξ}t}{\hat{s}}^{\′}\left(\mathrm{\ξ}\right){\widehat{\mathrm{\ψ}}}^{*}\left(\frac{{\mathrm{\ω}}_{\mathrm{\ψ}}\mathrm{\ξ}}{\mathrm{\ω}}\right)d\mathrm{\ξ}\end{array}$

Wherein,It is the Fourier transform of s ' (t)；s^′+T () is the positive frequency part of signal s ' (t), its expression formula is：Wavelet functions of the ψ (t) for wavelet transformation,For the Fourier transform of ψ (t), condition is metSubscript * represents conjugate operation；It is the crest frequency of small echo；Wavelet function is just adopting logarithm State is distributed small echo, itsAnd ω_ψIt is expressed as follows：

$\widehat{\mathrm{\ψ}}\left(\mathrm{\ξ}\right)=e^{-{\left(2{\mathrm{\πf}}_0\ln \mathrm{\ξ}\right)}^2/2},{\mathrm{\ω}}_{\mathrm{\ψ}}=0$

Wherein, f₀It is the resolution parameter for weighing time and frequency resolution in conversion process, generally gives tacit consent to f₀=1；

Step C-2：Check whether wavelet transformation is optimal time-frequency representation, if it is not, then adopting adding window Fourier transformation G_s′(ω, T), adding window Fourier transform definition is as follows：

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

Wherein, g (t) is the window function of adding window Fourier transformation,For the Fourier transform of g (t), condition is met：Select Gaussian window as the window function of adding window Fourier transformation, its expression formula is：

$\hat{g}\left(\mathrm{\ξ}\right)=e^{-{\left(f_0\mathrm{\ξ}\right)}^2/2}\⇔g\left(t\right)=\frac{1}{\sqrt{2\mathrm{\π}}{f}_0}{e}^{-{(t/f_0)}^2/2}$

Step C-3：Find out all of ridge curve in the time-frequency representation of signal s ' (t)DefinitionIt is the ridge of h subharmonic Curve；

In sometime t_n, using following formula algorithm, find out h maximum point；

${\mathrm{\ω}}_p\left(t\right)=\underset{\mathrm{\ω}\∈\[{\mathrm{\ω}}_{-}\left(t_n\right),{\mathrm{\ω}}_{+}\left(t_n\right)\]}{\mathrm{argmax}}\left|H_{s^{\′}}(\mathrm{\ω},t)\right|$

N=1,2 in above formula ..., N；N is data length；H_s′(ω t) is become through above-mentioned Fourier transformation or adding window Fourier Time-frequency representation after changing, as W_s′(ω, t) or G_s′(ω,t)；

By H_s′(ω, t) in the ridge point line at all moment found out, constitute h vallate curves

Step C-4：Using ridge curve Reconstruction h order harmonic componentsWherein A^(h)(t)、It is its amplitude, phase place and frequency respectively；Calculating process is as follows：

If the time-frequency representation of signal s ' (t) adopts Fourier transformation, h order harmonic components x^(h)T () is obtained by formula (1)：

If the time-frequency representation of signal s ' (t) adopts adding window Fourier transformation, h order harmonic components x^(h)T () is obtained by formula (2) Arrive：

Wherein,WithRespectively Fourier transformation and adding window Fourier transformation are by produced by parabola interpolation The amendment that discretization affects；

Step C-5：The method of inspection is substituted with noise immunity and determines effective harmonic component；The method differentiates to extract using alternate data The true and false of the harmonic component for going out, filters out all of real harmonic component, and when continuous three harmonic components are judged as Fictitious time stops catabolic process；Comprise the following steps that：

(1) calculate identification statistic D of a certain harmonic component for extracting₀(α_A,α_ν)；

Amplitude A of each harmonic component for extracting^(h)(t) and frequency ν^(h)T the degree of order of () composes entropy with whichWithQuantitatively to weigh, wherein,WithRespectively A^(h)(t) and ν^(h)The Fourier transform of (t), identification Statistic D is defined as follows：

$D({\mathrm{\α}}_A,{\mathrm{\α}}_{\mathrm{\ν}})\≡{\mathrm{\α}}_{A}Q\[{\hat{A}}^{\left(h\right)}\left(\mathrm{\ξ}\right)\]+{\mathrm{\α}}_{\mathrm{\ν}}Q\[{\widehat{\mathrm{\ν}}}^{\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.$

Wherein, α_AAnd α_vRespectivelyWithWeight coefficient；

(2) N is created for signal s ' (t)_sIndividual Fourier transform alternate data, its production method is：

$y\left(t\right)=\frac{1}{2\mathrm{\π}}\∫e^{-i\mathrm{\ξ}t}\left|{\hat{s}}^{\′}\left(\mathrm{\ξ}\right)\right|e^{{\mathrm{i\φ}}_{\mathrm{\ξ}}}d\mathrm{\ξ}$

Wherein, φ_ξObey [0,2 π) on be uniformly distributed, each φ_ξOne Fourier transform alternate data of correspondence；

(3) time-frequency representation corresponding with each alternate data is calculated, and therefrom extracts each harmonic component respectively, calculated The identification statistic of each alternate dataDefining significance level index is：

$significance-level=\frac{N_{D_s>D_0}}{N_s}$

In formula,To meet D_s>D₀Alternate data number；Assume to create N_sIndividual alternate data and by significance level Setup measures are p, i.e. at least N_s× p alternate data meets D_s>D₀Just think that the component is not noise, so as to continue to decompose Process；Parameter (α of the detection using three groups of different values_A,α_ν), i.e., calculate respectively D (1,1), D (0,1) and D (1, value 0), As long as wherein at least one value does not meet null hypothesises, then it is assumed that meet D_s>D₀；

(4) calculate the comprehensive measurement value of the degree of association between harmonic wave

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

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_{\mathrm{\&nu;}}^{\left(h\right)}\&equiv;\exp \{-\frac{\sqrt{<{\&lsqb;{\mathrm{\&nu;}}^{\left(h\right)}\left(t\right)-{\mathrm{h\&nu;}}^{\left(1\right)}\left(t\right)\&rsqb;}^2>}}{<{\mathrm{\&nu;}}^{\left(h\right)}\left(t\right)>}\}.$

In formula, w_A,w_φ,w_vRepresentWeights；Acquiescence uses ρ^(h)≡ρ^(h)(1,1,0) it is consistent with phase place for amplitude Property the equal weights of distribution, and do not distribute weights to frequency invariance；

(5) in order to reduce the false judgment to true harmonic components, the threshold value for defining comprehensive measurement value is：

${\mathrm{\&rho;}}_{min}={0.5}^{w_A+w_{\mathrm{\&phi;}}}$

(6) when the comprehensive measurement value index of a harmonic component meets ρ^(h)≥ρ_minAnd significance level index During significance-level >=p=95%, then it is assumed that this harmonic component, by inspection, is real harmonic component；If Can not be substituted by noise immunity and be checked, then be stopped nonlinear model and decompose；

Step C-6：All of real harmonic component is added and constitutes a nonlinear model component c₁′(t)；

Step C-7：The nonlinear model component is deducted from signal s ' (t), repeat step C-1 obtains all of to step C-6 Nonlinear model component c_i′(t)；

Step D：To the nonlinear model component c for obtaining_i' (t) integrations obtain b_i(t)；

Step E：To b_iT () carries out nonlinear model decomposition with nonlinear model decomposition method again；

Step E-1：By the process of step C-1 to step C-6 to each b_iT () carries out nonlinear model decomposition；

Step E-2：Extract each b_iNonlinear model of the first nonlinear model component of (t) as primary signal s (t) Component c_i(t)；Finally obtain the nonlinear model decomposition result of original signal：

$s\left(t\right)=\underset{i}{\mathrm{\&Sigma;}}{c}_i\left(t\right)+n\left(t\right)$

Step F：Primary signal s (t) to obtaining carries out analysis of spectrum, extracts harmonic signal；

Step F-1：To nonlinear model component c_iT () does Hilbert transform and produces quadrature component

Step F-2：Construction complex signal z_iT (), its expression formula is：

$z_i\left(t\right)=c_i\left(t\right)+j{\tilde{c}}_i\left(t\right)$

Step F-3：By complex signal z_iT () is converted into polar form, try to achieve instantaneous envelope a_i(t) and instantaneous frequency ω_i(t)；z_i T the polar form of () is：

$z_i\left(t\right)=a_i\left(t\right)e^{{\mathrm{j\&phi;}}_i\left(t\right)}$

Wherein, instantaneous envelope a_i(t) and instantaneous phase φ_iT () is expressed as follows：

$a_i\left(t\right)=\sqrt{{c_i}^2\left(t\right)+{{\tilde{c}}_i}^2\left(t\right)}$

${\mathrm{\&phi;}}_i\left(t\right)=\arctan \&lsqb;\frac{{\tilde{c}}_i\left(t\right)}{c_i\left(t\right)}\&rsqb;$

Its instantaneous frequency ω_iT () is instantaneous phase φ_iT the derivative of () to the time, is expressed as：

${\mathrm{\&omega;}}_i\left(t\right)=\frac{{\mathrm{d\&phi;}}_i\left(t\right)}{dt}$

Step F-4：Seek all nonlinear model component c_i(t) hilbert spectrum H (ω, t)；

Hilbert spectrum is defined as：

$H(\mathrm{\&omega;},t)=\mathrm{Re}\&lsqb;\underset{i=1}{\overset{n}{\&Sigma;}}{a}_i\left(t\right)\exp (j\&Integral;{\mathrm{\&omega;}}_i(t\left)dt\right)\&rsqb;$

Step F-5：Hilbert marginal spectrum is sought, the spectral peak of signal harmonic component is extracted；

The expression formula of Hilbert marginal spectrum is：

$h\left(\mathrm{\&omega;}\right)=\underset{0}{\overset{T}{\&Integral;}}H(\mathrm{\&omega;},t)dt$

Wherein, T represents the sampling time, and its computing formula is：T=N/f_s。