<jupyter_start><jupyter_text>Oscillation PresenceTopic: neural oscillations are variably present, and therefore have to be explicitly detected. IssueNeural activity contains both periodic (oscillatory) and aperiodic activity, which contributes power across all frequencies. Neural oscillations (periodic activity) are variable in their occurrence. Due to aperiodic activity, there will always be power at any given frequency.  SolutionInvestigations of neural oscillations should include a _detection_ step, in which the presence of oscillatory activity is explicitly validated. This can be done be with spectral parameterization, as the presence of a frequency specific oscillation can be detected as a peak in the power spectrum.  Related WorkThe following papers include discussion of this topic:- [Are EEGs mainly rhythmic? Assessment of periodicity in wide-band time series](https://doi.org/10.1016/S0306-4522(03)00208-2), Bullock et al, 2003- [Scaling Brain Size, Keeping Timing: Evolutionary Preservation of Brain Rhythms](https://doi.org/10.1016/j.neuron.2013.10.002), Buzsaki, Logothetis, Singer, 2013- [Parameterizing Neural Power Spectra into Periodic and Aperiodic Components](https://doi.org/10.1038/s41593-020-00744-x), Donoghue et al, 2020 Related Discussion / ExamplesThe following examples discuss related topics:- A discusion on[signal rhythmicity](https://fooof-tools.github.io/fooof/auto_motivations/concepts/plot_DoYouEvenOscillate.html),and representation of different types of signals in the frequency domain. - A discussion on [filtering aperiodic activity](https://fooof-tools.github.io/fooof/auto_motivations/concepts/plot_IfYouFilterTheyWillCome.html), which includes demonstrations of filtering aperiodic activity.<jupyter_code>import numpy as np
import matplotlib.pyplot as plt

from fooof import FOOOF
from fooof.bands import Bands
from fooof.utils import trim_spectrum
from fooof.plts.annotate import plot_annotated_model
from fooof.plts.spectra import plot_spectrum, plot_spectrum_shading

from neurodsp.filt import filter_signal
from neurodsp.spectral import compute_spectrum
from neurodsp.sim import sim_oscillation, sim_powerlaw
from neurodsp.plts.time_series import plot_time_series
from neurodsp.utils import create_times, set_random_seed
# Add path and import custom code for the project
import sys
sys.path.append('code')
from plts import style_psd, savefig
from settings import BANDS, BAND_COLORS<jupyter_output><empty_output><jupyter_text>Settings<jupyter_code>import seaborn as sns
sns.set_context('poster')
# Set random seed
set_random_seed(808)
# Define general simulation settings
n_seconds = 20
fs = 1000

times = create_times(n_seconds, fs)
# Define parameters for the simulations
exp = -1.5
cf = 10
ap_filt = (1.5, 75)
# Set frequency ranges of interest
plt_range = [2, 30]
filt_xlim = [5, 8]
# Plot settings
plt_kwargs = {'xlabel' : '', 'ylabel' : '', 'lw' : 3}
shade_colors = [BAND_COLORS[band] for band in BANDS.labels]
# Set whether to save out figures
SAVE_FIG = False<jupyter_output><empty_output><jupyter_text>Aperiodic Signal Have Power Across all FrequenciesFor an initial example, let's consider the dirac delta - a single value of 1 in a series of zeros.<jupyter_code># Simulate a delta function
n_points = len(times)
dirac_sig = np.zeros([n_points])
dirac_sig[round(n_points/2)] = 1
# Plot the time series of the delta signal
plot_time_series(times, dirac_sig, lw=4)
plt.axis('off');
savefig(SAVE_FIG, '01-dirac_ts')
# Compute a power spectrum of the Dirac delta
freqs, powers = compute_spectrum(dirac_sig, fs)
plot_spectrum(*trim_spectrum(freqs, powers, [1.5, 30]), 
              color='black', lw=6, log_powers=False)

plt.xlim(plt_range)
style_psd(plt.gca())
savefig(SAVE_FIG, '01-dirac_psd')<jupyter_output><empty_output><jupyter_text>Rhythmic signals have power at specific frequenciesNext, lets examine a pure sinusoid.<jupyter_code># Simulate a sinusoid
osc_sig = sim_oscillation(n_seconds, fs, 10)
# Plot the time series of the sinusoidal signal
plot_time_series(times, osc_sig, lw=3, xlim=[2, 3])
plt.axis('off');
savefig(SAVE_FIG, '01-sine-ts')
# Compute a power spectrum of the Dirac delta
freqs, powers = compute_spectrum(osc_sig, fs, nperseg=2*fs)
plot_spectrum(*trim_spectrum(freqs, powers, [1, 50]), 
              color='black', lw=6, log_powers=False)

plt.xlim(plt_range)
style_psd(plt.gca())

savefig(SAVE_FIG, '01-sine-psd')<jupyter_output><empty_output><jupyter_text>Simulate Aperiodic Time SeriesFirst, let's simulate an example time series. For this example, we want to emphasize that measures may not necessarily reflect oscillatory activity. To do so, we will simulate a signal with neurally plausible aperiodic activity, but with _no_ oscillations. In this scenario, we can examine measurements that reflect only the aperiodic activity.<jupyter_code># Simulate aperiodic time series
ap_sig = sim_powerlaw(n_seconds, fs, exponent=exp, f_range=ap_filt)
# Plot the simulated time series
_, ax = plt.subplots(figsize=(7, 2))
plot_time_series(times, ap_sig, xlim=(10, 13), alpha=0.75, ax=ax)
plt.axis('off')
savefig(SAVE_FIG, '01-ts')
# Zoom in and plot a segment of the simulated signal
st, en = 10.2, 11
_, ax = plt.subplots(figsize=(8, 2.5))
plot_time_series(times, ap_sig, ax=ax, xlim=[st, en], **plt_kwargs)
ax.set_xticks([]); ax.set_yticks([]);
savefig(SAVE_FIG, '01-ts-segment')<jupyter_output><empty_output><jupyter_text>Compute Power SpectrumNext, we can compute the power spectrum of our measured signal. Since we simulated a signal with aperiodic activity, note that we expect the power spectrum to reflect this 1/f property of the data.<jupyter_code># Compute power spectrum of the simulated signal
freqs, powers = trim_spectrum(*compute_spectrum(ap_sig, fs, nperseg=fs), plt_range)
# Plot the power spectrum of the simulated signal
plot_spectrum_shading(freqs, powers, lw=5, color='black', log_powers=True,
                      shades=BANDS.definitions, shade_colors=shade_colors)
plt.xlim(plt_range)
style_psd(plt.gca())
savefig(SAVE_FIG, '01-ap_psd_col')<jupyter_output><empty_output><jupyter_text>In the spectrum above we do see the expected 1/f nature of the power spectrum. Notably, due to the 1/f there is power across all frequencies, despite there being no _oscillatory_ power.   Examine Filtered RangesNext, let's examine what the data looks like when filtered in band-specific frequency ranges.<jupyter_code># Create a plot of the data filtered into different frequency ranges
_, axes = plt.subplots(BANDS.n_bands, 1, figsize=(8, 9))
for ax, color, (label, f_range) in zip(axes, shade_colors, BANDS):

    band_sig = filter_signal(ap_sig, fs, 'bandpass', f_range)
    plot_time_series(times, band_sig, ax=ax,
                     xlim=filt_xlim, ylim=(-1.25, 1.25), **plt_kwargs)
    ax.axvspan(filt_xlim[0], filt_xlim[1], alpha=0.2, color=color)
    ax.set_xticks([]); ax.set_yticks([]); 
    
plt.subplots_adjust(hspace=0.05)
savefig(SAVE_FIG, '01-ap_ts_bands')<jupyter_output><empty_output><jupyter_text>Note in the above filtered data, that there is activity in all bands, and it even seems to have interesting dynamics. This is because the aperiodic activity has activity at all frequencies, that can be extracted with filters. In addition, the seeming sinusoidality of the filtered traces stems from the filters themselves. However, these traces do not reflect true oscillatory activity, and interpreting it as such would be wrong.   Compare to a signal with an oscillation<jupyter_code># Simulate a time series with an oscillation
osc_sig = sim_oscillation(n_seconds, fs, freq=cf, variance=0.5)

# Combine the oscillation with the aperiodic signal from before
sig_comb = ap_sig + osc_sig
# Plot the simulated time series
plot_time_series(times, sig_comb, xlim=(10, 13), lw=2.5, alpha=0.75)
plt.axis('off')
savefig(SAVE_FIG, '01-comb_ts')
# Compute power spectrum of the simulated signal
freqs, powers = trim_spectrum(*compute_spectrum(sig_comb, fs, nperseg=1.5*fs), plt_range)
# Plot the power spectrum of the simulated signal
plot_spectrum_shading(freqs, powers, lw=5, color='black', log_powers=True,
                      shades=BANDS.definitions, shade_colors=shade_colors)
plt.xlim(plt_range)
style_psd(plt.gca())
savefig(SAVE_FIG, '01-comb_psd_col')<jupyter_output><empty_output><jupyter_text>Peak detectionIn order to detect an oscillation, we can examine the the power spectrum for a peak of frequency-specific activity. Here, we will use [spectral parameterization](https://fooof-tools.github.io/)to detect peaks in the power spectrum.<jupyter_code># Initialize an object for spectral parameterization
fm = FOOOF(verbose=False)

# Fit the spectral parameterization model
fm.fit(freqs, powers)
# Set the plot color for the shaded peak to be the alpha color
from fooof.plts.settings import PLT_COLORS
PLT_COLORS['periodic'] = shade_colors[2]
# Plot the power spectrum, with the model, highlighting the detected peak
plot_annotated_model(fm, annotate_peaks=False, annotate_aperiodic=False)
plt.gca().get_legend().remove()
plt.xlim(plt_range)
style_psd(plt.gca())
savefig(SAVE_FIG, '01-comb_psd_peak')<jupyter_output><empty_output><jupyter_text>Having detected a peak, we can restrict our analyses to frequency bands with verified oscillations.<jupyter_code># Plot the filtered range for the band with a detected peak
_, ax = plt.subplots(figsize=(8, 3))
band_sig = filter_signal(sig_comb, fs, 'bandpass', BANDS.alpha)
plot_time_series(times, band_sig, ax=ax,
                 xlim=filt_xlim, ylim=(-1.5, 1.5), **plt_kwargs)
ax.set_xticks([]); ax.set_yticks([]); 
ax.axvspan(filt_xlim[0], filt_xlim[1], alpha=0.2, color=shade_colors[2])
savefig(SAVE_FIG, '01-alpha_ts_bands')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09360213947747377, 'fraction_numerical': 0.0202633203044641, 'mean_word_length': 6.005043227665706, 'pattern_counts': {'":': 0, '<': 41, '<?xml version=': 0, '>': 41, 'https://': 6, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '280323', 'license': 'MIT', 'n_tokens_mistral': 3337, 'n_tokens_neox': 3075, 'n_words': 990, 'path': '01-OscillationPresence.ipynb', 'repo_name': 'OscillationMethods/OscillationMethods', 'size': '9722'}