Plot Spectrogram

Example of how to generate a spectrogram using bmtool¶

By Gregory Glickert

First we read in our data and process it into an LFP

In [1]:

from bmtool.analysis.lfp import load_ecp_to_xarray, cwt_spectrogram_xarray, fit_fooof, get_windowed_data, ecp_to_lfp, plot_spectrogram
import numpy as np
import matplotlib.pyplot as plt

ecp = load_ecp_to_xarray('../network_to_analysis/long/ecp.h5')
lfp = ecp_to_lfp(ecp)

from bmtool.analysis.lfp import load_ecp_to_xarray, cwt_spectrogram_xarray, fit_fooof, get_windowed_data, ecp_to_lfp, plot_spectrogram import numpy as np import matplotlib.pyplot as plt ecp = load_ecp_to_xarray('../network_to_analysis/long/ecp.h5') lfp = ecp_to_lfp(ecp)

calculate spectrogram¶

Then we need to calculate the spectrogram. This will just make an array of the data and then we will process it further to make our plot

In [2]:

downsample_freq = 200 # if you want to further down sample
tseg = 0.5  # time segment length for PSD
axis = lfp.dims.index('time') # dimension index of time axis in x

lfp_sxx = cwt_spectrogram_xarray(lfp, lfp.fs, axis=axis, downsample_fs = downsample_freq,
                                        channel_coords={'channel_id': lfp.channel_id}, freq_range=(1 / tseg, np.inf))

downsample_freq = 200 # if you want to further down sample tseg = 0.5 # time segment length for PSD axis = lfp.dims.index('time') # dimension index of time axis in x lfp_sxx = cwt_spectrogram_xarray(lfp, lfp.fs, axis=axis, downsample_fs = downsample_freq, channel_coords={'channel_id': lfp.channel_id}, freq_range=(1 / tseg, np.inf))

Plot spectrogram¶

Whole simulation spectrogram with the aperiodic removed can be done like this. It may be good to check the aperiodic fit of the model before completely trusting the plot. This can be done by setting the plot argument in the fit_fooof function to True.

In [3]:

fooof_params = dict(aperiodic_mode='fixed', freq_range=(1,100), peak_width_limits=100.)
plt_range = [2., 100.]
clr_freq_range = None
log_power = 'dB'

fig, ax = plt.subplots(1,1,figsize=(10,6))
sxx = lfp_sxx.sel(channel_id=0)
sxx_tot = sxx.PSD.mean(dim='time')

fooof_results, _ = fit_fooof(sxx_tot.frequency.values, sxx_tot.values, **fooof_params,plot=False)
_ = plot_spectrogram(sxx, remove_aperiodic=fooof_results, log_power=log_power,
                        plt_range=plt_range, clr_freq_range=clr_freq_range, pad=0.01, ax=ax)

fooof_params = dict(aperiodic_mode='fixed', freq_range=(1,100), peak_width_limits=100.) plt_range = [2., 100.] clr_freq_range = None log_power = 'dB' fig, ax = plt.subplots(1,1,figsize=(10,6)) sxx = lfp_sxx.sel(channel_id=0) sxx_tot = sxx.PSD.mean(dim='time') fooof_results, _ = fit_fooof(sxx_tot.frequency.values, sxx_tot.values, **fooof_params,plot=False) _ = plot_spectrogram(sxx, remove_aperiodic=fooof_results, log_power=log_power, plt_range=plt_range, clr_freq_range=clr_freq_range, pad=0.01, ax=ax)

Average spectrogram and plot¶

For this simulation be analyzed there were 9 pulses given. It may be easier to see the pulse response if we average the times together

In [4]:

# in this example there are 9 pulses given for 1 second with a 500ms break before the next pulse. The pulse starts 1 second into the simulation 
on_time = 1 # how long in seconds is our stimulus
off_time = 0.5 # how long is the break between stim is 
t_start = 1 # when starts in sec
stimulus_count = 9

times_to_avg = []
for i in range(stimulus_count):
    end_time = t_start+on_time+off_time
    times_to_avg.append([t_start,end_time])
    t_start = end_time 

# make list into np array 
times_to_avg = np.array(times_to_avg)

_, _, lfp_sxx_avg = get_windowed_data(lfp_sxx.PSD, times_to_avg, {0: np.arange(times_to_avg.shape[0])})
lfp_sxx_avg = lfp_sxx_avg[0].mean_.sel(unique_cycle=0).to_dataset(name='PSD')

fig, ax = plt.subplots(1,1,figsize=(10,6))
sxx = lfp_sxx_avg.sel(channel_id=0)
sxx_tot = sxx.PSD.mean(dim='time')
fooof_results, _ = fit_fooof(sxx_tot.frequency.values, sxx_tot.values, **fooof_params,plot=False)
_ = plot_spectrogram(sxx, remove_aperiodic=fooof_results, log_power=log_power,
                        plt_range=plt_range, clr_freq_range=clr_freq_range, pad=0.01, ax=ax)

# in this example there are 9 pulses given for 1 second with a 500ms break before the next pulse. The pulse starts 1 second into the simulation on_time = 1 # how long in seconds is our stimulus off_time = 0.5 # how long is the break between stim is t_start = 1 # when starts in sec stimulus_count = 9 times_to_avg = [] for i in range(stimulus_count): end_time = t_start+on_time+off_time times_to_avg.append([t_start,end_time]) t_start = end_time # make list into np array times_to_avg = np.array(times_to_avg) _, _, lfp_sxx_avg = get_windowed_data(lfp_sxx.PSD, times_to_avg, {0: np.arange(times_to_avg.shape[0])}) lfp_sxx_avg = lfp_sxx_avg[0].mean_.sel(unique_cycle=0).to_dataset(name='PSD') fig, ax = plt.subplots(1,1,figsize=(10,6)) sxx = lfp_sxx_avg.sel(channel_id=0) sxx_tot = sxx.PSD.mean(dim='time') fooof_results, _ = fit_fooof(sxx_tot.frequency.values, sxx_tot.values, **fooof_params,plot=False) _ = plot_spectrogram(sxx, remove_aperiodic=fooof_results, log_power=log_power, plt_range=plt_range, clr_freq_range=clr_freq_range, pad=0.01, ax=ax)

/home/gjgpb9/miniconda3/envs/bmtk/lib/python3.8/site-packages/numpy/lib/nanfunctions.py:1879: RuntimeWarning: Degrees of freedom <= 0 for slice.
  var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof,