Dynamic Bayesian Multitaper Spectral Analysis

TL;DR

This paper introduces a Bayesian state-space approach to spectral analysis that improves upon traditional sliding window methods by providing more accurate, high-resolution spectrograms with confidence intervals, especially in noisy conditions.

Contribution

It combines multitaper spectral analysis with Bayesian state-space modeling and EM algorithms to better capture spectral dynamics and provide confidence measures.

Findings

Enhanced spectral estimates with higher resolution.

Robustness to noise demonstrated on EEG and electric network data.

Theoretical bias-variance analysis shows performance improvements.

Abstract

Spectral analysis using overlapping sliding windows is among the most widely used techniques in analyzing non-stationary time series. Although sliding window analysis is convenient to implement, the resulting estimates are sensitive to the window length and overlap size. In addition, it undermines the dynamics of the time series as the estimate associated to each window uses only the data within. Finally, the overlap between consecutive windows hinders a precise statistical assessment. In this paper, we address these shortcomings by explicitly modeling the spectral dynamics through integrating the multitaper method with state-space models in a Bayesian estimation framework. The underlying states pertaining to the eigen-spectral quantities arising in multitaper analysis are estimated using instances of the Expectation-Maximization algorithm, and are used to construct spectrograms and their respective confidence intervals. We propose two spectral estimators that are robust to noise and are able to capture spectral dynamics at high spectrotemporal resolution. We provide theoretical analysis of the bias-variance trade-off, which establishes performance gains over the standard overlapping multitaper method. We apply our algorithms to synthetic data as well as real data from human EEG and electric network frequency recordings, the results of which validate our theoretical analysis.

Figures

Equations

$$y_t={\int }_{-\frac{1}{2}}^{\frac{1}{2}}{e}^{i2\pi ft}dz(f),$$

$$y(f)={\int }_{-\frac{1}{2}}^{\frac{1}{2}}\frac{\sin W\pi (f-\zeta )}{\sin \pi (f-\zeta )}{e}^{-i2\pi (f-\zeta )\frac{W-1}{2}}dz(\zeta ).$$

$$U^{(k)}(f):=\sum _{l=0}^{W-1}{u}_l^{(k)}{e}^{-i2\pi fl}.$$

$${\hat{S}}^{(\mathsf{mt})}(f):=\frac{1}{K}\sum _{k=1}^K|{\hat{x}}^{(k)}(f){|}^2,$$

$$y_t={\int }_{-1/2}^{1/2}{x}_t(f)e^{j2\pi ft}df,$$

$${\tilde{\mathbf{y}}}_n={\mathbf{F}}_n{\mathbf{x}}_n+{\mathbf{v}}_n,$$

$${\tilde{\mathbf{y}}}_n^{(k)}={\mathbf{F}}_n{\mathbf{x}}_n^{(k)}+{\mathbf{v}}_n^{(k)},$$

$${\mathbf{x}}_n^{(k)}={\alpha }^{(k)}{\mathbf{x}}_{n-1}^{(k)}+{\mathbf{w}}_n^{(k)},$$

$$\frac{{\hat{S}}^{(k)}(f)}{S(f)}\sim \frac{{\chi }_2^2}{2},0<f<1/2,$$

$${\psi }^{(k)}(f)=s^{(k)}(f)+{\varphi }^{(k)}(f),$$

$$p(\varphi )=\frac{1}{2}\exp \left(\varphi -\frac{1}{2}\exp (\varphi )\right),$$

$${\mathbf{\psi }}_n^{(k)}={\mathbf{s}}_n^{(k)}+{\mathbf{\varphi }}_n^{(k)},$$

$${\mathbf{s}}_n^{(k)}={\theta }^{(k)}{\mathbf{s}}_{n-1}^{(k)}+{\mathbf{e}}_n^{(k)},$$

$$\underset{{\mathbf{x}}_1^{(k)},{\mathbf{x}}_2^{(k)},\cdots ,{\mathbf{x}}_N^{(k)}}{\min }$$

$$\displaystyle\!\!\!\!(\mathbf{x}_{n}^{(k)}-\alpha\mathbf{x}_{n-1}^{(k)})^{H}{\mathbf{Q}^{(k)}_{n}}^{-1}(\mathbf{x}_{n}^{(k)}-\alpha\mathbf{x}_{n-1}^{(k)})\bigg{]},$$

$$\underset{{\mathbf{s}}_1^{(k)},{\mathbf{s}}_2^{(k)},\cdots ,{\mathbf{s}}_N^{(k)}}{\min }$$

$$+$$

$$\log L(\alpha ,{\mathbf{Q}}_{1:N})$$

$$\displaystyle\!\!\!\!\!\!\!\!\!\!\!\!\!+(\mathbf{x}_{n}-\alpha\mathbf{x}_{n-1})^{H}\mathbf{Q}_{n}^{-1}(\mathbf{x}_{n}-\alpha\mathbf{x}_{n-1})\bigg{]}+{c},$$

$${\mathbf{x}}_{n|n-1}={\alpha }^{[l]}{\mathbf{x}}_{n-1|n-1}$$

$${\mathbf{\Sigma }}_{n|n-1}={{\alpha }^{[l]}}^2{\mathbf{\Sigma }}_{n-1|n-1}+{\mathbf{Q}}^{[l]}$$

$${\mathbf{K}}_n={\mathbf{\Sigma }}_{n|n-1}{\mathbf{F}}_n^H({\mathbf{F}}_n{\mathbf{\Sigma }}_{n|n-1}{\mathbf{F}}_n^H+{\sigma }^2\mathbf{I}{)}^{-1}$$

$${\mathbf{x}}_{n|n}={\mathbf{x}}_{n|n-1}+{\mathbf{K}}_n({\tilde{\mathbf{y}}}_n-{\mathbf{F}}_n{\mathbf{x}}_{n|n-1})$$

$${\mathbf{\Sigma }}_{n|n}={\mathbf{\Sigma }}_{n|n-1}-{\mathbf{K}}_n{\mathbf{F}}_n{\mathbf{\Sigma }}_{n|n-1}$$

$${\mathbf{B}}_n={\alpha }^{[l]}{\mathbf{\Sigma }}_{n|n}{\mathbf{\Sigma }}_{n+1|n}^{-1}$$

$${\mathbf{x}}_{n|N}={\mathbf{x}}_{n|n}+{\mathbf{B}}_n({\mathbf{x}}_{n+1|N}-{\mathbf{x}}_{n+1|n})$$

$${\mathbf{\Sigma }}_{n|N}={\mathbf{\Sigma }}_{n|n}+{\mathbf{B}}_n({\mathbf{\Sigma }}_{n+1|N}-{\mathbf{\Sigma }}_{n+1|n}){\mathbf{B}}_n^H$$

$${\mathbf{\Sigma }}_{n,n-1|N}={\mathbf{B}}_{n-1}{\mathbf{\Sigma }}_{n|N}$$

$$\displaystyle\leavevmode\resizebox{390.25534pt}{}{$\displaystyle\alpha^{[l+1]}=\frac{\sum_{n=2}^{N}\text{Tr}(\bm{\Sigma}_{n,n-1|N}{\mathbf{Q}^{[l]}}^{-1})+\mathbf{x}_{n-1|N}^{H}\mathbf{\mathbf{Q}^{[l]}}^{-1}\mathbf{x}_{n|N}}{\sum_{n=2}^{N}\text{Tr}(\bm{\Sigma}_{n-1|N}{\mathbf{Q}^{[l]}}^{-1})+\mathbf{x}_{n-1|N}^{H}{\mathbf{Q}^{[l]}}^{-1}\mathbf{x}_{n-1|N}}$},$$

$$\displaystyle\qquad\qquad\qquad\leavevmode\resizebox{303.53267pt}{}{$\displaystyle-\alpha^{[l+1]}(\mathbf{x}_{n-1|N}\mathbf{x}_{n|N}^{H}\!+\!\mathbf{x}_{n|N}\mathbf{x}_{n-1|N}^{H}\!+\!2\bm{\Sigma}_{n,n-1|N})\big{]}$}.$$

$$\displaystyle\qquad\qquad\quad\leavevmode\resizebox{294.85841pt}{}{$\displaystyle\!\!\!\!\!\!\!\!\!-\!\alpha^{[l+1]}(\mathbf{x}_{n-1|N}\mathbf{x}_{n|N}^{H}\!+\!\mathbf{x}_{n|N}\mathbf{x}_{n-1|N}^{H}\!+\!2\bm{\Sigma}_{n,n-1|N})\Big{]}$}.$$

$${\hat{D}}_n(f_j)=\frac{1}{K}\sum _{k=1}^K{\left|{\left({\hat{\mathbf{x}}}_n^{(k)}\right)}_j\right|}^2,$$

$$p(\varphi )=\frac{1}{2^{\nu }\Gamma (\nu )}\exp \left(\nu \varphi -\frac{1}{2}\exp (\varphi )\right).$$

$$\displaystyle+J\bigg{(}\nu\log 2+\log\Gamma(\nu)\bigg{)}+(\mathbf{s}_{n}\!-\theta\mathbf{s}_{n-1})^{H}\mathbf{R}_{n}^{-1}(\mathbf{s}_{n}\!-\theta\mathbf{s}_{n-1})$$

$$\displaystyle\qquad\qquad\qquad\qquad\quad\quad+\log|\det\mathbf{R}_{n}|\bigg{]}+{c},$$

Full text

I Introduction

Spectral analysis techniques are among the most important tools for extracting information from time series data recorded from naturally occurring processes. Examples include speech [2], images [3], electroencephalography (EEG) [4], oceanography [5], climatic time series [6] and seismic data [7]. Due to the exploratory nature of most of these applications, non-parametric techniques based on Fourier methods and wavelets are among the most widely used. In particular, the multitaper (MT) method excels among the available non-parametric techniques due to both its simplicity and control over the bias-variance trade-off via bandwidth adjustment [8, 9, 10].

Most existing spectral analysis techniques assume that the time series is stationary. In many applications of interest, however, the energy of the various oscillatory components in the data exhibits dynamic behavior. Extensions of stationary time series analysis to these non-stationary processes have led to ‘time-varying’ spectral descriptions such as the Wigner-Ville distribution [11, 12], the evolutionary spectra and its generalizations [13, 14], and the time-frequency operator symbol formulation [15] (See [16] for a detailed review). A popular approach to estimating such time-varying spectra is to subdivide the data into overlapping windows or segments and estimate the spectrum locally for each window using various Fourier or wavelet-based methods[17], assuming the underlying process is quasi-stationary, i.e., the spectrum changes slowly with time. Thereby, the so-called spectrogram analysis is obtained by using sliding windows with overlap in order to capture non-stationarity.

Although sliding window processing is widely used due to its fast implementation, it has several major drawbacks. First, the window length and extent of overlap are subjective choices and can drastically change the overall attribute of the spectrogram if chosen poorly. Second, given that the estimate associated to a given window is obtained by only the data within, it ignores the common dynamic trends shared across multiple windows, and thereby fails to fully capture the degree of smoothness inherent in the signal. Instead, the smoothness of the estimates is enforced by the amount of overlap between adjacent windows. Third, although techniques such as the MT analysis are able to mitigate the variabilities arising from finite data duration or the so-called ‘sampling’ noise by averaging over multiple tapers, their spectral resolution degrades when applied to data within small windows due to the increase in the Rayleigh resolution [18]. In addition, they do not have a mechanism in place to suppress the additive measurement noise that commonly contaminates empirical observations. Fourth, the overlap between adjacent windows hinders a precise statistical assessment of the estimates, such as constructing confidence intervals due to the high dependence of estimates across windows. To address this issue, statistical corrections for multiple comparisons need to be employed [19], which in turn limit the resulting test powers when multiple windows are involved.

In recent years, several alternative approaches to non-stationary spectral analysis have been proposed, such as the empirical mode decomposition (EMD) [20], synchrosqueezed wavelet transform [21, 22], time-frequency reassignment [23], time-frequency ARMA models [24], and spectrotemporal pursuit [25]. These techniques aim at decomposing the data into a small number of smooth oscillatory components in order to produce spectral representations that are smooth in time but sparse or structured in frequency. Although they produce spectral estimates that are highly localized in the time-frequency plane, they require certain assumptions on the data to hold. For example, EMD analysis assumes the signal to be deterministic and does not take into account the effect of observation noise [20]. Other methods assume that the underlying spectrotemporal components pertain to certain structures such as amplitude-modulated narrowband mixtures [21, 22], sparsity [25] or chirp-like dynamics [23]. In addition, they lack a statistical characterization of the estimates. Finally, although these sophisticated methods provide spectrotemporal resolution improvements, they do not yield implementations as simple as those of the sliding window-based spectral estimators.

In this paper, we address the above-mentioned shortcomings of sliding window multitaper estimators by resorting to state-space modeling. State-space models provide a flexible and natural framework for analyzing systems that evolve with time [26, 27, 28, 29], and have been previously used for parametric [30, 24] and non-parametric [25] spectral estimation. The novelty of our approach is in the integration of techniques from MT analysis and state-space modeling in a Bayesian estimation framework. To this end, we construct state-space models in which the underlying states pertain to the eigen-spectral quantities, such as the empirical eigen-coefficients and eigen-spectra arising in MT analysis. We employ state dynamics that capture the evolution of these quantities, coupled with observation models that reflect the effect of measurement and sampling noise. We then utilize Expectation-Maximization (EM) to find the maximum a posteriori (MAP) estimate of the states given the observed data to construct our spectral estimators as well as statistical confidence intervals.

We provide theoretical analysis of the bias-variance trade-off, which reveals two major features of our proposed framework: 1) our methodology inherits the control mechanism of the bias-variance trade-off from the MT framework by means of changing the design bandwidth parameters [10], and 2) our algorithms enjoy the optimal data combining and denoising features of Bayesian filtering and smoothing. In addition, due to the simplicity and wide usage of Bayesian filtering and smoothing algorithms, our algorithms are nearly as simple to implement as the sliding window-based spectrograms. To further demonstrate the performance of our algorithms, we apply them to synthetic as well as real data including human EEG recordings during sleep and electric network frequency data from audio recordings. Application of our proposed estimators to these data provides spectrotemporal features that are significantly denoised, are smooth in time, and enjoy high spectral resolution, thereby corroborating our theoretical results.

The rest of the paper is organized as follows: In Section II, we present the preliminaries and problem formulation. In Section III, we develop our proposed estimators. Application of our estimators to synthetic and real data are given in Section IV, followed by our theoretical analysis in Section V. Finally, our concluding remarks are presented in Section VI.

II Preliminaries and Problem Formulation

II-A Non-stationary Processes and Time-Varying Spectrum

Consider a finite realization of $T$ samples from a discrete-time non-stationary process $y_{t},t=1,2,\cdots T$, obtained via sampling a continuous-time signal above Nyquist rate. We assume that the non-stationary process $y_{t}$ is harmonizable so that it admits a Cramér representation[31] of the form:

[TABLE]

where $dz(f)$ is the generalized Fourier transform of the process. This process has a covariance function of the form:

[TABLE]

where $\gamma_{L}(f_{1},f_{2})\!\!:=\!\!\mathbbm{E}[dz(f_{1})dz^{*}(f_{2})]$ is referred to as the generalized spectral density or the Loève spectrum [32]. Due to the difficulty in extracting physically-plausible spectrotemporal information from the two-dimensional function $\gamma_{L}(f_{1},f_{2})$, other forms of spectrotemporal characterization that are two-dimensional functions over time and frequency have gained popularity [17]. To this end, by defining the coordinate rotations $t:=(t_{1}+t_{2})/2$, $\tau:=t_{1}-t_{2}$, $f:=(f_{1}+f_{2})/2$, and $g:=f_{1}-f_{2}$ and by substituting in the definition of the covariance function in (2), we obtain:

[TABLE]

where $f$ and $g$ are referred to as the ordinary and non-stationary frequencies, respectively[32], and $\gamma(g,f)dfdg:=\gamma_{L}(f_{1},f_{2})df_{1}df_{2}$ is the Loève spectrum in the rotated coordinates. To obtain one such two-dimensional spectral density representation over time and frequency, we define:

[TABLE]

which coincides with the expected value of the Wigner-Ville distribution [11]. The ‘time-varying’ spectral representation $D(t,f)$ captures the spectral information of the data as a function of time, and thus provides a useful framework for analyzing non-stationary time series. However, estimating $D(t,f)$ from finite samples of the process is challenging, considering that the expectation needs to be replaced by time averages which may smooth out the time-varying features of the signal [12].

In order to address this challenge, certain additional assumptions need to be imposed on the underlying process, which as a matter of fact restrict the extent of temporal or spectral variations the signal can exhibit. In this regard, two such popular assumptions are posed in terms of the quasi-stationary and underspread properties. In order to define these properties quantitatively, let us first define the Expected Ambiguity Function (EAF) [16, 14] as:

[TABLE]

A signal whose EAF has a limited spread along $g$ (i.e., negligible spectral correlation) is called quasi-stationary. If the signal is concentrated around the origin with respect to both $\tau$ and $g$ (i.e., small spectral and temporal correlation), it is called underspread. There exists a large body of work on estimating time-varying spectra under the quasi-stationarity assumption, such as short-time periodograms, pseudo-Wigner estimators [12], and estimates of the evolutionary spectra [13]. More recent methods such as the Generalized Evolutionary Spectra (GES) estimators [14] and time-frequency auto-regressive moving-average (TFARMA) [24] estimators rely on the underspread property of the underlying signals. We refer the reader to [33, Chapter 10] for a detailed discussion of these properties.

It is noteworthy that both assumptions are fairly general, encompass a broad range of naturally-occurring processes, and have resulted in successful applications in real life problems. We will next discuss one of the widely used methods for estimating time-varying spectra under the quasi-stationary assumption, which extends MT spectral analysis beyond second-order stationary processes and has gained popularity in exploratory studies of naturally occurring processes [32].

II-B The Sliding Window MT Spectral Analysis

One of the popular non-parametric techniques for estimating the ‘time-varying’ spectral representation $D(t,f)$ is achieved by subdividing the data into overlapping windows or segments and estimating the spectrum for each window independently using the MT method [32]. This method naturally and intuitively extends the popular non-parametric MT method to the non-stationary scenario under the assumption of quasi-stationarity, which enables one to treat the time series within segments (locally) as approximately second-order stationary [34]. The resulting spectrotemporal representations are smoothed version of the Wigner-Ville distribution [16] and are referred to as spectrogram. In what follows, we briefly describe the MT spectrogram method, since it preludes the rest of our treatment.

The MT method is an extension of single-taper spectral analysis, where the data is element-wise multiplied by a taper prior to forming the spectral representation to mitigate spectral leakage [9, 8]. In the MT method, spectral representation is computed as the average of several such single-taper PSDs, where the tapers are orthogonal to each other and exhibits good leakage properties. This can be achieved by using the discrete prolate spheroidal sequences (dpss) or Slepian sequences [35], due to their orthogonality and optimal leakage properties.

Another viewpoint of the MT method with this particular choice of data tapers is the decomposition of the spectral representation of the process over a set of orthogonal basis functions. Indeed, these basis functions originate from an approximate solution to the integral equation expressing the projection of $dz(f)$ onto the Fourier transform of the data:

[TABLE]

where $W$ is window length, i.e., number of samples, and $dz(\zeta)$ is an orthogonal increment process. This integral equation can be approximated using a local expansion of the increment process over an interval $[-B,B]$, for some small design band-width $B$, in the space spanned by the eigenfunctions of the Dirichlet kernel $\frac{\sin W\pi f}{\sin\pi f}$ [9, 8].

These eigenfunctions are known as the prolate spheroidal wave functions (PSWFs), which are a set of doubly-orthogonal functions over $[-B,B]$ and $[-{1}/{2},{1}/{2}]$, with time-domain representations given by the dpss sequences. Let $u_{l}^{(k)}$ be the $l$th sample of the $k$th dpss sequence, for a given bandwidth $B$ and window length $W$. The $k$th PSWF is then defined as:

[TABLE]

Choosing $K\leqslant\lfloor 2WB\rfloor-1$ dpss having eigenvalues close to $1$ as data tapers, the MT spectral estimate can be calculated as follows:

[TABLE]

where ${\widehat{x}^{(k)}(f)}:=\sum_{l=0}^{W-1}e^{-i2\pi fl}u_{l}^{(k)}y_{l}$ for $k=1,2,\cdots,K$ are called the ‘eigen-coefficients’. The ‘eigen-spectra’, $\widehat{S}^{(k)}(f):=|{\widehat{x}^{(k)}(f)}|^{2}$ can be viewed as the expansion coefficients of the decomposition.

To estimate time-varying spectra under the MT framework, sliding windows with overlap are used to enforce temporal smoothness and increase robustness of the estimates [36, 32], resulting in MT spectrogram estimates. Although this ‘overlapping’ MT procedure overcomes frequency leakage issues and produces consistent estimates, subjective choices of the window length and the degree of overlap can change the overall appearance of the spectrogram drastically when poor choices of these parameters are used. In addition, these estimates lack precise statistical inference procedures, such as hypothesis testing, due to the statistical dependence induced by the overlaps. The objective of this work is to overcome these limitations by directly modeling and estimating the evolution of the process without committing to overlapping sliding windows, while achieving fast and efficient implementations. Before presenting our proposed solutions, we give a brief overview of other existing approaches in the literature in order to put our contributions in context.

II-C Motivation and Connection to Existing Literature

Our goal is to overcome the foregoing challenges faced by the sliding window MT spectrogram analysis in estimating the expected value of the Wigner-Ville distribution (See Eq. (3)) under the quasi-stationarity assumption. In addition, our work can be viewed in the context of the spectrogram approximation to the evolutionary spectra [13]. The evolutionary spectra is obtained by considering an expansion of $y_{t}$ over the set of complex sinusoids as

[TABLE]

with uncorrelated time-varying expansion coefficients, i.e., $\displaystyle\mathbbm{E}[x_{t}(f_{1})x_{t}^{*}(f_{2})]=D(t,f_{1})\delta(f_{1}-f_{2})$. In standard MT spectrogram analysis, the extent of overlap between consecutive segments dictates the amount of temporal smoothness in the estimates. Our approach is to avoid the usage of overlapping windows by modeling and estimating the dependence of the spectra across windows using state-space models, while retaining the favorable leakage properties of the MT analysis. As will be revealed in the subsequent sections, in the same vein as the sliding window multitaper analysis our methods pertain to the class of spectrogram estimates, which are viewed as smoothed versions of the Wigner-Ville spectrum [16].

Due to the underlying quasi-stationarity assumption, i.e., negligible spectral correlation, the domain of applicability of our methods might be narrower than the more general non-stationary spectral analysis methods such as GES and Weyl spectral estimation and TFARMA modeling; however, our methods admit simple and efficient implementations, which makes them attractive for exploratory applications in which sliding window processing is widely used with subjective and ad hoc choices of design parameters. In this context, the novelty of our contributions lies in:

Capturing the evolution of the spectra across windows by modeling the dynamics of certain eigen-spectral quantities arising in MT analysis (e.g., spectral eigen-coefficients and eigen-spectra);

• 2)

Addressing the additive measurement noise and multiplicative sampling noise, which severely distort the spectrograms obtained by the multitaper framework; and

• 3)

Constructing a framework for precise statistical assessment of the estimates, by addressing the dependency among windows using a Bayesian formulation.

As it will be evident in Section V, the use of state-space models in the context of MT analysis results in adaptive weighting of the estimates of the eigen-coefficients or eigen-spectra across windows, thanks to the optimal data combining feature of Bayesian smoothing. These adaptive weights depend on the common dynamic trends shared across windows and hence result in capturing the degree of smoothness inherent in the signal, while producing estimates robust against uncertainties due to observation noise and limited data.

It is noteworthy that the use of state-space models in our work is significantly different from those used in parametric non-stationary spectral analysis methods such as TFARMA modeling [24]. In the TFARMA formulation, time delays and frequency shifts are used to model the non-stationary dynamics of the process in a physically intuitive way. These state-space models therefore determine the functional form of the resulting spectral estimates in closed form in terms of the finite set of ARMA coefficients. The state-space models used in our work, however, do not determine the functional form of the spectral estimates at each window, and rather control the temporal smoothness of the eigen-spectral quantities via forming a regularization mechanism in the underlying Bayesian estimation framework (See Section II-E). In our approach, we indeed estimate the spectrogram at a given number of frequency bins in each window, which scales with the total number of samples, in the same vein as sliding window MT spectrogram.

In light of the above, our algorithms belong to the class of semi-parametric estimation methods, as the underlying model is a hybrid of parametric unobservable state evolution process and a non-parametric data generating process [37]. In Section IV-A, we will compare our proposed semi-parametric methodology with both non-parametric and parametric techniques, namely the MT spectrogram analysis and the Time-Frequency Autoregressive (TFAR) modeling technique [24].

II-D Problem Formulation

Assume, without loss of generality, that an arbitrary window of length $W$ is chosen so that for some integer $N$, $NW=T$ and let \mathbf{y}_{n}=\big{[}y_{(n-1)W+1},y_{(n-1)W+2},\cdots,y_{nW}\big{]}^{\top} for $n=1,2\cdots N$, denotes the data in the $n$th window. This way, the entire data is divided into $N$ non-overlapping segments of length $W$ each. To this end, we invoke quasi-stationarity assumption by modeling $y_{t}$ to be stationary within each of these segments of length $W$. With this assumption, motivated by the major sources of uncertainty in spectral estimation, i.e., measurement noise and sampling noise, we formulate two state-space frameworks in the following subsections.

II-D1 Mitigating the Measurement Noise

Suppose that $\widetilde{y}_{t}$ is the noise corrupted observation obtained from the true signal $y_{t}$, i.e., $\widetilde{y}_{t}=y_{t}+v_{t}$, where $(v_{t})_{t=1}^{T}$ is an i.i.d. zero-mean Gaussian noise sequence with fixed variance $\sigma^{2}$. By discretizing the representation in (1) at a frequency spacing of $2\pi/J$ with $J$ an integer, at any arbitrary window $n$, we have

[TABLE]

where $\mathbf{F}_{n}$ is a matrix with elements (\mathbf{F}_{n})_{l,j}:=\exp\big{(}i2\pi(((n-1)W+l)\frac{j-1}{J})\big{)} for $l=1,2,\cdots,W$ and $j=1,2,,\cdots J$; \widetilde{\mathbf{y}}_{n}:=\big{[}\widetilde{y}_{(n-1)W+1},\widetilde{y}_{(n-1)W+2},\cdots,\widetilde{y}_{nW}\big{]}^{\top} is the noisy observation of the true signal $\mathbf{y}_{n}$; $x_{n}(f)$ and \mathbf{x}_{n}:=\big{[}x_{n}(0),x_{n}(2\pi\frac{1}{J}),\cdots,x_{n}(2\pi\frac{J-1}{J})\big{]}^{\top} denote the orthogonal increment process and its discretized version, respectively at window $n$ and \mathbf{v}_{n}=\big{[}v_{(n-1)W+1},v_{(n-1)W+2},\cdots,v_{nW}\big{]}^{\top} is zero-mean Gaussian noise with covariance $\text{Cov}\{\mathbf{v}_{i},\mathbf{v}_{j}\}=\sigma^{2}\mathbf{I}\delta_{i,j}$ .

Let \mathbf{u}^{(k)}\!:=\!\big{[}u_{1}^{(k)},u_{2}^{(k)},\cdots,u_{W}^{(k)}\big{]}^{\top} denotes the $k$th dpss taper and $\widetilde{\mathbf{y}}_{n}^{(k)}:=\mathbf{u}^{(k)}\odot\widetilde{\mathbf{y}}_{n}$, where $\odot$ denotes element-wise multiplication. Let $x_{n}^{(k)}(f)$ and $\mathbf{x}_{n}^{(k)}:=\left[x_{n}^{(k)}(0),x_{n}^{(k)}(2\pi\frac{1}{J}),\cdots,x_{n}^{(k)}(2\pi\frac{J-1}{J})\right]^{\top}$ denote the $k$th spectral eigen-coefficient of $\mathbf{y}_{n}$ and its discretized version, respectively, for $k=1,2,\cdots,K$. Then, following (7) we consider the following spectrotemporal representation of the tapered data segments:

[TABLE]

where $\mathbf{v}^{(k)}_{n}$ is the contribution of $\mathbf{v}_{n}$ to the $k$th tapered data, assumed to be independent of $\mathbf{x}_{1:n-1}^{(k)}$, and identically distributed according to a zero-mean Gaussian distribution with covariance $\text{Cov}\{\mathbf{v}^{(k)}_{i},\mathbf{v}^{(k)}_{j}\}={\sigma^{(k)}}^{2}\mathbf{I}\delta_{i,j}$. We view $\widetilde{\mathbf{y}}_{n}^{(k)}$ as a noisy observation corresponding to the true eigen-coefficient $\mathbf{x}_{n}^{(k)}$, which provides a linear Gaussian forward model for the observation process.

In order to capture the evolution of the spectrum and hence systematically enforce temporal smoothness, we impose a stochastic continuity constraint on the eigen-coefficients $(\mathbf{x}_{n}^{(k)})_{n=1}^{N}$ for $k=1,2,\cdots K$, using a first-order difference equation:

[TABLE]

where $0\leqslant\alpha^{(k)}<1$, and $\mathbf{w}^{(k)}_{n}$ is independent of $\mathbf{x}_{1:n-1}^{(k)}$ and assumed to be independently distributed according to a zero-mean Gaussian distribution with diagonal covariance $\text{Cov}\{\mathbf{w}^{(k)}_{i},\mathbf{w}^{(k)}_{j}\}=\mathbf{Q}^{(k)}_{i}\delta_{i,j}$. Under this assumption, the discrete-time process, $(\mathbf{x}_{n}^{(k)})_{n=1}^{N}$ forms a jointly Gaussian random process with independent increments, while the process itself is statistically dependent. An estimate of the unobserved states (true eigen-coefficients) from the observations (tapered data) under this model suppresses the measurement noise and captures the state dynamics.

II-D2 Mitigating the Sampling Noise

Suppose the additive measurement noise is negligible, i.e. $v_{t}\approxeq 0$. For now, consider only a single window of length $W$. It is known that when the spectrum does not rapidly vary over the chosen design bandwidth $B$, the eigen-spectra are approximately uncorrelated and the following approximation holds for $k=1,2\cdots,K$ [32, 9]:

[TABLE]

where $\widehat{S}^{(k)}(f)$ and $S(f)$ are the tapered estimate and true PSD, respectively. In other words, the empirical eigen-spectra of the process can be thought of as the true spectra corrupted by a multiplicative noise, due to sampling and having access to only a single realization of the process. We refer to this uncertainty induced by sampling as sampling noise. By defining $\psi^{(k)}(f):=\log{\widehat{S}^{(k)}(f)}+\log 2$ and $s^{(k)}(f):=\log{S(f)}$, we can transform the multiplicative effect of the sampling noise in Eq. (10) to the following additive forward model[38]:

[TABLE]

where $\phi^{(k)}(f)$ is a log-chi-square distributed random variable, capturing the uncertainty due to sampling noise. It can be shown that $\phi^{(k)}(f)$ has a density given by:

[TABLE]

which belongs to the family of log-Gamma distributions, including the Gumbel and Bramwell-Holdsworth-Pinton distributions common in extreme value statistics [39].

In order to incorporate this observation model in our dynamic framework, we define the state vector, \mathbf{s}_{n}^{(k)}\!:=\!\Big{[}s_{n}^{(k)}(0),s_{n}^{(k)}\Big{(}\textstyle 2\pi\frac{1}{J}\Big{)},\cdots,s_{n}^{(k)}\Big{(}\textstyle 2\pi\frac{J-1}{J}\Big{)}\Big{]}^{\top}; the observation vector, \bm{\psi}_{n}^{(k)}\!:=\!\Big{[}\psi_{n}^{(k)}(0),\psi_{n}^{(k)}\Big{(}\textstyle 2\pi\frac{1}{J}\Big{)},\cdots,\psi_{n}^{(k)}\Big{(}\textstyle 2\pi\frac{J-1}{J}\Big{)}\Big{]}^{\top} and the observation noise vector, \bm{\phi}_{n}^{(k)}:=\Big{[}\phi_{n}^{(k)}(0),\phi_{n}^{(k)}\Big{(}\textstyle 2\pi\frac{1}{J}\Big{)},\cdots,\phi_{n}^{(k)}\Big{(}\textstyle 2\pi\frac{J-1}{J}\Big{)}\Big{]}^{\top}. Then, the forward model at window $n$ can be stated as:

[TABLE]

where each element of $\bm{\phi}_{n}^{(k)}$ is log-chi-square distributed. Similar to the preceding model, we impose a stochastic continuity constraint over the logarithm of the eigen-spectra as follows:

[TABLE]

where $0\leqslant\theta^{(k)}<1$, and $\mathbf{e}_{n}^{(k)}$ is assumed to be a zero-mean Gaussian vector independent of $\mathbf{s}_{1:n-1}^{(k)}$ and with a diagonal covariance $\text{Cov}\{\mathbf{e}^{(k)}_{i},\mathbf{e}^{(k)}_{j}\}=\mathbf{R}^{(k)}_{i}\delta_{i,j}$. Note that the logarithm function maps the range of the eigen-spectra in $[0,\infty)$ to $(-\infty,\infty)$ which makes the Gaussian state evolution plausible. An estimate of the unobserved states (logarithm of the true spectra) from the observations (logarithm of the empirical eigen-spectra) under this model suppresses the sampling noise and captures the state dynamics.

In summary, through these models we project the data of each short window onto the functional space spanned by the PSWFs and impose stochastic continuity constraints ((9) and (14)) on these projections (eigen-coefficients or eigen-spectra) in order to recover spectral representations that are smooth in time and robust against measurement or sampling noise. Fig. 1 provides a visual illustration of the proposed modeling paradigm.

II-E The Inverse Problem

We formulate the spectral estimation problem as one of Bayesian estimation, in which the Bayesian risk/loss function, fully determined by the posterior density of $(\mathbf{x}_{n}^{(k)})_{n=1,k=1}^{N,K}$ (resp. $(\mathbf{s}_{n}^{(k)})_{n=1,k=1}^{N,K}$) given the observations $(\widetilde{\mathbf{y}}_{n}^{(k)})_{n=1,k=1}^{N,K}$ (resp. $(\bm{\psi}_{n}^{(k)})_{n=1,k=1}^{N,K}$) is minimized. We first consider the forward model of (8), which provides the observed data likelihood given the states. Under the state-space model of Eq. (9), the $k$th eigen-coefficient can be estimated by solving the following maximum a posteriori (MAP) problem:

[TABLE]

for $k=1,2,\cdots,K$. Similarly, in the second state space framework (13), the eigen-spectra can be obtained by solving another MAP estimation problem:

[TABLE]

for $k=1,2,\cdots,K$, where $\mathbf{1}_{J}$ is the vector of all ones of length $J$. We call the MAP estimation problems in (II-E) and (II-E) the Dynamic Bayesian Multitaper (DBMT) and the log-DBMT estimation problems, respectively. Similarly, the respective spectrogram estimates will be denoted by the DBMT and log-DBMT estimates.

Equation (II-E) is a strictly convex function of $\mathbf{x}^{(k)}_{n}\in\mathbb{C}^{W}$ and $\mathbf{Q}^{(k)}_{n}\in\mathbb{S}^{W}_{++}$ for $n=1,2\cdots,N$, which can be solved using standard optimization techniques. However, these techniques do not scale well with the data length $N$. A careful examination of the log-posterior reveals a block tri-diagonal structure of the Hessian, which can be used to develop efficient recursive solutions that exploit the temporal structure of the problem. A similar argument holds for the optimization problem in (II-E). However, the parameters of these state-space models need to be estimated from the data. In the next section, we show how the EM algorithm can be used to both estimate the parameters and states efficiently from the optimization problems (II-E) and (II-E).

III Fast Recursive Solutions via the EM Algorithm

In order to solve the MAP problem in (II-E), we need to find the parameters $\mathbf{Q}_{n}^{(k)}\in\mathbb{S}^{W}_{++}$ and $\alpha^{(k)}\in(0,1]$ for $n=1,2\cdots,N$ and $k=1,2,\cdots,K$. Similarly $\mathbf{R}_{n}^{(k)}\in\mathbb{S}^{W}_{++}$ and $\theta^{(k)}\in(0,1]$ need to be estimated for the problem in (II-E). If the underlying states were known, one could further maximize the log-posterior with respect to the parameters. This observation can be formalized in the EM framework [40, 28, 27]. To avoid notational complexity, we drop the dependence of the various variables on the taper index $k$ in the rest of this subsection.

III-A The DBMT Spectrum Estimation Algorithm

By treating $\mathbf{x}_{n},n=1,2,\cdots,N$ as the hidden variables and $\alpha,\mathbf{Q}_{n},n=1,2,\cdots,N$ as the unknown parameters to be estimated, we can write the complete log-likelihood as:

[TABLE]

where $c$ represents the terms that do not depend on $\alpha$, $(\mathbf{Q}_{n})_{n=1}^{N}$ or $(\mathbf{x}_{n})_{n=1}^{N}$. For simplicity of exposition, we assume that $\mathbf{Q}_{n}=\mathbf{Q}$ for $n=1,2,\cdots,N$. The forthcoming treatment can be extended to the general case with little modification. Also, note that $\sigma^{2}$ can be absorbed in $\mathbf{Q}$, and thus is assumed to be known. At the $l$th iteration, we have:

III-A1 E-Step

Given $\alpha^{[l]},\mathbf{Q}^{[l]}$, for $n=1,2,\cdots,N$, the expectations, $\mathbf{x}_{n|N}:=\mathbbm{E}[\mathbf{x}_{n}|\widetilde{\mathbf{y}}_{1:N},\alpha^{[l]},\mathbf{Q}^{[l]}],$ $\bm{\Sigma}_{n|N}:=\mathbbm{E}[(\mathbf{x}_{n}-\mathbf{x}_{n|N})(\mathbf{x}_{n}-\mathbf{x}_{n|N})^{H}|\widetilde{\mathbf{y}}_{1:N},\alpha^{[l]},\mathbf{Q}^{[l]}],$ $\bm{\Sigma}_{n,n-1|N}:=\mathbbm{E}[(\mathbf{x}_{n}\!-\!\mathbf{x}_{n|N})(\mathbf{x}_{n-1}\!-\!\mathbf{x}_{n-1|N})^{H}|\widetilde{\mathbf{y}}_{1:N},\alpha^{[l]},\mathbf{Q}^{[l]}],$ can be calculated using the Fixed Interval Smoother (FIS) [41] (lines 4 and 5) and the state-space covariance smoothing algorithm [42] (line 6). These expectations can be used to compute the expectation of the complete data log-likelihood \mathbbm{E}\big{[}\log L(\alpha,\mathbf{Q})|\widetilde{\mathbf{y}}_{1:N},\alpha^{[l]},\mathbf{Q}^{[l]}\big{]}.

III-A2 M-Step

The parameters for subsequent iterations, $\alpha^{[l+1]}$ and $\mathbf{Q}^{[l+1]}$ can be obtained by maximizing the expectation of (III-A). Although this expectation is convex in $\alpha$ and $\mathbf{Q}$ individually, it is not a convex function of both. Hence, we perform cyclic iterative updates for $\alpha^{[l+1]}$ and $\mathbf{Q}^{[l+1]}$ given by:

[TABLE]

and

[TABLE]

These iterations can be performed until convergence to a possibly local maximum. However, with even one such update, the overall algorithm forms a majorization-minimization (MM) procedure, generalizing the EM procedure and enjoying from similar convergence properties [43]. One possible implementation of this iterative procedure is described in Algorithm 1. Once the DBMT estimates of all the $K$ eigen-coefficients $\widehat{\mathbf{x}}_{n}^{(k)}$ are obtained, for $n=1,2,\cdots,N$ and $k=1,2,\cdots,K$, the DBMT spectrum estimate is constructed similar to (5):

[TABLE]

where $f_{j}:=\frac{2\pi(j-1)}{J}$ for $j=1,2,\cdots,J$ and $n=1,2,\cdots,N$. Confidence intervals can be computed by mapping the Gaussian confidence intervals for $\widehat{\mathbf{x}}_{n}^{(k)}$’s to the final DBMT estimate.

III-B The log-DBMT Spectrum Estimation Algorithm

We utilize a similar iterative procedure based on the EM algorithm to find the log-DBMT spectrum estimate. As before, we treat $\mathbf{s}_{n},n=1,2,\cdots,N$ as hidden variables and $\theta,\mathbf{R}_{n},n=1,2,\cdots,N$ as the unknown parameters to be estimated. In order to give more flexibility to the observation model, we consider the observation noise to be distributed as log-chi-square with degrees of freedom $2\nu$, for some positive integer $\nu$ to be estimated. The density of each element of $\bm{\phi}_{n}^{(k)}$ is then given by:

[TABLE]

We can express the complete data log-likelihood as:

[TABLE]

where $c$ represents the terms that do not depend on $\nu$, $\theta$, $(\mathbf{R}_{n})_{n=1}^{N}$ or $(\mathbf{s}_{n})_{n=1}^{N}$. Again assuming $\mathbf{R}_{n}=\mathbf{R}$ for all $n=1,2,\cdots,N$ for simplicity, the following EM algorithm can be constructed:

III-B1 E-step

Computation of the conditional expectation of the log-likelihood in (III-B) requires evaluating $\mathbbm{E}[\mathbf{s}_{n}|\bm{\psi}_{1:N},\mathbf{R}^{[l]},\theta^{[l]},\nu^{[l]}]$ and $\mathbbm{E}[\exp(-\mathbf{s}_{n})|\bm{\psi}_{1:N},\mathbf{R}^{[l]},\theta^{[l]},\nu^{[l]}]$ for $n=1,2,\cdots,N$. Unlike the DBMT estimation problem, the forward model in this case is non-Gaussian, and hence we cannot apply the Kalman filter and FIS to find the state expectations. To compute the conditional expectation, the distribution of $\mathbf{s}_{n}|\bm{\psi}_{1:n},\mathbf{R}^{[l]},\theta^{[l]},\nu^{[l]}$ or its samples are required [29]. Computation of the distribution $\mathbf{s}_{n}|\bm{\psi}_{1:n},\mathbf{R}^{[l]},\theta^{[l]},\nu^{[l]}$ involves intractable integrals and sampling from the distribution using numerical methods such as Metropolis-Hastings is not computationally efficient, especially for long data, given that it has to be carried out at every iteration.

Since the posterior distribution is unimodal and a deviation from the Gaussian posterior, we approximate the distribution of $\mathbf{s}_{n}|\bm{\psi}_{1:n},\mathbf{R}^{[l]},\theta^{[l]},\nu^{[l]}$ as a Gaussian distribution by matching its mean and covariance matrix to the log-posterior in Eq. (III-B). To this end, the mean is approximated by the mode of $f_{\mathbf{s}_{n}|\bm{\psi}_{1:n},\mathbf{R}^{[l]},\theta^{[l]},\nu^{[l]}}$ and the covariance is set to the inverse of the negative Hessian of the log-likelihood in (III-B) [44, 28]. Under this approximation, computing $\mathbbm{E}[\exp(-\mathbf{s}_{n})|\bm{\psi}_{1:n},\mathbf{R}^{[l]},\theta^{[l]},\nu^{[l]}]$ is also facilitated thanks to the closed-form moment generating function of $\mathbf{z}\sim\mathcal{N}(\bm{\mu},\bm{\Sigma})$:

[TABLE]

Similar to the case of DBMT, we can exploit the block tri-diagonal structure of the Hessian in (II-E) to carry out the E-step efficiently using forward filtering and backward smoothing.

III-B2 M-step

Once the conditional expectation of the log-likelihood in (III-B) given $\bm{\psi}_{1:n}$, $\mathbf{R}^{[l]}$, $\theta^{[l]}$, $\nu^{[l]}$ is available, we can update $\mathbf{R}^{[l+1]}$ and $\theta^{[l+1]}$ using similar closed form equations as in (III-A2). But updating $\nu^{[l+1]}$ by maximizing the conditional expectation of the log likelihood in (III-B) wrt. $\nu^{[l+1]}$ requires solving following nonlinear equation:

[TABLE]

where $\digamma(\cdot)$ is the digamma function. We can use Newton’s method to solve this equation up to a given precision. An implementation of the log-DBMT is given by Algorithm 2. Note that unlike the DBMT algorithm which pertains to a Gaussian observation model, the forward filtering step to compute $\mathbf{s}_{n|n}$ is nonlinear, and standard techniques such as Newton’s method can be used to solve for $\mathbf{s}_{n|n}$. We use the log-DBMT algorithm to find all the $K$ estimates of true log-spectra and construct the log-DBMT estimate as:

[TABLE]

where $f_{j}:=\frac{2\pi(j-1)}{J}$ for $j=1,2,\cdots,J$ and $n=1,2,\cdots,N$. Again, confidence intervals can be computed by mapping the Gaussian confidence intervals for $\widehat{\mathbf{s}}_{n}^{(k)}$’s to the final log-DBMT estimate.

III-C Parameter Selection

The window length $W$, design bandwidth $B$, and the number of tapers $K$ need to be carefully chosen. Since both proposed algorithms are motivated by the standard overlapping MT method, we use the same guidelines for choosing these parameters [32, 36]. The window length $W$ is determined based on the expected rate of change of the PSD (given domain-specific knowledge) in order to make sure that the quasi-stationarity assumption holds. The design bandwidth $B$ is chosen small enough to be able to resolve the dominant frequency components in the data, while being large enough to keep the time-bandwidth product $\rho:=WB\geqslant 1$. The number of tapers $K$ is then chosen as $K\leqslant\lfloor 2\rho\rfloor-1$ [32].

IV Application to Synthetic and Real Data

Before presenting our theoretical analysis, we examine the performance of DBMT and log-DBMT spectrogram estimators on synthetic data, and then demonstrate their utility in two real world data applications, namely spectral analysis of human EEG during sleep and Electric Network Frequency signal detection.

IV-A Application to Synthetic Data

The synthetic data consists of the linear combination of two amplitude-modulated and frequency-modulated processes with high dynamic range (i.e., high-Q). The amplitude-modulated component $y^{(1)}_{t}$ is generated through modulating an $\text{AR}(6)$ process tuned around $11~{}\text{Hz}$ by a cosine at a low frequency $f_{0}=0.02~{}\text{Hz}$. The frequency-modulated component $y^{(2)}_{t}$ is a realization of an $\text{ARMA}(6,4)$ with varying pole loci. To this end, the process has a pair of 3rd order poles at $\omega_{t}:=2\pi f_{t}$ and $-\omega_{t}$, where $f_{t}$ increases from $5~{}\text{Hz}$, starting at $t=0$, every $\sim 26~{}\text{s}$ by increments of $0.48~{}\text{Hz}$, to achieve frequency modulation. In summary, the noisy observations are given by:

[TABLE]

where $v_{t}$ is a white Gaussian noise process and $\sigma$ is chosen to achieve an SNR of $30~{}\text{dB}$. The process is truncated at $600~{}\text{s}$ to be used for spectrogram analysis. Figure 2 shows a $12$ second sample window of the process.

In addition to the standard overlapping MT, we present comparison to the TFAR method as an example of parametric state-space modeling approaches to non-stationary spectral analysis. This method is known to be well suited to processes whose time-varying spectra exhibits sharp peaks, i.e., signals consisting of several narrow-band components [24]. The TFAR model is defined by the input-output relation

[TABLE]

where $e_{t}$ is a stationary white noise process with unit variance, and $(a_{m,l})_{m=1,l=-L_{A}}^{M_{A},L_{A}}$ and $(b_{0,l})_{l=-L_{B}}^{L_{B}}$ are the autoregressive (AR) and zero-delay moving average (MA) parameters, respectively. The integers $M_{A}$ and $L_{A}$ are respectively the delay and Doppler model orders of the AR component and $L_{B}$ denotes the Doppler model order of the zero-delay MA component. The AR and MA parameters are estimated by solving the time-frequency Yule-Walker equations, from which the evolutionary spectra can be constructed (See the methods described in [24] for more details).

Figure 3 shows the true as well as estimated spectrograms by the standard overlapping MT, DBMT, log-DBMT and the TFAR estimators. Each row consists of three panels: the left panel shows the entire spectrogram; the middle panel shows a zoomed-in spectrotemporal region marked by the dashed box in the left panel; and the right panel shows the PSD along with confidence interval (CI) in gray hull, at a selected time point marked by a dashed vertical line in the middle panel. Note that for the standard MT estimates, the CIs are constructed assuming a $\chi_{2K}^{2}$ distribution of the estimates around the true values [32, 9], whereas for DBMT and log-DBMT estimate by mapping the Gaussian confidence intervals for eigen-coefficients or eigen-spectra to the final estimates. We were not able to evaluate the CIs for the TFAR estimates, since to the best of our knowledge we are not aware of any method to do so. Figure 3A shows true spectrogram of the synthetic process, in which the existence of both amplitude and frequency modulations makes the spectrogram estimation a challenging problem.

Fig. 3B shows the standard overlapping MT spectrogram estimate. We used windows of length $6~{}\text{s}$ and the first $3$ tapers corresponding to a time-bandwidth product of $3$ and $50\%$ overlap to compute the estimates (note that the same window length, tapers and time-bandwidth product are used for the DBMT and log-DBMT estimators). Although the standard MT spectrogram captures the dynamic evolution of both components, it is blurred by the background noise and picks up spectral artifacts (i.e., vertical lines) due to window overlap, frequency mixing, and sampling noise. Fig. 3C demonstrates how the DBMT spectrogram estimate overcomes these deficiencies of the overlapping MT spectrogram: the spectrotemporal localization is sharper and smoother across time, artifacts due to overlapping between windows are vanished, and frequency mixing is further mitigated. By comparing the right panel of the second and third rows, two important observations can be made: first, the DBMT captures the true dynamic range of the original noiseless PSD, while the standard MT estimate fails to do so. Second, the CIs in Fig. 3C as compared to 3B are wider when the signal is weak (e.g., near $5~{}\text{Hz}$) and tighter when the signal is strong (e.g., near $11~{}\text{Hz}$). The latter observation highlights the importance of the model-based confidence intervals in interpreting the denoised estimates of DBMT: while the most likely estimate (i.e., the mean) captures the true dynamic range of the noiseless PSD, the estimator does not preclude cases in which the noise floor of $-40~{}\text{dB}$ is part of the true signal, while showing high confidence in detecting the spectral content of the true signal that abides by the modeled dynamics.

Next, Fig. 3D shows the log-DBMT spectrogram estimate, which shares the artifact rejection feature of the DBMT spectrogram. However, the log-DBMT estimate is smoother than both the standard overlapping MT and DBMT spectrograms in time as well as in frequency (see the zoomed-in middle panels), due to its sampling noise mitigation feature (by design). In addition, similar to the standard MT estimate, the log-DBMT estimator is not able to denoise the estimate by removing the observation noise. Though, the confidence intervals of the log-DBMT PSD estimate are tighter than those of the standard overlapping MT estimate due to averaging across multiple windows via Bayesian filtering/smoothing. As we will show in Section V, these qualitative observations can be established by our theoretical analysis.

Finally, Fig. 3E shows the TFAR spectral estimate. The model orders are chosen as $M_{A}=20$, $L_{A}=25$, $L_{B}=25$, large enough to allow the parametric model to achieve high time-frequency resolution. As it is shown in the right panel, the TFAR method provides the smoothest estimate along the frequency axis. However, it is not successful in capturing the true dynamic range of the signal due to spectral leakage (See middle and right panels). In addition, it is contaminated by similar vertical frequency artifacts as in the case of the MT spectrogram (See Figure 3A).

In the spirit of easing reproducibility, we have deposited a MATLAB implementation of these algorithms on the open source repository GitHub [45], which generates Figure 3.

IV-B Application to EEG data

To illustrate the utility of our proposed spectrogram estimators, we apply them to human EEG data recorded during sleep. In the interest of space, in the remainder of this section, we only present comparisons with the MT spectrogram as a non-parametric benchmark. The EEG data set is available online as part of the SHHS Polysomnography Database (https://www.physionet.org/pn3/shhpsgdb/). The data is $900~{}\text{s}$ long during stage 2 sleep, and sampled at $250~{}\text{Hz}$. During stage 2 sleep, the EEG is known to manifest delta waves ($0-4~{}\text{Hz}$) and sleep spindles (transient wave packets with frequency $12-14~{}\text{Hz}$) [46, 47]. Accurate localization of these spectrotemporal features has significant applications in studying sleep disorders and cognitive function [46]. Since the transient spindles occur at a time scale of seconds, we choose a window length of $2.25~{}\text{s}$ for all algorithms (with $50\%$ overlap for the standard overlapping MT estimate). We also chose a time-bandwidth product of $2.25$ for all algorithms, in order to keep the frequency resolution at $2~{}\text{Hz}$. Figs. 4A, B and C show the MT, DBMT and log-DBMT spectrogram estimates, respectively, with a similar presentational structure as in Fig. 3. As the middle panels reveal, the overlapping MT estimate is not able to clearly distinguish the delta waves and sleep spindles due to high background noise. The DBMT estimate shown in Fig. 4B, however, provides a significantly denoised spectrogram, in which the delta waves and sleep spindles are visually separable. The log-DBMT estimator shown in Fig. 4C provides significant spectrotemporal smoothing, and despite not fully reducing the background noise, provides a clear separation of the delta waves and spindles (see the PSD in the right panel). Similar to the analysis of synthetic data, the same observations regarding the confidence intervals of the estimators can be made.

IV-C Application to ENF data

Finally, we examine the performance of our proposed algorithms in tracking the Electrical Network Frequency (ENF) signals from audio recordings. The ENF signal corresponds to the supply frequency of the power distribution network which is embedded in audio recordings [48, 49]. The instantaneous values of this time-varying frequency and its harmonics form the ENF signal. The ability to detect and track the spectrotemporal dynamics of ENF signals embedded in audio recordings has shown to be crucial in data forensics applications [48].

Fig. 5A shows the spectrogram estimates around the sixth harmonic of the nominal $60~{}\text{Hz}$ ENF signal (data from [49]). We used $1000~{}\text{s}$ of audio recordings, and constructed spectrograms with windows of length $5~{}\text{s}$ and using the first $3$ tapers corresponding to a time-bandwidth product of $3$ for all three methods (with $25\%$ overlap for the overlapping MT estimate). The two dominant components around the sixth ENF harmonic exhibit temporal dynamics, but are hard to distinguish from the noisy background. Fig. 5B shows the DBMT spectrogram, in which the background noise is significantly suppressed, yielding a crisp and temporally smooth estimate of the ENF dynamics. The log-DBMT estimate is shown in Fig. 5C, which provides higher spectrotemporal smoothness than the standard MT estimate. Although the log-DBMT shows smaller variability in the estimates (middle and right panels), the gain is not as striking as in the cases of synthetic data and EEG analysis, due to the usage of longer windows which mitigates the sampling noise for all algorithms. Similar observations as in the previous two cases regarding the statistical confidence intervals can be made, which highlight the advantage of modeling the spectrotemporal dynamics in spectrogram estimation.

V Theoretical Analysis

V-A Filter Bank Interpretation

In order to characterize the spectral properties of any non-parametric spectrum estimator, the tapers applied to the data need to be carefully inspected. In the MT framework, the dpss sequences are used as tapers, which are known to produce negligible side-lobes in the frequency domain [8, 9]. The DBMT and log-DBMT algorithms also use the dpss tapers to alleviate the problem of frequency leakage. However, because of the stochastic continuity constraint we introduced, the estimate associated to any given window is now a function of the data in all the windows. Therefore, the theoretical properties of the MT method do not readily apply to our estimators.

To characterize the statistical properties of our estimates, we first need take a detour from the usual analysis of spectrum estimation techniques. In what follows, we mainly focus on the DBMT algorithm for the sake of presentation. By virtue of the FIS procedure under the assumptions that: 1) the window length $W$ is an integer multiple of $J$, the number of discrete frequencies, so that $\mathbf{F}_{n}=\mathbf{F}_{1},\forall n$, and 2) the state noise covariance matrices are time-invariant, i.e., $\mathbf{Q}_{n}=\mathbf{Q},\forall n$, one obtains the following expansion of $\mathbf{x}_{n|N}^{(k)}$ in terms of the observed data [25]:

[TABLE]

In other words, the DBMT algorithm maps the entire data $\widetilde{\mathbf{y}}:=[\widetilde{y}_{1},\widetilde{y}_{2},\cdots,\widetilde{y}_{T}]^{\top}$ to the vector of coefficients $\widehat{\mathbf{X}}^{(k)}$ according to [25]:

[TABLE]

where $\mathbf{F}$ and $\mathbf{U}^{(k)}$ are block-diagonal matrices with $\mathbf{F}_{1}$ and $\mathbf{U}_{k}:=\text{diag}[\mathbf{u}^{(k)}]$ as the diagonal blocks, respectively, and $\mathbf{G}$ is a weighting matrix which depends only on $\mathbf{Q}_{\infty}=\lim_{l\rightarrow\infty}\mathbf{Q}^{[l]}$, $\alpha_{\infty}=\lim_{l\rightarrow\infty}\alpha^{[l]}$, and window length, $W$. The rows of $\mathbf{G}^{(k)}\mathbf{F}^{H}\mathbf{U}^{(k)}$ form a filter bank whose output is equivalent to the time-frequency representation.

In order to continue our analysis, we make two assumptions common in the analysis of adaptive filters[50, 51]. First, we assume that the parameter estimates $\mathbf{Q}_{\infty}$ and $\alpha_{\infty}$ are close enough to the true values of $\mathbf{Q}$ and $\alpha$, and therefore replace them by $\mathbf{Q}$ and $\alpha$, i.e., as if the true parameters were known. Note that we have discarded the dependence of $\mathbf{Q}$ and $\alpha$ on $k$ for lucidity of analysis. Second, noting that $\alpha(\mathbf{I}-\mathbf{K}_{m}\mathbf{F}_{m})=\alpha\bm{\Sigma}_{m|m}\bm{\Sigma}_{m|m-1}^{-1}$ and $\mathbf{B}_{m}=\alpha\mathbf{\Sigma}_{m|m}\mathbf{\Sigma}_{m+1|m}^{-1}$ and that in steady state we have $\bm{\Sigma}_{m|m}:=\bm{\Sigma}_{\infty}$ and $\bm{\Sigma}_{m|m-1}=\alpha^{2}\bm{\Sigma}_{\infty}+\mathbf{Q}$, Eq. (V-A) can be approximated by:

[TABLE]

for $1\ll n\ll N$, where $\mathbf{\Lambda}=\alpha\mathbf{\Sigma}_{\infty}(\alpha^{2}\mathbf{\Sigma}_{\infty}+\mathbf{Q})^{-1}$, and \mathbf{\Gamma}=(\alpha^{2}\mathbf{\Sigma}_{\infty}+\mathbf{Q})\big{[}\mathbf{I}-rW\big{(}(\alpha^{2}\mathbf{\Sigma}_{\infty}+\mathbf{Q})^{-1}+rW\mathbf{I}\big{)}^{-1}\big{]}. That is, for values of $n$ far from the data boundaries, the weighting matrix is equivalent to a weighted set of dpss tapers in matrix form acting on all the data windows with an exponential decay with respect to the $n$th window.

As an example, the equivalents filters of the DBMT estimator corresponding to the first taper for the $11$Hz and $9$Hz frequencies around $300$ secs, from the synthetic data example are shown in Fig. 6. They are also compared to the equivalent filters corresponding to first taper of standard MT method in the frequency domain. As apparent from Fig. 6, the weighting matrix sets the gain of these filters in an adaptive fashion across all windows, unlike the standard MT method which only uses the data in window $n$. In addition, the filter corresponding to frequency of $9$Hz, which is negligible in the data, is highly attenuated, resulting in significant noise suppression. In this sense, the proposed estimation method can be identified as a data-driven denoising method for constructing time-frequency representations given noisy time series data. Next, we will characterize the performance of the DBMT estimator in terms of bias-variance trade-off.

V-B Bias and Variance Analysis

We first consider the implication of the stochastic continuity constraint of Eq. (9) on the evolution of the orthogonal increment processes governing the time series data. We first assume that the parameters $\alpha^{(k)}=\alpha$, for all $k=1,2,\cdots,K$. Suppose that the data in window $n$ has a Cramér representation with an orthogonal increment process $dz_{n}(f)$, $n=1,2,\cdots,N$. Then, one way to achieve the stochastic continuity of Eq. (9) is to assume:

[TABLE]

where $d\epsilon_{n}(f)$ is a Gaussian orthogonal increment process, independent of $dz_{n}(f)$. In the forthcoming analysis we also assume the locally stationarity condition, i.e., the generalized Fourier transform of the process remains stationary within each window. This assumption is common in the analysis of non-parametric spectral estimators [36, 32]. Finally, we assume a scaling of $K,N,W\rightarrow\infty$, $B\rightarrow 0$, $BW\rightarrow\rho$, for some constant $\rho$ [52]. The following theorems characterize the bias and variance of the DBMT estimator:

Theorem 1**.**

Suppose that the locally stationary process $y_{t}$ is governed by orthogonal increment processes evolving according to the dynamics $dz_{n+1}(f)=\alpha dz_{n}(f)+d\epsilon_{n}(f)$, with $\alpha<1$, where the noise process $d\epsilon_{n}(f),\forall f\in(-1/2,1/2]$ is a zero-mean Gaussian increment process with variance $q(f)>0$, independent of $dz_{n}(f)$. If the process is corrupted by additive zero-mean Gaussian noise with variance $\sigma^{2}$, then for $f\in\{f_{1},f_{2},\cdots,f_{J}\}$, the DBMT estimate satisfies:

[TABLE]

where $\lambda_{k}$ is the eigenvalue associated with the $k$th PSWF, $D(f):=q(f)/(1-\alpha)$, and $\kappa_{n}(f),\mu_{n}(f)$ are functions of $\alpha$ and $q(f)$ and explicitly given in the proof.

Theorem 2**.**

Under the assumptions of Theorem 1, the variance of the DBMT estimate $\widehat{D}_{n}(f)$ satisfies:

[TABLE]

The proofs of Theorems 1 and 2 integrate the treatment of [52] with the structure of the FIS estimates, and are presented in Appendix A. In order to illustrate the implications of these theorems, several remarks are in order:

Remark 1. The function $\kappa_{n}(f)$ controls the trade-off between bias and variance: for values of $\kappa_{n}(f)<1$, the bound on the variance decreases while the bias bound increases, and for $\kappa_{n}(f)\approx 1$, all the terms in the bias bound become negligible, while the variance bound increases. The function $\mu_{n}(f)$, on the other hand, reflects observation noise suppression in both the bias and variance. Note that these upper bounds are tight and achieved for a signal with flat spectrum.

Remark 2. The bias and variance bounds of [52] for the standard MT method can be recovered by setting $\kappa_{n}(f)=1$, $\mu_{n}(f)=1$, and $\sigma^{2}=0$ in the results of Theorems 1 and 2, i.e., in the absence of signal dynamics and measurement noise. For the DBMT estimator, signal and measurement noise variances, respectively contribute to the bias/variance upper bounds in different fashions through $\kappa_{n}(f)$ and $\mu_{n}(f)$, due to the distinction of the signal and measurement noise in our state-space model. In contrast, in the standard MT method, possible measurement noise is treated in the same way as the true data, and hence both the signal and noise variances have equal contributions in the estimator bias/variance.

The functions $\kappa_{n}(f)$ and $\mu_{n}(f)$ do not have closed-form expressions with respect to the state-space parameters $\alpha$, $\sigma^{2}$ and $q(f)$. In order to illustrate the roles of $\kappa_{n}(f)$ and $\mu_{n}(f)$, we consider the scenario under which the upper bounds on the bias and variance are achieved, i.e., $q(f)$ being independent of $f$, and hence $\kappa_{n}(f)=\kappa_{n}$ and $\mu_{n}(f)=\mu_{n}$, $\forall f$. In this scenario, even though the dependence of $\mu_{n}$ and $\kappa_{n}$ on the state-space parameters are quite involved, it is possible to obtain upper and lower bounds on $\kappa_{n}$ and $\mu_{n}$. As it is shown in Proposition 3 in Appendix B, the main parameters determining the behavior of $\kappa_{n}$ and $\mu_{n}$ are $q_{n}/\sigma$ (i.e., the SNR) and $\alpha$ (i.e., temporal signal dependence). Here, we present a numerical example for clarification. Fig. 7A and B shows the plot of $\mu_{n}$ vs. $\alpha$ and $\kappa_{n}$ vs. $\alpha$ for $n=50$ and $q/\sigma^{2}=10$. It is apparent that $\mu_{n}$ increases with $\alpha$ and does not exceed $1$. The fact that $\mu_{n}<1$ implies that the DBMT estimator achieves a higher noise suppression compared to the standard MT method. This fact agrees with the noise suppression performances observed in Section IV.

Fig. 7B shows the plot of $\kappa_{n}$ vs. $\alpha$, which exhibits a similar increasing trend, but eventually exceeds $1$. This result implies that with a careful choice $\alpha$, it is possible to achieve $\kappa_{n}<1$, and hence obtain lower variance than that of the standard MT estimate. Fig. 8 illustrates this statement by showing the value of $\alpha$ for which $\kappa_{n}\approx 1$ vs. $q_{n}/\sigma^{2}$. For models with high temporal dependence (i.e., $\alpha$ close to $1$), it is possible to achieve $\kappa_{n}<1$ and hence reduce the estimator variance, due to the increase in the weight of data pooled from adjacent windows, even for small values of $q_{n}/\sigma^{2}$ (i.e., low SNR). However, this reduction in variance comes with the cost of increasing the bias. On the contrary when the data across windows have low temporal dependence (i.e., $\alpha\ll 1$), it is only possible to achieve a reduction in variance for high values of $q_{n}/\sigma^{2}$ (i.e., high SNR). This is due to the fact that at low SNR with low temporal dependence, pooling data from adjacent windows is not beneficial in reducing the variance in a particular window, and indeed can result in higher bias.

Remark 3. Even though the parameter $\alpha$ is estimated in a data-driven fashion, it can be viewed as a tuning parameter controlling the bias-variance trade-off, given the foregoing discussion. For a given SNR, fixing $\alpha$ at a small value can help reduce the variance but with a cost of increasing bias, and vice versa. In light of this observation, Fig. 8 can be thought of as a guideline for choosing $\alpha$ to achieve $\kappa_{n}\approx 1$, so that the estimator is nearly unbiased, and achieves a lower variance than that of the standard overlapping MT estimator due to the noise suppression virtue of the state-space model. Although we focused on the case of flat spectrum in the foregoing discussion, it is possible to numerically compute these trade-off curves for more general cases, given the general expressions for $\kappa_{n}(f)$ and $\mu_{n}(f)$ given in Appendix A.

Remark 4. Extending Theorems 1 and 2 to the log-DBMT algorithm is not straightforward, due to the high nonlinearity of the underlying state-space models. However, under the common Gaussian approximation of the log-posterior density, the well-known variance reduction property of the fixed interval smoother [41, 51] carries over to the estimator of $\log S_{n}(f)$. That is, the variance of the log-DBMT estimate of $\log S_{n}(f)$ obtained using the state-space model is lower than that of the standard MT estimate which only uses the data within window $n$. This fact agrees with our earlier observation in Section IV-A regarding the tightening of the confidence intervals for log-DBMT as compared to the overlapping MT estimates.

VI Concluding Remarks

Spectral analysis of non-stationary time series data poses serious challenges for classical non-parametric techniques, in which temporal smoothness of the spectral representations are implicitly captured using sliding windows with overlap. This widely-practiced approach tends to ignore the inherent smoothness of the data and is not robust against measurement/sampling noise. In this paper, we address these issues and provide an alternative to the sliding window spectrogram analysis paradigm. We propose two semi-parametric spectrogram estimators, namely the DBMT and log-DBMT estimators, by integrating techniques from MT analysis and Bayesian estimation. To this end, we explicitly model the temporal dynamics of the spectrum using a state-space model over the spectral features obtained by multitapering. Therefore our algorithms inherit the optimality features of both Bayesian estimators and MT analysis.

Our algorithms admit efficient and simple implementations, thanks to the Expectation-Maximization algorithm and the well-known fixed interval state-space smoothing procedure. Unlike existing approaches, our algorithms require no a priori assumptions about the structure of the spectral representation and operate in a fully data-driven fashion. While both algorithms yield spectral estimates that are continuous in time, by design the DBMT algorithm significantly suppresses the measurement noise in forming the spectrogram and the log-DBMT algorithm mitigates sampling noise due to small observation length. We establish the performance gains provided by our algorithms through theoretical analysis of the bias-variance trade-off, as well as application to synthetic and real data from human EEG recording during sleep and ENF signals from audio recordings.

Appendix A Proofs of Theorems 1 and 2

Recall from Eq. (29) that the $k$th eigen-coefficient estimate at window $n$ can be written in terms of the observed data as:

[TABLE]

where $\widetilde{\mathbf{y}}_{s}^{(k)}=\mathbf{u}^{(k)}\odot\widetilde{\mathbf{y}}_{s}=\mathbf{U}^{(k)}\widetilde{\mathbf{y}}_{s}$. Given that $\mathbf{Q}$ is a diagonal matrix with elements $q(f_{j})$, $j=1,2,\cdots,J$, it can be shown that $\mathbf{\Sigma}_{\infty}$, $\mathbf{\Lambda}$ and $\mathbf{\Gamma}$ are also diagonal matrices. Denoting the elements of $\mathbf{\Sigma}_{\infty}$, $\mathbf{\Gamma}$ and $\mathbf{\Lambda}$, respectively by $\tau(f)$, $\eta(f)$ and $\gamma(f)$, for $f\in\{f_{1},f_{2},\cdots,f_{J}\}$, we have $\eta(f)=\frac{\alpha^{2}\tau(f)+q(f)}{1+rW(\alpha^{2}\tau(f)+q(f))}$.

A-A Proof of Theorem 1

First note that Eq. (30) implies that $\mathbbm{E}[{dz_{n}(f)dz_{n+t}^{*}(f^{\prime})}]=\alpha^{t}D(f)\delta(f-f^{\prime})dfdf^{\prime}$. By invoking the Cramér representation, the covariance of the data tapered by the $k$th and $l$th dpss sequences can be expressed as:

[TABLE]

From Eqs. (20) and (31), we get:

[TABLE]

Taking the expectation of both sides and after some simplification, one arrives at:

[TABLE]

Using the orthogonality of the PSWFs as in [52], and using the fact that $D(\beta)\leqslant\sup_{f}D(f),\forall\beta$, we get:

[TABLE]

[TABLE]

$\text{and }\mu_{n}(f):=\eta(f)^{2}\sum_{s=1}^{N}\gamma(f)^{2|s-n|},$ for $f\in\{f_{1},f_{2},\cdots,f_{J}\}$. Using the triangle inequality, the bound of Theorem 1 on \big{|}\mathbbm{E}[\widehat{D}_{n|N}(f)]-D(f)\big{|} follows. $\blacksquare$

A-B Proof of Theorem 2

Using the notation of Appendix A, we have:

[TABLE]

Note that we have omitted the integral limits, as they are understood to be same as in (1) henceforth. After summing over all tapers and rearranging the summations, the first expression within the brackets in (A-B) becomes:

[TABLE]

Let

[TABLE]

Using the Schwarz inequality, the integral in Eq. (A-B) can be bounded by:

[TABLE]

Using bounds on the convolutions of PSWFs from [52], and upper bounding $\mathbbm{A}(n,m,\beta)$ by $\sup_{f}\{\kappa_{n}(f)D(f)+\mu_{n}(f)\sigma^{2}\}$, the statement of the theorem on the variance of the DBMT estimate $\widehat{D}_{n}(f)$ follows. $\blacksquare$

Appendix B Characterization of $\kappa_{n}(f)$: bounds and parameter dependence

B-A Lower and Upper Bounds on $\kappa_{n}(f)$

Consider the scenario where $q(f)=q$, i.e., flat spectrum. Then, the dependent of $\gamma(f)$ and $\kappa_{n}(f)$ on $f$ is suppressed. We have the following bound on $\kappa_{n}$:

Proposition 3**.**

For $0<\gamma,\alpha<1$, the quantity $\kappa_{n}$ can be bounded as:

[TABLE]

where $\displaystyle T_{0}:=\frac{1+\gamma^{2}-\gamma^{2n}-\gamma^{2(N-n+1)}}{1-\gamma^{2}}$.

Proof:

To get an upper bound on $\kappa_{n}$, we rewrite the expression defining $\kappa_{n}$ as:

[TABLE]

Now, let us define $T_{t}:=\sum_{s=1}^{N}\gamma^{|s-n|}\gamma^{|s-t-n|}$. Then it can be verified that:

[TABLE]

and

[TABLE]

Also, we have:

[TABLE]

Similarly, we have:

[TABLE]

which along with (B-A) leads to the claimed upper bound. For the lower bound, we use the fact that

[TABLE]

which implies $\sum_{t=0}^{N-1}\alpha^{t}T_{t}\geqslant\frac{1-(\alpha\gamma)^{N}}{1-\alpha\gamma}T_{0}$ and $\sum_{t=-N+1}^{0}\alpha^{|t|}T_{t}\geqslant\frac{1-(\alpha\gamma)^{N}}{1-\alpha\gamma}T_{0}$. Using the latter lower bounds for $\kappa_{n}$ yield the claimed lower bound. ∎

B-B Relation between $\kappa_{n}$ and $\mathbf{Q}$

The expressions for $\kappa_{n}(f)$ and $\mu_{n}(f)$ in Appendix A depend on $\gamma(f)$ and $\alpha$. But $\gamma(f)$ itself depends on $q(f)$ and $\alpha$, and it is not straightforward to give a closed-form expression of $\gamma(f)$ merely in terms of $q(f)$ and $\alpha$, since it requires computation of the filtered error covariance matrix $\mathbf{\Sigma}_{n|n}$ given $\mathbf{Q}$. Again, by invoking the stead-state approximation, and defining $\mathbf{\Sigma}_{n|n-1}=:\mathbf{\Sigma},\forall n=1,2,\cdots,N$, the matrix $\bm{\Sigma}$ can be obtained by solving the following algebraic Riccati equation:

[TABLE]

and thereby the steady-state error covariance matrix $\bm{\Sigma}_{n|n}=:\bm{\Sigma}_{\infty}$ is given by $\mathbf{\Sigma}_{\infty}=\frac{1}{\alpha^{2}}(\mathbf{\Sigma}-\mathbf{Q})$.

Although this procedure can be carried out numerically, in general it is not possible to solve the Riccati equation to get a closed-form expression for arbitrary $\mathbf{Q}$. In order to illustrate the explicit dependent of $\kappa_{n}(f)$ on the state-space parameters, we consider the case of flat spectrum where $\mathbf{Q}=q\mathbf{I}$. In this case, it can be shown that $\mathbf{\Sigma}=\zeta\mathbf{I}$ for some $\zeta>0$ and the matrix equation (45) reduces to a simpler scalar equation for $\zeta$, given by:

[TABLE]

Following simplification, a quadratic equation for $\zeta$ results, and since $\zeta\geqslant 0$, the positive solution for $\zeta$ is given by:

[TABLE]

Then, $\gamma$ can be computed as \gamma=\frac{1}{\alpha}\big{(}1-\frac{q/\sigma^{2}}{\zeta/\sigma^{2}}\big{)}. Using these expressions for $\gamma$ and $\alpha$, the functions $\kappa_{n}$ and $\mu_{n}$ can be computed. Fig. 7B shows the upper and lower bounds on $\kappa_{n}$ evaluated for $n=50$ and different values of $\alpha$ for $q/\sigma^{2}=10$. The upper and lower bounds are simple functions of $\alpha$ and $q/\sigma^{2}$ and can be used to further inspect the performance trade-offs of the DBMT algorithm with respect to the state-space model parameters.

Acknowledgment

We would like to thank Adi Hajj-Ahmad and Min Wu for providing us the ENF data from [49].