Schedule for: 22w2230 - Multitaper Spectral Analysis (Online)

"In this talk, I trace the development of "multitaper spectrum analysis'' from stumbling over Slepian and Pollak's 1961 paper during my first year as an undergraduate, through my initial, circa 1973, multitaper code, to the present. I had a few conference papers prior to the 1982 paper but some aspects of data analysis, such as autocorrelations and nonstationarity were only given a line or two.

Timelimited (or spacelimited) functions cannot be simultaneously bandlimited (in frequency). Yet the finite precision of measurement and computation unavoidably bandlimits our observation and modeling scientific data, and we often only have access to, or are only interested in, a study area that is temporally or spatially bounded. In the geosciences we may be interested in spectrally modeling a time series defined only on a certain interval, or we may want to characterize a specific geographical area observed using an effectively bandlimited measurement device. Here, we give a theoretical overview of one approach to this "concentration" problem, as proposed for time series by Slepian and coworkers in the 1960s, and made famous by David Thomson in the 1980s and beyond. We show how this framework leads to practical algorithms and statistically performant methods for the analysis of signals (linear inverse theory) and their power spectra (quadratic spectral estimators) in one and two dimensions, and for scalar and vectorial signals defined on the surface of a unit sphere. In particular, we develop procedures to invert for planetary potential fields from vector observations collected within a spatially bounded region at varying satellite altitude. Our method relies on the construction of spatiospectrally localized bases of functions that mitigate the noise amplification caused by downward continuation while balancing the conflicting demands for spatial concentration and spectral limitation. Practical examples will be drawn from the analysis of gravitational potential data, vectorial measurements of internal and external planetary fields, and other applications.

In the mid 1990s, Jeffrey Park of Yale University and I developed the “Multitaper Method-Singular Value Decomposition” or “MTM-SVD” as an approach for detecting and reconstruct the characteristics of narrowband spatiotemporal signals immersed in colored noise. I will review the approach and its application to climate observations, paleoclimate data and climate model simulations to the detection of a putative low-frequency “Atlantic Multidecadal Oscillation” or “AMO” signal in the climate record.

In 2014, a 13-element broadband array was installed at Pinyon Flat Observatory to test installation options and characteristics for the EarthScope USArray TA project in preparation to deployment in Alaska. In this array, each element uses a broadband posthole sensor (STS-5 or Trillium 120PH) and either a Setra 278 barometer or an Episensor strong motion accelerometer. Each posthole is about 3 meters deep in cased holes emplaced in decomposed granite by a light weight portable drill developed for the TA using a hammer bit. Continuous data from has been collected at 1, 40, and 200 sps on the broadband sensors. The array has an aperture of 1 km with a minimum station spacing of 40 meters. This array presents a unique opportunity to study the noise and coherence characteristics of the seismic signals from this new generation of broadband seismic installation technique. We will present results from multitaper spectra from seismic noise as well as from local, regional, and teleseismic events recorded on the array. We will complement the spectra results by presenting multitaper coherence analysis between various combinations of sensors.

(This research is joint with Shuhan Tang and Yunzhang Zhu, The Ohio State University.) Spectral estimation provides key insights into the frequency domain characteristics of a time series. Naive nonparametric estimates of the spectral density, such as the periodogram, are inconsistent, and the more advanced lag window or multitaper estimators are often still too noisy. We propose an L1 penalized quasi-likelihood Whittle framework based on multitaper spectral estimates which performs semiparametric spectral estimation for regularly sampled univariate stationary time series. Our new approach circumvents the problematic Gaussianity assumption required by least square approaches and achieves sparsity for a wide variety of basis functions. We present an alternating direction method of multipliers (ADMM) algorithm to efficiently solve the optimization problem, and develop universal threshold and generalized information criterion (GIC) strategies for efficient tuning parameter selection that outperform cross-validation methods. Theoretically, a fast convergence rate for the proposed spectral estimator is established. We demonstrate the utility of our methodology on simulated series and to the spectral analysis of electroencephalogram (EEG) data.

The SS seismic phase is a shear wave that reflects at Earth’s surface at the source-receiver mid-point. SS precursors reflect at the underside of sharp internal interfaces below the raypath midpoint; deeper interfaces lead to larger precursory time advances. SS precursors can target features in the mantle transition zone (MTZ) far from both earthquake sources and seismic receivers. To resolve interfaces from scattered waves in the time-domain, current SS precursor techniques require tightly bandpassed signals (e.g., 0.02-0.05 Hz), limiting resolution. Higher frequency content would allow for the detection of finer layering in and around the MTZ. Additionally, current methods often produce side lobes around large signals that hamper the identification of small-scale features and low-velocity zones surrounding the MTZ. We benchmark a new SS precursor deconvolution technique, based on multiple-taper correlation (MTC), on full-waveform synthetic seismograms computed via AxiSem3D for a 1-D reference Earth model. MTC deconvolution of SS precursors can increase the frequency cutoff up to ~0.5 Hz, thereby increasing vertical resolution to ~10 km. The low-frequency bandpass cutoff can be lowered or eliminated, reducing the side lobes of large pulses. We apply our MTC SS precursor deconvolution to ~7000 seismograms recorded at sensors of the Global Seismographic Network with source-receiver mid-points in the North-Central Pacific Ocean. We resolve a distinct SS precursor phase from above the MTZ that had previously been inferred from an asymmetry in sidelobe amplitudes, demonstrating it to be an interface (<10-km thickness), rather than a thick wavespeed gradient.

Eigenvectors of the operator that truncates to a neighborhood of the origin then bandlimits to a neighborhood of zero in the spectral domain are investigated in the context of finite-dimensional Boolean hypercubes. The eigenvectors of these spatio--spectral limiting operators can be regarded as analogues of Slepian vectors---finite-dimensional versions of zero order prolate spheroidal wave functions. Unlike the linear case, in the hypercube setting the composition of truncation in the spatial and spectral domains does not commute with what might be considered the natural analogue of the differential operator whose eigenfunctions are the prolate functions. Thus a different approach is needed to identify and compute the eigenvectors of this composition. Such an approach is outlined here. Potential application to spectrum estimation of so-called (local) dyadic processes in a manner analogous to Thomson's multitaper method is outlined. Extension of this approach to discrete tori is also outlined. This talk is based on join work with Jeffrey A. Hogan.

This study has two goals: (1) measure solid Earth diurnal and semi-diurnal tidal constituents using multitaper spectral analysis, and (2) calculate the time delay of the Earth’s solid tidal response between two land-based geographic locations. For these goals, high-precision tidal signals were collected from borehole strainmeters at Plate Boundary Observatory (PBO) stations in Yellowstone, WY, USA and Sequim, WA, USA. Signals from synchronized microstrain time series sampled at 5 minute intervals from 1 January 2011 to 31 December 2015 were evaluated with power spectra and harmonic analysis with statistical F-testing. Along with the usual major constituents predicted for the Earth tides, at Yellowstone, WY non-tidal components related to episodic lake seiches are present, and at Sequim, WA lines are detected at frequencies associated with ψ1 and ϕ1 tides. The K1 and K2 frequencies at the two stations indicate delays in excess of those predicted from a solid Earth tidal model for the two stations due to ocean tidal effects.

Although multitaper spectral analysis of second-order stationary (s.o.s.) time-series is well established, some of the key characteristics of neural data such as non-stationarity and non-linearity of neuronal activity cannot be readily captured by existing techniques. For instance, the classical sliding window multitaper method used in analyzing non-stationary time series may undermine the underlying dynamics and lacks a precise statistical characterization of the estimates. Furthermore, existing methods for multitaper analysis of neuronal spiking observations require estimation of latent time-series that underlie spiking activity from the spiking statistics. The time-domain smoothing in the estimation procedure inadvertently introduces distortion in the subsequent spectral-domain analysis. Here, we address these issues by explicitly modeling the spectral dynamics using a state-space model where the underlying states pertain to the eigen-spectra of a quasi-stationary process. Once these states are estimated within a Bayesian framework, the spectrograms and their respective confidence intervals can be readily constructed. Next, we propose a multitaper variant tailored for neural spiking data. Assuming an s.o.s. latent process, we construct auxiliary spiking statistics akin to tapered data following point process theory, from which the corresponding eigen-spectra can be directly inferred. Finally, we develop a unified framework to infer multitaper semi-stationary spectra of neuronal ensembles from multivariate spiking observations, by combining the aforementioned approaches. We demonstrate the performance gains of our proposed methods over existing techniques through extensive simulations and application to experimentally recorded neuronal data, as well as theoretical analyses of the bias-variance trade-off.

While multitaper analysis stands on some substantial theoretical foundations, current literature does little to address the following questions: why does it appear so robust, and to what extent are the goodness-of-fit metrics good decent proxies for model performance? The early 20th Century saw the construction of a time-series model class amenable to spectral analysis for processes containing a nonrandom, harmonic signal buried in stationary, Gaussian noise. Since the turn of the millennium, the analysis of complex-valued stochastic processes has remained an active field in signal processing to meet the demands of high-quality inference with the advent of high-quality big data. Considerable attention has been paid to cyclostationary processes, which have themselves been used to augment the model class amenable to spectral analysis. Current literature posits that periodicity of a cyclostationary process is accurately explained by a small fraction of the Slepian eigencoefficient processes of the discrete-Fourier-transform process. Moreover, these spectral processes are assumed complex-Gaussian. Justification for this has typically been through goodness-of-fit metrics, or through citations without a comprehensive illustration. This talk introduces a model class containing both cyclostationary and stationary processes that are amenable to multitaper spectral analysis, revealing how the treatment of cyclostationary rather stationary noise requires a revised distribution theory to justify goodness of fit: from the choice of a suitable dimension and reliability of the complex-Gaussian assumption; to what random spectral quantities are being considered whose parameters are inferred; and to establishing which physical mechanisms can give rise to members of the model class.

This work determines how to define a discretely implemented Fourier transform when analysing an observed spatial point process. To develop this transform we answer four questions; first what is the natural definition of a Fourier transform, and what are its spectral moments, second we calculate fourth order moments of the Fourier transform using Campbell's theorem. Third we determine how to implement tapering, an important component for spectral analysis of other stochastic processes. Fourth we answer the question of how to produce an isotropic representation of the Fourier transform of the process. This determines the basic spectral properties of an observed spatial point process.

Maximum likelihood estimation of time series and spatial parameters is often intractable with large datasets and with non-trivial parametric covariance structures. In this talk I discuss the de-biased Whittle Likelihood - a method we recently developed that can estimate time series and spatial parameters for massive datasets using Fast Fourier transforms. The procedure is related to, but distinct from, the standard and well-known Whittle Likelihood - and I will make these distinctions clear in the talk. The key distinction is a de-biasing procedure to remove well known sources of bias from the standard Whittle likelihood. In addition, the de-biased Whittle likelihood can be generalised to simultaneously account for missing and multivariate data, as well as irregular sampling boundaries in higher-dimension random fields. We also show how tapering can be naturally incorporated into the framework which may lead to variance reduction of the parameter estimates in some instances. We have found numerous application benefits, and I will highlight results from our most recent project on estimating parameters of ocean wave models from significant wave height observations. The time series paper is available in Biometrika, and the spatial paper will shortly appear in Series B of The Journal of the Royal Statistical Society.

Thomson's multi-taper estimator has been extended to more general acquisition geometries with great success. Except in a few cases that exhibit useful symmetries, the calculation of the corresponding Slepian tapers is ill-posed and standard eigen-solvers may fail. We argue that this is not an obstacle to the computation of the multi-taper: while each tapered periodogram is impacted by the mentioned instability, their average is largely insensitive to it. This observation helps explain why common practice with irregular acquisition geometries yields remarkable accurate results, and also suggests an alternative method to compute the multi-taper based on block eigen-decomposition. Second, we derive risk rates for covariance estimation of stationary Gaussian processes with C2 densities, both for individual frequencies or in the spectral norm. When applied to the one-dimensional case, these show that Thomson's classical estimator achieves the optimal minimax rate. This analysis also contributes to the discussion on the optimal relation between the desired bandwidth and the number of tapers. Refs. [1] Abreu, L.D., Romero, J. L. MSE estimates for multitaper spectral estimation and off-grid compressive sensing. IEEE Trans. Inf. Theory. 63(12):7770–7776, 2017. [2] Andén, J., Romero, J. L. Multitaper estimation on arbitrary domains. SIAM J. Imaging Sci. 13, no. 3, 1565–1594, 2020. [3] Romero, J. L., Speckbacher, M. Spectral-norm risk rates for multi-taper estimation of Gaussian processes. J. Nonparametr. Stat. 34(2):448–464, (2022).

The remarkable work of D. Slepian, H. Landau and H. Pollak back in the 1960s at Bell Labs has had remarkable impact in many areas of mathematics, probably far beyond their dreams. I will describe a few of these searches that trace back their origin to their work. The word serendipity appears twice in D. Slepian's 1974 von Neuman lecture and recent developments prove him right once again.

This paper analyzes three 60 d sections of geomagnetic data from Honolulu Observatory during 2001−2, showing the ubiquitous presence of narrowband, very statistically significant, high Q features in multitaper power spectra, along with pervasive nonstationarity as measured by the frequency offset coherence over 400−4000 μHz (or 2500−250 s period). This behavior is nearly identical in the H and Z components of the geomagnetic field, and more subdued in the much weaker D component. The peak frequencies correlate well with the optically−measured frequencies of solar p−modes, and the raw Qs are defined by the resolution bandwidths of the estimates, with values ranging from hundreds to thousands. Further, spectral peaks are consistently coherent across frequency due to nonstationarity, and frequently exhibit cyclostationarity at offset frequencies of ±0.5 cycles per day. None of these characteristics are consistent with internal magnetospheric processes. A mixture central/noncentral chi square model was fit to raw spectral estimates in an attempt to model narrowband, high Q, quasi−deterministic modes embedded in a stochastic background. This model yielded noncentral fractions of 0.13 (1000−2000 μHz), 0.24 (1500−2500 μHz), 0.35 (2000−3000 μHz), 0.30 (2500−3500 μHz) and 0.17 (3000−4000 μHz). These values suggest that up to 35% of the power in the geomagnetic field in the 1000−4000 μHz band averaged over 60 days is forced by solar p−modes.

Elastic anisotropy in Earth’s interior arises from shear deformation associated with plate tectonics, mountain-building, aligned cracks or melt pockets, or fine-layered sediments. Many researchers have used the birefringence of P-to-S converted waves from the Moho discontinuity at the base of Earth’s crust to constrain its anisotropy. However, this practice ignores the substantial influence that anisotropy has on the initial amplitude of the converted wave. Theory predicts that P-to-S (and S-to-P) conversions to (and from) the two independent shear polarizations have a preferred “handedness” that varies with the approach angle (back-azimuth) of an earthquake wave. We deconvolve scattered waves with a coordinate transformation of the S-polarizations in the frequency domain, and form time-domain receiver functions (RFs) after an inverse DFT. RFs estimate the timing, polarity and amplitude of P-to-S (or S-to-P) converted waves caused by horizontal discontinuities in elastic properties. We project seismic signals onto the wrong “handedness” to estimate signal-generated noise. Using multitaper correlation (MTC), we find large amplitude variations with back-azimuth, consistent with P-wave and S-wave anisotropy. We confirm that variations are largest for crustal anisotropy with a tilted axis of symmetry, a geometry that is often neglected in birefringence interpretations. S-to-P scattering is not typically studied for anisotropy, partly because the SV and SH polarizations are often combined in earthquake signals and scatter differently. An improved simultaneous deconvolution has potential to test the hypothesis that the seismic lithosphere-asthenosphere boundary (LAB) is caused by a transition in anisotropic layering at the base of Earth’s tectonic plates.

Stars oscillate similar to musical instruments but at millions of closely-spaced frequencies. By estimating the spectrum underlying the time-series of a star, we can analyze the frequency structure of its oscillations (also called modes) and precisely estimate its age. However, using a biased and inconsistent spectrum estimator can cause false mode detections or imprecise mode frequency estimates, leading to problems with age estimation. The Lomb-Scargle periodogram, the most commonly used spectrum estimator in astronomy, faces these issues. We thus introduce multitaper spectral estimation in the frequency analysis of oscillation modes, a method known to tackle the problems of spectral leakage and noise. While multitaper spectral analysis is well-developed and widely used for regularly sampled data, its mapping to irregular sampling is limited. We develop a new multitaper method, mtNUFFT, for accurately analyzing irregularly-sampled asteroseismic and other astronomical data. We also introduce powerful features of multitapering (like the F-test for detecting purely periodic signals) into our new method. We are currently developing an open-source Python package for this method that generally applies to time-domain astronomy and other fields. In this talk, I will describe our method and its application to the age estimation of stars. I will also show that the new F-test helps detect exoplanets orbiting around stars. Finally, I will discuss the promising implications of my work for understanding the history of galaxies.

We confirm a fact discovered by Thomson [Spectrum estimation and harmonic analysis, Proc. IEEE, 1982]: assuming bandwidth W and N time domain observations, the average of the square of the first K = 2NW Slepian functions approaches, as K grows, an ideal band-pass kernel for the interval [−W, W]. We provide an analytic proof of this fact and measure the corresponding rate of convergence in the L1 norm. This validates a heuristic approximation used to control the MSE of the multitaper estimator.

The English word "spectrum'' was coined by Newton in 1671 and, following the translation of Fourier's book, the concept and uses of spectra was extensively discussed by Kelvin and Stokes. This work was all for stationary processes, that is, cases where the statistics did not evolve with time. Spectra for non--stationary process were published by Loève in 1945 and by Koenig et al. in 1946. Koenig's "spectrograph'' was used immediately and extensively in applications such as underwater sound where changes in the spectrum should be detected rapidly. In contrast, Loève two-dimensional spectra had been described in terms of expected values and, prior to multitapers, was almost impossible to estimate as a general tool. However, because there are almost no truly stationary processes, Loève's ideas should be more generally used. Turning to multitaper, my 1982 paper showed an example of a Loève spectrum in Figure 30, but despite being computable, these spectra remained difficult to interpret. In this talk I discuss some ideas for making these estimates into useful analysis tools.

Time series regression is perhaps the most obvious approach to expressing a response time series as a function of one or more predictor series. While a linear filter model may feature less easily interpretable coefficients, it's often better suited to incorporating potential temporal trends inherent to the data - trends the traditional regression model would ignore in assuming its residuals are uncorrelated. The filter model can be represented by a complex regression on the frequency-domain counterparts of the original series, offering insight into what coherency may be present in the model. Within the multitaper framework, these objects take the form of eigencoefficients related via a transfer function of frequency. Multitaper Transfer Function Estimates (MTFEs) offer an alluring topic of research in that their distributional behaviours offer promising signal detection techniques. Under certain conditions, the MTFE's modulus and variance respectively provide test statistics for detecting signals in the response and predictor, with the former (T1) demonstrating less sensitivity to false positives than the MSC, and the latter (T2) performing more robustly to frequency modulation than the Harmonic F-test. The MTFE isn't frequency stationary; this diminishes the effectiveness of T2 when restricted to a single realization. The MTFE's phase distribution, however, may be able to expose the same information as T2 without succumbing to this phenomenon, and is specifically sensitive to signals coherent between the predictor and response. This presentation hence explores the MTFE's phase as a signal detection tool for models whose underlying noise is assumed to be Gaussian and weakly stationary.

Chave recently proposed an estimator for multitaper spectral density where the time series contains missing values. In this article, we generalize this technique to a multitaper estimator of coherence and phase and show that one can also obtain bootstrapped confidence intervals. We give three examples, the first of which is a toy example in which the true coherence is known. In the second example we show that the multitaper missing-data coherence estimator computed on real data with a single gap comprising 11% of the data outperforms the Daniell-smoothed coherence estimator where there are no gaps. Finally, vertical wind speed time series collected by a Doppler lidar is punctuated by gaps at 15-minute intervals, which makes it impossible to observe convective boundary layer processes which typically evolve on time scales less than one hour. Spectrum and coherence estimation using missing-data methods shows boundary layer meteorological features at the time scales of interest.