Skip to content

Analysis API Reference

This page provides API reference documentation for the Analysis module.

The Analysis module provides tools for processing and analyzing simulation results from BMTK models, including spike data and LFP/ECP data.

Spikes

The spikes module provides functions for loading and analyzing spike data from simulations.

bmtool.analysis.spikes.load_spikes_to_df(spike_file, network_name, sort=True, config=None, groupby='pop_name')

Load spike data from an HDF5 file into a pandas DataFrame.

Parameters:

Name Type Description Default
spike_file str

Path to the HDF5 file containing spike data

 required
network_name str

The name of the network within the HDF5 file from which to load spike data

 required
sort bool

Whether to sort the DataFrame by 'timestamps' (default: True)

 True
config str

Path to configuration file to label the cell type of each spike (default: None)

 None
groupby Union[str, List[str]]

The column(s) to group by (default: 'pop_name')

 'pop_name'

Returns:

Type Description
DataFrame

A pandas DataFrame containing 'node_ids' and 'timestamps' columns from the spike data, with additional columns if a config file is provided

Examples:

>>> df = load_spikes_to_df("spikes.h5", "cortex")
>>> df = load_spikes_to_df("spikes.h5", "cortex", config="config.json", groupby=["pop_name", "model_type"])
Source code in bmtool/analysis/spikes.py 
14
15
16
17
18
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
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
def load_spikes_to_df(spike_file: str, network_name: str, sort: bool = True, config: str = None, groupby: Union[str, List[str]] = 'pop_name') -> pd.DataFrame:
    """
    Load spike data from an HDF5 file into a pandas DataFrame.

    Parameters
    ----------
    spike_file : str
        Path to the HDF5 file containing spike data
    network_name : str
        The name of the network within the HDF5 file from which to load spike data
    sort : bool, optional
        Whether to sort the DataFrame by 'timestamps' (default: True)
    config : str, optional
        Path to configuration file to label the cell type of each spike (default: None)
    groupby : Union[str, List[str]], optional
        The column(s) to group by (default: 'pop_name')

    Returns
    -------
    pd.DataFrame
        A pandas DataFrame containing 'node_ids' and 'timestamps' columns from the spike data,
        with additional columns if a config file is provided

    Examples
    --------
    >>> df = load_spikes_to_df("spikes.h5", "cortex")
    >>> df = load_spikes_to_df("spikes.h5", "cortex", config="config.json", groupby=["pop_name", "model_type"])
    """
    with h5py.File(spike_file) as f:
        spikes_df = pd.DataFrame({
            'node_ids': f['spikes'][network_name]['node_ids'],
            'timestamps': f['spikes'][network_name]['timestamps']
        })

        if sort:
            spikes_df.sort_values(by='timestamps', inplace=True, ignore_index=True)

        if config:
            nodes = load_nodes_from_config(config)
            nodes = nodes[network_name]

            # Convert single string to a list for uniform handling
            if isinstance(groupby, str):
                groupby = [groupby]

            # Ensure all requested columns exist
            missing_cols = [col for col in groupby if col not in nodes.columns]
            if missing_cols:
                raise KeyError(f"Columns {missing_cols} not found in nodes DataFrame.")

            spikes_df = spikes_df.merge(nodes[groupby], left_on='node_ids', right_index=True, how='left')

    return spikes_df

bmtool.analysis.spikes.compute_firing_rate_stats(df, groupby='pop_name', start_time=None, stop_time=None)

Computes the firing rates of individual nodes and the mean and standard deviation of firing rates per group.

Args: df (pd.DataFrame): Dataframe containing spike timestamps and node IDs. groupby (str or list of str, optional): Column(s) to group by (e.g., 'pop_name' or ['pop_name', 'layer']). start_time (float, optional): Start time for the analysis window. Defaults to the minimum timestamp in the data. stop_time (float, optional): Stop time for the analysis window. Defaults to the maximum timestamp in the data.

Returns: Tuple[pd.DataFrame, pd.DataFrame]: - The first DataFrame (pop_stats) contains the mean and standard deviation of firing rates per group. - The second DataFrame (individual_stats) contains the firing rate of each individual node.

Source code in bmtool/analysis/spikes.py 
 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
 99
100
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
136
def compute_firing_rate_stats(df: pd.DataFrame, groupby: Union[str, List[str]] = "pop_name", start_time: float = None, stop_time: float = None) -> Tuple[pd.DataFrame, pd.DataFrame]:
    """
    Computes the firing rates of individual nodes and the mean and standard deviation of firing rates per group.

    Args:
        df (pd.DataFrame): Dataframe containing spike timestamps and node IDs.
        groupby (str or list of str, optional): Column(s) to group by (e.g., 'pop_name' or ['pop_name', 'layer']).
        start_time (float, optional): Start time for the analysis window. Defaults to the minimum timestamp in the data.
        stop_time (float, optional): Stop time for the analysis window. Defaults to the maximum timestamp in the data.

    Returns:
        Tuple[pd.DataFrame, pd.DataFrame]: 
            - The first DataFrame (`pop_stats`) contains the mean and standard deviation of firing rates per group.
            - The second DataFrame (`individual_stats`) contains the firing rate of each individual node.
    """

    # Ensure groupby is a list
    if isinstance(groupby, str):
        groupby = [groupby]

    # Ensure all columns exist in the dataframe
    for col in groupby:
        if col not in df.columns:
            raise ValueError(f"Column '{col}' not found in dataframe.")

    # Filter dataframe based on start/stop time
    if start_time is not None:
        df = df[df["timestamps"] >= start_time]
    if stop_time is not None:
        df = df[df["timestamps"] <= stop_time]

    # Compute total duration for firing rate calculation
    if start_time is None:
        min_time = df["timestamps"].min()
    else:
        min_time = start_time

    if stop_time is None: 
        max_time = df["timestamps"].max()
    else:
        max_time = stop_time

    duration = max_time - min_time  # Avoid division by zero

    if duration <= 0:
        raise ValueError("Invalid time window: Stop time must be greater than start time.")

    # Compute firing rate for each node
    import pandas as pd

    # Compute spike counts per node
    spike_counts = df["node_ids"].value_counts().reset_index()
    spike_counts.columns = ["node_ids", "spike_count"]  # Rename columns

    # Merge with original dataframe to get corresponding labels (e.g., 'pop_name')
    spike_counts = spike_counts.merge(df[["node_ids"] + groupby].drop_duplicates(), on="node_ids", how="left")

    # Compute firing rate
    spike_counts["firing_rate"] = spike_counts["spike_count"] / duration * 1000 # scale to Hz
    indivdual_stats = spike_counts

    # Compute mean and standard deviation per group
    pop_stats = spike_counts.groupby(groupby)["firing_rate"].agg(["mean", "std"]).reset_index()

    # Rename columns
    pop_stats.rename(columns={"mean": "firing_rate_mean", "std": "firing_rate_std"}, inplace=True)

    return pop_stats,indivdual_stats

bmtool.analysis.spikes._pop_spike_rate(spike_times, time=None, time_points=None, frequency=False)

Calculate the spike count or frequency histogram over specified time intervals.

Parameters:

Name Type Description Default
spike_times Union[ndarray, list]

Array or list of spike times in milliseconds

 required
time Optional[Tuple[float, float, float]]

Tuple specifying (start, stop, step) in milliseconds. Used to create evenly spaced time points if time_points is not provided. Default is None.

 None
time_points Optional[Union[ndarray, list]]

Array or list of specific time points for binning. If provided, time is ignored. Default is None.

 None
frequency bool

If True, returns spike frequency in Hz; otherwise, returns spike count. Default is False.

 False

Returns:

Type Description
ndarray

Array of spike counts or frequencies, depending on the frequency flag.

Raises:

Type Description
ValueError

If both time and time_points are None.
Source code in bmtool/analysis/spikes.py 
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
def _pop_spike_rate(spike_times: Union[np.ndarray, list], time: Optional[Tuple[float, float, float]] = None, 
                   time_points: Optional[Union[np.ndarray, list]] = None, frequency: bool = False) -> np.ndarray:
    """
    Calculate the spike count or frequency histogram over specified time intervals.

    Parameters
    ----------
    spike_times : Union[np.ndarray, list]
        Array or list of spike times in milliseconds
    time : Optional[Tuple[float, float, float]], optional
        Tuple specifying (start, stop, step) in milliseconds. Used to create evenly spaced time points 
        if `time_points` is not provided. Default is None.
    time_points : Optional[Union[np.ndarray, list]], optional
        Array or list of specific time points for binning. If provided, `time` is ignored. Default is None.
    frequency : bool, optional
        If True, returns spike frequency in Hz; otherwise, returns spike count. Default is False.

    Returns
    -------
    np.ndarray
        Array of spike counts or frequencies, depending on the `frequency` flag.

    Raises
    ------
    ValueError
        If both `time` and `time_points` are None.
    """
    if time_points is None:
        if time is None:
            raise ValueError("Either `time` or `time_points` must be provided.")
        time_points = np.arange(*time)
        dt = time[2]
    else:
        time_points = np.asarray(time_points).ravel()
        dt = (time_points[-1] - time_points[0]) / (time_points.size - 1)

    bins = np.append(time_points, time_points[-1] + dt)
    spike_rate, _ = np.histogram(np.asarray(spike_times), bins)

    if frequency:
        spike_rate = 1000 / dt * spike_rate

    return spike_rate

bmtool.analysis.spikes.get_population_spike_rate(spike_data, fs=400.0, t_start=0, t_stop=None, config=None, network_name=None, save=False, save_path=None, normalize=False)

Calculate the population spike rate for each population in the given spike data, with an option to normalize.

Parameters:

Name Type Description Default
spike_data DataFrame

A DataFrame containing spike data with columns 'pop_name', 'timestamps', and 'node_ids'

 required
fs float

Sampling frequency in Hz, which determines the time bin size for calculating the spike rate (default: 400.0)

 400.0
t_start float

Start time (in milliseconds) for spike rate calculation (default: 0)

 0
t_stop Optional[float]

Stop time (in milliseconds) for spike rate calculation. If None, defaults to the maximum timestamp in the data

 None
config Optional[str]

Path to a configuration file containing node information, used to determine the correct number of nodes per population. If None, node count is estimated from unique node spikes (default: None)

 None
network_name Optional[str]

Name of the network used in the configuration file, allowing selection of nodes for that network. Required if config is provided (default: None)

 None
save bool

Whether to save the calculated population spike rate to a file (default: False)

 False
save_path Optional[str]

Directory path where the file should be saved if save is True (default: None)

 None
normalize bool

Whether to normalize the spike rates for each population to a range of [0, 1] (default: False)

 False

Returns:

Type Description
Dict[str, ndarray]

A dictionary where keys are population names, and values are arrays representing the spike rate over time for each population. If normalize is True, each population's spike rate is scaled to [0, 1].

Raises:

Type Description
ValueError

If save is True but save_path is not provided.
Notes
  • If config is None, the function assumes all cells in each population have fired at least once; otherwise, the node count may be inaccurate.
  • If normalization is enabled, each population's spike rate is scaled using Min-Max normalization based on its own minimum and maximum values.
Source code in bmtool/analysis/spikes.py 
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
244
245
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
276
277
278
279
280
281
282
283
284
285
286
287
def get_population_spike_rate(spike_data: pd.DataFrame, fs: float = 400.0, t_start: float = 0, t_stop: Optional[float] = None, 
                              config: Optional[str] = None, network_name: Optional[str] = None,
                              save: bool = False, save_path: Optional[str] = None,
                              normalize: bool = False) -> Dict[str, np.ndarray]:
    """
    Calculate the population spike rate for each population in the given spike data, with an option to normalize.

    Parameters
    ----------
    spike_data : pd.DataFrame
        A DataFrame containing spike data with columns 'pop_name', 'timestamps', and 'node_ids'
    fs : float, optional
        Sampling frequency in Hz, which determines the time bin size for calculating the spike rate (default: 400.0)
    t_start : float, optional
        Start time (in milliseconds) for spike rate calculation (default: 0)
    t_stop : Optional[float], optional
        Stop time (in milliseconds) for spike rate calculation. If None, defaults to the maximum timestamp in the data
    config : Optional[str], optional
        Path to a configuration file containing node information, used to determine the correct number of nodes per population.
        If None, node count is estimated from unique node spikes (default: None)
    network_name : Optional[str], optional
        Name of the network used in the configuration file, allowing selection of nodes for that network.
        Required if `config` is provided (default: None)
    save : bool, optional
        Whether to save the calculated population spike rate to a file (default: False)
    save_path : Optional[str], optional
        Directory path where the file should be saved if `save` is True (default: None)
    normalize : bool, optional
        Whether to normalize the spike rates for each population to a range of [0, 1] (default: False)

    Returns
    -------
    Dict[str, np.ndarray]
        A dictionary where keys are population names, and values are arrays representing the spike rate over time for each population.
        If `normalize` is True, each population's spike rate is scaled to [0, 1].

    Raises
    ------
    ValueError
        If `save` is True but `save_path` is not provided.

    Notes
    -----
    - If `config` is None, the function assumes all cells in each population have fired at least once; otherwise, the node count may be inaccurate.
    - If normalization is enabled, each population's spike rate is scaled using Min-Max normalization based on its own minimum and maximum values.
    """
    pop_spikes = {}
    node_number = {}

    if config is None:
        print("Note: Node number is obtained by counting unique node spikes in the network.\nIf the network did not run for a sufficient duration, and not all cells fired, this count might be incorrect.")
        print("You can provide a config to calculate the correct amount of nodes!")

    if config:
        if not network_name:
            print("Grabbing first network; specify a network name to ensure correct node population is selected.")

    for pop_name in spike_data['pop_name'].unique():
        ps = spike_data[spike_data['pop_name'] == pop_name]

        if config:
            nodes = load_nodes_from_config(config)
            if network_name:
                nodes = nodes[network_name]
            else:
                nodes = list(nodes.values())[0] if nodes else {}
            nodes = nodes[nodes['pop_name'] == pop_name]
            node_number[pop_name] = nodes.index.nunique()
        else:
            node_number[pop_name] = ps['node_ids'].nunique()

        if t_stop is None:
            t_stop = spike_data['timestamps'].max()

        filtered_spikes = spike_data[
            (spike_data['pop_name'] == pop_name) & 
            (spike_data['timestamps'] > t_start) & 
            (spike_data['timestamps'] < t_stop)
        ]
        pop_spikes[pop_name] = filtered_spikes

    time = np.array([t_start, t_stop, 1000 / fs])
    pop_rspk = {p: _pop_spike_rate(spk['timestamps'], time) for p, spk in pop_spikes.items()}
    spike_rate = {p: fs / node_number[p] * pop_rspk[p] for p in pop_rspk}

    # Normalize each spike rate series if normalize=True
    if normalize:
        spike_rate = {p: (sr - sr.min()) / (sr.max() - sr.min()) for p, sr in spike_rate.items()}

    if save:
        if save_path is None:
            raise ValueError("save_path must be provided if save is True.")

        os.makedirs(save_path, exist_ok=True)

        save_file = os.path.join(save_path, 'spike_rate.h5')
        with h5py.File(save_file, 'w') as f:
            f.create_dataset('time', data=time)
            grp = f.create_group('populations')
            for p, rspk in spike_rate.items():
                pop_grp = grp.create_group(p)
                pop_grp.create_dataset('data', data=rspk)

    return spike_rate

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 
16
17
18
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
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 
49
50
51
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
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 
 93
 94
 95
 96
 97
 98
 99
100
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
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 
130
131
132
133
134
135
136
137
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
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., 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(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 
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
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.)  # Convert back from log
    res_fit = np.insert(full_fit - ap_fit, 0, 0.)
    ap_fit = np.insert(ap_fit, 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 
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
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 
325
326
327
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
354
355
356
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 
359
360
361
362
363
364
365
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 
550
551
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
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 
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
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

Entrainment Analysis

The entrainment module provides tools for analyzing the entrainment of spikes and lfp

bmtool.analysis.entrainment.calculate_signal_signal_plv(signal1, signal2, fs, freq_of_interest=None, filter_method='wavelet', lowcut=None, highcut=None, bandwidth=2.0)

Calculate Phase Locking Value (PLV) between two signals using wavelet or Hilbert method.

Parameters:

Name Type Description Default
signal1 ndarray

First input signal (1D array)

 required
signal2 ndarray

Second input signal (1D array, same length as signal1)

 required
fs float

Sampling frequency in Hz

 required
freq_of_interest float

Desired frequency for wavelet PLV calculation, required if filter_method='wavelet'

 None
filter_method str

Method to use for filtering, 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: 2.0)

 2.0

Returns:

Type Description
ndarray

Phase Locking Value (1D array)
Source code in bmtool/analysis/entrainment.py 
16
17
18
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
50
51
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
def calculate_signal_signal_plv(signal1: np.ndarray, signal2: np.ndarray, fs: float, freq_of_interest: float = None, 
                  filter_method: str = 'wavelet', lowcut: float = None, highcut: float = None, 
                  bandwidth: float = 2.0) -> np.ndarray:
    """
    Calculate Phase Locking Value (PLV) between two signals using wavelet or Hilbert method.

    Parameters
    ----------
    signal1 : np.ndarray
        First input signal (1D array)
    signal2 : np.ndarray
        Second input signal (1D array, same length as signal1)
    fs : float
        Sampling frequency in Hz
    freq_of_interest : float, optional
        Desired frequency for wavelet PLV calculation, required if filter_method='wavelet'
    filter_method : str, optional
        Method to use for filtering, 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: 2.0)

    Returns
    -------
    np.ndarray
        Phase Locking Value (1D array)
    """
    if len(signal1) != len(signal2):
        raise ValueError("Input signals must have the same length.")

    if filter_method == 'wavelet':
        if freq_of_interest is None:
            raise ValueError("freq_of_interest must be provided for the wavelet method.")

        # Apply CWT to both signals
        theta1 = wavelet_filter(x=signal1, freq=freq_of_interest, fs=fs, bandwidth=bandwidth)
        theta2 = wavelet_filter(x=signal2, freq=freq_of_interest, fs=fs, bandwidth=bandwidth)

    elif filter_method == 'butter':
        if lowcut is None or highcut is None:
            print("Lowcut and/or highcut were not defined, signal will not be filtered and will just take Hilbert transform for PLV calculation")

        if lowcut and highcut:
            # Bandpass filter and get the analytic signal using the Hilbert transform
            filtered_signal1 = butter_bandpass_filter(data=signal1, lowcut=lowcut, highcut=highcut, fs=fs)
            filtered_signal2 = butter_bandpass_filter(data=signal2, lowcut=lowcut, highcut=highcut, fs=fs)
            # Get phase using the Hilbert transform
            theta1 = signal.hilbert(filtered_signal1)
            theta2 = signal.hilbert(filtered_signal2)
        else:
            # Get phase using the Hilbert transform without filtering
            theta1 = signal.hilbert(signal1)
            theta2 = signal.hilbert(signal2)

    else:
        raise ValueError("Invalid method. Choose 'wavelet' or 'butter'.")

    # Calculate phase difference
    phase_diff = np.angle(theta1) - np.angle(theta2)

    # Calculate PLV from standard equation from Measuring phase synchrony in brain signals(1999)
    plv = np.abs(np.mean(np.exp(1j * phase_diff), axis=-1))

    return plv

bmtool.analysis.entrainment.calculate_spike_lfp_plv(spike_times=None, lfp_data=None, spike_fs=None, lfp_fs=None, filter_method='butter', freq_of_interest=None, lowcut=None, highcut=None, bandwidth=2.0, filtered_lfp_phase=None)

Calculate spike-lfp unbiased phase locking value

Parameters:

Name Type Description Default
spike_times ndarray

Array of spike times

 None
lfp_data ndarray

Local field potential time series data. Not required if filtered_lfp_phase is provided.

 None
spike_fs float

Sampling frequency in Hz of the spike times, only needed if spike times and LFP have different sampling rates

 None
lfp_fs float

Sampling frequency in Hz of the LFP data

 None
filter_method str

Method to use for filtering, either 'wavelet' or 'butter' (default: 'butter')

 'butter'
freq_of_interest float

Desired frequency for wavelet phase extraction, required if filter_method='wavelet'

 None
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: 2.0)

 2.0
filtered_lfp_phase ndarray

Pre-computed instantaneous phase of the filtered LFP. If provided, the function will skip the filtering step.

 None

Returns:

Type Description
float

Phase Locking Value (unbiased)
Source code in bmtool/analysis/entrainment.py 
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
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
136
137
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
def calculate_spike_lfp_plv(spike_times: np.ndarray = None, lfp_data: np.ndarray = None, spike_fs: float = None,
                   lfp_fs: float = None, filter_method: str = 'butter', freq_of_interest: float = None,
                   lowcut: float = None, highcut: float = None, bandwidth: float = 2.0,
                   filtered_lfp_phase: np.ndarray = None) -> float:
    """
    Calculate spike-lfp unbiased phase locking value 

    Parameters
    ----------
    spike_times : np.ndarray
        Array of spike times
    lfp_data : np.ndarray
        Local field potential time series data. Not required if filtered_lfp_phase is provided.
    spike_fs : float, optional
        Sampling frequency in Hz of the spike times, only needed if spike times and LFP have different sampling rates
    lfp_fs : float
        Sampling frequency in Hz of the LFP data
    filter_method : str, optional
        Method to use for filtering, either 'wavelet' or 'butter' (default: 'butter')
    freq_of_interest : float, optional
        Desired frequency for wavelet phase extraction, required if filter_method='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: 2.0)
    filtered_lfp_phase : np.ndarray, optional
        Pre-computed instantaneous phase of the filtered LFP. If provided, the function will skip the filtering step.

    Returns
    -------
    float
        Phase Locking Value (unbiased)
    """

    if spike_fs is None:
        spike_fs = lfp_fs
    # Convert spike times to sample indices
    spike_times_seconds = spike_times / spike_fs

    # Then convert from seconds to samples at the new sampling rate
    spike_indices = np.round(spike_times_seconds * lfp_fs).astype(int)

    # Filter indices to ensure they're within bounds of the LFP signal
    if filtered_lfp_phase is not None:
        valid_indices = [idx for idx in spike_indices if 0 <= idx < len(filtered_lfp_phase)]
    else:
        valid_indices = [idx for idx in spike_indices if 0 <= idx < len(lfp_data)]

    if len(valid_indices) <= 1:
        return 0

    # Get instantaneous phase
    if filtered_lfp_phase is None:
        instantaneous_phase = get_lfp_phase(lfp_data=lfp_data, filter_method=filter_method, 
                                           freq_of_interest=freq_of_interest, lowcut=lowcut, 
                                           highcut=highcut, bandwidth=bandwidth, fs=lfp_fs)
    else:
        instantaneous_phase = filtered_lfp_phase

    # Get phases at spike times
    spike_phases = instantaneous_phase[valid_indices]

    # Number of spikes
    N = len(spike_phases)

    # Convert phases to unit vectors in the complex plane
    unit_vectors = np.exp(1j * spike_phases)

    # Sum of all unit vectors (resultant vector)
    resultant_vector = np.sum(unit_vectors)

    # Calculate plv^2 * N
    plv2n = (resultant_vector * resultant_vector.conjugate()).real / N  # plv^2 * N
    plv = (plv2n / N) ** 0.5
    ppc = (plv2n - 1) / (N - 1)  # ppc = (plv^2 * N - 1) / (N - 1)
    plv_unbiased = np.fmax(ppc, 0.) ** 0.5  # ensure non-negative

    return plv_unbiased

bmtool.analysis.entrainment.calculate_ppc(spike_times=None, lfp_data=None, spike_fs=None, lfp_fs=None, filter_method='wavelet', freq_of_interest=None, lowcut=None, highcut=None, bandwidth=2.0, ppc_method='numpy', filtered_lfp_phase=None)

Calculate Pairwise Phase Consistency (PPC) between spike times and LFP signal. Based on https://www.sciencedirect.com/science/article/pii/S1053811910000959

Parameters:

Name Type Description Default
spike_times ndarray

Array of spike times

 None
lfp_data ndarray

Local field potential time series data. Not required if filtered_lfp_phase is provided.

 None
spike_fs float

Sampling frequency in Hz of the spike times, only needed if spike times and LFP have different sampling rates

 None
lfp_fs float

Sampling frequency in Hz of the LFP data

 None
filter_method str

Method to use for filtering, either 'wavelet' or 'butter' (default: 'wavelet')

 'wavelet'
freq_of_interest float

Desired frequency for wavelet phase extraction, required if filter_method='wavelet'

 None
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: 2.0)

 2.0
ppc_method str

Algorithm to use for PPC calculation: 'numpy', 'numba', or 'gpu' (default: 'numpy')

 'numpy'
filtered_lfp_phase ndarray

Pre-computed instantaneous phase of the filtered LFP. If provided, the function will skip the filtering step.

 None

Returns:

Type Description
float

Pairwise Phase Consistency value
Source code in bmtool/analysis/entrainment.py 
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
244
245
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
276
277
278
279
280
281
282
283
284
285
286
287
def calculate_ppc(spike_times: np.ndarray = None, lfp_data: np.ndarray = None, spike_fs: float = None,
                  lfp_fs: float = None, filter_method: str = 'wavelet', freq_of_interest: float = None,
                  lowcut: float = None, highcut: float = None, bandwidth: float = 2.0, 
                  ppc_method: str = 'numpy', filtered_lfp_phase: np.ndarray = None) -> float:
    """
    Calculate Pairwise Phase Consistency (PPC) between spike times and LFP signal.
    Based on https://www.sciencedirect.com/science/article/pii/S1053811910000959

    Parameters
    ----------
    spike_times : np.ndarray
        Array of spike times
    lfp_data : np.ndarray
        Local field potential time series data. Not required if filtered_lfp_phase is provided.
    spike_fs : float, optional
        Sampling frequency in Hz of the spike times, only needed if spike times and LFP have different sampling rates
    lfp_fs : float
        Sampling frequency in Hz of the LFP data
    filter_method : str, optional
        Method to use for filtering, either 'wavelet' or 'butter' (default: 'wavelet')
    freq_of_interest : float, optional
        Desired frequency for wavelet phase extraction, required if filter_method='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: 2.0)
    ppc_method : str, optional
        Algorithm to use for PPC calculation: 'numpy', 'numba', or 'gpu' (default: 'numpy')
    filtered_lfp_phase : np.ndarray, optional
        Pre-computed instantaneous phase of the filtered LFP. If provided, the function will skip the filtering step.

    Returns
    -------
    float
        Pairwise Phase Consistency value
    """
    if spike_fs is None:
        spike_fs = lfp_fs
    # Convert spike times to sample indices
    spike_times_seconds = spike_times / spike_fs

    # Then convert from seconds to samples at the new sampling rate
    spike_indices = np.round(spike_times_seconds * lfp_fs).astype(int)

    # Filter indices to ensure they're within bounds of the LFP signal
    if filtered_lfp_phase is not None:
        valid_indices = [idx for idx in spike_indices if 0 <= idx < len(filtered_lfp_phase)]
    else:
        valid_indices = [idx for idx in spike_indices if 0 <= idx < len(lfp_data)]

    if len(valid_indices) <= 1:
        return 0

    # Get instantaneous phase
    if filtered_lfp_phase is None:
        instantaneous_phase = get_lfp_phase(lfp_data=lfp_data, filter_method=filter_method, 
                                           freq_of_interest=freq_of_interest, lowcut=lowcut, 
                                           highcut=highcut, bandwidth=bandwidth, fs=lfp_fs)
    else:
        instantaneous_phase = filtered_lfp_phase

    # Get phases at spike times
    spike_phases = instantaneous_phase[valid_indices]

    n_spikes = len(spike_phases)

    # Calculate PPC (Pairwise Phase Consistency)
    if n_spikes <= 1:
        return 0

    # Explicit calculation of pairwise phase consistency
    # Vectorized computation for efficiency
    if ppc_method == 'numpy':
        i, j = np.triu_indices(n_spikes, k=1)
        phase_diff = spike_phases[i] - spike_phases[j]
        sum_cos_diff = np.sum(np.cos(phase_diff))
        ppc = ((2 / (n_spikes * (n_spikes - 1))) * sum_cos_diff)
    elif ppc_method == 'numba':
        ppc = _ppc_parallel_numba(spike_phases)
    elif ppc_method == 'gpu':
        ppc = _ppc_gpu(spike_phases)
    else:
        raise ValueError("Please use a supported ppc method currently that is numpy, numba or gpu")
    return ppc

bmtool.analysis.entrainment.calculate_ppc2(spike_times=None, lfp_data=None, spike_fs=None, lfp_fs=None, filter_method='wavelet', freq_of_interest=None, lowcut=None, highcut=None, bandwidth=2.0, filtered_lfp_phase=None)

-----------------------------------------------------------------------------

PPC2 Calculation (Vinck et al., 2010)

-----------------------------------------------------------------------------

Equation(Original):

PPC = (2 / (n * (n - 1))) * sum(cos(φ_i - φ_j) for all i < j)

Optimized Formula (Algebraically Equivalent):

PPC = (|sum(e^(i*φ_j))|^2 - n) / (n * (n - 1))

-----------------------------------------------------------------------------

Parameters:

Name Type Description Default
spike_times ndarray

Array of spike times

 None
lfp_data ndarray

Local field potential time series data. Not required if filtered_lfp_phase is provided.

 None
spike_fs float

Sampling frequency in Hz of the spike times, only needed if spike times and LFP have different sampling rates

 None
lfp_fs float

Sampling frequency in Hz of the LFP data

 None
filter_method str

Method to use for filtering, either 'wavelet' or 'butter' (default: 'wavelet')

 'wavelet'
freq_of_interest float

Desired frequency for wavelet phase extraction, required if filter_method='wavelet'

 None
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: 2.0)

 2.0
filtered_lfp_phase ndarray

Pre-computed instantaneous phase of the filtered LFP. If provided, the function will skip the filtering step.

 None

Returns:

Type Description
float

Pairwise Phase Consistency 2 (PPC2) value
Source code in bmtool/analysis/entrainment.py 
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
326
327
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
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
def calculate_ppc2(spike_times: np.ndarray = None, lfp_data: np.ndarray = None, spike_fs: float = None,
                  lfp_fs: float = None, filter_method: str = 'wavelet', freq_of_interest: float = None,
                  lowcut: float = None, highcut: float = None, bandwidth: float = 2.0,
                  filtered_lfp_phase: np.ndarray = None) -> float:
    """
    # -----------------------------------------------------------------------------
    # PPC2 Calculation (Vinck et al., 2010) 
    # -----------------------------------------------------------------------------
    # Equation(Original):
    #   PPC = (2 / (n * (n - 1))) * sum(cos(φ_i - φ_j) for all i < j)
    # Optimized Formula (Algebraically Equivalent):
    #   PPC = (|sum(e^(i*φ_j))|^2 - n) / (n * (n - 1))
    # -----------------------------------------------------------------------------

    Parameters
    ----------
    spike_times : np.ndarray
        Array of spike times
    lfp_data : np.ndarray
        Local field potential time series data. Not required if filtered_lfp_phase is provided.
    spike_fs : float, optional
        Sampling frequency in Hz of the spike times, only needed if spike times and LFP have different sampling rates
    lfp_fs : float
        Sampling frequency in Hz of the LFP data
    filter_method : str, optional
        Method to use for filtering, either 'wavelet' or 'butter' (default: 'wavelet')
    freq_of_interest : float, optional
        Desired frequency for wavelet phase extraction, required if filter_method='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: 2.0)
    filtered_lfp_phase : np.ndarray, optional
        Pre-computed instantaneous phase of the filtered LFP. If provided, the function will skip the filtering step.

    Returns
    -------
    float
        Pairwise Phase Consistency 2 (PPC2) value
    """

    if spike_fs is None:
        spike_fs = lfp_fs
    # Convert spike times to sample indices
    spike_times_seconds = spike_times / spike_fs

    # Then convert from seconds to samples at the new sampling rate
    spike_indices = np.round(spike_times_seconds * lfp_fs).astype(int)

    # Filter indices to ensure they're within bounds of the LFP signal
    if filtered_lfp_phase is not None:
        valid_indices = [idx for idx in spike_indices if 0 <= idx < len(filtered_lfp_phase)]
    else:
        valid_indices = [idx for idx in spike_indices if 0 <= idx < len(lfp_data)]

    if len(valid_indices) <= 1:
        return 0

    # Get instantaneous phase
    if filtered_lfp_phase is None:
        instantaneous_phase = get_lfp_phase(lfp_data=lfp_data, filter_method=filter_method, 
                                           freq_of_interest=freq_of_interest, lowcut=lowcut, 
                                           highcut=highcut, bandwidth=bandwidth, fs=lfp_fs)
    else:
        instantaneous_phase = filtered_lfp_phase

    # Get phases at spike times
    spike_phases = instantaneous_phase[valid_indices]

    # Calculate PPC2 according to Vinck et al. (2010), Equation 6
    n = len(spike_phases)

    if n <= 1:
        return 0

    # Convert phases to unit vectors in the complex plane
    unit_vectors = np.exp(1j * spike_phases)

    # Calculate the resultant vector
    resultant_vector = np.sum(unit_vectors)

    # PPC2 = (|∑(e^(i*φ_j))|² - n) / (n * (n - 1))
    ppc2 = (np.abs(resultant_vector)**2 - n) / (n * (n - 1))

    return ppc2