For Loops and Trial Averaged Frequency Spectra

Course Catalog

LFP Analysis in Python

Frequency Analysis

Authors

Dr. Atle Rimehaug | Dr. Nicholas Del Grosso

Download Materials

In this notebook, you will learn how to apply frequency analysis methods to real LFP and spike data.

Setup

Import Libraries

import xarray as xr
from pathlib import Path

import numpy as np
import matplotlib.pyplot as plt

from scipy import signal
from nitime.algorithms import multi_taper_psd

Download Dataset

import owncloud

owncloud.Client.from_public_link('https://uni-bonn.sciebo.de/s/M56DCkpyCJ9FRoi', folder_password="ibots"
).get_file('/', 'dataset_session_754312389.nc')

True

Section 1: Calculating Trial Averaged Frequency Spectra

To reduce the prominence of noise in our results, we often calculate the average across multiple trials. When your data from each trial or channel is of the same length, you can simply apply the functions to compute the power spectra directly to the data and you will get spectra calculated for each trial and channel (or cells, in the case of spike data) organized the same way as your original data.

Code	Description
`freqs, power_density = signal.welch(data, fs = sampling_frequency)`	Calculate the power spectral density using Welch’s method
`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.
`np.mean(data, axis = (dim_num1, dim_num2))` or `data.mean(axis = dim_num1, dim_num2)`	Calculate the average of the data across both the `dim_num1` and the `dim_num2` dimensions of the array.
`plt.plot(x, y)`	Plot y values against x values.
`plt.semilogx(x, y)`	Plot y values against x values with a logarithmic x-axis and linear y-axis.
`plt.semilogy(x, y)`	Plot y values against x values with a logarithmic y-axis and a linear x-axis.
`plt.loglog(x, y)`	Plot y values against x values with logarithmic x- and y-axis.
`plt.xlim([some_lower_xlim, some_upper_xlim])`	Limit the x-axis in a plot to be between the numbers `some_lower_xlim` and `some_upper_xlim`.
`np.arange(start, stop, step)`	Makes an array with all the integers between two values, with a given spacing

Load dataset

Run cells below to load dataset and extract relevant data

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

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

# extract data

lfp_V1 = dataset['lfp_V1']

channels_depth = dataset['channel_depth']
sampling_frequency_lfp = dataset.sampling_frequency_lfp # Hz
time_lfp = lfp_V1.time_lfp.values # ms

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

sampling_frequency_spikes = dataset.sampling_frequency_spikes # Hz
time_spikes = dataset['time_spikes'] # ms

lfp_V1

<xarray.DataArray 'lfp_V1' (channel_depth: 23, trial_nr: 75, time_lfp: 1875)> Size: 26MB
array([[[-8.51569421e-06, -1.15631923e-05, -3.76924291e-07, ...,
         -2.17654367e-06,  1.68119712e-07, -6.94200500e-06],
        [ 2.81771330e-08,  6.29403676e-06,  9.98029865e-06, ...,
          1.01304291e-05,  5.93295341e-06,  1.72969475e-05],
        [-1.27271427e-06, -3.68783112e-06, -4.46340366e-06, ...,
          3.65666426e-06,  2.96670978e-06,  1.24568597e-06],
        ...,
        [ 4.70565253e-06,  1.94317598e-06,  1.20455675e-06, ...,
          4.07503496e-06,  7.65054254e-06,  3.91673831e-06],
        [ 8.38144661e-07,  2.28123812e-06,  3.70896321e-06, ...,
          1.34013793e-05,  1.15462596e-05,  8.42411093e-06],
        [ 7.32479652e-06,  1.07711181e-06,  7.00885884e-06, ...,
          1.92065863e-05,  1.92065863e-05,  1.92065863e-05]],
       [[ 3.91952787e-06,  1.70056905e-06,  5.97275690e-06, ...,
          4.06974274e-06,  4.04094749e-06, -4.68055730e-06],
        [-4.10311367e-06,  5.00856269e-06,  6.20500275e-06, ...,
          8.75868311e-06,  5.69604080e-06,  8.57744575e-06],
        [ 1.09646506e-06,  1.62070262e-06,  1.10953040e-07, ...,
          1.01990671e-05, -8.36754580e-06, -1.43387447e-05],
...
        [-6.46738552e-06, -1.19701456e-05, -2.41882630e-05, ...,
         -2.31175915e-05, -2.37374262e-05, -7.53926446e-06],
        [ 2.91850049e-05,  2.81713398e-05,  5.17618757e-05, ...,
          4.44764295e-05,  4.08778688e-05,  3.83233192e-05],
        [-1.80642354e-05, -1.06068608e-05, -3.71177368e-06, ...,
          1.48063260e-05,  1.48063260e-05,  1.48063260e-05]],
       [[ 4.42062699e-05,  5.99208645e-05,  7.60618853e-05, ...,
          9.85508013e-05,  8.87923195e-05,  8.02998878e-05],
        [-6.30930699e-06,  2.14882295e-05,  5.40961203e-06, ...,
         -6.18517313e-05, -8.49961217e-05, -7.07025406e-05],
        [ 1.32347041e-04,  9.66830285e-05,  9.77152893e-05, ...,
         -3.68393584e-05, -2.08663734e-05, -1.27901811e-05],
        ...,
        [ 7.59530363e-07, -9.00525049e-06, -1.95382678e-05, ...,
         -2.78150664e-05, -1.46349189e-05,  3.65395950e-06],
        [ 3.22381374e-05,  3.90163375e-05,  6.39709157e-05, ...,
          5.17422904e-05,  4.94353600e-05,  3.90770304e-05],
        [-2.50190016e-05, -1.82145705e-05, -1.97245214e-05, ...,
          6.13148062e-06,  6.13148062e-06,  6.13148062e-06]]])
Coordinates:
  * channel_depth  (channel_depth) int64 184B 0 -40 -80 -120 ... -800 -840 -880
  * trial_nr       (trial_nr) int32 300B 0 1 2 3 4 5 6 ... 68 69 70 71 72 73 74
  * time_lfp       (time_lfp) float64 15kB -1e+03 -999.2 -998.4 ... 498.4 499.2

Example: Calculate the power spectrum for all trials and all cells containing spike data from LM with Welch’s method. Make a 2D colormap of the trial averaged power spectrum with the freqs on the x-axis and cell numbers in LM (cell_nrs_LM) on the y-axis.

freqs, psd_spikes_LM_trials = signal.welch(spike_counts_LM.values, sampling_frequency_spikes)

cell_nrs_LM = np.arange(spike_counts_LM.shape[0])
plt.pcolormesh(freqs, cell_nrs_LM, psd_spikes_LM_trials.mean(axis=1), cmap = 'plasma', shading='gouraud')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Cell nr.')
plt.xlim([0,100]);

Exercise: Calculate the power spectrum for all trials and all cells containing spike data from V1 with Welch’s method. Make a 2D colormap of the trial averaged power spectrum with the freqs on the x-axis and cell numbers in V1 (cell_nrs_V1) on the y-axis.

Hint: To focus on a subset of frequencies, for example the frequencies below 100 Hz, you can limit the x-axis to focus on them using plt.xlim([lower_limit, upper_limit]) (insert numbers for lower_limit and upper_limit).

Solution

freqs, psd_spikes_V1_trials = signal.welch(spike_counts_V1.values, sampling_frequency_spikes)

cell_nrs_V1 = np.arange(spike_counts_V1.shape[0])
plt.pcolormesh(freqs, cell_nrs_V1, psd_spikes_V1_trials.mean(axis=1), cmap = 'plasma', shading='gouraud')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Cell nr.')
plt.xlim([0,100]);

Exercise: Calculate the power spectrum for all trials and all cells containing spike data from LGN with Welch’s method. make a 2D colormap of the trial averaged power spectrum with the freqs on the x-axis and cell numbers on the y-axis.

Do you notice anything different with this plot compared to the equivalent plots for V1 and LM spikes? Is there a frequency where the power is notably higher for some cells?

Solution

freqs, psd_spikes_LGN_trials = signal.welch(spike_counts_LGN.values, sampling_frequency_spikes)

cell_nrs_LGN = np.arange(spike_counts_LGN.shape[0])
plt.pcolormesh(freqs, cell_nrs_LGN, psd_spikes_LGN_trials.mean(axis=1), cmap = 'plasma', shading='gouraud')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Cell nr.')
plt.xlim([0,100]);

Exercise: Calculate the power spectrum for all trials and channels containing LFP data from V1 with Welch’s method. Compute the trial averaged power spectrum and make a 2D colormap with freqs on the x-axis and channels_depth on the y-axis.

NB: Make sure to use the sampling frequency for the LFP (sampling_frequency_lfp) in the function for Welch’s method.

Solution

freqs, psd_lfp_V1 = signal.welch(lfp_V1.values, sampling_frequency_lfp)

plt.pcolormesh(freqs, channels_depth, psd_lfp_V1.mean(axis=1), cmap = 'plasma', shading='gouraud')
plt.xlabel('Frequency (Hz)')
plt.ylabel(r'Depth ($\mu$m)')
plt.xlim([0,100]);

Exercise: From the LFP spectrum computed in the exercise above, calculate the average power spectrum for V1 LFP across both channels and trials. Make a 1D plot of that average against frequencies.

Solution

psd_lfp_trial_avg = np.mean(psd_lfp_V1, axis=(0,1))
plt.plot(freqs, psd_lfp_trial_avg);

Exercise: It’s often useful (and common practice) to plot frequency spectra on logarithmic axes rather than on linear axes. Create the same plot as above, but with a logarithmic x-axis (and a linear y-axis).

Solution

plt.semilogx(freqs, psd_lfp_trial_avg)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power spectral density');

Exercise: Create the same plot as above, but with a logarithmic y-axis (and a linear x-axis).

At approximately which frequencies are there bumps in the spectrum?

Solution

plt.loglog(freqs, psd_lfp_trial_avg)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power spectral density');

Exercise: Create the same plot as above, but with both axes logarithmic.

Solution

plt.loglog(freqs, psd_lfp_trial_avg)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power spectral density');

Demo

The spectra of single spikes (delta pulse) is a flat line.

fake_spikes = np.zeros(1000)

fake_spikes[100] = 1

plt.plot(fake_spikes)

freqs, psd = signal.welch(fake_spikes, 1000)

plt.plot(freqs, psd)
plt.title('Power spectrum of single spike (delta pulse)')
plt.xlabel('Frequency (Hz)')

The power spectrum for a poisson process is also a flat, albeit noisier, line.

poisson_process = np.random.poisson(size=(100,1000))

freqs, psd = signal.welch(poisson_process, 1000)

plt.plot(freqs, psd.mean(axis=0))
plt.title('Power spectrum of poisson process')
plt.xlabel('Frequency (Hz)')

The power spectra of a delta impulse and a poisson process are useful to have in mind when we assess the power spectrum of spike trains. If the power spectrum of a spike train is a flat and noisy line, it indicates that the neuron fires as a poisson process.

Example: Calculate the power spectrum for V1 spikes and average across both cells and trials. Make a 1D plot of the average power spectral density against frequencies.

freqs, psd_spikes_V1_trials = signal.welch(spike_counts_V1.values, sampling_frequency_spikes)

plt.plot(freqs, psd_spikes_V1_trials.mean(axis=(0,1)))
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power spectral density');

Exercise: Calculate the power spectrum for LGN spikes and average across both cells and trials. Make a 1D plot of the power spectral density against frequencies.

Is there a correspondence between any of the frequencies where you find peaks in the power spectrum for LGN spikes and where you find peaks in the spectrum for V1 LFP?

Solution

freqs, psd_spikes_LGN_trials = signal.welch(spike_counts_LGN.values, sampling_frequency_spikes)

plt.plot(freqs, psd_spikes_LGN_trials.mean(axis=(0,1)))
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power spectral density')
plt.xlim(0,100);

Section 2: Iteration with For Loops

There are different kinds of loops in programming languages that are used for different purposes. When you know when to start and stop an iteration, for loops are the most appropriate. For example, when you want to start at the first trial and stop at the last trial, or start at the 5th trial and stop at the 10th trial, you can specify that with a for loop. Going through the 5th to the 10th trial would have the following structure in Python code:

for itrial in range(4, 10):
    # do some analysis of the data in the trial with index "itrial"

Code	Description
`for letter in ["a", "b", "c"]:` `print(letter)`	Print every `letter` in the list `["a", "b", "c"]`
`for i in range(5):` `print(i)`	Print every number `i` between 0 and 5 (not including 5)
`for i in range(2, 7):` `print(i)`	Print every number `i` between 2 and 7 (not including 7)
`for i in range(1, 10, 2):` `print(i)`	Print eververy second number `i` between 1 and 10 (not including 10)
`for x_i in x:` `print(x_i)`	Print every element `x_i` in the collection `x`
`for x_i in x:` `for y_i in y`: `print(x_i, y_i)`	Print every element `x_i` and `y_i` in the collections `x` and `y`
`new_list = []` or `new_list = list()`	Initialize an empty list named `new_list`
`new_list.append(element)`	Append an element to the list `new_list`

Exercise: Write a for loop that prints the numbers 0 to 9.

Solution

for num in range(9):
    print(num)

Exercise: Write a for loop that prints the numbers 5 to 10.

Solution

for i in range(5,10):
    print(i)

Exercise: Write a for loop that counts from 0 to 100 in steps of 10.

Solution

for j in range(0,100,10):
    print(j)

Exercise: Use a loop to print all fruits and colors, for example banana is yellow".

Solution

fruits = ["banana", "orange", "cherry"]
colors= ["yellow", "orange", "green"];

for fruit in fruits:
    for color in colors:
        print(fruit+' is '+color)

banana is yellow
banana is orange
banana is green
orange is yellow
orange is orange
orange is green
cherry is yellow
cherry is orange
cherry is green

Example: Write a loop that calculates the square of each element in numbers and stores the result in a new list called squared_numbers. Display the squared_numbers list at the end to check that it is correct.

numbers = [4, 6, 12, 50]

squared_numbers = []

for i in range(len(numbers)):
    num_squared = numbers[i]**2
    squared_numbers.append(num_squared)

squared_numbers

[16, 36, 144, 2500]

Exercise: Write a loop that multiplies each element in the list numbers by 4 and store the results in a new list called multiplied_numbers. Display multiplied_numbers to check that it is correct.

Solution

numbers = [4, 6, 12, 50]

multiplied_numbers = []
for num in numbers:
    mult_number = 4*num

    multiplied_numbers.append(mult_number)

multiplied_numbers

[16, 24, 48, 200]

Exercise: Write a loop that print out "Run experiment for animal n" for each animal id in the list below.

Solution

animal_ids = ["m23", "m59", "m145"]

for animal_id in animal_ids:
    print('Run experiment for animal '+animal_id)

Run experiment for animal m23
Run experiment for animal m59
Run experiment for animal m145

Exercise: Use a for loop within another for loop to print "This is block b" for each block in blocks and
"This is trial t" for each t in trials.

Solution

blocks = [1, 2, 3]
trials = [1, 2, 3, 4, 5, 6, 7, 8 , 9 ,10]

for block in blocks:
    print(f'This is block {block}')
    for trial in trials:
        print(f'This is trial {trial}')

This is block 1
This is trial 1
This is trial 2
This is trial 3
This is trial 4
This is trial 5
This is trial 6
This is trial 7
This is trial 8
This is trial 9
This is trial 10
This is block 2
This is trial 1
This is trial 2
This is trial 3
This is trial 4
This is trial 5
This is trial 6
This is trial 7
This is trial 8
This is trial 9
This is trial 10
This is block 3
This is trial 1
This is trial 2
This is trial 3
This is trial 4
This is trial 5
This is trial 6
This is trial 7
This is trial 8
This is trial 9
This is trial 10

Section 3: Analysis Inside For Loops

Sometimes, your data (for example from different trials) isn’t all of the same length, so you can’t apply the functions to do your analysis directly to the data with all trials together. Or you might want to skip trials that don’t meet certain conditions from your analysis. In those cases, you will often need to use for loops and carry out your analysis inside each loop. Here, you’ll get some practice doing your frequency analysis inside a loop and storing the results.

Code	Description
`freqs, psd = signal.welch(sig, fs = sampling_frequency)`	Calculate the power spectral density for a signal `sig` using Welch’s method
`plt.plot(x, y)`	Plot y values against x values.
`np.shape(data)` or `data.shape`	Gives you the dimensions of the data. The output for `data` containing 20 channels on the first dimension, 50 trials, and 1000 time points would be (20,50,1000). `data.shape[2]` would give you the number of the second dimension (1000 in this example).
`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.
`lfp_trace = lfp_V1[ichan, itrial]`	Get the lfp trace from the channel `ichan` in trial `itrial`.
`new_list = []` or `new_list = list()`	Initialize an empty list named `new_list`
`new_list.append(element)`	Append an element to the list `new_list`

Exercises

Example: Loop through the first five trials and use Welch’s method to calculate the power spectrum for the LFP trace from the channel at index 8 in each trial. Store the results in a list. Calculate the trial averaged frequency spectrum and plot it against the frequencies (freqs).

Note: the frequencies for which the spectrum is calculated will be identical for all trials and do not need to be stored in a list. The frequencies returned in the last iteration can be used.

psd_lfp_V1 = []

# select which channel to get the lfp trace from
ichan = 8

# loop through the first 5 trials
for itrial in range(0,5,1):

    lfp_trace = lfp_V1[ichan, itrial]

    # calculate the power spectrum for this lfp trace
    freqs, psd_lfp_trial = signal.welch(lfp_trace, fs = sampling_frequency_lfp)
    
    # store the psd calculated for lfp trace in this trial
    psd_lfp_V1.append(psd_lfp_trial)

psd_lfp_trial_avg = np.mean(psd_lfp_V1, axis = 0)

plt.loglog(freqs, psd_lfp_trial_avg)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power spectral density');

Exercise: Loop through trials 5 to 10 and calculate the frequency spectrum for the LFP trace from channel at index 20 in each trial. Store the results in a list. Calculate the trial averaged frequency spectrum and plot it against the frequencies (freqs).

Solution

psd_lfp_V1 = []
ichan = 20
for itrial in range(5,10,1):
    lfp_trace = lfp_V1[ichan, itrial]
    freqs, psd_lfp_trial = signal.welch(lfp_trace, fs = sampling_frequency_lfp)
    
    psd_lfp_V1.append(psd_lfp_trial)

psd_lfp_trial_avg = np.mean(psd_lfp_V1, axis = 0)

plt.loglog(freqs, psd_lfp_trial_avg)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power spectral density')

Exercise: Loop through all (75) trials, calculate the frequency spectrum for the LFP trace from channel at index 8 in each trial, and store the results in a list. Calculate the trial averaged frequency spectrum and plot it against the frequencies (freqs).

Hint: Instead of writing 75 in the range of the loop, you can use the np.shape(data)[some_number] or data.shape[some_number] to get the number of trials in the array containing lfp traces from all trials and channels. Trials are on dimension 1.

Solution

psd_lfp_V1 = []
ichan = 8
for itrial in range(lfp_V1.shape[1]):
    lfp_trace = lfp_V1[ichan, itrial]
    freqs, psd_trial = signal.welch(lfp_trace, fs = sampling_frequency_lfp)
    
    psd_lfp_V1.append(psd_trial)

psd_lfp_trial_avg = np.mean(psd_lfp_V1, axis = 0)

plt.loglog(freqs, psd_lfp_trial_avg)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power spectral density')

Example: Loop through all trials and calculate the frequency spectrum for the LGN spikes of the cell at index 5. Store the results in a list and calculate and plot the trial averaged frequency spectrum of the spikes of this cell.

icell = 5

psd_spikes_LGN_all_trials = []
for itrial in range(spike_counts_LGN.shape[1]):
    spikes = spike_counts_LGN[icell, itrial]
    freqs, psd = signal.welch(spikes, sampling_frequency_spikes)

    psd_spikes_LGN_all_trials.append(psd)

psd_spikes_LGN_trial_avg = np.mean(psd_spikes_LGN_all_trials, axis = 0)

plt.plot(freqs, psd_spikes_LGN_trial_avg)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power spectral density');

Exercise: Loop through all trials and calculate the frequency spectrum for the V1 spikes of the cell at index 5. Store the results in a list and calculate and plot the trial averaged frequency spectrum of the spikes of this cell.

Solution

icell = 5

psd_spikes_V1_all_trials = []
for itrial in range(spike_counts_V1.shape[1]):
    spikes = spike_counts_V1[icell, itrial]

    freqs, psd = signal.welch(spikes, sampling_frequency_spikes)

    psd_spikes_V1_all_trials.append(psd)

psd_V1_trial_avg = np.mean(psd_spikes_V1_all_trials, axis = 0)

plt.plot(freqs, psd_V1_trial_avg)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power spectral density')

Exercise: Loop through all cells and calculate the frequency spectrum for the LGN spikes of each cell in the first trial (itrial = 0) and store the results in a list. Calculate and plot the population averaged frequency spectrum of the spikes from all cells in this single trial.

Solution

itrial = 0

psd_spikes_all_units = []
for icell in range(spike_counts_LGN.shape[0]):
    spikes = spike_counts_LGN[icell, itrial]
    freqs, psd = signal.welch(spikes, sampling_frequency_spikes)

    psd_spikes_all_units.append(psd)

psd_population_avg = np.mean(psd_spikes_all_units, axis = 0)

plt.plot(freqs, psd_population_avg)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power spectral density')

Section 4: Nested For Loops for Multidimensional Data

If you need to go through both different channels and trials or both different cells and trials for an analysis, you will need to use a loop within a loop, or a nested for loop.

Code	Description
`for x_i in x:` `for y_i in y`: `print(x_i, y_i)`	Print every element `x_i` and `y_i` in the collections `x` and `y`
`new_list = []` or `new_list = list()`	Initialize an empty list named `new_list`
`new_list.append(element)`	Append an element to the list `new_list`
`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.

Exercises

Example: Loop through both the cells and the trials and calculate the power spectrum for the LGN spike count of each cell in each trial. Store the resulting power spectra for each cell and trial in lists. Calculate the trial average and make a 2D colormap plot of the results with frequencies on the x-axis and cells on the y-axis.

psd_spikes_LGN = []
for icell in range(spike_counts_LGN.shape[0]):
    psd_spikes_LGN_all_trials = []
    for itrial in range(spike_counts_LGN.shape[1]):
        spikes = spike_counts_LGN[icell, itrial]
        freqs, psd = signal.welch(spikes, sampling_frequency_spikes)

        psd_spikes_LGN_all_trials.append(psd)

    psd_spikes_LGN.append(psd_spikes_LGN_all_trials)

psd_spikes_LGN_trial_avg = np.mean(psd_spikes_LGN, axis = 1)

cell_nrs_LGN = np.arange(spike_counts_LGN.shape[0])
plt.pcolormesh(freqs, cell_nrs_LGN, psd_spikes_LGN_trial_avg, cmap='plasma', shading='gouraud')
plt.xlim([0,100])
plt.xlabel('Frequency (Hz)')
plt.ylabel('Cell nr.')
plt.title('Trial averaged power spectrum LGN spikes');

Exercise: Loop through both the cells and the trials and calculate the power spectrum for the LM spike count (spike_counts_LM) of each cell in each trial. Store the resulting power spectra for each cell and trial in lists. Calculate the trial average and make a 2D colormap plot of the results with frequencies on the x-axis and cells on the y-axis.

Solution

psd_spikes_LM = []
for icell in range(spike_counts_LM.shape[0]):
    psd_spikes_LM_all_trials = []
    for itrial in range(spike_counts_LM.shape[1]):
        spikes = spike_counts_LM[icell, itrial]
        freqs, psd = signal.welch(spikes, sampling_frequency_spikes)

        psd_spikes_LM_all_trials.append(psd)

    psd_spikes_LM.append(psd_spikes_LM_all_trials)

psd_spikes_LM_trial_avg = np.mean(psd_spikes_LM, axis = 1)

cell_nrs_LM = np.arange(spike_counts_LM.shape[0])
plt.pcolormesh(freqs, cell_nrs_LM, psd_spikes_LM_trial_avg, cmap='plasma', shading='gouraud')
plt.xlim([0,100])
plt.xlabel('Frequency (Hz)')
plt.ylabel('Cell nr.')
plt.title('Trial averaged power spectrum')

Exercise: Loop through both the cells and the trials and calculate the power spectrum for the V1 spike count (spike_counts_V1) of each cell in each trial. Store the resulting power spectra for each cell and trial in lists. Calculate the trial average and make a 2D colormap plot of the results with frequencies on the x-axis and cells on the y-axis.

Solution

psd_spikes_LGN = []
for icell in range(spike_counts_LGN.shape[0]):
    psd_spikes_LGN_all_trials = []
    for itrial in range(spike_counts_LGN.shape[1]):
        spikes = spike_counts_LGN[icell, itrial]
        freqs, psd = signal.welch(spikes, sampling_frequency_spikes)

        psd_spikes_LGN_all_trials.append(psd)

    psd_spikes_LGN.append(psd_spikes_LGN_all_trials)

psd_spikes_LGN_trial_avg = np.mean(psd_spikes_LGN, axis = 1)

cell_nrs_LGN = np.arange(spike_counts_LGN.shape[0])
plt.pcolormesh(freqs, cell_nrs_LGN, psd_spikes_LGN_trial_avg, cmap='plasma', shading='gouraud')
plt.xlim([0,100])
plt.xlabel('Frequency (Hz)')
plt.ylabel('Cell nr.')
plt.title('Trial averaged power spectrum')

Exercise: Loop through both the channels and the trials and calculate the frequency spectrum for all LFP traces. Store the resulting power spectra. Calculate the trial average, and make a 2D colormap plot of the results with frequencies on the x-axis and cells on the y-axis.

NB: Make sure to use the sampling frequency for the LFP (sampling_frequency_lfp) in the function for Welch’s method.

Solution

psd_lfp_V1_chans = []
for ichan in range(lfp_V1.shape[0]):
    psd_single_chan_all_trials = []
    for itrial in range(lfp_V1.shape[1]):
        lfp_trace = lfp_V1[ichan, itrial]
        freqs, psd = signal.welch(lfp_trace, fs = sampling_frequency_lfp)
        
        psd_single_chan_all_trials.append(psd)
        
    psd_lfp_V1_chans.append(psd_single_chan_all_trials)

psd_lfp_trial_avg_chans = np.mean(psd_lfp_V1_chans, axis = 1)

plt.pcolormesh(freqs, channels_depth, psd_lfp_trial_avg_chans, cmap='plasma', shading='gouraud')
plt.xlim([0,100])
plt.xlabel('Frequency (Hz)')
plt.ylabel(r'Depth ($\mu$m)')