
On the HilbertHuang Transform Theoretical DevelopmentsSemion Kizhner, Karin Blank, Thomas Flatley, Norden E. Huang, David Petrick and Phyllis Hestnes National Aeronautics and Space Administration Goddard Space Flight Center Greenbelt Road, Greenbelt MD, 20771 3012861294 Semion.Kizhner1@nasa.gov AbstractOne of the main heritage tools used in scientific and engineering data spectrum analysis is the Fourier Integral Transform and its high performance digital equivalent – the Fast Fourier Transform (FFT). Both carry strong apriori assumptions about the source data, such as linearity, of being stationary, and of satisfying the Dirichlet conditions. A recent development at the National Aeronautics and Space Administration (NASA) Goddard Space Flight Center (GSFC), known as the HilbertHuang Transform (HHT), proposes a novel approach to the solution for the nonlinear class of spectrum analysis problems. Using aposteriori data processing based on the Empirical Mode Decomposition (EMD) sifting process (algorithm), followed by the normalized Hilbert Transform of the decomposition data, the HHT allows spectrum analysis of nonlinear and nonstationary data. The EMD sifting process results in a nonconstrained decomposition of a source real value data vector into a finite set of Intrinsic Mode Functions (IMF). These functions form a near orthogonal adaptive basis, a basis that is derived from the data. The IMFs can be further analyzed for spectrum interpretation by the classical Hilbert Transform. A new engineering spectrum analysis tool using HHT has been developed at NASA GSFC, the HHT Data Processing System (HHTDPS). As the HHTDPS has been successfully used and commercialized, new applications post additional questions about the theoretical basis behind the HHT and EMD algorithms. Why is the fastest changing component of a composite signal being sifted out first in the EMD sifting process? Why does the EMD sifting process seemingly converge and why does it converge rapidly? Does an IMF have a distinctive structure? Why are the IMFs near orthogonal? We address these questions and develop the initial theoretical background for the HHT. This will contribute to the developments of new HHT processing options, such as realtime and 2D processing using Field Programmable Array (FPGA) computational resources, enhanced HHT synthesis, and broaden the scope of HHT applications for signal processing. ^ One of the main heritage tools used in scientific and engineering data spectrum analysis is the Fast Fourier Transform that carries strong apriori assumptions about the source data, such as linearity and of being stationary. The data must originate from a periodic waveform and also satisfy the Dirichlet conditions of having a finite number of discontinuities and extremas, and be integrable in any sequence of time intervals with length of the period T [8]. Heritage spectrum analysis methods use a fixed basis in their transforms while EMD derives its basis adaptively from the data itself. The EMD sifting process [1][3] is a novel algorithm for digital signal processing of nonlinear and nonstationary data. Given an arbitrary input vector of rational numbers, the EMD algorithm invariably sifts out IMF components of different time scales with the fastest varying component being sifted out first. This is accomplished by a process known as “sifting” which is repeatedly applied to the signal until it converges on criteria that defines an IMF. In addition, the EMD process itself converges and produces a finite number of IMFs. All this has been observed in all tests conducted so far, in other words, empirically. However, it is not that obvious why the EMD algorithm behaves this way and this paper addresses these questions by using the following methodology. 1.2 MethodologyWe formulated the EMD algorithm by an indepth procedural description of the EMD algorithm, as it is implemented in the up to date HHTDPS Release 1.4. The EMD algorithm, in addition to the original publications and patents [1][3], was presented now quiet a few times, for example in [4][5]. However, an overview of it is given here once more, but for its implementation in HHTDPS, and being also formulated with different procedural details for the purpose to use it in the following research of the theoretical bases of why the EMD empirical algorithm works. We then formally stated the research problem as the need of understanding and developing the theoretical foundations of the EMD algorithm empirical behavior. This is followed by examination of a few analogies and intuitive examples of signals containing artificially created fast and slow varying components. We also consider a general case for the EMD sifting process and establish the theoretical foundations of the EMD algorithm’s sifting sequence of scales, IMF structure and EMD theoretical convergence. We accomplish this by providing a hypothesis for why the fastest scale is sifted out first and report two hypotheses about the sifting process’s symmetric pair invariance and resulting IMF structure, and its theoretical rapid convergence in diminishing amplitude regions. We then propose a few new applications based on this research, as well as define the areas of future research work on the HHT. ^ The EMD algorithm’s empirical behavior is determined by its builtin definitions and criterias as well as by the user’s supplied run configuration vector Ψ as described in more details in Section 2.1 Ψ = Ψ(Dt, m, k, p, …) (1) Ψ is composed of the sampling time interval Dt (used after EMD in spectrum analysis) and empirical parameters supplied by a user when running HHTDPS. Among these parameters are m, the maximum number of IMFs to strive for, k is the maximum allowable number of EMD sifts for one IMF, and p, which enables a user to select the “beyondtheenvelopeendpoints” prediction algorithm option. The following formulation of the Empirical Mode Decomposition algorithm is describing the EMD implementation in the latest HHTDPS Release 1.4. ^ The EMD algorithm is based on a few distinct procedural steps (a)  (h) as follows: EMD Algorithm Entry Point User specifies EMD Configuration Vector Ψ = Ψ(Dt, m, k, p, …) and input signal s(t) sampling vector – a vector of signal energy values s(t_{i}), where t_{i} is an equally spaced sequence of time arguments for i[1,…, N]. InitializationFor the first iteration of this algorithm beginning in Step (a) an initialization takes place. This results in the “data parent” p(t_{i}) vector for the EMD process iteration loop becoming the original signal s(t_{i}). The input signal values s(t_{i}) are assigned to residue R(t_{i}) (which becomes the parent signal p(t_{i}) in Step (a) and p(t_{i}) is later reset to running residue r(t_{i}) in Step (g)) for usage in Step (b)). The IMF storage array is cleared to zeros. j:=0, where j is the number of IMFs found R(t_{i}) := s(t_{i}), i Î (1,…,N) IMF := 0 (a) EMD IMF Registry and Completion Criteria Check Entry Point IMF Registry R(t_{i}) := R(t_{i}) – IMF_{j}(t_{i}), i Î (1,…,N) ^ R(t_{i}) is checked against the Completion Criteria where E(R(t_{i})) is the number of local extremas in R(t_{i}) and m is the user supplied parameter in configuration vector (1) ( (E(R(t_{i}) <= 3) OR (j = m) ) (2) If the above Completion Criteria is satisfied the EMD algorithm completes with the last signal residue R(t_{i}) becoming the process residue from which an IMF couldn’t be made. ^ . If the above Completion Criteria is not met p(t_{i}) := R(t_{i}) #Sifts:=0 and the EMD process continues with the following iteration loop (b) EMD Sifting Process Iterative Loop Entry Point The first step of the algorithm is to search the parent dataset and find the local maximum and minimum values, also known as extrema. These points form the set E(p(t)) which is divided into two sets, p^{max} and p^{min}, containing the respective local maximum and minimum extrema values. (c) The two subsets for both extrema point types p^{max} and p^{min} are splined using piecewise cubic splines, comprising the upper and lower boundary of the extrema subsets’ envelope boundaries u(t) and d(t). The use of the cubic spline provides a relatively slow changing background median M(t) cubic spline (as described below) against which local fast variations of the signal become prominent and which, in turn, forms the basis for the most variant signal being sifted out first by the EMD sifting process. But first the extrema sets are extended a bit beyond [t_{1}, t_{N}] to always maintain data over the original argument interval, as described next. (d) The sets u(t), d(t) are extended a bit beyond the data original argument interval [t_{1}, t_{N}] with a few predicted maxima and minima extrema points on both ends of the original argument interval [t_{1}, t_{N}]. There are many prediction algorithms, meaning there is no accepted one. The extended extrema sets are then splined beyond original argument interval [t_{1}, t_{N}] using predicted extrema points and then resampled on [t_{1}, t_{N}]. This original argument interval is always maintained throughout the HHTDPS. (e) As stated above in (c), splines u(t), d(t) are then resampled for sampling time arguments t_{i}, where i Î [1,…,N], and the discrete median vector M(t_{i}) is computed as M(t_{i}) = (u(t_{i}) + d(t_{i}) ) / 2 i Î (1,…,N) (3) This definition of median is essential, because it assures the sifting process convergence. When the median computation is iterated ktimes in the sifting process, the divisor becomes 2^{k}, which rapidly becomes a large number. The convergence theoretics are based on this number and are described in more details in Section 3. (f) The difference r(t_{i}) = p(t_{i}) – M(t_{i}) is formed and the sifting residue r(t_{i}) is checked against the IMF criteria, as follows. (g) The IMF Criteria Check The difference between the number E(r(t_{i})) of extrema points in r(t_{i}), and the number Z(r(t_{i})) of its zerocrossings or the change in two adjacent extrema points’ magnitude sign, may not vary by more than 1. Additionally, each IMF must have at least 3 extrema points ( (E(r(t_{i}) > 3) && (E(r(t_{i}),) – Z(r(t_{i})) <= 1) ) OR (#Sifts=k) i[1,…,N] (4) An r(t_{i}), which satisfies the IMF Criteria, is called an IMF. Parameter k is user defined in configuration vector (1). The r(t_{i}) is treated as an IMF if the limit of sifts k is reached. If the r(t_{i}) satisfies the IMF Criteria it is stored during this EMD algorithm verification of the IMF Criteria, the number of IMFs j is incremented by 1 and control is passed to the EMD IMF Registry Entry Point and Exit Criteria Check (a) where the IMF is subtracted (registered) from the signal residue R(t) to obtain the EMD next sifting parent signal p(t). This signal residue becomes the new input parent to the EMD sifting process iteration step. Note, that at no time does a median M(t) become an input to the EMD sifting process iterative loop. If the sifting residue r(t_{i}) does not satisfy the IMF Criteria p(t_{i}) := r(t_{i}) The number of sifts #Sifts is incremented by 1 and control is passed to EMD algorithm iteration loop entry point (b). (h) EMD Algorithm End. ^ Naturally, because of the EMD algorithm construct described above, the sum of all IMFs and the last signal residue R(t) (which is counted towards the number of IMFs) synthesize the original input signal s(t) s(t) = ∑IMF_{l} + R(t), where 1<=l<=m1 (5) The set of IMFs, which is derived from the data, comprises the signal s(t) nearorthogonal adaptive basis and is used for the following signal timespectrum analysis. Additional details concerning the EMD can be found in references [4][5]. With the EMD algorithm described above, the research problem is to understand why it works this way and try to develop the theoretical fundamentals of this algorithm. ^ The fastest changing component of a composite signal is invariably being sifted out first in the Empirical Mode Decomposition algorithm and the first hypothesis in the development of the EMD algorithm theory is related to the question: “why?” This empirical fact is significant in the respect that it provides the first insight into how the EMD algorithm works. Also the EMD algorithm’s computational performance depends on the number of extrema in the input data set and taking out first all the extrema for the fastest varying component greatly contributes to the EMD algorithm’s next IMF sifting steps computational performance. Hypothesis 1: Assuming theoretical convergence of the EMD sifting process, the fastest scale is being sifted out first, because the composite signal s(t) extremas’ envelope median is approximating the slower variance signal in presence of a fast varying component. It is implied in this Section that the EMD sifting process converges theoretically in all cases. The EMD sifting process rapid theoretical convergence hypotheses are formulated in this paper and will be reported in future papers. In order to prove this hypothesis we are first considering two analogies to the EMD sifting process  one from optical physics and the other from electrical and electronics engineering disciplines. These analogies provide initial insight into why the signal fastest changing component is being sifted out first by the EMD sifting process. We then consider three examples of the EMD sifting for a few artificially created signals comprising of fast and slow varying components. These analogies and intuitive examples allow us to gain an initial insight into the mechanism of why the fastest scale gets sifted out first. We then examine the EMD sifting of an arbitrary signal and prove the hypothesis for a general case. ^ The first analogy has to do with light. Light spectrum analysis, which has played such an important part in astronomical work, is essentially a method of ascertaining the nature of a remote celestial body by a process of sifting or analyzing into different components the light received from them [9]. It was first clearly established by Newton that ordinary white light is comprised of waves of different frequencies. The prism sifts out different colors. This known empirical phenomenon has a wellestablished theoretical basis in that the prism bends light of different wavelength by a different degree ρ, resulting in a distinctive order of output light spectrum (colors), with higher frequency components appearing “ first” in the prism spectrum. It results in a definitive sequence of colors at the output of a prism. ^ The second analogy has to do with and electrical RCchain (circuitry). This circuitry comprises a resistance (R) and capacitor (C), with the RCchain characteristic constant τ = RC as depicted in Picture 1. This circuitry is often used to smooth (filter out) fast variable fluctuations of a constant input voltage U_{Input}. The functioning of the RCchain as a filter depends on its user selected hardware parameter, the characteristic constant τ. Fluctuation R U_{Input } C_{ }U_{Output} Picture 1. RCchain The RCchain’s transfer function K(ω) [7] can be described as the ratio of the output voltage to input voltage, with both voltages considered as a function of frequency ω: K(ω) = U(ω)_{Output} _{ }/ U_{Input}(ω) = 1 / √(1 + ω^{2 }C^{2 }R^{2}) or K(ω) = 1 / √(1 + ω^{2 }τ^{2}) (6) For a near constant input voltage (ω=0, absence of fluctuations) the capacitor impedance X_{C} = 1/ωC = 1/0*C = ∞ and K(0) = 1. As ω → ∞ increases, the capacitor impedance to it is decreasing, Xc → 0, and the amplitude of the output voltage variable component factor K(ω) is approaching 0, as shown on the graph in Picture 2. K(ω) 1 0.7071 0 ω_{B} ω Picture 2. RCchain Transfer Function If we arbitrarily select the output voltage amplitude as a fraction of the input voltage, for example, as K(ω) = 0.707 = 1/√2 (this threshold value of 1/√2 is selected to get the following convenient computations and is actually the root mean square (rms) power of a sinusoid voltage function over its period), then the corresponding frequency band ω_{B} or filter upper boundary can be evaluated from the transfer function as follows: 1 / √2 = 1 / √ (1 + ω^{2}_{B }τ^{2}) or (1 + ω^{2}_{B }τ^{2}) = 2 → ω^{2}_{B }τ^{2} = 1 → ω_{B }τ = 1 or ω_{B} = 1 / τ Frequencies higher than ω_{B} are filtered out (sifted out) by this RCchain filter and the stability of the output is increasing with the increase in the RCchain’s constant τ. ^ The EMD sifting process is computationally analogous to the “hardware” phenomenon of a prism sifting out white light into a distinctive sequence of components of different frequencies, from highest frequency in violet (790 Tera Hertz  Thz) to blue, cyan, green, yellow, orange and to lower frequency in red (480 Thz). If the bottom face of the prism is horizontal, the light beam entering from the bottom of the prism is not deviated. When the light leaves the prism at the inclined face and enters the air, it is refracted, and the beam is deflected to the right; the deflection is larger at shorter, bluer, wavelengths. The EMD sifting process is also analogous to the RCchain that is sifting out the input voltage fast varying fluctuations of frequencies exceeding the less varying or a constant input voltage, as determined by the RCchain characteristic constant τ. Analogous to these two processes, with their characteristic parameters ρ and τ, the EMD sifting process is defined both by its implementation definitions and criterions in the HHTDPS, as well as by its user supplied empirical runtime configuration vector Ψ. The EMD sifts out the input signal’s components in a definitive sequence of scales, analogous to that of a prism sifting out the light waveform into a sequence from shortest to longest wavelength, and performs it in a way similar to that of an RCchain filtering out the frequency band around the slow signal component or the signal intrinsic median level (s), as further elaborated below. ^ Intuitively, the EMD sifting process can be initially rationalized in the following way. Visualize a composite signal comprising two components of which one is a highly variable component and the other is a relatively slow varying trendlike component. It is obvious that the fast varying component will be observable as a curve following in the shape of the slow varying component in the composite signal, similar to the multitudes of small inlets and bay fractals in the broad outline of an ocean shore, when observed from afar at sea, or from Space. In other words, the slow varying component can be interpreted as an approximate median of the composite signal’s envelope built on fast varying component extrema points with maybe a few standout extrema points (accidental large amplitude spikes) belonging to the slow varying component. The signal extremas envelope’s upper and lower boundaries are defined by the EMD HHTDPS implementation as two piecewise cubic splines, connecting the signal’s adjacent maxima or minima points correspondingly. The piecewise curves are also smoothly connected, meaning their slopes must be contiguous and be of the same value at juncture points. The EMD sifting process computes the composite signal envelope upper and lower boundaries’ median as the algebraic sum of the envelope upper and lower curves’ data at sampling time points divided by 2, and subtracts the resulting median from the parent signal to arrive at the next “sifted out” residue component candidate, as defined above in the EMD sifting algorithm overview (Section 2 equation (3)). Intuitively, the EMD essentially subtracts the slow varying component from the composite signal, yielding the signal’s fastest varying component. This sifting residue is then checked against the IMF criteria and the process is repeated when necessary until the fastest varying IMF is found or made from the input composite signal. Before proving the proposed Hypothesis 1 and elaborating on the theoretical details of the EMD sifting process, let us first consider three intuitive examples, using a few artificially constructed composite signals for which the workings and the results of the EMD sifting process are known in advance due to the way in which the input signals were constructed. ^ The purpose of the following three examples is to demonstrate how the EMD sifting process works for data when the results of the sifting process are known intuitively and a priori. We know these results beforehand, because of the artificial way the input data in these examples was constructed. We know that the envelope of a fast varying sinusoid signal is obviously composed of two straight lines and its median is a straight line too. If we include in the signal a slower varying signal as a straight line the median intuitively will be approximated by this straight line. Considering then the low varying component a straight line (which approximates the envelope median), makes the composite signal sifting process’ results intuitively clear, because the subtraction of the straight line median results in the fast varying sinusoid component. Following this line we constructed three composite signals with the purpose of observing how they are processed by the EMD. 3.5.1 Example 1Consider a composite signal comprised of a known constant signal s1=0.5 (nonvariant signal presented geometrically by a horizontal straight line) and a known relatively fast (as compared to s1) varying signal s2 that is a 5 Hz sinusoid with amplitude in the range [1.0, 1.0] s(t) = s1 + s2 = 0.5 + 1.0*cos(2*pi*f*t) The EMD of this signal results in sifting out first the 2 Hz sinusoid followed by the straight line s1=0.5 residue. 3.5.2 Example 2We consider now a composite signal comprised of a known and relatively slow varying signal s1(t) = k*t representing an inclined straight line with slope k and a second known signal that is a 5 Hz sinusoid with amplitude 1.0, identical to s2(t) in Example 1: s(t) = 1*t + 1.0*cos(2*pi*5*t) The EMD of this signal results in sifting out first the 5 Hz sinusoid followed by the line s1=t residue. 3.5.3 Example 3Finally, we consider the third intuitive example, a composite signal s comprising four components s1, s2, s3 and bias b1. Namely these are a 1 Hz, 2 Hz and 50 Hz sinusoids (cos function) and a constant bias with amplitude 1.0: s1 = b1 + 1.0*cos(2*pi*1*t); s2 = 2.0*cos(2*pi*2*t)); s3 = 2.0*cos(2*pi*50*t); s = s1 + s2 + s3 The EMD of this signal results in sifting out first the 50 Hz sinusoid followed by the 2 Hz sinusoid and residue equal to component s1. ^ In a general case, we now consider an arbitrary composite signal s(t). We will prove the Hypothesis 1 for an arbitrary signal by demonstrating the reason why the fastest varying component in an arbitrary signal s(t) is being sifted out first by the EMD sifting process. We are assuming the sifting process’ theoretical convergence. The theoretical convergence of the EMD sifting process is proved in the following Section 4. The only notrestrictive convention we make here is that we view an arbitrary signal being comprised of two components – a relatively fast or high variance component s_{h}(t) and a slower or low variance component s_{l}(t), similar to signals in the analogies and intuitive examples considered above in Sections 2.12.5 s(t) = s_{h}(t) + s_{l}(t) (7) This convention about the signal s(t) structure, again, is based on the two analogies we had described above, and is notrestrictive in any sense. There is no mystery of what s_{h}(t) may be. For example, the locations of the initial set of extremas of s(t) coincide, in general, with the location of the s_{h}(t) extremas, except an occasional spike rooted in the remaining components. The set of extremas of the parent signal characterizes the fastest scale at the time of the sifting process. ^ We will proceed with two instances of a general case signal s(t) represented in equation (7), before taking on the case of an arbitrary signal EMD processing directly. We will first consider an instance of in which signal has extrema symmetry in its fast varying component. We will then follow it up by a more general linear approximation case of the lower varying component, and complete the proof of Hypothesis 1. We will then formulate Hypothesis 2, 3. ^ We are only assuming in this instance of a general signal that the fast varying component’s amplitude is locally symmetric or that two adjacent maxima and minima are approximately equal in magnitude but different in direction or sign in the offset against the slow varying component. In other words, we assume that on an interval of interest tÎ[t1, t1 + dt] we have (max)s_{h}(t)  (min)s_{h}(t) ≈ 0 (8) The signals in the above Examples 13 were of kind. For such an instance of general symmetry we can approximate the signal maxima and minima envelope spline values s^{max}(t) and s^{min}(t) as s^{max}(t) ≈ s_{l}(t) + (max)s_{h}(t) s^{min}(t) ≈ s_{l}(t)  (min)s_{h}(t) for any t on interval [t1, t1 + dt]. The signal envelope’s median M(t) on [t, t + dt] is following along the signal slow varying component and can then be expressed as M(t) ≈ (s^{max}(t) + s^{min}(t) ) / 2 ≈ s_{l}(t) and s(t) – M(t) = s(t) – s_{l}(t) = s_{h}(t) → IMF. We rationalized in this instance (omitting the details of the validity of affects of proximity dependence on to fast scale definition in 9’) that the fast varying component’s amplitude local symmetry results in the otherwise arbitrary fast varying component are being sifted out first in the EMD sifting process. This is also a strong consequence of Hypothesis 2 proven below in Section 3. Namely, fast varying component’s adjacent extrema of opposite type are the first to be recognized as invariant pairs when they exist or they will be “polished” into local symmetry in a few iterative sifting steps, with the remaining parts of the signal to converge rapidly to a nearzero amplitude curve. The “polished off” part of the fast varying component together with nearzero convergence areas form the fastest scale IMF. When this IMF is subtracted from the parent signal the new remainder comprises the slower varying components of the composite signal. ^ In order to explain why the fastest scale is being sifted our first priority is to consider the approximation of the slow varying signal component s_{l}(t) on a small interval tÎ[t1, t1+ dt] by a straight line passing through the interval two end points (linear approximation is here actually the definition of a slow varying component) { (t1, s_{l}(t1)), (t1+ dt, s_{l}(t1+ dt) ) } or (s_{l}(t) – s_{l}(t1) ) / (tt1) = (s_{l}(t1+ dt) – s_{l}(t1) ) / dt = k or s_{l}(t) ≈ (s_{l}(t + dt) – s_{l}(t) ) / dt = s_{l}(t1) + k*(t – t1) or s_{l}(t) ≈ s_{l}(t1) + k*(t – t1) Slope k is negligibly small because of the selection of s_{l}(t) as the slow varying component and thus the straight line is being almost parallel to the Ot axis. A fast varying oscillation component combined with a line segment must satisfy above equality (9). This leads to a median point M(t) for any tÎ[t1, t1 + dt] being computed as M(t) = ( ( (s_{l}(t1) + k*(t – t1) ) + (max)s_{h}(t) ) + ( (s_{l}(t1) + k*(t – t1) ) + (min)s_{h}(t) ) ) / 2 = s_{l}(t1) + k*(t – t1) ≈ s_{l}(t1) ≈ s_{l}(t) or that the resulting median point is being situated approximately on the slower varying component s_{l}(t). When this is subtracted from the signal we get the fastest varying component of the signal s_{h}(t) → IMF. We proved that the fast varying component is being sifted out first in the EMD sifting process when the arbitrary slow varying component can be approximated by a piecewise straight line. ^ The selection of a piecewise straight line approximation for the signal envelope boundaries u(t), d(t), except for very small time intervals dt used in the above conceptual proof, is unproductive for an EMD implementation. Instead the EMD uses the piecewise cubic spline interpolation. Surprisingly, the above examined linear approximations of the low varying component of a composite signal on small intervals of the argument prove useful, when considering the implications of piecewise cubic splines selection for envelope boundaries. In the algorithm, we deal with discrete data and their envelopes being interpolated by EMD using a piecewise cubic spline function for adjacent extrema points of the same type, either maxima or minima. To distinguish these extrema data points among all signal data s(t) we will call them x(t). On an interval of two adjacent extrema of the same type [x_{i}, x_{j}], the function x(t) can be expressed by a cubic polynomial of variable t belonging to interval [0, 1] x(t) = a + b*t + c*t^{2} + d*t^{3} (9) For small values of t (in the neighborhood of t=0) the second and third terms in x(t) become very small and can be ignored. In general, x(t) can also be expressed as a sum of two terms, yielding x(t) = (a + bt) + (c*t^{2} + d*t^{3}) While the upper and lower envelopes represented by a piecewise cubic spline connecting maxima and minima extrema points on an interval [t1=t_{i}, (t1 + Δt)= t_{j}] can be conceptually viewed as a line segment (a + bt) or slow varying signal component modulated with a fast varying component (c*t^{2} + d*t^{3}). The same considerations hold for the extrema envelope lower boundary approximation by piecewise cubic splines. It is obvious that all three approximation linear segments – for the upper envelope, the lower envelope and the median are straight lines that are locally equidistant, with the median line following the slower varying signal as a trend. For a fast varying component the interval [t1=t_{i}, (t1 + dt)=t_{j}] between two adjacent arguments is very small and the linear part of the piecewise cubic spline is the envelopes’ and their median’ approximation on the entire interval and the proof for the general case signal is similar to one provided above in Section 3.6.1. However, two adjacent extrema points, as indicated by their two indices, of the same type x_{i}, x_{i+1=j} may be widely separated among all signal data points and this case requires further considerations as follows. The coefficients {a, b, c, d} for the piece of a cubic spline passing through point pairs x_{i} and x_{i+1} for parameter arguments t_{i}=0 and t_{i+1}=1 and 0<=i<=N1 can be computationally determined from the four conditions for these two obvious end points and two more controversially selected slopes of x(t) at interval ends, namely k_{i} and k_{i+1}: x_{i} = a (at given argument t=0 and signal data x_{i}) x_{i+1} = a + b + c + d (at given argument t=1 and signal data x_{i +1}) k_{i} = b (the first derivative of x(t) b + 2ct + 3dt^{2} at t=0) k_{i+1} = b + 2c + 3d (the first derivative of x(t) b + 2ct + 3dt^{2} at t=1) The solution for the cubic spline polynomial coefficients {a, b, c, d} in terms of the two known extrema data point magnitudes xi and x_{i+1}, and EMD process selected finite slopes k_{i}, k_{i+1}, are: a = x_{i} b = k_{i} (10)_{ } c = 3(x_{i+1} – x_{i}) – 2k_{i} – k_{i+1} d = 2(x_{i} – x_{i+1} ) + k_{i} + k_{i+1} For example, if we consider the two slopes being zero or k_{i} = k_{i+1} = 0 (as they should be for s(t) at any extrema point argument), then the solution determines the following polynomial x (t) = x_{i} + 3(x_{i+1} – x_{i})t^{2} + 2(x_{i} – x_{i+1})t^{3} (11) or x(t) = x_{i} + (x_{i+1} – x_{i})t^{2} * (3  2t) (12) If x_{i} ≈ x_{i+1}, then the cubic spline polynomial (12) becomes a straight line x(t) = αx_{i}. If the envelope becomes very symmetric and extrema of the same kind are close in value, then the cubic spline becomes very close to straight line segments x(t) ≈ x_{i} and this is resulting in fast varying component adjacent extrema (pair of a consecutive maxima and minima) to be selected as first invariants. In the general case there are N points and N1 piecewise splines. If we denote a piecewise spline that begins at point x_{i} as Y_{i}(t), we have N1 polynomials and need to determine {b_{i}, c_{i}, d_{i,}, k_{i}} for 0<=i<=N2 or 4(N1) unknowns. For this we will need 4(N1) equations. System of equations (10) is not valid for the last polynomial and we need more conditions imposed on the polynomials. These are obtained by requiring the polynomial’s continuity at joint points, the continuity of their first and second derivatives, as well as second derivatives being zero at the two data end points. This completes the specifications of the full system of linear equations for {a, b, c, d} for N1 cubic splines which has been proven to have a unique solution [10]. It is sufficient for our task to just know that this system has a solution and that the polynomials present a “natural” curve smoothly connecting the control data points. It is not obvious but the median M(t) of two piecewise cubic splines u(t) and d(t), as defined in Section 2 equation (3). is also a piecewise cubic spline. Indeed, if we set u(t) = a1 + b1t +c1t^{2} + d1t^{3} d(t) = a2 + b2t +c2t^{2} + d2t^{3} then M(t) = (u(t) – d(t))/2 = (a1 + a2)/2 + ((b1 + b2)/2)t + ((c1 + c2)/2)t^{2} + ((d1 + d2)/2)t^{3} and function M(t) is obviously contiguous and differentiable and its first derivative is also continues and differentiable with the resulting second derivative being contiguous. In other words, function M(t) is also a smooth function. In the proximity of an extrema point M(t) is close to the line of which the other end is the next opposite type extrema, making this pair an invariant for the following sifting steps. In the beginning of the EMD sifting process all the extremas obviously belong to the fastest varying component. Smooth connection of them by a piecewise cubic spline results in a piecewise cubic spline median which is also smooth and as a consequence the fast varying adjacent extrema pairs are near symmetric or quickly molded (being made) into being symmetric and thus become first invariant pairs in the sifting process. This results in the fastest varying component of an arbitrary signal s(t) to be sifted out first by the EMD sifting process. Finally, what is important is to revisit the concept of piecewise cubic spline smoothness by evaluating the polynomials’ global maximum G on [x_{i}, x_{i+1}] G = max(s{x_{i}, x_{i+1}, x(t_{1}^{0}), x(t_{2}^{0})}) (13) where t_{1}^{0}, t_{2}^{0} are the two roots of the cubic spline polynomial (9) first derivative b + 2ct + 3dt^{2} = 0 (14) It is important that the piecewise cubic splines u(t), d(t) control points are the data extrema points {x_{i} }. This ensures that G is very close to max(s_{1}, s_{2}, …,s_{N})=max{x_{i} }. It is also, for example, important to work with nearzero slopes for which G is close to max{x_{i}, x_{i+1}}, assuring nearlinear behavior of the cubic spline between two sparse extrema points of the same type. For example, if k_{i} – k_{i+1} = ε and x_{i+1} – x_{i} = ε, where ε is small Then c and d are also small and x(t) = a + bt G in (13) can also be evaluated using the law of the mean from Calculus [6], [11]. The first derivative of the cubic spline(14) is a parabola varying from b_{i} to 2c_{i}+3d_{i} (for 0 and 1 end parameter values at t=t_{i}). The second derivative of the cubic spline is 2c+6dt and is a straight line while it can only grow or only decrease on [x_{i}, x_{i+1}] and it is definitely differentiable on the open interval (x_{i}, x_{i+1}). The law of the mean states: Let x be a function of t which is continuous at each point of the closed interval A<=t<=B, and let it have a derivative at each point of the open interval A x(t_{i+1})  x(t_{i}) = (t_{i+1}  t_{i}) * x^{’}(θ) The point θ is such that the derivative at θ gives a slope equal to the slope of the straight line connecting x(t_{i+1}) and x(t_{i}) and which is the straight line around which the fastest varying component is locally weaving around on interval [x_{i}, x_{i+1}]. That the locally fastest varying component is being sifted out first in such a case has been already proven above. This concludes the initial proof of Hypothesis 1 that for an arbitrary signal the fastest scale component is sifted out first by the EMD sifting process, and as a consequence, the remaining sifting becomes computationally less difficult. ^ Other issues, such as the theory of local symmetry and IMF structure, the theory of the EMD rapid convergence, HHT signal synthesis, and ways of breaking up the signal for faster EMD processing, are left for a future paper in which we will prove the following Hypotheses 2 –3: Hypothesis 2: The EMD sifting process preserves an intermediate locally symmetric zerocrossing pair of extrema points with interleaved regions of diminishing amplitudes yielding an IMF with a definitive structure. Hypothesis 3: The EMD sifting process’ rapid convergence is of order O(1/2^{k}). This is a consequence of the EMD envelope control points definition as sets of extremas of the same type, its interpolation by piecewise cubic spline whose control points are the data extremas, and envelope’s median construction as an arithmetic median the envelops’ sum divided by 2. ConclusionsWe have reported the initial theoretical proof of why the fastest changing component of a composite signal is being sifted out first in the Empirical Mode Decomposition sifting process as implemented in the HHTDPS. We have also provided the two hypotheses for the theoretical explanation of why the EMD algorithm converges and converges rapidly while using cubic splines for signal envelope interpolation. As part of this research we developed a few practical techniques for cutting a large input data set into smaller files, which facilitates faster HHTDPS processing of large sound files. We have also developed in parallel with this investigation a few related to this research applications. For example, we have developed one of the first techniques for signal synthesis using the Hilbert Spectrum, by “painting” a recognizable 2D pattern in Hilbert Spectrum image and then tracing back (direct inversion) the painted subset pixels to their origin in the IMFs, and deleting them from the IMFs. This then allows reconstructing the portion of the input signal (its synthesis) from the modified IMFs that depicts the extracted feature being successfully removed. These developments will be presented in follow up papers. Future work includes research on the IMFs basis orthogonality, the affects of normalization on the Hilbert Transform and resulting instantaneous frequency, as well as further research in the 2D signal processing domain and handling of intermittency. AcknowledgementsThis work was performed by the authors as a team of NASA GSFC Applied Engineering and Technology Directorate (AETD) researchers and was funded by the AETD 2005 Research and Technology Development of Core Capabilities Grant (Advanced Communications R&TD 2005), as well as by the NASA Technology Transfer Accelerated Partnerships (TTAP) program. The assistance of Michael A. Johnson/AETD Code 560 Chief Technologist and Nona Cheeks/TTAP Code 504 Chief in sponsoring and encouraging this work is greatly appreciated and acknowledged. References[1] The empirical mode decomposition and the Hilbert spectrum for nonlinear and non stationary time series analyses by Norden E. Huang et al. Proc. R Soc. London. (1998) 454, 903995 [2] A new view of nonlinear water waves: The Hilbert Spectrum by Norden E. Huang at al Annual. Rev. Fluid Mech. 1999, 31 417457 [3] Norden E. Huang Patent #6,311,130 B1 Patent Date: Oct. 30, 2001 Patent Name: Computer Implemented Empirical Mode Decomposition Method, Apparatus, and Article of Manufacture for TwoDimensional Signals [4] Semion Kizhner, Thomas P. Flatley, Dr. Norden E. Huang, Karin Blank, Evette Conwell “On the HilbertHuang Transform Data Processing System Development Phase 1”, 2003 MAPLD Conference Washington D.C., September 911, 2003 [5] Semion Kizhner, Thomas P. Flatley, Dr. Norden E. Huang, Karin Blank, Evette Conwell “On the HilbertHuang Transform Data Processing System Development”, 2004 IEEE Aerospace Conference Proceedings, Big Sky Montana, March 613, 2004 [6] Angus E. Taylor, Advanced Calculus, Blaisdell Publishing Company, 1955, p.483 [7] Yakov Z. Tsypkin, “Sampling Systems Theory and its Applications Volume 1”, The Macmillan Company, New York 1964 [8] Robert W. Ramirez, “The FFT Fundamentals and Concepts”, PrenticeHall, 1985 [9] Arthur Berry, “A Short History of Astronomy”, Dover Publication 1960, par. 299 [10] Teukolsky, Vettering, Flannery “Numerical Recipes in C The art of Scientific Computing 2d Edition”, Cambridge University Press 1999 [11] Beyer, “CRC Standard Mathematical Tables 26^{th} Edition”, CRC Press, Inc. 
On the HilbertHuang Transform Data Processing System Development  An Introduction to HilbertHuang Transform: a plea for Adaptive Data Analysis Norden E. Huang. Research Center for Adaptive Data Analysis. National Central University Часть Сущность преобразования. Обработка и анализ данных Норден E. Хуанг. Исследовательский центр адаптивного анализа данных. Национальный Центральный Университет  
S. P. Shen Hilbert—Huang Transform and Its Applications Результаты emdhsa не имеют ложных гармоник (результатов наложения свойств линейности на нелинейные системы) и не ограничиваются...  Фон: • Экспертиза анализа Fourier ...  
Zhaohua Wu1 и Norden E. Huang Она применяется в широком диапазоне приложений для того, чтобы выделить сигналы из данных, сгенерированных в шумных нелинейных и...  Документы 1. /Theoretical Endings.doc  
Документы 1. /TRANSFORm.doc  Документы 1. /transform.doc  
Документы 1. /transform.doc  The results of theoretical research Vladimir Fedorovich Vlasov Science, Physics, Head of Research Department 