[ENH] PSD fitting based on AR form

[ryanhammonds]

Adds spectral fitting based on the PSD form of an AR(p) model:

$$ P(f) = \frac{\sigma^{2}}{|1 - \sum\limits_{k=1}^p \varphi_{k} e\strut^{-j 2 \pi f k}|^{2}} $$

Given a high enough AR order, a wide range spectral shapes can be fit. This gives us a simulator for free, since we can take the $\varphi$'s and a white noise input to generate a signal with the desired spectral shape. Here is a double knee example:

import numpy as np
import matplotlib.pyplot as plt
from timescales.autoreg import ARPSD, ar_spectrum

# Simulate spectral form
fs = 2000
freqs = np.linspace(.1, 1000, 2000)
k = np.arange(1, 2)
exp = np.exp(-2j * np.pi * np.outer(freqs, k) / fs).T
powers = ar_spectrum(exp, 0.99, 1) + ar_spectrum(exp, 0.7, 50.)

# Fit PSD with AR form
order = 100
ar = ARPSD(order, fs, ar_bounds=(0., 1.))
ar.fit(freqs, powers)

# Simulate timeseries
sig = ar.simulate(200, fs)

After fitting, using ar.plot():

Then we can take the signal from ar.simulate() and confirm the spectrum of the simulated signal matches the desired shape:

Oscillations also work. The fit spectrum (a higher order may give more accurate results):

And the corresponding simulated signal based on the learned ar coefficients.

There seems to be a wide range of spectral shapes that high order AR models can approximate. It turns out the form of P(f) is a universal approximator since the denominator is a polynomial + Weierstrass approximation theorem. This PR gives a simulation method for arbitrary spectral shapes.

[codecov-commenter]

Welcome to Codecov 🎉

Once you merge this PR into your default branch, you're all set! Codecov will compare coverage reports and display results in all future pull requests.

Thanks for integrating Codecov - We've got you covered ☂️