Time-resolved Frequency Analysis

Courses

Authors

Dr. Atle Rimehaug | Dr. Nicholas Del Grosso

Sometimes, it’s not enough to know that an oscillation at a specific frequency is present in your signal, you also need to know when the oscillation occurred. For that, we have to do a time-resolved frequency analysis of the signal. That means that you compute a power spectrum that shows how the power of different frequencies in your signal varies over time. In this notebook, you will learn about methods for computing time-resolved frequency spectra of signals.

Setup

Import Libraries

import pywt
import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path
import xarray as xr

Utility Functions

Run the cell below, but feel free to ignore this code — it’s some utility functions we’ll use later for convenience.

def calc_wavelet_spectrum(sig, sampling_frequency, number_of_freqs, bandwidth, center_frequency):
    '''
    Calculates the time-resolved power spectrum of a signal using Morlet wavelets

    Args:
        sig: signal that the wavelet spectrum is calculated for
        sampling_frequency: sampling frequency of signal
        number_of_freqs: number of frequencies the wavelet spectrum is calculated for
        bandwidth: the bandwidth of the Morlet wavelet
        center_frequency: the central frequency of the Morlet wavelet
    Returns:
        freqs: frequencies for which the spectral values are calculated
        spectrum_amp: spectral amplitudes
    '''
    # the bandwidth and the central frequency sets which Morlet wavelet is used
    wavelet = "cmor"+str(bandwidth)+"-"+str(center_frequency)

    # widths of the Morlet wavelet determines which frequencies the wavelet spectrum is calculated for
    widths = np.logspace(start=3, stop=10, num=number_of_freqs, base=2)

    sampling_period = 1/sampling_frequency
    spectrum, freqs = pywt.cwt(sig, widths, wavelet, sampling_period=sampling_period)

    # take absolute value to get real part and therefore amplitude of spectrum
    spectrum_amp = np.abs(spectrum)

    return freqs, spectrum_amp

def plot_spectrum_with_coi(time, freqs, spectrum, t_coi, apply_shading=True) -> None:
    '''
    Plots wavelet spectrum together with cone of influence

    Args:
        time: time points of spectrum
        freqs: frequencies the spectrum is calculated for
        spectrum: wavelet spectrum
        t_coi: time-points of cone of influence
        apply_shading: whether to apply shading to spectrum plot
    '''

    # plot wavelet spectrum
    if apply_shading:
        plt.pcolormesh(time, freqs, spectrum, cmap='plasma', shading='gouraud')
    else:
        plt.pcolormesh(time, freqs, spectrum, cmap='plasma')
    plt.ylim([0,100])

    # plot cone of influence
    plt.plot(time[0]+t_coi, np.arange(freqs.size)[::-1], 'g-', label = 'Cone of influence')
    plt.plot(time[-1]-t_coi, np.arange(freqs.size)[::-1],'g-')

    plt.xlabel('Time from stimulus onset (ms)')
    plt.ylabel('Frequency (Hz)')
    plt.title('Time-resolved frequency spectrum with\ncone of influence')
    plt.legend();

class utils:
    calc_wavelet_spectrum = calc_wavelet_spectrum
    plot_spectrum_with_coi = plot_spectrum_with_coi

Download dataset

import owncloud

owncloud.Client.from_public_link('https://uni-bonn.sciebo.de/s/rY4YCPXmkpespKG').get_file('/', 'dataset_session_754312389.nc')

True

Section 1: Time-resolved frequency analysis with Morlet wavelets

One of the most important and frequently used methods for computing time-resolved frequency spectra is to convolve the signal with Morlet wavelets. Morlet wavelets are wavelets with specific frequencies, and the process of convolving them with the signal can be seen as a matching process. You check how well the frequency of the signal matches with the frequency of the Morlet wavelet at a given time point. The better the match is, the higher the power will be for that frequency at that time point in the power spectrum. The plot below illustrates a Morlet wavelet whose frequency matches the frequency of the signal.

Code	Description
`freqs, spectrum_amp = utils.calc_wavelet_spectrum(sig, sampling_frequency)`	Calculates the wavelet spectrum `spectrum` for a signal `sig` with a given `sampling_frequency`. `number_of_freqs` sets the number of frequencies the spectrum is calculated for. `bandwidth` and `center_frequency` sets the parameters for the Morlet wavelet.
`plt.pcolormesh(x, y, C, cmap = 'plasma', shading = 'gouraud')`	Make a 2D colormap of values in a 2D array (C) against x and y values. The optional parameters `cmap` and `shading` defines the colormap and the interpolation applied to the plot, respectively.

Load dataset

Run the cell below to load the dataset.

Exercises

loadpath = 'dataset_session_754312389.nc'

dataset_path = Path(loadpath)

dataset = xr.load_dataset(dataset_path)
dataset

<xarray.Dataset> Size: 225MB Dimensions: (channel_depth: 23, trial_nr: 75, time_lfp: 1875, time_csd: 1875, time_whole_rec_lfp: 373124, time_whole_rec_ecp: 373124, unit_id_LGN: 27, time_spikes: 1500, unit_id_V1: 91, unit_id_LM: 13, stimulus_start_times: 75, stimulus_stop_times: 75, channel_id: 23) Coordinates: (12/13)

channel_depth (channel_depth) int64 184B 0 -40 -80 … -840 -880
trial_nr (trial_nr) int32 300B 0 1 2 3 4 5 … 70 71 72 73 74
time_lfp (time_lfp) float64 15kB -1e+03 -999.2 … 499.2
time_csd (time_csd) float64 15kB -1e+03 -999.2 … 499.2
time_whole_rec_lfp (time_whole_rec_lfp) float64 3MB 1.286e+06 … 1….
time_whole_rec_ecp (time_whole_rec_ecp) float64 3MB 1.286e+06 … 1…. … …
time_spikes (time_spikes) float64 12kB -999.5 -998.5 … 499.5
unit_id_V1 (unit_id_V1) int32 364B 951795075 … 951798053
unit_id_LM (unit_id_LM) int32 52B 951791074 … 951792163
stimulus_start_times (stimulus_start_times) float64 600B 1.29e+06 … …
stimulus_stop_times (stimulus_stop_times) float64 600B 1.29e+06 … 1…
channel_id (channel_id) int32 92B 850144538 … 850144362 Data variables: lfp_V1 (channel_depth, trial_nr, time_lfp) float64 26MB … csd_V1 (channel_depth, trial_nr, time_csd) float64 26MB … lfp_whole_recording_V1 (channel_depth, time_whole_rec_lfp) float64 69MB … ecp_whole_recording_V1 (channel_depth, time_whole_rec_ecp) float64 69MB … spike_counts_LGN (unit_id_LGN, trial_nr, time_spikes) int16 6MB 0 … spike_counts_V1 (unit_id_V1, trial_nr, time_spikes) int16 20MB 0 … spike_counts_LM (unit_id_LM, trial_nr, time_spikes) int16 3MB 0 …. pupil_width (trial_nr) float64 600B 39.12 40.53 … 46.57 44.31 run_speed (trial_nr) float64 600B 1.155 1.597 … 1.251 1.711 Attributes: time_unit: millisecond lfp_unit: Volt channel_depth_unit: micrometer note_channel_depth: Measured in distance from electrode closest t… sampling_frequency_lfp: 1250 sampling_frequency_spikes: 1000 sampling_frequency_unit: Hz stimulus_onset: 1000 stimulus_offset: 1250

Example: Calculate the time-resolved power spectrum for the constructed signal below using a Morlet wavelet. Plot the resulting spectrum using the pcolormesh function.

# provided

# Generate signal
sampling_frequency = 1000
durationA = 0.5
time_vecA = np.linspace(0,durationA,int(sampling_frequency*durationA)) # seconds
frequencyA = 5 # Hz
sineA = np.sin(2*np.pi*frequencyA*time_vecA)

durationB = 0.5
time_vecB = np.linspace(durationA, durationA+durationB, int(sampling_frequency*durationB)) # seconds
frequencyB = 15 # Hz
sineB = np.sin(2*np.pi*frequencyB*time_vecB)

time_vec = np.linspace(0,durationA+durationB,int(sampling_frequency*(durationA+durationB)))
sig = np.concatenate([sineA,sineB])

plt.figure()
plt.plot(time_vec,sig)
plt.title(f'Sinusoid signal that changes\nfrequency from {frequencyA} Hz to {frequencyB} Hz');

Solution

# solution

number_of_freqs = 10
bandwidth = 1.5
center_frequency = 1.0
freqs, spectrum = utils.calc_wavelet_spectrum(sig, sampling_frequency, number_of_freqs, bandwidth, center_frequency)

plt.pcolormesh(time_vec, freqs, spectrum, cmap = 'plasma', shading='gouraud')
plt.xlabel('Time (s)')
plt.ylabel('Frequencies (Hz)')
plt.title('Wavelet spectrum of signal');

Exercise: Calculate the time-resolved power spectrum for the same constructed sinusoid below using a Morlet wavelet. But this time change the number_of_freqs to 100. This parameter specifices how many frequencies to calculate the spectrum for. Plot the result using the pcolormesh function.

What happens to the frequency resolution when you compare to the plot above?

Solution

number_of_freqs = 100
bandwidth = 1.5
center_frequency = 1.0
freqs, spectrum = utils.calc_wavelet_spectrum(sig, sampling_frequency, number_of_freqs, bandwidth, center_frequency)

plt.pcolormesh(time_vec, freqs, spectrum, cmap = 'plasma', shading='gouraud')
plt.xlabel('Time (s)')
plt.ylabel('Frequencies (Hz)')
plt.title('Wavelet spectrum of signal');

The Morlet wavelet comes in many different shapes and sizes. There are two parameters that determine which kind of Morlet wavelet is used in the calculation of the spectrum: the bandwidth and the center frequency of the wavelet. The bandwidth determines how broad the wavelet is and the center frequency determines how many oscillations it has. Wavelets with different bandwidths and center frequencies are illustrated in the figure below. These two parameters determine the frequency- and time-resolution of the resulting power spectrum.

Exercise: Calculate the time-resolved power spectrum for the same constructed sinusoid below using a Morlet wavelet. Keep the number_of_frequencies parameter at 100, but change the center_frequency parameter to 0.5.

What happens to the frequency resolution now?

Solution

number_of_freqs = 100
bandwidth = 1.5
center_frequency = 0.5
freqs, spectrum = utils.calc_wavelet_spectrum(sig, sampling_frequency, number_of_freqs, bandwidth, center_frequency)

plt.pcolormesh(time_vec, freqs, spectrum, cmap = 'plasma', shading='gouraud')
plt.xlabel('Time (s)')
plt.ylabel('Frequencies (Hz)')
plt.title('Wavelet spectrum of signal');

Exercise: Change the center_frequency parameter to 2.0. Calculate and plot the time-resolved spectrum for the same constructed sinusoid below using a Morlet wavelet.

What happens to the frequency resolution? What happens to the time resolution?

Solution

number_of_freqs = 100
bandwidth = 1.5
center_frequency = 2.0
freqs, spectrum = utils.calc_wavelet_spectrum(sig, sampling_frequency, number_of_freqs, bandwidth, center_frequency)

plt.pcolormesh(time_vec, freqs, spectrum, cmap = 'plasma', shading='gouraud')
plt.xlabel('Time (s)')
plt.ylabel('Frequencies (Hz)')
plt.title('Wavelet spectrum of signal');

Exercise: Set the center_frequency parameter back to 1.0, but change the bandwidth parameter to 0.5. Calculate and plot the time-resolved spectrum for the same constructed sinusoid again.

What happens to the frequency resolution this time?

Solution

number_of_freqs = 100
bandwidth = 0.5
center_frequency = 1.0
freqs, spectrum = utils.calc_wavelet_spectrum(sig, sampling_frequency, number_of_freqs, bandwidth, center_frequency)

plt.pcolormesh(time_vec, freqs, spectrum, cmap = 'plasma', shading='gouraud')
plt.xlabel('Time (s)')
plt.ylabel('Frequencies (Hz)')
plt.title('Wavelet spectrum of signal');

Exercise: Change the bandwidth of the wavelet to 2.5. Calculate and plot the spectrum using pcolormesh.

What happens to the frequency resolution compared to the last plot?

Solution

number_of_freqs = 100
bandwidth = 2.5
center_frequency = 1.0
freqs, spectrum = utils.calc_wavelet_spectrum(sig, sampling_frequency, number_of_freqs, bandwidth, center_frequency)

plt.pcolormesh(time_vec, freqs, spectrum, cmap = 'plasma', shading='gouraud')
plt.xlabel('Time (s)')
plt.ylabel('Frequencies (Hz)')
plt.title('Wavelet spectrum of signal');

Comment: The exercises and plots above should demonstrate that if you increase the time resolution, you decrease the frequency resolution in your spectrum. And vice versa, if you increase the frequency resolution, you decrease the time resolution of your spectrum. There’s a fundamental trade-off between resolution in frequency and resolution in time.

The different parameters used for the Morlet wavelet determine where you strike the balance between frequency and time resolution. The frequency resolution increases and the time resolution goes down with greater bandwidth or greater center frequency.

Section 2: Time-resolved frequency analysis with morlet wavelets on spikes and LFP data

In this section, you will use Morlet wavelets to calculate the time-resolved spectra for spikes and LFP data.

Code	Description
`freqs, spectrum_amp = utils.calc_wavelet_spectrum(sig, sampling_frequency)`	Calculates the wavelet spectrum `spectrum` for a signal `sig` with a given `sampling_frequency`. `number_of_freqs` sets the number of frequencies the spectrum is calculated for. `bandwidth` and `center_frequency` sets the parameters for the Morlet wavelet.
`np.mean(data, axis = (dim_num))` or `data.mean(axis = dim_num)`	Calculate the average of the data across the `dim_num` dimension of the array.
`plt.pcolormesh(x, y, C, cmap = 'plasma', shading = 'gouraud')`	Make a 2D colormap of values in a 2D array (C) against x and y values. The optional parameters `cmap` and `shading` defines the colormap and the interpolation applied to the plot, respectively.

Run the cell below to extract the LFP data from the dataset.

Exercises

lfp_V1 = dataset['lfp_V1']
lfp_whole_rec = dataset['lfp_whole_recording_V1']
channels_depth = dataset['channel_depth']

sampling_frequency_lfp = dataset.sampling_frequency_lfp # Hz
time_lfp = lfp_V1.time_lfp.values # ms

We will stick to the parameters for the Morlet wavelet used in the first example in this section. That is, bandwidth = 1.5, and center_frequency = 1.0. Example: Calculate the time-resolved power spectrum of an LFP trace from the channel on index 8 in trial 0 using a morlet wavelet. Plot the result using pcolormesh.

ichan = 8
itrial = 0

number_of_freqs = 100
bandwidth = 1.5
center_frequency = 1.0
freqs, spectrum = utils.calc_wavelet_spectrum(lfp_V1[ichan, itrial], sampling_frequency_lfp, number_of_freqs, bandwidth, center_frequency)

# plot result using matplotlib's pcolormesh
plt.pcolormesh(time_lfp, freqs, spectrum, cmap = 'plasma', shading='gouraud')
plt.xlabel('Time from stimulus onset (ms)')
plt.ylabel('Frequencies (Hz)')
plt.title('Power spectrum')

Exercise: Calculate the time-resolved power spectrum of an LFP trace from the channel on index 2 in trial 0 using a Morlet wavelet. Plot the result using pcolormesh.

Is there a noticable increase in power for some frequencies at any time during the trial?

Solution

ichan = 2
itrial = 0

number_of_freqs = 100
bandwidth = 1.5
center_frequency = 1.0
freqs, spectrum = utils.calc_wavelet_spectrum(lfp_V1[ichan, itrial], sampling_frequency_lfp, number_of_freqs, bandwidth, center_frequency)

plt.pcolormesh(time_lfp, freqs, spectrum, cmap = 'plasma', shading='gouraud')
plt.xlabel('Time from stimulus onset (ms)')
plt.ylabel('Frequencies (Hz)')

Exercise: Calculate the spectrum of an LFP from the channel on index 8 in all trials (you don’t need to use a for loop). Calculate the trial-averaged spectrum and plot the result.

Tip: You might want to decrease the number_of_freqs parameter if you feel that it takes too long to calculate the spectrum.

Solution

ichan = 8

number_of_freqs = 50
bandwidth = 1.5
center_frequency = 1.0
freqs, spectrum = utils.calc_wavelet_spectrum(lfp_V1[ichan], sampling_frequency_lfp, number_of_freqs, bandwidth, center_frequency)

spectrum_trial_avg = np.mean(spectrum, axis = 1)

plt.pcolormesh(time_lfp, freqs, spectrum_trial_avg, cmap = 'plasma', shading='gouraud')
plt.xlabel('Time from stimulus onset (ms)')
plt.ylabel('Frequencies (Hz)')

Run the cell below to extract the spike data.

spike_counts_LGN = dataset['spike_counts_LGN']
spike_counts_V1 = dataset['spike_counts_V1']
spike_counts_LM = dataset['spike_counts_LM']

time_spikes = spike_counts_LM.time_spikes
sampling_frequency_spikes = dataset.sampling_frequency_spikes

Exercise: Calculate the time-resolved power spectrum for the trial-averaged firing rate of the LM cell at index 6. Plot the result using pcolormesh.

Is there a notable increase in power at any time point during the trial?

NB: You have to change the sampling frequency parameter from sampling_frequency_lfp to sampling_frequency_spikes in the calc_wavelet_spectrum function. And the time in the pcolormesh function also needs to be changed to time_spikes.

Solution

icell = 6

number_of_freqs = 50
bandwidth = 1.5
center_frequency = 1.0
freqs, spectrum = utils.calc_wavelet_spectrum(spike_counts_LM[icell].mean('trial_nr'), sampling_frequency_spikes, number_of_freqs, bandwidth, center_frequency)

plt.pcolormesh(time_spikes, freqs, spectrum, cmap = 'plasma', shading='gouraud')
plt.xlabel('Time from stimulus onset (ms)')
plt.ylabel('Frequencies (Hz)')

Comment: For time-resolved spectra of spikes, you will probably get more useful spectra if you average across trials (or cells) to get the firing rate first and then calculate the spectrum, like in the exercise above. However, the fact that the spike trains won’t all be in phase will mean that the spectral estimate isn’t entirely accurate.

Section 3: Cone of influence

When the signal is convolved with a wavelet, the wavelet will partly be outside of the signal window close to the edges (the start and end point of your signal). This will lead to inaccurate spectral estimates in these regions. Which parts of your spectral estimate may be influenced by such edge-effect artifacts can be calculated with the formula for a “cone of influence”. Regions of your spectrum outside this cone can’t be trusted and should be ignored when examining your time-resolved power spectrum. The formula for the cone of influence is given as:

t_{COI} = \frac{\sigma}{2\pi}\frac{\sqrt{2}}{f}

where $\sigma$ is the bandwidth of the Morlet wavelet and $f$ is frequencies for which the spectrum is calculated.

Code	Description
`freqs, spectrum_amp = utils.calc_wavelet_spectrum(sig, sampling_frequency)`	Utility function to calculate the wavelet spectrum `spectrum` for a signal `sig` with a given `sampling_frequency`. `number_of_freqs` sets the number of frequencies the spectrum is calculated for. `bandwidth` and `center_frequency` sets the parameters for the Morlet wavelet.
`t_coi = (bandwidth) / (2np.pi) np.sqrt(2) / freqs`	Calculate the cone of influence (COI) for a spectrum computed with a Morlet wavelet with a given bandwidth.
`utils.plot_spectrum_with_coi(time_lfp, freqs, spectrum, t_coi, apply_shading = True)`	Makes a plot of the wavelet spectrum with the cone of influence superimposed.

Exercises

Example: Calculate the cone of influence (COI) for the Morlet wavelet spectrum computed for a single channel, single trial LFP trace (ichan = 2, itrial = 10). Plot the COI on top of the the wavelet spectrum. You can use the utility functions to calculate the spectrum and to make the plot.

ichan = 2
itrial = 10

number_of_freqs = 50
bandwidth = 1.5
center_frequency = 1.0
freqs, spectrum = utils.calc_wavelet_spectrum(lfp_V1[ichan, itrial], sampling_frequency_lfp, number_of_freqs, bandwidth, center_frequency)

t_coi = (bandwidth) / (2*np.pi) * np.sqrt(2) / freqs * 1E3

# plot spectrum and cone of influence
utils.plot_spectrum_with_coi(time_lfp, freqs, spectrum, t_coi)

Exercise: Change the bandwidth parameter to 3.5 and calculate of the spectrum for the same LFP trace. Calculate and plot the COI on top of the spectrum below. What happens to the COI?

Solution

ichan = 2
itrial = 10


number_of_freqs = 50
bandwidth = 3.5
center_frequency = 1.0
freqs, spectrum = utils.calc_wavelet_spectrum(lfp_V1[ichan, itrial], sampling_frequency_lfp, number_of_freqs, bandwidth, center_frequency)

t_coi = (bandwidth) / (2*np.pi) * np.sqrt(2) / freqs * 1E3

# plot spectrum and cone of influence
utils.plot_spectrum_with_coi(time_lfp, freqs, spectrum, t_coi, apply_shading = True)

One way to avoid having to disregard a significant part of the spectrum because it’s outside of the COI is to calculate it for a larger time window.

Exercise: Below, the spectrum is calculated with bandwidth=3.5 for an LFP trace for a quarter of the whole recording period instead of just a single trial. Calculate and plot the COI for this spectrum. What happens to the COI in this case? Is it a large part of the spectrum that you have to disregard because it’s outside of the COI?

Hint: It’s probably best to set the apply_shading parameter to False in the plot_spectrum_with_coi function to make the plotting go faster.

ichan = 2
lfp_trace = lfp_whole_rec[ichan, :int(lfp_whole_rec.shape[1]/4)]
time_lfp_part_of_rec = lfp_whole_rec.time_whole_rec_lfp.values[:len(lfp_trace)]

number_of_freqs = 20
bandwidth = 3.5
center_frequency = 1.0
freqs, spectrum = utils.calc_wavelet_spectrum(lfp_trace, sampling_frequency_lfp, number_of_freqs, bandwidth, center_frequency)

Solution

# solution 

t_coi = (bandwidth) / (2*np.pi) * np.sqrt(2) / freqs * 1E3

# plot spectrum and cone of influence
utils.plot_spectrum_with_coi(time_lfp_part_of_rec, freqs, spectrum, t_coi, apply_shading = False)

Section 4: Reorganizing spectrum calculated for whole period into trials.

If you want to look at the spectrum in individual trials or the trial average, but you have calculated the spectrum for a whole recording period to limit the influence of edge effects, you have to reorganize the calculated spectrum into trials again. To do that, you can use information about the time points at which the stimuli are presented and logical indexing.

Code	Description
`mask = np.logical_and(array > start, array < stop)`	Create logical mask where both conditions are met. The array, e.g. a time vector, is bigger than the starting point and smaller than the stopping point.

Run cell below to load meta data on stimulus presentation start times and LFP from whole stimulus presentation period

Exercises

stim_start_times = dataset['stimulus_start_times'].values

ichan = 2
lfp_trace = lfp_whole_rec[ichan]
time_lfp_whole_rec = lfp_whole_rec.time_whole_rec_lfp.values

Example: Get LFP from 100 ms before the first stimulus presentation to 250 ms after onset using logical indexing. Plot the LFP in this selected window.

time_before_onset = 100
time_after_onset = 250

# get time of first stimulus presentation
idx_stim_start = 0
time_stim_onset = stim_start_times[idx_stim_start]

# create boolean mask for time window around stimulus onset
mask_stim_window = np.logical_and(time_lfp_whole_rec > (time_stim_onset-time_before_onset), 
                                  time_lfp_whole_rec < (time_stim_onset+time_after_onset))

# get time points in window
time_lfp_window = time_lfp_whole_rec[mask_stim_window]
time_lfp_window -= (time_lfp_window.min()+time_before_onset)

# get lfp in time window
lfp_window = lfp_whole_rec[:, mask_stim_window]

# plot lfp in selected time window
plt.pcolormesh(time_lfp_window, channels_depth, lfp_window, cmap='bwr', shading='gouraud')
plt.xlabel('Time from stimulus onset (ms)')
plt.ylabel(r'Depth ($\mu$m)');

Exercise: Get LFP from 500 ms before the first stimulus presentation to 500 ms after onset of the 5th stimulus presentation using logical indexing. Plot the LFP in this selected window.

Solution

time_before_onset = 500
time_after_onset = 500

# get time of selected stimulus presentation
idx_stim_start = 4
time_stim_onset = stim_start_times[idx_stim_start]

# create boolean mask for time window around stimulus onset
mask_stim_window = np.logical_and(time_lfp_whole_rec > (time_stim_onset-time_before_onset), 
                                  time_lfp_whole_rec < (time_stim_onset+time_after_onset))

# get time points in window
time_lfp_window = time_lfp_whole_rec[mask_stim_window]
time_lfp_window -= (time_lfp_window.min()+time_before_onset)

# get lfp in time window
lfp_window = lfp_whole_rec[:, mask_stim_window]

# plot lfp in selected time window
plt.pcolormesh(time_lfp_window, channels_depth, lfp_window, cmap='bwr', shading='gouraud')
plt.xlabel('Time from stimulus onset (ms)')
plt.ylabel(r'Depth ($\mu$m)');

Exercise: Reorganize the array containing LFP from the whole recording into an array organized by trials with 250 ms before and after stimulus onset by using the stimulus start time information.

Hint: Loop through all trials, get the lfp 250 ms before and after stimulus onset in each trial, and store it in a variable that contains the LFP in this window for all stimulus presentations. Plot the trial averaged LFP in this window.

Solution

time_before_onset = 250
time_after_onset = 250

n_trials = lfp_V1.shape[1]
lfp_windows = []
for idx_stim_start in range(n_trials-1): # 1 is subtracted just because the data was cut too early in the last trial
    # get time of selected stimulus presentation
    time_stim_onset = stim_start_times[idx_stim_start]

    # create boolean mask for time window around stimulus onset
    mask_stim_window = np.logical_and(time_lfp_whole_rec > time_stim_onset-time_before_onset, 
                                      time_lfp_whole_rec < time_stim_onset+time_after_onset)

    # get lfp in time window and store in list
    lfp_windows.append(lfp_whole_rec[:, mask_stim_window])

# get time points in window
time_lfp_window = time_lfp_whole_rec[mask_stim_window]
time_lfp_window -= (time_lfp_window.min()+time_before_onset)

# calculate trial averaged lfp
lfp_window_trial_avg = np.mean(lfp_windows, axis = 0)

plt.pcolormesh(time_lfp_window, channels_depth, lfp_window_trial_avg, cmap='bwr', shading='gouraud')
plt.xlabel('Time from stimulus onset (ms)')
plt.ylabel(r'Depth ($\mu$m)');

Exercise: Reorganize the spectrum calculated for the whole recording below into trials by using the stimulus start time information. Recipe: Loop through all trials, get the lfp 250 ms before and after stimulus onset in each trial, and store it in a variable that contains the spectrum in this window for all stimulus presentations. Plot the trial averaged spectrum in this window.

Note: Calculating the spectrum for the whole recording period will take some time.

# provided
number_of_freqs = 20
bandwidth = 3.5
center_frequency = 1.0
freqs, spectrum = utils.calc_wavelet_spectrum(lfp_trace, sampling_frequency_lfp, number_of_freqs, bandwidth, center_frequency)

Solution

# solution

n_trials = lfp_V1.shape[1]

# initialize empty list to store lfp from each trial
spectrum_windows = []
for idx_stim_start in range(n_trials-1): # 1 is subtracted because the data was cut too early in the last trial
    
    # get time of selected stimulus presentation
    time_stim_onset = stim_start_times[idx_stim_start]

    # create boolean mask for time window around stimulus onset
    mask_stim_window = np.logical_and(time_lfp_whole_rec > time_stim_onset-250, time_lfp_whole_rec < time_stim_onset+250)

    # get spectrum in time window and store in list
    spectrum_windows.append(spectrum[:, mask_stim_window])

# get time points in window
time_lfp_window = time_lfp_whole_rec[mask_stim_window]
time_lfp_window -= (time_lfp_window.min()+250)

# calculate trial averaged spectrum
spectrum_window_trial_avg = np.mean(spectrum_windows, axis = 0)

plt.pcolormesh(time_lfp_window, freqs, spectrum_window_trial_avg, cmap='inferno', shading='gouraud')
plt.xlabel('Time from stimulus onset (ms)')
plt.ylabel(r'Depth ($\mu$m)');