Harmonic connectivity

Harmonic connectivity#

Cross-channel coupling, cross-frequency, and connectivity matrices. Everything that compares two-or-more signals lives here.

Quick start#

Peak-based connectivity (between extracted peak lists, one scalar per electrode pair) — uses the legacy peak-extraction pipeline:

from biotuner.harmonic_connectivity import harmonic_connectivity

hc = harmonic_connectivity(sf=1000, data=data_array, peaks_function="EMD",
                           precision=0.5, min_freq=2, max_freq=30, n_peaks=5)
H_mat  = hc.compute_harm_connectivity(metric="harmsim")           # legacy H
PC_mat = hc.compute_peak_phase_coupling_connectivity(coupling_metric="nm_plv")
R_mat  = hc.compute_peak_resonance_connectivity(combine="product")

Spectrum-based cross-channel resonance (per-frequency H/PC/R between two signals, with the same swappable kernels/metrics/combine rules as biotuner.resonance):

from biotuner.harmonic_connectivity import compute_cross_resonance

cross = compute_cross_resonance(sig1, sig2, sf=1000)
# cross.resonance_spectrum["1to2"]   — asymmetric (sig1 at i, sig2 at j)
# cross.resonance_spectrum["2to1"]   — transposed
# cross.resonance_spectrum["all"]    — symmetrized average
# cross.factors["H"][...], cross.factors["PC"][...]
# cross.summaries["H"/"PC"/"R"]      — complexity dict per spectrum

Connectivity matrices (loop compute_cross_resonance over all electrode pairs):

M = hc.compute_cross_resonance_connectivity(
    factor="R", flavor="all", aggregate="peak_to_median",
)

Statistical inference via surrogate-z-scoring (the principled way to separate true cross-channel phase coupling from broadband-power artifacts):

obs, z, p = hc.compute_cross_resonance_connectivity_zscore(
    surrogate_kind="iaaft", n_surrogates=200,
)

Sister modules#

  • biotuner.harmonic_spectrum — single-signal H spectrum.

  • biotuner.resonance — single-signal H × PC = R pipeline, including the registries this module dispatches against (PAIRWISE_COUPLING_METRICS, HARMONIC_KERNELS, RATIO_KERNELS, COMBINE_RULES).

Public API#

The class:

  • harmonic_connectivity — peak-based + new spectrum-based methods:

    • .compute_harm_connectivity(metric=...) — legacy peak-based H

    • .compute_peak_phase_coupling_connectivity(coupling_metric=...) — registry-based PC

    • .compute_peak_resonance_connectivity(combine=...) — H × PC peak scalar

    • .compute_cross_resonance_connectivity(factor=, flavor=, aggregate=) — spectrum-based matrix

    • .compute_cross_resonance_connectivity_zscore(surrogate_kind=, n_surrogates=) — surrogate-normalized inference

    • .compute_harmonic_spectrum_connectivity(...) — legacy cross-spectrum DataFrame

Spectrum-based cross-channel orchestrator:

  • compute_cross_resonance() — single-call cross-channel R(f) on a signal pair

  • CrossResonanceResult — output dataclass with 3 reducer flavors ('1to2' / '2to1' / 'all') per factor. The full N×N matrices are accessible as result.harmonicity_matrix and result.phase_coupling_matrix when the config has return_intermediates=True.

Legacy shim:

  • compute_cross_spectrum_harmonicity() — delegates to compute_cross_resonance with explicit legacy preset; bit-exact reproduction guaranteed for paper workflows.

Standalone coupling utilities (kept for backward compat):

  • wPLI_crossfreq(), wPLI_multiband()

  • cross_frequency_rrci(), rhythmic_ratio_coupling_imaginary(), compute_rhythmic_ratio()

  • n_m_phase_locking()

  • compute_mutual_information(), MI_spectral()

IMF-based utilities (different abstraction layer):

  • HilbertHuang1D_nopeaks(), EMD_time_resolved_harmonicity()

  • temporal_correlation_fdr()

biotuner.harmonic_connectivity — cross-channel coupling, cross-frequency, and connectivity matrices.

Module type: Object

This module is the cross-channel counterpart to biotuner.harmonic_spectrum (single-signal H) and biotuner.resonance (single-signal R = H · PC). Everything that compares TWO-OR-MORE signals lives here.

Quick start#

Peak-based connectivity (between extracted peak lists, one scalar per electrode pair) — uses the legacy peak-extraction pipeline:

from biotuner.harmonic_connectivity import harmonic_connectivity
hc = harmonic_connectivity(sf=1000, data=data_array, peaks_function="EMD",
                            precision=0.5, min_freq=2, max_freq=30, n_peaks=5)
H_mat = hc.compute_harm_connectivity(metric="harmsim")       # legacy H
PC_mat = hc.compute_peak_phase_coupling_connectivity(coupling_metric="nm_plv")
R_mat = hc.compute_peak_resonance_connectivity(combine="product")

Spectrum-based cross-channel resonance (per-frequency H/PC/R between two signals, with the same swappable kernels/metrics/combine rules as biotuner.resonance):

from biotuner.harmonic_connectivity import compute_cross_resonance
from biotuner.resonance import ResonanceConfig

cross = compute_cross_resonance(sig1, sig2, sf=1000)
# cross.resonance_spectrum["1to2"]  — asymmetric (sig1 at i, sig2 at j)
# cross.resonance_spectrum["2to1"]  — transposed
# cross.resonance_spectrum["all"]   — symmetrized average
# cross.factors["H"][...], cross.factors["PC"][...]
# cross.summaries["H"/"PC"/"R"]     — complexity dict per spectrum

Connectivity matrices (loop compute_cross_resonance over all electrode pairs):

M = hc.compute_cross_resonance_connectivity(
    factor="R", flavor="all", aggregate="peak_to_median",
)

Statistical inference via surrogate-z-scoring (the principled way to separate true cross-channel phase coupling from broadband-power artifacts):

obs, z, p = hc.compute_cross_resonance_connectivity_zscore(
    surrogate_kind="iaaft", n_surrogates=200,
)

Sister modules#

  • biotuner.harmonic_spectrum — single-signal H spectrum.

  • biotuner.resonance — single-signal H × PC = R pipeline, including the registries this module dispatches against (PAIRWISE_COUPLING_METRICS, HARMONIC_KERNELS, RATIO_KERNELS, COMBINE_RULES).

class harmonic_connectivity(sf=None, data=None, peaks_function='EMD', precision=0.1, n_harm=10, harm_function='mult', min_freq=2, max_freq=80, n_peaks=5)[source]#

Bases: object

Class used to compute harmonicity metrics between pairs of sensors.

Parameters:

  • sf (int) – sampling frequency (in Hz)

  • data (2D array (elec, numDataPoints)) – Electrodes x Time series to analyse.

  • peaks_function (str, default=’EMD’) – See compute_biotuner class for details.

  • precision (float, default=0.1) – Precision of the peaks (in Hz) When HH1D_max is used, bins are in log scale.

  • n_harm (int, default=10) – Set the number of harmonics to compute in harmonic_fit function

  • harm_function (str, default=’mult’) – Computes harmonics from iterative multiplication (x, 2x, 3x, …nx) or division (x, x/2, x/3, …x/n). Choose between ‘mult’ and ‘div’

  • min_freq (float, default=2) – Minimum frequency (in Hz) to consider for peak extraction.

  • max_freq (float, default=80) – Maximum frequency (in Hz) to consider for peak extraction.

  • n_peaks (int, default=5) – Number of peaks to extract per frequency band.

compute_harm_connectivity(metric='harmsim', delta_lim=20, save=False, savename='_', graph=True, FREQ_BANDS=None, max_denom_rrci=16)[source]#

Computes the harmonic connectivity matrix between electrodes.

Parameters:

  • metric (str, optional) – The metric to use for computing harmonic connectivity. Default is ‘harmsim’.

    Possible values are:

    • ‘harmsim’:

      computes the harmonic similarity between each pair of peaks from the two electrodes. It calculates the ratio between each pair of peaks and computes the mean harmonic similarity.

    • ‘euler’:

      computes the Euler’s Gradus Suavitatis on the concatenated peaks of the two electrodes.

    • ‘harm_fit’:

      computes the number of common harmonics between each pair of peaks from the two electrodes. It evaluates the harmonic fit between each peak pair and counts the number of common harmonics.

    • ‘subharm_tension’:

      computes the tension between subharmonics of two electrodes. It evaluates the tension between subharmonics of the two electrodes by comparing the subharmonics and their ratios.

    • ‘RRCi’:

      computes the Rhythmic Ratio Coupling with Imaginary Component (RRCi) metric between each pair of peaks from the two electrodes, using a bandwidth of 2 Hz and a max_denom of 16. This metric calculates the imaginary part of the complex phase differences between two filtered signals.

    • ‘wPLI_crossfreq’:

      computes the weighted Phase Lag Index (wPLI) for cross-frequency coupling between each pair of peaks from the two electrodes. The wPLI measures the phase synchronization between two signals, with a value close to 0 indicating no synchronization and a value close to 1 indicating perfect synchronization.

    • ‘wPLI_multiband’:

      computes the weighted Phase Lag Index (wPLI) for multiple frequency bands between the two electrodes. It calculates wPLI for each frequency band and returns an array of wPLI values for the defined frequency bands.

  • delta_lim (int, optional) – The delta limit for the subharmonic tension metric. Default is 20.

  • save (bool, optional) – Whether to save the connectivity matrix. Default is False.

  • savename (str, optional) – The name to use when saving the connectivity matrix. Default is ‘_’.

  • graph (bool, optional) – Whether to display a heatmap of the connectivity matrix. Default is True.

  • FREQ_BANDS (list, optional) – The frequency bands to use for the computation of the wPLI_multiband metric. Default is None. If None, the following frequency bands will be used: [2, 3.55], [3.55, 7.15], [7.15, 14.3], [14.3, 28.55], [28.55, 49.4].

  • max_denom_rrci (int, optional) – The maximum denominator to use for the computation of the RRCi metric. Default is 16.

Returns:

conn_matrix (numpy.ndarray) – The harmonic connectivity matrix between electrodes.

compute_peak_phase_coupling_connectivity(coupling_metric='nm_plv', coupling_metric_params=None, ratio_kernel='binary', ratio_kernel_params=None, bandwidth=3.0, aggregate='max', FREQ_BANDS=None, graph=True, save=False, savename='_')[source]#

Peak-based phase-coupling connectivity matrix.

For each electrode pair: extract peaks from both signals, for every peak pair (f_i, f_j) determine the n:m harmonic ratio via the chosen ratio kernel, bandpass-filter both signals at their peak frequencies, Hilbert-transform to obtain analytic signals, and compute the chosen pairwise coupling metric. Aggregate over peak pairs to a scalar.

Parameters:

  • coupling_metric (str, default=’nm_plv’) – Name of any registered pairwise coupling metric. Options include: 'nm_plv', 'nm_pli', 'nm_wpli', 'nm_rrci', 'nm_plv_canonical', 'nm_wpli_complex'. See biotuner.resonance.coupling for definitions.

  • coupling_metric_params (dict, optional) – Extra kwargs passed to the metric function.

  • ratio_kernel (str, default=’binary’) – Name of registered ratio kernel. 'binary' is the legacy n:m gate.

  • ratio_kernel_params (dict, optional) – Kwargs for the ratio kernel; defaults to {'max_nm': 3, 'tolerance': 0.05, 'fallback_to_1_1': True}.

  • bandwidth (float, default=3.0) – Hz width of the bandpass filter applied around each peak. Must be wide enough to capture the oscillation’s instantaneous-frequency excursions: a too-narrow window (e.g. 1 Hz) corrupts the Hilbert phase on non-stationary signals and collapses all coupling toward 0. ~3 Hz is a robust default; widen further for strongly drifting rhythms.

  • aggregate ({‘mean’, ‘max’, ‘sum’}, default=’max’) – How to reduce the per-peak-pair scalars to one number per electrode pair. 'max' reports the strongest n:m coupling across peak pairs and is recommended — 'mean' dilutes a genuine pairwise lock among the (often spurious) other peak-pair combinations.

  • FREQ_BANDS (list of [low, high], optional) – Override the default frequency-band partition used by peak extraction.

  • graph (bool, default=True) – If True, displays a heatmap of the result.

  • save, savename (passthrough to plotting.)

Returns:

matrix (ndarray (n_elec, n_elec)) – Pairwise scalar coupling values.

compute_peak_resonance_connectivity(harm_metric='harmsim', coupling_metric='nm_plv', combine='product', coupling_metric_params=None, ratio_kernel='binary', ratio_kernel_params=None, bandwidth=3.0, coupling_aggregate='max', delta_lim=20, FREQ_BANDS=None, graph=True, save=False, savename='_')[source]#

Peak-based resonance connectivity: H × PC per electrode pair.

Combines a peak-based harmonicity scalar (from compute_harm_connectivity()) and a peak-based phase-coupling scalar (from compute_peak_phase_coupling_connectivity()) via a registered combine rule. Default combine='product' gives the standard H · PC resonance, analogous to the per-frequency R(f) = H(f) · PC(f) in the single-channel framework.

Parameters:

  • harm_metric (str, default=’harmsim’) – Harmonic-similarity metric for the H scalar (see compute_harm_connectivity()).

  • coupling_metric (str, default=’nm_plv’) – Phase-coupling metric for the PC scalar (see compute_peak_phase_coupling_connectivity()).

  • combine (str, default=’product’) – Name in biotuner.resonance.registry.COMBINE_RULES. Options: 'product' (legacy), 'geomean', 'harmmean', 'min', 'weighted_log'.

  • Other parameters are forwarded to the two underlying methods.

Returns:

matrix (ndarray (n_elec, n_elec))

compute_cross_resonance_connectivity(config=None, factor='R', flavor='all', aggregate='peak_to_median', graph=True, save=False, savename='_', store_full_results=False)[source]#

Cross-resonance connectivity matrix: loop compute_cross_resonance over every electrode pair and reduce each result to a scalar.

Parameters:

  • config (ResonanceConfig, optional) – Forwarded to compute_cross_resonance(). Defaults to the refined cross-channel config (joint PC + n:m kernel).

  • factor ({‘H’, ‘PC’, ‘R’}, default ‘R’) – Which factor’s spectrum to reduce to a scalar.

  • flavor ({‘1to2’, ‘2to1’, ‘all’}, default ‘all’) – Which of the 3 reducer flavors to use.

  • aggregate (str, default 'peak_to_median') – How to reduce the per-bin spectrum to a single number per electrode pair. Simple aggregates inherit the joint-probability reducer’s broadband-power bias (any channel with broadband content ranks high); normalized aggregates measure peak structure ABOVE the spectrum’s noise floor:

    'max' — peak value. Caveat: rewards broadband-power channels

    (e.g. pink noise) because the joint-probability reducer gives large weighted sums everywhere the partner channel has power.

    'mean' — average. Dilutes peaks; same broadband bias. 'sum' — total area. Scales with bandwidth. 'peak' — value at the prominence-detected peak frequency.

    Mildly better but still inherits H’s broadband bias.

    'peak_to_median' — log10(max / median). Log-scale

    peak-to-floor ratio; robust to broadband channels. DEFAULT.

    'peak_over_median' — (max - median) / (max + median).

    DEPRECATED — saturates at 1.0 for any spectrum with a sharp peak relative to its baseline; prefer peak_to_median.

    'spectral_concentration' — fraction of energy in the top-3

    bins. Sharp peaks → high; diffuse spectra → low.

    'peak_z' — z-score of the dominant peak relative to the

    off-peak distribution (all detected peaks masked out, so multi-peak spectra are NOT penalized). Often inflates multi-peak spectra; use peak_to_median for the cleanest broadband-vs-focal discrimination.

  • graph (bool) – Heatmap if True.

  • store_full_results (bool, default False) – If True, also stores self.cross_resonance_results as a list of (i, j, CrossResonanceResult) for downstream graph analysis. Costly for large electrode counts (stores n_elec² ResonanceResults).

Returns:

matrix (ndarray (n_elec, n_elec)) – The scalar connectivity matrix.

Notes

There is no single scalar aggregate of R(f) that cleanly separates true phase-coupling from broadband-power overlap — H’s joint p1·p2 reducer rewards any partner power overlap, so a focal signal paired with broadband noise will produce sharp peaks in R(f) at the focal channel’s frequency. For paper-quality discrimination of true cross-channel phase coupling vs broadband artifact, use compute_cross_resonance_connectivity_zscore() with surrogate_kind='iaaft' — IAAFT surrogates preserve per-channel PSD while destroying cross-channel phase, so a high z-score signals genuine phase coupling above the broadband null.

The default aggregate='peak_to_median' was chosen empirically on the 6-channel validation dataset as the cleanest scalar discriminator: locked alpha pair > 1:2 harmonic pair > drifting alpha > pink-noise pair, in the expected order. peak_z has a tendency to inflate multi-peak spectra and is offered as an alternative.

compute_cross_resonance_connectivity_zscore(config=None, factor='R', flavor='all', aggregate='max', surrogate_kind='iaaft', n_surrogates=50, rng_seed=42, graph=True)[source]#

Surrogate-normalized z-scored cross-resonance connectivity matrix.

For each surrogate, the multichannel data is replaced by independently surrogated copies of each channel, the n_elec × n_elec connectivity matrix is computed, and the surrogate mean/std is collected per cell. The observed matrix is then z-scored against the surrogate distribution.

Parameters:

  • surrogate_kind ({‘phase_randomize’, ‘iaaft’, ‘time_shuffle’}) – Which surrogate generator to use (see biotuner.resonance.nulls). Default 'iaaft' — tightest cross-channel null per Fig 30.

  • n_surrogates (int, default 50) – Number of surrogate matrices.

  • Other args forwarded to :meth:`compute_cross_resonance_connectivity`.

Returns:

  • observed (ndarray (n_elec, n_elec))

  • z_matrix (ndarray (n_elec, n_elec))

  • p_matrix (ndarray (n_elec, n_elec) one-sided empirical p-value)

compute_IMF_correlation(nIMFs=5, freq_range=(1, 60), precision=0.5, delta_lim=50)[source]#

Compute the correlation, coherence, and peak frequency between each pair of IMFs for each pair of electrodes.

Returns:

df (pd.DataFrame) – DataFrame with columns: elec1, elec2, imf1, imf2, pearson, coherence, and peak_freq.

compute_time_resolved_harm_connectivity(sf, nIMFs, metric='harmsim', delta_lim=50)[source]#

Computes the time-resolved harmonic connectivity matrix between electrodes, which is a harmonic connectivity matrix for each intrinsic mode function (IMF), and each time point.

Parameters:

  • data (numpy.ndarray) – Input data with shape (num_electrodes, numDataPoints)

  • sf (int) – Sampling frequency

  • nIMFs (int) – Number of intrinsic mode functions (IMFs) to consider.

  • metric (str, optional) – The metric to use for computing harmonic connectivity. Default is ‘harmsim’.

  • delta_lim (int, optional) – The delta limit for the subharmonic tension metric. Default is 20.

Returns:

connectivity_matrices (numpy.ndarray) – Time-resolved harmonic connectivity matrices with shape (IMFs, numDataPoints, electrodes, electrodes).

Notes

!!! This method is very computationally expensive and can take a long time to run. !!!

transitional_connectivity(data=None, sf=None, mode='win_overlap', overlap=10, delta_lim=20, graph=False, n_trans_harm=3)[source]#

This function calculates the transitional connectivity among electrodes, utilizing concepts of transitional harmony and temporal correlation. It does so by employing a two-step approach:

  • Transitional Harmony Calculation: For every electrode, the transitional harmony is determined first.

    This measure reflects the patterns or sequences of signal transitions over time that each electrode experiences.

  • Temporal Correlation with FDR Correction: After obtaining transitional harmonies, the function evaluates the temporal correlation

    between these harmonies for each pair of electrodes. Temporal correlation quantifies how similar the timing and pattern of transitional harmonies are between each pair of electrodes.

The entire process takes into account multiple comparisons, utilizing False Discovery Rate (FDR) correction to reduce the likelihood of false positives. This correction is crucial in maintaining the validity and robustness of the results, especially when dealing with a large number of comparisons, as is the case with numerous electrodes.

Parameters:

  • data (numpy.ndarray) – Multichannel EEG data with shape (n_electrodes, n_timepoints).

  • sf (float) – Sampling frequency of the EEG data in Hz.

  • mode (str, optional, default=’win_overlap’) – The mode to compute the transitional harmony. ‘win_overlap’ computes the transitional harmony using a sliding window with overlap.

  • overlap (int, optional, default=10) – The percentage of overlap between consecutive windows when computing the transitional harmony. Default is 10.

  • delta_lim (int, optional, default=20) – The minimum delta value for the computation of subharmonic tension, in ms. Default is 20.

  • graph (bool, optional, default=False) – If True, it will plot the graph of the transitional harmony. Default is False.

  • n_trans_harm (int, optional, default=3) – Number of transitional harmonics to compute. Default is 3.

Returns:

  • conn_mat (numpy.ndarray) – Connectivity matrix of shape (n_electrodes, n_electrodes) representing the transitional connectivity between electrodes.

  • pval_mat (numpy.ndarray) – P-value matrix of shape (n_electrodes, n_electrodes) with FDR-corrected p-values for the computed connectivity values.

plot_conn_matrix(conn_matrix=None, node_names=None, n_lines=50)[source]#

Plots a connectivity matrix in a circle plot.

Parameters:

  • conn_matrix (ndarray, optional) – The connectivity matrix to plot. If None, uses the object’s own attribute conn_matrix. If conn_matrix is still None, a ValueError will be raised.

  • node_names (list, optional) – The labels for the nodes. If None, a range object will be converted to a list of string values for node names.

  • n_lines (int, default=50) – The number of lines to draw in the plot. Default is 50.

Returns:

fig (matplotlib.figure.Figure) – A matplotlib figure object.

Raises:

ValueError – Raised when no connectivity matrix is found.

Notes

This method uses the function plot_connectivity_circle to plot the connectivity matrix.

compute_harmonic_spectrum_connectivity(sf=None, data=None, precision=0.5, fmin=None, fmax=None, noverlap=1, power_law_remove=False, n_peaks=5, metric='harmsim', n_harms=10, delta_lim=0.1, min_notes=2, plot=False, smoothness=1, smoothness_harm=1, phase_mode=None, save_fig=False, savename='harmonic_spectrum_connectivity.png')[source]#

Computes the harmonic spectrum connectivity between pairs of electrodes. For more details on the harmonic spectrum, see the function compute_cross_spectrum_harmonicity().

Parameters:

  • sf (float, optional) – Sampling frequency of the data. If not provided, uses the object’s own attribute sf.

  • data (ndarray, optional) – 2D data array of shape (n_electrodes, n_datapoints). If not provided, uses the object’s own attribute data.

  • precision (float, default=0.5) – Precision of the frequency axis of the cross-spectrum.

  • fmin (float, optional) – Minimum frequency for computation. If not provided, uses the smallest possible frequency (i.e., zero).

  • fmax (float, optional) – Maximum frequency for computation. If not provided, uses the Nyquist frequency.

  • noverlap (int, default=1) – Argument passed to the scipy.signal.stft() function.

  • power_law_remove (bool, default=False) – If True, removes the power-law noise from the cross-spectrum.

  • n_peaks (int, default=5) – Number of peaks to derive from the harmonic spectrum.

  • metric (str, default=’harmsim’) – Name of the metric to be used in the computation of harmonicity. Choose between:

    • ‘harmsim’

    • ‘subharm_tension’

  • n_harms (int, default=10) – Number of harmonics to consider in the computation.

  • delta_lim (float, default=0.1) – Threshold for the delta limit used in subharmonic tension computation.

  • min_notes (int, default=2) – Minimum number of notes required for a harmonic pattern to be considered valid, when subharmonic tension is used as the metric.

  • plot (bool, default=False) – If True, plots the results.

  • smoothness (int, default=1) – Smoothness factor of the power spectrum. When smoothness=1, no smoothing is applied.

  • smoothness_harm (int, default=1) – Smoothness factor of the harmonic spectrum. When smoothness_harm=1, no smoothing is applied.

  • phase_mode (str, optional) – If set to ‘weighted’, the phase coupling is weighted by the power of associated frequencies.

  • save_fig (bool, default=False) – If True, saves the resulting plot as a .png file.

  • savename (str, default=’harmonic_spectrum_connectivity.png’) – Name of the .png file to save if save_fig is True.

Returns:

output (pandas.DataFrame) – DataFrame containing the results of the harmonic spectrum connectivity computation.

Notes

The harmonic spectrum connectivity is computed between each pair of electrodes. This is achieved by computing the cross-spectrum harmonicity for each pair of electrodes, and storing the results in a DataFrame. The DataFrame includes indices of the electrode pairs, alongside the results of the harmonicity computation.

get_harm_spectrum_metric_matrix(metric)[source]#

Method to retrieve a matrix of harmonic spectrum metric values between pairs of electrodes.

Parameters:

metric (str) – The specific metric from the harmonic spectrum connectivity data to create the matrix from.

Returns:

matrix (pandas.DataFrame or None) – A pivot table with ‘elec1’ and ‘elec2’ as indices and the provided ‘metric’ as values. If ‘harmonic_spectrum_connectivity’ is None or ‘metric’ is not found in ‘harmonic_spectrum_connectivity’, this method will print an error message and return None.

Notes

This method checks if ‘harmonic_spectrum_connectivity’ exists in the object. If it does, it further checks if ‘metric’ exists in the DataFrame’s columns. If both checks pass, a pivot table is created using ‘elec1’ and ‘elec2’ as indices and ‘metric’ as values.

plot_IMF_correlation_matrix(df, elec1, elec2, variable='pearson', savepath=None)[source]#

Plot a heatmap of correlations between IMFs for specified electrodes.

Parameters:

  • df (pd.DataFrame) – DataFrame containing the correlation data with columns: elec1, elec2, imf1, imf2, and pearson.

  • elec1 (int) – The first electrode of interest.

  • elec2 (int) – The second electrode of interest.

Returns:

None

compute_cross_resonance(signal1: ndarray, signal2: ndarray, sf: float, config: ResonanceConfig | None = None) → CrossResonanceResult[source]#

Cross-channel resonance spectrum H × PC = R between two signals.

The cross-channel analog of biotuner.resonance.compute_resonance(), producing three reducer flavors per factor (1→2 asymmetric, 2→1 asymmetric, symmetrized ‘all’) and a corresponding resonance spectrum for each. With the default ResonanceConfig this reproduces the legacy compute_cross_spectrum_harmonicity() numerics within atol=1e-5 on the snapshot regression set.

Parameters:

  • signal1, signal2 (1-D arrays) – Time-domain signals to compare.

  • sf (float) – Sampling frequency (Hz).

  • config (ResonanceConfig, optional) – If None, the legacy-default config is used (matches the historical compute_cross_spectrum_harmonicity behavior).

Returns:

CrossResonanceResult

class CrossResonanceResult(freqs: ~numpy.ndarray, resonance_spectrum: ~typing.Dict[str, ~numpy.ndarray] = <factory>, factors: ~typing.Dict[str, ~typing.Dict[str, ~numpy.ndarray]] = <factory>, summaries: ~typing.Dict[str, ~typing.Dict[str, ~typing.Any]] = <factory>, peaks: ~typing.Dict[str, ~numpy.ndarray] = <factory>, config: ~biotuner.resonance.orchestrator.ResonanceConfig | None = None, intermediates: ~typing.Dict[str, ~typing.Any] | None = None)[source]#

Bases: object

Output of compute_cross_resonance(). Mirrors biotuner.resonance.ResonanceResult but with three reducer flavors per factor (asymmetric 1→2, asymmetric 2→1, symmetrized ‘all’), reflecting the directionality of cross-channel weighting:

'1to2' — H[i] = (p1[i] * Σⱼ p2[j] * S[i,j]) / (2T) (signal1 at i, signal2 at j) '2to1' — H[i] = (p2[i] * Σⱼ p1[j] * S[i,j]) / (2T) (transposed) 'all' — H[i] = (H_1to2[i] + H_2to1[i]) / 2 (symmetrized)

where T = sum(psd1_clean) + sum(psd2_clean) (legacy normalization preserved for bit-exact reproduction of compute_cross_spectrum_harmonicity).

freqs: ndarray#
resonance_spectrum: Dict[str, ndarray]#
factors: Dict[str, Dict[str, ndarray]]#
summaries: Dict[str, Dict[str, Any]]#
peaks: Dict[str, ndarray]#
config: ResonanceConfig | None = None#
intermediates: Dict[str, Any] | None = None#
property harmonicity_matrix: ndarray#

The N×N cross-channel harmonic-similarity matrix S[i, j].

Symmetric in (i, j). Requires ResonanceConfig(return_intermediates=True).

property phase_coupling_matrix: ndarray#

The N×N cross-channel phase-coupling matrix Φ[i, j].

Each cell is W[i, j] · metric(channel-1 at i, channel-2 at j, n, m). NOT symmetric — Φ[i, j] and Φ[j, i] correspond to different directional weightings. Requires ResonanceConfig(return_intermediates=True).

compute_cross_spectrum_harmonicity(signal1, signal2, precision_hz, fmin=None, fmax=None, noverlap=1, fs=44100, power_law_remove=False, n_peaks=5, metric='harmsim', n_harms=10, delta_lim=0.1, min_notes=2, plot=False, smoothness=1, smoothness_harm=1, phase_mode=None, save_fig=False, save_name='harmonic_spectrum_connectivity.png')[source]#

Compute the cross-spectrum harmonicity between two signals. This function is useful for analyzing the interaction between frequency components of two different signals.

Parameters:

  • signal1, signal2 (array_like) – Input signals.

  • precision_hz (float) – The precision in Hz for the short time Fourier transform.

  • fmin, fmax (float or None) – Minimum and maximum frequency for computing the power spectral density.

  • noverlap (int) – Number of points to overlap between segments for computing the power spectral density.

  • fs (int) – The sampling frequency of the signals.

  • power_law_remove (bool) – If True, power-law noise is removed from the power spectral densities of the signals.

  • n_peaks (int) – The number of peaks to identify in the spectra.

  • metric ({“harmsim”, “subharm_tension”}) – The metric to use for computing the harmonicity between the signals.

  • n_harms (int) – Number of harmonics to consider for computing subharmonic tension.

  • delta_lim (float) – Delta limit for computing subharmonic tension.

  • min_notes (int) – Minimum number of notes for computing subharmonic tension.

  • plot (bool) – If True, the function will plot the cross harmonic and phase coupling spectrum.

  • smoothness (int) – Smoothness of the Fourier transform. Higher values result in a smoother transform.

  • smoothness_harm (int) – Smoothness of the harmonic spectrum. Higher values result in a smoother spectrum.

  • phase_mode ({None, “weighted”}) – If “weighted”, the phase coupling is computed as a weighted sum. If None, it is computed as a simple average.

  • save_fig (bool) – If True, the plot is saved as a png file.

  • save_name (str) – Name of the png file if save_fig is True.

Returns:

DataFrame – A pandas DataFrame with columns representing various computed metrics of the signals, such as spectral flatness, spectral entropy, etc. The DataFrame also contains columns for peak frequencies in the spectra and their harmonic similarities.

wPLI_crossfreq(signal1, signal2, peak1, peak2, sf)[source]#

Computes the weighted phase lag index (wPLI) between two signals in specified frequency bands, centered around provided peak frequencies.

Parameters:

  • signal1 (ndarray) – First signal in time series.

  • signal2 (ndarray) – Second signal in time series.

  • peak1 (float) – Center of the frequency band for the first signal.

  • peak2 (float) – Center of the frequency band for the second signal.

  • sf (float) – Sampling frequency.

Returns:

wPLI (float) – Weighted phase lag index between two signals in specified frequency bands.

Examples

>>> signal1 = np.random.normal(0, 1, 5000)
>>> signal2 = np.random.normal(0, 1, 5000)
>>> peak1 = 10.0
>>> peak2 = 20.0
>>> sf = 100.0
>>> wPLI = wPLI_crossfreq(signal1, signal2, peak1, peak2, sf)
wPLI_multiband(signal1, signal2, freq_bands, sf)[source]#

Computes the weighted phase lag index (wPLI) between two signals for multiple frequency bands.

Parameters:

  • signal1 (ndarray) – First signal in time series.

  • signal2 (ndarray) – Second signal in time series.

  • freq_bands (list of tuple) – List of frequency bands. Each band is represented as a tuple (lowcut, highcut).

  • sf (float) – Sampling frequency.

Returns:

wPLI_values (list of float) – List of wPLI values for each frequency band.

Examples

>>> signal1 = np.random.normal(0, 1, 5000)
>>> signal2 = np.random.normal(0, 1, 5000)
>>> freq_bands = [(8, 12), (13, 30), (30, 70)]
>>> sf = 100.0
>>> wPLI_values = wPLI_multiband(signal1, signal2, freq_bands, sf)
cross_frequency_rrci(signal1, signal2, sfreq, freq_peak1, freq_peak2, bandwidth=1, max_denom=50)[source]#

Computes the Rhythmic Ratio Coupling Index (RRCI) between two signals for cross-frequency.

The function first calculates the rhythmic ratio between two peak frequencies. It then filters the two input signals around a frequency band centered on their respective peak frequencies. Finally, it calculates the rhythmic ratio coupling index (RRCI) for these filtered signals. The RRCI is a measure of how much the rhythms of the two signals, in terms of their phase information, are coupled across different frequencies. In other words, it provides a measure of phase-to-phase coupling across these frequencies.

Parameters:

  • signal1 (ndarray) – First signal in time series.

  • signal2 (ndarray) – Second signal in time series.

  • sfreq (float) – Sampling frequency.

  • freq_peak1 (float) – Peak frequency for the first signal.

  • freq_peak2 (float) – Peak frequency for the second signal.

  • bandwidth (float, default=1) – Frequency bandwidth for filtering the signals.

  • max_denom (int, default=50) – The maximum denominator for the rhythmic ratio.

Returns:

rrci (float) – Rhythmic Ratio Coupling Index between the two signals.

Examples

>>> signal1 = np.random.normal(0, 1, 5000)
>>> signal2 = np.random.normal(0, 1, 5000)
>>> sfreq = 100.0
>>> rrci = cross_frequency_rrci(signal1, signal2, sfreq, 10, 20, 2, 4)
n_m_phase_locking(signal1, signal2, n, m, sf)[source]#

Calculate n:m phase locking between two signals.

Parameters:

  • signal1 – First time series.

  • signal2 – Second time series.

  • n – Multiplicative factor for the first signal.

  • m – Multiplicative factor for the second signal.

  • sf – Sampling frequency of the signals.

Returns:

Mean resultant length (Rn:m) indicating the phase locking value.

rhythmic_ratio_coupling_imaginary(signal1, signal2, rhythmic_ratio, sfreq, freq_band1, freq_band2)[source]#

Computes the Imaginary part of Rhythmic Ratio Coupling between two signals. From the paper : “On cross-frequency phase-phase coupling between theta and gamma oscillations in the hippocampus”

Parameters:

  • signal1 (ndarray) – First signal in time series.

  • signal2 (ndarray) – Second signal in time series.

  • rhythmic_ratio (tuple) – Rhythmic ratio (numerator, denominator).

  • sfreq (float) – Sampling frequency.

  • freq_band1 (tuple) – Frequency band for the first signal (lowcut, highcut).

  • freq_band2 (tuple) – Frequency band for the second signal (lowcut, highcut).

Returns:

imaginary_part (float) – Imaginary part of the Rhythmic Ratio Coupling.

compute_rhythmic_ratio(freq1, freq2, max_denom)[source]#
HilbertHuang1D_nopeaks(data, sf, graph=False, nIMFs=5, min_freq=1, max_freq=45, precision=0.5, bin_spread='log', smooth_sigma=None)[source]#

The Hilbert-Huang transform provides a description of how the energy or power within a signal is distributed across frequency. The distributions are based on the instantaneous frequency and amplitude of a signal.

Parameters:

  • data (array (numDataPoints,)) – Single time series.

  • sf (int) – Sampling frequency.

  • graph (bool, default=False) – Defines if a graph is generated.

  • nIMFs (int, default=5) – Number of intrinsic mode functions (IMFs) to keep when Empirical Mode Decomposition (EMD) is computed.

  • min_freq (float, default=1) – Minimum frequency to consider.

  • max_freq (float, default=80) – Maximum frequency to consider.

  • precision (float, default=0.1) – Value in Hertz corresponding to the minimal step between two frequency bins.

  • bin_spread (str, default=’log’) –

    • ‘linear’

    • ‘log’

  • smooth_sigma (float, default=None) – Sigma value for gaussian smoothing.

Returns:

  • IF (array (numDataPoints,nIMFs)) – instantaneous frequencies associated with each IMF.

  • IP (array (numDataPoints,nIMFs)) – instantaneous power associated with each IMF.

  • IA (array (numDataPoints,nIMFs)) – instantaneous amplitude associated with each IMF.

  • spec (array (nIMFs, nbins)) – Power associated with all bins for each IMF

  • bins (array (nIMFs, nbins)) – Frequency bins for each IMF

EMD_time_resolved_harmonicity(time_series1, time_series2, sf, nIMFs=5, method='harmsim')[source]#

Computes the harmonicity between the instantaneous frequencies (IF) for each point in time between all pairs of corresponding intrinsic mode functions (IMFs).

Parameters:

  • time_series1 (array (numDataPoints,)) – First time series.

  • time_series2 (array (numDataPoints,)) – Second time series.

  • sf (int) – Sampling frequency.

  • nIMFs (int, default=5) – Number of intrinsic mode functions (IMFs) to consider.

Returns:

harmonicity (array (numDataPoints, nIMFs)) – Harmonicity values for each pair of corresponding IMFs.

temporal_correlation_fdr(data)[source]#

Compute the temporal correlation for each pair of electrodes and output a connectivity matrix and a matrix of FDR-corrected p-values. This function is used in the computation of the transitional harmony connectivity.

Parameters:

data (array-like) – An array of shape (electrodes, samples) containing the electrode recordings.

Returns:

  • connectivity_matrix (ndarray) – A connectivity matrix of shape (electrodes, electrodes) with the temporal correlation for each pair of electrodes.

  • fdr_corrected_pvals (ndarray) – A matrix of FDR-corrected p-values of shape (electrodes, electrodes).

Notes

This function calculates the temporal correlation for each pair of electrodes, creating a connectivity matrix. Simultaneously, it calculates a matrix of p-values and corrects for multiple comparisons using the False Discovery Rate (FDR) method.

Contents