Resonance

Modular, strategy-registry-based resonance framework. Builds per-frequency H × PC = R spectra from a signal via swappable harmonic kernels, ratio gates, phase estimators, pairwise coupling metrics, and combine rules — with surrogate normalization and a separate path for higher-order coupling metrics (Phase 3).

Public entry points:

• biotuner.resonance.compute_resonance()

• biotuner.resonance.ResonanceConfig

• biotuner.resonance.ResonanceResult

• biotuner.resonance.HigherOrderResult

• biotuner.resonance.with_surrogate_null()

biotuner.resonance — modular, strategy-registry-based resonance framework.

This package builds the per-frequency resonance spectrum R(f) = H(f) · PC(f) from a single signal via swappable harmonic kernels, ratio kernels, phase estimators, pairwise coupling metrics, persistence (Q) measures, and combine rules — with optional surrogate normalization.

Quick start

The recommended path uses default config (joint-probability PC + n:m ratio gating + IAAFT-friendly cross-channel hooks already wired):

from biotuner.resonance import compute_resonance, ResonanceConfig

# Default config — recommended for new analyses
result = compute_resonance(signal, sf=1000)
# result.factors["H"], result.factors["PC"]
# result.resonance_spectrum  — H · PC
# result.summaries["H"/"PC"/"R"]  — complexity dict per spectrum
# result.peaks  — prominence-detected peak freqs

# To reproduce legacy compute_global_harmonicity bit-exactly:
cfg = ResonanceConfig(psd_normalization="minmax_prob", ...)
result = compute_resonance(signal, sf=1000, config=cfg)

Sister modules

• biotuner.harmonic_spectrum — narrow H-only entry point (compute_harmonic_spectrum).

• biotuner.harmonic_connectivity — cross-channel API: compute_cross_resonance for two signals, harmonic_connectivity(...).compute_cross_resonance_connectivity() for N-channel matrices, compute_cross_resonance_connectivity_zscore() for surrogate-normalized statistical inference.

Public API

• compute_resonance() — main entry point

• ResonanceConfig — all swappable knobs

• ResonanceResult — output dataclass

• HigherOrderResult — for Phase 3 higher-order metrics

• with_surrogate_null() — surrogate-z-scored variant of compute_resonance

• results_to_dataframe() — pack multiple results into a DataFrame for the plotters

• harmonic_spectrum_plot_* — comparison plots over collections of results

Internals

Strategy catalogs live in biotuner.resonance.registry. To add a new kernel, metric, or combine rule, call the appropriate register_* function from your own module; it’ll be discoverable by name via ResonanceConfig.

list_strategies(verbose: bool = True) → Dict[str, Dict[str, Callable]]

Discover the strategies currently registered in the resonance package.

Returns a dict keyed by registry name ('HARMONIC_KERNELS', 'PAIRWISE_COUPLING_METRICS', etc.). When verbose=True (the default), also prints a friendly summary.

Use this as the first stop when starting a new analysis:

>>> from biotuner.resonance import list_strategies
>>> list_strategies()
HARMONIC_KERNELS (2)
  - harmsim
  - subharm_tension
RATIO_KERNELS (1)
  - binary
...

All names returned here are valid for the corresponding ResonanceConfig field (e.g. ResonanceConfig(harmonic_kernel=...)).

class ResonanceConfig(psd_method: ~typing.Literal['welch', 'multitaper'] = 'welch', remove_aperiodic: bool = True, psd_normalization: ~typing.Literal['prob', 'minmax_prob', 'none'] = 'minmax_prob', harmonic_kernel: str = 'harmsim', harmonic_kernel_params: ~typing.Dict[str, ~typing.Any] = <factory>, ratio_kernel: str = 'binary', ratio_kernel_params: ~typing.Dict[str, ~typing.Any] = <factory>, phase_estimator: str = 'stft', phase_estimator_params: ~typing.Dict[str, ~typing.Any] = <factory>, coupling_metric: str = 'nm_plv', coupling_metric_params: ~typing.Dict[str, ~typing.Any] = <factory>, pac_convention: ~typing.Literal['row', 'col', 'symmetrize'] = 'row', higher_order_coupling: str | None = None, higher_order_params: ~typing.Dict[str, ~typing.Any] = <factory>, project_higher_order_to_bins: bool = False, projection_op: ~typing.Literal['max', 'sum', 'complexity_weighted_sum'] = 'complexity_weighted_sum', persistence: str | None = None, persistence_params: ~typing.Dict[str, ~typing.Any] = <factory>, combine: str = 'product', combine_params: ~typing.Dict[str, ~typing.Any] = <factory>, alpha_self: float = 1.0, alpha_partner: float = 1.0, legacy_self_pair_subtract: bool = True, gaussian_smooth_sigma: float = 1.0, detrend: bool = False, rescale_factors_after_detrend: bool = True, precision_hz: float = 1.0, fmin: float = 1.0, fmax: float = 30.0, noverlap: int = 1, smoothness: float = 1.0, n_peaks: int = 5, normalize: bool = True, bandwidth_correction: bool = False, null_model: ~typing.Dict[str, ~typing.Any] | None = None, cross_pc_reducer: ~typing.Literal['count', 'joint', 'joint_2T_count'] = 'joint', cross_use_ratio_kernel: bool = True, return_intermediates: bool = False)

Bases: object

Configuration for compute_resonance(). Plan §4.1.

Default values reproduce the legacy compute_global_harmonicity pipeline so that the snapshot regression test (tests/resonance/test_snapshot.py) passes. New code should override these to opt in to cleaner numerics (e.g. switch psd_normalization to 'prob' and legacy_self_pair_subtract to False).

psd_method: Literal['welch', 'multitaper'] = 'welch'

remove_aperiodic: bool = True

psd_normalization: Literal['prob', 'minmax_prob', 'none'] = 'minmax_prob'

harmonic_kernel: str = 'harmsim'

harmonic_kernel_params: Dict[str, Any]

ratio_kernel: str = 'binary'

ratio_kernel_params: Dict[str, Any]

phase_estimator: str = 'stft'

phase_estimator_params: Dict[str, Any]

coupling_metric: str = 'nm_plv'

coupling_metric_params: Dict[str, Any]

pac_convention: Literal['row', 'col', 'symmetrize'] = 'row'

higher_order_coupling: str | None = None

higher_order_params: Dict[str, Any]

project_higher_order_to_bins: bool = False

projection_op: Literal['max', 'sum', 'complexity_weighted_sum'] = 'complexity_weighted_sum'

persistence: str | None = None

persistence_params: Dict[str, Any]

combine: str = 'product'

combine_params: Dict[str, Any]

alpha_self: float = 1.0

alpha_partner: float = 1.0

legacy_self_pair_subtract: bool = True

gaussian_smooth_sigma: float = 1.0

detrend: bool = False

rescale_factors_after_detrend: bool = True

precision_hz: float = 1.0

fmin: float = 1.0

fmax: float = 30.0

noverlap: int = 1

smoothness: float = 1.0

n_peaks: int = 5

normalize: bool = True

bandwidth_correction: bool = False

null_model: Dict[str, Any] | None = None

cross_pc_reducer: Literal['count', 'joint', 'joint_2T_count'] = 'joint'

cross_use_ratio_kernel: bool = True

return_intermediates: bool = False

class ResonanceResult(freqs: ~numpy.ndarray, resonance_spectrum: ~numpy.ndarray, resonance_spectrum_z: ~numpy.ndarray | None = None, surrogate_mean: ~numpy.ndarray | None = None, surrogate_std: ~numpy.ndarray | None = None, factors: ~typing.Dict[str, ~numpy.ndarray] = <factory>, summaries: ~typing.Dict[str, ~typing.Any] = <factory>, config: ~biotuner.resonance.orchestrator.ResonanceConfig | None = None, peaks: ~typing.Dict[str, ~numpy.ndarray] | None = None, higher_order: ~biotuner.resonance.orchestrator.HigherOrderResult | None = None, participation_spectrum: ~numpy.ndarray | None = None, intermediates: ~typing.Dict[str, ~typing.Any] | None = None)

Bases: object

Plan §4.3. Output of compute_resonance().

freqs: ndarray

resonance_spectrum: ndarray

resonance_spectrum_z: ndarray | None = None

surrogate_mean: ndarray | None = None

surrogate_std: ndarray | None = None

factors: Dict[str, ndarray]

summaries: Dict[str, Any]

config: ResonanceConfig | None = None

peaks: Dict[str, ndarray] | None = None

higher_order: HigherOrderResult | None = None

participation_spectrum: ndarray | None = None

intermediates: Dict[str, Any] | None = None

class HigherOrderResult(method: str, triplets: ~typing.List[~typing.Tuple] | None = None, polyrhythms: ~typing.List[~typing.Tuple] | None = None, coupled_pairs: ~typing.List[~typing.Tuple] | None = None, gplv: float | None = None, singular_vectors: ~typing.Tuple | None = None, summaries: ~typing.Dict[str, float] = <factory>)

Bases: object

Plan §4.4. Populated only when config.higher_order_coupling is set.

method: str

triplets: List[Tuple] | None = None

polyrhythms: List[Tuple] | None = None

coupled_pairs: List[Tuple] | None = None

gplv: float | None = None

singular_vectors: Tuple | None = None

summaries: Dict[str, float]

compute_resonance(signal: ndarray, sf: float, config: ResonanceConfig | None = None, freqs: ndarray | None = None) → ResonanceResult

Compute the resonance spectrum + factor breakdown for a 1-D signal.

With default ResonanceConfig this reproduces the legacy compute_global_harmonicity numerics within atol=1e-6 on the snapshot regression set (tests/resonance/snapshots/).

Returns:

ResonanceResult – freqs, resonance_spectrum, factors={‘H’, ‘PC’}, summaries (added by callers), config, peaks={‘H’, ‘PC’, ‘R’}, optional surrogate fields and higher_order.

with_surrogate_null(signal: ndarray, sf: float, config, *, surr_type: str = 'IAAFT', n: int = 200, correction: Literal['zscore', 'pvalue', 'both'] = 'zscore', parallel: bool = True, rng_seed: int | None = None)

Compute resonance on signal and on n surrogates; z-score the spectrum.

Returns a ResonanceResult with resonance_spectrum_z, surrogate_mean, surrogate_std populated (and summaries['p_value_spectrum'] if correction includes pvalue).

Parameters:

• signal (1-D ndarray)

• sf (sampling frequency (Hz))

• config (ResonanceConfig (the null_model field is ignored to avoid recursion))

• surr_type (surrogate type understood by biotuner.surrogates.generate_surrogate)

• n (number of surrogates)

• correction (‘zscore’ | ‘pvalue’ | ‘both’)

• parallel (if True and joblib is importable, parallelize across surrogates)

• rng_seed (optional seed for reproducibility)

harmonic_spectrum_plot_trial_corr(df_all, df_all_rnd, label1='Brain Signals', label2='Random Signals')

Per-trial correlation between harmonicity and phase-coupling spectra.

Expects df_all and df_all_rnd to have a ‘trial’ column and ‘harmonicity’ / ‘phase_coupling’ array-valued columns.

harmonic_spectrum_plot_freq_corr(df1, df2, mean_phase_coupling=False, label1='Brain Signals', label2='Random Signals', fmin=2, fmax=30, xlim=None)

Per-frequency-bin correlation between harmonicity and phase-coupling.

harmonic_spectrum_plot_avg_corr(df1, df2, label1='Brain Signals', label2='Random Signals')

Scatter mean-harmonicity vs mean-phase-coupling for two groups.

results_to_dataframe(results)

Convert a list of ResonanceResult into the legacy DataFrame format.

Each result becomes one row with columns ‘harmonicity’, ‘phase_coupling’, ‘resonance’, and ‘trial’ (0-indexed).

biotuner.resonance.orchestrator — main entry point for resonance computation.

The orchestrator dispatches strategies named in ResonanceConfig against the registries in biotuner.resonance.registry, runs the per-bin pipeline (harmonic kernel → ratio kernel → phase coupling → combine), and returns a ResonanceResult.

Output-arity invariant (plan §4.4): the per-bin resonance spectrum is built ONLY from pairwise-or-lower factors. Higher-order coupling metrics (bplv / mplv / cf_plm / gpla) run on a separate code path and attach as ResonanceResult.higher_order.

class ResonanceConfig(psd_method: ~typing.Literal['welch', 'multitaper'] = 'welch', remove_aperiodic: bool = True, psd_normalization: ~typing.Literal['prob', 'minmax_prob', 'none'] = 'minmax_prob', harmonic_kernel: str = 'harmsim', harmonic_kernel_params: ~typing.Dict[str, ~typing.Any] = <factory>, ratio_kernel: str = 'binary', ratio_kernel_params: ~typing.Dict[str, ~typing.Any] = <factory>, phase_estimator: str = 'stft', phase_estimator_params: ~typing.Dict[str, ~typing.Any] = <factory>, coupling_metric: str = 'nm_plv', coupling_metric_params: ~typing.Dict[str, ~typing.Any] = <factory>, pac_convention: ~typing.Literal['row', 'col', 'symmetrize'] = 'row', higher_order_coupling: str | None = None, higher_order_params: ~typing.Dict[str, ~typing.Any] = <factory>, project_higher_order_to_bins: bool = False, projection_op: ~typing.Literal['max', 'sum', 'complexity_weighted_sum'] = 'complexity_weighted_sum', persistence: str | None = None, persistence_params: ~typing.Dict[str, ~typing.Any] = <factory>, combine: str = 'product', combine_params: ~typing.Dict[str, ~typing.Any] = <factory>, alpha_self: float = 1.0, alpha_partner: float = 1.0, legacy_self_pair_subtract: bool = True, gaussian_smooth_sigma: float = 1.0, detrend: bool = False, rescale_factors_after_detrend: bool = True, precision_hz: float = 1.0, fmin: float = 1.0, fmax: float = 30.0, noverlap: int = 1, smoothness: float = 1.0, n_peaks: int = 5, normalize: bool = True, bandwidth_correction: bool = False, null_model: ~typing.Dict[str, ~typing.Any] | None = None, cross_pc_reducer: ~typing.Literal['count', 'joint', 'joint_2T_count'] = 'joint', cross_use_ratio_kernel: bool = True, return_intermediates: bool = False)

Bases: object

Configuration for compute_resonance(). Plan §4.1.

Default values reproduce the legacy compute_global_harmonicity pipeline so that the snapshot regression test (tests/resonance/test_snapshot.py) passes. New code should override these to opt in to cleaner numerics (e.g. switch psd_normalization to 'prob' and legacy_self_pair_subtract to False).

psd_method: Literal['welch', 'multitaper'] = 'welch'

remove_aperiodic: bool = True

psd_normalization: Literal['prob', 'minmax_prob', 'none'] = 'minmax_prob'

harmonic_kernel: str = 'harmsim'

harmonic_kernel_params: Dict[str, Any]

ratio_kernel: str = 'binary'

ratio_kernel_params: Dict[str, Any]

phase_estimator: str = 'stft'

phase_estimator_params: Dict[str, Any]

coupling_metric: str = 'nm_plv'

coupling_metric_params: Dict[str, Any]

pac_convention: Literal['row', 'col', 'symmetrize'] = 'row'

higher_order_coupling: str | None = None

higher_order_params: Dict[str, Any]

project_higher_order_to_bins: bool = False

projection_op: Literal['max', 'sum', 'complexity_weighted_sum'] = 'complexity_weighted_sum'

persistence: str | None = None

persistence_params: Dict[str, Any]

combine: str = 'product'

combine_params: Dict[str, Any]

alpha_self: float = 1.0

alpha_partner: float = 1.0

legacy_self_pair_subtract: bool = True

gaussian_smooth_sigma: float = 1.0

detrend: bool = False

rescale_factors_after_detrend: bool = True

precision_hz: float = 1.0

fmin: float = 1.0

fmax: float = 30.0

noverlap: int = 1

smoothness: float = 1.0

n_peaks: int = 5

normalize: bool = True

bandwidth_correction: bool = False

null_model: Dict[str, Any] | None = None

cross_pc_reducer: Literal['count', 'joint', 'joint_2T_count'] = 'joint'

cross_use_ratio_kernel: bool = True

return_intermediates: bool = False

class HigherOrderResult(method: str, triplets: ~typing.List[~typing.Tuple] | None = None, polyrhythms: ~typing.List[~typing.Tuple] | None = None, coupled_pairs: ~typing.List[~typing.Tuple] | None = None, gplv: float | None = None, singular_vectors: ~typing.Tuple | None = None, summaries: ~typing.Dict[str, float] = <factory>)

Bases: object

Plan §4.4. Populated only when config.higher_order_coupling is set.

method: str

triplets: List[Tuple] | None = None

polyrhythms: List[Tuple] | None = None

coupled_pairs: List[Tuple] | None = None

gplv: float | None = None

singular_vectors: Tuple | None = None

summaries: Dict[str, float]

class ResonanceResult(freqs: ~numpy.ndarray, resonance_spectrum: ~numpy.ndarray, resonance_spectrum_z: ~numpy.ndarray | None = None, surrogate_mean: ~numpy.ndarray | None = None, surrogate_std: ~numpy.ndarray | None = None, factors: ~typing.Dict[str, ~numpy.ndarray] = <factory>, summaries: ~typing.Dict[str, ~typing.Any] = <factory>, config: ~biotuner.resonance.orchestrator.ResonanceConfig | None = None, peaks: ~typing.Dict[str, ~numpy.ndarray] | None = None, higher_order: ~biotuner.resonance.orchestrator.HigherOrderResult | None = None, participation_spectrum: ~numpy.ndarray | None = None, intermediates: ~typing.Dict[str, ~typing.Any] | None = None)

Bases: object

Plan §4.3. Output of compute_resonance().

freqs: ndarray

resonance_spectrum: ndarray

resonance_spectrum_z: ndarray | None = None

surrogate_mean: ndarray | None = None

surrogate_std: ndarray | None = None

factors: Dict[str, ndarray]

summaries: Dict[str, Any]

config: ResonanceConfig | None = None

peaks: Dict[str, ndarray] | None = None

higher_order: HigherOrderResult | None = None

participation_spectrum: ndarray | None = None

intermediates: Dict[str, Any] | None = None

compute_resonance(signal: ndarray, sf: float, config: ResonanceConfig | None = None, freqs: ndarray | None = None) → ResonanceResult

Compute the resonance spectrum + factor breakdown for a 1-D signal.

With default ResonanceConfig this reproduces the legacy compute_global_harmonicity numerics within atol=1e-6 on the snapshot regression set (tests/resonance/snapshots/).

Returns:

ResonanceResult – freqs, resonance_spectrum, factors={‘H’, ‘PC’}, summaries (added by callers), config, peaks={‘H’, ‘PC’, ‘R’}, optional surrogate fields and higher_order.

biotuner.resonance.kernels_harmonic — harmonic similarity kernels.

Each kernel takes two frequency arrays and returns an (len(freqs_i), len(freqs_j)) similarity matrix. Phase 1 ships harmsim and subharm_tension (bit-equivalent ports of the legacy biotuner.harmonic_spectrum.harmonicity_matrices() formulas). Phase 2 adds sethares, stolzenburg, harmonic_entropy; Phase 3 adds hopf and lorentzian.

References:

Dyad similarity / harmsim: biotuner native; see biotuner.metrics.dyad_similarity(). Subharmonic tension: see biotuner.metrics.compute_subharmonic_tension().

kernel_harmsim(freqs_i: ndarray, freqs_j: ndarray, *, diagonal: float | None = None, **_unused) → ndarray

Harmonic similarity matrix: dyad_similarity(f_i / f_j).

Bit-equivalent to the legacy biotuner.harmonic_spectrum.harmonicity_matrices() branch for metric='harmsim'. f_j == 0 entries are zero; otherwise the raw (un-reduced) ratio f_i / f_j is passed to dyad_similarity (which performs its own reduction internally).

kernel_subharm_tension(freqs_i: ndarray, freqs_j: ndarray, *, n_harms: int = 10, delta_lim: float = 20, min_notes: int = 2, diagonal: float | None = None, **_unused) → ndarray

1 - subharmonic_tension for each (f_i, f_j) pair.

Bit-equivalent to the legacy harmonicity_matrices branch for metric='subharm_tension'.

biotuner.resonance.kernels_ratio — soft (and binary) n:m ratio gates.

Ratio kernels return a weight in [0, 1] for each (f_i, f_j) pair that gates the phase-coupling matrix. Phase 1 ships the binary 5% legacy gate (with 1:1 fallback); Phase 2 adds the soft Arnold-tongue and Stern-Brocot variants from plan §5.2.

Each kernel returns (W, N, M) where:

W : (Ni, Nj) weights N, M : (Ni, Nj) int arrays giving the best (n, m) pair per cell

Reference for soft Arnold-tongue (Phase 2): Pikovsky, Rosenblum, Kurths 2001 Synchronization Ch. 7.

binary_nm_kernel(freqs_i: ndarray, freqs_j: ndarray, *, max_nm: int = 3, tolerance: float = 0.05, fallback_to_1_1: bool = True, **_unused)

Binary n:m gate (legacy): returns 1 for the best n:m match within tolerance, 0 otherwise.

Mirrors the behavior of legacy biotuner.harmonic_spectrum.get_harmonic_ratio(): iterate over (n, m) with 1 <= n, m <= max_nm, pick the one minimizing |ratio - m/n| / (m/n) if any falls below tolerance. If no match is found and fallback_to_1_1 is True, the pair gets (W=1, N=1, M=1); otherwise (W=0, N=0, M=0).

Returns:

(W, N, M) (tuple of (Ni, Nj) ndarrays)

biotuner.resonance.phase_estimators — instantaneous phase extraction.

Phase estimators produce a (n_freqs, n_times) phase matrix from a 1-D signal. Phase 1 ships stft_phase (the legacy STFT-bin phase). Phase 2 adds hilbert_bandpass_phase and morlet_wavelet_phase from plan §5.3.

Note: STFT-bin phase is not true oscillation phase (it’s affected by leakage). New code should prefer Hilbert-bandpass once available; STFT is retained for backward compatibility only.

stft_phase(signal: ndarray, sf: float, *, precision_hz: float, noverlap: int = 10, smoothness: float = 1, **_unused) → ndarray

STFT-bin phase. Bit-equivalent to legacy compute_phase_values.

Parameters:

• signal (1-D array)

• sf (sampling frequency (Hz))

• precision_hz (frequency precision (Hz); nperseg = int(sf / precision_hz).)

• noverlap (STFT overlap (samples))

• smoothness (divides nperseg; smoothness=1 means nperseg as computed.)

Returns:

ndarray (n_freqs, n_times) – np.angle(Zxx) from scipy.signal.stft.

biotuner.resonance.coupling — phase coupling metrics and per-bin reducers.

Pairwise coupling metrics (arity 2) build a Φ[i,j] matrix that feeds the per-bin phase-coupling spectrum. Higher-order metrics (arity 3, K, survey, state) go in plan Phase 3 and run on a separate code path.

Phase 1 ships four pairwise-symmetric variants on the n:m complex-exponential mean, each emphasizing a different aspect of phase coherence:

nm_plv — |<exp(i*(n*φᵢ - m*φⱼ))>|

Classical phase-locking value. Maximum when phase difference is constant across time. Sensitive to volume conduction (0-lag bias). Ref: Tass et al. 1998 PRL 81:3291.

nm_pli — |<sign(Im(exp(i*(n*φᵢ - m*φⱼ))))>|

Phase-Lag Index. Counts only the SIGN of the imaginary part; zero when phases are exactly aligned or anti-aligned. Robust to instantaneous (0-lag) common reference but discards magnitude information. Ref: Stam, Nolte, Daffertshofer 2007 Hum Brain Mapp 28:1178.

nm_wpli — |<|Im(X)| * sign(Im(X))>| / <|Im(X)|>

Weighted Phase-Lag Index. Same 0-lag robustness as PLI but weights by magnitude of the imaginary part — less variance-biased and more sensitive than PLI on noisy data. Ref: Vinck et al. 2011 NeuroImage 55:1548.

nm_rrci — |Im(<exp(i*(n*φᵢ - m*φⱼ))>)|

Rhythmic Ratio Coupling, imaginary part. Like PLV but discards the real part of the mean exponential, isolating non-zero-lag coupling. Ref: Scheffer-Teixeira & Tort 2016 eLife 5:e20515.

Higher-order metrics (bplv triplet, mplv N-ary, cf_plm survey, gpla state) land in Phase 3 and are NOT valid for ResonanceConfig.coupling_metric — see plan §4.4 (arity contract).

nm_plv(phase_i: ndarray, phase_j: ndarray, n: int, m: int) → float

n:m phase-locking value (Tass 1998): | <exp(i*(n*φᵢ - m*φⱼ))> |.

nm_pli(phase_i: ndarray, phase_j: ndarray, n: int, m: int) → float

n:m phase-lag index (Stam 2007): |<sign(Im(exp(i*(n*φᵢ - m*φⱼ))))>|.

Volume-conduction robust: returns 0 for perfectly synchronous (0-lag) or anti-synchronous (π-lag) phase differences. Counts only sign of Im, so discards magnitude.

nm_wpli(phase_i: ndarray, phase_j: ndarray, n: int, m: int) → float

n:m weighted phase-lag index (Vinck 2011): |<|Im(X)| * sign(Im(X))>| / <|Im(X)|> where X = exp(i*(n*φᵢ - m*φⱼ)).

Same 0-lag robustness as PLI but uses Im magnitude as weight, reducing variance bias and improving noise sensitivity.

nm_rrci(phase_i: ndarray, phase_j: ndarray, n: int, m: int) → float

n:m Rhythmic Ratio Coupling, imaginary (Scheffer-Teixeira & Tort 2016): |Im(<exp(i*(n*φᵢ - m*φⱼ))>)|.

Like PLV but discards the real part of the mean exponential, isolating out-of-phase (non-zero-lag) coupling.

nm_wpli_complex(analytic_i: ndarray, analytic_j: ndarray, n: int = 1, m: int = 1, epsilon: float = 1e-10) → float

Amplitude-weighted n:m wPLI on complex analytic signals.

Computes:

|⟨Im(X * conj(Y))⟩| / ⟨|Im(X * conj(Y))|⟩

where X = |a_i| · exp(i·n·φ_i), Y = |a_j| · exp(i·m·φ_j), and the analytic signals carry both amplitude and phase. For n = m = 1 this is simply |⟨Im(a_i · conj(a_j))⟩| / ⟨|Im(a_i · conj(a_j))|⟩.

This matches the cross-spectrum formula used in legacy compute_cross_spectrum_harmonicity, where analytic_* are STFT coefficients (Zxx) at the relevant frequency bins.

Differs from nm_wpli() in that the latter discards amplitude information (uses only sin(Δφ)), while this variant weights by the instantaneous magnitudes |a_i| · |a_j| — more sensitive when joint high-amplitude epochs carry the coupling signal.

Reference: Vinck et al. 2011 NeuroImage 55:1548 (wPLI); applied to STFT cross-spectrum coefficients.

nm_plv_canonical(phase_i: ndarray, phase_j: ndarray, n: int, m: int) → float

n:m PLV with the Tass 1998 convention.

Note on convention

The legacy biotuner ratio kernel (get_harmonic_ratio / binary_nm_kernel) returns (n, m) such that freq_j / freq_i ≈ m / n. The standard Tass et al. 1998 convention is freq_j / freq_i = n / m, with PLV defined as |<exp(i(n*φᵢ - m*φⱼ))>| — so when n*f_i = m*f_j (a true n:m mode-lock), the phase difference is constant in time and PLV = 1.

The legacy nm_plv applies n*φᵢ - m*φⱼ with the swapped (n, m), which yields a non-stationary phase difference even for perfectly locked harmonics — it measures STFT-phase-progression coherence rather than true n:m phase locking. Bit-exact reproduction of the legacy snapshot preserves this behavior.

This nm_plv_canonical variant internally swaps (n, m) to recover the Tass convention, so it correctly returns 1.0 for perfectly locked harmonic pairs. Use this for new analyses where standard n:m PLV semantics matter.

build_pairwise_coupling_matrix(phase: ndarray, freqs: ndarray, ratio_kernel_fn, ratio_kernel_params: dict, metric_fn=None) → ndarray

Build the (n_freqs, n_freqs) symmetric coupling matrix.

For each upper-triangular pair (i, j), consults ratio_kernel_fn to determine the best (n, m); when the ratio kernel returns W=0, the pair gets coupling=0. With the legacy binary kernel and fallback_to_1_1=True, every pair gets a value (n:m if matched, else 1:1).

Parameters:

• phase ((n_freqs, n_times) phase time series)

• freqs ((n_freqs,) frequency grid)

• ratio_kernel_fn (callable returning (W, N, M) arrays)

• ratio_kernel_params (dict passed to ratio_kernel_fn)

• metric_fn (scalar pairwise metric f(phase_i, phase_j, n, m) -> float.) – Default is nm_plv() for backward compatibility.

build_nm_plv_matrix(phase: ndarray, freqs: ndarray, ratio_kernel_fn, ratio_kernel_params: dict) → ndarray

Backwards-compat alias for build_pairwise_coupling_matrix() with metric_fn=nm_plv.

reduce_matrix_to_spectrum(matrix: ndarray, psd_prob: ndarray, *, normalize: bool = True, legacy_self_pair_subtract: bool = True, alpha_self: float = 1.0, alpha_partner: float = 1.0) → ndarray

Reduce an N×N similarity/coupling matrix to a length-N per-bin spectrum.

Two reduction modes:

legacy_self_pair_subtract=True (DEFAULT, matches legacy compute_phase_spectrum and compute_harmonic_power):

v[i] = p_i * Σ_j (M[i,j] * p_j) - M[i,i] * p_i^2

legacy_self_pair_subtract=False (plan §A.4 — recommended for new code):

v[i] = p_i^alpha_self * Σ_{j ≠ i} (M[i,j] * p_j^alpha_partner) Off-diagonal mask, no subtraction artifact.

Parameters:

• matrix ((N, N) similarity/coupling matrix)

• psd_prob ((N,) probability weights summing to 1)

• normalize (if False, mirrors legacy normalize=False branch (compute_phase_spectrum) – uses an asymmetric formula; compute_harmonic_power uses a row sum of M*p_i*p_j)

• legacy_self_pair_subtract (reproduce legacy diagonal-subtraction quirk)

• alpha_self, alpha_partner (only used in non-legacy mode)

biotuner.resonance.combine — combine rules for stacking factors into resonance.

Each rule takes a list of length-N arrays (the factors: H, PC, optionally Q) and returns a single length-N resonance spectrum. The legacy default is product; plan §5.7 adds geomean, harmmean, min, weighted_log for downstream experimentation.

product(factors, **_unused)

np.prod(factors, axis=0). Legacy default.

geomean(factors, **_unused)

Geometric mean: same zeros as product, but linear-scale interpretation.

harmmean(factors, **_unused)

Harmonic mean: penalizes asymmetry between factors.

min_combine(factors, **_unused)

np.min(factors, axis=0). ‘Bottleneck’ reading.

weighted_log(factors, weights=None, **_unused)

Generalized geometric mean: exp(Σ w_i * log(f_i)).

biotuner.resonance.nulls — surrogate-based null normalization.

Wraps biotuner.surrogates.generate_surrogate() to z-score (and/or p-value) the resonance spectrum and scalar summaries against a surrogate distribution. Off by default; opt in via ResonanceConfig.null_model.

For cross-channel analyses, three null generators are exposed for use with biotuner.harmonic_connectivity.compute_cross_resonance():

phase_randomize_surrogate(signal)

Fourier phase randomization — preserves PSD exactly, destroys phase structure. Permissive null (any signal that has the same PSD passes).

iaaft_surrogate(signal)

Iterated Amplitude-Adjusted Fourier Transform (Schreiber & Schmitz 1996). Preserves BOTH the PSD and the amplitude distribution. Tighter than plain phase randomization — rejects more spurious cross-channel coupling.

time_shuffle_surrogate(signal)

Block-shuffle the time-domain signal. Destroys temporal phase relationships while preserving local PSD structure. Strict null for cross-channel coherence testing.

Plan §5.6 and Appendix A.5.

with_surrogate_null(signal: ndarray, sf: float, config, *, surr_type: str = 'IAAFT', n: int = 200, correction: Literal['zscore', 'pvalue', 'both'] = 'zscore', parallel: bool = True, rng_seed: int | None = None)

Compute resonance on signal and on n surrogates; z-score the spectrum.

Returns a ResonanceResult with resonance_spectrum_z, surrogate_mean, surrogate_std populated (and summaries['p_value_spectrum'] if correction includes pvalue).

Parameters:

• signal (1-D ndarray)

• sf (sampling frequency (Hz))

• config (ResonanceConfig (the null_model field is ignored to avoid recursion))

• surr_type (surrogate type understood by biotuner.surrogates.generate_surrogate)

• n (number of surrogates)

• correction (‘zscore’ | ‘pvalue’ | ‘both’)

• parallel (if True and joblib is importable, parallelize across surrogates)

• rng_seed (optional seed for reproducibility)

phase_randomize_surrogate(signal: ndarray, rng: Generator) → ndarray

Fourier phase randomization. Preserves PSD exactly, destroys phase structure. Most permissive null for cross-channel coupling tests.

iaaft_surrogate(signal: ndarray, rng: Generator, n_iter: int = 100, tol: float = 1e-06) → ndarray

Iterated Amplitude-Adjusted Fourier Transform (Schreiber & Schmitz 1996, PRL 77:635). Preserves BOTH the power spectrum AND the empirical amplitude distribution of the original signal — a tighter null than plain phase randomization.

Algorithm:

1. Start from a random shuffle of the signal.

2. Iterate: (a) replace the FFT magnitudes with the original PSD’s,

   2. rank-match back to the original amplitude distribution.

3. Stop when amplitude distribution converges or n_iter reached.

time_shuffle_surrogate(signal: ndarray, rng: Generator, block_size: int | None = None) → ndarray

Block-shuffle surrogate. Cuts the signal into blocks of length block_size (default len(signal)//20) and randomly reorders them. Destroys long-range temporal phase relationships while preserving local PSD structure within blocks. Strictest cross-channel null.

biotuner.resonance.registry — name → callable lookup for resonance axes.

Each axis (harmonic kernel, ratio kernel, phase estimator, coupling metric, persistence, combine rule, surrogate type) is a dict-based registry mapping short string names to callables. The orchestrator dispatches by name; users configure pipelines by setting strings in ResonanceConfig.

Coupling metrics carry an explicit arity tag (pairwise_symmetric, pairwise_asymmetric, triplet, nary, survey, state) that the orchestrator enforces — only pairwise metrics feed the per-bin resonance spectrum. Higher-order metrics run on a separate code path and are not valid values for ResonanceConfig.coupling_metric.

register_harmonic_kernel(name: str, fn: Callable) → None

register_ratio_kernel(name: str, fn: Callable) → None

register_phase_estimator(name: str, fn: Callable) → None

register_coupling_metric(name: str, fn: Callable, arity: str, input_type: str = 'phase') → None

Register a coupling metric with its arity tag and input-type tag.

arity must be one of:

‘pairwise_symmetric’, ‘pairwise_asymmetric’ (valid for coupling_metric) ‘triplet’, ‘nary’, ‘survey’, ‘state’ (valid for higher_order_coupling)

input_type must be one of:

‘phase’ — fn(phase_i, phase_j, n, m) on real-valued phase angles ‘analytic’ — fn(analytic_i, analytic_j, n, m) on complex analytic signals

register_persistence(name: str, fn: Callable) → None

register_combine_rule(name: str, fn: Callable) → None

register_surrogate_type(name: str, fn: Callable) → None

is_pairwise(metric_name: str) → bool

list_strategies(verbose: bool = True) → Dict[str, Dict[str, Callable]]

Discover the strategies currently registered in the resonance package.

Returns a dict keyed by registry name ('HARMONIC_KERNELS', 'PAIRWISE_COUPLING_METRICS', etc.). When verbose=True (the default), also prints a friendly summary.

Use this as the first stop when starting a new analysis:

>>> from biotuner.resonance import list_strategies
>>> list_strategies()
HARMONIC_KERNELS (2)
  - harmsim
  - subharm_tension
RATIO_KERNELS (1)
  - binary
...

All names returned here are valid for the corresponding ResonanceConfig field (e.g. ResonanceConfig(harmonic_kernel=...)).

biotuner.resonance.plots — comparison plots over collections of resonance results.

These three plotters consume DataFrames containing the full H/PC/R triple — the output of stacking multiple compute_resonance() (or legacy compute_global_harmonicity) calls. They were moved here from biotuner.harmonic_spectrum because they visualize the FULL framework, not just the harmonicity factor.

Expected DataFrame columns (any of the following per row):

• ‘harmonicity’ : 1-D ndarray, the H(f) spectrum

• ‘phase_coupling’ : 1-D ndarray, the PC(f) spectrum

• ‘trial’ : trial index (for plot_trial_corr)

A convenience constructor that builds the expected DataFrame from a list of ResonanceResult objects lives at results_to_dataframe().

results_to_dataframe(results)

Convert a list of ResonanceResult into the legacy DataFrame format.

Each result becomes one row with columns ‘harmonicity’, ‘phase_coupling’, ‘resonance’, and ‘trial’ (0-indexed).

harmonic_spectrum_plot_trial_corr(df_all, df_all_rnd, label1='Brain Signals', label2='Random Signals')

Per-trial correlation between harmonicity and phase-coupling spectra.

Expects df_all and df_all_rnd to have a ‘trial’ column and ‘harmonicity’ / ‘phase_coupling’ array-valued columns.

harmonic_spectrum_plot_freq_corr(df1, df2, mean_phase_coupling=False, label1='Brain Signals', label2='Random Signals', fmin=2, fmax=30, xlim=None)

Per-frequency-bin correlation between harmonicity and phase-coupling.

harmonic_spectrum_plot_avg_corr(df1, df2, label1='Brain Signals', label2='Random Signals')

Scatter mean-harmonicity vs mean-phase-coupling for two groups.