## Abstract

We introduce a new algorithm, the intrinsic time-scale decomposition (ITD), for efficient and precise time–frequency–energy (TFE) analysis of signals. The ITD method overcomes many of the limitations of both classical (e.g. Fourier transform or wavelet transform based) and more recent (empirical mode decomposition based) approaches to TFE analysis of signals that are nonlinear and/or non-stationary in nature. The ITD method decomposes a signal into (i) a sum of proper rotation components, for which instantaneous frequency and amplitude are well defined, and (ii) a monotonic trend. The decomposition preserves precise temporal information regarding signal critical points and riding waves, with a temporal resolution equal to the time-scale of extrema occurrence in the input signal. We also demonstrate how the ITD enables application of single-wave analysis and how this, in turn, leads to a powerful new class of real-time signal filters, which extract and utilize the inherent instantaneous amplitude and frequency/phase information in combination with other relevant morphological features.

## 1. Introduction

Signal and other data analyses are an essential part of pure and applied research in all disciplines. Real-world data are often non-stationary, noisy and obtained from underlying nonlinear systems or processes about which dynamical information is unknown or limited. However, the majority of tools available for quantitative analysis of such signals are better suited for the analysis of stationary signals, with little or no noise, obtained from linear systems.

The evolution of methods for time–frequency–energy (TFE) analysis of complex signals of arbitrary origin has been guided by attempts to:

carefully balance the trade-offs between precise temporal information and precise frequency information, obtaining temporal information on the time-scales determined by the frequencies present in the signal itself,

avoid

*a priori*assumptions about the content/morphology of the signal being analysed (e.g. make the decomposition ‘basis free’), andperform the analysis in an efficient and rapid manner (ideally, online and in real-time).

The Fourier transform has long been the dominant TFE analysis tool for all types of signals and data, despite its inherent assumption of signal stationarity. Attempts to address this limitation in analysis of non-stationary signals typically involve the application of the Fourier transform in moving windows of demeaned and/or tapered data or the application of linear, finite (FIR) or infinite impulse response filters in order to extract the desired TFE information (e.g. Oppenheim & Schafer 1989). However, it is well known that these approaches force a difficult choice between: (i) good frequency resolution in the TFE information (by using long windows/filters) at the expense of temporal resolution, or (ii) good temporal resolution (by using short windows/filters) at the expense of frequency resolution. In addition, although efficient, the fast Fourier transform uses a rectangular partitioning of the time–frequency plane on which the information is obtained, so that the temporal information regarding high-frequency fluctuations in the signal is blurred over the entire window. Additional drawbacks of Fourier-based analysis include the introduction of spurious harmonics and the necessity to use a wide frequency spectrum in order to represent brief transient changes (e.g. impulses or steps) or other data that are non-stationary in time.

Over the past couple of decades, alternative approaches for obtaining TFE information based upon the wavelet transform have become popular (e.g. Strang 1989; Haar 1910; Daubechies 1992; Strang & Nguyen 1996), owing to this transform's dyadic partitioning of the time–frequency plane, which allows a better temporal localization of high-frequency signal changes and better frequency information for low-frequency changes. The primary strength of wavelets over Fourier analysis lies in their ability, via dilation and translation of basis elements with compact support, to provide TFE information that more closely match the time-scales at which the quantified signal changes occur. Figure 1*a–d* illustrates the performance and some limitations of these popular approaches in extracting TFE information from a sample brain signal, which serves as model real-world signal that may exhibit nonlinear, non-stationary, high-dimensional and noisy characteristics.

A drawback inherent to both the Fourier and the wavelet analyses as well as most other approaches to TFE analysis, such as those based on the Wigner–Ville distribution (e.g. Cohen 1995) or the evolutionary spectrum (Priestly 1965), is that these methods represent the signal in terms of basis functions that are defined *a priori* and the performance of the analysis method is closely tied to how well the morphology of signals under study is represented by the selected basis elements. For example, sound waves produced by an orchestra and consisting of superimposed pure tones with characteristic harmonics and fixed durations (quarter notes, sixteenth notes, etc.) are well represented via Fourier analysis (e.g. sheet music), whereas more complex non-stationary signals, such as ambient room noise, meteorologic signals, financial signals or brain waves, make *a priori* construction of suitable basis elements difficult, if not impossible, and also preclude accurate temporal segmentation with prespecified windows, as the precise times and time-scales of signal changes are typically unknown in advance.

An alternative approach to TFE analysis that made progress towards addressing these difficulties is known as the empirical mode decomposition (EMD; Huang *et al*. 1998, 2003). The EMD is an algorithm for decomposing physical signals, including those that may be nonlinear or non-stationary, into a collection of intrinsic mode functions (IMFs), which the method's developers claim are indicative of intrinsic oscillatory modes in the physical signals. The EMD method obtains this decomposition by defining a baseline trend signal for a given window of an input signal, such that when this baseline is subtracted from the input, the resulting residual signal (which they call an IMF) is representative of the highest relative frequency riding waves present in the input. The baseline is constructed as the mean value of the upper and lower envelopes of the input signal, with these upper and lower envelopes defined by cubic spline interpolation of the input signal's local maxima and local minima, respectively.

The EMD procedure for obtaining the baseline is intended to produce a residual IMF that is well behaved, with well-defined instantaneous amplitude and frequency. More precisely, the EMD's objective is to produce an IMF that is a ‘proper rotation’, which is a signal that has strictly positive values at all local maxima and strictly negative values at all local minima. Once the first residual IMF function has been obtained, it is subtracted from the input signal and the process is repeated on the resulting lower frequency baseline signal. This process is iterated, producing a multilevel decomposition of the original or raw input signal into IMFs of successively lower relative frequency, until all the riding waves are removed and the decomposition is completed, leaving a trend signal. The Hilbert transform can then be applied to the IMF components to extract instantaneous amplitude and frequency information.

While the introduction of the EMD constitutes a conceptual advance in TFE analysis for non-stationary and nonlinear signals, it has several practical limitations that reduce its practical utility and may cause inaccuracies in depicting signal dynamics. The EMD-defined baseline often fails to produce a residual that has the important property of being a proper rotation. In an effort to work around this inability to consistently generate proper rotation components, Huang *et al*. (1998) developed a ‘sifting’ process that is applied in an iterative manner to generate a sequence of signal baseline candidates until one is found with a proper rotation residual. This sifting process has two stated goals: (i) to separate out high-frequency, small-amplitude waves, which are ‘riding’ atop, or superimposed on, larger amplitude, lower frequency waves, and (ii) to smooth out uneven amplitudes in the IMF being extracted. However, these goals are often conflicting for non-stationary signals, since riding waves may be transient in nature and/or highly variable in amplitude, and smoothing out the uneven amplitudes via sifting can prevent faithful extraction of these waves. Moreover, repetitive sifting causes smearing of TFE information across different decomposition levels and an intra-level smoothing of TFE information, which is unlikely to reflect the intrinsic characteristics of the signal under analysis. Despite this effort, the sifting process still often fails to produce proper rotations and thus requires a stopping criterion1 in the EMD procedure to keep the process bounded and limit the amount of computational energy expended.

The requisite sifting also forces the EMD to be applied in a window-based manner, with data-dependent (and thus uncertain) computational expense and loss of TFE information on time-scales longer than the window length. Each sifting then causes an inward propagation of window boundary or edge effects, hindering its ability to extract meaningful TFE information. In addition, the fact that the EMD's baseline is exclusively determined by extrema, and ignores all other critical signal information, contributes to mislocalization of temporal information and produces spurious phase shifts and distortion in the components it extracts. Its use of splines to obtain the baseline also causes the EMD to inherit all the well-known difficulties associated with this form of interpolation. For example, overshoots and undershoots of the interpolating cubic splines generate spurious extrema, and may shift or exaggerate the existing ones. This results in the EMD's production of IMFs with critical points, which typically differ from those of the original signal and do not retain the important TFE information that these points naturally possess.

In this paper, we present a new method for TFE analysis, called the intrinsic time-scale decomposition (ITD), which is specifically formulated for application to nonlinear or non-stationary signals of arbitrary origin and obtained from complex systems with underlying dynamics that change on multiple time-scales simultaneously. The ITD overcomes the limitations of EMD listed earlier, as well as those previously mentioned and associated with more classical approaches such as Fourier and wavelets. In particular, the ITD provides:

efficient signal decomposition into ‘proper rotation’ components, for which instantaneous frequency and amplitude are well defined, along with the underlying monotonic signal trend, without the need for laborious and ineffective sifting or splines,

precise temporal information regarding instantaneous frequency and amplitude of component signals with a temporal resolution equal to the time-scale of occurrence of extrema in the input signal, and

a

*new class of real-time signal filters*that utilize the newly available instantaneous amplitude and frequency/phase information together with additional features and morphology information obtained via single-wave analysis. Moreover, the resulting feature extraction and feature-based filtering can be performed in real-time and is easily adapted to extract feature signals of interest at the time-scales on which they naturally occur, while preserving their morphology and relative phases.

## 2. The intrinsic time-scale decomposition

Given a signal, *X*_{t}, we define an operator, , which extracts a baseline signal from *X*_{t} in a manner that causes the residual to be a proper rotation. More specifically, *X*_{t} can be decomposed as(2.1)where is the baseline signal and is a proper rotation.

Suppose {*X*_{t}, *t*≥0} is a real-valued signal, and let {*τ*_{k}, *k*=1, 2, …} denote the local extrema of *X*_{t}, and for convenience define *τ*_{0}=0. In the case of intervals on which *X*_{t} is constant, but which contain extrema due to neighbouring signal fluctuations, *τ*_{k} is chosen as the right endpoint of the interval. To simplify notation, let *X*_{k} and *L*_{k} denote *X*(*τ*_{k}) and *L*(*τ*_{k}), respectively.

Suppose that *L*_{t} and *H*_{t} have been defined on [0, *τ*_{k}] and that *X*_{t} is available for *t*∈[0, *τ*_{k+2}]. We can then define a (piece-wise linear) baseline-extracting operator, , on the interval (*τ*_{k}, *τ*_{k+1}] between successive extrema as follows:(2.2)where(2.3)and 0<*α*<1 is typically fixed with *α*=1/2. We construct the baseline signal, *L*_{t}, in this manner in order to maintain the monotonicity of *X*_{t} between extrema, while at the same time remaining inside an ‘envelope’ generated by some wave riding atop this baseline. The extrema are interpreted as evidence of some proper rotation, riding wave to be extracted. The baseline is constructed as a linearly transformed contraction of the original signal in order to make the residual function monotonic between extrema, a necessity for proper rotations. The approach also enables information ‘intrinsic’ to the original signal to be passed down to the baseline and residual components. We found that other attempts at baseline construction that were not based on use of the input signal, inevitably failed to produce proper rotation residuals.

After defining the baseline signal according to equations (2.2) and (2.3), we are able to define the residual, proper-rotation-extracting operator, , as(2.4)

Since

*X*_{t}is monotonic on (*τ*_{k},*τ*_{k+1}], and*L*_{t}and*H*_{t}are linear transformations of*X*_{t}on this interval, it follows that*L*_{t}and*H*_{t}are also monotonic on (*τ*_{k},*τ*_{k+1}]. The extrema of*H*_{t}, i.e. {*τ*_{k},*k*≥1}, coincide with extrema of the input signal,*X*_{t}. The extrema of*X*_{t}are either extrema for*L*_{t}or inflection points.The manner of defining these components of

*X*_{t}on inter-extrema intervals allows the operation to be performed in a recursive manner, interval by interval, with the arrival of each new local extrema. The process is illustrated in figure 2 and allows the method to operate very efficiently and in real-time. Thus, the method is applicable to both online or offline signal analysis.The components are defined at the same time points as the original signal, whether

*X*_{t}is a function of a continuous variable,*t*, or a discrete (sampled) variable, and the method*does not*require the raw data to be evenly sampled in time. We have attempted to present the description using a notation suitable for either a discrete- or a continuous-time framework.Although

*α*is most naturally and typically chosen to be 0.5 (as in figure 2), any choice of*α*in the interval (0, 1) produces a proper rotation component. As we will show, the parameter*α*acts as a gain control that linearly scales the amplitude of the proper rotation component extracted at each application.It is important to note that the proper rotation components produced by the ITD method are not unique and the decomposition is nonlinear in the sense that the decomposition of a sum of two signals need not produce components equal to the sum of the components obtained from decomposing each signal individually.

The extracted baseline signal,

*L*_{t}, has the same smoothness/differentiability as the input signal,*X*_{t}, at values of*t*between extrema. At the extrema, successive baseline components will be continuous and differentiable, but in general will not be twice differentiable. This introduction of non-smoothness in successive baselines must be kept in mind in interpreting instantaneous frequency measurements obtained from analysis of the resulting components.To initialize the decomposition on the interval [0,

*τ*_{1}], one may consider the first point of the signal as an extremum (i.e. define*τ*_{0}=0), and define . Other alternatives, such as those that extend the function to values of*t*<0 via reflections, can also be used. The effect of this initialization on the resulting decomposition is confined to the interval [0,*τ*_{2}]. If the ITD is used to analyse a fixed window of data, a similar procedure may be used to construct the decomposition beyond the last measured extrema in the interval, should this be desired. In practice, if window-based analysis is used, then resulting TFE information obtained from these beginning and ending ‘partial waves’ should be interpreted with caution.

*H*_{t} *is a proper rotation on the interval* [*τ*_{1}, *τ*_{N}] *for any N*≥1.

Without loss of generality, assume that *X*_{t} has a local maximum at *τ*_{k}, and adjacent local minima at *τ*_{k−1} and *τ*_{k+1}. It suffices to show that *H*_{t} has a local maximum at *τ*_{k} and is monotonic on [*τ*_{k−1},*τ*_{k}] and on [*τ*_{k},*τ*_{k+1}]. Now(2.5)where , ∀*k*.

*X*_{t} monotonic on [*τ*_{k−1}, *τ*_{k}] and [*τ*_{k}, *τ*_{k+1}]⇒*H*_{t} monotonic on [*τ*_{k−1}, *τ*_{k}] and [*τ*_{k},*τ*_{k+1}].

By construction, *L*_{k−1}>*X*_{k−1}, *L*_{k+1}>*X*_{k+1} and *L*_{k}<*X*_{k}⇒*H*_{k−1}=*X*_{k−1}−*L*_{k−1}<0, *H*_{k}=*X*_{k}−*L*_{k}>0 and *H*_{k+1}=*X*_{k+1}−*L*_{k+1}<0. It follows that *H*_{t} has a local maximum at *t*=*τ*_{k}. ▪

Once the input signal has been decomposed into a baseline and a proper rotation component, with the latter representing the highest relative frequency ‘riding wave’ present in the input signal, the process can be re-applied using the baseline signals as input. This procedure can be iterated until a monotonic baseline signal (i.e. the ‘trend’) is obtained. This decomposes the raw signal into a sequence of proper rotations of successively decreasing instantaneous frequency at each subsequent level of the decomposition (figure 3). More precisely,(2.6)where *X*_{t} is the (*k*+1)st level proper rotation and *X*_{t} is either the monotonic trend or the lowest frequency baseline extracted if the decomposition is stopped before the monotonic trend is reached. Here, we have made use of the fact that (by definition) and used the telescoping sum formula .

### (a) Instantaneous time–frequency–energy information

Once the ITD method has been applied to decompose the input signal into a set of proper rotation components and a residual signal, one may now analyse the proper rotations to extract instantaneous amplitude, phase and frequency information. The conventional approach to obtain instantaneous amplitude, *A*_{t}, instantaneous phase, *θ*_{t} and instantaneous frequency, *f*_{t}, information from a proper rotation signal, *R*_{t}, is based on the application of the Hilbert transform, *h*[.], as follows:(2.7)where , and ‘P.V.’ denotes the Cauchy principal value of the integral. However, in practice, this Hilbert transform-based approach to instantaneous TFE analysis has some difficulties for rapid online analysis that stems from its edge effects, non-recursive nature and the occasional appearance of negative frequencies, which should not occur for proper rotations. Therefore, we introduce two useful alternatives for determining instantaneous amplitude, phase and frequency that do not use the Hilbert transform. These alternatives are ‘wave based’ (i.e. they define the instantaneous TFE information in a piece-wise manner, on each time-interval between successive up-crossings of a proper rotation, and based only upon information about the single wave of the proper rotation occurring during that period) and, unlike the Hilbert transform, guarantee a monotonically increasing phase angle when applied to proper rotation components. The two approaches obtain their respective phase angles, and , from the signal *x*_{t} for one full wave (i.e. the portion of the signal between successive zero up-crossings) at a time, as follows:(2.8)and(2.9)where *A*_{1}>0 and *A*_{2}>0 are the respective amplitudes of the positive and negative half-waves between successive zero up-crossings, *t*_{1} and *t*_{5}, and where *t*_{2} is the (earliest) time of the maxima (*A*_{1}) on the positive half-wave; *t*_{3} is the zero down-crossing time and *t*_{4} is the (earliest) time of the minima (−*A*_{2}) on the negative half-wave. The second approach is similar to the first, but provides an approximation alternative that is more computationally efficient (i.e. avoiding arcsine evaluation) by linearly interpolating phase between zero crossings and extrema. In either approach, the instantaneous amplitude is piece-wise constant and determined by the extrema values of the proper rotations between zero crossings, as follows:(2.10)

These relatively simple definitions ensure that the phase angle is 0 at every zero up-crossing,

*π*/2 at the local maxima of the wave,*π*at each zero down-crossing and 3*π*/2 at each local minima.The fact that the set of proper rotation components produced ITD method is not unique implies that the same holds for the resulting TFE information obtained from this approach. As with any TFE analysis of non-stationary signals, one must take care in interpreting the derived information.

Incomplete waves at the beginning and/or end of a proper rotation component, i.e. data prior to the first zero up-crossing or after the last zero up-crossing, may be treated according to these definitions over the appropriate sub-interval of [

*t*_{1},*t*_{5}].In online applications of the ITD, the evolution of the phase angle over any monotonic segment of a proper rotation can be computed between times of successive extrema as soon as the segment has been obtained from the ITD decomposition. This is a significant advantage over the Hilbert transform-based approaches, which require longer windows of data and not suited for application to single waves or short oscillations.

The smoothness/differentiability of the instantaneous phase signal, which directly impacts the extraction of instantaneous frequency information, clearly coincides with the smoothness of the proper rotation component signal,

*x*_{t}, between extrema.Using the notation to explicitly denote the dependence of

*H*_{t}in equation (2.4) on the parameter*α*, it is straightforward to derive the identity . This implies that the parameter*α*simply acts as a linear gain on the amplitude of the proper rotation component, and that the instantaneous frequency of the proper rotation is independent of the choice of*α*. Using a value of*α*near 0 extracts a very small proper rotation from the input signal, leaving more activity of the same frequency in the residual to be extracted in subsequent iterations. Alternatively, larger choices of*α*(i.e. above 0.5) decrease the number of components in the resulting decomposition, since more of the highest frequency activity is extracted at each iteration. Therefore, while the choice does not affect the frequencies of a proper rotation obtained in a single step of the decomposition, varying*α*will change the total, iterative decomposition of the original signal and thus it alters the representation of TFE information obtained.

Figure 4 illustrates the instantaneous frequencies computed for each of the ITD-generated proper rotation components shown in figure 3, through decomposition of the signal of figure 1*a*. Figure 4*a* shows the instantaneous frequencies obtained using the Hilbert transform, while figure 4*b* shows the corresponding instantaneous frequencies computed using the alternative phase angle definition, , given in equation (2.8). In both the plots, instantaneous frequency was obtained by differentiating the instantaneous phase angles using an 11 coefficient least-squares FIR differentiating filter (Frei *et al*. 1999). The instantaneous frequencies obtained from are clearly stabler, whereas the Hilbert transform-based output is prone to large spikes. These spikes are due to 2*π* jumps in the Hilbert transform generated phase angles, which occur during the unwrapping procedure when imposing the constraint that the phase angles be monotonically increasing, a necessity since the phases are derived from proper rotation signals. While removing this constraint eliminates these jumps, it is not an improvement since the phase angles are then allowed to (improperly) decrease at certain points, resulting in unwanted negative frequencies. These problems are eliminated using the alternative phase angle definitions we propose. The instantaneous amplitudes *A*_{1} and *A*_{2}, for each half-wave of each proper rotation component are shown in figure 4*c* and are similar to that obtained via the Hilbert transform (not shown). Figure 4*d* presents a revised rendering of the instantaneous frequency information shown in figure 4*b*, but only colours the instantaneous frequencies that correspond to half-waves with amplitude exceeding 0.1 mV, leaving the others coloured grey. This practice aids in visualization of the most significant waves present and provides an alternative view of TFE information to that shown in figure 1*e*.

## 3. Time–frequency–energy analysis of non-stationary signals

Once an input signal has been decomposed into proper rotation components and a low-frequency baseline trend signal, various morphologic and other features of waves or segments of the proper rotations can be analysed and used to reconstruct novel and useful ‘filtered’ signals. Given a proper rotation signal, we shall use the term ‘half-wave’ to indicate a segment of the signal between successive zero crossings and ‘monotonic segment’ to indicate a segment of the signal between successive extrema. The ITD method produces proper rotations, one monotonic segment at a time, and their precise duration, starting and ending values are automatically available as they are constructed. These features alone can be used to measure wavespeed (i.e. the inverse of twice the half-wave duration) and (signed) amplitude of individual half-waves (figure 3). More detailed morphological features can be easily extracted (with some additional computational expense) by examining such characteristics as moments (e.g. skewness and kurtosis) or deviation from a ‘pure tone’ (i.e. a portion of a sine wave with constant frequency equal to the half-wave or segment's wavespeed). Additional features, such as those that may be derived from a wave's latency from a fiducial time marker or its relation to corresponding waves on other input channels being analysed, can also be used in constructing useful filtered components.

Figure 1*e* shows an ITD-based TFE distribution for the signal of figure 1*a*, illustrating the significant improvement in time–frequency localization of energy provided by the ITD as compared to Fourier and wavelet analyses for this non-stationary signal. The displayed line segments' start and end points correspond *exactly* to the start and end times (abscissa) and wavespeeds (ordinate) of corresponding half-waves in ITD-produced proper rotation components. The colouring of each line segment is based upon the amplitude of the corresponding wave (red indicates greater amplitude).

Decomposition of brain wave signals into seizure and non-seizure components.

The ITD is ideally suited for the analysis of non-stationary biological signals such as those generated by the human brain cortex. A particular application of interest involves use of the ITD as a filter to decompose the electrical brain signals of a subject with epilepsy into two components, which may be conceptualized as a ‘seizure component’ and a ‘non-seizure component’ (Osorio *et al*. 1998). Figure 5 shows how the method can be used to accomplish this goal.

The raw signal (figure 5*a*) consists of 10 s of electrocorticogram data containing a seizure that begins 5 s into the recording. After the ITD is applied and the TFE distribution is computed, the absolute amplitude and wavespeed of each individual half-wave are obtained. Figure 5*b*,*c* shows the one-dimensional (marginal) densities of these two features along with discriminators (solid vertical lines) that were selected for purposes of illustrating the capabilities of this method to separate seizure from non-seizure signal activity. Figure 5*d* shows the (absolute amplitude and wavespeed) ordered pairs for each half wave. Two-dimensional feature discriminators are applied, separating pairs with wavespeed less than 20 Hz, between 20 and 50 Hz and greater than 50 Hz, and amplitudes above and below 0.01 mV. Feature vector pairs (i.e. points in figure 5*d*) falling into certain specific regions of the two-dimensional feature range are classified (and coloured) according to the region to which they belong. Figure 5*e* displays three filtered components constructed using only those segments with points in specific regions (corresponding to three of the regions in figure 5*d*). Component one, C_{1}, is the output of a filter that contains the signal reconstructed using only those waves that have absolute amplitude exceeding 0.01 mV and wavespeed between 20 and 50 Hz. Component two, C_{2}, is the output of a filter that contains the signal reconstructed using only those waves that have absolute amplitude less than 0.01 mV and wavespeed between 20 and 50 Hz. Component three, C_{3}, is the output of a filter that contains the signal reconstructed using only those waves that have absolute amplitude exceeding 0.01 mV and wavespeed below 20 Hz. A portion of the output obtained by combining some of these components (i.e. C_{1}+C_{2}+C_{3} and C_{1}+C_{2}) are shown in figure 5*f*, along with the output of a traditional, linear, [20, 50] Hz band-pass digital filter provided for comparison. The extrema-preserving properties of this new type of filter are evident upon inspection of figure 5*f*, and even more clearly illustrated through inspection of the residual signals obtained by subtracting the filter outputs from the raw signal (figure 5*g*). The improved separation between abnormal and normal signal components afforded by the ITD-based filter, as compared to the Fourier-based linear band-pass filter is readily apparent.

In practice, ITD-based filters may be designed by specifying *a priori* feature constraints (as in example 3.1) or can be trained by using available signal information and associated, ITD-derived, features. For example, one may compare the multidimensional probability densities of several features obtained from available recordings for one signal state (e.g. non-seizure) to those obtained for another state (e.g. seizure) and use maximum-likelihood discriminators to produce a feature-vector look-up table that can be updated with time and prospectively applied to decompose the incoming raw signal into components representative of the respective states and, if desired, a residual which differs from the past waveform characteristics observed in either state. This concept enables efficient real-time adaptive filtering useful for state change detection, as well as noise reduction or removal.

Decomposition of a Lorenz signal by ITD-filtering.

In this example, we illustrate the use of the ITD together with single-wave analysis in constructing wave-based feature-filtered signal components. Figure 6*a* illustrates 15 s of the first coordinate trajectory obtained by a simulation of the well-known nonlinear chaotic Lorenz system. Figure 6*b* shows the wavespeed and signed amplitude features of each ITD-generated half-wave, plotted as ordered pairs of feature values. Analysis of these two features, together with half-wave kurtosis (to differentiate between points that overlap in the plot) allows identification of distinct clusters and corresponding half-wave classification (each plotted in different colours) using these features. Reconstruction of components corresponding to the displayed feature clusters is shown in figure 6*c*, along with the original signal (top), illustrating the filter's capability to extract components consisting of half-waves in the original signal of similar morphology while precisely preserving their temporal position. Note that, for example, the sum of the first two components (green and red) may be used to count impending transitions between the two lobes of the chaotic attractor.

## 4. Comparison between the intrinsic time-scale and empirical mode decomposition methods

The ITD method overcomes essentially all the limitations of EMD described earlier. In particular, the ITD only involves *O*(*n*) computations, operates in real-time, and although it can be applied to windows of data, it naturally operates on inter-extrema segments of data, with a time-delay equal to the inter-extrema interval, which is determined by the frequency content of the signal being decomposed. When applied to a window of data, the unavoidable edge effects that stem from uncertainty regarding wave morphology in the intervals before the first and after the last signal extrema do not propagate inward as they can with the EMD, but rather are temporally constrained to margins defined by the first and last two extrema in the window. The ITD does not use splines and thus avoids the ‘difficulties’ they introduce. By construction, the ITD method guarantees proper rotation components2 and precisely preserves temporal information regarding TFE information and all the signal critical points. Finally, and perhaps most importantly, the sifting procedure introduced by Huang *et al*. (1998) in an attempt to produce proper rotation components with the EMD, but which also causes smearing and smoothing of TFE information, is not needed in the ITD. However, for heuristic purposes, we investigated its application within the ITD framework and determined that, as with the EMD, it has the general effect of smoothing variations in amplitude between adjacent extrema. Exhaustive sifting within a window of data results in a component that is purely frequency modulated with constant amplitude. Unlike the EMD, after each ITD sifting iteration, the extrema points remain identical to those of the input signal.

Figure 7 presents a comparison between ITD and EMD methods in decomposing (i) a non-stationary test signal (figure 7*a*,*b*), and (ii) a brain wave signal containing an epileptic seizure (figure 7*c–e*). In both examples of figure 7, the IMF components produced by the EMD method failed to be proper rotations. Instead, the sifting process was terminated when the extracted IMFs satisfied the stopping criteria prescribed in Huang *et al*. (1998). Removing this stopping rule and applying the EMD to the brain wave signal of figure 7*c* required 56 393 siftings before the *very first* resulting IMF was produced that satisfied the definition of a proper rotation. Moreover, *each* such sifting took more computations/processing time than the ITD-produced proper rotation, and the latter was obtained in only one sifting-free step. The decomposition of the test signal in figure 7*a* illustrates the detrimental edge effects of the EMD along with their inward propagation with each successive component and the improved accuracy in trend extraction afforded by the ITD. The lower frequency components extracted by the EMD also contain oscillations that do not seem to be present in the raw signal, stemming from the step-type drop and rise in the non-stationary signal. In the brain wave signal decomposition, the signal changes associated with the epileptic seizure (between 300 and 350 s into the segment) induce propagating oscillations in the EMD components that precede and follow the seizure changes. The trend signals (magnified for ease of inspection) illustrate the ITD's ability to detect small signal changes in baseline during the seizure, which go undetected by the EMD method.

## 5. Conclusion

The ITD method represents a significant advance in TFE analysis and nonlinear filtering that is especially useful for analysis of non-stationary or nonlinear signals. It improves upon both conventional methods of TFE analysis (e.g. Fourier and wavelet-based approaches) and the EMD, in accuracy and efficiency, accomplishing in a single step that which the EMD may require numerous iterations to achieve. The ITD method provides a powerful method of real-time signal analysis and decomposition/filtering well suited for applications that involve non-stationary signals.

## Acknowledgments

This research was supported in part by NIH/NINDS grants nos. 1R01NS046602-01 and 1R43NS39240-01.

## Footnotes

↵These criteria may be, for example, (i) when the resulting IMF becomes a proper rotation, (ii) when a Cauchy criteria is satisfied, or (iii) when some predetermined number of iterations has been exhausted without obtaining a proper rotation.

↵Numerically, the resulting HF components are proper rotations up to machine precision, beyond which a computer is unable to differentiate differences in extrema.

- Received December 2, 2005.
- Accepted July 11, 2006.

- © 2006 The Royal Society