[ENH] - Simulation of Central Frequency in the Time Domain

[elybrand]

Overview

This PR builds a new combined aperiodic and periodic simulation which allows the user to simulate an aperiodic component with a specified power law exponent and a single oscillatory component with specified central frequency, bandwidth, and relative height above the aperiodic component in the log-power spectrum.

Here is a minimal working example with the analytical power spectrum superimposed for comparison.

import numpy as np
import matplotlib.pyplot as plt
from neurodsp.utils.data import create_times
from neurodsp.plts import plot_power_spectra, plot_time_series
from neurodsp.sim import sim_central_freq
from neurodsp.spectral import compute_spectrum

np.random.seed(0)

n_seconds = 10
fs = 10**3
chi = -2
central_freq = 20
bw = 5
ht = 2
sig = sim_central_freq(n_seconds, fs, chi, central_freq, bw, ht)
freqs, psd = compute_spectrum(sig, fs, )#nperseg=fs*n_seconds)
times = create_times(n_seconds, fs)

# Build a reference power law to visually track if the psd has exponent chi.
psd_baseline = psd[1]/(freqs[1]**chi)*np.array([10**(ht * np.exp(-(freq - central_freq)**2/(2*bw**2))) * freq**(chi) for freq in freqs])

fig, axes = plt.subplots(2,1, figsize=(15,10))
axes[0].set_title("Aperiodic Noise with Central Frequency")
plot_time_series(times, sig, ax=axes[0])
plot_power_spectra(freqs, [psd, psd_baseline], ax=axes[1], labels=["Empirical PSD", "Reference PSD"])
plt.tight_layout()

Comments

The way the code is written is a bit limiting. For example, it doesn't allow the presence of multiple central frequencies. The code is easily adaptable to handle this case though. In general, given a signal in the time domain, this code adds sinusoids at the relevant frequencies to add a ``bump'' in the log power spectrum to the inputted signal. To get multiple central frequencies, you'd just do a few recursive calls on aperiodic noise with different central frequencies, bandwidths, and relative heights.

@ryanhammonds @TomDonoghue Do you think I should rewrite this code so that it takes as input a general time series and then adds sinusoids to augment its power spectrum, or do you think that should be a separate function and we can keep this sim as is?

The Math Behind The Code

[ryanhammonds]

This is really cool! I fit the psd from your example with fooof to compare chi, central_freq, bw, and ht. Everything matched up as expected, except for a sign flip in the chi (negative here, but fooof reports a positive). Anyways, I'll leave a proper review for this soon.

[TomDonoghue]

This looks awesome Eric - I'm looking forward to trying it out when I get a chance!

The different sign of 1/f is expected based on different conventions between FOOOF & NDSP. Sorta annoying, but sorta grand-fathered in at this point, and consistent with everything else.

[elybrand]

@rdgao just wanted to tag you in to let you look at the math if you're so inclined.

[rdgao]

holy moly that's a lot of math

[ryanhammonds]

Do you think I should rewrite this code so that it takes as input a general time series and then adds sinusoids to augment its power spectrum, or do you think that should be a separate function and we can keep this sim as is?

I pushed an update so any aperiodic function may be used. However, it defaults to sim_powerlaw with an exponent of -2. I also fixed a doc typo, trimmed line lengths for pep8, etc.

[ryanhammonds]

@TomDonoghue I addressed all of your comments, renamed the func/variables, added array support, ran pylint to resolve multi-line white space issues. I also added an accuracy test to make sure a peak was close to the center freq: np.argmax(spectrum-spectrum_ap).