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

For nonstationary time series consisting ofmultiple time-varying frequency (TVF) componentswhere 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.

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 timefrequency structures.We give two examples illustrating the techniques.

𝐺-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.
Jiang et al. [7] define the generalized instantaneous period (GIP) for the -stationary process to be (; , ) =  −1 (()+)−, where  is the period of the dual process, and the generalized instantaneous frequency (GIF) is 1/(; , ).They showed that the GIF is approximately proportional to   ().So, for the -stationary process where () = ln(), it follows that GIF ∝ 1/, and GIP ∝ ; that is, the generalized instantaneous period is changing linearly with time.For the ()-stationary process where () = (  − 1)/, the GIF changes like a polynomial function, that is, like  −1 .For the linear chirp stationary process where () =  2 + , the GIF changes linearly with time, that is, like 2 + .
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 timefrequency 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.

𝐺−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,  = (   − 1)/Δ   − , are   = (( + )Δ   + 1) 1/ where (Δ   + 1) 1/ 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, (  ),  = 1, 2, . .., 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  ∈ [0, 400].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  1 () and  2 () where  1 () = cos(0.0002 2 /2) and  2 () = cos(0.0008 2 /2).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  1 () and  2 () 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  = 2.1 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  = 0.025 and  = 0.075.Step 2 involves filtering the dual data, which we accomplish by using low-pass and high-pass Butterworth filters with a cutoff frequency   = 0.05.(The -function butterworthT is available at the website listed above.) The low and high frequency dual components are shown in Figures 2(c 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 ( = 0 to  = 400).For example,   = 0.1 could be used to separate the signals from  = 0 to about  = 150, but beyond  = 300 both TVF components are associated with frequencies below   = 0.1.
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 timevarying frequency structures.

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 () = ( 2.1 − 1)/2.1 = .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 2(); that is, the instantaneous frequency is   () [11].Let  = () be a monotonic timedeformation function.Then in the dual space, the phase function becomes 2( −1 ()), 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 2 1 () and the other 2 2 (), the instantaneous frequencies of these two components are   1 () and   2 ().In the dual space, they are   1 ()/  () =   1 ( −1 ())/   ( −1 ()) and   2 ()/  () =   2 ( −1 ())/  ( −1 ()).In Example 3, we obtained a single time transformation, () = ( 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   () > 0).Note that (4) says that     () stays between the upper and lower curves shown in the Wigner-Ville plot.From (3) we see that   1 ()/  () and   2 ()/  () 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 welldefined time-deformation function  = () and a constant   , such that   1 () ≤     () ≤   2 (), for  ∈ (, ), then these two components can be extracted from the original data.Example 4 illustrates these comments.200 <  < 250 is similar to the frequency behavior in   () for 0 <  < 50.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  = 1.7 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  = 0.13.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 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   = 0.13 and where () is the () transformation used in Example 4 with  = 1.7 and offset 3; that is, The solid curves are the upper and lower Wigner-Ville curves shown in Figure 4(b) which represent   1 () and   2 ().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   1 ()/  () and   2 ()/  (), respectively.The straight dashed line corresponds to   = 0.13.
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.

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   (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  ≤ 100.
(4) The frequencies associated with the components intersect at two places: to create the actual signal,   , which is the sum of the components; that is,   =  (1)   +  (2)   +  (3)   ,  = 1, . . ., 400.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  ≤ 100, 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 (2) 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  = 0.15 in the Wigner-Ville plot).Note that although the segment of the original data to which we applied the time deformation (1 ≤  ≤ 300) 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 AR() 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  = 250 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 (101 ≤  ≤ 400) and corresponding Wigner-Ville plot are shown in Figures 10(a 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  Step 3 (Extract Component 3).The third component is obtained by subtracting the two recovered components from the original data; that is, X(3)  =   − X(1)  − X(2)  , for  = 1, . . ., 400 where X()  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.

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.

Figure 1 :
Figure 1: (a) A realization from the linear chirp model in Example 3; (b) the Wigner-Ville time-frequency plot for the data in (a); (c) lower frequency TVF component,  1 (); (d) Wigner-Ville plot for the signal in (c); (e) higher frequency TVF component,  2 (); (f) Wigner-Ville plot for the signal in (e).

Figure 2 :
Figure 2: (a) The dual data obtained from Figure 1(a) using GWS software; (b) Wigner-Ville time-frequency plot for (a); (c) recovered low frequency component from (a); (d) recovered high frequency component from (a).

Figure 4 :
Figure 4: (a) A realization of the process in Example 4; (b) the Wigner-Ville time-frequency distribution of the data; (c), (d) the two components.

Figure 5 :
Figure 5: (a) Data on transformed time axis; (b) the Wigner-Ville time-frequency distribution of (a); (c) data in (a) after low-pass filter; (d) data in (a) after high-pass filter.

Figure 7 :
Figure 7: (a) Sketch of upper and lower Wigner-Ville curves for Example 4 showing   () for   = 0.13 and () given in (6); (b) sketch of Wigner-Ville curves in dual space; (c) sketch of upper and lower Wigner-Ville curves for Example 4 showing   () for   = 0.13 and () given in (7); (d) sketch of Wigner-Ville curves in dual space.

3 Figure 8 :
Figure 8: (a) A realization of the process in Example 5; (b) the Wigner-Ville time-frequency distribution of the data; (c), (d), and (e) the three components.

chirp component 1 Figure 9 :
Figure 9: (a) The dual data (1 : 300) after the first time deformation; (b) the Wigner-Ville time-frequency distribution of the dual; (c) the extracted stationary dual component; (d) the recovered first linear chirp component.

( 1 )
(a), and the Wigner-Ville plot of this realization is shown in Figure 8(b).From Figure 8(b) we observe the following.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 1 ≤  ≤ 250; (b) component 2 exists for approximately  ≥ 100; (c) component 3 exists for approximately  ≥ 200.

2 Figure 10 :
Figure 10: (a) The data (101 : 400) with the first filtered component subtracted; (b) the Wigner-Ville time-frequency distribution of the data; (c) the dual data; (d) the Wigner-Ville time-frequency distribution of the dual; (e) the extracted stationary dual component; (f) the recovered -stationary component.
Figures 8(b)-8(c) show the component signals actually usedto create the actual signal,   , which is the sum of the components; that is,   = (1)   + (2)   + (3)   ,  = 1, . . ., 400.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.
(a).Using this method, a () model with  = 0 is selected with offset 69.Applying this time transformation to the data for 101 ≤  ≤ 400, we obtain the dual data shown in Figure10(c) and associated Wigner-Ville plot in Figure10(d).From the Wigner-Ville plot in Figure10(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  = 0.1 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 lowpass Butterworth filter to extract the dual for component 2 for the range 341-400.Forecasting and backcasting are used as inStep 1 for recovering the component 2 dual for the range 281-340.The resulting extracted signal associated with the dual Figure 11: (a), (c), and (e): recovered components for Example 5 and (b), (d), and (e): the corresponding Wigner-Ville plots.