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.

Args: 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'. Defaults to True. config (str, optional): Will label the cell type of each spike. groupby (str or list of str, optional): The column(s) to group by. Defaults to 'pop_name'.

Returns: pd.DataFrame: A pandas DataFrame containing 'node_ids' and 'timestamps' columns from the spike data.

Example: df = load_spikes_to_df("spikes.h5", "cortex")

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
def load_spikes_to_df(spike_file: str, network_name: str, sort: bool = True, config: str = None, groupby: str = 'pop_name') -> pd.DataFrame:
    """
    Load spike data from an HDF5 file into a pandas DataFrame.

    Args:
        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'. Defaults to True.
        config (str, optional): Will label the cell type of each spike.
        groupby (str or list of str, optional): The column(s) to group by. Defaults to 'pop_name'.

    Returns:
        pd.DataFrame: A pandas DataFrame containing 'node_ids' and 'timestamps' columns from the spike data.

    Example:
        df = load_spikes_to_df("spikes.h5", "cortex")
    """
    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 
 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
 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
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, frequeny=False)

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

Args: 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. frequeny (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 frequeny flag.

Raises: ValueError: If both time and time_points are None.

Source code in bmtool/analysis/spikes.py 
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
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, frequeny: bool = False) -> np.ndarray:
    """
    Calculate the spike count or frequency histogram over specified time intervals.

    Args:
        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.
        frequeny (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 `frequeny` 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 frequeny:
        spike_rate = 1000 / dt * spike_rate

    return spike_rate

bmtool.analysis.spikes.get_population_spike_rate(spikes, 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.

Args: spikes (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 is 400. t_start (float, optional): Start time (in milliseconds) for spike rate calculation. Default is 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 is 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 is None. save (bool, optional): Whether to save the calculated population spike rate to a file. Default is False. save_path (Optional[str], optional): Directory path where the file should be saved if save is True. If save is True and save_path is None, a ValueError is raised. normalize (bool, optional): Whether to normalize the spike rates for each population to a range of [0, 1]. Default is 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.

Source code in bmtool/analysis/spikes.py 
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
244
245
246
247
248
249
250
251
252
253
254
def get_population_spike_rate(spikes: 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.

    Args:
        spikes (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 is 400.
        t_start (float, optional): Start time (in milliseconds) for spike rate calculation. Default is 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 is 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 is None.
        save (bool, optional): Whether to save the calculated population spike rate to a file. Default is False.
        save_path (Optional[str], optional): Directory path where the file should be saved if `save` is True. If `save` is True and `save_path` is None, a ValueError is raised.
        normalize (bool, optional): Whether to normalize the spike rates for each population to a range of [0, 1]. Default is 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 spikes['pop_name'].unique():
        ps = spikes[spikes['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 = spikes['timestamps'].max()

        filtered_spikes = spikes[
            (spikes['pop_name'] == pop_name) & 
            (spikes['timestamps'] > t_start) & 
            (spikes['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 
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
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 
 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
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 
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
220
221
222
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 
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
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 
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
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)

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

Source code in bmtool/analysis/lfp.py 
279
280
281
282
283
284
285
286
def wavelet_filter(x: np.ndarray, freq: float, fs: float, bandwidth: float = 1.0, axis: int = -1) -> np.ndarray:
    """
    Compute the Continuous Wavelet Transform (CWT) for a specified frequency using a complex Morlet wavelet.
    """
    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 
289
290
291
292
293
294
295
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.calculate_signal_signal_plv(x1, x2, fs, freq_of_interest=None, method='wavelet', lowcut=None, highcut=None, bandwidth=2.0)

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

Parameters: - x1, x2: Input signals (1D arrays, same length) - fs: Sampling frequency - freq_of_interest: Desired frequency for wavelet PLV calculation - method: 'wavelet' or 'hilbert' to choose the PLV calculation method - lowcut, highcut: Cutoff frequencies for the Hilbert method - bandwidth: Bandwidth parameter for the wavelet

Returns: - plv: Phase Locking Value (1D array)

Source code in bmtool/analysis/lfp.py 
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
def calculate_signal_signal_plv(x1: np.ndarray, x2: np.ndarray, fs: float, freq_of_interest: float = None, 
                  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:
    - x1, x2: Input signals (1D arrays, same length)
    - fs: Sampling frequency
    - freq_of_interest: Desired frequency for wavelet PLV calculation
    - method: 'wavelet' or 'hilbert' to choose the PLV calculation method
    - lowcut, highcut: Cutoff frequencies for the Hilbert method
    - bandwidth: Bandwidth parameter for the wavelet

    Returns:
    - plv: Phase Locking Value (1D array)
    """
    if len(x1) != len(x2):
        raise ValueError("Input signals must have the same length.")

    if 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=x1, freq=freq_of_interest, fs=fs, bandwidth=bandwidth)
        theta2 = wavelet_filter(x=x2, freq=freq_of_interest, fs=fs, bandwidth=bandwidth)

    elif method == 'hilbert':
        if lowcut is None or highcut is None:
            print("Lowcut and or highcut were not definded, signal will not be filter and just take hilbert transform for plv calc")

        if lowcut and highcut:
            # Bandpass filter and get the analytic signal using the Hilbert transform
            x1 = butter_bandpass_filter(x1, lowcut, highcut, fs)
            x2 = butter_bandpass_filter(x2, lowcut, highcut, fs)

        # Get phase using the Hilbert transform
        theta1 = signal.hilbert(x1)
        theta2 = signal.hilbert(x2)

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

    # 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.lfp.calculate_spike_lfp_plv(spike_times=None, lfp_signal=None, spike_fs=None, lfp_fs=None, method='hilbert', freq_of_interest=None, lowcut=None, highcut=None, bandwidth=2.0)

Calculate spike-lfp phase locking value Based on https://www.sciencedirect.com/science/article/pii/S1053811910000959

Parameters: - spike_times: Array of spike times - lfp_signal: Local field potential time series - spike_fs: Sampling frequency in Hz of the spike times only needed if spikes times and lfp has different fs - lfp_fs : Sampling frequency in Hz of the LFP - method: 'wavelet' or 'hilbert' to choose the phase extraction method - freq_of_interest: Desired frequency for wavelet phase extraction - lowcut, highcut: Cutoff frequencies for bandpass filtering (Hilbert method) - bandwidth: Bandwidth parameter for the wavelet

Returns: - ppc1: Phase-Phase Coupling value - spike_phases: Phases at spike times

Source code in bmtool/analysis/lfp.py 
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
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
def calculate_spike_lfp_plv(spike_times: np.ndarray = None, lfp_signal: np.ndarray = None, spike_fs : float = None,
                   lfp_fs: float = None, method: str = 'hilbert', freq_of_interest: float = None,
                   lowcut: float = None, highcut: float = None,
                   bandwidth: float = 2.0) -> tuple:
    """
    Calculate spike-lfp phase locking value Based on https://www.sciencedirect.com/science/article/pii/S1053811910000959

    Parameters:
    - spike_times: Array of spike times 
    - lfp_signal: Local field potential time series
    - spike_fs: Sampling frequency in Hz of the spike times only needed if spikes times and lfp has different fs
    - lfp_fs : Sampling frequency in Hz of the LFP
    - method: 'wavelet' or 'hilbert' to choose the phase extraction method
    - freq_of_interest: Desired frequency for wavelet phase extraction
    - lowcut, highcut: Cutoff frequencies for bandpass filtering (Hilbert method)
    - bandwidth: Bandwidth parameter for the wavelet

    Returns:
    - ppc1: Phase-Phase Coupling value
    - spike_phases: Phases at spike times
    """

    if spike_fs == 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
    valid_indices = [idx for idx in spike_indices if 0 <= idx < len(lfp_signal)]
    if len(valid_indices) <= 1:
        return 0, np.array([])

    # Extract phase using the specified method
    if method == 'wavelet':
        if freq_of_interest is None:
            raise ValueError("freq_of_interest must be provided for the wavelet method.")

        # Apply CWT to extract phase at the frequency of interest
        lfp_complex = wavelet_filter(x=lfp_signal, freq=freq_of_interest, fs=lfp_fs, bandwidth=bandwidth)
        instantaneous_phase = np.angle(lfp_complex)

    elif method == 'hilbert':
        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 PPC1 calculation")
            filtered_lfp = lfp_signal
        else:
            # Bandpass filter the signal
            filtered_lfp = butter_bandpass_filter(lfp_signal, lowcut, highcut, lfp_fs)

        # Get phase using the Hilbert transform
        analytic_signal = signal.hilbert(filtered_lfp)
        instantaneous_phase = np.angle(analytic_signal)

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

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

    # Calculate PPC1
    n = len(spike_phases)

    # 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)

    # Plv is the squared length of the resultant vector divided by n²
    plv = (np.abs(resultant_vector) ** 2) / (n ** 2)

    return plv

bmtool.analysis.lfp.calculate_ppc(spike_times=None, lfp_signal=None, spike_fs=None, lfp_fs=None, method='hilbert', freq_of_interest=None, lowcut=None, highcut=None, bandwidth=2.0, ppc_method='numpy')

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

Parameters: - spike_times: Array of spike times - lfp_signal: Local field potential time series - spike_fs: Sampling frequency in Hz of the spike times only needed if spikes times and lfp has different fs - lfp_fs: Sampling frequency in Hz of the LFP - method: 'wavelet' or 'hilbert' to choose the phase extraction method - freq_of_interest: Desired frequency for wavelet phase extraction - lowcut, highcut: Cutoff frequencies for bandpass filtering (Hilbert method) - bandwidth: Bandwidth parameter for the wavelet - ppc_method: which algo to use for PPC calculate can be numpy, numba or gpu

Returns: - ppc: Pairwise Phase Consistency value

Source code in bmtool/analysis/lfp.py 
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
550
551
552
553
554
555
556
557
558
559
560
561
562
def calculate_ppc(spike_times: np.ndarray = None, lfp_signal: np.ndarray = None, spike_fs: float = None,
                  lfp_fs: float = None, method: str = 'hilbert', freq_of_interest: float = None,
                  lowcut: float = None, highcut: float = None,
                  bandwidth: float = 2.0,ppc_method: str = 'numpy') -> tuple:
    """
    Calculate Pairwise Phase Consistency (PPC) between spike times and LFP signal.
    Based on https://www.sciencedirect.com/science/article/pii/S1053811910000959

    Parameters:
    - spike_times: Array of spike times 
    - lfp_signal: Local field potential time series
    - spike_fs: Sampling frequency in Hz of the spike times only needed if spikes times and lfp has different fs
    - lfp_fs: Sampling frequency in Hz of the LFP
    - method: 'wavelet' or 'hilbert' to choose the phase extraction method
    - freq_of_interest: Desired frequency for wavelet phase extraction
    - lowcut, highcut: Cutoff frequencies for bandpass filtering (Hilbert method)
    - bandwidth: Bandwidth parameter for the wavelet
    - ppc_method: which algo to use for PPC calculate can be numpy, numba or gpu

    Returns:
    - ppc: 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
    valid_indices = [idx for idx in spike_indices if 0 <= idx < len(lfp_signal)]
    if len(valid_indices) <= 1:
        return 0, np.array([])

    # Extract phase using the specified method
    if method == 'wavelet':
        if freq_of_interest is None:
            raise ValueError("freq_of_interest must be provided for the wavelet method.")

        # Apply CWT to extract phase at the frequency of interest
        lfp_complex = wavelet_filter(x=lfp_signal, freq=freq_of_interest, fs=lfp_fs, bandwidth=bandwidth)
        instantaneous_phase = np.angle(lfp_complex)

    elif method == 'hilbert':
        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 PPC calculation")
            filtered_lfp = lfp_signal
        else:
            # Bandpass filter the signal
            filtered_lfp = butter_bandpass_filter(lfp_signal, lowcut, highcut, lfp_fs)

        # Get phase using the Hilbert transform
        analytic_signal = signal.hilbert(filtered_lfp)
        instantaneous_phase = np.angle(analytic_signal)

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

    # 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, spike_phases

    # Explicit calculation of pairwise phase consistency
    sum_cos_diff = 0.0

    # # Σᵢ Σⱼ₍ᵢ₊₁₎ f(θᵢ, θⱼ)
    # for i in range(n_spikes - 1):  # For each spike i
    #     for j in range(i + 1, n_spikes):  # For each spike j > i
    #         # Calculate the phase difference between spikes i and j
    #         phase_diff = spike_phases[i] - spike_phases[j]

    #         #f(θᵢ, θⱼ) = cos(θᵢ)cos(θⱼ) + sin(θᵢ)sin(θⱼ) can become #f(θᵢ, θⱼ) = cos(θᵢ - θⱼ)
    #         cos_diff = np.cos(phase_diff)

    #         # Add to the sum
    #         sum_cos_diff += cos_diff

    # # Calculate PPC according to the equation
    # # PPC = (2 / (N(N-1))) * Σᵢ Σⱼ₍ᵢ₊₁₎ f(θᵢ, θⱼ)
    # ppc = ((2 / (n_spikes * (n_spikes - 1))) * sum_cos_diff)

    # same as above (i think) but with vectorized computation and memory fixes so it wont take forever to run.
    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 ExceptionType("Please use a supported ppc method currently that is numpy, numba or gpu")
    return ppc

bmtool.analysis.lfp.calculate_ppc2(spike_times=None, lfp_signal=None, spike_fs=None, lfp_fs=None, method='hilbert', freq_of_interest=None, lowcut=None, highcut=None, bandwidth=2.0)

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

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: Array of spike times - lfp_signal: Local field potential time series - spike_fs: Sampling frequency in Hz of the spike times only needed if spikes times and lfp has different fs - lfp_fs: Sampling frequency in Hz of the LFP - method: 'wavelet' or 'hilbert' to choose the phase extraction method - freq_of_interest: Desired frequency for wavelet phase extraction - lowcut, highcut: Cutoff frequencies for bandpass filtering (Hilbert method) - bandwidth: Bandwidth parameter for the wavelet

Returns: - ppc2: Pairwise Phase Consistency 2 value - spike_phases: Phases at spike times

Source code in bmtool/analysis/lfp.py 
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
647
648
649
def calculate_ppc2(spike_times: np.ndarray = None, lfp_signal: np.ndarray = None, spike_fs: float = None,
                  lfp_fs: float = None, method: str = 'hilbert', freq_of_interest: float = None,
                  lowcut: float = None, highcut: float = None,
                  bandwidth: float = 2.0) -> tuple:
    """
    # -----------------------------------------------------------------------------
    # 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: Array of spike times 
    - lfp_signal: Local field potential time series
    - spike_fs: Sampling frequency in Hz of the spike times only needed if spikes times and lfp has different fs
    - lfp_fs: Sampling frequency in Hz of the LFP
    - method: 'wavelet' or 'hilbert' to choose the phase extraction method
    - freq_of_interest: Desired frequency for wavelet phase extraction
    - lowcut, highcut: Cutoff frequencies for bandpass filtering (Hilbert method)
    - bandwidth: Bandwidth parameter for the wavelet

    Returns:
    - ppc2: Pairwise Phase Consistency 2 value
    - spike_phases: Phases at spike times
    """

    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
    valid_indices = [idx for idx in spike_indices if 0 <= idx < len(lfp_signal)]
    if len(valid_indices) <= 1:
        return 0, np.array([])

    # Extract phase using the specified method
    if method == 'wavelet':
        if freq_of_interest is None:
            raise ValueError("freq_of_interest must be provided for the wavelet method.")

        # Apply CWT to extract phase at the frequency of interest
        lfp_complex = wavelet_filter(x=lfp_signal, freq=freq_of_interest, fs=lfp_fs, bandwidth=bandwidth)
        instantaneous_phase = np.angle(lfp_complex)

    elif method == 'hilbert':
        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 PPC2 calculation")
            filtered_lfp = lfp_signal
        else:
            # Bandpass filter the signal
            filtered_lfp = butter_bandpass_filter(lfp_signal, lowcut, highcut, lfp_fs)

        # Get phase using the Hilbert transform
        analytic_signal = signal.hilbert(filtered_lfp)
        instantaneous_phase = np.angle(analytic_signal)

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

    # 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, spike_phases

    # 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

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 
732
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
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 
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
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.plot_spectrogram(sxx_xarray, remove_aperiodic=None, log_power=False, plt_range=None, clr_freq_range=None, pad=0.03, ax=None)

Plot spectrogram. Determine color limits using value in frequency band clr_freq_range

Source code in bmtool/analysis/lfp.py 
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
def plot_spectrogram(sxx_xarray, remove_aperiodic=None, log_power=False,
                     plt_range=None, clr_freq_range=None, pad=0.03, ax=None):
    """Plot spectrogram. Determine color limits using value in frequency band clr_freq_range"""
    sxx = sxx_xarray.PSD.values.copy()
    t = sxx_xarray.time.values.copy()
    f = sxx_xarray.frequency.values.copy()

    cbar_label = 'PSD' if remove_aperiodic is None else 'PSD Residual'
    if log_power:
        with np.errstate(divide='ignore'):
            sxx = np.log10(sxx)
        cbar_label += ' dB' if log_power == 'dB' else ' log(power)'

    if remove_aperiodic is not None:
        f1_idx = 0 if f[0] else 1
        ap_fit = gen_aperiodic(f[f1_idx:], remove_aperiodic.aperiodic_params)
        sxx[f1_idx:, :] -= (ap_fit if log_power else 10 ** ap_fit)[:, None]
        sxx[:f1_idx, :] = 0.

    if log_power == 'dB':
        sxx *= 10

    if ax is None:
        _, ax = plt.subplots(1, 1)
    plt_range = np.array(f[-1]) if plt_range is None else np.array(plt_range)
    if plt_range.size == 1:
        plt_range = [f[0 if f[0] else 1] if log_power else 0., plt_range.item()]
    f_idx = (f >= plt_range[0]) & (f <= plt_range[1])
    if clr_freq_range is None:
        vmin, vmax = None, None
    else:
        c_idx = (f >= clr_freq_range[0]) & (f <= clr_freq_range[1])
        vmin, vmax = sxx[c_idx, :].min(), sxx[c_idx, :].max()

    f = f[f_idx]
    pcm = ax.pcolormesh(t, f, sxx[f_idx, :], shading='gouraud', vmin=vmin, vmax=vmax)
    if 'cone_of_influence_frequency' in sxx_xarray:
        coif = sxx_xarray.cone_of_influence_frequency
        ax.plot(t, coif)
        ax.fill_between(t, coif, step='mid', alpha=0.2)
    ax.set_xlim(t[0], t[-1])
    #ax.set_xlim(t[0],0.2)
    ax.set_ylim(f[0], f[-1])
    plt.colorbar(mappable=pcm, ax=ax, label=cbar_label, pad=pad)
    ax.set_xlabel('Time (sec)')
    ax.set_ylabel('Frequency (Hz)')
    return sxx