Hybrid Adaptive Parametric Frequency Analysis

Introduction: High quality spectra are very useful in anesthesia related procedures where Electroencephalogram (EEG) frequency content has been shown to drastically help track different brain states. A recent work (Konstantinidis & Brown, 2019 [1]) introduced the Gaussian Hybrid Autoregressive Model as a parametric method to generate smooth, very high resolution spectrograms of non-stationary EEG data of humans under propofol. Objective: In this paper, we extend the model proposed in [1] to incorporate non-Gaussian state noise. Methods: A Monte Carlo Markov Chain (MCMC) filtering procedure on a self-organizing state-space model is presented. Results: We test the extended model on EEGs from human patients under propofol, ketamine and sevoflurane and illustrate the advantages over its Gaussian counterpart. Conclusion: The suitability of the proposed method for online use, in combination with its ability to smoothly track frequency changes in human EEG signals under the most common anesthetics, suggests that it can be used for real-time brain state tracking. Such online use can facilitate the design of more precise closed loop systems for automatic control of brain states under general anesthesia.


Introduction
Non-stationary time-series frequency analysis traditionally attracts the attention of the scientific community. This interest is due to its enormous practical benefits that arise from the fact that most signals with a temporal representation encountered in real-life applications are characterized by time-varying statistics. An example of a real-life application is the use of EEG to track brain states under anesthesia. EEG has been extensively used to characterize brain states under various anesthetics [2], [3], [4] and to automatically adjust drug infusion rates using closed loop control systems using instantaneous brain states as control signals [5]. The issue of non-stationarity is overcome by assumptions about local stationarity properties over short time intervals and application of techniques suitable only for stationary signals. Even though such assumptions # This work was supported by the The Bertarelli Program in Translational Neuroscience and Neuroengineering (https://bertarelli.hms.harvard.edu/) * Corresponding author. Email: kritonkonstantinidis@gmail.com can have practical advantages, they are not always sufficient and can be a source of noise. Current techniques used to estimate frequency spectra of non-stationary signals include various non-parametric methods, such as the multitaper approach developed in [6] and variants thereof. Non-parametric approaches cannot yield spectrograms with both high temporal and frequency resolution simultaneously. Improving time resolution by using a shorter, assumed stationary, data segment results in lower frequency resolution and vice versa. Parametric methods (e.g., autoregressive models) are usually based on time-varying linear predictive models and provide an elegant alternative. In [1], a hybrid autoregressive model is used in a state-space form with Gaussian state noise and Kalman Filtering is employed to fit the model. Hybridity is shown to yield smoother EEG spectrograms than non-parametric and discrete parametric methods. In [7], Cauchy state noise is used and is shown to provide advantages over the use of Gaussian state noise in efficiently tracking EEG frequency changes during Event Related Synchronization (ERS). The scale of the Cauchy distribution is assumed identical for all parameters and is jointly estimated with the state (autoregressive parameters) using a sequential importance resampling algorithm. Motivated by the work of [7], this paper extends the hybrid approach of [1] and develops an algorithm to fit hybrid autoregressive models using Cauchy state noise by jointly estimating a different scale for each autoregressive parameter. Finally, it explores the performance of the hybrid algorithm on different EEG datasets of humans under anesthesia and compares it to that obtained by other widely used non-parametric and discrete parametric methods.

Hybrid Cauchy Model
In [7], the performance of the filter is better than the Gaussian one in identifying abrupt frequency changes but the use of Cauchy noise can introduce some disturbances in more stationary regimes of simulated data. Motivated by the approach in [7], a bootstrap hybrid algorithm is developed, that attempts to improve the tracking performance of the Gaussian filter while mitigating outliers in more stationary parts, by taking advantage of hybridity.

Self-Organizing State-Space Formulation.
State Equation: where w j (t) ∼ C(0, q j (t)) is a continuous Cauchy process with 0 mean and scale parameter q j (t) and j (t) ∼ N (0, R ) is a continuous Gaussian process with 0 mean and variance R v k ∼ N (0, R(t)) is a discrete scalar white Gaussian observation noise with 0 mean and variance R(t) x(t) = a 1 (t), a 2 (t), . . . , a p (t), ln(q 1 (t)), ln(q 2 (t)), Let T denote the total number of observations in the data, indexed by k, N denote the total number of particles used for the sequential importance resampling, indexed by i, j index the autoregressive parameters a 1 (t), a 2 (t), . . . , a p (t) and f s denote the sampling frequency.

Discretization of the state equation.
To start the algorithm, it is necessary to discretize equation (1).
For a discretization step ∆, (1) can be written as . The Ito integral in equation (3) can be discretized as , for the first p dimensions where Q is the diagonal matrix with entries q j (t) and ∆w(t) ∼ N (0, Q q ∆), where Q q is the diagonal matrix with entries j (t) for the next p dimensions. The Self-Organizing state-space equations can now be written as: Log-likelihood of the model. Let θ denote the scale parameters q j . The log-likelihood that the observed data are generated by the model is defined as k denotes the particle i at time k. Using (7), (6) can be approximated as k is the unnormalized weight of particle i at time k. 2 EAI Endorsed Transactions on Bioengineering and Bioinformatics Online First Hybrid Adaptive Parametric Frequency Analysis

Model Selection, Spectral Estimation & Roughness Metric
Model selection is performed using the well-known Akaike Information Criterion (AIC), defined as AIC(p) = 2p − 2L ML , where p is the order of the autoregressive model and L ML is the maximized marginal log-likelihood for that model. The order that yields the lowest AIC will be chosen. L ML is approximated by (8). Model selection was performed in an initial burn-in period of the anesthesia procedure (10 seconds for the experiments presented).
Having calculated the optimal estimates for the autoregressive parameters, the spectral density at frequency f and time t is calculated using: Where a k (t) is the k-th estimated autoregressive coefficient at time t.
Finally, we define a roughness metric of a coefficient a k (t) between 2 time points t 1 and t 2 , R k = [8]. The lower the metric, the smoother the temporal evolution of the autoregressive coefficients.

Results
Results on non-parametric, discrete Gaussian (DC) and hybrid Gaussian (HG) and Cauchy (HC) models will be discussed. Discrete Cauchy models lead to even noisier spectrograms than the Gaussian models, as expected due to the nature of the heavy-tailed Cauchy distribution. Spectrograms are shown for the discrete Gaussian and hybrid Cauchy models.

Simulation Results
To validate the extended model and to obtain an insight on how hybridity and use of Cauchy state noise enhances smoothness and instantaneous frequency estimation accuracy, it was fit along with its Gaussian counterparts to a simulated noisy sinusoidal wave with stepwise frequency modulation, in order to simulate both abrupt and stationary regimes.
An artificial sinusoidal signal of T = 60sec, z k = A k sin(ω k k f s ) + v k was generated at a typical EEG sampling rate f s = 250Hz. v k is added white Gaussian noise with zero mean and unit variance. Amplitude evolves as EAI Endorsed Transactions on Bioengineering and Bioinformatics Online First Let c k|N denote the roots of the characteristic polynomial of the autoregression, which can be written as c k|N = r k|N e −iω k|N , where r k|N is the modulus and ω k|N is the phase of each root c k|N . The dominant frequency is proportional to the phase of the complex roots: 2π . Applying this procedure for a discrete and hybrid models on the simulated data, using an autoregressive model of order p = 2, R = 1 and Q = 10 −3 I, B = 100, Q q = 10 −4 I the following estimations are obtained (Figure 1). A low autoregressive order was used for increased tractability of the poles. The values for the remaining parameters were set empirically.
All models are able to track the true frequency evolution, but the hybrid ones result in much less fluctuation around the true instantaneous value. Hybrid Cauchy Model does not suffer from the lag of the Hybrid Gaussian model at the times of abrupt frequency changes. The use of Cauchy model in combination with the adaptive estimation of the scales allows for very high adaptability. Mean Square Errors (MSE) are shown  Table 1. Values shown are averages and standard deviations for 100 runs.

Spectrograms of EEG under anesthetics
Spectrograms for propofol, ketamine and sevoflurane were calculated using the hybrid and discrete parametric methods, as well as widely used non-parametric methods including the periodogram, the multitaper spectrogram, the state-space periodogram and the recently developed state-space multitaper spectrogram [9]. For the multitapers, 3 tapers were used at a 2Hz spectral resolution. The time window of assumed stationarity was 2 seconds. For all Cauchy spectrograms, B = 10 bootstrap samples and N = 500 particles were used. The performance of the proposed model was compared to that of discrete methods that have been used before for EEG spectrogram generation under anesthesia [10] and to that of the non-parametric approaches mentioned above. All EEG recordings are part of deidentified data collected from patients at Massachusetts General Hospital (MGH), as a part of a MGH Human Research Committee approved protocol. Spectrograms of human EEG under propofol, ketamine and sevoflurane are discussed.
Propofol. In the case of propofol, non-parametric spectrograms ( Figure 2) clearly show an α frequency band in addition to the slow oscillations, as expected [2].
The periodogram and the multitaper spectrogram are quite noisy. On the other hand, their state-space counterparts seem to underfit as it is very implausible that the α band is evolving as a completely straight line in an actual EEG experiment.
Autoregressive models of selected order p = 14 were used to calculate the parametric spectrograms. In the purely discrete case (Figure 3), the resulting spectrogram seems able to identify the subtle changes in frequency but is not able to reduce the underlying 4 EAI Endorsed Transactions on Bioengineering and Bioinformatics Online First Hybrid Adaptive Parametric Frequency Analysis

Ketamine.
A β-γ frequency band is apparent in all non-parametric spectrograms ( Figure 5) with peaks at around 30Hz and 50Hz, as expected [2]. The periodogram and multitaper spectrogram seem to identify changes across the bands but they suffer from high noise. On the other hand, as it was the case with propofol, state-space periodogram and multitaper spectrogram oversimplify the frequency evolution and yield a constantly evolving β-γ band.
Order p = 11 was selected for the parametric spectrograms. The discrete one ( Figure 6) is quite noisy, while the Cauchy hybrid spectrogram (Figure 7) provides a better result.  The ability to efficiently track high frequency events due to ketamine induced spikes in the EEG is a strong validation for the proposed model.

Sevoflurane.
Non-parametric spectrograms of EEG under sevoflurane follow (Figure 8). The interchange between the merge and separation of the α and slow oscillations is observed throughout the whole recording, as expected [2].
Order p = 18 was selected for the parametric spectrograms. The hybrid model ( Figure 10) provides again the best result by efficiently smoothing the discrete spectrogram (Figure 9).
The conclusion drawn by visual inspection for smoother spectrograms of the hybrid models is corroborated by the roughness metrics of Table 2, as defined in (10). The reported final value for each model is the average of the metrics for each individual coefficient taking part in the model. 5 EAI Endorsed Transactions on Bioengineering and Bioinformatics Online First Konstantinidis, Brown   (4.08 ± 0.034) · 10 5 (4.74 ± 0.14) · 10 5 Ketamine (8.36 ± 0.18) · 10 4 (1.03 ± 0.083) · 10 5 Sevoflurane (9.01 ± 0.09) · 10 5 (1.03 ± 0.011) · 10 6 As expected, the hybrid Cauchy model yields less smooth coefficients than the hybrid Gaussian but much smoother than the discrete Gaussian. We believe that the hybrid Cauchy model estimates the instantaneous frequency with higher accuracy, balancing the tendency for oversmoothing of the hybrid approach by the heavier tales of the Cauchy distribution and thus the proposal of larger variations of the autoregressive parameters.
Apart from the better performance in simulated data, log-likelihood values of Table 3 corroborate the aforementioned hypothesis. Averages and standard deviations are shown for 100 runs. Cauchy models yield higher likelihood values in all three cases.
Finally, the Ljung-Box scores for residual autocorrelation are shown in Table 4. In accordance with [7], the scores are much higher than the cutoff value for 5 lags (11.07), indicating strong residual autocorrelation. This is a common phenomenon when fitting models to

Online use
All algorithms described in this paper can be a great tool for offline analysis of the EEGs. On the other hand, we are ultimately interested in developing real-time control systems and therefore the practicability of the algorithms for online use should be investigated. In this section, we report the time elapsed for the algorithms to be run on the 10-minute dataset of EEG under ketamine used to generate spectrograms of Figures 7-9. Run on a 2.9 GHz Intel Core i5 with an 8GB RAM the following execution times were obtained for the hybrid and discrete models ( Table 5). The algorithm was run five times for each case to establish an indicative average and standard deviation. The frequency resolution was 4 index/Hz. Autoregressive model selection was performed before the execution of the filtering procedure and is not included in the reported times. Model order selection can be determined through offline analysis of the EEG, or through the use of a short burn-in period in the very beginning of the anesthesia procedure. All reported times are much less than the duration of the recording (10 minutes) and therefore the algorithms are suitable for online use (the time required to process each new incoming observation is less than the sampling rate).
Longer execution times are observed in the hybrid models in comparison to the discrete versions. This is expected since the integration of the differential equations for the state vector and the state covariance matrix, that is the main idea of this approach, takes more time than just "jumping" from a value at time t k to the updated value at t k+1 . Furthermore, very high 6 EAI Endorsed Transactions on Bioengineering and Bioinformatics Online First temporal resolution of the parametric models leads to an inherent higher computational burden, that the nonparametric models avoid due to the binning procedure. Finally, regarding the Cauchy model, the use of the MCMC approach in combination with the hybridity and thus the need to sample from the discretized model (∆ = f s 50 in this paper) leads to the longest execution times. On the other hand, this method is shown to provide the most accurate results in the simulated data and the highest likelihood values in human EEG datasets. Thus, it would be beneficial to adapt the algorithm for higher computational efficiency [11].

Discussion
The hybrid filter stands at the border between the underfitting regime of the state-space multitaper periodogram and the overfitting regime of the multitaper spectrogram. In relation to the purely discrete version, the hybrid filter is better able to capture the smoothness in the data, while at the same time this ability does not prevent it from identifying sharp frequency changes. To stress the importance of the hybrid model, we would like to emphasize that lowering the covariance matrix of the discrete model as a measure towards smoothness of EEG spectrograms of patients under propofol was not an efficient approach, as shown in [1]. Lower quality spectrograms with induced artifacts due to a low covariance matrix were also observed for EEGs under ketamine and sevoflurane.
In this hybrid approach, no prior knowledge about the evolution of the parameters in the time intervals between the observations is assumed. As a result, the autoregressive coefficients are assumed to follow a random walk. The proposed model can become even more potent in calculating frequency spectra if some prior knowledge about the parameters is available that will guide a reasonable choice of the state transition matrix. On the other hand, predicting a priori how the parameters are expected to evolve is not a trivial issue and demands a fair amount of research and experience with EEG datasets and effects of anesthetics on the brain. Nevertheless, it is expected that incorporating such prior information can enhance the model's performance even more, since by integrating the state equation, the model will be able to predict values for the parameters even in the time intervals between the observations where no new information is available.
Concerning the two different hybrid models, our results agree with [7] and indicate that the use of Cauchy noise is also advantageous when calculating spectra of EEG under anesthesia, in addition to EEG spectral analysis for Event Related Desynchronization that is discussed in [7]. The use of the self-organizing model allows the scales of the Cauchy distribution to be adjusted to higher values during abrupt frequency changes for better tracking. On the other hand, fitting of Cauchy models is more computationally expensive due to the need for Monte Carlo simulations. Gaussian assumptions lead to convenient closed form solutions and allow the use of Kalman Filter that is perfectly suited for online use and do not impose the high computational burden of Monte Carlo simulations, which, even though can be alleviated by parallelization, requires substantial hardware. For this reason, the Gaussian [1] or the Cauchy model presented in this paper are available to allow for optimal trade-off between desirable accuracy and use of resources one wants to allocate.

Conclusion
Overall, the proposed method is an accurate tool for frequency analysis of non-stationary data. It provides a complete framework that incorporates robust autoregressive model selection and an efficient way to calculate instantaneous frequency estimates. To the best of our knowledge, no other method has shown similar, smooth yet accurate, results for frequency analysis of EEG under anesthetics. The suitability of the proposed model for online use, coupled with its ability to provide smoother spectrograms, suggests that it can be used as a monitoring tool for real time brain state tracking under general anesthesia that will further potentiate the design of closed loop systems for automatic and precise control of brain states. Future work will test the accuracy of brain-state classification when using our model for feature extraction. Finally, note that even though we have focused on anesthesia EEGs, the work presented in this paper can be applied to any non-stationary time-series.