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

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 𝐴(𝑓)=∫^∞_−∞𝑎(𝑡)𝑒^{−2𝜋𝑖𝑓𝑡}𝑑𝑡 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 𝑔(𝑡)=ln(𝑡),𝑡∈(0,∞), 𝑋(𝑡) is called an 𝑀-stationary process ((Gray et al. (1988)), [6]);(3)when 𝑔(𝑡)=(𝑡^𝜆−1)/𝜆,𝑡∈(0,∞) (a Box-Cox transformation), 𝑋(𝑡) is called a 𝐺(𝜆)-stationary process [7];(4)when 𝑔(𝑡)=𝑎𝑡²+𝑏𝑡,𝑡∈(0,∞) with 𝑎>0 and 𝑏≥0, 𝑋(𝑡) 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 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 𝑔(𝑡)=𝑎𝑡²+𝑏𝑡, 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 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, 𝑘=(𝑡^𝜆_𝑘−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 𝑌(𝑡)=cos0.0002𝜋𝑡²2+cos0.0008𝜋𝑡²2, (1) 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 𝑌₁(𝑡) and 𝑌₂(𝑡) where 𝑌₁(𝑡)=cos(0.0002𝜋𝑡²/2) and 𝑌₂(𝑡)=cos(0.0008𝜋𝑡²/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 𝑌₁(𝑡) and 𝑌₂(𝑡) display the lower frequency and higher frequency TVF behavior (resp.) shown in Figure 1(b).

(a)  TVF signal

(b) Wigner-Ville plot for (a)

(c) Low frequency component

(d) Wigner-Ville plot for (c)

(e) High frequency component

(f) Wigner-Ville plot for (e)

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).

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.)

(a)  Data after time transformation

(b) Wigner-Ville plot for (a)

(c) Low-pass filtered data

(d) Recovered high frequency TVF signal

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).

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 𝑓=0.025, while the Wigner-Ville plot for Figure 2(d) would show a single horizontal line at about 𝑓=0.075. 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.

(a)   Recovered low frequency TVF signal

(b) Recovered high frequency TVF signal

(c) Wigner-Ville plot for (a)

(d) Wigner-Ville plot for (b)

Figure 3 

(a) Recovered low frequency TVF component; (b) recovered high frequency TVF component; (c) Wigner-Ville plot for (a); (d) Wigner-Ville plot for (b).

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 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 𝑔(𝑡)=(𝑡^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 time-deformation function. Then in the dual space, the phase function becomes 2𝜋𝜑(𝑔⁻¹(𝑢)), and the instantaneous frequency is 𝑑𝜑(𝑔⁻¹(𝑢))/𝑑𝑢=𝜑^′(𝑔⁻¹(𝑢))(𝑑𝑔⁻¹(𝑢)/𝑑𝑢)=𝜑^′(𝑔⁻¹(𝑢))(1/𝑔^′(𝑡))=𝜑^′(𝑡)/𝑔^′(𝑡)=𝜑^′(𝑔⁻¹(𝑢))/𝑔^′(𝑔⁻¹(𝑢)). 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𝜋𝜑₁(𝑡) and the other 2𝜋𝜑₂(𝑡), 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, 𝑔(𝑡)=(𝑡^2.1−1)/2.1=𝑢, for which 𝜑^′₁(𝑡)/𝑔^′(𝑡) and 𝜑^′₂(𝑡)/𝑔^′(𝑡) were both approximately constant. Visually, from Figure 2(b) we see that 𝜑^′₁(𝑡)/𝑔^′(𝑡)≈0.09 and 𝜑^′₂(𝑡)/𝑔^′(𝑡)≈0.02 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 𝑔^′(𝑡)>0 for all 𝑡∈(𝑎,𝑏), then 𝜑^′₁𝑔⁻¹(𝑢)𝑔^′𝑔⁻¹(𝑢≥𝜑)^′2𝑔⁻¹(𝑢)𝑔^′𝑔⁻¹(𝑢.) (2) 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 𝜑^′₁(𝑡)𝑔^′(𝑡)≤𝑓_𝑐≤𝜑^′₂(𝑡)𝑔^′(𝑡) (3) or 𝜑^′₁(𝑡)≤𝑓_𝑐𝑔^′(𝑡)≤𝜑^′₂(𝑡) (4) 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 𝜑^′₁(𝑡)/𝑔^′(𝑡) 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, 𝑋(𝑡)=cos2𝜋−𝑡+50100−3³+0.15𝑡+cos2𝜋𝑡+50100−1³+0.14𝑡=𝑋_𝐿(𝑡)+𝑋_𝐻(𝑡) (5) for 0<𝑡<250. 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 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.

(a) TVF signal

(b) Wigner-Ville plot for (a)

(c) Low frequency TVF component

(d) High frequency TVF component

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.

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) 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.

(a) Data after time transformation

(b) Wigner-Ville plot for (a)

(c) Low-pass filtered data from (a)

(d) High pass filtered data from (a)

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.

(a)   Recovered low frequency TVF signal

(b) Recovered high frequency TVF signal

(c) Wigner-Ville plot for (a)

(d) Wigner-Ville plot for (b)

Figure 6 

(a) Figure 5(c) plotted on original time scale; (b) Figure 5(d) plotted on original time scale; (c) Wigner-Ville plot for data in (a); (d) Wigner-Ville plot for data in (b).

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, 𝑔(𝑡)=(𝑡+3)^1.7−1.1.7 (6) 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 𝑓_𝑐=0.13.

(a)  𝜆=1.7

(b) Dual plot for (a)

(c) 𝜆=2

(d) Dual plot for (c)

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.

In Figures 7(c) and 7(d) we show the corresponding plots associated with making a time transformation 𝑔(𝑡)=(𝑡+68)²−12. (7) 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 1≤𝑡≤250;(b)component 2 exists for approximately 𝑡≥100;(c)component 3 exists for approximately 𝑡≥200.(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:(a)frequencies associated with components 1 and 2 intersect at about 𝑡=150;(b)frequencies associated with components 2 and 3 intersect at about 𝑡=275.
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, 𝑋_𝑡=𝑋_𝑡⁽¹⁾+𝑋_𝑡⁽²⁾+𝑋_𝑡⁽³⁾,𝑡=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)  Data for Example 5

(b) Wigner-Ville plot for (a)

(c) Component 1

(d) Component 2

(e) Component 3

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.

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.

(a)  Dual data for Example 5, Step 1

(b) Wigner-Ville plot for (a)

(c) Extracted dual component 1

(d) Recovered chirp component 1

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.

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) and 10(b). We note the following:(i)only TVF components 2 and 3 remain;(ii)component 2 exists by itself for 101≤𝑡≤200.
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 𝜆=0 is selected with offset 69. Applying this time transformation to the data for 101≤𝑡≤400, 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 𝑓=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 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.

(a)  Data (101–400) with First comp. subtracted

(b) Wigner-Ville plot for (a)

(c) Dual data (101–400)

(d) Wigner-Ville plot for (c)

(e) Extracted dual component 2

(f) Recovered component 2

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.

Step 3 (Extract Component 3). The third component is obtained by subtracting the two recovered components from the original data; that is, 𝑋_𝑡⁽³⁾=𝑋_𝑡−𝑋_𝑡⁽¹⁾−𝑋_𝑡⁽²⁾, for 𝑡=1,…,400 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.

(a)  Recovered first component (linear chirp)

(b) Wigner-Ville plot for (a)

(c) Recovered second component (𝑀-stationary)

(d) Wigner-Ville plot for (c)

(e) Recovered third component (linear chirp)

(f) Wigner-Ville plot for (e)

Figure 11 

(a), (c), and (e): recovered components for Example 5 and (b), (d), and (e): the corresponding Wigner-Ville plots.

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.