memspectrum

Core module of memspectrum. Implements class MESA, which computes the spectrum of a given time-series using the Burg algorithm, plus some helpers.

class MESA(filename=None, *args, **kwargs)

Class the implement the reproduces the Maximum Entropy Spectrum of a given time-series.

compute_autocorrelation(dt, normalize=True, scipy_convention=False)

Compute the autocovariance \(C(\tau)\) of the data based on the autoregressive coefficients. The autocovariance is defined as:

\[C(\tau) = E_t[(x_t - \mu)(x_{t+\tau} -\mu)]\]

where \(\mu\) is the mean value of the timeseries. It amounts to the inverse Fourier transform of the PSD.

If option scipy_convention is set, we center the autocorrelation, following correlate() conventions.

\[C(\tau) = E_t[(x_t - \mu)(x_{t-\tau+T} -\mu)]\]

where \(T\) is the segment length

Parameters:

• dt (float) – Sampling interval

• normalize (bool) – Whether the autocovariance should be normalized s.t. \(C(\tau = 0) = 1\)

Returns:

autocorr – Autocorrelation of the model

Return type:

ndarray

entropy_rate(dt)

Compute the entropy gain \(\Delta H\) for a given power spectrum:

\[\Delta H = \int_{- Ny}^{Ny}\log S(f) df\]

It is described in Eq (84) in Papoulis or Eq (10) in Martini et al.

Parameters:

dt (float) – Sampling interval

Returns:

rate – Entropy rate

Return type:

float

forecast(data, length, number_of_simulations=1, P=None, include_data=False, seed=None, verbose=False)

Forecasting on an observed process for a total number of points given by length. It computes number_of_simulations realization of the forecast time series. This method can only be used if a_k coefficients have been computed already. Use solve method before forecasting.

Parameters:

• data (ndarray) – data for the spectrum calculation

• length ('np.int') – Number of future points to be predicted

• number_of_simulations ('np.int') – Total number of simulations of the process

• P (float) – Variance of white noise for the autoregressive process. Default is None and uses the estimate obtained with Burg’s algorithm.

• include_data (bool) – Whether to prepend to the output the input time series

• seed (int) – Seed for the random generator. If is None, no initialization of the seed is done and the authomatic numpy setting is used.

• verbose (bool) – Whether to print the status of the forecasting

Returns:

predictions – Array containing the forecasted points for every simulation of the process (Shape (number_of_simulations, length))

Return type:

ndarray

generate_noise(T, sampling_rate=1.0, fmin=None, fmax=None, N_series=1)

Generate some noise starting from the internal \(AR(p)\) model. The noise generated has the same features as the data given in input to the mesa.

Parameters:

• T (float) – Length (in seconds) of the signal to generate

• sampling_rate (float) – Sampling rate

• fmin (float) – Minimum frequency in the signal (if None, is equal to zero)

• fmax (float) – Maximum frequency in the signal (if None, Nyquist frequency is used: f_Ny = 0.5*sampling_rate)

• N_series (float) – Number of time series to generate

Returns:

• times (ndarray) – Time grid at which the noise time series is evaluated at

• times_series (ndarray) – Time series (shape (N_series, sampling_rate*T) )

• frequencies (ndarray) – Frequency grid at which the noise frequency series is evaluated at

• frequency_series (ndarray) – Frequency series from which the times_series is computed (shape (N_series, sampling_rate*T) )

• psd (ndarray) – PSD of the model, including both positive and negative frequencies

get_p()

Returns the order of the autoregressive process that defines the PSD.

Returns:

p – Order of the autoregressive process that define the PSD

Return type:

int

load(filename)

Load the class from a given file. The file shall be the same format produced by save().

Parameters:

filename (str) – Name of the file to load the data from

logL(data, dt)

Compute the log likelihood given the current spectrum data smust be in the time domain.

Parameters:

• data (ndarray) – Stretch of data to compute the likelihood of

• dt (float) – Sampling interval

Returns:

white_data – Stretch of (trimmed) whitened data

Return type:

ndarray

save(filename)

Save the class to file (if the spectral density analysis is already performed). The output file can be used to load the class with method load() File is a 1D array with the format: [P, N, mu, a_k]. The header holds the shapes for each array

Parameters:

filename (str) – Name of the file to save the data at

save_spectrum(filename, dt, frequencies=None)

Saves the power spectral density computed by the model to a txt file. The PSD is evaluated on a user given grid of frequencies. If None, a standard grid is used (as computed by np.fft.fftfreq). The spectrum is saved as a 2D array: [f, PSD(f)]

Parameters:

• filename (str) – Name of the file to save the PSD at

• dt (float) – Sampling rate for the time series

• frequencies (ndarray) – Frequencies to compute the PSD at. If None, a standard grid will be computed

solve(data, m=None, optimisation_method='FPE', method='Fast', regularisation=1.0, early_stop=True, verbose=False)

Computes the power spectral density of the given data data by solving the Levinson recursion for the autoregressive coefficients a_k and the noise variance P. It uses the standard Burg method recursive and a Faster (but less stable) version. Default is ‘Fast’.

Parameters:

• data (ndarray) – One dimensional array with the time series data for the spectrum calculation

• m (int) – Maximum number of recursions for the computation of the power spectral density. Maximum autoregressive order is p = m-1 Default is None, that means m = 2N / log(2N)

• optimisation_method ('str') – Method used to select the best recursive order. The order is chosen minimizing the corresponding method. Available methods are “FPE”, “OBD”, “CAT”, “AIC”, “Fixed”. If optimisation_method is “Fixed”, the autoregressive order is always set to m, without looking for a minimum. Deafult is “FPE”.

• method ('str') – Can be “standard” or “Fast”. Selects the algorithm used to compute the power spectral density with. Default is “Fast”

• regularisation (float) – Implement Tikhonov regularisation. Should be a number slightly larger than one. Default is 1, which means no regularisation

• early_stop (bool) – Default is True. Breaks the iteration if there is no new global maximum after 100 iterations. Recommended for every optimisation method but CAT. Has no effect with “Fixed” loss function.

• verbose (bool) – Whether to print the status of the mesa solution

Returns:

• P (float) – Variance of white noise for the associated autoregressive process

• a_k (ndarray) – The coefficient used to compute the power spectral density

• optimization (ndarray) – The values of the chosen optimisation_method at every iteration

spectrum(dt=1.0, frequencies=None, onesided=False)

Computes the power spectral density of the model. Default returns power spectral density and frequency array automatically computed by sampling theory. It can also be computed on a user-given frequency grid passing proper frequency array

Parameters:

• dt (float) – Sampling interval

• frequencies (ndarray) – (positive) frequencies to evaluate the spectrum at If None, a equally spaced frequency grid is used (and returned)

• onesided (bool) – Whether the one sided PSD (only positive frequencies) shall be returned. It has effect only if a frequency array is not given

Returns:

• frequencies (ndarray) – Frequencies at which power spectral density is evaluated (only returned if frequencies = None)

• spectrum (ndarray) – Power spectral density

whiten(data, trim=None, zero_phase=False)

Whiten the data by convolving with the autoregressive coeffiecients a_k

Parameters:

• data (ndarray) – Stretch of data to whiten

• trim (int) – Number of points to remove at the begining and at the end of the stretch of data to remove edge effects. If None, trim is set to be equal to the autoregressive order.

Returns:

white_data – Stretch of (trimmed) whitened data

Return type:

ndarray

GenerateTimeSeries