[go: up one dir, main page]

CN107748734B - Analytic-empirical mode decomposition method - Google Patents

Analytic-empirical mode decomposition method Download PDF

Info

Publication number
CN107748734B
CN107748734B CN201711049787.0A CN201711049787A CN107748734B CN 107748734 B CN107748734 B CN 107748734B CN 201711049787 A CN201711049787 A CN 201711049787A CN 107748734 B CN107748734 B CN 107748734B
Authority
CN
China
Prior art keywords
imf
mode decomposition
decomposition
signal
mode
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
Application number
CN201711049787.0A
Other languages
Chinese (zh)
Other versions
CN107748734A (en
Inventor
陈中柘
毛自荐
颜小刚
朱圆
王霄
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
University of Electronic Science and Technology of China
Original Assignee
University of Electronic Science and Technology of China
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by University of Electronic Science and Technology of China filed Critical University of Electronic Science and Technology of China
Priority to CN201711049787.0A priority Critical patent/CN107748734B/en
Publication of CN107748734A publication Critical patent/CN107748734A/en
Application granted granted Critical
Publication of CN107748734B publication Critical patent/CN107748734B/en
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
    • G06F17/141—Discrete Fourier transforms
    • G06F17/142—Fast Fourier transforms, e.g. using a Cooley-Tukey type algorithm
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Mathematical Physics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • Discrete Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • Other Investigation Or Analysis Of Materials By Electrical Means (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)

Abstract

The invention discloses an analytic-empirical mode decomposition method, which comprises the following steps: s1, screening the original signals through an empirical mode decomposition method to obtain n IMF components and 1 residual amount; and S2, processing the IMF components through an analytic mode decomposition algorithm, and then reconstructing the mode according to a certain rule. The invention provides a screening termination criterion based on an effective data segment aiming at the problem of the screening termination criterion of empirical mode decomposition, provides an analysis-empirical mode decomposition method aiming at the mode confusion phenomenon, carries out analysis mode decomposition on the modes with mode confusion, and finally reconstructs the modes. The analytic-empirical mode decomposition method can effectively inhibit mode confusion, and has short decomposition time and higher method execution efficiency.

Description

Analytic-empirical mode decomposition method
Technical Field
The invention relates to an analytic-empirical mode decomposition method.
Background
Empirical Mode Decomposition (EMD) is a self-adaptive signal processing method which is widely applied, and aiming at the problem of screening termination criteria in the implementation process of empirical mode decomposition, the screening termination criteria based on an effective data segment is provided and compared with the traditional screening termination criteria, the new screening criteria avoid the influence of end point effect, and the accuracy of empirical mode decomposition is improved. Aiming at the mode confusion phenomenon, an analytic-empirical mode decomposition method is provided, the modes which have mode confusion in the empirical mode decomposition are subjected to analytic mode decomposition, and the modes are reconstructed, so that the mode confusion of each order of modes is obviously smaller than the empirical mode decomposition. Compared with the existing main improvement method of empirical mode decomposition, namely the ensemble average empirical mode decomposition method, the method provided by the invention can effectively inhibit mode aliasing as the ensemble average empirical mode decomposition method. However, the method has more obvious effect of inhibiting mode confusion and has higher execution efficiency.
The Empirical Mode Decomposition (EMD) method is based on the local characteristic time scale of a signal, and can decompose a complex signal into the sum of a limited number of Intrinsic Mode functions (IMF for short), each Intrinsic Mode Function represents an Intrinsic vibration Mode of an original signal, and the stable Intrinsic Mode functions are analyzed to extract the characteristic information of the signal.
Empirical mode decomposition has attracted considerable attention from the advent of numerous scholars. The method is widely applied to the fields of economics, biomedicine, engineering science and the like, but has some defects, mainly including Mode aliasing (Mode Mixing), End Effects (End Effects) and Stop conditions (Stop criterion) of Intrinsic Mode Function (IMF) criterion, overcladding and undercladding phenomena and the like. The empirical mode decomposition has the defects that the further popularization and the use of the empirical mode decomposition are restricted.
The condition for judging the intrinsic mode function in the empirical mode decomposition method is that the average value of the upper envelope line and the lower envelope line of the intrinsic mode function is zero, but in the factual decomposition process, the average value of the envelope lines of the intrinsic mode function components obtained by the empirical mode decomposition cannot be zero due to the interference of cubic spline curve interpolation, the endpoint effect, the influence of the adopted frequency and the like. Therefore, it is necessary to define a stopping criterion of the decomposition process feasible in engineering, i.e. the intrinsic mode function criterion problem of empirical mode decomposition is also called screening stopping criterion problem.
The Mode aliasing refers to that one IMF (intrinsic Mode function) includes a characteristic time scale with a great difference, or similar characteristic time scales are distributed in different IMFs, so that adjacent 2 IMF waveforms are aliased, mutually affected and difficult to recognize.
The mode confusion is an important defect of the EMD, once the mode confusion is generated, the physical meaning of the finally obtained intrinsic mode function component is not clear, the accuracy of signal decomposition is affected, and the application of the mode confusion in engineering is seriously hindered.
In order to solve the problem of mode confusion of empirical mode decomposition, scholars propose some improved methods. Most notably, Huang's proposed Ensemble Empirical Mode Decomposition (EEMD) is the best known method. The method utilizes the statistical characteristic that Gaussian white noise has uniform frequency distribution, so that signals added with the Gaussian white noise have continuity on different scales, and the problem of mode aliasing is solved on a certain scale. However, there are problems such as the number of times of adding white noise and the selection of the magnitude of the added white noise are highly subjective and the method is poorly adaptive. In addition, when the number and amplitude of the added white noise are reasonably selected, the mode confusion of low frequency can be artificially caused while the mode confusion of high frequency is suppressed. Moreover, the overall local mean decomposition method has a complex algorithm and a long program running time.
Disclosure of Invention
The invention aims to overcome the defects of the prior art and provide a method for resolving modes of modes with mode confusion and finally reconstructing the modes. The analytic-empirical mode decomposition method can effectively inhibit mode confusion, has short decomposition time and higher execution efficiency.
The purpose of the invention is realized by the following technical scheme: an analytical-empirical mode decomposition method, comprising the steps of:
s1, screening the original signals through an empirical mode decomposition method to obtain n IMF components and 1 residual amount;
and S2, processing the IMF components through an analytic mode decomposition algorithm, and reconstructing the mode.
Further, the step S1 includes the following sub-steps:
s11, determining all maximum points and minimum points of original signal x (t), fitting upper and lower envelope curves of the maximum points by cubic spline curve, ensuring the interval between the upper and lower envelope curves of signal x (t), calculating local mean value of the upper and lower envelope curves, and recording as m1；
S12, calculating the difference h between the signal x (t) and the local mean value1(t)：
h1(t)＝x(t)-m1 (1)；
S13, judgment h1(t) whether or not a condition for establishing the IMF can be satisfied; if not, h is1(t) repeating steps S11 and S12 as a new original signal until h1(t) satisfying a condition that IMF is established;
s14, order C1(t)＝h1(t)，C1(t) is the first IMF component obtained by decomposition;
s15, mixing C1(t) is separated from the signal x (t) and the resulting residual signal is noted: r is1(t)＝x(t)-C1(t)；
S16, mixing1(t) repeating steps S11-S15 as new original signal, and performing n finite cycles to obtain n IMF components and 1 residual rn(t); namely:
Figure GDA0003160195680000031
the original signal x (t) is:
Figure GDA0003160195680000032
wherein, Ci(t) denotes the i-th intrinsic mode function component, rn(t) is called residual function.
Further, in step S16, the cycle number n is calculated by using a screening termination criterion based on the valid data segment, and the specific operation method includes: when intrinsic mode function Cn(t) stopping the cycle when the following conditions are met:
Figure GDA0003160195680000033
in the formula (I), the compound is shown in the specification,
Figure GDA0003160195680000034
is represented by CnAnd (t) in the effective section of the upper and lower envelope mean curves, delta represents a deviation threshold value, and the value range is 0.01-0.1.
Further, the step S2 includes the following sub-steps:
s21, sequentially recording the n intrinsic mode function components obtained by calculation in the step S1 as IMF1(t)～IMFn(t)；
For IMF1(t) performing fast Fourier transform on the component to obtain an amplitude spectrum thereof, and determining a boundary frequency f from the amplitude spectrumb1(ii) a For IMF1(t) performing analytical mode decomposition from the boundary frequency fb1The decomposition yields two signals c1(t) and
Figure GDA0003160195680000035
c1(t) represents IMF1(ii) a correction component of (t),
Figure GDA0003160195680000036
is c1(t) a residual component; namely:
Figure GDA0003160195680000037
s22, mixing
Figure GDA0003160195680000038
Addition to IMF2(t) obtaining an updated IMF2(t) as
Figure GDA0003160195680000039
To pair
Figure GDA00031601956800000310
Performing fast Fourier transform to obtain
Figure GDA00031601956800000311
And determining the boundary frequency f from the amplitude spectrumb2(ii) a To pair
Figure GDA00031601956800000312
Performing analytic mode decomposition with the boundary frequency fb2To obtain c2(t) and
Figure GDA00031601956800000313
wherein, c2(t) is
Figure GDA00031601956800000314
The correction component of (a) is,
Figure GDA00031601956800000315
is c2(t) a residual component; namely:
Figure GDA00031601956800000316
s23, removing IMFn(t) repeating the operations of steps S21 and S22 for all intrinsic mode function components except for (t);
s24 for IMFn(t) performing the following operations: will be provided with
Figure GDA00031601956800000317
Addition to IMFnIn (t), obtaining IMF* n(t) for IMF* n(t) performing fast Fourier transform and plotting the amplitude spectrum thereof, determining the demarcation frequency fbk(ii) a For IMF* n(t) performing analytical mode decomposition from the boundary frequency fbkObtain a signal cn(t) and
Figure GDA0003160195680000041
namely:
Figure GDA0003160195680000042
will be provided with
Figure GDA0003160195680000043
Added to the residual function rn(t) obtaining the final residual error, denoted as vn(t)。
The invention has the beneficial effects that:
1. the invention provides an analysis-empirical mode decomposition method aiming at the mode confusion phenomenon existing in empirical mode decomposition, and the analysis mode decomposition is carried out on the modes with mode confusion, and finally the modes are reconstructed. The analytic-empirical mode decomposition method can effectively inhibit mode confusion, and has short decomposition time and higher method execution efficiency.
2. Aiming at the problem of the inherent modal function screening termination criterion of empirical mode decomposition, the invention provides the effective data segment-based screening termination criterion, compared with the traditional screening termination criterion, the result of empirical mode decomposition obtained by using the screening termination criterion of the invention is more accurate.
Drawings
FIG. 1 is a time domain diagram of a simulated signal x (t) and two component components used in a verification test (one);
FIG. 2 is a decomposition result diagram of a verification test (one) employing a cycle termination criterion based on valid data segments;
FIG. 3 is a graph of the decomposition results of the cycle termination criteria proposed by Huang in the validation test (one);
FIG. 4 is a time domain plot of a vibration signal for verifying a gear tooth break fault used in test (II);
FIG. 5 is a diagram of a first intrinsic mode function component of vibration signal analysis-empirical mode decomposition for a gear tooth breakage fault;
FIG. 6 is a time domain waveform of the simulated signal x (t) used in the verification test (III);
FIG. 7 is a graph of empirical mode decomposition results of simulation signals from a verification test (III);
FIG. 8 is a graph of the simulated signal ensemble-averaged empirical mode decomposition results of the verification test (III);
FIG. 9 is a graph of simulated signal analysis-empirical mode decomposition results for verification test (III);
FIG. 10 is a time domain plot of a rotor fault signal used in the validation test (IV);
FIG. 11 shows a partial mean decomposition of a rotor fault signal;
fig. 12 shows the result of the rotor fault signal analysis-empirical mode decomposition.
Detailed Description
The technical scheme of the invention is further explained by combining the attached drawings.
The invention discloses an analytic-empirical mode decomposition method, which comprises the following steps:
s1, screening the original signals through an empirical mode decomposition method to obtain n IMF components and 1 residual amount; empirical mode decomposition is a "sieving" process; screening out the component with minimum extreme value time characteristic scale, screening out the component with larger characteristic scale until the component with minimum extreme value time characteristic scale is screened outThe component with the largest characteristic dimension is screened out. That is, empirical mode decomposition yields an intrinsic mode function component C1，C2The average frequency of … … is gradually reduced. The method specifically comprises the following substeps:
s11, determining all maximum points and minimum points of original signal x (t), fitting upper and lower envelope curves of the maximum points by cubic spline curve, ensuring the interval between the upper and lower envelope curves of signal x (t), calculating local mean value of the upper and lower envelope curves, and recording as m1；
S12, calculating the difference h between the signal x (t) and the local mean value1(t)：
h1(t)＝x(t)-m1 (1)；
S13, judgment h1(t) whether or not a condition for establishing the IMF can be satisfied; if not, h is1(t) repeating steps S11 and S12 as a new original signal until h1(t) satisfying a condition that IMF is established;
s14, order C1(t)＝h1(t)，C1(t) is the first IMF component obtained by decomposition;
s15, mixing C1(t) is separated from the signal x (t) and the resulting residual signal is noted: r is1(t)＝x(t)-C1(t)；
S16, mixing1(t) repeating steps S11-S15 as new original signal, and performing n finite cycles to obtain n IMF components and 1 residual rn(t); namely:
Figure GDA0003160195680000051
the original signal x (t) is:
Figure GDA0003160195680000052
wherein, Ci(t) denotes the i-th intrinsic mode function component, rn(t) is called residual function.
The above-mentioned loop process of step S16 is when rn(t) is oneThe loop terminates when the function is monotonic, but in the actual decomposition process when r isnWhen the (t) satisfies the monotonic function condition, the number of cycles is often very large, and the number of intrinsic mode function components finally obtained by decomposition will be large. From a large number of experiments of empirical mode decomposition, it is found that most of the final components obtained by empirical mode decomposition are false intrinsic mode function components, and have no substantial physical significance.
The stopping criterion of the decomposition process proposed by Huang is achieved by limiting the size of the standard deviation between 2 adjacent IMF components, by an iteration threshold SdTo control the number of siftings. SdIs defined as:
Figure GDA0003160195680000061
in the formula: t is the time length of the signal; h isk(t)，hk-1(t) is the two adjacent IMF components in the decomposition process; when S isdAnd stopping the circulation when the range of 0.2-0.3 is satisfied.
It is very important to determine a reasonable iteration threshold, if the threshold is too small, the calculation amount will be greatly increased, and finally, the obtained intrinsic mode function component is meaningless. If the threshold is too large, it will be difficult to satisfy the decomposition result of the intrinsic mode function condition.
The screening termination criteria defined by Huang are in most cases more effective. It also presents some problems. As can be seen from formula (3), when h iskWhen (t) is 0, the resulting iteration threshold will be an indeterminate value if it is clearly inappropriate to use it as a screening criterion. In addition, the screening termination criterion does not fully consider the influence of the endpoint effect, and the final decomposition result has a large error after the distorted endpoint data is used. In many cases, the end-point effect of empirical mode decomposition is relatively significant, and even if some method is taken to suppress the end-point effect, the end-point effect is still present.
For SdIn the step S16, the loop is calculated by using the screening termination criteria based on the valid data segmentsThe number n, the valid segment, is the remaining segment excluding the two-end-point distorted segment. The influence of the endpoint effect on the screening criterion is fully considered in the criterion, and only the data of the effective data segment is used as the basis for calculating the iteration threshold, so that the influence of the empirical mode decomposition endpoint effect is avoided. The specific operation method in the step comprises the following steps: when intrinsic mode function Cn(t) stopping the cycle when the following conditions are met:
Figure GDA0003160195680000062
in the formula (I), the compound is shown in the specification,
Figure GDA0003160195680000063
is represented by CnAnd (t) in the effective section of the upper and lower envelope mean curves, delta represents a deviation threshold value, and the value range is 0.01-0.1.
Verification test (one); the application of the cycle termination criteria of the present invention in the empirical mode decomposition method is illustrated by a simulation signal and compared to the proposed criteria of n.e. huang. The expression of the simulation signal x (t) is:
x(t)＝x1(t)+x2(t),t∈[0,1] (5)
x1(t)＝(1+0.5sin(2π*4t))*2cos(2π*80t),x2(t)＝sin(2π*24t) (6)
the simulation signal x (t) is composed of x1(t) and x2(t) the two components. The simulated signal and the two component components are shown in figure 1. The results of empirical mode decomposition using the loop termination criteria proposed by the present invention are shown in fig. 2. The results of empirical mode decomposition using the termination criteria proposed by Huang are shown in FIG. 3. Comparing FIGS. 1, 2 and 3, it can be seen that the component C obtained by using the loop termination criteria of the present invention1Component imf derived from the termination criteria set forth in Huang1Are all relatively close to the ideal result x1(t)，C2And imf2Are all relatively close to the ideal result x2(t) which illustrates that empirical mode decomposition will fit this in the original signal regardless of which screening termination criterion is usedTwo components x1(t) and x2(t) accurately isolated. However, C1Curve ratio imf near the two end points of (A)1Is closer to the ideal curve x1(t)，C2Curve ratio imf near the two end points of (A)2Is closer to the ideal curve x2(t), which shows that the results are more accurate than those of empirical mode decomposition using the loop termination criteria proposed by the present invention.
Verification test (ii): FIG. 4 shows a gear vibration signal upon a broken tooth fault. The number of teeth of the gear is z-75, the modulus m-2, and the frequency f-n2And 60 is approximately equal to 13.6Hz, the vibration signal is decomposed by adopting an empirical mode decomposition method, and the screening termination condition of the empirical mode decomposition method adopts the criterion provided by the invention. The final result of the decomposition is 5 intrinsic mode function components imf1,imf2,…imf5。imf1Containing a large amount of information about gear failure, imf1Analysis was performed, imf1Is shown in fig. 5. From the 1 st intrinsic mode function component imf1The obvious modulation characteristic can be seen in the waveform of (1), the period T of the modulation wave is 0.074s, the frequency is approximately 13.6Hz, and the frequency is exactly equal to the rotating frequency of the fault gear. The empirical mode method based on the cycle termination criterion of the effective data segment provided by the invention can be used for accurately resolving the fault information of the actual gear vibration signal.
S2, processing the IMF components through an analytic mode decomposition algorithm, and reconstructing a mode; aiming at the mode confusion phenomenon of empirical mode decomposition, the analytic-empirical mode decomposition method is provided, the analytic mode decomposition method is used for decomposing the modes of the empirical mode decomposition method with mode confusion, and then the modes are reconstructed, so that the mode confusion phenomenon is reduced to a great extent.
An Analytical Mode Decomposition (AMD) method can separate harmonic components of each frequency band from a signal, and the principle of the method is as follows:
signal x0(t) from the dividing frequency fbThe decomposition is carried out in two parts: a fast varying signal xf(t) and a slowly varying signal xs(t)：
x0(t)＝xf(t)+xs(t) (7)
The corresponding Fourier transforms are respectively denoted as Xf(ω)、Xs(ω) that do not overlap in frequency band. Respectively to cos (2 pi f)bt)·x0(t) and sin (2 π f)bt)·x0(t) performing Hilbert transform to obtain:
H[cos(2πfbt)·x0(t)]＝xs(t)·H[cos(2πfbt)]+cos(2πfbt)·H[xf(t)] (8)
H[sin(2πfbt)·x0(t)]＝xs(t)·H[sin(2πfbt)]+sin(2πfbt)·H[xf(t)] (9)
deriving x from the above formulaf(t) and xsThe calculation method of (t) is as follows:
xs(t)＝sin(2πfbt)·H[cos(2πfbt)·x0(t)]-cos(2πfbt)·H[sin(2πfbt)·x0(t)] (10)
xf(t)＝x0(t)-xs(t) (11)。
the step S2 specifically includes the following sub-steps:
s21, sequentially recording the n intrinsic mode function components obtained by calculation in the step S1 as IMF1(t)～IMFn(t)；
For IMF1(t) performing fast Fourier transform on the component to obtain an amplitude spectrum thereof, and determining a boundary frequency f from the amplitude spectrumb1(ii) a For IMF1(t) performing analytical mode decomposition from the boundary frequency fb1The decomposition yields two signals c1(t) and
Figure GDA0003160195680000081
c1(t) represents IMF1(ii) a correction component of (t),
Figure GDA0003160195680000082
is c1(t) a residual component; namely:
Figure GDA0003160195680000083
s22, mixing
Figure GDA0003160195680000084
Addition to IMF2(t) obtaining an updated IMF2(t) as
Figure GDA0003160195680000085
To pair
Figure GDA0003160195680000086
Performing fast Fourier transform to obtain
Figure GDA0003160195680000087
And determining the boundary frequency f from the amplitude spectrumb2(ii) a To pair
Figure GDA0003160195680000088
Performing analytic mode decomposition with the boundary frequency fb2To obtain c2(t) and
Figure GDA0003160195680000089
wherein, c2(t) is
Figure GDA00031601956800000810
The correction component of (a) is,
Figure GDA00031601956800000811
is c2(t) a residual component; namely:
Figure GDA00031601956800000812
s23, removing IMFn(t) repeating the operations of steps S21 and S22 for all intrinsic mode function components except for (t);
s24 for IMFn(t) carrying out the following operationsThe method comprises the following steps: will be provided with
Figure GDA00031601956800000813
Addition to IMFnIn (t), obtaining IMF* n(t) for IMF* n(t) performing fast Fourier transform and plotting the amplitude spectrum thereof, determining the demarcation frequency fbk(ii) a For IMF* n(t) performing analytical mode decomposition from the boundary frequency fbkObtain a signal cn(t) and
Figure GDA00031601956800000814
namely:
Figure GDA00031601956800000815
will be provided with
Figure GDA00031601956800000816
Added to the residual function rn(t) obtaining the final residual error, denoted as vn(t)。
Verification test (iii): in order to prove the effectiveness and the advantages of the analytic-empirical mode decomposition method, the overall average empirical mode decomposition method and the analytic-empirical mode decomposition method are adopted to decompose the simulation signals respectively, and the decomposition results are compared.
The simulation signals x (t) are:
x(t)＝x1(t)+x2(t)+x3(t)；t∈[0,1] (15)
Figure GDA00031601956800000817
x2(t)＝2cos(60πt)；x3(t)＝2.5cos(28πt)；t∈[0,1] (17)
the simulated signal is shown in fig. 6. The results of the empirical mode decomposition are shown in fig. 7. As can be seen from FIG. 7, the resulting modal aliasing is more severe, the first oneMode C1(t) aliasing the high frequency intermittent signal with a 30Hz sinusoidal signal; second mode C2(t), the 30Hz sinusoidal signal is aliased with the 14Hz sinusoidal signal; the third and fourth modes are spurious intrinsic mode function components. The ensemble average empirical mode decomposition method has an ensemble average number N of 30, a noise amplitude of 0.01 times the standard deviation of the original signal, and the result of the ensemble average empirical mode decomposition is imf1、imf2、imf3、imf4As shown in fig. 8. The results of the analytical-empirical mode decomposition are shown in fig. 9. During the program run, the run times (6 run times) for the three methods were recorded, as shown in table 1. In table 1, EMD represents the empirical mode decomposition method, EEMD represents the ensemble average empirical mode decomposition method, and AEMD represents the analysis-empirical mode decomposition method.
TABLE 1 three methods run time (units s)
Method For the first time For the second time The third time Fourth time Fifth time The sixth time
EMD 1.3988 1.3860 1.3930 1.3846 1.3820 1.3940
EEMD 41.4875 44.8383 44.1800 42.8480 43.2321 44.024
AEMD 1.4027 1.3832 1.3843 1.3934 1.3880 1.3848
By comparing the decomposition results of the three methods, it can be found that both the ensemble empirical mode decomposition and the analytical empirical mode decomposition can suppress mode confusion, but the analytical empirical mode decomposition effect is better than that of the ensemble empirical mode decomposition method. By comparing the decomposition times of the three methods, it can be seen that the runtime of the analytical-empirical mode decomposition method is significantly lower than that of the ensemble average empirical mode decomposition method, and is similar to that of the empirical mode method. Therefore, the analytic-empirical mode decomposition has certain advantages.
Verification test (iv): fig. 10 is a vibration signal of a rotor rub fault. The rotating speed of the rotor is controlled by an input motor, and the rotating speed of the rotor is set to be 300 r/min. And the displacement sensor is horizontally arranged and used for measuring a vibration signal of the rotor with local rubbing fault. The sensor sampling frequency was set to 8 kHz. The vibration signal was decomposed using a modified method of empirical mode decomposition, local mean, and the results are shown in fig. 11. The results of the decomposition using the analytical-empirical mode decomposition of the present invention are shown in fig. 12. Comparing fig. 11 and fig. 12, it can be seen that the result mode confusion of the local mean method is more serious, and the mode confusion of the analysis empirical mode decomposition result is weaker, which indicates that the analysis empirical mode can suppress the mode confusion, and has a certain practical value for analyzing the rotor vibration signal.
It will be appreciated by those of ordinary skill in the art that the embodiments described herein are intended to assist the reader in understanding the principles of the invention and are to be construed as being without limitation to such specifically recited embodiments and examples. Those skilled in the art can make various other specific changes and combinations based on the teachings of the present invention without departing from the spirit of the invention, and these changes and combinations are within the scope of the invention.

Claims (2)

1.一种解析-经验模态分解方法，其特征在于，包括以下步骤：1. an analytical-empirical mode decomposition method, is characterized in that, comprises the following steps: S1、通过经验模态分解方法对原始信号进行筛分，得到n个IMF分量以及1个残余量；原始信号为发生断齿故障的齿轮振动信号或转子碰摩故障的振动信号；包括以下子步骤：S1. The original signal is sieved by the empirical mode decomposition method to obtain n IMF components and 1 residual quantity; the original signal is the vibration signal of the gear broken tooth fault or the vibration signal of the rotor rubbing fault; including the following sub-steps : S11、确定原始信号x(t)的所有极大值点与极小值点，利用三次样条曲线分别拟合出极值点的上、下包络线，并确保信号x(t)在上下包络线之间，然后计算出上下包络线的局部均值，记作m1；S11. Determine all the maximum points and minimum points of the original signal x(t), use the cubic spline curve to fit the upper and lower envelopes of the extreme points respectively, and ensure that the signal x(t) is up and down between the envelopes, and then calculate the local mean of the upper and lower envelopes, denoted as m 1 ; S12、计算出信号x(t)与局部均值的差值h1(t)：S12. Calculate the difference h 1 (t) between the signal x(t) and the local mean: h1(t)＝x(t)-m1 (1)；h 1 (t)=x(t)-m 1 (1); S13、判断h1(t)能否满足IMF成立的条件；若不满足，则将h1(t)作为新的原始信号，重复步骤S11和S12，直到h1(t)满足IMF成立的条件；S13. Determine whether h 1 (t) can satisfy the conditions for the establishment of the IMF; if not, take h 1 (t) as a new original signal, and repeat steps S11 and S12 until h 1 (t) satisfies the conditions for the establishment of the IMF ; S14、令C1(t)＝h1(t)，C1(t)即为分解得到的第一个IMF分量；S14. Let C 1 (t)=h 1 (t), and C 1 (t) is the first IMF component obtained by decomposition; S15、将C1(t)从信号x(t)中分离出来，得到的剩余信号记为：r1(t)＝x(t)-C1(t)；S15. Separate C 1 (t) from the signal x(t), and the obtained residual signal is denoted as: r 1 (t)=x(t)-C 1 (t); S16、将r1(t)作为新的原始信号，重复步骤S11～S15，进行n次有限的循环，得到n个IMF分量以及1个残余量rn(t)；即：S16, take r 1 (t) as a new original signal, repeat steps S11 to S15, and perform n limited cycles to obtain n IMF components and 1 residual r n (t); namely:
Figure FDA0003132473480000011
Figure FDA0003132473480000011
 原始信号x(t)即为：
Figure FDA0003132473480000012
其中，Ci(t)表示第i个内禀模态函数分量，rn(t)称为残差函数；The original signal x(t) is:
Figure FDA0003132473480000012
Among them, C i (t) represents the i-th intrinsic mode function component, and r n (t) is called the residual function; S2、通过解析模态分解算法对IMF分量进行处理，重构模态；包括以下子步骤：S2. Process the IMF components through an analytical modal decomposition algorithm to reconstruct the modal; including the following sub-steps: S21、将步骤S1计算得到的n个内禀模态函数分量依次记为IMF1(t)～IMFn(t)；S21, the n intrinsic mode function components calculated in step S1 are sequentially recorded as IMF 1 (t) ~ IMF n (t); 对IMF1(t)分量做快速傅里叶变换，获得其幅值谱，由幅值谱图确定分界频率fb1；对IMF1(t)进行解析模式分解，由分界频率fb1分解得到两个信号c1(t)和
Figure FDA0003132473480000013
c1(t)表示IMF1(t)的校正分量，
Figure FDA0003132473480000014
是c1(t)的剩余分量；即：Perform fast Fourier transform on the IMF 1 (t) component to obtain its amplitude spectrum, and determine the boundary frequency f b1 from the amplitude spectrum; perform analytical mode decomposition on IMF 1 (t), and decompose the boundary frequency f b1 to obtain two signals c 1 (t) and
Figure FDA0003132473480000013
c 1 (t) represents the correction component of IMF 1 (t),
Figure FDA0003132473480000014
is the remaining component of c 1 (t); that is:
Figure FDA0003132473480000021
Figure FDA0003132473480000021
 S22、将
Figure FDA0003132473480000022
添加到IMF2(t)中，得到更新后的IMF2(t)，记作IMF2 *(t)；对IMF2 *(t)进行快速傅里叶变换，得到IMF2 *(t)的幅值谱，并由幅值谱图确定分界频率fb2；对IMF2 *(t)进行解析模式分解，由分界频率fb2得到c2(t)和
Figure FDA0003132473480000023
其中，c2(t)就是IMF2 *(t)的校正分量，
Figure FDA0003132473480000024
是c2(t)的剩余分量；即：S22, will
Figure FDA0003132473480000022
Add to IMF 2 (t) to obtain the updated IMF 2 (t), denoted as IMF 2 * (t); perform fast Fourier transform on IMF 2 * (t) to obtain IMF 2 * (t) amplitude spectrum, and determine the boundary frequency f b2 from the amplitude spectrum; perform analytical mode decomposition on IMF 2 * (t), and obtain c 2 (t) and
Figure FDA0003132473480000023
where c 2 (t) is the correction component of IMF 2 * (t),
Figure FDA0003132473480000024
is the remaining component of c 2 (t); that is:
Figure FDA0003132473480000025
Figure FDA0003132473480000025
 S23、对除IMFn(t)以外的所有内禀模态函数分量重复步骤S21的和S22的操作；S23, repeat the operations of step S21 and S22 for all intrinsic mode function components except IMF n (t); S24、对于IMFn(t)进行以下操作：将
Figure FDA0003132473480000026
添加到IMFn(t)中，得IMF* n(t)，对IMF* n(t)做快速傅里叶变换并画出其幅值谱，确定分界频率fbk；对IMF* n(t)进行解析模式分解，由分界频率fbk得到信号cn(t)和
Figure FDA0003132473480000027
即：S24. Perform the following operations for IMF n (t): set the
Figure FDA0003132473480000026
Add to IMF n (t), get IMF * n (t), perform fast Fourier transform on IMF * n (t) and draw its amplitude spectrum, determine the boundary frequency f bk ; for IMF * n (t ) for analytical mode decomposition, and the signals c n ( t ) and
Figure FDA0003132473480000027
which is:
Figure FDA0003132473480000028
Figure FDA0003132473480000028
 将
Figure FDA0003132473480000029
添加到残差函数rn(t)中，得到最终的残差，记作vn(t)。Will
Figure FDA0003132473480000029
Add to the residual function r n (t) to get the final residual, denoted as v n (t). 2.根据权利要求1所述的一种解析-经验模态分解方法，其特征在于，所述步骤S16中采用基于有效数据段的筛分终止准则来计算循环次数n，具体操作方法为：当内禀模态函数Cn(t)满足下列条件时，停止循环：2. a kind of analytical-empirical mode decomposition method according to claim 1, is characterized in that, adopts the sieving termination criterion based on valid data segment in described step S16 to calculate the number of cycles n, and the concrete operation method is: when The loop stops when the intrinsic modal function C n (t) satisfies the following conditions:
Figure FDA00031324734800000210
Figure FDA00031324734800000210
 式中，
Figure FDA00031324734800000211
表示Cn(t)的上下包络均值曲线的有效段，δ表示偏差阈值，取值范围为0.01～0.1。In the formula,
Figure FDA00031324734800000211
Indicates the effective segment of the upper and lower envelope mean value curves of C n (t), and δ represents the deviation threshold, which ranges from 0.01 to 0.1.
CN201711049787.0A 2017-10-31 2017-10-31 Analytic-empirical mode decomposition method Expired - Fee Related CN107748734B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201711049787.0A CN107748734B (en) 2017-10-31 2017-10-31 Analytic-empirical mode decomposition method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201711049787.0A CN107748734B (en) 2017-10-31 2017-10-31 Analytic-empirical mode decomposition method

Publications (2)

Publication Number Publication Date
CN107748734A CN107748734A (en) 2018-03-02
CN107748734B true CN107748734B (en) 2021-08-27

Family

ID=61252821

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201711049787.0A Expired - Fee Related CN107748734B (en) 2017-10-31 2017-10-31 Analytic-empirical mode decomposition method

Country Status (1)

Country Link
CN (1) CN107748734B (en)

Families Citing this family (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108508268A (en) * 2018-04-02 2018-09-07 山西大学 A kind of novel digital type flickermeter decomposed based on time-frequency domain
CN108596215B (en) * 2018-04-04 2020-11-17 中国航发湖南动力机械研究所 Multi-modal signal analysis and separation method, device, equipment and storage medium
CN109567743A (en) * 2018-10-24 2019-04-05 加康康健有限公司 A kind of signal reconstruction method, apparatus, terminal equipment and storage medium based on EMD
CN109582913B (en) * 2018-11-20 2023-08-22 国家海洋局第一海洋研究所 Method for empirical mode decomposition screening iteration process termination criterion
CN109859174B (en) * 2019-01-09 2020-09-25 东莞理工学院 An OLED defect detection method based on empirical mode decomposition and regression model
CN109799535B (en) * 2019-03-14 2020-10-30 中船海洋探测技术研究院有限公司 Filtering method for full-tensor magnetic gradient positioning detection data
CN110096673B (en) * 2019-04-29 2023-03-14 河北工业大学 EMD improvement method suitable for signal decomposition
CN110441063B (en) * 2019-06-12 2021-04-27 祝思宁 Method for monitoring and diagnosing cracks of large high-speed rotor shaft
CN111537989B (en) * 2020-03-25 2022-07-15 中国电子科技集团公司第二十九研究所 Method for extracting signal micro Doppler modulation component based on empirical mode decomposition
CN111553308A (en) * 2020-05-11 2020-08-18 成都亿科康德电气有限公司 Reconstruction method of partial discharge signal of power transformer
CN114936375B (en) * 2022-06-02 2024-06-04 北京理工大学 Asymmetric integrated optical encryption system based on two-dimensional empirical mode decomposition

Citations (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102778357A (en) * 2012-08-15 2012-11-14 重庆大学 Mechanical failure feature extracting method based on optimal parameter ensemble empirical mode decomposition (EEMD)
CN103617684A (en) * 2013-12-12 2014-03-05 威海北洋电气集团股份有限公司 Interference type optical fiber perimeter vibration intrusion recognition algorithm
CN103870694A (en) * 2014-03-18 2014-06-18 江苏大学 Empirical mode decomposition denoising method based on revised wavelet threshold value
CN104007315A (en) * 2014-04-29 2014-08-27 昆明理工大学 Improved empirical mode decomposition processing method
CN104006961A (en) * 2014-04-29 2014-08-27 北京工业大学 Cycloid bevel gear fault diagnosis method based on empirical mode decomposition and cepstrum
CN105893773A (en) * 2016-04-20 2016-08-24 辽宁工业大学 Vibration signal frequency characteristic extraction method based on EEMD technology, CMF technology and WPT technology
CN105938542A (en) * 2016-03-16 2016-09-14 南京大学 Empirical-mode-decomposition-based noise reduction method for bridge strain signal
CN106019102A (en) * 2016-06-27 2016-10-12 国网北京市电力公司 Signal de-noising method and apparatus
CN106610918A (en) * 2015-10-22 2017-05-03 中央大学 Empirical mode decomposition method and system for adaptive binary, conjugate mask grid
CN106844293A (en) * 2017-02-03 2017-06-13 中国铁道科学研究院 A kind of adaptive decoupling method of modal overlap problem in empirical mode decomposition
CN106874549A (en) * 2017-01-10 2017-06-20 西安理工大学 A kind of discrete distribution parabolic equation method in the arrowband of high-precision forecast ASF
CN107064752A (en) * 2017-03-22 2017-08-18 北京航空航天大学 A kind of distinguished number of aviation fault electric arc detection
CN107255563A (en) * 2017-06-28 2017-10-17 石家庄铁道大学 Realize gear-box mixed fault signal blind source separation method

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
TWI454248B (en) * 2008-09-23 2014-10-01 Ind Tech Res Inst Method of multi-dimensional empirical mode decomposition for image morphology
TW201217993A (en) * 2010-10-20 2012-05-01 Huafan University employing operation on decomposed matrices to reduce operation amount for single matrix per unit time for light-weighting matrix operation process in simpler operation circuit

Patent Citations (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102778357A (en) * 2012-08-15 2012-11-14 重庆大学 Mechanical failure feature extracting method based on optimal parameter ensemble empirical mode decomposition (EEMD)
CN103617684A (en) * 2013-12-12 2014-03-05 威海北洋电气集团股份有限公司 Interference type optical fiber perimeter vibration intrusion recognition algorithm
CN103870694A (en) * 2014-03-18 2014-06-18 江苏大学 Empirical mode decomposition denoising method based on revised wavelet threshold value
CN104007315A (en) * 2014-04-29 2014-08-27 昆明理工大学 Improved empirical mode decomposition processing method
CN104006961A (en) * 2014-04-29 2014-08-27 北京工业大学 Cycloid bevel gear fault diagnosis method based on empirical mode decomposition and cepstrum
CN106610918A (en) * 2015-10-22 2017-05-03 中央大学 Empirical mode decomposition method and system for adaptive binary, conjugate mask grid
CN105938542A (en) * 2016-03-16 2016-09-14 南京大学 Empirical-mode-decomposition-based noise reduction method for bridge strain signal
CN105893773A (en) * 2016-04-20 2016-08-24 辽宁工业大学 Vibration signal frequency characteristic extraction method based on EEMD technology, CMF technology and WPT technology
CN106019102A (en) * 2016-06-27 2016-10-12 国网北京市电力公司 Signal de-noising method and apparatus
CN106874549A (en) * 2017-01-10 2017-06-20 西安理工大学 A kind of discrete distribution parabolic equation method in the arrowband of high-precision forecast ASF
CN106844293A (en) * 2017-02-03 2017-06-13 中国铁道科学研究院 A kind of adaptive decoupling method of modal overlap problem in empirical mode decomposition
CN107064752A (en) * 2017-03-22 2017-08-18 北京航空航天大学 A kind of distinguished number of aviation fault electric arc detection
CN107255563A (en) * 2017-06-28 2017-10-17 石家庄铁道大学 Realize gear-box mixed fault signal blind source separation method

Non-Patent Citations (9)

* Cited by examiner, † Cited by third party
Title
A new machine learning approach to select adaptive IMFs of EMD;Uddin Md Burhan 等;《2016 2nd International Conference on Electrical, Computer & Telecommunication Engineering (ICECTE)》;20170316;1-4 *
Adaptive Empirical Mode Decomposition for Bearing Fault Detection;Van Tuan Do 等;《Journal of Mechanical Engineering》;20160321;第62卷(第5期);281-290 *
Gear Fault Diagnosis Based on Improved EMD Method and the Energy Operator Demodulation Approach;Han, J. 等;《Journal of Changsha University of Science & Technology(Natural Science)》;20150228(第12期);66-71 *
Vibration Analysis as a Diagnosis Tool for Health Monitoring of Industrial Machines;Peiming Shi 等;《Shock and Vibration》;20160209;1-11 *
基于EEMD和Hilbert变换的齿轮箱故障诊断;林近山;《机械传动》;20100515;第34卷(第05期);62-64 *
基于EMD与分批估计的动态称量快速融合方法;蒋婷 等;《仪器仪表学报》;20150615;第36卷(第06期);1406-1414 *
基于形态奇异值分解和经验模态分解的滚动轴承故障特征提取方法;汤宝平 等;《机械工程学报》;20100305;第46卷(第05期);37-42 *
自适应总体平均经验模式分解及其在行星齿轮箱故障检测中的应用;雷亚国 等;《机械工程学报》;20131114;第50卷(第03期);64-70 *
非平稳信号分析的广义解析模态分解方法;郑近德 等;《电子学报》;20160615;第44卷(第06期);1458-1464 *

Also Published As

Publication number Publication date
CN107748734A (en) 2018-03-02

Similar Documents

Publication Publication Date Title
CN107748734B (en) Analytic-empirical mode decomposition method
CN112446323B (en) HHT harmonic analysis method based on improved EMD modal aliasing and endpoint effect
CN103941091A (en) Power system HHT harmonious wave detection method based on improved EMD end point effect
CN107024718A (en) Poststack earthquake fluid Forecasting Methodology based on CEEMD SPWVD Time-frequency Spectrum Analysis
CN111222088B (en) An Improved Flat-Top Self-Convolution Window Weighted Power Harmonic Amplitude Estimation Method
CN110095650A (en) The complicated harmonic detecting analysis method of four spectral line interpolation FFTs based on five Rife-Vincent (I) windows
CN109783767B (en) Self-adaptive selection method for short-time Fourier transform window length
CN108764073A (en) A kind of acceleration filter of combination spectrum energy form fitting is made an uproar and integration method
CN103424183A (en) Method for eliminating abnormal interference on detection for mechanical vibration signals
CN102254103A (en) Method for determining screening time in self-adaptive ensemble empirical mode decomposition (EEMD)
CN105067025A (en) Method for utilizing monostable system stochastic resonance effect to detect weak signals
Li et al. A UV-visible absorption spectrum denoising method based on EEMD and an improved universal threshold filter
CN104597502A (en) Novel petroleum seismic exploration data noise reduction method
CN103674235B (en) Single-frequency alarm sound feature detection method based on short-time Fourier transform
CN105509771A (en) Signal noise reduction method of engine lubricating oil metal particle on-line monitor
CN110287853B (en) Transient signal denoising method based on wavelet decomposition
CN103488882B (en) A kind of hydraulic structure safety automation monitoring system noise-reduction method
Yan et al. Adaptive linear chirplet synchroextracting transform for time-frequency feature extraction of non-stationary signals
CN111596358B (en) Multiple suppression method, device and system
CN104569581A (en) Multi-level set and single-cycle estimation method of power grid frequency measuring
CN105678049A (en) End effect suppression method based on improved SVR continuation
CN106909787B (en) A kind of spark plug gap prediction technique and device
CN108680782A (en) The voltage flicker parameter detection method decomposed based on extreme point symmetric pattern
CN105466710A (en) Local mean value decomposition end point effect improvement method based on frequency domain similarity
CN110703089A (en) Wavelet threshold denoising method for low-frequency oscillation Prony analysis

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
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20210827