Resonance

Resonance#

Modular, strategy-registry-based resonance framework. Builds the per-frequency H × PC = R spectrum from a single 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).

Quick start#

from biotuner.resonance import compute_resonance, ResonanceConfig

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

Discovery: see what kernels / coupling metrics / combine rules are available by calling list_strategies(). All names returned are valid for the corresponding ResonanceConfig field.

Sister modules#

biotuner.harmonic_spectrum — narrow H-only entry point (biotuner.harmonic_spectrum.compute_harmonic_spectrum()).
biotuner.harmonic_connectivity — cross-channel APIs: compute_cross_resonance for two signals, harmonic_connectivity(...).compute_cross_resonance_connectivity() for N-channel matrices, and ...connectivity_zscore() for surrogate-normalized statistical inference.

Public API#

Main entry points:

compute_resonance() — main entry point
ResonanceConfig — all swappable knobs (see return_intermediates to keep the full N×N matrices on the result)
ResonanceResult — output dataclass; see “Result views” below
HigherOrderResult — for Phase 3 higher-order metrics
with_surrogate_null() — surrogate-z-scored variant of compute_resonance. Default surr_type='IAAFT' (iterated AAFT: preserves PSD + amplitude distribution); also accepts 'phase_randomize', 'time_shuffle', and the generate_surrogate types ('AAFT', 'phase', 'shuffle', colored). Populates per-factor z-scores — see factor_z below.

Result views (on ResonanceResult):

result.factors["H" | "PC"] and result.resonance_spectrum — reduced 1-D spectra, length n_freqs
result.factor_z["H" | "PC" | "R"] — per-frequency surrogate z-scores (populated by with_surrogate_null()), with matching result.factor_surrogate_mean / result.factor_surrogate_std. Use factor_z["PC"] for phase-coupling inference: R is harmonicity-dominated (PSD-driven), so factor_z["R"] is largely blind to phase coupling under a PSD-preserving null. result.resonance_spectrum_z mirrors factor_z["R"].
result.harmonicity_matrix and result.phase_coupling_matrix — full N×N matrices S[i, j] and Φ[i, j] (need ResonanceConfig(return_intermediates=True))
result.summaries["H" | "PC" | "R"] — scalar metrics per spectrum: avg, max, flatness, entropy, spread, higuchi, peaks, peaks_avg, peak_harmsim, peak_harmsim_avg, peak_harmsim_max (one-to-one with legacy compute_global_harmonicity columns)
result.peaks["H" | "PC" | "R"] — peak frequencies per spectrum

Discovery + plotting helpers:

list_strategies() — print/return the registered strategies across all 8 registries
results_to_dataframe() — pack multiple results into a wide pandas DataFrame; flattens summaries into legacy-named columns (harm_flatness, phase_peaks_avg, res_peak_harmsim_max, …) by default; pass flatten_summaries=False for the narrow 4-column form
harmonic_spectrum_plot_avg_corr, harmonic_spectrum_plot_trial_corr, harmonic_spectrum_plot_freq_corr — 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 and listed by list_strategies().

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#

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

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 (2)
  - binary
  - fraction
...

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)[source]#

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, factor_z: ~typing.Dict[str, ~numpy.ndarray] | None = None, factor_surrogate_mean: ~typing.Dict[str, ~numpy.ndarray] | None = None, factor_surrogate_std: ~typing.Dict[str, ~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)[source]#

Bases: object

Plan §4.3. Output of compute_resonance().

Quick reference for getting different views of the analysis:

Reduced 1-D spectra (length n_freqs)
result.factors["H"] — harmonicity spectrum result.factors["PC"] — phase coupling spectrum result.resonance_spectrum — R = combine(H, PC)
Full 2-D matrices (shape n_freqs × n_freqs)
result.harmonicity_matrix — S[i, j] harmonic similarity result.phase_coupling_matrix — Φ[i, j] phase coupling These need config.return_intermediates=True — the matrix properties raise a clear error otherwise.
Scalar metrics per spectrum
result.summaries["H" | "PC" | "R"] → dict with avg, max, peaks, peaks_avg, flatness, entropy, spread, higuchi, peak_harmsim, peak_harmsim_avg, peak_harmsim_max. Match the legacy compute_global_harmonicity columns one-to-one.
Peak frequencies per spectrum
result.peaks["H" | "PC" | "R"] → ndarray of peak frequencies (the same array that lives at summaries[...]["peaks"]).

For a flat pandas DataFrame of multiple results, see biotuner.resonance.results_to_dataframe().

freqs: ndarray#

resonance_spectrum: ndarray#

resonance_spectrum_z: ndarray | None = None#

surrogate_mean: ndarray | None = None#

surrogate_std: ndarray | None = None#

factor_z: Dict[str, ndarray] | None = None#

factor_surrogate_mean: Dict[str, ndarray] | None = None#

factor_surrogate_std: Dict[str, 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#

property harmonicity_matrix: ndarray#

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

Requires ResonanceConfig(return_intermediates=True). Off-diagonal cells encode harmonic similarity between frequency bins i and j (high where the ratio is musically simple). Row-summing S * p[j] and scaling by p[i] gives the reduced factors["H"] spectrum.

property phase_coupling_matrix: ndarray#

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

Requires ResonanceConfig(return_intermediates=True). Each cell is W[i, j] · metric(phase_i, phase_j, n, m) where (n, m) comes from the ratio kernel and metric from the coupling metric. Row-summing gives the reduced factors["PC"] spectrum.

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>)[source]#

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[source]#

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)[source]#

Compute resonance on signal and on n surrogates; z-score every factor.

Returns a ResonanceResult with, for the resonance spectrum R: resonance_spectrum_z, surrogate_mean, surrogate_std — and, for ALL three factors, factor_z / factor_surrogate_mean / factor_surrogate_std dicts keyed "H"/"PC"/"R".

Because R is H-dominated (H is PSD-driven), R_z is largely blind to phase coupling under a PSD-preserving null; factor_z["PC"] is the correct detector for n:m phase coupling. See the resonance_paper validation suite.

Parameters:

signal (1-D ndarray)
sf (sampling frequency (Hz))
config (ResonanceConfig (the null_model field is ignored to avoid recursion))
surr_type ('IAAFT' (default; iterated AAFT, preserves PSD + amplitude) – distribution), 'phase_randomize', 'time_shuffle', or any type accepted by biotuner.surrogates.generate_surrogate() ('AAFT', 'TFT', 'phase', 'shuffle', 'white'/'pink'/'brown'/'blue').
n (number of surrogates)
correction (‘zscore’ | ‘pvalue’ | ‘both’ (p-values added per factor when) – pvalue is requested: summaries['p_value_H'|'p_value_PC'|'p_value_spectrum'])
parallel (if True and joblib is importable, parallelize across surrogates)
rng_seed (optional seed for reproducibility)

nm_intertrial_plv(phase_epochs_i: ndarray, phase_epochs_j: ndarray, n: int = 1, m: int = 1) → float[source]#

Inter-trial n:m phase-locking value across epochs (Tass convention).

Inputs are 2-D epoched phase arrays (n_epochs, n_times) for the two frequencies. For each epoch the time-averaged n:m relative-phase resultant < exp(i*(m*φ_i - n*φ_j)) >_t is formed (note the m/n swap, matching nm_plv_canonical()); the inter-trial PLV is the magnitude of the mean of those per-epoch resultants across epochs:

ITC = | (1/E) Σ_e  < exp(i*(m*φ_i[e] - n*φ_j[e])) >_t |

This is the correct estimator when coupling is consistent across trials but the absolute phase resets between trials — the regime where a continuous single-trial PLV cannot distinguish genuine coupling from a stationary process (see resonance_paper Study 1B / Study 5). It is provided as a standalone utility because it requires an epoch dimension that the single-trial orchestrator does not model; it is therefore NOT registered as a PAIRWISE_COUPLING_METRIC (those receive 1-D phase from the orchestrator).

Parameters:

phase_epochs_i, phase_epochs_j (ndarray (n_epochs, n_times)) – Instantaneous phase per epoch at the two frequencies (e.g. from the hilbert phase estimator applied per epoch).
n, m (int) – n:m ratio. For 1:1 this is the standard inter-trial coherence of the phase difference.

Returns:

float in [0, 1] — 1 = perfectly trial-consistent n (m relative phase.)

detect_nm_coupling(a, b, sf, freq_pairs, *, metrics=('nm_plv', 'nm_rho_entropy', 'nm_phase_mi'), bandwidth=3.0, max_denom=16, n_surrogates=99, seed=0, scope_guard=True)[source]#

Detect n:m phase coupling between signals a and b at the given frequency pairs.

Parameters:

a, b (1-D arrays — the two signals (ideally INDEPENDENT sources; see scope note).)
sf (sampling frequency (Hz).)
freq_pairs (list of (f_a, f_b) — the component frequencies to test (e.g. (10, 15) for 2:3).)
metrics (panel of technique names (subset of PANEL). Default = plv + rho_entropy + phase_mi.)
bandwidth (Hz width of the band-pass around each frequency (>=3 keeps Hilbert phase stable).)
n_surrogates (number of IAAFT-of-b surrogates for the null.)
scope_guard (if True, warn when a and b are highly correlated (within-signal / shared source).)

Returns:

dict with ``results`` (one entry per freq pair (resolved n:m, and per-metric value / surrogate z /)
rank-p) and warning (the scope-guard message, or None).

nm_multipliers(f_a, f_b, max_denom=16)[source]#

Correct (Tass) integer multipliers (n, m) for frequencies f_a, f_b: n*f_a = m*f_b.

n multiplies phi_a, m multiplies phi_b. For f_b/f_a = p/q (lowest terms) this returns (n=p, m=q) so that n*phi_a - m*phi_b is stationary for a genuine lock.

harmonic_spectrum_plot_trial_corr(df_all, df_all_rnd, label1='Brain Signals', label2='Random Signals')[source]#

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)[source]#: Per-frequency-bin correlation between harmonicity and phase-coupling.

harmonic_spectrum_plot_avg_corr(df1, df2, label1='Brain Signals', label2='Random Signals')[source]#: Scatter mean-harmonicity vs mean-phase-coupling for two groups.

results_to_dataframe(results, *, flatten_summaries: bool = True)[source]#

Convert a list of ResonanceResult into a wide DataFrame.

By default, expands result.summaries into one column per {spectrum}_{metric} pair so the returned frame matches the pre-refactor compute_global_harmonicity column layout.

Columns always present (one row per result):

trial              — 0-indexed
freqs              — 1-D ndarray, the frequency grid
harmonicity        — 1-D ndarray, the H(f) spectrum
phase_coupling     — 1-D ndarray, the PC(f) spectrum
resonance          — 1-D ndarray, the R(f) spectrum

Additional columns when flatten_summaries=True (default) — for each spectrum s ∈ {'harm', 'phase', 'res'} and each metric exposed by ResonanceResult.summaries:

{s}_avg                {s}_max               {s}_peaks
{s}_peak_indices       {s}_peaks_avg         {s}_flatness
{s}_entropy            {s}_spread            {s}_higuchi
{s}_peak_harmsim       {s}_peak_harmsim_avg  {s}_peak_harmsim_max

The s prefix follows the legacy DataFrame convention (harm_*, phase_*, res_*) rather than the result’s internal keys (H, PC, R).

Parameters:

results (list of ResonanceResult)
flatten_summaries (bool, default=True) – If False, return only the bare 5 columns (back-compat).

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.

Bases: object

Configuration for compute_resonance(). Plan §4.1.

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#

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

Bases: object

Plan §4.3. Output of compute_resonance().

Quick reference for getting different views of the analysis:

Reduced 1-D spectra (length n_freqs)
result.factors["H"] — harmonicity spectrum result.factors["PC"] — phase coupling spectrum result.resonance_spectrum — R = combine(H, PC)
Full 2-D matrices (shape n_freqs × n_freqs)
result.harmonicity_matrix — S[i, j] harmonic similarity result.phase_coupling_matrix — Φ[i, j] phase coupling These need config.return_intermediates=True — the matrix properties raise a clear error otherwise.
Scalar metrics per spectrum
result.summaries["H" | "PC" | "R"] → dict with avg, max, peaks, peaks_avg, flatness, entropy, spread, higuchi, peak_harmsim, peak_harmsim_avg, peak_harmsim_max. Match the legacy compute_global_harmonicity columns one-to-one.
Peak frequencies per spectrum
result.peaks["H" | "PC" | "R"] → ndarray of peak frequencies (the same array that lives at summaries[...]["peaks"]).

For a flat pandas DataFrame of multiple results, see biotuner.resonance.results_to_dataframe().

freqs: ndarray#

resonance_spectrum: ndarray#

resonance_spectrum_z: ndarray | None = None#

surrogate_mean: ndarray | None = None#

surrogate_std: ndarray | None = None#

factor_z: Dict[str, ndarray] | None = None#

factor_surrogate_mean: Dict[str, ndarray] | None = None#

factor_surrogate_std: Dict[str, 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#

property harmonicity_matrix: ndarray#

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

property phase_coupling_matrix: ndarray#

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

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

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[source]#

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[source]#

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 — n:m ratio gates for phase coupling.

Ratio kernels decide, for each frequency pair (f_i, f_j), which integer ratio (n, m) to test for phase coupling and how much weight to give the result. The orchestrator’s PC pipeline calls these to drive build_pairwise_coupling_matrix.

Each kernel returns (W, N, M) where:: W : (Ni, Nj) weights in [0, 1] N, M : (Ni, Nj) int arrays giving the (n, m) pair to test per cell

Convention (shared by all registered kernels)#

All ratio kernels return (n, m) such that ratio = f_j / f_i ≈ m / n. This is the legacy biotuner convention, preserved across all kernels for consistency. To get a mathematically correct n:m phase-locking test (Tass 1998 form n·φ_i − m·φ_j = const), ALWAYS pair these kernels with the coupling metric nm_plv_canonical (which swaps internally). Pairing with the raw nm_plv keeps the legacy behavior preserved for bit-exact paper reproduction.

Quick rule#

New analyses: coupling_metric='nm_plv_canonical' with any ratio kernel.
Paper reproduction: coupling_metric='nm_plv' with ratio_kernel='binary'.

Registered kernels#

fraction
For each pair, computes Fraction(f_j / f_i).limit_denominator(max_denom) to get the EXACT closest rational. Works for any frequency pair — e.g. (10 Hz, 17 Hz) gives (n=10, m=17), testing the actual 10:17 mode-lock. Weight W = exp(-beta * log2(n*m)) penalizes high-order ratios.

binary (DEFAULT)
Legacy gate (preserves bit-exact reproduction): tries (n, m) pairs with 1 ≤ n, m ≤ max_nm (default 3) and picks the best match within tolerance. Returns W=1 if any match found, W=0 otherwise (or W=1 at (1,1) if fallback_to_1_1=True). Misses coupling at any ratio outside the small preset table — use ‘fraction’ for new analyses.

Phase 2 will add ‘arnold_tongue’ (soft Gaussian membership of Arnold tongues, Pikovsky-Rosenblum-Kurths 2001) and ‘stern_brocot’ (depth-weighted complexity).

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

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)

fraction_kernel(freqs_i: ndarray, freqs_j: ndarray, *, max_denom: int = 16, beta: float = 1.0, **_unused)[source]#

For each pair, picks (n, m) as the exact closest rational of f_j/f_i.

For ANY frequency pair this returns a meaningful (n, m) — the legacy binary_nm_kernel could only handle ratios up to max_nm=3.

Convention#

Returns (n, m) such that ratio = f_j / f_i ≈ m / n — the SAME convention as binary_nm_kernel(). This is the legacy biotuner convention; to get a mathematically correct n:m phase-locking test, pair this kernel with coupling_metric='nm_plv_canonical' (which swaps internally to apply the Tass 1998 convention).

The weight W penalizes high-order ratios via Tenney height:

W[i, j] = exp(-beta * log2(n * m))

Examples (max_denom=16, beta=1.0):: (10, 20) → (n=1, m=2), W = 0.368 octave (10, 15) → (n=2, m=3), W = 0.075 perfect fifth (10, 17) → (n=10, m=17), W ≈ 6e-4 complex, but exact test (10, 14.14) → (n=12, m=17), W ≈ 5e-4 closest rational to √2

The recommended pairing for new analyses is:

ResonanceConfig(
    ratio_kernel='fraction',
    coupling_metric='nm_plv_canonical',   # not just 'nm_plv'!
)

Parameters:

freqs_i, freqs_j (1-D arrays)
max_denom (int, default=16) – Maximum denominator passed to Fraction.limit_denominator. Larger values give more exact ratios; smaller values force simpler approximations. The legacy cross_frequency_rrci used 16 as well.
beta (float, default=1.0) – Tenney-height complexity penalty exponent. beta=0 gives W=1 everywhere; beta=1 matches typical musicological assumptions.

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, with row ``i`` aligned to ``freqs[i]`` (the analysis frequency grid passed by the orchestrator). Two estimators are registered:

stft (default)
STFT-bin phase. Bit-equivalent to legacy compute_phase_values for the raw transform, but now the returned rows are selected to match the analysis freqs grid. (Historically the full 0..Nyquist STFT grid was indexed with the [fmin, fmax]-clipped freqs, so phase[i] was off by fmin — that bug corrupted every phase-coupling entry. Passing freqs here fixes it.) STFT-bin phase is still leakage-limited and is not true oscillation phase.

hilbert
Hilbert analytic phase of the signal band-pass-filtered around each analysis frequency. This is true instantaneous phase and recovers n:m phase locking that the STFT-bin estimator misses (validated in resonance_paper Study 5). Aligned to freqs by construction.

Both accept and ignore extra kwargs so the orchestrator can pass a common set.

stft_phase(signal: ndarray, sf: float, *, precision_hz: float, noverlap: int = 10, smoothness: float = 1, freqs: ndarray | None = None, **_unused) → ndarray[source]#

STFT-bin phase, with rows aligned to freqs when provided.

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.)
freqs (analysis frequency grid. When given, the returned phase rows are the) – STFT bins nearest each freqs value, so phase[i] is the phase at freqs[i]. When None, the full 0..Nyquist grid is returned (legacy).

Returns:

ndarray (n_freqs, n_times) of np.angle(Zxx).

hilbert_bandpass_phase(signal: ndarray, sf: float, *, freqs: ndarray, bandwidth: float = 2.0, filter_order: int = 4, **_unused) → ndarray[source]#

Hilbert analytic phase of the signal band-passed around each freqs[i].

Returns true instantaneous phase at sample resolution, aligned to freqs by construction. For a band centered at f the passband is [f - bandwidth/2, f + bandwidth/2] (clamped to (0, Nyquist)).

Returns:

ndarray (n_freqs, n_times) of unwrapped-then-wrapped instantaneous phase
angles (np.angle of the analytic signal).

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[source]#: n:m phase-locking value (Tass 1998): | <exp(i*(n*φᵢ - m*φⱼ))> |.

nm_pli(phase_i: ndarray, phase_j: ndarray, n: int, m: int) → float[source]#

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[source]#

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[source]#

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[source]#

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[source]#: 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.

nm_intertrial_plv(phase_epochs_i: ndarray, phase_epochs_j: ndarray, n: int = 1, m: int = 1) → float[source]#

Inter-trial n:m phase-locking value across epochs (Tass convention).

ITC = | (1/E) Σ_e  < exp(i*(m*φ_i[e] - n*φ_j[e])) >_t |

Parameters:

phase_epochs_i, phase_epochs_j (ndarray (n_epochs, n_times)) – Instantaneous phase per epoch at the two frequencies (e.g. from the hilbert phase estimator applied per epoch).
n, m (int) – n:m ratio. For 1:1 this is the standard inter-trial coherence of the phase difference.

Returns:

float in [0, 1] — 1 = perfectly trial-consistent n (m relative phase.)

nm_rho_entropy(phase_i, phase_j, n, m, nbins: int = 18)[source]#

Tass 1998 n:m entropy synchronization index in [0, 1].

rho = (Hmax - H)/Hmax where H is the Shannon entropy of the histogram of psi = n*phi_i - m*phi_j (wrapped) and Hmax = ln(nbins). Reads ANY departure of psi from uniformity (all moments), so it detects multimodal n:m locks the PLV family misses. Positively biased at small N — use a surrogate z. Ref: Tass et al. 1998 PRL 81:3291.

nm_conditional_prob(phase_i, phase_j, n, m, nbins: int = 18, min_count: int = 5)[source]#

Tass 1998 n:m conditional-probability index in [0, 1].

Bins n*phi_i and takes the mean resultant of m*phi_j within each bin. Captures dependence of one phase given the other, but each bin’s resultant is again a first moment — so it shares the PLV family’s blindness to within-bin antipodal structure. Ref: Tass et al. 1998 PRL 81:3291.

nm_phase_mi(phase_i, phase_j, n, m, nbins: int = 16)[source]#

Normalized mutual information in [0, 1] between n*phi_i and m*phi_j.

Model-free, all-moment statistical dependence — the most general n:m detector; recovers any multimodal/nonlinear lock. Most data-hungry; use a surrogate z. Ref: Palus 1997 Phys Lett A 235:341.

build_pairwise_coupling_matrix(phase: ndarray, freqs: ndarray, ratio_kernel_fn, ratio_kernel_params: dict, metric_fn=None) → ndarray[source]#

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[source]#: 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[source]#

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)

nm_pli_canonical(arg_i, arg_j, n, m, **kw)#: Canonical (Tass-convention) variant of nm_pli() — swaps (n, m) so the genuine n:m mode-lock is tested with the legacy ratio kernels.

nm_wpli_canonical(arg_i, arg_j, n, m, **kw)#: Canonical (Tass-convention) variant of nm_wpli() — swaps (n, m) so the genuine n:m mode-lock is tested with the legacy ratio kernels.

nm_rrci_canonical(arg_i, arg_j, n, m, **kw)#: Canonical (Tass-convention) variant of nm_rrci() — swaps (n, m) so the genuine n:m mode-lock is tested with the legacy ratio kernels.

nm_wpli_complex_canonical(arg_i, arg_j, n, m, **kw)#: Canonical (Tass-convention) variant of nm_wpli_complex() — swaps (n, m) so the genuine n:m mode-lock is tested with the legacy ratio kernels.

nm_rho_entropy_canonical(arg_i, arg_j, n, m, **kw)#: Canonical (Tass-convention) variant of nm_rho_entropy() — swaps (n, m) so the genuine n:m mode-lock is tested with the legacy ratio kernels.

nm_conditional_prob_canonical(arg_i, arg_j, n, m, **kw)#: Canonical (Tass-convention) variant of nm_conditional_prob() — swaps (n, m) so the genuine n:m mode-lock is tested with the legacy ratio kernels.

nm_phase_mi_canonical(arg_i, arg_j, n, m, **kw)#: Canonical (Tass-convention) variant of nm_phase_mi() — swaps (n, m) so the genuine n:m mode-lock is tested with the legacy ratio kernels.

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)[source]#: np.prod(factors, axis=0). Legacy default.

geomean(factors, **_unused)[source]#: Geometric mean: same zeros as product, but linear-scale interpretation.

harmmean(factors, **_unused)[source]#: Harmonic mean: penalizes asymmetry between factors.

min_combine(factors, **_unused)[source]#: np.min(factors, axis=0). ‘Bottleneck’ reading.

weighted_log(factors, weights=None, **_unused)[source]#: 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.

Compute resonance on signal and on n surrogates; z-score every factor.

Parameters:

signal (1-D ndarray)
sf (sampling frequency (Hz))
config (ResonanceConfig (the null_model field is ignored to avoid recursion))
surr_type ('IAAFT' (default; iterated AAFT, preserves PSD + amplitude) – distribution), 'phase_randomize', 'time_shuffle', or any type accepted by biotuner.surrogates.generate_surrogate() ('AAFT', 'TFT', 'phase', 'shuffle', 'white'/'pink'/'brown'/'blue').
n (number of surrogates)
correction (‘zscore’ | ‘pvalue’ | ‘both’ (p-values added per factor when) – pvalue is requested: summaries['p_value_H'|'p_value_PC'|'p_value_spectrum'])
parallel (if True and joblib is importable, parallelize across surrogates)
rng_seed (optional seed for reproducibility)

phase_randomize_surrogate(signal: ndarray, rng: Generator) → ndarray[source]#: 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[source]#

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:

Start from a random shuffle of the signal.
Iterate: (a) replace the FFT magnitudes with the original PSD’s,
1. rank-match back to the original amplitude distribution.
Stop when amplitude distribution converges or n_iter reached.

time_shuffle_surrogate(signal: ndarray, rng: Generator, block_size: int | None = None) → ndarray[source]#: 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[source]#

register_ratio_kernel(name: str, fn: Callable) → None[source]#

register_phase_estimator(name: str, fn: Callable) → None[source]#

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

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[source]#

register_combine_rule(name: str, fn: Callable) → None[source]#

register_surrogate_type(name: str, fn: Callable) → None[source]#

is_pairwise(metric_name: str) → bool[source]#

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

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 (2)
  - binary
  - fraction
...

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, *, flatten_summaries: bool = True)[source]#

Convert a list of ResonanceResult into a wide DataFrame.

By default, expands result.summaries into one column per {spectrum}_{metric} pair so the returned frame matches the pre-refactor compute_global_harmonicity column layout.

Columns always present (one row per result):

trial              — 0-indexed
freqs              — 1-D ndarray, the frequency grid
harmonicity        — 1-D ndarray, the H(f) spectrum
phase_coupling     — 1-D ndarray, the PC(f) spectrum
resonance          — 1-D ndarray, the R(f) spectrum

Additional columns when flatten_summaries=True (default) — for each spectrum s ∈ {'harm', 'phase', 'res'} and each metric exposed by ResonanceResult.summaries:

{s}_avg                {s}_max               {s}_peaks
{s}_peak_indices       {s}_peaks_avg         {s}_flatness
{s}_entropy            {s}_spread            {s}_higuchi
{s}_peak_harmsim       {s}_peak_harmsim_avg  {s}_peak_harmsim_max

The s prefix follows the legacy DataFrame convention (harm_*, phase_*, res_*) rather than the result’s internal keys (H, PC, R).

Parameters:

results (list of ResonanceResult)
flatten_summaries (bool, default=True) – If False, return only the bare 5 columns (back-compat).

harmonic_spectrum_plot_trial_corr(df_all, df_all_rnd, label1='Brain Signals', label2='Random Signals')[source]#

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)[source]#: Per-frequency-bin correlation between harmonicity and phase-coupling.

harmonic_spectrum_plot_avg_corr(df1, df2, label1='Brain Signals', label2='Random Signals')[source]#: Scatter mean-harmonicity vs mean-phase-coupling for two groups.

Resonance

Contents

Resonance#

Quick start#

Sister modules#

Public API#

Internals#

Quick start#

Sister modules#

Public API#

Internals#

Convention (shared by all registered kernels)#

Quick rule#

Registered kernels#

Convention#

Note on convention#