noda

crosspeak.noda

MovingWindowResult

Bases: NamedTuple

Correlation maps from a moving-window 2DCOS analysis.

Attributes:

Name	Type	Description
sync	ndarray	Synchronous matrices, shape (n_windows, n_wavenumbers, n_wavenumbers).
asyn	ndarray	Asynchronous matrices, same shape.
centres	ndarray	Perturbation value at the middle of each window, shape (n_windows,).

Source code in src/crosspeak/noda.py

class MovingWindowResult(NamedTuple):
    """Correlation maps from a moving-window 2DCOS analysis.

    Attributes
    ----------
    sync
        Synchronous matrices, shape (n_windows, n_wavenumbers, n_wavenumbers).
    asyn
        Asynchronous matrices, same shape.
    centres
        Perturbation value at the middle of each window, shape (n_windows,).
    """

    sync: np.ndarray
    asyn: np.ndarray
    centres: np.ndarray

AutopeakResult

Bases: NamedTuple

Autopeaks found on a synchronous diagonal.

Attributes:

Name	Type	Description
positions	ndarray	Wavenumbers at which autopeaks were detected, in the order of the input wavenumber axis.
intensities	ndarray	Autopower values (diagonal of the synchronous matrix) at those positions.

Source code in src/crosspeak/noda.py

class AutopeakResult(NamedTuple):
    """Autopeaks found on a synchronous diagonal.

    Attributes
    ----------
    positions
        Wavenumbers at which autopeaks were detected, in the order of the
        input wavenumber axis.
    intensities
        Autopower values (diagonal of the synchronous matrix) at those
        positions.
    """

    positions: np.ndarray
    intensities: np.ndarray

hilbert_noda_matrix(n)

The Hilbert-Noda matrix for n perturbation points.

The discrete kernel that turns the asynchronous calculation into a matrix product: with N in hand, Ψ = Y^T @ N @ Y / (m - 1). Off the diagonal N_jk = 1 / (π(k - j)); the diagonal is zero, since a point has no phase lag against itself.

The sign convention matters. With N positive in the upper triangle, a positive Ψ(ν₁, ν₂) alongside a positive Φ means ν₁ leads ν₂ under Noda's rules. Transpose N and every ordering flips — so if this ever gets swapped out, check the sign against a sequence you already know.

Parameters:

Name	Type	Description	Default
n		Number of perturbation points (m in the usual notation), at least 2.	required

Returns:

Type	Description
ndarray	The (n, n) matrix, antisymmetric with a zero diagonal.

Source code in src/crosspeak/noda.py

def hilbert_noda_matrix(n):
    """The Hilbert-Noda matrix for `n` perturbation points.

    The discrete kernel that turns the asynchronous calculation into a matrix
    product: with N in hand, Ψ = Y^T @ N @ Y / (m - 1). Off the diagonal
    N_jk = 1 / (π(k - j)); the diagonal is zero, since a point has no phase
    lag against itself.

    The sign convention matters. With N positive in the upper triangle, a
    positive Ψ(ν₁, ν₂) alongside a positive Φ means ν₁ leads ν₂ under Noda's
    rules. Transpose N and every ordering flips — so if this ever gets swapped
    out, check the sign against a sequence you already know.

    Parameters
    ----------
    n
        Number of perturbation points (m in the usual notation), at least 2.

    Returns
    -------
    np.ndarray
        The (n, n) matrix, antisymmetric with a zero diagonal.
    """
    if not isinstance(n, (int, np.integer)):
        raise TypeError(f"n must be an integer, got {type(n).__name__}")
    if n < 2:
        raise ValueError(f"need at least 2 perturbation points, got {n}")

    indices = np.arange(n)
    diff = indices[None, :] - indices[:, None]

    matrix = np.zeros((n, n))
    off_diagonal = diff != 0
    matrix[off_diagonal] = 1.0 / (np.pi * diff[off_diagonal])

    return matrix

synchronous(series)

Synchronous 2DCOS matrix Φ of a spectral series.

Φ(ν₁, ν₂) = Y^T @ Y / (m - 1) over the mean-centred intensities, which is just the sample covariance across the perturbation. Large positive values mean two wavenumbers rise and fall together; negative means they move in opposition. The diagonal is the autopower spectrum — where the perturbation does the most work.

Parameters:

Name	Type	Description	Default
series		A SpectralSeries with at least two perturbation points.	required

Returns:

Type	Description
ndarray	The (n, n) synchronous matrix, symmetric by construction.

Source code in src/crosspeak/noda.py

def synchronous(series):
    """Synchronous 2DCOS matrix Φ of a spectral series.

    Φ(ν₁, ν₂) = Y^T @ Y / (m - 1) over the mean-centred intensities, which is
    just the sample covariance across the perturbation. Large positive values
    mean two wavenumbers rise and fall together; negative means they move in
    opposition. The diagonal is the autopower spectrum — where the perturbation
    does the most work.

    Parameters
    ----------
    series
        A `SpectralSeries` with at least two perturbation points.

    Returns
    -------
    np.ndarray
        The (n, n) synchronous matrix, symmetric by construction.
    """
    if not isinstance(series, SpectralSeries):
        raise TypeError(f"expected SpectralSeries, got {type(series).__name__}")

    m = series.n_perturbations
    Y = series.intensities - series.intensities.mean(axis=0)
    return (Y.T @ Y) / (m - 1)

asynchronous(series)

Asynchronous 2DCOS matrix Ψ of a spectral series.

Ψ(ν₁, ν₂) = Y^T @ N @ Y / (m - 1), with N the Hilbert-Noda matrix. Where Φ says which bands move together, Ψ says in what order: read against the sign of Φ, a nonzero Ψ resolves which of two overlapping features responds first. Antisymmetric, with a zero diagonal.

Same mean-centring and conventions as synchronous; see hilbert_noda_matrix for the sign that fixes the direction of the ordering.

Parameters:

Name	Type	Description	Default
series		A SpectralSeries with at least two perturbation points.	required

Returns:

Type	Description
ndarray	The (n, n) asynchronous matrix.

Source code in src/crosspeak/noda.py

def asynchronous(series):
    """Asynchronous 2DCOS matrix Ψ of a spectral series.

    Ψ(ν₁, ν₂) = Y^T @ N @ Y / (m - 1), with N the Hilbert-Noda matrix. Where Φ
    says which bands move together, Ψ says in what order: read against the sign
    of Φ, a nonzero Ψ resolves which of two overlapping features responds
    first. Antisymmetric, with a zero diagonal.

    Same mean-centring and conventions as `synchronous`; see
    `hilbert_noda_matrix` for the sign that fixes the direction of the ordering.

    Parameters
    ----------
    series
        A `SpectralSeries` with at least two perturbation points.

    Returns
    -------
    np.ndarray
        The (n, n) asynchronous matrix.
    """
    if not isinstance(series, SpectralSeries):
        raise TypeError(f"expected SpectralSeries, got {type(series).__name__}")

    m = series.n_perturbations
    Y = series.intensities - series.intensities.mean(axis=0)
    N = hilbert_noda_matrix(m)
    return (Y.T @ N @ Y) / (m - 1)

moving_window(series, *, window_size=5)

Moving-window 2DCOS: correlation over a sliding perturbation window.

Full-range 2DCOS gives one pair of maps summarising the whole series, which hides when things happened — a system that reorganises in two stages produces a single map with both stages blended into it. Moving-window analysis slides a short window along the perturbation axis and correlates within each position, so features can be seen switching on and off as the perturbation advances.

Each window is mean-centred against its own local mean rather than the mean of the whole series. That is what makes the result local.

The window is constant in index, not in perturbation value. For an unevenly sampled series — water contents of 0, 1, 2, 5, 10, 15, 20, 30, say — a five-point window spans 0 to 10 at one end and 10 to 30 at the other. This is the standard formulation, but it means centres should be read alongside the maps rather than assumed evenly spaced.

Short windows localise sharply but leave few points to correlate, and the asynchronous map degrades first. Windows of 5 to 11 points are usual.

Parameters:

Name	Type	Description	Default
series	SpectralSeries	A SpectralSeries. Crop to the region of interest first: the result holds two (n, n) matrices per window position, so an uncropped spectrum can run to tens of megabytes per window.	required
window_size	int	Perturbation points per window. Must be odd, at least 3, and no larger than the series. Default 5. Noda writes this as 2m+1, so 5 is m=2.	5

Returns:

Type	Description
MovingWindowResult	Named tuple of (sync, asyn, centres), with n_windows = n_perturbations - window_size + 1.

Raises:

Type	Description
ValueError	If window_size is even, below 3, or larger than the series.

Source code in src/crosspeak/noda.py

def moving_window(
    series: SpectralSeries,
    *,
    window_size: int = 5,
) -> MovingWindowResult:
    """Moving-window 2DCOS: correlation over a sliding perturbation window.

    Full-range 2DCOS gives one pair of maps summarising the whole series, which
    hides *when* things happened — a system that reorganises in two stages
    produces a single map with both stages blended into it. Moving-window
    analysis slides a short window along the perturbation axis and correlates
    within each position, so features can be seen switching on and off as the
    perturbation advances.

    Each window is mean-centred against its own local mean rather than the mean
    of the whole series. That is what makes the result local.

    The window is constant in *index*, not in perturbation value. For an
    unevenly sampled series — water contents of 0, 1, 2, 5, 10, 15, 20, 30, say
    — a five-point window spans 0 to 10 at one end and 10 to 30 at the other.
    This is the standard formulation, but it means `centres` should be read
    alongside the maps rather than assumed evenly spaced.

    Short windows localise sharply but leave few points to correlate, and the
    asynchronous map degrades first. Windows of 5 to 11 points are usual.

    Parameters
    ----------
    series
        A `SpectralSeries`. Crop to the region of interest first: the result
        holds two (n, n) matrices per window position, so an uncropped spectrum
        can run to tens of megabytes per window.
    window_size
        Perturbation points per window. Must be odd, at least 3, and no larger
        than the series. Default 5. Noda writes this as 2m+1, so 5 is m=2.

    Returns
    -------
    MovingWindowResult
        Named tuple of (sync, asyn, centres), with
        `n_windows = n_perturbations - window_size + 1`.

    Raises
    ------
    ValueError
        If `window_size` is even, below 3, or larger than the series.
    """
    if not isinstance(series, SpectralSeries):
        raise TypeError(f"expected SpectralSeries, got {type(series).__name__}")
    if not isinstance(window_size, (int, np.integer)):
        raise TypeError(f"window_size must be an integer, got {type(window_size).__name__}")
    if window_size % 2 == 0:
        raise ValueError(f"window_size must be odd, got {window_size}")
    if window_size < 3:
        raise ValueError(f"window_size must be at least 3, got {window_size}")

    m = series.n_perturbations
    if window_size > m:
        raise ValueError(
            f"window_size ({window_size}) exceeds the number of perturbation points ({m})"
        )

    n_windows = m - window_size + 1
    n_wavenumbers = series.n_wavenumbers
    half = window_size // 2

    sync = np.empty((n_windows, n_wavenumbers, n_wavenumbers))
    asyn = np.empty((n_windows, n_wavenumbers, n_wavenumbers))

    for i in range(n_windows):
        window = SpectralSeries(
            wavenumbers=series.wavenumbers,
            perturbations=series.perturbations[i : i + window_size],
            intensities=series.intensities[i : i + window_size],
            name=series.name,
        )
        sync[i] = synchronous(window)
        asyn[i] = asynchronous(window)

    return MovingWindowResult(
        sync=sync,
        asyn=asyn,
        centres=series.perturbations[half : half + n_windows],
    )

find_autopeaks(phi, wavenumbers, *, prominence_frac=0.02, **kwargs)

Detect autopeaks on the diagonal of a synchronous 2DCOS matrix.

The diagonal of Phi is the autopower spectrum — the variance of each wavenumber across the perturbation series. Local maxima of this spectrum are Noda's autopeaks and identify the wavenumbers that change most strongly under the perturbation.

Parameters:

Name	Type	Description	Default
phi	ndarray	Synchronous 2DCOS matrix, shape (n, n), as returned by synchronous().	required
wavenumbers	ndarray	Wavenumber axis, shape (n,), strictly monotonic.	required
prominence_frac	float	Minimum peak prominence, as a fraction of the diagonal maximum. Default 0.02. Set to 0 to disable the prominence filter.	0.02
**kwargs		Additional keyword arguments passed through to scipy.signal.find_peaks (e.g. height, width, distance). Passing prominence directly overrides the value derived from prominence_frac.	{}

Returns:

Type	Description
AutopeakResult	Named tuple of (positions, intensities). Empty arrays if none found.

Source code in src/crosspeak/noda.py

def find_autopeaks(
    phi: np.ndarray,
    wavenumbers: np.ndarray,
    *,
    prominence_frac: float = 0.02,
    **kwargs,
) -> AutopeakResult:
    """Detect autopeaks on the diagonal of a synchronous 2DCOS matrix.

    The diagonal of Phi is the autopower spectrum — the variance of each
    wavenumber across the perturbation series. Local maxima of this spectrum
    are Noda's autopeaks and identify the wavenumbers that change most
    strongly under the perturbation.

    Parameters
    ----------
    phi
        Synchronous 2DCOS matrix, shape (n, n), as returned by `synchronous()`.
    wavenumbers
        Wavenumber axis, shape (n,), strictly monotonic.
    prominence_frac
        Minimum peak prominence, as a fraction of the diagonal maximum.
        Default 0.02. Set to 0 to disable the prominence filter.
    **kwargs
        Additional keyword arguments passed through to `scipy.signal.find_peaks`
        (e.g. `height`, `width`, `distance`). Passing `prominence` directly
        overrides the value derived from `prominence_frac`.

    Returns
    -------
    AutopeakResult
        Named tuple of (positions, intensities). Empty arrays if none found.
    """
    phi = np.asarray(phi)
    wavenumbers = np.asarray(wavenumbers)

    if phi.ndim != 2:
        raise ValueError(f"phi must be 2D, got shape {phi.shape}")
    if phi.shape[0] != phi.shape[1]:
        raise ValueError(f"phi must be square, got shape {phi.shape}")
    if wavenumbers.ndim != 1:
        raise ValueError(f"wavenumbers must be 1D, got shape {wavenumbers.shape}")
    if wavenumbers.size != phi.shape[0]:
        raise ValueError(
            f"wavenumbers size ({wavenumbers.size}) does not match phi shape ({phi.shape})"
        )

    diagonal = np.diag(phi)

    # User-supplied absolute prominence wins over prominence_frac
    kwargs.setdefault("prominence", prominence_frac * diagonal.max())

    indices, _ = find_peaks(diagonal, **kwargs)

    return AutopeakResult(
        positions=wavenumbers[indices],
        intensities=diagonal[indices],
    )

synchronous_hetero(series_1, series_2)

Heterospectral synchronous 2DCOS matrix.

Computes Φ(ν₁, ν₂) = Y_1^T @ Y_2 / (m - 1), where Y_1 and Y_2 are the mean-centred intensities of the two series. The two series must share the same perturbation array (rows are paired by index).

Parameters:

Name	Type	Description	Default
series_1	SpectralSeries	SpectralSeries supplying ν₁ — maps to matrix axis 0 (rows), and to the y-axis when the result is plotted.	required
series_2	SpectralSeries	SpectralSeries supplying ν₂ — maps to matrix axis 1 (cols), and to the x-axis when the result is plotted.	required

Returns:

Type	Description
ndarray	Shape (n_1, n_2). Not symmetric in general.

Raises:

Type	Description
ValueError	If the two series have different perturbation arrays.

Source code in src/crosspeak/noda.py

def synchronous_hetero(
    series_1: SpectralSeries,
    series_2: SpectralSeries,
) -> np.ndarray:
    """Heterospectral synchronous 2DCOS matrix.

    Computes Φ(ν₁, ν₂) = Y_1^T @ Y_2 / (m - 1), where Y_1 and Y_2 are the
    mean-centred intensities of the two series. The two series must share
    the same perturbation array (rows are paired by index).

    Parameters
    ----------
    series_1
        SpectralSeries supplying ν₁ — maps to matrix axis 0 (rows), and to
        the y-axis when the result is plotted.
    series_2
        SpectralSeries supplying ν₂ — maps to matrix axis 1 (cols), and to
        the x-axis when the result is plotted.

    Returns
    -------
    np.ndarray
        Shape `(n_1, n_2)`. Not symmetric in general.

    Raises
    ------
    ValueError
        If the two series have different perturbation arrays.
    """
    _check_paired_perturbations(series_1, series_2)
    m = series_1.n_perturbations
    Y_1 = series_1.intensities - series_1.intensities.mean(axis=0)
    Y_2 = series_2.intensities - series_2.intensities.mean(axis=0)
    return (Y_1.T @ Y_2) / (m - 1)

asynchronous_hetero(series_1, series_2)

Heterospectral asynchronous 2DCOS matrix.

Computes Ψ(ν₁, ν₂) = Y_1^T @ N @ Y_2 / (m - 1), where N is the Hilbert-Noda matrix of shape (m, m).

See synchronous_hetero for axis conventions.

Raises:

Type	Description
ValueError	If the two series have different perturbation arrays.

Source code in src/crosspeak/noda.py

def asynchronous_hetero(
    series_1: SpectralSeries,
    series_2: SpectralSeries,
) -> np.ndarray:
    """Heterospectral asynchronous 2DCOS matrix.

    Computes Ψ(ν₁, ν₂) = Y_1^T @ N @ Y_2 / (m - 1), where N is the
    Hilbert-Noda matrix of shape (m, m).

    See `synchronous_hetero` for axis conventions.

    Raises
    ------
    ValueError
        If the two series have different perturbation arrays.
    """
    _check_paired_perturbations(series_1, series_2)
    m = series_1.n_perturbations
    Y_1 = series_1.intensities - series_1.intensities.mean(axis=0)
    Y_2 = series_2.intensities - series_2.intensities.mean(axis=0)
    N = hilbert_noda_matrix(m)
    return (Y_1.T @ N @ Y_2) / (m - 1)