Harmonic spectrum#

biotuner.harmonic_spectrum — harmonicity-field computations from PSDs.

Module type: Functions

Resonance values, phase coherence, harmonic power, harmonicity matrices, and peak-finding utilities operating on power spectral densities.

compute_global_harmonicity(signal, precision_hz, fmin=1, fmax=30, noverlap=1, fs=1000, power_law_remove=True, n_peaks=5, metric='harmsim', n_harms=10, delta_lim=20, min_notes=2, plot=False, smoothness=1, smoothness_harm=1, save=False, savename='', phase_mode=None, normalize=True, bandwidth_correction=False, detrend_harmonicity=False, return_fig=False)[source]#

Compute global harmonicity, phase coupling, and resonance spectrum from a signal.

This function computes the Power Spectral Density (PSD) of the signal, applies power law removal if required, calculates the phase matrix, computes dyad similarities and phase couplings for each pairs of frequencies, calculates harmonicity, phase coupling and resonance spectrum, identifies spectral peaks, and returns a dataframe summarizing these metrics.

Parameters:

signal (array_like) – 1-D input signal.
precision_hz (float) – Frequency precision for computing spectra.
fmin (float, default=1) – Minimum frequency to consider in the spectral analysis.
fmax (float, default=30) – Maximum frequency to consider in the spectral analysis.
noverlap (int, default=1) – Number of points of overlap between segments for PSD computation.
fs (int, default=1000) – Sampling frequency.
power_law_remove (bool, default=False) – If True, applies power law removal to the PSD.
n_peaks (int, default=5) – Number of spectral peaks to identify.
metric (str, default=’harmsim’) – Method for computing dyad similarity. Options are:

‘harmsim’ : Harmonic similarity ‘subharm_tension’ : Subharmonic tension
n_harms (int, default=10) – Number of harmonics to consider in dyad similarity computation.
delta_lim (float, default=20) – Limit in ms used when metric is ‘subharm_tension’.
min_notes (int, default=2) – Minimum number of notes for dyad similarity computation.
plot (bool, default=False) – If True, plots the resulting spectra.
smoothness (int, default=1) – Smoothing factor applied to the PSD before computing spectra.
smoothness_harm (int, default=1) – Smoothing factor applied to harmonicity values.
save (bool, default=False) – If True, saves the plot as a .png file.
savename (str, default=’’) – Name for the saved plot file.
phase_mode (str, default=None) – Method for weighting phase coupling computation. Options are ‘weighted’ and ‘None’.
normalize (bool, default=True) – If True, normalize the harmonicity and phase coupling values by dividing by the total power.
bandwidth_correction (bool, default=False) – If True, apply correction for the bias that lower frequencies have more potential harmonic partners within the analysis bandwidth. This corrects the systematic advantage of low frequencies.
detrend_harmonicity (bool, default=False) – If True, remove linear trend from harmonicity, phase coupling, and resonance spectra. This eliminates systematic frequency-dependent bias (e.g., bandwidth effects) while preserving relative structure. Applied after bandwidth_correction if both are enabled.

Returns:

df (DataFrame) – A DataFrame containing computed harmonicity, phase coupling, and resonance values, spectral flatness, entropy, Higushi Fractal Dimension, and spectral spread for each of the three spectra (harmonicity, phase coupling, resonance). Also includes average values and maximum values for these metrics, peak frequencies for each spectrum, and ‘harmsim’ values for peak frequencies.
harmonicity_matrix (ndarray) – The harmonicity matrix of the signal, which corresponds to the power x power comodulogram weighted by dyad similarity.

compute_phase_values(signal, precision_hz, fs=1000, noverlap=10, smoothness=1)[source]#

Compute the phase matrix of a signal using the Short-Time Fourier Transform (STFT).

Parameters:

signal (ndarray) – Input signal.
precision_hz (float) – Frequency precision.
fs (int) – The sampling frequency of the signal.
noverlap (int) – The number of points of overlap between blocks.
smoothness (float) – A parameter for smoothing the signal.

Returns:

ndarray – The phase matrix of the signal.

harmonicity_matrices(freqs, metric='harmsim', n_harms=5, delta_lim=150, min_notes=2)[source]#

Compute harmonicity matrix of frequencies.

Parameters:

freqs (ndarray) – Array of frequencies.
metric (str, optional) – The metric to compute dyad similarity. Default is ‘harmsim’.
n_harms (int, optional) – The number of harmonics. Default is 5.
delta_lim (int, optional) – The delta limit. Default is 150.
min_notes (int, optional) – The minimum number of notes. Default is 2.

Returns:

tuple of ndarrays – The harmonicity and the phase coupling matrix.

get_harmonic_ratio(freq_i, freq_j, max_n=3, max_m=3, tolerance=0.05)[source]#

Identify if two frequencies are at an n:m harmonic ratio.

Parameters:

freq_i (float) – First frequency.
freq_j (float) – Second frequency.
max_n (int, default=3) – Maximum value for n in n:m ratio.
max_m (int, default=3) – Maximum value for m in n:m ratio.
tolerance (float, default=0.05) – Relative tolerance for ratio matching (5%).

Returns:

tuple or None – (n, m) if frequencies are harmonically related, None otherwise.

compute_nm_plv(phase_i, phase_j, n, m)[source]#

Compute Phase Locking Value for n:m phase locking. Tests if n*φᵢ - m*φⱼ is constant over time.

Parameters:

phase_i (ndarray) – Phase time series for frequency i.
phase_j (ndarray) – Phase time series for frequency j.
n (int) – Multiplier for phase i.
m (int) – Multiplier for phase j.

Returns:

float – PLV value [0, 1] indicating strength of n:m phase coupling.

PLV_comod(phase, freqs, psd_clean)[source]#

Compute pure phase coupling matrix using n:m phase locking detection.

This function measures phase locking independent of harmonic relationships. The n:m ratio detection determines which phase relationship to test, but the output PLV is not weighted by harmonic similarity. This prevents double-counting when combining with harmonicity in the resonance spectrum.

Parameters:

phase (ndarray) – Phase matrix from STFT, shape (n_freqs, n_times).
freqs (ndarray) – Array of frequencies.
psd_clean (ndarray) – Cleaned power spectral density (not used but kept for API consistency).

Returns:

ndarray – Phase coupling matrix (n_freqs, n_freqs) containing pure PLV values [0, 1].

Notes

The algorithm: - For frequencies at detectable n:m ratios: computes n:m PLV - For other frequencies: computes standard 1:1 PLV - No harmonic weighting applied (kept independent for resonance calculation) - PLV measures phase synchronization strength regardless of harmonic structure

count_theoretical_harmonic_partners(freq, fmin, fmax, max_ratio=5)[source]#

Count theoretical maximum number of harmonic partners within bandwidth.

This function counts how many frequencies at simple n:m ratios (up to max_ratio) would fall within the analysis bandwidth [fmin, fmax] for a given frequency. Lower frequencies have more potential partners, creating systematic bias.

Parameters:

freq (float) – The frequency to analyze.
fmin (float) – Minimum frequency of analysis bandwidth.
fmax (float) – Maximum frequency of analysis bandwidth.
max_ratio (int, default=5) – Maximum value for n and m in n:m ratios to consider.

Returns:

int – Number of theoretical harmonic partners within bandwidth.

compute_harmonic_power(freqs, dyad_similarities, psd_clean, normalize=True, bandwidth_correction=False, fmin=1, fmax=30)[source]#

Compute harmonicity as probability-weighted average of dyad similarities.

The harmonicity of frequency i represents the expected harmonic similarity with other frequencies, weighted by the joint probability of observing each frequency pair. This formulation is scale-invariant and not sensitive to spectral shape or total power distribution.

Parameters:

freqs (ndarray) – Array of frequencies.
dyad_similarities (ndarray) – Dyad similarities matrix.
psd_clean (ndarray) – The cleaned Power Spectral Density (PSD).
normalize (bool, default=True) – If True, compute harmonicity as weighted average (recommended). If False, compute as sum of weighted harmonic contributions.
bandwidth_correction (bool, default=False) – If True, apply correction for bandwidth bias where lower frequencies have more potential harmonic partners within the analysis range.
fmin (float, default=1) – Minimum frequency of analysis bandwidth (used for bandwidth correction).
fmax (float, default=30) – Maximum frequency of analysis bandwidth (used for bandwidth correction).

Returns:

tuple of ndarrays – The harmonicity values and the harmonicity matrix.

Notes

The weighting uses psd_prob[i] * psd_prob[j] which requires both frequencies to have high power for high harmonicity contribution, preserving the conceptual model that harmonic relationships need both frequencies to be present with significant power.

Bandwidth correction addresses the systematic bias that lower frequencies have more of their harmonic series within typical analysis windows, leading to artificially inflated harmonicity scores.

compute_phase_spectrum(freqs, phase_coupling_matrix, psd_clean, normalize=True)[source]#

Compute phase coupling spectrum as probability-weighted average.

The phase coupling value for each frequency represents the expected phase coupling strength with other frequencies, weighted by the joint probability of observing each frequency pair.

Parameters:

freqs (ndarray) – Array of frequencies.
phase_coupling_matrix (ndarray) – The phase coupling matrix (n_freqs, n_freqs).
psd_clean (ndarray) – The cleaned Power Spectral Density (PSD).
normalize (bool, default=True) – If True, compute as weighted average (recommended). If False, compute as sum of weighted contributions.

Returns:

ndarray – Phase coupling values for each frequency.

Notes

Uses the same probability-weighted averaging approach as compute_harmonic_power to ensure consistency and scale-invariance.

compute_resonance_values(harmonicity_values, phase_coupling_values)[source]#

Compute resonance values from harmonicity and phase coupling values.

Resonance is the product of harmonicity and phase coupling, representing frequencies where both harmonic structure AND phase synchronization are present.

This function preserves the natural scaling of the input spectra rather than normalizing to [0,1], which would destroy the relationship between harmonicity and phase coupling.

Parameters:

harmonicity_values (ndarray) – Harmonicity values for each frequency.
phase_coupling_values (ndarray) – Phase coupling values for each frequency.

Returns:

ndarray – Resonance values for each frequency.

Notes

The product is computed directly without normalization to preserve the independence of harmonicity and phase coupling measures. Both inputs are already probability-weighted averages, so their product is meaningful.

find_spectral_peaks(values, freqs, n_peaks, prominence_threshold=0.5)[source]#

Identify the prominent spectral peaks in a frequency spectrum.

This function uses the peak prominence to select the most notable peaks, and returns their frequencies and indices. Prominence is a measure of how much a peak stands out due to its intrinsic height and its location relative to other peaks.

Parameters:

values (array_like) – 1-D array of values for the frequency spectrum.
freqs (array_like) – 1-D array of frequencies corresponding to the values in ‘values’.
n_peaks (int) – The number of top prominent peaks to return.
prominence_threshold (float, default=0.5) – The minimum prominence a peak must have to be considered notable. Peaks with a prominence less than this value will be ignored.

Returns:

peak_frequencies (ndarray) – Frequencies of the ‘n_peaks’ most prominent peaks.
prominent_peaks (ndarray) – Indices in ‘values’ and ‘freqs’ of the ‘n_peaks’ most prominent peaks.

Harmonic spectrum

Contents

Harmonic spectrum#