Skip to content

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 
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
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 
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
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 
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
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 
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
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 
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
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.get_fooof_freq_band(gaussian_params, freq_range, width_limit=(0.0, np.inf), top_n_peaks=1, bandwidth_n_sigma=1.5)

Get frequency band of the top N peaks in the FOOOF results within a given band of interest.

Parameters:

Name Type Description Default
gaussian_params NDArray[float]

Gaussian parameters from FOOOF results, FOOOFResults.gaussian_params.

 required
freq_range tuple[float, float]

Frequency band of interest

 required
width_limit tuple[float, float]

Width limit of the peaks in terms of the standard deviation of the Gaussian parameters.

 (0.0, inf)
top_n_peaks int

Number of top peaks to include in the band.

 1
bandwidth_n_sigma float

Multiplier of sigma of the Gaussian parameters to define the bandwidth of the peak.

 1.5

Returns:

Name Type Description
band tuple[float, float]

Combined frequency band of the top N peaks within the given band of interest. If no peaks are found within the given band of interest, return (np.nan, np.nan).
peak_inds array_like of bool

Boolean array of the peaks within the given band of interest.
Source code in bmtool/analysis/lfp.py 
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
def get_fooof_freq_band(
    gaussian_params,
    freq_range,
    width_limit=(0., np.inf),
    top_n_peaks=1,
    bandwidth_n_sigma=1.5
):
    """Get frequency band of the top N peaks in the FOOOF results within a given band of interest.

    Parameters
    ----------
    gaussian_params : NDArray[float]
        Gaussian parameters from FOOOF results, FOOOFResults.gaussian_params.
    freq_range : tuple[float, float]
        Frequency band of interest
    width_limit : tuple[float, float]
        Width limit of the peaks in terms of the standard deviation of the Gaussian parameters.
    top_n_peaks : int
        Number of top peaks to include in the band.
    bandwidth_n_sigma : float
        Multiplier of sigma of the Gaussian parameters to define the bandwidth of the peak.

    Returns
    -------
    band : tuple[float, float]
        Combined frequency band of the top N peaks within the given band of interest.
        If no peaks are found within the given band of interest, return (np.nan, np.nan).
    peak_inds : array_like of bool
        Boolean array of the peaks within the given band of interest.
    """
    # find peaks within the given band of interest
    peak_inds = (gaussian_params[:, 0] >= freq_range[0]) & (gaussian_params[:, 0] <= freq_range[1])
    peak_inds = peak_inds & (gaussian_params[:, 2] >= width_limit[0]) & (gaussian_params[:, 2] <= width_limit[1])
    peak_inds = np.nonzero(peak_inds)[0]
    if peak_inds.size == 0:
        return (np.nan, np.nan), peak_inds

    # find top peaks by height
    top_n_peaks = max(top_n_peaks, 1)  # at least one peak
    idx = np.argsort(gaussian_params[peak_inds, 1])[::-1][:top_n_peaks]
    peak_inds = peak_inds[idx]

    # get the combined frequency band of the top peaks
    band_peaks = gaussian_params[peak_inds, :]
    band_widths = bandwidth_n_sigma * band_peaks[:, 2]  # one-sided bandwidth
    band_freqs = np.fmax(band_peaks[:, [0]] + np.outer(band_widths, [-1, 1]), 0.)
    band = (band_freqs[:, 0].min(), band_freqs[:, 1].max())  # combined frequency band
    return band, peak_inds

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 
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
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 
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
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 
450
451
452
453
454
455
456
457
458
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.get_lfp_power(lfp_data, freq_of_interest, fs, filter_method='wavelet', lowcut=None, highcut=None, bandwidth=1.0)

Compute the power of the raw LFP signal in a specified frequency band, preserving xarray structure if input is xarray.

Parameters:

Name Type Description Default
lfp_data ndarray or DataArray

Raw local field potential (LFP) time series data

 required
freq_of_interest float

Center frequency (Hz) for wavelet filtering method

 required
fs float

Sampling frequency (Hz) of the input data

 required
filter_method str

Filtering method to use, either 'wavelet' or 'butter' (default: 'wavelet')

 'wavelet'
lowcut float

Lower frequency bound (Hz) for butterworth bandpass filter, required if filter_method='butter'

 None
highcut float

Upper frequency bound (Hz) for butterworth bandpass filter, required if filter_method='butter'

 None
bandwidth float

Bandwidth parameter for wavelet filter when method='wavelet' (default: 1.0)

 1.0

Returns:

Type Description
ndarray or DataArray

Power of the filtered signal (magnitude squared) with same structure as input
Notes
  • The 'wavelet' method uses a complex Morlet wavelet centered at the specified frequency
  • The 'butter' method uses a Butterworth bandpass filter with the specified cutoff frequencies
  • When using the 'butter' method, both lowcut and highcut must be provided
  • If input is an xarray DataArray, the output will preserve the same structure with coordinates
Source code in bmtool/analysis/lfp.py 
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
def get_lfp_power(
    lfp_data,
    freq_of_interest: float,
    fs: float,
    filter_method: str = "wavelet",
    lowcut: float = None,
    highcut: float = None,
    bandwidth: float = 1.0,
):
    """
    Compute the power of the raw LFP signal in a specified frequency band,
    preserving xarray structure if input is xarray.

    Parameters
    ----------
    lfp_data : np.ndarray or xr.DataArray
        Raw local field potential (LFP) time series data
    freq_of_interest : float
        Center frequency (Hz) for wavelet filtering method
    fs : float
        Sampling frequency (Hz) of the input data
    filter_method : str, optional
        Filtering method to use, either 'wavelet' or 'butter' (default: 'wavelet')
    lowcut : float, optional
        Lower frequency bound (Hz) for butterworth bandpass filter, required if filter_method='butter'
    highcut : float, optional
        Upper frequency bound (Hz) for butterworth bandpass filter, required if filter_method='butter'
    bandwidth : float, optional
        Bandwidth parameter for wavelet filter when method='wavelet' (default: 1.0)

    Returns
    -------
    np.ndarray or xr.DataArray
        Power of the filtered signal (magnitude squared) with same structure as input

    Notes
    -----
    - The 'wavelet' method uses a complex Morlet wavelet centered at the specified frequency
    - The 'butter' method uses a Butterworth bandpass filter with the specified cutoff frequencies
    - When using the 'butter' method, both lowcut and highcut must be provided
    - If input is an xarray DataArray, the output will preserve the same structure with coordinates
    """
    import xarray as xr

    # Check if input is xarray
    is_xarray = isinstance(lfp_data, xr.DataArray)

    if is_xarray:
        # Get the raw data from xarray
        raw_data = lfp_data.values
        # Check if 'fs' attribute exists in the xarray and override if necessary
        if "fs" in lfp_data.attrs and fs is None:
            fs = lfp_data.attrs["fs"]
    else:
        raw_data = lfp_data

    if filter_method == "wavelet":
        filtered_signal = wavelet_filter(raw_data, freq_of_interest, fs, bandwidth)
    elif filter_method == "butter":
        if lowcut is None or highcut is None:
            raise ValueError(
                "Both lowcut and highcut must be specified when using 'butter' method."
            )
        filtered_signal = butter_bandpass_filter(raw_data, lowcut, highcut, fs)
    else:
        raise ValueError("Invalid method. Choose 'wavelet' or 'butter'.")

    # Calculate power (magnitude squared of filtered signal)
    power = np.abs(filtered_signal) ** 2

    # If the input was an xarray, return an xarray with the same coordinates
    if is_xarray:
        power_xarray = xr.DataArray(
            power,
            coords=lfp_data.coords,
            dims=lfp_data.dims,
            attrs={
                **lfp_data.attrs,
                "filter_method": filter_method,
                "frequency_of_interest": freq_of_interest,
                "bandwidth": bandwidth,
                "lowcut": lowcut,
                "highcut": highcut,
                "power_type": "magnitude_squared",
            },
        )
        return power_xarray

    return power

bmtool.analysis.lfp.get_lfp_phase(lfp_data, freq_of_interest, fs, filter_method='wavelet', lowcut=None, highcut=None, bandwidth=1.0)

Calculate the phase of the filtered signal, preserving xarray structure if input is xarray.

Parameters:

Name Type Description Default
lfp_data ndarray or DataArray

Input LFP data

 required
freq_of_interest float

Frequency of interest (Hz)

 required
fs float

Sampling frequency (Hz)

 required
filter_method str

Method for filtering the signal ('wavelet' or 'butter')

 'wavelet'
bandwidth float

Bandwidth parameter for wavelet filter when method='wavelet' (default: 1.0)

 1.0
lowcut float

Low cutoff frequency for Butterworth filter when method='butter'

 None
highcut float

High cutoff frequency for Butterworth filter when method='butter'

 None

Returns:

Type Description
ndarray or DataArray

Phase of the filtered signal with same structure as input
Notes
  • The 'wavelet' method uses a complex Morlet wavelet centered at the specified frequency
  • The 'butter' method uses a Butterworth bandpass filter with the specified cutoff frequencies followed by Hilbert transform to extract the phase
  • When using the 'butter' method, both lowcut and highcut must be provided
  • If input is an xarray DataArray, the output will preserve the same structure with coordinates
Source code in bmtool/analysis/lfp.py 
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
def get_lfp_phase(
    lfp_data,
    freq_of_interest: float,
    fs: float,
    filter_method: str = "wavelet",
    lowcut: float = None,
    highcut: float = None,
    bandwidth: float = 1.0,
) -> np.ndarray:
    """
    Calculate the phase of the filtered signal, preserving xarray structure if input is xarray.

    Parameters
    ----------
    lfp_data : np.ndarray or xr.DataArray
        Input LFP data
    freq_of_interest : float
        Frequency of interest (Hz)
    fs : float
        Sampling frequency (Hz)
    filter_method : str, optional
        Method for filtering the signal ('wavelet' or 'butter')
    bandwidth : float, optional
        Bandwidth parameter for wavelet filter when method='wavelet' (default: 1.0)
    lowcut : float, optional
        Low cutoff frequency for Butterworth filter when method='butter'
    highcut : float, optional
        High cutoff frequency for Butterworth filter when method='butter'

    Returns
    -------
    np.ndarray or xr.DataArray
        Phase of the filtered signal with same structure as input

    Notes
    -----
    - The 'wavelet' method uses a complex Morlet wavelet centered at the specified frequency
    - The 'butter' method uses a Butterworth bandpass filter with the specified cutoff frequencies
      followed by Hilbert transform to extract the phase
    - When using the 'butter' method, both lowcut and highcut must be provided
    - If input is an xarray DataArray, the output will preserve the same structure with coordinates
    """
    import xarray as xr

    # Check if input is xarray
    is_xarray = isinstance(lfp_data, xr.DataArray)

    if is_xarray:
        # Get the raw data from xarray
        raw_data = lfp_data.values
        # Check if 'fs' attribute exists in the xarray and override if necessary
        if "fs" in lfp_data.attrs and fs is None:
            fs = lfp_data.attrs["fs"]
    else:
        raw_data = lfp_data

    if filter_method == "wavelet":
        if freq_of_interest is None:
            raise ValueError("freq_of_interest must be provided for the wavelet method.")
        # Wavelet filter returns complex values directly
        filtered_signal = wavelet_filter(raw_data, freq_of_interest, fs, bandwidth)
        # Phase is the angle of the complex signal
        phase = np.angle(filtered_signal)
    elif filter_method == "butter":
        if lowcut is None or highcut is None:
            raise ValueError(
                "Both lowcut and highcut must be specified when using 'butter' method."
            )
        # Butterworth filter returns real values
        filtered_signal = butter_bandpass_filter(raw_data, lowcut, highcut, fs)
        # Apply Hilbert transform to get analytic signal (complex)
        analytic_signal = signal.hilbert(filtered_signal)
        # Phase is the angle of the analytic signal
        phase = np.angle(analytic_signal)
    else:
        raise ValueError(f"Invalid method {filter_method}. Choose 'wavelet' or 'butter'.")

    # If the input was an xarray, return an xarray with the same coordinates
    if is_xarray:
        phase_xarray = xr.DataArray(
            phase,
            coords=lfp_data.coords,
            dims=lfp_data.dims,
            attrs={
                **lfp_data.attrs,
                "filter_method": filter_method,
                "freq_of_interest": freq_of_interest,
                "bandwidth": bandwidth,
                "lowcut": lowcut,
                "highcut": highcut,
            },
        )
        return phase_xarray

    return phase

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 
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
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 
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
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

bmtool.analysis.lfp.get_windowed_data(x, windows, win_grp_idx, dim='time', win_dim='cycle', win_coord=None, grp_dim='unique_cycle')

Apply functions of windowing to data x: DataArray windows: windows for windowed_xarray win_grp_idx: win_grp_idx for group_windows dim: dimension along which to divide win_dim: dimension for different windows win_coord: pandas Index object of win_dim coordinates grp_dim: dimension along which to stack average of window groups. If None or empty or False, do not calculate average. Return: data returned by three functions, windowed_xarray, group_windows, average_group_windows

Source code in bmtool/analysis/lfp.py 
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
def get_windowed_data(
    x, windows, win_grp_idx, dim="time", win_dim="cycle", win_coord=None, grp_dim="unique_cycle"
):
    """Apply functions of windowing to data
    x: DataArray
    windows: `windows` for `windowed_xarray`
    win_grp_idx: `win_grp_idx` for `group_windows`
    dim: dimension along which to divide
    win_dim: dimension for different windows
    win_coord: pandas Index object of `win_dim` coordinates
    grp_dim: dimension along which to stack average of window groups.
        If None or empty or False, do not calculate average.
    Return: data returned by three functions,
        `windowed_xarray`, `group_windows`, `average_group_windows`
    """
    x_win = windowed_xarray(x, windows, dim=dim, new_coord_name=win_dim, new_coord=win_coord)
    x_win_onff = group_windows(x_win, win_grp_idx, win_dim=win_dim)
    if grp_dim:
        x_win_avg = [average_group_windows(x, win_dim=win_dim, grp_dim=grp_dim) for x in x_win_onff]
    else:
        x_win_avg = None
    return x_win, x_win_onff, x_win_avg