Abstract:
The damping behavior of a physical system is evaluated by mapping of damping loss factors calculated from time-dependent amplitude and frequency functions to form a full time-frequency dependent damping spectrum. The time-dependent amplitude and frequency functions are formulated by calculations based on the intrinsic mode functions derived by empirical mode decomposition of the original time series dataset representing the system behavior. The amplitude and frequency functions for each intrinsic mode function are calculated from the polar representation of the Hilbert transform-based time-dependent complex function corresponding to the intrinsic mode function.

Description:
The present invention relates in general to mathematical analysis of time series data with respect to damping and decay relaxation effects. 
     BACKGROUND OF THE INVENTION 
     Damping effects govern the behavior of a wide variety of physical systems. Extracting damping characteristics from measured data is a practice utilized in scientific and engineering applications. For example, structural damping effects on how a bridge reacts during an earthquake or when aircraft experience turbulence is important in order to enable creation of better and safer designs. Damping parameter extraction algorithms have been developed using a variety of techniques and continue to be of considerable scientific and engineering importance. The present invention is a contribution to this technology area. 
     An empirical study usually involves analysis of time series obtained through sensors. A method for analysis of such test data exploiting the use of the Hilbert Transform, is explained in U.S. Pat. No. 5,983,162 to Huang, issued Nov. 9, 1999. Through computer implemented Empirical Mode Decomposition (EMD), such method decomposes a given dataset into a sum of Intrinsic Mode Functions (IMFs). The Hilbert Transform employed in the usual way specifies instantaneous amplitude and frequency behavior as functions of time for each IMF. The foregoing procedure enables time-frequency dependent signal-amplitude analysis (Hilbert-Huang Spectrum). 
     However, more detailed physical characteristics of time series datasets may be needed than that afforded by signal-amplitude analysis, such as transient decay relaxation rates, or damping characteristics of resonant mode families. Heretofore, only somewhat limited conventional analyses have been available to examine the damping characteristics of a system. It is therefore an important object of the present invention to provide a method that fully utilizes the benefits of an improved analytical technique realized by use of the EMD and Hilbert Transform in order to determine damping characteristics. 
     SUMMARY OF THE INVENTION 
     In accordance with the present invention the methodology disclosed in the aforementioned U.S. patent to Huang involving analysis of time series data, is expanded into applications such as time-frequency analysis of structural dynamics and structural shock responses involving development of calculation algorithms on time and frequency damping loss factors through which the damping characteristics of a system are evaluated. 
     As a first step, a time dependent decay rate function is defined for each empirical mode (k). Based on empirical mode decomposition (EMD), a dimensionless quantity as a function of time (t) is then formulated, corresponding to the critical damping ratio, which is one-half the damping loss factor η k (t) at structural resonance. This damping loss factor is then evaluated at given time and instantaneous damped frequency utilizing the Hilbert transform methodology described in the Huang patent hereinbefore referred to. Algorithmic calculation is utilized to quantify and map time-frequency dependent damping loss factors under a full damping spectrum. A root mean square time average of the damping loss factors is then formulated as a function of frequency and averaging time (T) to provide an analysis of the system damping characteristics based on the time series data obtained from single time measurements. 
    
    
     BRIEF DESCRIPTION OF DRAWING 
     A more complete appreciation of the invention and many of its attendant advantages will be readily appreciated as the same becomes better understood by reference to the following detailed description when considered in connection with the accompanying drawing wherein: 
     The drawing FIGURE is a flow chart diagram of a computer based method for analyzing damping behavior of a system utilizing the Hilbert Damping Spectrum, pursuant to the present invention. 
    
    
     DESCRIPTION OF PREFERRED EMBODIMENTS 
     Referring now to the drawing, it diagrams the algorithmic method of the present invention initiated by input of time series data  10  from measurements of a physical system being analyzed. Under a step  12  derived from the methodology disclosed in U.S. Pat. No. 5,983,162 to Huang, Empirical Mode Decomposition (EMD) is performed to decompose the time series data into a set of simple oscillatory functions denoted as  14 , defined as Intrinsic Mode functions (IMF) corresponding to a total of (n) intrinsic modes (k), n being sufficiently large that no further IMFs result from further EMD. Such (IMF) functions  14  must satisfy the conditions that the number of extrema must be equal to the number of zero crossings (or differ by one at most). Such components of the (IMF) functions  14  are obtained by a repeated application of an interactive procedure referred to as “sifting”. The original dataset is recovered by summing all of the intrinsic mode functions in accordance with the following equation:                  X        (   t   )       =         ∑     k   =   1     n            c   k          (   t   )         +     r   n         ,           (   A   )                                
     where c k (t) is the k th  IMF. X(t) is measured timed series data and r n  is a residual term. 
     The next step  16  of the diagrammed method involves use of a Hilbert transform to define time-dependent amplitudes a k (t) and frequencies ω k (t) for each intrinsic mode (k). For a given intrinsic mode function c k (t), there is a corresponding Hilbert transform d k (t) determined by the following equation:                    d   k          (   t   )       =       1   π        P          ∫     -   ∞     ∞                c   k          (   t   )         t   -   t               t             ,           (   B   )                                
     where P denotes the Cauchy principal value. A complex signal Z k (t) is defined by use of the k th  IMF function component c k (t) and its Hilbert Transform d k (t). Using such components, the physical motion is expressed in terms of amplitude and frequency as shown in the following equations:                    z   k          (   t   )       =         c   k          (   t   )       +       id   k          (   t   )           ;           (   C   )                     a   k          (   t   )       =             c   k          (   t   )       2     +         d   k          (   t   )       2           ;           (   D   )                     θ   k          (   t   )       =       tan     -   1                d   k          (   t   )           c   k          (   t   )             ;   and           (   E   )                     ω   k          (   t   )       =                         θ   k          (   t   )                         t         ,           (   F   )                                
     where a k (t) is an instantaneous amplitude defined by the complex mode function z k (t), θ k (t) is a monotone instantaneous phase function, and ω k (t) is the instantaneous frequency. Systematic extraction of the intrinsic mode functions by computer implemented EMD or sifting as hereinbefore referred to is performed as described in the Huang patent in order to derive the instantaneous phase and frequency functions and complete the diagrammed Hilbert Transform step  16  for each of the modes k=1, 2, 3 . . . n. 
     Pursuant to the present invention, the Hilbert Transform step  16  is performed so that the frequency function ω 0k (t) is defined by the following equation:                    ω     0      k            (   t   )       =             ω   k          (   t   )       2     +       (         γ   k          (   t   )       2     )     2           ,           (   G   )                                
     in which γ k (t) is defined by the following equation:                  γ   k          (   t   )       =           -   2                                         a   k          (   t   )              t             a   k          (   t   )         .             (   H   )                                
     Using such amplitude and frequency functions a k (t) and ω k (t) obtained upon completion of step  16 , the next step  18  is performed to calculate damping loss factor η k (t) under a Hilbert Damping Spectrum. Damping is thereby described for the k-th mode through step  18  by use of the following equation:                  η   k          (   t   )       =         -   2                                         a   k          (   t   )                         t               a   k          (   t   )                  [       ω   k          (   t   )       ]     2     +       [       1       a   k          (   t   )         ·                         a   k          (   t   )                         t         ]     2                     (   I   )                                
     The damping loss factor η k (t) is expressed independent of mode index (k) by performance of the next step  20  mapping the damping loss factor, Equation (I), for each mode on a time-frequency domain, whereby the loss factor η(ω,t) is extracted for all modes by expression in terms of both frequency (ω) and time (t). Step  20  accordingly involves representation of the damping loss factor η k (t) as a joint function of time and frequency in three dimensional (3D) space. When such 3D space mapping  20  is performed, θ k (t) is replaced by η(ω,t) for performance of a following step  22  pursuant to one embodiment of the invention resulting in the full damping spectrum reflected by the following equation:                  η        (     ω   ,   t     )       2     =           [         -   2       a        (     ω   ,   t     )         ·                       a        (     ω   ,   t     )                         t         ]     2         ω   2     +       [                         a        (     ω   ,   t     )                         t         a        (     ω   ,   t     )         ]     2         .             (   J   )                                
     As an alternative or additional embodiment, the damping loss factor output η(ω,t) of mapping step  20  may be utilized to create a loss factor (η) as a function of only one of either frequency (ω) or time (t) through steps  24  and  26 . Thus, a set of sample points is taken over a total sample time (T), a root mean square may be chosen to calculate the damping loss factor in accordance with the following two equations:                η        (     ω   ,   T     )       =         1   T            ∫   0   T              η        (     ω   ,   t     )       2             t                     (   K   )                                
     As T approaches infinity, the foregoing equation is isolated for simplification of loss factor calculation as a function of frequency (ω) under step  24  as follows:                η        (   ω   )       =       lim     T   -&gt;   ∞            η        (     ω   ,   T     )                 (   L   )                                
     If one chooses a selected frequency band B=(ω 1 , ω 2 ) to define a time dependent damping behavior, the damping loss factor as a function of time η(t,B) may be calculated under step  26  as follows:                  η        (     t   ,   B     )       =         1        B                 ∫     ω   1       ω   2                η        (     ω   ,   T     )       2             ω               ,           (   M   )                                
     where |B|=|ω 1 −ω 2 | is the chosen bandwidth. 
     Obviously, other modifications and variations of the present invention may be possible in light of the foregoing teachings. It is therefore to be understood that within the scope of the appended claims the invention may be practiced otherwise than as specifically described.