Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update acf <-> psd conversions #32

Open
wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Update acf <-> psd conversions #32

wants to merge 1 commit into from

Conversation

ryanhammonds
Copy link
Member

@ryanhammonds ryanhammonds commented Aug 27, 2024

Addresses #29 and updates the conversion functions. For ACF -> PSD:

import matplotlib.pyplot as plt
import numpy as np
from neurodsp.spectral import compute_spectrum
from statsmodels.tsa.stattools import acf
from timescales.sim import sim_exp_decay, sim_ar
from timescales.conversions import acf_to_psd, psd_to_acf

# Parameters
fs = 1000
phi = 0.95

# Target ACF
tau = -1/(np.log(phi) * fs)
lags = np.arange(5000)
corrs = sim_exp_decay(lags, fs, tau, 1.)

# Simulated sample
sig = sim_ar(1000, 1000, np.array([phi]))

# ACF to PSD
freqs, powers = acf_to_psd(corrs, fs)
freqs_welch, powers_welch = compute_spectrum(sig, fs, nperseg=5001)

plt.loglog(freqs_welch[1:], powers_welch[1:] / powers_welch[:20].mean(), label='Welch, from simulated signal')
plt.loglog(freqs[1:], powers[1:] / powers[1], label='abs(iff(acf(t)))^2, from theoretical acf')
plt.legend()

For PSD -> ACF:

# AR(1) PSD
freqs = np.linspace(0.01, fs//2, 20000)
powers = 1 / (1 - 2 * phi * np.cos(2 * np.pi * freqs * 1/fs) + phi**2)

# Convert PSD to ACF
lags_ifft, corrs_ifft = psd_to_acf(freqs, powers, fs)

plt.plot(lags[:100], corrs[:100], label='simulated acf')
plt.plot(lags_ifft[:50], corrs_ifft[:50], label='ifft(powers), from simulate AR(1) psd', ls='--')

These fft conversions agree on the PSD forms from on the AR(1) PSD (#31 has spectral fitting for AR(p)), and OU and AR(1) simulations. This deviates from the Lorentzian form in specparam. The specparam form can be derived from the AR(1) PSD using a couple approximations, e.g. taylor approximation of the cosine term above. These approximations result in bias as the knee frequency increases towards nyquist.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

1 participant