Multitaper.jl: A Julia package for frequency domain analysis of time series

Spectral analysis is widely used to infer physical mechanisms for underlying process dynamics from a realization of a stationary time series. The multitaper method is a nonparametric technique for estimating the power spectrum of a discrete time series that simultaneously controls bias and variance of the estimator by premultiplying the data by a set of orthogonal sequences—discrete prolate spheroidal sequences that are optimally concentrated in both time and frequency. These have the effect of stabilizing the bias and variance of the spectrum estimator (Thomson, 1982). While multitaper codes have been introduced in multiple languages, including Julia, the Multitaper.jl package that we present offers functionality beyond univariate and bivariate time series analysis and provides routines for a number of (selected) research-level topics not found in other packages.


Summary
Spectral analysis is widely used to infer physical mechanisms for underlying process dynamics from a realization of a stationary time series. The multitaper method is a nonparametric technique for estimating the power spectrum of a discrete time series that simultaneously controls bias and variance of the estimator by premultiplying the data by a set of orthogonal sequences-discrete prolate spheroidal sequences that are optimally concentrated in both time and frequency. These have the effect of stabilizing the bias and variance of the spectrum estimator (Thomson, 1982). While multitaper codes have been introduced in multiple languages, including Julia, the Multitaper.jl package that we present offers functionality beyond univariate and bivariate time series analysis and provides routines for a number of (selected) research-level topics not found in other packages.

Statement of need
Multitaper.jl is a Julia package for spectrum analysis of time series, multivariate time series, and spatial or space-time processes. Multitaper.jl was designed to be useful to researchers in diverse fields, including geophysics (climate, seismology, limnology, and stratigraphy), cognitive radio, space science (solar physics), speech processing, astronomy, and biomedicine. For example, a researcher might want to compute the multitaper spectrum of a time series so he or she can identify which periodic components contribute the most to signal variance, and do so with jackknifed error bounds on the oscillation amplitudes.
The high-level character of Julia allows for widely readable and extendible codes, while the low-level functionality provides speed and efficiency. The Multitaper.jl package provides a user-friendly implementation of many of the basic concepts such as spectrum analysis, Ftesting for harmonic analysis (Thomson, 1982), coherence and phase, jackknifed variance estimates (Thomson & Chave, 1991), and complex demodulation (Thomson et al., 2007); more advanced techniques such as dual-frequency spectra, cepstrum, multitaper for time series containing gaps (Chave, 2019), T 2 tests for multiple line components (Thomson, 2011) implementations of higher-dimensional Slepian tapers on Cartesian domains (Simons & Wang, 2011) (Geoga et al., 2018); and others (Thomson & Haley, 2014) (Haley & Anitescu, 2017). In addition, we provide tutorial-style notebooks to allow accessibility to those new to these concepts or to Julia in general.
Multitaper.jl has been used in graduate courses to provide fast spectrum estimates of unequally sampled time series. Early versions of this code have also been used to compute statistics for research publications.

Other software
Multitaper spectrum analysis is implemented in Julia in the DSP.jl package, but it is limited to estimation of the spectrum. In the R programming language the R multitaper package gives an R-wrapped Fortran 77 implementation, which provides fast spectrum estimates, F-tests, jackknifing, coherences, and complex demdodulation (Rahim & Burr, 2013). A C-subroutine for multitaper methods was introduced in (Lees & Park, 1995), and the derivative Python implementation pymutt (Smith, 2011) is a wrapped version of the former. The Fortran 90 library of multitaper methods (Prieto et al., 2009) provides for spectrum analysis, F-testing, spectral reshaping, coherences, and quadratic inverse spectrum estimation. While wrapped versions of low-level codes run rapidly, they can be more difficult to extend by the research community. The Matlab Signal Processing Toolbox implements discrete prolate spheroidal sequences and multitaper spectrum estimation, while standalone contributions on the Matlab file exchange describe the extension to multitaper methods with gaps (Chave, 2019), higherdimensional spectrum estimation on Cartesian domains (Simons & Wang, 2011) and the sphere (Simons et al., 2006) and the freely available jlab codes (Lilly, 2019).

To contribute
We welcome input of any kind via bitbucket issues or by pull requests. Support requests can be directed to haley@anl.gov.