Adaptive State-Space Multitaper Spectral Estimation

Nonstationary time series with time-varying probability structures are ubiquitous. Some examples include speech and image recordings [1, 2], oceanographic and seismic signals [3], neural spike trains [4], and electroencephalogram (EEG) [5]. We are interested in analyzing the nonstationary data through the lens of the time-varying spectral dynamics, which yields valuable information on the underlying system. The traditional approach has been to segment the data into independent overlapping or non-overlapping windows, assuming local stationarity [6] within each window, and to apply Fourier or wavelet transform [7, 8].

Despite the popularity, the windowing approach suffers from spectral estimates of high variance within each window, due to finite window-length [9] and the restrictive independence assumption for different windows. The state-space multitaper (SSMT) framework [10] and its extensions [11, 12, 13, 14] have been proposed as the solutions, by positing a time-invariant latent state-space model in the time-frequency domain, with each state representing the Fourier coefficients in each window. Since the states are linked by stochastic continuity constraint across windows, the spectral estimates are not independent.

However, the use of time-invariant parameters limits the capacity of SSMT to track strong nonstationarity that commonly occurs in systems neuroscience and manifests as strong spectral fluctuations due to fluctuations in system properties. Although time-varying state-space model offers an obvious solution, traditional parameter estimation approaches based on expectation-maximization (EM) algorithm [15] makes it less practical for real-time applications - The parameters must be re-estimated with every incoming batch of data, and thus incurs high computational costs.

We propose a time-varying extension of SSMT with an efficient and adaptive parameter estimation scheme, termed an adaptive SSMT (ASSMT) method. Based on a novel nonstationarity metric derived from spectral fluctuation, ASSMT adaptively adjusts the state parameters and consequently Kalman gain. This allows ASSMT to switch between removing background noise and tracking spectral fluctuation in a data-driven manner. The estimation for time-varying parameters requires only single pass through the data, making it well suited for real-time applications. The paper is organized as follows. In Section II, we review the SSMT framework. In Section III, we formulate the ASSMT framework. In Section IV and V, the results and conclusion are presented.

We first review the SSMT algorithm in [10]. Consider a nonstationary time series $y_{t}$ sampled at frequency $f_{s}$ as

where $x_{t}$ is a locally stationary latent Gaussian process [16] with the measurement noise $\varepsilon_{t}\sim\mathcal{N}(0,\sigma_{\varepsilon}^{2})$. Leveraging the local stationary property, we divide these signals into $K$ nonoverlapping stationary intervals of $J$ samples, such that $T=KJ$. The segmented vectors for interval $k$ are denoted as $Y_{k},X_{k},\varepsilon_{k}\in\mathbb{R}^{J}$, with the $j$th element as $Y_{k,j}=y_{J(k-1)+j}$ for $k=1,\ldots,K$ and $j=1,\ldots,J$.

To perform time-frequency analysis, we introduce the latent $Z_{k}=(Z_{k,1},\ldots,Z_{k,J})\in\mathbb{C}^{J}$, where $Z_{k,j}$ is a complex Gaussian variable with the magnitude corresponding to the power at the normalized frequency $\omega_{j}=2\pi(j-1)/J$ and interval $k$ [17]. To model the evolution of spectra across the windows, we assume that $\{Z_{k}\}_{k=1}^{K}$ follow a random walk prior

where the state noise $v_{k}$ is a complex Gaussian process with a diagonal covariance $I(\sigma^{2}_{v,j})$, i.e., $v_{k}\sim\mathcal{CN}(0,I(\sigma^{2}_{v,j}))$. The prior encodes two important properties of $\{Z_{k}\}_{k=1}^{K}$. First, the state variance $\sigma_{v,j}^{2}$ controls the smoothness of the process, with a large value indicative of non-smooth or fluctuating process. Second, $Z_{k,j}$ and $Z_{k,j^{\prime}}$ for $j\neq j^{\prime}$ are independent a priori.

We can link the time-frequency process $\{Z_{k}\}_{k=1}^{K}$ with time-domain observation $\{Y_{k}\}_{k=1}^{K}$, using the Fourier transform matrix $F\in\mathbb{C}^{J\times J}$ with $F_{j,l}=J^{-1/2}\exp(-i2\pi(l-1)j/J)$ and the inverse Fourier matrix $W\in\mathbb{C}^{J\times J}$, such that $FW=I$,

where $\varepsilon_{k}^{F}\sim\mathcal{CN}(0,I(\sigma_{\varepsilon}^{2}))$, since $FI(\sigma_{\varepsilon}^{2})W=I(\sigma_{\varepsilon}^{2})$.

Along with the state-space model, we incorporate the data tapers to further reduce the variance of the spectral estimates. Specifically, we use $M$ Slepian tapers, leading to the multitaper (MT) method that optimally balances the bias-variance trade-off via bandwidth adjustment [18, 19]. This essentially produces $M$ independent set of state-space models,

where $v_{k}^{(m)}\sim\mathcal{CN}(0,I(\sigma^{2,(m)}_{v,j}))$, $\varepsilon_{k}^{(m),F}\sim\mathcal{CN}(0,I(\sigma_{\varepsilon}^{2,(m)}))$, $Y_{k}^{(m)}$ denotes the $m^{\text{th}}$ Slepian taper applied to $Y_{k}$, $Y_{k}^{(m),F}$ denotes Fourier transform of $Y_{k}^{(m)}$, i.e., $Y_{k}^{(m),F}=FY_{k}^{(m)}$, and $Z^{(m)}_{k}$ represents the $m$th spectral eigen-coefficient of $Z_{k}$. The variants of SSMT modify Eqs. (4) or (5) [11, 12, 13, 14].

We now focus on $\omega_{j}$ for simplicity. Based on Eqs. (4) and (5), we derive a Kalman filter algorithm for the estimation of $Z^{(m)}_{k,j}$ using Kalman gain $C^{(m)}_{k,j}$

where the notation $k|s$ denotes the estimate on interval $k$ given the data observed up to interval $s$ and $C^{(m)}_{k,j}$ is given as

The spectrogram estimate at frequency $\omega_{j}$ on interval $k$ is

The parameters $\{\sigma_{v,j}^{2,(m)}\}_{j,m=1}^{J,M}$ and $\{\sigma_{\varepsilon}^{2,(m)}\}_{m=1}^{M}$ are estimated with EM algorithm [15]. SSMT tracks nonstationarity in the data with the time-invariant parameters $\sigma^{2,(m)}_{v,j}$ and $\sigma_{\varepsilon}^{2,(m)}$, since $Z^{(m)}_{k|k,j}\neq Z^{(m)}_{k-1|k-1,j}$ and thus $f^{\text{SSMT}}_{k}(\omega_{j})\neq f^{\text{SSMT}}_{k-1}(\omega_{j})$.

Although SSMT can model slowly time-varying spectral dynamics, it is restrictive for highly fluctuating spectral dynamics. Such strong nonstationarity (or fluctuation) is common in EEG with external intervention to the brain or with the change of the brain state during anesthesia/sleep [20]. Fig. 1 shows a snippet of human anesthesia EEG in which SSMT (Fig. 1(b)) cannot track the apparent dynamics shown by the MT approach (Fig. 1(a)). However, the proposed ASSMT (Fig. 1(c)) is able to capture the spectral dynamics and remove non-relevant background activity.

The failure of SSMT is due to the time-invariant parameters, as the Kalman gain quickly converges to a steady-state value $C_{\infty,j}^{(m)}$ [21]. Denoting the observation prediction error as $\Delta Y^{(m)}_{k,j}=Y^{(m),F}_{k,j}-Z_{k-1|k-1,j}^{(m),F}$ and the latent estimate change as $\Delta Z_{k,j}^{(m)}=Z^{(m)}_{k|k,j}-Z^{(m)}_{k-1|k-1,j}$, we can express Eq. (6) as $\Delta Z_{k,j}^{(m)}=C^{(m)}_{\infty,j}\Delta Y_{k,j}^{(m)}$. If $C_{\infty,j}^{(m)}$ is low (SSMT for Fig. 1(b)), $Z^{(m)}_{k|k,j}$ fails to reflect fluctuation in $Y_{k,j}^{(m),F}$. Consequently, $f_{k}^{\text{SSMT}}$ cannot reliably track the spectral dynamics.

To resolve this issue, adaptive SSMT (ASSMT) posits a time-varying state-space model, to allow the adaptive change of $\sigma_{v,j}^{2,(m)}$. Specifically, we replace $v_{k}^{(m)}\sim\mathcal{CN}(0,I(\sigma^{2,(m)}_{v,j}))$ in Eq. (5) with $v^{(m)}_{k,j}\sim\mathcal{CN}(0,I(\sigma^{2,(m)}_{v,k,j}))$, to indicate the state variance’s dependence on time. This prevents Kalman gain from converging to a steady-state value, allowing the algorithm to produce more flexible spectrogram estimates. $\sigma_{\varepsilon}^{2,(m)}$ remains constant across windows because we assume stationary background noise.

With the modified generative model, we now address when and how ASSMT tracks varying degrees of nonstationarity. We first quantify the notion of nonstationarity, and then propose an adaptive parameter estimation approach.

We apply ASSMT to two datasets: 1) nonstationary simulated data and 2) EEG data from a patient under propofol anesthesia. We compare the spectrogram estimates between MT, SSMT, and ASSMT. All spectrograms are in the dB scale.

We introduced an adaptive state-space multitaper (ASSMT) framework for estimating spectral dynamics in the nonstationary time series. By relaxing the time-invariant parameters assumption and proposing an adaptive parameter estimation scheme, we demonstrated that ASSMT was able to capture strong power fluctuations more reliably compared to SSMT.

We apply ASSMT to EEG recording of another volunteer receiving propofol anesthesia. All parameter settings of SSMT and ASSMT are same as those used in the main paper.