iBOTS Learn: Filtering Methods in Analysis of LFP

Filtering Methods in Analysis of LFP

Authors

Dr. Atle Rimehaug | Dr. Nicholas Del Grosso

Download Materials

In this notebook, you will apply filtering methods to real LFP data and combine the results of this filtering with other analyses.

Lowpass filters - keep part of signal below a certain frequency.
Highpass filters - keep part of signal above a certain frequency.
Bandpass filters - keep part of signal inside a certain frequency range.
Bandstop filters - keep part of signal outside a certain frequency range.

Credit: allaboutcircuits.com

Setup

Import Libraries

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

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 style_plot(title, xlabel, ylabel):
    '''
    Styles plots displaying filters. Adds labels, grids and shows legends

    Args:
        title: title of plot
        xlabel: label x-axis
        ylabel: label y-axis
    '''
    
    plt.title(title)
    plt.xlabel(xlabel)
    plt.ylabel(ylabel)
    plt.margins(0, 0.1)
    plt.grid(which='both', axis='both')
    plt.legend();

class utils:
    style_plot = style_plot

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: Filter to obtain LFP from recorded extracellular potential

Lowpass filters are used to obtain LFP from the recorded extracellular potential but there’s no universally agreed upon convention on which cutoff frequency to use. Most labs will use a cutoff frequency between 100 Hz and 500 Hz. If you use higher cutoff frequencies, you are more likely to include multi-unit activity (MUA) - that is, spikes - in your LFP. And spikes is noise if you want to study LFP alone. However, if you use a low cutoff frequency, you may lose important information in your LFP. Your choice of a cutoff frequency will ultimately depend on your analysis, but 300 Hz can be seen as a reasonable default value. It can also be a good idea to experiment with different cutoff values to see how that may affect your analysis.

Code	Description
`b, a = signal.butter(filter_order, cutoff_frequency, btype = 'lowpass', fs = sampling_frequency)`	Create the coefficients `b` and `a` to be used in the lowpass filtering of the signal.
`w, h = signal.freqz(b, a, fs=sampling_frequency)`	Compute the frequency response of the filter with coefficients `b` and `a`. `h` is the frequency response (the filter’s magnitude at a given frequency), `w` is the frequencies the response is calculated for.
`lfp_trace = signal.filtfilt(b, a, ecp_trace))`	Apply the constructed filter to the `ecp` trace from a single channel using the coefficients `b` and `a`.
`lfp = signal.filtfilt(b, a, ecp, axis = time_dim))`	Apply the constructed filter to an array containing the `ecp` recorded from all channels using the coefficients `b` and `a`. When the ecp data is multidimensional, the axis the filter is applied to needs to be specified to be along the time dimension of the data (default is `axis = -1`).
`plt.plot(x, y)`	Plot y values against x values.
`plt.imshow(data, aspect = 'auto', cmap = 'selected_colormap')`	Plots a multidimensional array as an image. `aspect` sets the aspect ratio of the axis. `cmap` sets the colormap.

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

ecp_whole_rec = dataset['ecp_whole_recording_V1']
time_ecp = ecp_whole_rec.time_whole_rec_ecp
time_lfp = time_ecp

sampling_frequency_lfp = dataset.sampling_frequency_lfp

Exercise: Construct and plot a lowpass butterworth filter with a cutoff frequency at 300 Hz. Apply it to the ECP trace from channel index 2 in the recording and plot both the original and the filtered trace together in a new plot. The filtered trace consitutes the LFP.

Solution

filter_order = 5 
cutoff_frequency = 300 # Hz

# creating butter filter
b, a = signal.butter(filter_order, cutoff_frequency, btype = 'low', fs = sampling_frequency_lfp)

w, h = signal.freqz(b, a, fs=sampling_frequency_lfp)

plt.plot(w, np.abs(h), label='Filter')

plt.title('Butterworth filter frequency response')
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude')
plt.margins(0, 0.1)
plt.grid(which='both', axis='both')
plt.axvline(cutoff_frequency, color='green', label='Cutoff frequency')
plt.legend();

ichan = 2
ecp_single_chan = ecp_whole_rec[ichan]

filter_order = 5 #
cutoff_frequency = 300 # Hz

# creating butter filter
b, a = signal.butter(filter_order, cutoff_frequency, btype = 'lowpass', fs = sampling_frequency_lfp)

lfp_single_chan = signal.filtfilt(b, a, ecp_single_chan)

plt.plot(time_ecp, ecp_single_chan, label = 'ECP trace')
plt.plot(time_ecp, lfp_single_chan, label = 'LFP trace')
plt.legend()
plt.xlabel('Time (s)')
plt.ylabel('Amplitude (V)');

Exercise: Construct and apply a lowpass butterworth filter (you don’t need to plot it) with a cutoff frequency at 100 Hz to the ECP trace from channel index 2 in the recording. Plot the original and the filtered signals together.

Solution

ichan = 2
ecp_single_chan = ecp_whole_rec[ichan]

filter_order = 5
cutoff_frequency = 100 # Hz

# creating butter filter
b, a = signal.butter(filter_order, cutoff_frequency, btype = 'low', fs = sampling_frequency_lfp)

lfp_single_chan = signal.filtfilt(b, a, ecp_single_chan)

plt.plot(time_ecp, ecp_single_chan, label = 'ECP trace')
plt.plot(time_ecp, lfp_single_chan, label = 'LFP trace')
plt.legend()
plt.xlabel('Time (s)')
plt.ylabel('Amplitude (V)');

Demo:

It’s a bit easier to see the effects of applying the lowpass filter if we limit the x-axis to focus on a smaller time period, like in the plot below.

ichan = 2
ecp_single_chan = ecp_whole_rec[ichan]

filter_order = 5
cutoff_frequency = 100 # Hz

# creating butter filter
b, a = signal.butter(filter_order, cutoff_frequency, btype = 'low', fs = sampling_frequency_lfp)

lfp_single_chan = signal.filtfilt(b, a, ecp_single_chan)

plt.plot(time_ecp, ecp_single_chan, label = 'ECP trace')
plt.plot(time_ecp, lfp_single_chan, label = 'LFP trace')
plt.legend()
plt.xlabel('Time (s)')
plt.ylabel('Amplitude (V)')
plt.xlim(1.3*1E6,1.301*1E6);

Exercise: Construct and apply a butterworth filter with cutoff frequency 300 Hz to LFP from all channels. Hint: By specifying the axis in the data to which the filter is applied in the filtfilt function, you don’t need to select the trace from a single channel.

Solution

filter_order = 5
cutoff_frequency = 100 # Hz

# creating butter filter
b, a = signal.butter(filter_order, cutoff_frequency, btype = 'low', fs = sampling_frequency_lfp)

lfp = signal.filtfilt(b, a, ecp_whole_rec, axis=-1)

Exercise: Construct and apply a highpass butterworth filter with a cutoff frequency at 600 Hz to the ECP. Plot single channel traces from the ECP and the MUA (the resulting highpass-filtered signal) together.

This high-frequency part of the signal should mainly contain the multi-unit activity (MUA) - the spikes - in the extracellular potential. Is the amplitude of the resulting MUA signal what you would expect when you compare it to the amplitude of the LFP above?

Hint: Change the btype parameter to create a highpass filter.

Solution

filter_order = 5
cutoff_frequency = 600 # Hz
sampling_frequency = 1250 # Hz

# creating butter filter
b, a = signal.butter(filter_order, cutoff_frequency, btype = 'highpass', fs = sampling_frequency_lfp)

mua_whole_rec = signal.filtfilt(b, a, ecp_whole_rec, axis=-1)

ichan = 2
ecp_single_chan = ecp_whole_rec[ichan]
mua_single_chan = mua_whole_rec[ichan]

plt.plot(time_ecp, ecp_single_chan, label = 'Original signal')
plt.plot(time_ecp, mua_single_chan, label = 'Filtered signal')
plt.legend()
plt.xlabel('Time (s)')
plt.ylabel('Amplitude (V)');

Section 2: Apply a bandpass filter to extract signal frequencies within specific range.

Code	Description
`b, a = ss.butter(filter_order, [start_freq_band, end_freq_band], btype = 'bandpass', fs = sampling_frequency)`	Create the coefficients `b` and `a` to be used in the bandpass (or bandstop) filtering of the signal. To create a bandstop filter instead of a bandpass filter, simply set `btype = 'bandstop'`
`w, h = ss.freqz(b, a, fs=sampling_frequency)`	Compute the frequency response of the filter with coefficients `b` and `a`. `h` is the frequency response (the filter’s magnitude at a given frequency), `w` is the frequencies at which the response is calculated.
`lfp_trace = ss.filtfilt(b, a, ecp_trace))`	Apply the constructed filter to the ecp signal at a specific channel using the coefficients `b` and `a`.
`plt.plot(x, y)`	Plot y values against x values.
`plt.imshow(data, aspect = 'auto', cmap = 'selected_colormap')`	Plots a multidimensional array as an image. `aspect` sets the aspect ratio of the axis. `cmap` sets the colormap.

Load LFP from whole recording

Exercises

lfp_whole_rec = dataset['lfp_whole_recording_V1']
time_lfp = lfp_whole_rec.time_whole_rec_lfp
time_lfp = time_ecp

Exercise: Construct, plot, and apply a bandpass butterworth filter with a low cutoff at 250 Hz and a high cutoff at 500 Hz to the whole LFP data. Plot a single channel trace from the original and the filtered signal together in one plot.

Solution

filter_order = 5 
sampling_frequency = 1250 # Hz
lowcut = 250 # Hz
highcut = 500 # Hz

# creating butter filter
b, a = signal.butter(filter_order, [lowcut, highcut], btype = 'band', fs = sampling_frequency_lfp)

w, h = signal.freqz(b, a, fs=sampling_frequency) 

plt.plot(w, np.abs(h), label='Filter')
plt.axvline(lowcut, color='red', label='Lower band cutoff')
plt.axvline(highcut, color='green', label='Upper band cutoff')

utils.style_plot(title='Butterworth filter frequency response', xlabel = 'Frequency [Hz]', ylabel='Amplitude')

lfp_bandpass_filtered = signal.filtfilt(b, a, lfp_whole_rec, axis = -1)

ichan = 5
lfp_single_chan = lfp_whole_rec[ichan]
lfp_bandpass_single_chan = lfp_bandpass_filtered[ichan]
plt.plot(time_lfp, lfp_single_chan, label = 'Original signal')
plt.plot(time_lfp, lfp_bandpass_single_chan , label = 'Filtered signal')
plt.legend()
plt.xlabel('Time (s)')
plt.ylabel('Amplitude (V)');

Exercise: Construct and apply a bandpass butterworth filter with a low cutoff at 50 Hz and a high cutoff at 70 Hz to the whole LFP data. Plot a single channel trace from the original and the filtered signal together in one plot.

Solution

filter_order = 5
sampling_frequency = 1250 # Hz
lowcut = 50 # Hz
highcut = 70 # Hz

# creating butter filter
b, a = signal.butter(filter_order, [lowcut, highcut], btype = 'bandpass', fs = sampling_frequency_lfp)

bandpass_filtered_lfp = signal.filtfilt(b, a, lfp_whole_rec, axis = -1)

ichan = 5
lfp_single_chan = lfp_whole_rec[ichan]
bp_lfp_single_chan = bandpass_filtered_lfp[ichan]
plt.plot(time_lfp, lfp_single_chan, label = 'Original signal')
plt.plot(time_lfp, bp_lfp_single_chan, label = 'Filtered signal')
plt.legend()
plt.xlabel('Time (s)');

Exercise: Use the code below that gets the delta band LFP (1-4 Hz) to construct and apply filters to get the theta (4-8 Hz), alpha (8-13 Hz), beta (13-30 Hz), and gamma band (30-120 Hz) LFP. Comment: filter_order = 5 creates a bandpass filter that’s too sharp for the lower, narrower frequency bands (delta, alpha), it will result in NaNs. A lower filter order should be used for the low and narrow frequency bands. A higher filter_order is fine for the higher, broader frequency bands.

# provided

filter_order = 2

''' Filter for delta band'''
lowcut = 4 # Hz
highcut = 100 # Hz

# creating butter filter
b, a = signal.butter(filter_order, [lowcut, highcut], btype = 'bandpass', fs = sampling_frequency_lfp)

lfp_delta = signal.filtfilt(b, a, ecp_whole_rec, axis = -1)

Solution

''' Filter for theta band'''
filter_order = 2
lowcut = 4 # Hz
highcut = 8 # Hz

# creating butter filter
b, a = signal.butter(filter_order, [lowcut, highcut], btype = 'bandpass', fs = sampling_frequency_lfp)

lfp_theta = signal.filtfilt(b, a, ecp_whole_rec, axis = -1)

''' Filter for alpha band'''
filter_order = 2
lowcut = 8 # Hz
highcut = 13 # Hz

# creating butter filter
b, a = signal.butter(filter_order, [lowcut, highcut], btype = 'bandpass', fs = sampling_frequency_lfp)

lfp_alpha = signal.filtfilt(b, a, ecp_whole_rec, axis = -1)

''' Filter for beta band'''
filter_order = 5
lowcut = 13 # Hz
highcut = 30 # Hz

# creating butter filter
b, a = signal.butter(filter_order, [lowcut, highcut], btype = 'bandpass', fs = sampling_frequency_lfp)

lfp_beta = signal.filtfilt(b, a, ecp_whole_rec, axis = -1)

''' Filter for gamma band'''
filter_order = 5
lowcut = 30 # Hz
highcut = 120 # Hz

# creating butter filter
b, a = signal.butter(filter_order, [lowcut, highcut], btype = 'bandpass', fs = sampling_frequency_lfp)

lfp_gamma = signal.filtfilt(b, a, ecp_whole_rec, axis = -1)

Section 3: Relate frequencies to behavior

LFP oscillations in certain frequency bands have been shown to reflect both neural phenomena and behavior (and if they don’t reflect a behavior or a neural phenomenon, they reflect an artifact of the experiment). In this section you’ll see how filtering can be used to uncover associations between the power of an oscillation and behavioral data like running speed and pupil size.

Code	Description
`b, a = signal.butter(filter_order, [start_freq_band, end_freq_band], btype = 'bandpass', fs = sampling_frequency)`	Create the coefficients `b` and `a` to be used in the lowpass filtering of the signal.
`lfp = signal.filtfilt(b, a, ecp, axis = time_dim))`	Apply the constructed filter to an array containing the ecp recorded from all channels using the coefficients `b` and `a`. When the ecp data is multidimensional, the axis at which the filter is applied needs to be specified to be along the time dimension of the data.
`plt.plot(x, y)`	Plot y values against x values.
`fig, ax = plt.subplots(nrows, ncols, figsize = (width, height))`	Make a 2x2 grid of subplots with 2 rows and 2 columns. `ax` can be indexed to refer to a specific subplot, f.ex. the subplot on the first row and second column would be `ax[0,1]`.
`ax_twin = ax1.twinx()`	Makes an extra axis in a subplot (the subplot of `ax1`) such that data with different units can be plotted together in the same plot.
`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.

Load behavioral data: Pupil width and running speed

Exercises

# Extract pupil width
avg_pupil_width = dataset['pupil_width']

# Extract running speed
avg_run_speed = dataset['run_speed']

Load LFP data organized by trials

# extract LFP trial data

lfp_V1 = dataset['lfp_V1']

Example: Calculate the average¹ absolute amplitude of the LFP in each trial. Plot the average LFP amplitude from channel index 5 together with the average running speed in each trial. ¹The average over the whole time period in a trial, not the trial average.

lfp_avg_abs_amp = np.mean(np.abs(lfp_V1), axis = -1)

fig, ax = plt.subplots()

ax_run = ax
ax_run.plot(avg_run_speed, label = 'Running speed')

ax_lfp = ax.twinx()

ichan = 5
ax_lfp.plot(lfp_avg_abs_amp[ichan], color = 'g', label = 'LFP avg. absolute amplitude')

ax.set_xlabel('Trial #')
ax_run.set_ylabel('Running speed (a.u.)')
ax_lfp.set_ylabel('LFP (V)')

fig.legend();

Exercise: Plot the average absolute amplitude of the LFP from channel index 5 together with the average pupil width in each trial.

Solution

lfp_avg_abs_amp = np.mean(np.abs(lfp_V1), axis = -1)

fig, ax = plt.subplots()

ax_run = ax
ax_run.plot(avg_pupil_width, label = 'Pupil width')

ax_lfp = ax.twinx()

ichan = 5
ax_lfp.plot(lfp_avg_abs_amp[ichan], color = 'g', label = 'LFP avg. absolute amplitude')

ax.set_xlabel('Trial #')
ax_run.set_ylabel('Running speed (a.u.)')
ax_lfp.set_ylabel('LFP (V)')

fig.legend();

Exercise: Construct and apply a bandpass filter with a low cutoff at 50 Hz and a high cutoff at 70 to get the narrow gamma range LFP.

Solution

filter_order = 5 #
lowcut = 50 # Hz
highcut = 70 # Hz

# creating butter filter
b, a = signal.butter(filter_order, [lowcut, highcut], btype = 'bandpass', fs = sampling_frequency_lfp)
lfp_narrow_gamma = signal.filtfilt(b, a, lfp_V1, axis = -1)

Exercise: Plot the average amplitude of the filtered gamma range LFP together with the average pupil size in each trial.

Compare the overlap between pupil size and narrow gamma-band amplitude LFP in this plot to the overlap between pupil size and the total LFP amplitude in the exercise above. What do you notice?

Solution

fig, ax = plt.subplots()

ax_pupil = ax
ax_pupil.plot(avg_pupil_width, label = 'Pupil width')

ax_lfp = ax.twinx()

lfp_narrow_gamma_abs_amp = np.mean(np.abs(lfp_narrow_gamma[ichan]), axis = 1)

ax_lfp.plot(lfp_narrow_gamma_abs_amp, color = 'g', label = 'LFP avg. absolute amplitude')
ax.set_xlabel('Trial #')
ax_pupil.set_ylabel('Pupil width (a.u.)')
ax_lfp.set_ylabel('LFP (V)')

fig.legend()

Demo: Run the cell below the average running speed together with the average pupil size in each trial.

fig, ax = plt.subplots()

ax_run = ax
ax_run.plot(avg_run_speed, label = 'Running speed')

ax_pupil = ax.twinx()
ax_pupil.plot(avg_pupil_width, 'r', label = 'Pupil size')

ax.set_xlabel('Trial #')
ax_run.set_ylabel('Running speed (a.u.)')
ax_pupil.set_ylabel('Pupil size (a.u.)')

fig.legend()

Section 4: Apply bandstop (notch) filters to remove frequencies and to eliminate artifacts.

Overly consistent frequencies across a recording period - stripes across a whole spectrogram in all channels, for example - may indicate an artifact. In the figure below you can see an example of spectrograms computed for LFP that has such stripes. Those artifacts can be removed with a bandstop filter (also called a notch filter - a notch filter is simply a bandstop filter with a narrow range.)

Note: The illustration plot is from a different animal than the one used in this workshop.

Code	Description
`b, a = signal.butter(filter_order, [start_freq_band, end_freq_band], btype = 'bandpass', fs = sampling_frequency)`	Create the coefficients `b` and `a` to be used in the lowpass filtering of the signal.
`w, h = signal.freqz(b, a, fs=sampling_frequency)`	Compute the frequency response of the filter with coefficients `b` and `a`. `h` is the frequency response (the filter’s magnitude at a given frequency), `w` is the frequencies at which the response is calculated.
`lfp_trace = signal.filtfilt(b, a, ecp_trace))`	Apply the constructed filter to the ecp signal at a specific channel using the coefficients `b` and `a`.
`lfp_filtered = signal.filtfilt(b, a, lfp, axis = time_dim))`	Apply the constructed filter to an array containing the LFP recorded from all channels using the coefficients `b` and `a`. When the LFP data is multidimensional, the axis the filter is applied to needs to be specified to be along the time dimension of the data.
`plt.plot(x, y)`	Plot y values against x values.

Exercise: Construct, plot, and apply a bandstop butterworth filter with a low cutoff at 50 Hz and a high cutoff at 60 Hz to the LFP.

Solution

filter_order = 5
lowcut = 50 # Hz
highcut = 60 # Hz

# creating butter filter
b, a = signal.butter(filter_order, [lowcut, highcut], btype = 'bandstop', fs = sampling_frequency_lfp)

w, h = signal.freqz(b, a, fs=sampling_frequency_lfp) # freqz because digital filter

plt.plot(w, np.abs(h), label='Filter')

plt.title('Butterworth filter frequency response')
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude')
plt.margins(0, 0.1)
plt.grid(which='both', axis='both')
plt.axvline(lowcut, color='red', label='Lower band cutoff') # cutoff frequency
plt.axvline(highcut, color='green', label='Upper band cutoff') # cutoff frequency
plt.legend();

lfp_bandstopped = signal.filtfilt(b, a, lfp_whole_rec, axis = -1)

ichan = 5
lfp_single_chan = lfp_whole_rec[ichan]
bandstopped_lfp_single_chan = lfp_bandstopped[ichan]
plt.plot(time_lfp, lfp_single_chan, label = 'Original signal')
plt.plot(time_lfp, bandstopped_lfp_single_chan, label = 'Filtered signal', ls = '--')
plt.legend()
plt.xlabel('Time (s)')

Exercise: Construct and apply a bandstop butterworth filter with a low cutoff at 380 Hz and a high cutoff at 420 Hz to the LFP. Plot a single channel trace from the original LFP and the bandstop filtered LFP together.

Solution

filter_order = 5
lowcut = 300 # Hz
highcut = 350 # Hz
sampling_frequency = 1250

# creating butter filter
b, a = signal.butter(filter_order, [lowcut, highcut], btype = 'bandstop', fs = sampling_frequency)

lfp_bandstopped = signal.filtfilt(b, a, lfp_whole_rec, axis = -1)

ichan = 5
lfp_single_chan = lfp_whole_rec[ichan]
bandstopped_lfp_single_chan = lfp_bandstopped[ichan]
plt.plot(time_lfp, lfp_single_chan, label = 'Original signal')
plt.plot(time_lfp, bandstopped_lfp_single_chan, label = 'Filtered signal')
plt.legend()
plt.xlabel('Time (s)')

Exercise: In the figure provided in the introduction to this section, there’s a highly consistent and prominent horizontal stripe in the frequency spectrum at a little below 400 Hz in all channels. This is a candidate for an artifact. Remove this (potential) artifact from the LFP using a bandstop filter.

Solution

filter_order = 5
lowcut = 380 # Hz
highcut = 400 # Hz
sampling_frequency = 1250

# creating butter filter
b, a = signal.butter(filter_order, [lowcut, highcut], btype = 'bandstop', fs = sampling_frequency)

lfp_bandstopped = signal.filtfilt(b, a, lfp, axis = -1)