CN106546818B - A kind of harmonic signal detection method based on differential nonlinearity Mode Decomposition - Google Patents
A kind of harmonic signal detection method based on differential nonlinearity Mode Decomposition Download PDFInfo
- Publication number
- CN106546818B CN106546818B CN201610916571.9A CN201610916571A CN106546818B CN 106546818 B CN106546818 B CN 106546818B CN 201610916571 A CN201610916571 A CN 201610916571A CN 106546818 B CN106546818 B CN 106546818B
- Authority
- CN
- China
- Prior art keywords
- signal
- nonlinear model
- harmonic
- component
- follows
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
- G01R23/00—Arrangements for measuring frequencies; Arrangements for analysing frequency spectra
- G01R23/16—Spectrum analysis; Fourier analysis
Abstract
The invention discloses a kind of harmonic signal detection methods based on differential nonlinearity Mode Decomposition, difference algorithm is first applied to original signal by this new harmonic detecting method, then obtains a series of significant nonlinear model components using nonlinear model decomposition algorithm；Nonlinear model decomposition is carried out again after integrating to each nonlinear model component, and extracts nonlinear model component of first nonlinear model component as original signal；Finally, extracting the various harmonic signals of signal with Hilbert marginal spectrum.Harmonic signal detection method based on differential nonlinearity Mode Decomposition inherits the noise robustness characteristic of nonlinear model decomposition method, and is able to suppress the interference of chaos and noise signal, advantageous in terms of extracting higher harmonic components by a small margin.
Description
Technical field
The present invention relates to Power estimations and signal processing applications field, more particularly to one kind to be based on differential nonlinearity Mode Decomposition
Harmonic signal detection method.
Background technique
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
In the clutter etc. that ocean surface generates.Signal is detected and analyzed with certain difficulty under Chaotic Background.
Fourier transformation has the advantages that simple and calculating is efficient, is the analysis most common algorithm of harmonic wave.But it exists
Three major defects: aliasing, fence effect and spectrum leakage are not suitable for analysis non-stationary signal.Based on Radial Basis Function neural
The method of network can also be used for the harmonic amplitude of detection 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 transformation to each IMF, are closed
In the temporal frequency Joint Distribution of original signal.EMD is widely used in non-linear and Non-stationary Signal Analysis.But in reality
EMD will appear mode mixing phenomenon in, show two aspects: first is that it is poor to contain scale in the same IMF component
Different biggish signal component；Second is that the signal component of the same scale has appeared in different IMF.It is asked to solve this
Topic plus white Gaussian noise and applies EMD to signal, then, gathers empirical mode decomposition (Ensemble Empirical
Mode Decomposition, EEMD) it comes into being.Difference empirical mode decomposition (Differential Empirical Mode
Decomposition, DEMD) first then EMD is applied to original signal calculus of differences, EMD can be extracted and EEMD can not
The higher harmonic components by a small margin of separation.However, three of the above method all haves the defects that noise-sensitive.
Signal decomposition can be nonlinear model significant on a series of physical point by nonlinear model decomposition method (NMD)
Amount, while eliminating noise.Relative to past method, NMD has fabulous make an uproar because parameter can be adaptive selected
Sound robust property.But performance of the NMD when harmonic signal by a small margin detects not is very well, especially to deposit in Chaotic Background
In case, performance is general.
Summary of the invention
The purpose of the present invention is to provide a kind of harmonic signal detection method based on differential nonlinearity 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 sensitive to chaos and noise signal and nonlinear model decomposition method performance when small signal harmonic detects
Bad disadvantage.
To achieve the above object, the technical solution adopted by the present invention are as follows: a kind of based on the humorous of differential nonlinearity Mode Decomposition
Wave signal detecting method, comprising the following steps:
Step A: prepare the original signal s (t) of processing to be detected, sample rate fs, data length N；
Step B: calculus of differences is carried out to original 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 s ' (t) is calculateds′(ω, t), wavelet transformation is defined as:
Wherein,It is the Fourier transform of s ' (t)；s′+(t) be signal s ' (t) positive frequency part, expression formula
Are as follows:ψ (t) is the wavelet function of wavelet transformation,For the Fourier transform of ψ (t), meet conditionSubscript * indicates conjugate operation；It is the crest frequency of small echo；Wavelet function is using logarithm just
State is distributed small echo,And ωψIt is expressed as follows:
Wherein, f0It is the resolution parameter for weighing time and frequency resolution in conversion process, usually default f0=1；
Step C-2: checking whether wavelet transformation is best time-frequency representation, if it is not, then using adding window Fourier transformation Gs′
(ω, t), adding window Fourier transform definition are as follows:
Wherein, g (t) is the window function of adding window Fourier transformation,For the Fourier transform of g (t), meet condition:Select Gaussian window as the window function of adding window Fourier transformation, expression formula are as follows:
Step C-3: ridge curve all in the time-frequency representation of signal s ' (t) is found outHere it definesIt is h
The ridge curve of subharmonic；
In sometime tn, using following formula algorithm, h maximum point can be found out；
N=1,2 in above formula ..., N；N is data length；Hs′(ω, t) is by above-mentioned wavelet transformation or adding window Fourier
Transformed time-frequency representation, as Ws′(ω, t) or Gs′(ω,t)；ω-(tn) and ω+(tn) it is respectively tnMoment maximum value is searched
The lower bound of rope frequency range and the upper bound；
By Hs′The ridge point line at all moment found out in (ω, t) may make up h vallate curve
Step C-4: ridge curve Reconstruction h order harmonic components are utilizedWherein A(h)(t)、WithIt is its amplitude, phase and frequency respectively；Calculating process is as follows:
If the time-frequency representation of signal s ' (t) uses wavelet transformation, h order harmonic components x(h)(t) it is obtained by formula (1):
If the time-frequency representation of signal s ' (t) uses adding window Fourier transformation, h order harmonic components x(h)(t) by formula
(2) it obtains:
Wherein,WithRespectively wavelet transformation and adding window Fourier transformation is as produced by parabola interpolation
Discretization influence amendment；
Step C-5: effective harmonic component is determined with the noise immunity substitution method of inspection；Noise immunity substitutes the method for inspection and utilizes
Alternate data identifies the true and false of the harmonic component extracted, filters out all true harmonic components, and work as continuous three
Harmonic component is judged as fictitious time and stops decomposable process；Specific step is as follows:
(1) the identification statistic D of a certain harmonic component extracted is calculated0(αA,αν)；
The amplitude A of each harmonic component extracted(h)(t) and frequency ν(h)(t) the degree of order can compose entropy with itWithQuantitatively to measure, whereinWithRespectively A(h)(t) and ν(h)(t) Fourier
Leaf transformation, identification statistic D are defined as follows:
Wherein, αAAnd αvRespectivelyWithWeight coefficient；
(2) N is created for signal s ' (t)sA Fourier transform alternate data, production method are as follows:
Wherein, φξObey [0,2 π) on be uniformly distributed, each φξA corresponding Fourier transform alternate data；
(3) time-frequency representation corresponding with each alternate data is calculated, and therefrom extracts each harmonic component respectively, is counted
Calculate the identification statistic of each alternate dataHere significance index is defined are as follows:
In formula,To meet Ds> D0Alternate data number；Assuming that creation NsA alternate data and will be significant
Property horizontal setup measures be p, i.e., at least Ns× p alternate data meets Ds> D0Just think that the component is not noise, thus after
Continuous decomposable process；Noise immunity substitution examines the parameter (α for using three groups of different valuesA,αν), that is, it calculates separately out D (1,1), D (0,
1) and the value of D (1,0), as long as wherein at least one value does not meet null hypothesis, then it is assumed that meet Ds> D0；
(4) the comprehensive measurement value of the degree of correlation between harmonic wave is calculated
Wherein,
In formula, ah=A(h)(t)/A(1)It (t) is constant, wA,wφ,wvIt representsWeight；Default uses ρ(h)
≡ρ(h)(1,1,0) equal weight is distributed for amplitude and phase consistency, and does not distribute weight to frequency invariance；
(5) in order to reduce the false judgment to true harmonic components, the threshold value of comprehensive measurement value is defined are as follows:
(6) when the comprehensive measurement value index of a harmonic component meets ρ(h)≥ρminAnd significance index
When significance-level >=p=95%, then it is assumed that this harmonic component is true harmonic component by examining；If
It cannot be substituted and be examined by noise immunity, then stop nonlinear model decomposition；
Step C-6: all true harmonic components are added and constitute a nonlinear model component c '1(t)；
Step C-7: subtracting the nonlinear model component from signal s ' (t), repeats step C-1 to step C-6, obtains institute
Some nonlinear model component c 'i(t)；
Step D: to obtained nonlinear model component c 'i(t) integral obtains bi(t)；
Step E: to bi(t) nonlinear model decomposition is carried out with nonlinear model decomposition method again；
Step E-1: by the process of step C-1 to step C-6 to each bi(t) nonlinear model decomposition is carried out；
Step E-2: each b is extractedi(t) first nonlinear model component is as the non-linear of original signal s (t)
Mode component ci(t)；Finally obtain the nonlinear model decomposition result of original signal:
Step F: spectrum analysis is carried out to obtained s (t), extracts harmonic signal；
Step F-1: to ci(t) it does Hilbert transform and generates quadrature component
Step F-2: construction complex signal zi(t), expression formula are as follows:
Step F-3: by complex signal zi(t) it is converted into polar form, acquires instantaneous envelope ai(t) and instantaneous frequency ωi
(t)；zi(t) polar form are as follows:
Wherein, instantaneous envelope ai(t) and instantaneous phase φi(t) it is expressed as follows:
Its instantaneous frequency ωiIt (t) is instantaneous phase φi(t) it to the derivative of time, indicates are as follows:
Step F-4: all nonlinear model component c are soughti(t) hilbert spectrum H (ω, t)；
Hilbert spectrum is defined as:
Wherein, M is the number of nonlinear model component；
Step F-5: seeking Hilbert marginal spectrum, extracts the spectral peak of signal harmonic component；
The expression formula of Hilbert marginal spectrum are as follows:
Wherein, T indicates the sampling time, its calculation formula is: T=N/fs。
Compared with prior art, the present invention its remarkable advantage is:
1. differential nonlinearity Mode Decomposition method can more effectively inhibit interference caused by noise and chaotic signal；
2. differential nonlinearity Mode Decomposition method can extract wherein amplitude for the signal with multiple harmonic components
Lesser harmonic component.
Detailed description of the invention
Fig. 1 is the time domain waveform of chaotic signal d (t) and original signal s (t)；
The Hilbert marginal spectrum that Fig. 2 is obtained when being SNR=25dB using differential nonlinearity Mode Decomposition method；
The Hilbert marginal spectrum that Fig. 3 is obtained when being SNR=25dB using nonlinear model decomposition method；
The Hilbert marginal spectrum that Fig. 4 is obtained when being SNR=25dB using empirical mode decomposition method；
The Hilbert marginal spectrum obtained when Fig. 5 is SNR=25dB using set empirical mode decomposition method；
The Hilbert marginal spectrum that Fig. 6 is obtained when being SNR=25dB using difference empirical mode decomposition method.
Specific embodiment
Below by taking the original signal added with white Gaussian noise and Chaotic Background as an example, in conjunction with attached drawing, the present invention will be described in detail
Embodiment.
The present invention carries out calculus of differences to original signal s (t) first and obtains new signal, then carries out nonlinear model to it
Formula decomposes to obtain one group of nonlinear model component c 'i(t)；Obtained nonlinear model component is integrated to obtain bi(t)；Again
Nonlinear model decomposition is carried out to it, extracts bi(t) first nonlinear model component is as the non-linear of original signal
Mode component is denoted as ci(t)；Finally, using spectral analysis technology the field of signal processing the advantages of, to finally obtained non-linear
Mode component carries out Hilbert limit spectrum analysis, realizes the detection of harmonic signal.
Assuming that original signal s (t) are as follows:
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 fs=600Hz, data length N=650.
X (t) is generated by typical Lorenz chaos system, and system model is described as follows:
Wherein, σ=10, r=28, b=8/3, initial value x0=y0=z0=0.1.D (t) and s (t) time domain waveform such as Fig. 1
It is shown.
The present invention is that a kind of based on DNMD, (Differential Nonlinear Mode Decomposition, difference are non-
Linear model decompose) harmonic signal detection method, include the following steps:
Step A: prepare the original signal s (t) of processing to be detected, sample rate fs=600Hz, data length N=
650；
Step B: calculus of differences is carried out to original 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 s ' (t) is calculateds′(ω, t), wavelet transformation is defined as:
Wherein,It is the Fourier transform of s ' (t)；s′+(t) be signal s ' (t) positive frequency part, expression formula
Are as follows:ψ (t) is the wavelet function of wavelet transformation,For the Fourier transform of ψ (t), meet conditionSubscript * indicates conjugate operation；It is the crest frequency of small echo；Wavelet function is using logarithm just
State is distributed small echo,And ωψIt is expressed as follows:
Wherein, f0It is the resolution parameter for weighing time and frequency resolution in conversion process, usually default f0=1；
Step C-2: checking whether wavelet transformation is best time-frequency representation, if it is not, then using adding window Fourier transformation Gs′
(ω, t), adding window Fourier transform definition are as follows:
Wherein, g (t) is the window function of adding window Fourier transformation,For the Fourier transform of g (t), meet condition:Select Gaussian window as the window function of adding window Fourier transformation, expression formula are as follows:
Step C-3: ridge curve all in the time-frequency representation of signal s ' (t) is found outWhen so-called ridge curve is exactly
Frequency schemes the curve that upper some Local modulus maximas are formed by connecting, and each maximum point is just ridge point；Here it definesIt is h
The ridge curve of subharmonic；
In sometime tn, using following formula algorithm, h maximum point can be found out；
N=1,2 in above formula ..., N；N is data length；Hs′(ω, t) is by above-mentioned wavelet transformation or adding window Fourier
Transformed time-frequency representation, as Ws′(ω, t) or Gs′(ω,t)；ω_(tn) and ω+(tn) it is respectively rnMoment maximum value is searched
The lower bound of rope frequency range and the upper bound；
By Hs′The ridge point line at all moment found out in (ω, t) may be constructed h vallate curve
Step C-4: ridge curve Reconstruction h order harmonic components are utilizedWherein A(h)(t)、WithIt is its amplitude, phase and frequency respectively；Calculating process is as follows:
If the time-frequency representation of signal s ' (t) uses wavelet transformation, h order harmonic components x(h)(t) it is obtained by formula (1):
If the time-frequency representation of signal s ' (t) uses adding window Fourier transformation, h order harmonic components x(h)(t) by formula
(2) it obtains:
Wherein,WithRespectively wavelet transformation and adding window Fourier transformation is as produced by parabola interpolation
Discretization influence amendment；
Step C-5: effective harmonic component is determined with the noise immunity substitution method of inspection；Noise immunity substitutes the method for inspection and utilizes
Alternate data identifies the true and false of the harmonic component extracted, filters out all true harmonic components, and work as continuous three
Harmonic component is judged as fictitious time and stops decomposable process；Specific step is as follows:
(1) the identification statistic D of a certain harmonic component extracted is calculated0(αA,αν)；
The amplitude A of each harmonic component extracted(h)(t) and frequency ν(h)(t) the degree of order can compose entropy with itWithQuantitatively to measure, whereinWithRespectively A(h)(t) and ν(h)(t) Fourier
Leaf transformation, identification statistic D are defined as follows:
Wherein, αAAnd αvRespectivelyWithWeight coefficient；
(3) N is created for signal s ' (t)sA Fourier transform alternate data, production method are as follows:
Wherein, φξObey [0,2 π) on be uniformly distributed, each φξA corresponding Fourier transform alternate data；
(3) time-frequency representation corresponding with each alternate data is calculated, and therefrom extracts each harmonic component respectively, is counted
Calculate the identification statistic of each alternate dataHere significance index is defined are as follows:
In formula,To meet Ds> D0Alternate data number；Assuming that creation NsA alternate data and will be significant
Property horizontal setup measures be p, i.e., at least Ns× p alternate data meets Ds> D0Just think that the component is not noise, thus after
Continuous decomposable process；Noise immunity substitution examines the parameter (α for using three groups of different valuesA,αν), that is, it calculates separately out D (1,1), D (0,1)
With the value of D (1,0), as long as wherein at least one value does not meet null hypothesis, then it is assumed that meet Ds> D0；
(4) the comprehensive measurement value of the degree of correlation between harmonic wave is calculated
Wherein,
In formula, ah=A(h)(t)/A(1)It (t) is constant；wA,wφ,wvIt representsWeight；Herein, default
Use ρ(h)≡ρ(h)(1,1,0) equal weight is distributed for amplitude and phase consistency, and does not distribute weight to frequency invariance；
(5) in order to reduce the false judgment to true harmonic components, the threshold value of comprehensive measurement value is defined are as follows:
(6) when the comprehensive measurement value index of a harmonic component meets ρ(h)≥ρminAnd significance index
When significance-level >=p=95%, then it is assumed that this harmonic component is true harmonic component by examining；If
It cannot be substituted and be examined by noise immunity, then stop nonlinear model decomposition；
Step C-6: all true harmonic components are added and constitute a nonlinear model component c '1(t)；
Step C-7: subtracting the nonlinear model component from signal s ' (t), repeats step C-1 to step C-6, obtains institute
Some nonlinear model component c 'i(t)；
Step D: to obtained nonlinear model component c 'i(t) integral obtains bi(t)；
Step E: to bi(t) nonlinear model decomposition is carried out with nonlinear model decomposition method again；
Step E-1: by the process of step C-1 to step C-6 to each bi(t) nonlinear model decomposition is carried out；
Step E-2: each b is extractedi(t) first nonlinear model component is as the non-linear of original signal s (t)
Mode component ci(t)；Finally obtain the nonlinear model decomposition result of original signal:
Step F: spectrum analysis is carried out to obtained s (t), extracts harmonic signal；
Step F-1: to ci(t) it does Hilbert transform and generates quadrature component
Step F-2: construction complex signal zi(t), expression formula are as follows:
Step F-3: by complex signal zi(t) it is converted into polar form, acquires instantaneous envelope ai(t) and instantaneous frequency ωi
(t)；zi(t) polar form are as follows:
Wherein, instantaneous envelope ai(t) and instantaneous phase φi(t) it is expressed as follows:
Its instantaneous frequency ωiIt (t) is instantaneous phase φi(t) it to the derivative of time, indicates are as follows:
Step F-4: all nonlinear model component c are soughti(t) hilbert spectrum H (ω, t)；
Hilbert spectrum is defined as:
Wherein, M is the number of nonlinear model component；
Step F-5: seeking Hilbert marginal spectrum, extracts the spectral peak of signal harmonic component；
The expression formula of Hilbert marginal spectrum are as follows:
Wherein, T indicates the sampling time, its calculation formula is: T=N/fs。
Fig. 2 to Fig. 6 is followed successively by as Signal to Noise Ratio (SNR)=25dB, and differential nonlinearity 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 obtain
The Hilbert marginal spectrum arrived.From figure 2 it can be seen that differential nonlinearity Mode Decomposition method can clearly extract it is original
The fundamental frequency (30Hz) of signal and three times with 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
It measures (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, can not accurately extract other harmonic components of signal, only at third harmonic frequencies (90Hz)
Nearby there is a degree of protrusion.
For the performance for quantitatively comparing difference nonlinear model decomposition method Yu other methods, parameter Rk, it is defined asWherein, k is the number of harmonic wave, AkIt is the amplitude of corresponding harmonic frequency,It is to be decomposed point
The average value of the spectrum of amount.R of the every kind of method under different signal-to-noise ratiokValue is as shown in table 1.
R of the every kind of method of table 1 under different signal-to-noise ratiokValue
DNMD is differential nonlinearity Mode Decomposition method in table 1, and NMD is nonlinear model decomposition method, and EMD is Empirical Mode
State decomposition method, EEMD set empirical mode decomposition method and DEMD are difference empirical mode decomposition method, and SNR is signal-to-noise ratio.
From 1 column data of table it can be seen that as SNR=20dB, differential nonlinearity Mode Decomposition method can be extracted
The fundamental frequency (30Hz) of original signal, three times with quintuple harmonics frequency (90Hz, 150Hz), nonlinear model decomposition method is only
Fundametal compoment (30Hz) and third-harmonic component (90Hz) have been extracted, and empirical mode decomposition method, set empirical modal divide
Solution method and difference empirical mode decomposition method have only extracted fundametal compoment (30Hz).Although difference is non-in SNR=15dB
Linear model decomposition method is lost the ability for extracting quintuple harmonics (150Hz), but its R1=52.74dB and R3=28.30dB
It is all larger than the R of nonlinear model decomposition method1=51.90dB and R3=26.28dB.
In conclusion the above embodiments are merely illustrative of the technical solutions of the present invention, it is not intended to limit guarantor of the invention
Protect range.All within the spirits and principles of the present invention, any modification, equivalent substitution, improvement and etc. done should all cover
In scope of the presently claimed invention.
Claims (1)
1. a kind of harmonic signal detection method based on differential nonlinearity Mode Decomposition, it is characterised in that: the following steps are included:
Step A: prepare the original signal s (t) of processing to be detected, sample rate fs, data length N；
Step B: calculus of differences is carried out to original 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 s ' (t) is calculateds′(ω, t), wavelet transformation is defined as:
Wherein,It is the Fourier transform of s ' (t)；s′+(t) be signal s ' (t) positive frequency part, expression formula are as follows:ψ (t) is the wavelet function of wavelet transformation,For the Fourier transform of ψ (t), meet conditionSubscript * indicates conjugate operation；It is the crest frequency of small echo；Wavelet function is using logarithm just
State is distributed small echo,And ωψIt is expressed as follows:
Wherein, f0It is the resolution parameter for weighing time and frequency resolution in conversion process, usually default f0=1；
Step C-2: checking whether wavelet transformation is best time-frequency representation, if it is not, then using adding window Fourier transformation Gs′(ω,
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), meet condition:Select Gaussian window as the window function of adding window Fourier transformation, expression formula are as follows:
Step C-3: ridge curve all in the time-frequency representation of signal s ' (t) is found outDefinitionIt is h subharmonic
Ridge curve；
In sometime tn, using following formula algorithm, find out h maximum point；
N=1,2 in above formula ..., N；N is data length；Hs′(ω, t) is by above-mentioned wavelet transformation or adding window Fourier transformation
Time-frequency representation afterwards, as Ws′(ω, t) or Gs′(ω,t)；
By Hs′The ridge point line at all moment found out in (ω, t) constitutes h vallate curve
Step C-4: ridge curve Reconstruction h order harmonic components are utilizedWherein A(h)(t)、WithIt is its amplitude, phase and frequency respectively；Calculating process is as follows:
If the time-frequency representation of signal s ' (t) uses wavelet transformation, h order harmonic components x(h)(t) it is obtained by formula (1):
If the time-frequency representation of signal s ' (t) uses adding window Fourier transformation, h order harmonic components x(h)(t) it is obtained by formula (2)
It arrives:
Wherein,WithRespectively wavelet transformation and adding window Fourier transformation as caused by parabola interpolation from
The amendment that dispersion influences；
Step C-5: effective harmonic component is determined with the noise immunity substitution method of inspection；Noise immunity substitutes the method for inspection using substitution
Data identify the true and false of the harmonic component extracted, filter out all true harmonic components, and work as continuous three harmonic waves
Component is judged as fictitious time and stops decomposable process；Specific step is as follows:
(1) the identification statistic D of a certain harmonic component extracted is calculated0(αA,αν)；
The amplitude A of each harmonic component extracted(h)(t) and frequency ν(h)(t) the degree of order composes entropy with itWithQuantitatively to measure, whereinWithRespectively A(h)(t) and ν(h)(t) Fourier transform, identification
Statistic D is defined as follows:
Wherein, αAAnd αvRespectivelyWithWeight coefficient；
(2) N is created for signal s ' (t)sA Fourier transform alternate data, production method are as follows:
Wherein, φξObey [0,2 π) on be uniformly distributed, each φξA corresponding Fourier transform alternate data；
(3) time-frequency representation corresponding with each alternate data is calculated, and therefrom extracts each harmonic component respectively, is calculated
The identification statistic of each alternate dataDefine significance index are as follows:
In formula,To meet Ds> D0Alternate data number；Assuming that creation NsA alternate data and by conspicuousness water
Flat setup measures are p, i.e., at least Ns× p alternate data meets Ds> D0Just think that the component is not noise, to continue point
Solution preocess；Noise immunity substitution examines the parameter (α for using three groups of different valuesA,αν), that is, calculate separately out D (1,1), D (0,1) and
The value of D (1,0), as long as wherein at least one value does not meet null hypothesis, then it is assumed that meet Ds> D0；
(4) the comprehensive measurement value of the degree of correlation between harmonic wave is calculated
Wherein,
In formula, ah=A(h)(t)/A(1)(t), wA,wφ,wvIt representsWeight；Default uses ρ(h)≡ρ(h)(1,1,
0) equal weight is distributed for amplitude and phase consistency, and does not distribute weight to frequency invariance；
(5) in order to reduce the false judgment to true harmonic components, the threshold value of comprehensive measurement value is defined are as follows:
(6) when the comprehensive measurement value index of a harmonic component meets ρ(h)≥ρminAnd significance index
When significance-level >=p=95%, then it is assumed that this harmonic component is true harmonic component by examining；If
It cannot be substituted and be examined by noise immunity, then stop nonlinear model decomposition；
Step C-6: all true harmonic components are added and constitute a nonlinear model component c '1(t)；
Step C-7: subtracting the nonlinear model component from signal s ' (t), repeats step C-1 to step C-6, obtains all
Nonlinear model component c 'i(t)；
Step D: to obtained nonlinear model component c 'i(t) integral obtains bi(t)；
Step E: to bi(t) nonlinear model decomposition is carried out with nonlinear model decomposition method again；
Step E-1: by the process of step C-1 to step C-6 to each bi(t) nonlinear model decomposition is carried out；
Step E-2: each b is extractedi(t) nonlinear model of first nonlinear model component as original signal s (t)
Component ci(t)；Finally obtain the nonlinear model decomposition result of original signal:
Step F: spectrum analysis is carried out to obtained original signal s (t), extracts harmonic signal；
Step F-1: to nonlinear model component ci(t) it does Hilbert transform and generates quadrature component
Step F-2: construction complex signal zi(t), expression formula are as follows:
Step F-3: by complex signal zi(t) it is converted into polar form, acquires instantaneous envelope ai(t) and instantaneous frequency ωi(t)；zi
(t) polar form are as follows:
Wherein, instantaneous envelope ai(t) and instantaneous phase φi(t) it is expressed as follows:
Its instantaneous frequency ωiIt (t) is instantaneous phase φi(t) it to the derivative of time, indicates are as follows:
Step F-4: all nonlinear model component c are soughti(t) hilbert spectrum H (ω, t)；
Hilbert spectrum is defined as:
Step F-5: seeking Hilbert marginal spectrum, extracts the spectral peak of signal harmonic component；
The expression formula of Hilbert marginal spectrum are as follows:
Wherein, T indicates the sampling time, its calculation formula is: T=N/fs。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610916571.9A CN106546818B (en) | 2016-10-20 | 2016-10-20 | A kind of harmonic signal detection method based on differential nonlinearity Mode Decomposition |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610916571.9A CN106546818B (en) | 2016-10-20 | 2016-10-20 | A kind of harmonic signal detection method based on differential nonlinearity Mode Decomposition |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN106546818A CN106546818A (en) | 2017-03-29 |
| CN106546818B true CN106546818B (en) | 2019-04-09 |
Family
ID=58392002
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201610916571.9A Expired - Fee Related CN106546818B (en) | 2016-10-20 | 2016-10-20 | A kind of harmonic signal detection method based on differential nonlinearity Mode Decomposition |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN106546818B (en) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111397645A (en) * | 2020-04-06 | 2020-07-10 | 华中科技大学 | Phase difference decomposition and adjustment method and system |
Families Citing this family (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107219840A (en) * | 2017-05-05 | 2017-09-29 | 浙江理工大学 | Towards the regulating valve nonlinear characteristic detection method and system of Natural Gas Station |
| CN107229597A (en) * | 2017-05-31 | 2017-10-03 | 成都理工大学 | Synchronous extruding generalized S-transform signal Time-frequency Decomposition and reconstructing method |
| CN107403139B (en) * | 2017-07-01 | 2021-05-25 | 南京理工大学 | Urban rail train wheel flat scar fault detection method |
| CN107687941A (en) * | 2017-07-03 | 2018-02-13 | 昆明理工大学 | A kind of high-pressure diaphragm pump check valve Incipient Fault Diagnosis method based on analysis of vibration signal |
| CN108345289B (en) * | 2018-01-08 | 2020-03-31 | 浙江大学 | Industrial process stability detection method based on alternative data method |
| CN110363141B (en) * | 2019-07-15 | 2021-09-17 | 郑州大学 | Method for diagnosing a fault in a gas pressure regulator |
| CN113341226B (en) * | 2021-06-21 | 2022-04-29 | 合肥美的暖通设备有限公司 | Harmonic detection method, device, frequency converter and storage medium |
| CN117268483B (en) * | 2023-11-23 | 2024-02-23 | 青岛鼎信通讯科技有限公司 | Instantaneous flow metering method suitable for ultrasonic water meter |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6901353B1 (en) * | 2003-07-08 | 2005-05-31 | The United States Of America As Represented By The Administrator Of The National Aeronautics And Space Administration | Computing Instantaneous Frequency by normalizing Hilbert Transform |
| CN103163372A (en) * | 2013-03-26 | 2013-06-19 | 山西省电力公司长治供电分公司 | Method for analyzing harmonics of power system by adopting Hilbert-Huang transform (HHT) |
| CA2783089A1 (en) * | 2012-07-11 | 2014-01-11 | Farid Taheri | Damage detection in pipes and joint systems |
| CN103901273A (en) * | 2012-12-28 | 2014-07-02 | 白晓民 | Power harmonic detection method and power harmonic detection device |
| CN103941091A (en) * | 2014-04-25 | 2014-07-23 | 福州大学 | Power system HHT harmonious wave detection method based on improved EMD end point effect |
| CN105510711A (en) * | 2015-12-24 | 2016-04-20 | 合肥工业大学 | Empirical mode decomposition-based improved harmonic analysis method |
-
2016
- 2016-10-20 CN CN201610916571.9A patent/CN106546818B/en not_active Expired - Fee Related
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6901353B1 (en) * | 2003-07-08 | 2005-05-31 | The United States Of America As Represented By The Administrator Of The National Aeronautics And Space Administration | Computing Instantaneous Frequency by normalizing Hilbert Transform |
| CA2783089A1 (en) * | 2012-07-11 | 2014-01-11 | Farid Taheri | Damage detection in pipes and joint systems |
| CN103901273A (en) * | 2012-12-28 | 2014-07-02 | 白晓民 | Power harmonic detection method and power harmonic detection device |
| CN103163372A (en) * | 2013-03-26 | 2013-06-19 | 山西省电力公司长治供电分公司 | Method for analyzing harmonics of power system by adopting Hilbert-Huang transform (HHT) |
| CN103941091A (en) * | 2014-04-25 | 2014-07-23 | 福州大学 | Power system HHT harmonious wave detection method based on improved EMD end point effect |
| CN105510711A (en) * | 2015-12-24 | 2016-04-20 | 合肥工业大学 | Empirical mode decomposition-based improved harmonic analysis method |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111397645A (en) * | 2020-04-06 | 2020-07-10 | 华中科技大学 | Phase difference decomposition and adjustment method and system |
| CN111397645B (en) * | 2020-04-06 | 2020-12-18 | 华中科技大学 | Phase difference decomposition and adjustment method and system |
Also Published As
| Publication number | Publication date |
|---|---|
| CN106546818A (en) | 2017-03-29 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106546818B (en) | A kind of harmonic signal detection method based on differential nonlinearity Mode Decomposition | |
| CN112697887B (en) | Ultrasonic detection defect qualitative identification method based on neural network | |
| CN107340055B (en) | It is a kind of based on the random resonant weak signal detection method for more estimating fusion | |
| CN105486938B (en) | A kind of substation's mixed noise separation method | |
| CN108663605A (en) | Local discharge signal detection method based on Coupled Duffing oscillators | |
| CN108548957A (en) | The double-spectrum analysis method being combined based on circular modulating frequency spectrum and segmentation cross-correlation | |
| CN106771598B (en) | A kind of Adaptive spectra kurtosis signal processing method | |
| CN109001703A (en) | A kind of sea clutter denoising method based on the processing of wavelet packet multi-threshold | |
| CN109444519A (en) | Substation's noise source separation method towards complicated acoustic environment | |
| CN103905656B (en) | The detection method of residual echo and device | |
| Xu et al. | An adaptive spectrum segmentation method to optimize empirical wavelet transform for rolling bearings fault diagnosis | |
| Hua et al. | The methodology of modified frequency band envelope kurtosis for bearing fault diagnosis | |
| Ou et al. | Compound fault diagnosis of gearboxes based on GFT component extraction | |
| CN108761202A (en) | The harmonic detecting method that pole symmetric mode decomposition and Hilbert transform are combined | |
| CN106645947A (en) | Time-frequency analysis method based on nonlinear mode decomposition and adaptive optimal kernel | |
| Veeraiyan et al. | Frequency domain based approach for denoising of underwater acoustic signal using EMD | |
| Alkishriwo et al. | Signal separation in the Wigner distribution domain using fractional Fourier transform | |
| CN109460614A (en) | Signal time based on instant bandwidth-frequency decomposition method | |
| Dai et al. | Application of wavelet denoising and time-frequency domain feature extraction on data processing of modulated signals | |
| Ce et al. | An improved wavelet threshold de-noising method and its application | |
| Zheng et al. | Zero-Phase Filter-Based Adaptive Fourier Decomposition and Its Application to Fault Diagnosis of Rolling Bearing | |
| Yonghong et al. | Detection of small leakage from pipeline based on improved harmonic wavelet | |
| Xu et al. | Research on partial discharge de-noising for transformer based on synchro-squeezed continuous wavelet transform | |
| CN105044457B (en) | A kind of antimierophonic trend of harmonic detection method of power | |
| CN109632945A (en) | A kind of noise-reduction method suitable for Pulsed eddy current testing signal |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20190409 Termination date: 20201020 |
|
| CF01 | Termination of patent right due to non-payment of annual fee |