<figure><img src="image_1.png" /></figure>

HAL
open science

# An Improved Stationarity Test Based on Surrogates

Douglas Baptista de Souza, Jocelyn Chanussot, Anne-Catherine Favre, Pierre Borgnat

► To cite this version:

Douglas Baptista de Souza, Jocelyn Chanussot, Anne-Catherine Favre, Pierre Borgnat. An Improved Stationarity Test Based on Surrogates. IEEE Signal Processing Letters, 2019, 26 (10), pp.1431-1435. 10.1109/LSP.2019.2931150. [hal-02307460](https://hal.science/hal-02307460v1)

HAL Id: hal-02307460

https://hal.science/hal-02307460v1

Submitted on 4 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# An improved stationarity test based on surrogates

Douglas Baptista de Souza, Jocelyn Chanussot, Fellow, IEEE,
Anne-Catherine Favre, and Pierre Borgnat, Member, IEEE

**Abstract**—Over the last years, several stationarity tests have been proposed. One of these methods uses time-frequency representations and stationarized replicas of the signal (known as surrogates) for testing wide-sense stationarity. In this paper, we propose a procedure to improve the original surrogate test. The proposed methodology can be seen as a guideline on how the surrogate test can be improved. We show mathematically that the modified test should exhibit improved classification performance. Numerical simulations on synthetic and real-world signals are carried out to evaluate the modified test against competing ones.

**Index Terms**—Stationarity test, nonstationary signals, time-frequency analysis, surrogates.

## I. INTRODUCTION

TESTING stationarity is important in signal processing and other areas, such as economic and environmental sciences [1], [2]. Stationarity is usually defined as the temporal invariance of the statistical properties of the signal. In practice, researchers are mostly interested in testing for stationarity up to the second order [3]. A signal is called wide-sense stationary (WSS) if its first and second-order statistics are invariant in time [4]. The spectral content of a WSS signal is described by a time-invariant power spectrum density (PSD). Nonstationary signals, on the other hand, have time-dependent spectra that can be characterized by time-frequency (TF) representations [4], [5]. The advantages and disadvantages of different TF representations are well documented in the literature [4]–[6]. The various TF methods can be classified as parametric [7], [8] or nonparametric [9]–[13]. Since parametric methods rely on the adherence to the chosen models, which can be hard to assess in real-world datasets, nonparametric methods are often preferable when testing real-world signals [13]–[16].

Among the nonparametric TF methods, approaches based on the Wigner-Ville spectrum (WVS) and its modifications have been commonly used for analyzing nonstationary signals [17]–[20]. Several years ago, a test for wide-sense stationarity using the WVS was developed [21]. The idea of the test was to quantify the difference between local and global spectral

Manuscript received January xx, 2019; accepted xxxxxxxxx xx, 2019. Date of publication xxxxxxxxx xx, 2019; date of current version xxxxxxxxxxx xx 2019. The associate editor coordinating the review of this manuscript was xxxxxxxxxxxxxxxxxxxxxxxxx. (Corresponding author: Douglas Baptista de Souza.)

D. Baptista de Souza is with the Data Analytics Team at GE Renewable Energy, São Paulo, SP, Brazil (e-mail: douglas.souza@ge.com).

J. Chanussot is with the CNRS, GIPSA-LAB, Université Grenoble Alpes, 38000 Grenoble, France (e-mail: jocelyn.chanussot@gipsalab.grenoble-inp.fr).

A.-C. Favre is with the Institut des Géosciences de l'Environnement, Université Grenoble Alpes, 38000 Grenoble, France (e-mail: Anne-Catherine.FavrePugin@ense3.grenoble-inp.fr).

P. Borgnat is with the Laboratoire de Physique, ENS de Lyon, 69364 Lyon, France (e-mail: Pierre.Borgnat@ens-lyon.fr).

This letter has supplementary downloadable material available at http://ieeexplore.ieee.org, provided by the author.

features in the estimated WVS of the signal. Because of statistical fluctuations, these local and global features are not identical, even in case of stationarity. Hence, a hypothesis test was designed to test if the observed fluctuations were due to the supposed nonstationarity. To characterize the null hypothesis of stationarity, the authors proposed to use stationarized replicas of the signal known as surrogates [21]–[23].

The surrogate test works well for second-order nonstationary processes and signals modulated in amplitude or frequency [21]. On the other hand, the performances of the test decrease when testing signals with first-order or slowly-varying nonstationarities. To alleviate this problem, the test has been changed in [24] to improve the detection of first-order nonstationarities. However, the obtained performance gain was small and the applicability of the modified test was only verified empirically.

In this letter, we propose an improved surrogate test that is not dependent on the type of nonstationary signal being probed. This is an important point when testing real-world signals, whose underlying nonstationary characteristics are often unknown. The proposed technique can be seen as a generalization of [21] and [24]. More specifically, we first establish how the parameters of an improved version of the test should differ from the original ones. Then, we present a method to transform the original test parameters into the modified (desired) ones, provided that a time-varying feature with some special properties can be extracted from the TF representation of the signal. These properties are determined and we show mathematically and through simulations that the test can indeed be improved by the proposed modification.

This paper is organized as follows. In Section II, the original surrogate test is reviewed. Section III is devoted to the development of the modified framework. The experimental study and conclusions are shown in Sections IV and V, respectively.

## II. BACKGROUND ELEMENTS

The WVS of a given discrete-time signal $x(n)$ can be estimated by means of the multitaper spectrogram

$$S(n, f) = \frac{1}{K} \sum_{k=1}^{K} S_{h_k}(n, f) \quad (1)$$

obtained by averaging $k = 1, \dots, K$ spectrograms [5]

$$S_{h_k}(n, f) = \left| \sum_m x(n+m)h_k(m)e^{-j2\pi fm} \right|^2 \quad (2)$$

computed with K Hermite window functions given by

$$h_k(n) = e^{-n^2/2} H_k(n) / \sqrt{\pi^{1/2} 2^k k!} \quad (3)$$

where $H_k(n)$ are Hermite polynomials following the recursion $H_k(n) = 2nH_{k-1}(n) - 2(k-2)H_{k-2}(n)$ for $k \ge 2$ and the
initialization $H_0(n) = 1$ and $H_1(n) = 2n$ [21]. The Hermite windows are orthonormal and maximally concentrated in TF domain [25]. In [21], a free parameter $n_{h0}$ controls the length of $h_k(n)$ and the $n = 1, ..., M$ time points to compute $S(n, f)$.

For the stationarity test of [21], one evaluates if the local spectra $S(n, f)$ computed at these $1, ..., M$ time points are similar to the global average spectrum defined as

$$ \bar{S}(f) = \frac{1}{M} \sum_{n=1}^{M} S(n, f). \qquad (4) $$

For a given $n$, the similarity between (1) and (4) at frequency $f$ can be evaluated by a measure of local spectral mismatch

$$ d(n, f) = D[S(n, f), \bar{S}(f)]. \qquad (5) $$

To account for the combined effect in frequency of all individual values of $d(n, f)$, we integrate (5) over frequency

$$ c(n) = \sum_{f} d(n, f). \qquad (6) $$

Different distances taking form of (6) can be chosen, such as Kullback-Leibler, log-spectral and Itakura-Saito<fnref n="1" /> [21], [26]. If $x(n)$ is stationary, $S(n, f)$ reduces to the time-invariant PSD and the following approximation holds [4]:

$$ S(n, f) \approx \bar{S}(f), \forall n. \qquad (7) $$

Thus, $c(n) \approx 0$ if $x(n)$ is stationary. In practice, we always observe fluctuations in $S(n, f)$, which implies fluctuations in (6). The idea is to evaluate the significance of these fluctuations through stationarized versions of $x(n)$ known as surrogates. The surrogates are obtained by shuffling the phase of $x(n)$, which controls potential nonstationarities in the signal [27]. The phase of a surrogate is a resample of the original phase of $x(n)$. The surrogates are stationary and they have the same global average spectrum $\bar{S}(f)$ of $x(n)$ [21]. The test evaluates if the fluctuations observed in (6) are similar to those in the surrogate data. The test statistics $\Theta_1$ is defined as

$$ \Theta_1 = \mathrm{Var}[c(n)] \qquad (8) $$

which measures the fluctuations of $c(n)$. In (8), $\mathrm{Var}[\cdot]$ is the variance. To perform the hypothesis test, a sample with values of $\Theta_1$ under the null hypothesis of stationarity should be obtained. To do so, we first compute a collection of $J$ surrogates $s_1(n), ..., s_J(n)$ and their corresponding multitaper spectrograms $S_1(n, f), ..., S_J(n, f)$. Then, for each surrogate $s_j(n)$, the distances between local spectra and their global average are calculated as in (5), i.e., $d_j(n, f) = D[S_j(n, f), \bar{S}_j(f)]$, and then integrated over frequency to obtain

$$ c_j(n) = \sum_{f} d_j(n, f). \qquad (9) $$

This gives a collection of $J$ values of the test statistics under the null hypothesis of stationarity

$$ \Theta_0 = \{\mathrm{Var}[c_j(n)]\}, j = 1, ..., J. \qquad (10) $$

In practice, the variances in (8) and (10) are estimated by using the sample variance estimator [21]. It has been shown that the distribution of (10) can be modelled by a gamma probability density function (PDF) [28]. By fitting a gamma model to (10), a threshold $\gamma$ for stationarity can be derived given a prescribed false alarm rate. Taking into account these points, a one-sided hypothesis test<fnref n="2" /> for $x(n)$ can be built for testing the null hypothesis of stationarity ($H_0$) against the alternative one of nonstationarity ($H_1$)

$$ H(x) = \begin{cases} H_1 & \text{if } \Theta_1 > \gamma, \text{ "nonstationary"}, \\ H_0 & \text{if } \Theta_1 \le \gamma, \text{ "stationary"}. \end{cases} \qquad (11) $$

Finally, note that the test described in this section contains many steps (e.g., distance functions, surrogates, gamma modeling), leaving room for improvements in some parts. In the following section, we discuss how to improve this test by modifying some of its parameters.

### III. MODIFYING THE STATIONARITY TEST

#### A. How can the original test be improved?

Let $H'(x)$ represent a modified version of $H(x)$ in (11), which is obtained by changing parameters $\Theta_1$ and $\Theta_0$ as given in (8) and (10), respectively. Let $\Theta'_1$ and $\Theta'_0$ be the modified versions of these parameters. Here, we aim to find $\Theta'_1$ and $\Theta'_0$ which fulfill the following goals:

i) If $x(n)$ is stationary, $\Theta'_1 \approx \Theta_1$ and $\Theta'_0 \approx \Theta_0$.

ii) If $x(n)$ is nonstationary, $\Theta'_1 > \Theta_1$ and $\Theta'_0 \approx \Theta_0$.

Note that $\Theta'_0 \approx \Theta_0$ is a requirement of both goals i) and ii). Since the test threshold $\gamma$ in (11) is derived from the gamma PDF fitted to $\Theta_0$, such a requirement tells that $\gamma$ should remain approximately the same for $H'(x)$ (i.e., $\gamma' \approx \gamma$), be $x(n)$ stationary or nonstationary. Thus, if goal ii) holds, the nonstationarity detection rate of $H'(x)$ should be higher than the one of $H(x)$, since $\Theta'_1 > \Theta_1$ and $\gamma' \approx \gamma$. Conversely, if goal i) holds, the false positive rates of $H'(x)$ and $H(x)$ should be about the same, as $\Theta'_1 \approx \Theta_1$ and $\gamma' \approx \gamma$. Based on these criteria, if the two goals are fulfilled, $H'(x)$ should have better classification performances in comparison to $H(x)$. Having defined how the original test can be improved, we now present a way to obtain the modified parameters $\Theta'_1$ and $\Theta'_0$.

#### B. Obtaining $\Theta'_1$ and $\Theta'_0$

Let $\lambda(n)$ be a hypothetical time-varying feature extracted from $S(n, f)$, i.e., $\lambda(n) = f[S(n, f)]$. Consider one can choose any functional form $f(\cdot)$ for $\lambda(n)$, provided the chosen function (i.e., expression) exhibit the following properties:

P.1) If $x(n)$ is nonstationary, $\mathrm{Var}[\lambda(n)] > 0$ and $\mathrm{E}[\lambda(n)] \ge 1$.

P.2) If $x(n)$ is stationary, $\lambda(n) \approx 1$.

P.3) $\lambda(n)$ and individual values of $d(n, f)$ [see (5)] can be considered as approximately independent.

By using $\lambda(n)$, we propose to modify $\Theta_1$ and $\Theta_0$ as follows:

$$ \Theta'_1 = \mathrm{Var}[\lambda(n)c(n)] \qquad (12) $$

$$ \Theta'_0 = \{\mathrm{Var}[\lambda_j(n)c_j(n)]\}, j = 1, ..., J \qquad (13) $$

<footnote n="1">¹For these distances, for instance, (3) is given by $d(n, f) = [S(n, f) - \bar{S}(f)] \log[S(n, f)/\bar{S}(f)]$, $d(n, f) = S(n, f)/\bar{S}(f) - \log[S(n, f)/\bar{S}(f)] - 1$, and $d(n, f) = |\log[S(n, f)/\bar{S}(f)]|$, respectively.</footnote>

<footnote n="2">²Here, stationarity is being tested for the whole signal segment [21].</footnote>
where $\lambda_j(n)$ is the feature corresponding to the $j^{th}$ surrogate. It is assumed that the properties above hold for $\lambda_j(n)$. Thus, to obtain $\Theta'_1$ and $\Theta'_0$, we use the values of $\lambda(n)$ and $\lambda_j(n)$ as weights<fnref n="3" /> for $c(n)$ and $c_j(n)$, respectively. With (12) and (13) at hand, we show next how $\Theta'_1$ and $\Theta'_0$ fulfill goals i) and ii).

### C. How $\Theta'_1$ and $\Theta'_0$ fulfill goals i) and ii)

As pointed out in Section III-A, condition $\Theta'_0 \approx \Theta_0$ appears for goals i) and ii). We start by showing that such a condition holds. To do so, we rewrite (13) using $c_j(n)$ as given in (9)

$$ \Theta'_0 = \mathrm{Var} \left[ \sum_f y(n) d_j(n, f) \right]. \quad (14) $$

By using the variance of sum law [29], (14) is expanded to

$$ \Theta'_0 = \left\{ \sum_f \mathrm{Var} \left[ \lambda_j(n) d_j(n, f) \right] + 2 \sum_{f_i < f} \mathrm{Cov} \left[ \lambda_j(n) d_j(f_i, n), \lambda_j(n) d_j(f_i, n) \right] \right\}, j = 1, \dots, J, \quad (15) $$

where $\mathrm{Cov}[\cdot, \cdot]$ is the covariance. As $\lambda_j(n)$ and $d_j(n, f)$ are independent (see P.3), the first summand in (15) is equal to

$$ \mathrm{Var}[\lambda_j(n)d_j(n, f)] = \mathrm{Var}[d_j(n, f)]\mathrm{Var}[\lambda_j(n)] + \mathrm{Var}[\lambda_j(n)]\mathrm{E}^2[d_j(n, f)] + \mathrm{Var}[d_j(n, f)]\mathrm{E}^2[\lambda_j(n)] \quad (16) $$

due to the product property of the variance [29]. One can also show that the second summand in (15) can be rewritten as<fnref n="4" />

$$ \begin{aligned} \mathrm{Cov}[\lambda_j(n)d_j(f_i, n), \lambda_j(n)d_j(f_i, n)] = & \mathrm{Cov}[d_j(n, f_i), d_j(n, f_i)]\{\mathrm{Var}[\lambda_j(n)] + \mathrm{E}^2[\lambda_j(n)]\} \\ & + \mathrm{Var}[\lambda_j(n)]\mathrm{E}[d_j(n, f_i)]\mathrm{E}[d_j(n, f_i)]. \end{aligned} \quad (17) $$

By substituting (16) and (17) into (15), rearranging the resulting expression, and making use of (9), one gets the following expression for $\Theta'_0$:

$$ \Theta'_0 = \{\mathrm{Var}[\lambda_j(n)]\mathrm{E}^2[c_j(n)] + \mathrm{Var}[c_j(n)]\}\{\mathrm{Var}[\lambda_j(n)] + \mathrm{E}^2[\lambda_j(n)]\}, \quad j = 1, \dots, J. \quad (18) $$

Since all $j = 1, \dots, J$ surrogates are stationary by definition, $\lambda_j(n) \approx 1 \ \forall j$ (see property P.2), which implies that $\mathrm{Var}[\lambda_j(n)] \approx 0$ and $\mathrm{E}^2[\lambda_j(n)] \approx 1 \ \forall j$. Owing to this fact, it can be shown that (18) reduces to

$$ \Theta'_0 \approx \mathrm{Var}[c_j(n)], \quad j = 1, \dots, J. \quad (19) $$

Then, by comparing (19) with (10), it can be seen that

$$ \Theta'_0 \approx \Theta_0. \quad (20) $$

Note that (20) is valid if $x(n)$ stationary or nonstationary, requiring only properties P.2 and P.3, and the stationarity property of the surrogates to hold. Now, following the steps taken in (15)–(18) and using (6), one can rewrite (12) as

$$ \Theta'_1 = \mathrm{Var}[\lambda(n)]\mathrm{E}^2[c(n)] + \{\mathrm{Var}[\lambda(n)] + \mathrm{E}^2[\lambda(n)]\}\Theta_1. \quad (21) $$

If $x(n)$ is nonstationary, $\mathrm{Var}[\lambda(n)] > 0$ and $\mathrm{E}^2[\lambda(n)] \ge 1$ (see property P.1). Since all terms in (21) are positive, we have

$$ \Theta'_1 > \Theta_1. \quad (22) $$

As (20) and (22) hold true if $x(n)$ is nonstationary, we conclude that goal ii) is fulfilled. Now, if $x(n)$ is stationary, $\lambda \approx 1$, which yields $\mathrm{Var}[\lambda(n)] \approx 0$ and $\mathrm{E}^2[\lambda(n)] \approx 1$ (see property P.2). Then, (21) reduces to

$$ \Theta'_1 \approx \Theta_1 \quad (23) $$

and since (20) and (23) are verified in case $x(n)$ is stationary, goal i) is also achieved. Having shown that the modified test parameters in (12) and (13) fulfill goals i) and ii), we discuss in the next section the choice of $\lambda(n)$.

### D. Choosing a candidate $\lambda(n)$

The previous discussion provides general guidelines (for this work and possibly future ones) on how to improve the surrogate test by searching for a time-varying function $\lambda(n)$ that fulfills properties P.1 to P.3. Although there is no unique choice for $\lambda(n)$, note that the greater $\mathrm{Var}[\lambda(n)]$ is when $x(n)$ is nonstationary, the farther $\Theta'_1$ is from $\Theta_1$ [see (21)], and the more likely the test will (correctly) reject stationarity. By considering these points, we propose below a candidate $\lambda(n)$ that can be extracted from the TF representation of $x(n)$.

The center of gravity (or the first-order moment) of the local spectra (1) at the $n^{th}$ time point can be computed as

$$ f_m(n) = \frac{\sum_f f S(n, f)}{\sum_f S(n, f)}. \quad (24) $$

Due to (7), $f_m(n)$ should be approximately constant if $x(n)$ is stationary, but a time-varying function otherwise. Note that $f_m(n)$ can be seen as an approximation of the instantaneous frequency (IF) for monocomponent signals [4]. However, it is not the aim of this paper to interpret the physical meaning of $f_m(n)$, but to take advantage of the fact that (24) behaves differently in case of stationarity/nonstationarity, and captures the collective behavior of the spectra at time $n$. With $f_m(n)$ at hand, we propose to compute $\lambda(n)$ as

$$ \lambda(n) = 1 + \Delta f_m(n) \quad (25) $$

where $\Delta = \|f_m(n) - \bar{f}_m\|$ with $\| \cdot \|$ and $\bar{f}_m$ being the $l^2$ norm and the median value of $f_m(n)$, respectively. It can be shown that property P.1 holds for (25); since $f_m(n) > 0$, one has $\Delta \ge 0$, which implies that $\mathrm{E}[\lambda(n)] \ge 1$. Furthermore, if $x(n)$ is nonstationary $f_m(n)$ should fluctuate over time [and so $\lambda(n)$], yielding $\mathrm{Var}[\lambda(n)] > 0$. Note that, the wider the spread of $f_m(n)$, the larger the value of $\Delta$, and thus the greater $\mathrm{Var}[\lambda(n)]$. Property P.2 is also fulfilled by (25); since $f_m(n)$ should be approximately constant in case of stationarity one gets $\Delta \approx 0$, which yields $\lambda(n) \approx 1$. Finally, since (24) is computed only from the instantaneous spectrum at time $n$ (i.e., independently of the local spectra at other time points), while $d(n, f)$ requires the knowledge of the local spectra at all time instants [see (4) and (5)], one can consider that property P.3 should hold as well. By taking into account the procedure

<footnote n="3">³Note that setting $\lambda(n) = 1$ in (12) makes $\Theta'_1$ reduce to the original $\Theta'_1$ of [21]. In turn, if $\lambda(n) = \sum_f S(n, f)$ in (12), $\Theta'_1$ becomes the modified test statistic of [24]. Thus, the present method generalizes [21] and [24].</footnote>

<footnote n="4">⁴Steps to obtain (17) are shown in the supplementary material of this letter.</footnote>
<table><caption>TABLE I
PERCENTAGE OF "NONSTATIONARY" OUTCOMES OBTAINED BY TESTING 1000 MC SIMULATIONS OF THE SYNTHETIC SIGNALS.</caption><thead><tr><th rowspan="2">Test signals</th><th colspan="4">N = 250</th><th colspan="4">N = 500</th><th colspan="4">N = 750</th><th colspan="4">N = 1000</th></tr><tr><th>SUR<sub>new</sub></th><th>SUR<sub>ori</sub></th><th>SUR<sub>alt</sub></th><th>KPSS</th><th>SUR<sub>new</sub></th><th>SUR<sub>ori</sub></th><th>SUR<sub>alt</sub></th><th>KPSS</th><th>SUR<sub>new</sub></th><th>SUR<sub>ori</sub></th><th>SUR<sub>alt</sub></th><th>KPSS</th><th>SUR<sub>new</sub></th><th>SUR<sub>ori</sub></th><th>SUR<sub>alt</sub></th><th>KPSS</th></tr></thead><tbody><tr><td>Gauss. mean</td><td><b>12%</b></td><td>5.8%</td><td>4%</td><td>64.3%</td><td><b>16.2%</b></td><td>6.2%</td><td>5%</td><td>95.2%</td><td><b>18.5%</b></td><td>5.2%</td><td>5.3%</td><td>98.7%</td><td><b>22.4%</b></td><td>7.3%</td><td>8.6%</td><td>100%</td></tr><tr><td>Gauss. var</td><td><b>100%</b></td><td>100%</td><td>45%</td><td>6.1%</td><td><b>100%</b></td><td>100%</td><td>47.2%</td><td>6.3%</td><td><b>100%</b></td><td>100%</td><td>49.4%</td><td>6.9%</td><td><b>100%</b></td><td>100%</td><td>55.7%</td><td>7%</td></tr><tr><td>TVAR</td><td><b>32.4%</b></td><td>9%</td><td>14%</td><td>0.2%</td><td><b>50.3%</b></td><td>11.3%</td><td>22.2%</td><td>0.5%</td><td><b>65.4%</b></td><td>12%</td><td>24.8%</td><td>0.4%</td><td><b>80.2%</b></td><td>11.1%</td><td>17%</td><td>1.1%</td></tr><tr><td>UMP</td><td><b>69.3%</b></td><td>91%</td><td>0.4%</td><td>0.7%</td><td><b>99.8%</b></td><td>100%</td><td>0.8%</td><td>0.6%</td><td><b>100%</b></td><td>100%</td><td>1.1%</td><td>1.3%</td><td><b>100%</b></td><td>100%</td><td>2%</td><td>1.2%</td></tr><tr><td>WGN</td><td><b>6.7%</b></td><td>7.5%</td><td>2%</td><td>4.1%</td><td><b>6.1%</b></td><td>6.5%</td><td>0.1%</td><td>5.8%</td><td><b>1.4%</b></td><td>6%</td><td>1.5%</td><td>6.8%</td><td><b>5.4%</b></td><td>6.1%</td><td>1.9%</td><td>6.1%</td></tr><tr><td>AR(1)</td><td><b>6.5%</b></td><td>6.2%</td><td>2.5%</td><td>0.2%</td><td><b>5.9%</b></td><td>5.7%</td><td>2.5%</td><td>0.1%</td><td><b>5.2%</b></td><td>6.3%</td><td>2.8%</td><td>0%</td><td><b>5.1%</b></td><td>6.9%</td><td>2.7%</td><td>0%</td></tr></tbody></table>

discussed in Section III-C and $\lambda(n)$ proposed in this section, we now evaluate the performance of the modified test.

## IV. EXPERIMENTAL STUDY

Stationarity tests have been applied to uniformly modulated processes (UMP), time-varying autoregressive (TVAR) processes, and nonstationary Gaussian time series. These processes are often employed to model nonstationarities in practice. More details on these signals are given below. Stationary first-order autoregressive [AR(1)] and white Gaussian noise [WGN(0,1)] processes have been probed as well. The test signals have 250, 500, 700, and 1000 data points.

We have also tested a real-world dataset, which is the monthly average Central England Temperature (CET) time series [30]–[35]. The experimental study on the real-world data is presented in the supplementary material of this paper.

### A. Tested signals

1) UMP: In this paper, the following UMP has been tested:

$$x(n) = C n e^{\beta n/2} v_u(n) \qquad (26)$$

where $v_u(n) \sim U[-1, 1]$ is uniform random noise. The UMP in (26) has been used to model ground motion in earthquakes, with parameters $C$ and $\beta$ being the intensity of the ground acceleration, and the effective duration of the motion [36]. Here, we have set $\beta = 2/\rho N$ and $C = \beta/2e^{-1}$, with $\rho = 0.75$.

2) TVAR processes: The considered process is given as

$$x(n+1) = x(n) \sum_{m=0}^{q} a_m f_m(n) + \phi(n) \qquad (27)$$

where $\phi(n) \sim \text{WGN}(0, 1)$, $\{a_0, ..., a_q\}$ are constants, and $\{f_0(n), ..., f_q(n)\}$ are pre-determined basis functions chosen to characterize a given nonstationary behavior [37]. Here, we have considered $q = 2$, $a_0 = 1$, $a_1 = 2$, $a_2 = 1/2$, $f_0(n) = 1$, $f_1(n) = (n-1)/(N-2)$ and $f_2(n) = 3f_1^2(n)-1$ [38]. These basis functions are Legendre polynomials and have been used to model slowly-varying nonstationary behaviors [39].

3) Nonstationary Gaussian processes: These signals follow an AB model as in [40], [41], which means that their PDFs have a time-varying parameter $\theta(n)$ such as

$$\theta(n) = \begin{cases} \xi_1, & n = 1, ..., N/2, \\ \xi_2, & n = N/2 + 1, ..., N \end{cases} \qquad (28)$$

where $\xi_1$ and $\xi_2$ are the parameter values before and after change, respectively. Note that (28) describes a step change at $n = N/2$. The following Gaussian processes have been created using (28): a) signals with fixed mean ($\mu = 0$) and variance varying from $\xi_1 = 1$ to $\xi_1 = 4$, b) signals with fixed variance ($\sigma = 1$) and mean varying from $\xi_1 = 0$ to $\xi_1 = 4$.

### B. Test results

We have tested 1000 Monte Carlo (MC) simulations of the test signals with the new method and competing ones. The results are shown in Table I as percentage of cases each signal has been classified as “nonstationary”. In Table I, “SUR<sub>new</sub>” is the new surrogate test, “SUR<sub>ori</sub>” is the original surrogate test of [21], “SUR<sub>alt</sub>” is the alternative version of [21] proposed in [24], and “KPSS” is the classical KPSS test [1]. All tests have been performed with a significance level of 5%. For the surrogate-based method, the free parameter $n_{h0}$ has been set to $n_{h0} = 0.3$, the number of surrogates used has been $J = 100$, and (6) has been chosen as the log-spectral deviation [26].

As mentioned in Section I, the surrogate test works well for second-order nonstationary and amplitude-modulated signals (like the Gaussian processes with varying variance and the UMPs), but not for first-order and slowly nonstationary processes (like the Gaussian signals with varying mean and the TVAR processes). Detecting the latter signals is problematic for TF-based methods like the surrogate test. In addition to the fact that quadratic TF representations are better suited to characterize second-order evolutions [4], they often estimate poorly the spectral content at very low frequencies [5]. This is a problem for slowly nonstationary processes, whose spectral content is often concentrated at low frequency bands.

Note in Table I that, despite of the aforementioned limitations, the detection rates for the challenging signals (Gaussian mean and TVAR) at least doubled for SUR<sub>new</sub> in comparison to SUR<sub>ori</sub> and SUR<sub>alt</sub>. Nevertheless, SUR<sub>new</sub> has kept the same excellent performances as SUR<sub>ori</sub> for the best cases (Gaussian variance and UMP). The KPSS test, on the other hand, has performed best only for first-order nonstationary and AR(1) signals, which is reasonable given the way the KPSS test is designed [1]. By considering the tests as binary classifiers and computing the Accuracy (ACC) values, we obtain

$$\text{ACC}_{\text{new}} = 0.76, \quad \text{ACC}_{\text{ori}} = 0.66, \quad \text{ACC}_{\text{alt}} = 0.45, \quad \text{ACC}_{\text{KPSS}} = 0.48$$

which evaluate the overall classification performance of the tests. Note that the ACC of the new method is highest one.

## V. CONCLUDING REMARKS

This paper revisited the methodology for testing stationarity with surrogates. We studied a general weighting method to improve the test by transforming the original test parameters. A new time-varying feature with some special characteristics was designed by using the TF representation of the signal. We showed mathematically that the proposed procedure to modify the test can improve its classification performance. Simulations were carried out on synthetic and real-world signals to evaluate the performances of the modified test and competing ones.
REFERENCES

[1] D. Kwiatkowski, P. C. B. Phillips, P. Schmidt, and Y. Shin, "Testing the null hypothesis of stationarity against the alternative of a unit root: How sure are we that economic time series have a unit root," *J. Economet.*, vol. 54, pp. 159-178, Oct. 1992.

[2] H. Lins and T. Cohn, "Stationarity: Wanted dead or alive?" *J. Am. Water Resour. Assoc.*, vol. 47, no. 3, pp. 475-480, Jun. 2011.

[3] D. G. Manolakis, V. K. Ingle, and S. M. Kogon, *Statistical and Adaptive Signal Processing: Spectral estimation, Signal Modeling, Adaptive Filtering and Array Processing*. London: Artech House, 2005.

[4] P. Flandrin, *Time-Frequency/Time-Scale Analysis*, ser. Wavelet analysis and its applications. San Diego, CA: Academic Press, 1999.

[5] B. Boashash, *Time Frequency Signal Analysis and Processing: A Comprehensive Reference*, 1st ed. Oxford, UK: Elsevier, 2003.

[6] I. Shafi et al., "Techniques to obtain good resolution and concentrated time-frequency distributions: a review," *EURASIP J. Adv. Signal Process.*, vol. 2009, no. 1, pp. 1-43, Jun. 2009.

[7] Z. Leonowicz, *Parametric methods for time-frequency analysis of electric signals*. Poland: Faculty of Electrical Engineering, Wroclaw University of Technology, 2006.

[8] Y. Yang, Z. Peng, G. Meng, and W. Zhang, "Characterize highly oscillating frequency modulation using generalized warblet transform," *Mech. Syst. Signal Process.*, vol. 26, pp. 128-140, Jan. 2012.

[9] M. Bayram and R. Baraniuk, *Nonlinear and Nonstationary Signal Processing*. Cambridge Univ. Press, 2001, ch. Multiple Window Time-Varying Spectrum Estimation, pp. 292-316.

[10] M. Alam, M. Rahman, N. Parvin, and M. Sobhan, "Time-frequency representation of a signal through non-stationary multipath fading channel," in *Proc. Int. Conf. on Informatics, Electronics and Vision (ICIEV)*, Dhaka, Bangladesh, May 2012, pp. 1130-1135.

[11] A.-K. Seifert, A. M. Zoubir, and M. G. Amin, "Radar classification of human gait abnormality based on sum-of-harmonics analysis," in *Proc. IEEE Radar Conf.*, Oklahoma City, OK, Apr. 2018, pp. 940-945.

[12] P. Basu, D. Rudoy, and P. J. Wolfe, "A nonparametric test for stationarity based on local Fourier analysis," in *Proc. IEEE Int. Conf. Acoust., Speech, Signal Process. (ICASSP)*, Taiwan, Apr. 2009, pp. 3005-3008.

[13] D. Baptista de Souza, J. Chanussot, A.-C. Favre, and P. Borgnat, "A nonparametric test for slowly-varying nonstationarities," *Signal Process.*, vol. 143, pp. 241-252, Feb. 2018.

[14] Y. Xiang, J. Ding, and V. Tarokh, "Estimation of the evolutionary spectra with application to stationarity test," *IEEE Trans. on Signal Process.*, vol. 67, no. 5, pp. 1353-1365, Mar. 2019.

[15] M. A. Khan and J. W. Pierret, "A test for non-stationarity of synchrophasor measurements," in *Proc. IEEE Int. Conf. on Probabilistic Methods Applied to Power Systems (PMAPS)*, Boise, ID, Jun. 2018, pp. 1-6.

[16] Y. Xiang, J. Ding, and V. Tarokh, "Evolutionary spectra based on the multitaper method with application to stationarity test," in *Proc. IEEE Int. Conf. Acoust., Speech, Signal Process. (ICASSP)*, Calgary, Canada, Apr. 2018, pp. 3994-3998.

[17] G. Matz and F. Hlawatsch, "Wigner distributions (nearly) everywhere: time-frequency analysis of signals, systems, random processes, signal spaces and frames," *Signal Process.*, vol. 83, no. 7, pp. 1355-1378, Jul. 2003.

[18] R. Iqbal, T. Abhayapala, J. Ahmed, and T. Lamahewa, "Wigner-ville distribution of a type of non-stationary mobile rayleigh fading channels," in *Proc. IEEE Int. Multi-topic Conf. (INMIC)*, Islamabad, Pakistan, Dec. 2009, pp. 1-6.

[19] X. Chen and Z. Feng, "Iterative generalized time-frequency reassignment for planetary gearbox fault diagnosis under nonstationary conditions," *Mech. Syst. Signal Process.*, vol. 80, pp. 429-444, Dec. 2016.

[20] J. Bian et al., "A non-stationary mimo channel model for street corner scenarios considering velocity variations of the mobile station and scatterers," in *IEEE/CIC Int. Conf. on Comm. in China (ICCC)*, Qingdao, China, Oct. 2017, pp. 1-6.

[21] P. Borgnat, P. Flandrin, P. Honeine, C. Richard, and J. Xiao, "Testing stationarity with surrogates: A time-frequency approach," *IEEE Trans. Signal Process.*, vol. 58, no. 7, pp. 3459-3470, Jul. 2010.

[22] T. Schreiber and A. Schmitz, "Surrogate time series," *Physica D*, vol. 142, no. 3, pp. 346-382, 2000.

[23] J. Theiler, S. Eubank, A. Longtin, B. Goldrikian, and J. D. Farmer, "Testing for nonlinearity in time series: the method of surrogate data," *Physica D*, vol. 58, no. 1, pp. 77-94, Sep. 1992.

[24] D. Baptista de Souza, J. Chanussot, A.-C. Favre, and P. Borgnat, "A modified time-frequency method for testing wide-sense stationarity," in *Proc. IEEE Int. Conf. Acoust., Speech, Signal Process. (ICASSP)*, Kyoto, Japan, Mar. 2012, pp. 3409-3412.

[25] P. Flandrin, "Maximum signal energy concentration in a time-frequency domain," in *Proc. IEEE Int. Conf. Acoust., Speech, Signal Process. (ICASSP)*, vol. 4, New York, NY, Apr. 1988, pp. 2176-2179.

[26] M. Basseville, "Distances measures for signal processing and pattern recognition," *Signal Process.*, vol. 18, no. 4, pp. 349-369, Dec. 1989.

[27] J. Xiao, P. Borgnat, and P. Flandrin, "Testing stationarity with time-frequency surrogates," in *Proc. Eur. Signal Process. Conf. (EUSIPCO)*, Poznan, Poland, Sep. 2007, pp. 2020-2024.

[28] J. Xiao, "Contributions to nonstationary spectrum estimation and stationarity tests in the time-frequency plane," Ph.D. dissertation, École Normale Supérieure de Lyon, 2008.

[29] D. Blumenfeld, *Operations research calculations handbook*. Boca Raton, NJ: CRC Press LLC, 2001.

[30] Met Office Hadley Centre, "Central england temperature data," [Online]. Available: https://www.metoffice.gov.uk/hadobs/hadcet/data/download.html, Accessed on: Dec. 10, 2018.

[31] D. J. Thomson, "The seasons, global temperature, and precession," *Science*, vol. 268, no. 5207, pp. 59-68, Apr. 1995.

[32] T. Proietti, "Seasonal changes in central england temperatures," *J. R. Statist. Soc. A*, vol. 180, no. 3, pp. 769-791, May. 2017.

[33] T. Zhang and W. B. Wu, "Testing parametric assumptions of trends of nonstationary time series," *Biometrika*, vol. 98, no. 3, pp. 599-614, Jul. 2011.

[34] J. de Souza, S. M. D. Queirós, and A. M. Grimm, "Components of multifractality in the central england temperature anomaly series," *Chaos*, vol. 23, pp. 1-10, Jun. 2013.

[35] T. C. Benner, "Central england temperatures: Long-term variability and teleconnections," *Int. J. Climatol.*, vol. 19, pp. 391-403, Apr. 1999.

[36] A. Giaralis and P. D. Spanos, "Derivation of non-stationary stochastic process compatible with seismic response/design spectra," in *Proc. Int. Conf. on Comp. Stochastic Mechanics*, Rhodes, Greece, Jun. 2010, pp. 1-13.

[37] S. Kay, "A new nonstationarity detector," *IEEE Trans. Signal Process.*, vol. 56, no. 4, pp. 1440-1451, Apr. 2008.

[38] G. R. S. Reddy and R. Rao, "Performance analysis of basis functions in TVAR model," *Int. J. of Signal Processing, Image Processing and Pattern Recognition*, vol. 17, no. 3, pp. 317-338, Jun. 2014.

[39] Y. Li, Q. Liu, S.-R. Tan, and R. H. M. Chan, "High-resolution time-frequency analysis of EEG signals using multiscale radial basis functions," *Neurocomputing*, vol. 195, pp. 96-103, Jun. 2016.

[40] F. Lombard, "Rank test for changepoint problems," *Biometrika*, vol. 74, no. 3, pp. 615-624, Sep. 1987.

[41] J.-F. Quessy, A.-C. Favre, M. Saïd, and M. Champagne, "Statistical inference in Lombard's smooth-change model," *Environmetrics*, vol. 22, no. 7, pp. 882-893, Nov. 2011.