Orthogonal Subspace Projection¶

def osp(
    m: ArrayLike,
    resp: ArrayLike,
    resp_pxx: ArrayLike,
    f: ArrayLike,
    fs: float,
    min_resp_frequency: float = 0.1,
) -> OspResult:
    """Separate respiration-related and unrelated HRV modulation.

    Parameters
    ----------
    m : array_like
        Uniformly sampled, dimensionless HRV modulating signal.
    resp : array_like
        Respiration samples aligned with ``m`` on the same time grid.
    resp_pxx : array_like
        Finite, nonnegative respiration power spectral density.
    f : array_like
        Finite, strictly increasing frequency samples in hertz.
    fs : float
        Positive sampling frequency in hertz.
    min_resp_frequency : float, default=0.1
        Positive lower bound for the selected respiratory frequency.

    Returns
    -------
    OspResult
        Respiration-related modulation, orthogonal residual, and adaptive
        delayed-respiration model order. The component arrays align with
        ``m[delay - 1:]``.

    Raises
    ------
    TypeError
        If an input has an invalid numeric type.
    ValueError
        If vector shapes, spectrum values, frequencies, sampling parameters,
        or finite signal lengths violate the Biosiglib contract.

    Notes
    -----
    The Gram-matrix pseudoinverse uses the explicit Biosiglib binary64
    threshold, rather than NumPy's default pseudoinverse tolerance.

    Examples
    --------
    >>> result = osp(
    ...     [99, 1, 2, 3, 4, 5],
    ...     [1, 0, -1, 0, 1, 0],
    ...     [0, 0, 1],
    ...     [0, 0.5, 1],
    ...     1,
    ... )
    >>> result.delay
    2
    >>> np.allclose(result.m_resp + result.m_unrelated, [1, 2, 3, 4, 5])
    True
    """

    modulation = as_real_vector(m, name="m")
    respiration = as_real_vector(resp, name="resp")
    spectrum = as_real_vector(resp_pxx, name="resp_pxx")
    frequencies = as_real_vector(f, name="f")
    sampling_frequency = as_positive_real_scalar(fs, name="fs")
    minimum_frequency = as_positive_real_scalar(
        min_resp_frequency, name="min_resp_frequency"
    )

    _validate_spectrum(spectrum, frequencies)

    if modulation.size == 0 or respiration.size == 0:
        empty = np.asarray([], dtype=np.float64)
        return OspResult(empty.copy(), empty.copy(), empty.copy())
    if np.any(np.isnan(modulation)) or np.any(np.isnan(respiration)):
        empty = np.asarray([], dtype=np.float64)
        return OspResult(empty.copy(), empty.copy(), empty.copy())
    if np.any(np.isinf(modulation)) or np.any(np.isinf(respiration)):
        raise ValueError("m and resp must not contain infinite values")
    if modulation.size != respiration.size:
        raise ValueError("m and resp must have the same length")

    low_frequency, high_frequency = _occupied_power_limits(
        spectrum, frequencies
    )
    band_mask = (frequencies >= low_frequency) & (
        frequencies <= high_frequency
    )
    if not np.any(band_mask):
        band_mask = np.ones(frequencies.shape, dtype=bool)

    band_spectrum = spectrum[band_mask]
    band_frequencies = frequencies[band_mask]
    dominant_frequency = _dominant_frequency(
        band_spectrum, band_frequencies
    )
    dominant_frequency = max(dominant_frequency, minimum_frequency)
    delay = max(
        int(np.floor(2.0 * sampling_frequency / dominant_frequency + 0.5)),
        1,
    )

    if modulation.size < delay:
        return OspResult(np.nan, np.nan, delay)

    subspace = np.lib.stride_tricks.sliding_window_view(
        respiration, delay
    )
    gram = subspace.T @ subspace
    gram_inverse = _gram_pseudoinverse(gram)
    projection = subspace @ gram_inverse @ subspace.T
    delayed_modulation = modulation[delay - 1 :]
    m_resp = projection @ delayed_modulation
    m_unrelated = delayed_modulation - m_resp
    return OspResult(
        m_resp=np.asarray(m_resp, dtype=np.float64),
        m_unrelated=np.asarray(m_unrelated, dtype=np.float64),
        delay=delay,
    )