memspectrum 1.1.0

maximum entropy spectral analysis

memspectrum

Authors Alessandro Martini, Stefano Schmidt, Walter del Pozzo

email [e-mail]

Copyright Copyright (C) 2020 Alessandro Martini

Licence CC BY 4.0

Version 1.1.0

MAXIMUM ENTROPY ESTIMATION ALGORITHM FOR FAST PSD COMPUTATION

memspectrum is a package for the computation of power spectral densities of a time series. It implements a Fast verions of Burg method of Maximum Entropy Spectral Analysis. The method is fast and reliable and provides better performance than other standard methods.

The computation of the power spectral density requires solving the Levinson recursion for the forward prediction error coefficients a_k. The knowledge of such coefficients allows to characterize the observed process in terms of an autoregressive process of order p (AR(p)), being p + 1 the lenght of the a_k array. Together with a_k coefficients a P coefficient is estimated, and is to be interpreted as the variance of white noise component for the process. The computation of these quantities allow to perform high quality forecast for the time series. The estimate of the autoregressive order is via a second class, that implements several methods available in literature.

Usage of memspectrum

To get the PSD computed, the following steps are required

• Import the data

• Call MESA class passing data as argument

from memspectrum import MESA
m = MESA()

• Compute the coefficients via the solve() method: MANDATORY for further computations

m.solve(data)

• At this point you can compute the spectrum and forecast

m.spectrum(dt)
m.forecast(data)

Sinusoid example

To compute (and plot) the spectrum of a (noisy) sinusoidal signal:

from memspectrum import MESA
import numpy as np

Generating the data:

N, dt = 1000, .01  #Number of samples and sampling interval
time = np.arange(0, N) * dt
frequency = 2
data = np.sin(2 * np.pi * frequency * time) + np.random.normal(.4, size = 1000)
plt.plot(time, data, color = 'k')

Solving MESA is needed to compute PSD or forecast.

M = MESA()
M.solve(data)

The spectrum can be computed on sampling frequencies (automatically generated) or on some given interval

spectrum, frequencies = M.spectrum(dt)  #Computes on sampling frequencies
user_frequencies = np.linspace(1.5, 2.5)
user_spectrum = M.spectrum(dt, user_frequencies) #Computes on desired window

Plotting the two the following is obtained:

It can also be used to perform forecasting. For example, we consider the first 900 points of the data and try to infer the upcoming signal. 1000 simulations of 100 points are performed. Real observed data are compared with median estimate and 90% Credibility regions

M = MESA()
M.solve(data[:-100])
forecast = M.forecast(data[:-100], length = 100, number_of_simulations = 1000, include_data = False)
median = np.median(forecast, axis = 0) #Ensemble median
p5, p95 = np.percentile(forecast, (5, 95), axis = 0) #90% credibility boundaries

plt.plot(time[:-100], data[:-100], color = 'k')
plt.fill_between(time[-100:], p5, p95, color = 'b', alpha = .5, label = '90% Cr.')
plt.plot(time[-100:], data[-100:], color = 'k', linestyle = '-.', label = 'Observed data')
plt.plot(time[-100:], median, color = 'r', label = 'median estimate')

The forecast result is:

Generating data from PSD

memspectrum.generateTimeSeries provides a function that construct a time-series with a user-given power <<<<<<< HEAD spectral density. It can be called as :: ======= spectral density. It can be called as :: >>>>>>> af7984f4cf9ef4ad7b528e26a5fd54e6ccbfb4d4

from memspectrum.generateTimeSerie import generate_data f, psd = import wanted psd and frequency array time, time_series, frequency, frequency_series, psd = generate_data(f, psd, T, sampling_rate)

T represent the time length of the observation and sampling rate is equivalent to 1 / dt, with dt the sampling interval

Installation & documentation

To install the package:

pip install memspectrum

It requires numpy.

On the GitHub repository, a number of examples are available to the interested user:

• gwstrain.py: computes the PSD on a piece of gravitational waves data and perform some forecasting

• sunspots.py: using data from sunspots, it uses memspectrum to find an autoregressive process which describes them and forecast

• sound_MESA.py: given an input audio (wav) file reproducing the sound of a waterfall, it computes the PSD and generate a synthetic noise, resembling the original one.

For more advanced use or for more information, please refer to the code documentation:

import memspectrum
help(memspectrum)
help(memspectrum.<function_name>)

For full source code (and much more) see: https://github.com/martini-alessandro/Maximum-Entropy-Spectrum