LFP/ECP Analysis

The lfp module provides tools for analyzing local field potentials (LFP) and extracellular potentials (ECP).

bmtool.analysis.lfp.load_ecp_to_xarray(ecp_file, demean=False)

Load ECP data from an HDF5 file (BMTK sim) into an xarray DataArray.

Parameters:

ecp_file : str Path to the HDF5 file containing ECP data. demean : bool, optional If True, the mean of the data will be subtracted (default is False).

Returns:

xr.DataArray An xarray DataArray containing the ECP data, with time as one dimension and channel_id as another.

Source code in bmtool/analysis/lfp.py

def load_ecp_to_xarray(ecp_file: str, demean: bool = False) -> xr.DataArray:
    """
    Load ECP data from an HDF5 file (BMTK sim) into an xarray DataArray.

    Parameters:
    ----------
    ecp_file : str
        Path to the HDF5 file containing ECP data.
    demean : bool, optional
        If True, the mean of the data will be subtracted (default is False).

    Returns:
    -------
    xr.DataArray
        An xarray DataArray containing the ECP data, with time as one dimension
        and channel_id as another.
    """
    with h5py.File(ecp_file, "r") as f:
        ecp = xr.DataArray(
            f["ecp"]["data"][()].T,
            coords=dict(
                channel_id=f["ecp"]["channel_id"][()],
                time=np.arange(*f["ecp"]["time"]),  # ms
            ),
            attrs=dict(
                fs=1000 / f["ecp"]["time"][2]  # Hz
            ),
        )
    if demean:
        ecp -= ecp.mean(dim="time")
    return ecp

bmtool.analysis.lfp.ecp_to_lfp(ecp_data, cutoff=250, fs=10000, downsample_freq=1000)

Apply a low-pass Butterworth filter to an xarray DataArray and optionally downsample. This filters out the high end frequencies turning the ECP into a LFP

Parameters:

ecp_data : xr.DataArray The input data array containing LFP data with time as one dimension. cutoff : float The cutoff frequency for the low-pass filter in Hz (default is 250Hz). fs : float, optional The sampling frequency of the data (default is 10000 Hz). downsample_freq : float, optional The frequency to downsample to (default is 1000 Hz).

Returns:

xr.DataArray The filtered (and possibly downsampled) data as an xarray DataArray.

Source code in bmtool/analysis/lfp.py

def ecp_to_lfp(
    ecp_data: xr.DataArray, cutoff: float = 250, fs: float = 10000, downsample_freq: float = 1000
) -> xr.DataArray:
    """
    Apply a low-pass Butterworth filter to an xarray DataArray and optionally downsample.
    This filters out the high end frequencies turning the ECP into a LFP

    Parameters:
    ----------
    ecp_data : xr.DataArray
        The input data array containing LFP data with time as one dimension.
    cutoff : float
        The cutoff frequency for the low-pass filter in Hz (default is 250Hz).
    fs : float, optional
        The sampling frequency of the data (default is 10000 Hz).
    downsample_freq : float, optional
        The frequency to downsample to (default is 1000 Hz).

    Returns:
    -------
    xr.DataArray
        The filtered (and possibly downsampled) data as an xarray DataArray.
    """
    # Bandpass filter design
    nyq = 0.5 * fs
    cut = cutoff / nyq
    b, a = signal.butter(8, cut, btype="low", analog=False)

    # Initialize an array to hold filtered data
    filtered_data = xr.DataArray(
        np.zeros_like(ecp_data), coords=ecp_data.coords, dims=ecp_data.dims
    )

    # Apply the filter to each channel
    for channel in ecp_data.channel_id:
        filtered_data.loc[channel, :] = signal.filtfilt(
            b, a, ecp_data.sel(channel_id=channel).values
        )

    # Downsample the filtered data if a downsample frequency is provided
    if downsample_freq is not None:
        downsample_factor = int(fs / downsample_freq)
        filtered_data = filtered_data.isel(time=slice(None, None, downsample_factor))
        # Update the sampling frequency attribute
        filtered_data.attrs["fs"] = downsample_freq

    return filtered_data

bmtool.analysis.lfp.slice_time_series(data, time_ranges)

Slice the xarray DataArray based on provided time ranges. Can be used to get LFP during certain stimulus times

Parameters:

data : xr.DataArray The input xarray DataArray containing time-series data. time_ranges : tuple or list of tuples One or more tuples representing the (start, stop) time points for slicing. For example: (start, stop) or [(start1, stop1), (start2, stop2)]

Returns:

xr.DataArray A new xarray DataArray containing the concatenated slices.

Source code in bmtool/analysis/lfp.py

def slice_time_series(data: xr.DataArray, time_ranges: tuple) -> xr.DataArray:
    """
    Slice the xarray DataArray based on provided time ranges.
    Can be used to get LFP during certain stimulus times

    Parameters:
    ----------
    data : xr.DataArray
        The input xarray DataArray containing time-series data.
    time_ranges : tuple or list of tuples
        One or more tuples representing the (start, stop) time points for slicing.
        For example: (start, stop) or [(start1, stop1), (start2, stop2)]

    Returns:
    -------
    xr.DataArray
        A new xarray DataArray containing the concatenated slices.
    """
    # Ensure time_ranges is a list of tuples
    if isinstance(time_ranges, tuple) and len(time_ranges) == 2:
        time_ranges = [time_ranges]

    # List to hold sliced data
    slices = []

    # Slice the data for each time range
    for start, stop in time_ranges:
        sliced_data = data.sel(time=slice(start, stop))
        slices.append(sliced_data)

    # Concatenate all slices along the time dimension if more than one slice
    if len(slices) > 1:
        return xr.concat(slices, dim="time")
    else:
        return slices[0]

bmtool.analysis.lfp.fit_fooof(f, pxx, aperiodic_mode='fixed', dB_threshold=3.0, max_n_peaks=10, freq_range=None, peak_width_limits=None, report=False, plot=False, plt_log=False, plt_range=None, figsize=None, title=None)

Fit a FOOOF model to power spectral density data.

Parameters:

f : array-like Frequencies corresponding to the power spectral density data. pxx : array-like Power spectral density data to fit. aperiodic_mode : str, optional The mode for fitting aperiodic components ('fixed' or 'knee', default is 'fixed'). dB_threshold : float, optional Minimum peak height in dB (default is 3). max_n_peaks : int, optional Maximum number of peaks to fit (default is 10). freq_range : tuple, optional Frequency range to fit (default is None, which uses the full range). peak_width_limits : tuple, optional Limits on the width of peaks (default is None). report : bool, optional If True, will print fitting results (default is False). plot : bool, optional If True, will plot the fitting results (default is False). plt_log : bool, optional If True, use a logarithmic scale for the y-axis in plots (default is False). plt_range : tuple, optional Range for plotting (default is None). figsize : tuple, optional Size of the figure (default is None). title : str, optional Title for the plot (default is None).

Returns:

tuple A tuple containing the fitting results and the FOOOF model object.

Source code in bmtool/analysis/lfp.py

def fit_fooof(
    f: np.ndarray,
    pxx: np.ndarray,
    aperiodic_mode: str = "fixed",
    dB_threshold: float = 3.0,
    max_n_peaks: int = 10,
    freq_range: tuple = None,
    peak_width_limits: tuple = None,
    report: bool = False,
    plot: bool = False,
    plt_log: bool = False,
    plt_range: tuple = None,
    figsize: tuple = None,
    title: str = None,
) -> tuple:
    """
    Fit a FOOOF model to power spectral density data.

    Parameters:
    ----------
    f : array-like
        Frequencies corresponding to the power spectral density data.
    pxx : array-like
        Power spectral density data to fit.
    aperiodic_mode : str, optional
        The mode for fitting aperiodic components ('fixed' or 'knee', default is 'fixed').
    dB_threshold : float, optional
        Minimum peak height in dB (default is 3).
    max_n_peaks : int, optional
        Maximum number of peaks to fit (default is 10).
    freq_range : tuple, optional
        Frequency range to fit (default is None, which uses the full range).
    peak_width_limits : tuple, optional
        Limits on the width of peaks (default is None).
    report : bool, optional
        If True, will print fitting results (default is False).
    plot : bool, optional
        If True, will plot the fitting results (default is False).
    plt_log : bool, optional
        If True, use a logarithmic scale for the y-axis in plots (default is False).
    plt_range : tuple, optional
        Range for plotting (default is None).
    figsize : tuple, optional
        Size of the figure (default is None).
    title : str, optional
        Title for the plot (default is None).

    Returns:
    -------
    tuple
        A tuple containing the fitting results and the FOOOF model object.
    """
    if aperiodic_mode != "knee":
        aperiodic_mode = "fixed"

    def set_range(x, upper=f[-1]):
        x = np.array(upper) if x is None else np.array(x)
        return [f[2], x.item()] if x.size == 1 else x.tolist()

    freq_range = set_range(freq_range)
    peak_width_limits = set_range(peak_width_limits, np.inf)

    # Initialize a FOOOF object
    fm = FOOOF(
        peak_width_limits=peak_width_limits,
        min_peak_height=dB_threshold / 10,
        peak_threshold=0.0,
        max_n_peaks=max_n_peaks,
        aperiodic_mode=aperiodic_mode,
    )

    # Fit the model
    try:
        fm.fit(f, pxx, freq_range)
    except Exception as e:
        fl = np.linspace(f[0], f[-1], int((f[-1] - f[0]) / np.min(np.diff(f))) + 1)
        fm.fit(fl, np.interp(fl, f, pxx), freq_range)

    results = fm.get_results()

    if report:
        fm.print_results()
        if aperiodic_mode == "knee":
            ap_params = results.aperiodic_params
            if ap_params[1] <= 0:
                print(
                    "Negative value of knee parameter occurred. Suggestion: Fit without knee parameter."
                )
            knee_freq = np.abs(ap_params[1]) ** (1 / ap_params[2])
            print(f"Knee location: {knee_freq:.2f} Hz")

    if plot:
        plt_range = set_range(plt_range)
        fm.plot(ax=plt.gca(), plt_log=plt_log)
        plt.xlim(np.log10(plt_range) if plt_log else plt_range)
        # plt.ylim(-8, -5.5)
        if figsize:
            plt.gcf().set_size_inches(figsize)
        if title:
            plt.title(title)
        if is_notebook():
            pass
        else:
            plt.show()

    return results, fm

bmtool.analysis.lfp.generate_resd_from_fooof(fooof_model)

Generate residuals from a fitted FOOOF model.

Parameters:

fooof_model : FOOOF A fitted FOOOF model object.

Returns:

tuple A tuple containing the residual power spectral density and the aperiodic fit.

Source code in bmtool/analysis/lfp.py

def generate_resd_from_fooof(fooof_model: FOOOF) -> tuple:
    """
    Generate residuals from a fitted FOOOF model.

    Parameters:
    ----------
    fooof_model : FOOOF
        A fitted FOOOF model object.

    Returns:
    -------
    tuple
        A tuple containing the residual power spectral density and the aperiodic fit.
    """
    results = fooof_model.get_results()
    full_fit, _, ap_fit = gen_model(
        fooof_model.freqs[1:],
        results.aperiodic_params,
        results.gaussian_params,
        return_components=True,
    )

    full_fit, ap_fit = 10**full_fit, 10**ap_fit  # Convert back from log
    res_psd = np.insert(
        (10 ** fooof_model.power_spectrum[1:]) - ap_fit, 0, 0.0
    )  # Convert back from log
    res_fit = np.insert(full_fit - ap_fit, 0, 0.0)
    ap_fit = np.insert(ap_fit, 0, 0.0)

    return res_psd, ap_fit

bmtool.analysis.lfp.calculate_SNR(fooof_model, freq_band)

Calculate the signal-to-noise ratio (SNR) from a fitted FOOOF model.

Parameters:

fooof_model : FOOOF A fitted FOOOF model object. freq_band : tuple Frequency band (min, max) for SNR calculation.

Returns:

float The calculated SNR for the specified frequency band.

Source code in bmtool/analysis/lfp.py

def calculate_SNR(fooof_model: FOOOF, freq_band: tuple) -> float:
    """
    Calculate the signal-to-noise ratio (SNR) from a fitted FOOOF model.

    Parameters:
    ----------
    fooof_model : FOOOF
        A fitted FOOOF model object.
    freq_band : tuple
        Frequency band (min, max) for SNR calculation.

    Returns:
    -------
    float
        The calculated SNR for the specified frequency band.
    """
    periodic, ap = generate_resd_from_fooof(fooof_model)
    freq = fooof_model.freqs  # Get frequencies from model
    indices = (freq >= freq_band[0]) & (freq <= freq_band[1])  # Get only the band we care about
    band_periodic = periodic[indices]  # Filter based on band
    band_ap = ap[indices]  # Filter
    band_freq = freq[indices]  # Another filter
    periodic_power = np.trapz(band_periodic, band_freq)  # Integrate periodic power
    ap_power = np.trapz(band_ap, band_freq)  # Integrate aperiodic power
    normalized_power = periodic_power / ap_power  # Compute the SNR
    return normalized_power

bmtool.analysis.lfp.wavelet_filter(x, freq, fs, bandwidth=1.0, axis=-1, show_passband=False)

Compute the Continuous Wavelet Transform (CWT) for a specified frequency using a complex Morlet wavelet.

Parameters:

Name	Type	Description	Default
x	ndarray	Input signal	required
freq	float	Target frequency for the wavelet filter	required
fs	float	Sampling frequency of the signal	required
bandwidth	float	Bandwidth parameter of the wavelet filter (default is 1.0)	1.0
axis	int	Axis along which to compute the CWT (default is -1)	-1
show_passband	bool	If True, print the passband of the wavelet filter (default is False)	False

Returns:

Type	Description
ndarray	Continuous Wavelet Transform of the input signal

Source code in bmtool/analysis/lfp.py

def wavelet_filter(
    x: np.ndarray,
    freq: float,
    fs: float,
    bandwidth: float = 1.0,
    axis: int = -1,
    show_passband: bool = False,
) -> np.ndarray:
    """
    Compute the Continuous Wavelet Transform (CWT) for a specified frequency using a complex Morlet wavelet.

    Parameters
    ----------
    x : np.ndarray
        Input signal
    freq : float
        Target frequency for the wavelet filter
    fs : float
        Sampling frequency of the signal
    bandwidth : float, optional
        Bandwidth parameter of the wavelet filter (default is 1.0)
    axis : int, optional
        Axis along which to compute the CWT (default is -1)
    show_passband : bool, optional
        If True, print the passband of the wavelet filter (default is False)

    Returns
    -------
    np.ndarray
        Continuous Wavelet Transform of the input signal
    """
    if show_passband:
        lower_bound, upper_bound, passband_width = calculate_wavelet_passband(
            freq, bandwidth, threshold=0.3
        )  # kinda made up threshold gives the rough idea
        print(f"Wavelet filter at {freq:.1f} Hz Bandwidth: {bandwidth:.1f} Hz:")
        print(
            f"  Passband: {lower_bound:.1f} - {upper_bound:.1f} Hz (width: {passband_width:.1f} Hz)"
        )
    wavelet = "cmor" + str(2 * bandwidth**2) + "-1.0"
    scale = pywt.scale2frequency(wavelet, 1) * fs / freq
    x_a = pywt.cwt(x, [scale], wavelet=wavelet, axis=axis)[0][0]
    return x_a

bmtool.analysis.lfp.butter_bandpass_filter(data, lowcut, highcut, fs, order=5, axis=-1)

Apply a Butterworth bandpass filter to the input data.

Source code in bmtool/analysis/lfp.py

def butter_bandpass_filter(
    data: np.ndarray, lowcut: float, highcut: float, fs: float, order: int = 5, axis: int = -1
) -> np.ndarray:
    """
    Apply a Butterworth bandpass filter to the input data.
    """
    sos = signal.butter(order, [lowcut, highcut], fs=fs, btype="band", output="sos")
    x_a = signal.sosfiltfilt(sos, data, axis=axis)
    return x_a

bmtool.analysis.lfp.cwt_spectrogram(x, fs, nNotes=6, nOctaves=np.inf, freq_range=(0, np.inf), bandwidth=1.0, axis=-1, detrend=False, normalize=False)

Calculate spectrogram using continuous wavelet transform

Source code in bmtool/analysis/lfp.py

def cwt_spectrogram(
    x,
    fs,
    nNotes=6,
    nOctaves=np.inf,
    freq_range=(0, np.inf),
    bandwidth=1.0,
    axis=-1,
    detrend=False,
    normalize=False,
):
    """Calculate spectrogram using continuous wavelet transform"""
    x = np.asarray(x)
    N = x.shape[axis]
    times = np.arange(N) / fs
    # detrend and normalize
    if detrend:
        x = signal.detrend(x, axis=axis, type="linear")
    if normalize:
        x = x / x.std()
    # Define some parameters of our wavelet analysis.
    # range of scales (in time) that makes sense
    # min = 2 (Nyquist frequency)
    # max = np.floor(N/2)
    nOctaves = min(nOctaves, np.log2(2 * np.floor(N / 2)))
    scales = 2 ** np.arange(1, nOctaves, 1 / nNotes)
    # cwt and the frequencies used.
    # Use the complex morelet with bw=2*bandwidth^2 and center frequency of 1.0
    # bandwidth is sigma of the gaussian envelope
    wavelet = "cmor" + str(2 * bandwidth**2) + "-1.0"
    frequencies = pywt.scale2frequency(wavelet, scales) * fs
    scales = scales[(frequencies >= freq_range[0]) & (frequencies <= freq_range[1])]
    coef, frequencies = pywt.cwt(
        x, scales[::-1], wavelet=wavelet, sampling_period=1 / fs, axis=axis
    )
    power = np.real(coef * np.conj(coef))  # equivalent to power = np.abs(coef)**2
    # cone of influence in terms of wavelength
    coi = N / 2 - np.abs(np.arange(N) - (N - 1) / 2)
    # cone of influence in terms of frequency
    coif = COI_FREQ * fs / coi
    return power, times, frequencies, coif

bmtool.analysis.lfp.cwt_spectrogram_xarray(x, fs, time=None, axis=-1, downsample_fs=None, channel_coords=None, **cwt_kwargs)

Calculate spectrogram using continuous wavelet transform and return an xarray.Dataset x: input array fs: sampling frequency (Hz) axis: dimension index of time axis in x downsample_fs: downsample to the frequency if specified channel_coords: dictionary of {coordinate name: index} for channels cwt_kwargs: keyword arguments for cwt_spectrogram()

Source code in bmtool/analysis/lfp.py

def cwt_spectrogram_xarray(
    x, fs, time=None, axis=-1, downsample_fs=None, channel_coords=None, **cwt_kwargs
):
    """Calculate spectrogram using continuous wavelet transform and return an xarray.Dataset
    x: input array
    fs: sampling frequency (Hz)
    axis: dimension index of time axis in x
    downsample_fs: downsample to the frequency if specified
    channel_coords: dictionary of {coordinate name: index} for channels
    cwt_kwargs: keyword arguments for cwt_spectrogram()
    """
    x = np.asarray(x)
    T = x.shape[axis]  # number of time points
    t = np.arange(T) / fs if time is None else np.asarray(time)
    if downsample_fs is None or downsample_fs >= fs:
        downsample_fs = fs
        downsampled = x
    else:
        num = int(T * downsample_fs / fs)
        downsample_fs = num / T * fs
        downsampled, t = signal.resample(x, num=num, t=t, axis=axis)
    downsampled = np.moveaxis(downsampled, axis, -1)
    sxx, _, f, coif = cwt_spectrogram(downsampled, downsample_fs, **cwt_kwargs)
    sxx = np.moveaxis(sxx, 0, -2)  # shape (... , freq, time)
    if channel_coords is None:
        channel_coords = {f"dim_{i:d}": range(d) for i, d in enumerate(sxx.shape[:-2])}
    sxx = xr.DataArray(sxx, coords={**channel_coords, "frequency": f, "time": t}).to_dataset(
        name="PSD"
    )
    sxx.update(dict(cone_of_influence_frequency=xr.DataArray(coif, coords={"time": t})))
    return sxx