Journal of Probability and Statistics

Volume 2013 (2013), Article ID 569597, 15 pages

http://dx.doi.org/10.1155/2013/569597

## Using Time Deformation to Filter Nonstationary Time Series with Multiple Time-Frequency Structures

^{1}Statistics Genomics Unit, NIH/NIMH (National Institutes of Health/National Institute of Mental Health), Bethesda, MD 20892, USA^{2}Department of Statistical Science, Southern Methodist University, Dallas, TX 75205, USA

Received 13 August 2012; Accepted 6 February 2013

Academic Editor: Zhidong Bai

Copyright © 2013 Mengyuan Xu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

For nonstationary time series consisting of multiple time-varying frequency (TVF) components where the frequency of components overlaps in time, classical linear filters fail to extract components. The *G*-filter based on time deformation has been developed to extract components of multicomponent *G*-stationary processes. In this paper, we explore the wide application of the *G*-filter for filtering different types of nonstationary processes with multiple time-frequency structure. Simulation examples illustrate that the *G*-filter can be applied to filter a broad range of multicomponent nonstationary process where TVF components may in fact overlap in time.

#### 1. Introduction and Background

The traditional linear filter is defined as , where and and where and are the input and output processes. Papoulis [1] has shown that , where and denote the power spectra of the stationary input and output processes, and , and is the frequency response function. Based on this, certain filters (e.g., the Butterworth filter [2]) have been designed to filter or pass low frequency or high frequency components of the input process to the output process. Given a cutoff frequency, a low-pass (high-pass) Butterworth filter can extract components whose frequency content is below (above) this cutoff frequency.

In general, traditional linear filters, such as the Butterworth filter, are time invariant. They are designed to extract components from stationary processes where the frequency behavior of the signal does not change with time. However, for time series data with time-varying frequency behavior (TVF), these filters can produce very poor results because the time-invariant nature of the filters does not properly account for the time-varying frequency behavior of the data. That is, these filters do not properly adjust the cutoff frequency with time according to the frequency behavior of the data. See discussion and examples in Xu et al. [3].

In order to address the filtering problem for nonstationary data with TVF, Xu et al. [3] developed the -filter utilizing the time-deformation method by deforming the index set (time axis) of the time series data. The use of time-deformation or time-warping methods to define or analyze nonstationary data has been investigated by several researchers. For example, self-similar processes are obtained by applying the Lamperti operator on the time scale of a stationary process [4]. Gray and Zhang [5] base time deformation on a log transformation of the time axis. They refer to processes that are stationary on the log scale as -stationary processes and show that the resulting spectral representation is based on the Mellin transform. Flandrin et al. [4] point out that the log transformation of Gray and Zhang [5] is a special case of the Lamperti transform. Gray et al. [6] extend the -stationary model to analyze data collected at discrete time points. Jiang et al. [7] define a more general family of nonstationary processes called -stationary processes whose frequencies monotonically change with time. See also Woodward et al. [8, Chapter 13].

In Section 2, we give a brief review of -stationary processes and the -filter along with an illustrative example. In Section 3 we extend the applicability of the -filter to filter nonstationary time varying signals with varieties of time-frequency structures. We give two examples illustrating the techniques.

#### 2. -Stationary Processes and the -Filter

In this section, we define and discuss -stationary processes, we define the -filtering procedure given by Xu et al. [3], and, provide an example.

*Definition 1 (Dual Process). *A process is called a dual process of the process if there exists a mapping of which a well-defined inverse function, , exists, so that for . The mapping is called the time-deformation function.

*Definition 2 (-stationary Process). *A process is a -stationary process if and only if a stationary dual process exists under some time-deformation function .

Some special cases of -stationary processes are given in the following:(1)when , the -stationary process is simply the traditional weakly stationary process;(2)when , is called an -stationary process ((Gray et al. (1988)), [6]);(3)when (a Box-Cox transformation), is called a -stationary process [7];(4)when with and , is called the generalized linear chirp process [9].

Jiang et al. [7] define the generalized instantaneous period (GIP) for the -stationary process to be , where is the period of the dual process, and the generalized instantaneous frequency (GIF) is . They showed that the GIF is approximately proportional to . So, for the -stationary process where , it follows that GIF , and GIP ; that is, the generalized instantaneous period is changing linearly with time. For the -stationary process where , the GIF changes like a polynomial function, that is, like . For the linear chirp stationary process where , the GIF changes linearly with time, that is, like .

Given a realization, the frequency change measured by the GIF can be visually represented using time-frequency plots. In this paper we will use the nonparametric Wigner-Ville plots for this purpose. These plots display the time-frequency behavior in the data by computing inverse Fourier transforms of windowed versions of the sample autocovariance function (see [8, 10, 11]).

Xu et al. [3] define the -filter mathematically and give a straightforward procedure for -filtering a set of -stationary data with one or more time-varying frequency components.

##### 2.1. Filtering Procedure

*Step 1. *Fit a -stationary model to the data and transform the time axis to obtain a stationary dual process.

*Step 2. *Apply traditional time-invariant filtering methods to the dual data to extract components.

*Step 3. *Convert filtered dual components back into the original scale.

See Xu et al. [3] who also point out that the -filtering procedure described previously is an application of the principle of unitary equivalence systems [12].

In Example 3 we will show how to apply the -filter to filter linear chirp data. -stationary processes have an origin, and a given realization will begin at some offsets from the process origin. For a realization of length , our procedure is to obtain a dual realization (of length ) associated with a transformation . The values, , needed for obtaining the stationary dual based on the discrete Box-Cox time transformation, , are where is called the *realization offset* and is the sampling increment. Jiang et al. [7] employ a search routine to determine the values of and offset that produce the “most stationary” dual. For each set of and offset values considered in the search, the data, , , are approximated at the ’s using interpolation. By then indexing on , the dual realization associated with the given and offset is obtained, and for each combination of and offset, the resulting dual realization is checked for stationarity. Jiang et al. [7] suggest measuring stationarity by examining characteristics (e.g., sample autocorrelations) of the first and second halves of the transformed data. They employ a -statistic for measuring the difference between the sample autocorrelations of the two halves. This measure is based on the fact that the correlation structure under stationarity stays constant across the realization. The software, GWS, written in S+, can be used to perform this search procedure, and it is available from the authors at the website http://www.texasoft.com/atsa/. This software program can be used to fit a () model to a set of data, and it provides methods for spectral analysis, forecasting, and so forth. In the examples here involving analysis of TVF data, we use the GWS software package.

*Example 3. *We consider the linear chirp process
where . This is a linear chirp stationary process consisting of two chirp components. A realization of the process is shown in Figure 1(a) where it can be seen that frequency behavior is increasing (periods are shortening) in time. It also can be seen that there appear to be more than one TVF component. The Wigner-Ville plot in Figure 1(b) clearly shows two monotonically increasing frequency components in the data, and these frequency components visually appear to be increasing linearly with time. In Figures 1(c) and 1(e), we show the components and where and . The Wigner-Ville plots associated with these two component signals are shown in Figures 1(d) and 1(f) where it can be seen that and display the lower frequency and higher frequency TVF behavior (resp.) shown in Figure 1(b).

To perform Step 1 of the -filtering procedure, we use the method of Jiang et al. [7] (using GWS software) and find that a () model with and offset 19 is chosen as the transformation which most nearly transform the data in Figure 1(a) into a stationary dual. The resulting dual data are shown in Figure 2(a). Visual inspection suggests that the dual data have stable (nonchanging) frequency behavior with more than one frequency component. The Wigner-Ville plot for the dual data is shown in Figure 2(b) where it is seen that the data contain two stable frequency components at about and . Step 2 involves filtering the dual data, which we accomplish by using low-pass and high-pass Butterworth filters with a cutoff frequency . (The -function butterworthT is available at the website listed above.)

The low and high frequency dual components are shown in Figures 2(c) and 2(d), respectively. The Wigner-Ville plot (not shown) for the data in Figure 2(c) would have the appearance of a single horizontal line at about , while the Wigner-Ville plot for Figure 2(d) would show a single horizontal line at about . Step 3 in the -filtering procedure is accomplished by transforming the dual components in Figures 2(c) and 2(d) back to the original time scale. (This is accomplished using the -function reinterpol which is available at the website.) The recovered low and high frequency TVF components of the original data are shown in Figures 3(a) and 3(b). These closely resemble the actual components in Figures 1(c) and 1(e), respectively. The Wigner-Ville plots of the two filtered components are shown in Figures 3(c) and 3(d) where it can be seen that these plots indicate linearly decreasing frequency behavior very similar to that shown in Figures 1(d) and 1(f), respectively, for the true TVF components.

It is important to note that we could not have recovered the TVF components by application of the Butterworth filter directly to the data in Figure 1(a). To see this we notice from Figure 1(b) that there is no single cutoff frequency that can be used to separate the signals across the entire time span ( to ). For example, could be used to separate the signals from to about , but beyond both TVF components are associated with frequencies below .

For more details concerning the filtering method used in Example 3, see Xu et al. [3]. In the next section we illustrate the use of the time-transformation method to filter more general types of nonstationary processes with multiple time-varying frequency structures.

#### 3. Filtering More General Nonstationary Processes Using the -Filter

In Example 3, we showed how to use the -filter to filter a nonstationary processes, called a -stationary process, by first transforming the realization from a -stationary process into a realization from a stationary process using time deformation (i.e., transforming the time scale). Specifically, in Example 3 we used the time transformation . In Figures 2(a) and 2(b), it is clear that this transformation approximately stabilized the frequency structure essentially eliminating frequency change across the realization. In many cases, however, a single stationarizing time transformation does not exist. In this section, we will illustrate the extension of the techniques used in Example 3 to such cases.

First, we discuss how the time-deformation technique will affect the spectrum of the data. Suppose that a component has a phase function ; that is, the instantaneous frequency is [11]. Let be a monotonic time-deformation function. Then in the dual space, the phase function becomes , and the instantaneous frequency is . In order to stationarize a signal with one TVF, the procedure is to find a time transformation, , such that the instantaneous frequency, , is constant.

For a process consisting of two components, one with a phase function and the other , the instantaneous frequencies of these two components are and . In the dual space, they are and . In Example 3, we obtained a single time transformation, , for which and were both approximately constant. Visually, from Figure 2(b) we see that and are represented by the upper and lower horizontal lines, respectively, in the Wigner-Ville plot. This allowed us to select a cutoff frequency between 0.02 and 0.09 for purposes of separating the two components in the dual process.

More generally, it is important to note that regardless of whether and are constants, it follows that if and for all , then In other words, the component with a higher instantaneous frequency in the original time space still has a higher instantaneous frequency in the dual space. Suppose there exists a constant frequency, , such that or for all (assuming ). Note that (4) says that stays between the upper and lower curves shown in the Wigner-Ville plot. From (3) we see that and are depicted by the upper and lower Wigner-Ville curves in the dual space, and in this case a Butterworth filter with cutoff frequency should be able to separate these two components in the dual space. The important point is that this separation is possible even if the time transformation did not produce a stationary dual. Thus, if we can find a well-defined time-deformation function and a constant , such that , for , then these two components can be extracted from the original data. Example 4 illustrates these comments.

*Example 4. *Here, we consider a simulated data set which is the signal consisting of the sum of two quadratic chirps,
for . The data and the Wigner-Ville time-frequency distribution are displayed in Figures 4(a) and 4(b). The lower frequency TVF component, , is shown in Figure 4(c), while the higher frequency component, , is shown in Figure 4(d). By examining the Wigner-Ville plot in Figure 4(b) and the component plots in Figures 4(c) and 4(d), it can be seen that the frequency behavior in for is similar to the frequency behavior in for . Consequently, a simple Butterworth filter will not be able to separate the two components.

We apply the -filtering procedure in an attempt to recover the two TVF components. Step 1 involves using the procedure of Jiang et al. [7] in an attempt to stationarize the signal in Figure 4(a). However, the nature of the two quadratic chirp components is such that there is not a single time transformation, , that will stationarize the signal . Using the GWS software, a time transformation with and offset = 3 is selected as the transformation that most nearly stationarizes the data. Figures 5(a) and 5(b) show the “dual” data and the associated Wigner-Ville plot where it can be seen that the dual does not appear to be stationary. The key point, however, is that even though the transformation has not stationarized the data (the frequency lines in the Wigner-Ville plot are not straight horizontal lines), the upper and lower frequency components can be separated with a simple Butterworth filter at, for example, a cutoff frequency . Step 2 in the -filtering procedure can thus be performed, and the resulting low and high frequency components after filtering the dual data with a Butterworth filter (e.g., using -function butterworthT) are shown in Figures 5(c) and 5(d), respectively. We complete the example by applying Step 3 to plot the filtered dual components back on the original time scale. In Figures 6(a) and 6(b), we show the filtered data plotted on the original time scale (using -function reinterpol), and in Figures 6(c) and 6(d) we show the Wigner-Ville plots for these components. The Wigner-Ville plots in Figures 6(c) and 6(d) are good approximations to the lower and upper curves, respectively, in Figure 4(b) associated with the lower and upper components of the original signal. We comment that the poor behavior at the endpoints is a well-known property of filters such as the Butterworth, and no endpoint adjustments have been made here. Also, the loss of amplitude in the recovered higher frequency regions (see, e.g., Figure 6(d) compared to Figure 4(b)) is at least partially due to the linear interpolation that is used. Haney [13] has shown that the use of Fourier-based interpolation may be able to adjust for some of the signal loss.

Note that Example 4 illustrates the fact that if we can find a time transformation, , and frequency, , such that (4) holds, then the two time-varying frequencies can be separated. This is illustrated in Figure 7(a) where the dotted curve represents for and where is the transformation used in Example 4 with and offset 3; that is, The solid curves are the upper and lower Wigner-Ville curves shown in Figure 4(b) which represent and . It can be seen that for the time transformation in (6), the curve stays between the upper and lower Wigner-Ville curves. Figure 7(b) shows upper and lower Wigner-Ville curves in the associated dual space; that is, from (3) these are and , respectively. The straight dashed line corresponds to .

In Figures 7(c) and 7(d) we show the corresponding plots associated with making a time transformation It can be seen that either of the time transformation in (6) or (7) provides separation in the dual space sufficient for application of a filter such as the Butterworth.

#### 4. Extracting TVF Components Entering the Signal at Different Times

In this section we discuss methods for extracting components in more complicated TVF settings in which components may enter and/or exit the overall signal at different times. The procedure will be illustrated in Example 5. Suppose the high (low) frequency component in a TVF signal is -stationary over a certain time span, and we apply the appropriate time transformation to stationarize this -stationary component over that time span. Then by (2), on the transformed version of the time span being considered, this component will become stationary with a constant frequency that is still higher (lower) than the instantaneous frequency of the other component in the dual space. The key idea addressed in Example 5 is that in order to extract signals, it may be useful to segment the overall realization length into sections over which straightforward signal separations can be obtained using the techniques illustrated in Examples 3 and 4. We use these ideas in the following example.

*Example 5. *In this example, we simulate a situation similar to that studied by Papandreou-Suppappola and Suppappola [14]. The synthetic signal is shown in Figure 8(a), and the Wigner-Ville plot of this realization is shown in Figure 8(b). From Figure 8(b) we observe the following.(1)The overall signal consists of three TVF component signals. We will refer to these as Components 1, 2, and 3. Based on the Wigner-Ville plot we conclude that(a)component 1 exists for approximately ;(b)component 2 exists for approximately ;(c)component 3 exists for approximately .(2)Components 1 and 3 have Wigner-Ville curves that appear to be straight lines, suggesting that each of these is a linear chirp signal.(3)Component 1 exists in the signal by itself for about .(4)The frequencies associated with the components intersect at two places:(a)frequencies associated with components 1 and 2 intersect at about ;(b)frequencies associated with components 2 and 3 intersect at about .

Figures 8(b)-8(c) show the component signals actually used to create the actual signal, , which is the sum of the components; that is, . We see that our observations in conclusion 1 mentioned above are correct. Also, it should be noted that components 1 and 3 were generated as linear chirps (as was speculated in conclusion 2 mentioned above) and component 2 is -stationary, but this information will not be assumed to be known in the following analysis.

*Step 1 (Extract Component 1). *Since component 1 exists by itself for , we use the method of Jiang et al. [7] (using the GWS software) to find a model for stationarizing the signal composed of the first 100 points. A model is selected with offset 0. We apply this time transformation to the first 300 data values, and the resulting dual process is shown in Figure 9(a). The associated Wigner-Ville plot is shown in Figure 9(b) where it is seen that component 1 has now been stationarized (the associated frequency behavior has been stabilized, and it appears as a horizontal line at about in the Wigner-Ville plot). Note that although the segment of the original data to which we applied the time deformation () and the resulting dual data set both have 300 time values, the original time axis has been deformed. Consequently, the first 100 data values (involving only component 1) in the original time scale correspond to about the first 50 values in transformed time. By examining Figure 9(b) we see that the dual data are now split into three segments in the deformed time axis as follows: (i)the first segment (1–50) which contains only the transformed component 1 signal; (ii)the middle segment (51–100) where the frequencies intersect;(iii)the last segment (101–300) where the frequencies are well separated. In this segment the transformed component 1 signal is consistently higher frequency than the transformed component 2 and component 3 signals. In this range, a constant cutoff frequency (e.g., 0.12) can be used to separate component 1 from the other two components.

In order to extract the stationarized dual signal associated with component 1 we note that no extraction is necessary for the first 50 dual time values (the only signal present is the one associated with the first linear chirp component). As suggested in note (iii) above, extraction of component 1 in the range 101–300 can be accomplished using a high-pass Butterworth filter.

Extraction of the first linear chirp component in the range 51–100 is more difficult. The procedure we have selected is to fit models to the two extracted signals (i.e., 1–50 and 101–300) and then to use these models to forecast and backcast on the middle segment (51–100). The resulting extracted dual signal associated with component 1 is shown in Figure 9(c). Transforming this signal back to the original time scale gives Figure 9(d) which is similar to the true component shown in Figure 8(c). Note again that although the first dual component terminates at about a transformed time value of 210, the first component in the original time scale terminates at about as expected.

*Step 2 (Extract Component 2). *In the second step, we will extract component 2. We first subtract the recovered component 1 (shown in Figure 9(d)) from the original signal leaving only components 2 and 3. The resulting data () and corresponding Wigner-Ville plot are shown in Figures 10(a) and 10(b). We note the following:(i)only TVF components 2 and 3 remain;(ii)component 2 exists by itself for .

Based on note (ii) above, we apply the method of Jiang et al. [7] to find a stationarizing transformation for the data in time segment 101–200 in Figure 10(a). Using this method, a model with is selected with offset 69. Applying this time transformation to the data for , we obtain the dual data shown in Figure 10(c) and associated Wigner-Ville plot in Figure 10(d). From the Wigner-Ville plot in Figure 10(d) we note the following:(i)the second component seems to have been stationarized using this procedure (frequency behavior is represented by a horizontal line at about in Wigner-Ville plot);(ii)the first segment (101–250) contains only the stationarized component 2;(iii)transformed component 2 is consistently lower frequency than transformed component 3 for the range 251–280 and the signals could be separated by a constant cutoff frequency of about 0.075;(iv)the frequencies intersect in the range 281–340;(v)the transformed component 2 is consistently lower frequency than transformed component 3 for the range 341–400. A constant cutoff frequency of about 0.175 could be used to separate the signals in this range.

Based on these comments we apply a high-pass Butterworth filter to extract component 2 in the range 101–280 and a low-pass Butterworth filter to extract the dual for component 2 for the range 341–400. Forecasting and backcasting are used as in Step 1 for recovering the component 2 dual for the range 281–340. The resulting extracted signal associated with the dual for component 2 is shown in Figure 10(e), and Figure 10(f) shows this component transformed back to the original time scale.

*Step 3 (Extract Component 3). *The third component is obtained by subtracting the two recovered components from the original data; that is, , for where represents the recovered th component.

Figure 11 summarizes the results showing the three recovered components and associated Wigner-Ville plots which correspond well with those shown Figure 8.

#### 5. Conclusion

For nonstationary time series with multiple time-varying frequency structure, especially where the frequencies of components overlap over time, traditional linear filters are not able to successfully extract components. In order to address this problem, the -filter, a time-variant filter based on the time-deformation method, was studied by Xu et al. [3]. In that paper the -filter was applied to filter components from certain types of -stationary processes by transforming them to stationarity using an appropriate time transformation of the time scale.

In practical situations, most nonstationary data with multiple time-varying frequency structure cannot be transformed to stationary by applying any time deformation. In this paper, we showed that the -filtering procedure can be extended to cases in which the individual component signals may or may not be able to be fully stationarized using a time transformation. We first discussed the case (see Example 4) in which the frequency behaviors of the components of a process can be separated by a time transformation even though a stationarizing transformation is not available. This significantly extends the class of TVF processes that can be filtered using the -filter. Secondly, we show in Example 5 that the -filter can be used to filter signals with multiple -stationary components having varying points of entry into and exit from the signal.

#### References

- A. Papoulis,
*Probability, Random Variables, and Stochastic Processes*, McGraw-Hill, New York, NY, USA, 1984. View at MathSciNet - S. Butterworth, “On the theory of tilter amplifiers,”
*Experimental Wireless and the Wireless Engineer*, vol. 7, pp. 536–541, 1930. View at Google Scholar - M. Xu, K. B. Cohlmia, W. A. Woodward, and H. L. Gray, “G-filtering nonstationary time series,”
*Journal of Probability and Statistics*, vol. 2012, Article ID 738636, 15 pages, 2012. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - P. Flandrin, P. Borgnat, and P. O. Amblard, “From stationarity to self-similarity, and back: variations on the Lamperti transformation,” in
*Processes With Long-Range Correlations*, G. Raganjaran and M. Ding, Eds., pp. 88–117, Springer, New York, NY, USA, 2003. View at Google Scholar - H. L. Gray and N. F. Zhang, “On a class of nonstationary processes,”
*Journal of Time Series Analysis*, vol. 9, no. 2, pp. 133–154, 1988. View at Publisher · View at Google Scholar - H. L. Gray, Chu-P.C. Vijverberg, and W. A. Woodward, “Nonstationary data analysis by time deformation,”
*Communications in Statistics*, vol. 34, no. 1, pp. 163–192, 2005. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - H. Jiang, H. L. Gray, and W. A. Woodward, “Time-frequency analysis—G(
*λ*)-stationary processes,”*Computational Statistics & Data Analysis*, vol. 51, no. 3, pp. 1997–2028, 2006. View at Publisher · View at Google Scholar · View at MathSciNet - W. A. Woodward, H. L. Gray, and A. C. Elliott,
*Applied Time Series Analysis*, Chapman and Hall/CRC, Boca Raton, Fla, USA, 2012. View at MathSciNet - S. D. Robertson, H. L. Gray, and W. A. Woodward, “The generalized linear chirp process,”
*Journal of Statistical Planning and Inference*, vol. 140, no. 12, pp. 3676–3687, 2010. View at Publisher · View at Google Scholar · View at MathSciNet - W. Martin and P. Flandrin, “Wigner-Ville spectral analysis of nonstationary processes,”
*IEEE Transactions on Acoustics, Speech, and Signal Processing*, vol. 33, no. 6, pp. 1461–1470, 1985. View at Google Scholar · View at Scopus - B. Boashash,
*Time Frequency Analysis*, Elsevier, Oxford, UK, 2003. - R. G. Baraniuk and D. L. Jones, “Unitary equivalence: a new twist on signal processing,”
*IEEE Transactions on Signal Processing*, vol. 43, no. 10, pp. 2269–2282, 1995. View at Publisher · View at Google Scholar · View at Scopus - J. R. Haney,
*Analyzing time series with time-varying frequency behavior and conditional heteroskedasticity [Ph.D. thesis]*, Southern Methodist University, Department of Statistical Science, Dallas, Texas, USA, 2010. - A. Papandreou-Suppappola and S. B. Suppappola, “Analysis and classification of time-varying signals with multiple time-frequency structures,”
*IEEE Signal Processing Letters*, vol. 9, no. 3, pp. 92–95, 2002. View at Publisher · View at Google Scholar · View at Scopus