decompose.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Spectrogram decomposition
=========================
.. autosummary::
    :toctree: generated/

    decompose
    hpss
    nn_filter
"""

import numpy as np

import scipy.sparse
from scipy.ndimage import median_filter

import sklearn.decomposition

from . import core
from ._cache import cache
from . import segment
from . import util
from .util.exceptions import ParameterError
from typing import Any, Callable, List, Optional, Tuple, Union
from ._typing import _IntLike_co, _FloatLike_co

__all__ = ["decompose", "hpss", "nn_filter"]


def decompose(
    S: np.ndarray,
    *,
    n_components: Optional[int] = None,
    transformer: Optional[object] = None,
    sort: bool = False,
    fit: bool = True,
    **kwargs: Any,
) -> Tuple[np.ndarray, np.ndarray]:
    """Decompose a feature matrix.

    Given a spectrogram ``S``, produce a decomposition into ``components``
    and ``activations`` such that ``S ~= components.dot(activations)``.

    By default, this is done with with non-negative matrix factorization (NMF),
    but any `sklearn.decomposition`-type object will work.

    Parameters
    ----------
    S : np.ndarray [shape=(..., n_features, n_samples), dtype=float]
        The input feature matrix (e.g., magnitude spectrogram)

        If the input has multiple channels (leading dimensions), they will be automatically
        flattened prior to decomposition.

        If the input is multi-channel, channels and features are automatically flattened into
        a single axis before the decomposition.
        For example, a stereo input `S` with shape `(2, n_features, n_samples)` is
        automatically reshaped to `(2 * n_features, n_samples)`.

    n_components : int > 0 [scalar] or None
        number of desired components

        if None, then ``n_features`` components are used

    transformer : None or object
        If None, use `sklearn.decomposition.NMF`

        Otherwise, any object with a similar interface to NMF should work.
        ``transformer`` must follow the scikit-learn convention, where
        input data is ``(n_samples, n_features)``.

        `transformer.fit_transform()` will be run on ``S.T`` (not ``S``),
        the return value of which is stored (transposed) as ``activations``

        The components will be retrieved as ``transformer.components_.T``::

            S ~= np.dot(activations, transformer.components_).T

        or equivalently::

            S ~= np.dot(transformer.components_.T, activations.T)

    sort : bool
        If ``True``, components are sorted by ascending peak frequency.

        .. note:: If used with ``transformer``, sorting is applied to copies
            of the decomposition parameters, and not to ``transformer``
            internal parameters.

        .. warning:: If the input array has more than two dimensions
            (e.g., if it's a multi-channel spectrogram), then axis sorting
            is not supported and a `ParameterError` exception is raised.

    fit : bool
        If `True`, components are estimated from the input ``S``.

        If `False`, components are assumed to be pre-computed and stored
        in ``transformer``, and are not changed.

    **kwargs : Additional keyword arguments to the default transformer
        `sklearn.decomposition.NMF`

    Returns
    -------
    components: np.ndarray [shape=(..., n_features, n_components)]
        matrix of components (basis elements).
    activations: np.ndarray [shape=(n_components, n_samples)]
        transformed matrix/activation matrix

    Raises
    ------
    ParameterError
        if ``fit`` is False and no ``transformer`` object is provided.

        if the input array is multi-channel and ``sort=True`` is specified.

    See Also
    --------
    sklearn.decomposition : SciKit-Learn matrix decomposition modules

    Examples
    --------
    Decompose a magnitude spectrogram into 16 components with NMF

    >>> y, sr = librosa.load(librosa.ex('pistachio'), duration=5)
    >>> S = np.abs(librosa.stft(y))
    >>> comps, acts = librosa.decompose.decompose(S, n_components=16)

    Sort components by ascending peak frequency

    >>> comps, acts = librosa.decompose.decompose(S, n_components=16,
    ...                                           sort=True)

    Or with sparse dictionary learning

    >>> import sklearn.decomposition
    >>> T = sklearn.decomposition.MiniBatchDictionaryLearning(n_components=16)
    >>> scomps, sacts = librosa.decompose.decompose(S, transformer=T, sort=True)

    >>> import matplotlib.pyplot as plt
    >>> layout = [list(".AAAA"), list("BCCCC"), list(".DDDD")]
    >>> fig, ax = plt.subplot_mosaic(layout, constrained_layout=True)
    >>> librosa.display.specshow(librosa.amplitude_to_db(S, ref=np.max),
    ...                          y_axis='log', x_axis='time', ax=ax['A'])
    >>> ax['A'].set(title='Input spectrogram')
    >>> ax['A'].label_outer()
    >>> librosa.display.specshow(librosa.amplitude_to_db(comps,
    >>>                                                  ref=np.max),
    >>>                          y_axis='log', ax=ax['B'])
    >>> ax['B'].set(title='Components')
    >>> ax['B'].label_outer()
    >>> ax['B'].sharey(ax['A'])
    >>> librosa.display.specshow(acts, x_axis='time', ax=ax['C'], cmap='gray_r')
    >>> ax['C'].set(ylabel='Components', title='Activations')
    >>> ax['C'].sharex(ax['A'])
    >>> ax['C'].label_outer()
    >>> S_approx = comps.dot(acts)
    >>> img = librosa.display.specshow(librosa.amplitude_to_db(S_approx,
    >>>                                                        ref=np.max),
    >>>                                y_axis='log', x_axis='time', ax=ax['D'])
    >>> ax['D'].set(title='Reconstructed spectrogram')
    >>> ax['D'].sharex(ax['A'])
    >>> ax['D'].sharey(ax['A'])
    >>> ax['D'].label_outer()
    >>> fig.colorbar(img, ax=list(ax.values()), format="%+2.f dB")
    """
    # Do a swapaxes and unroll
    orig_shape = list(S.shape)

    if S.ndim > 2 and sort:
        raise ParameterError(
            "Parameter sort=True is unsupported for input with more than two dimensions"
        )

    # Transpose S and unroll feature dimensions
    # Use order='F' here to preserve the temporal ordering
    S = S.T.reshape((S.shape[-1], -1), order="F")

    if n_components is None:
        n_components = S.shape[-1]

    if transformer is None:
        if fit is False:
            raise ParameterError("fit must be True if transformer is None")

        transformer = sklearn.decomposition.NMF(n_components=n_components, **kwargs)

    # Suppressing type errors here because we don't want to overly restrict
    # the transformer object type
    activations: np.ndarray
    if fit:
        activations = transformer.fit_transform(S).T  # type: ignore
    else:
        activations = transformer.transform(S).T  # type: ignore

    components: np.ndarray = transformer.components_  # type: ignore
    component_shape = orig_shape[:-1] + [-1]
    # use order='F' here to preserve component ordering
    components = components.reshape(component_shape[::-1], order="F").T

    if sort:
        components, idx = util.axis_sort(components, index=True)
        activations = activations[idx]

    return components, activations


@cache(level=30)
def hpss(
    S: np.ndarray,
    *,
    kernel_size: Union[
        _IntLike_co, Tuple[_IntLike_co, _IntLike_co], List[_IntLike_co]
    ] = 31,
    power: float = 2.0,
    mask: bool = False,
    margin: Union[
        _FloatLike_co, Tuple[_FloatLike_co, _FloatLike_co], List[_FloatLike_co]
    ] = 1.0,
) -> Tuple[np.ndarray, np.ndarray]:
    """Median-filtering harmonic percussive source separation (HPSS).

    If ``margin = 1.0``, decomposes an input spectrogram ``S = H + P``
    where ``H`` contains the harmonic components,
    and ``P`` contains the percussive components.

    If ``margin > 1.0``, decomposes an input spectrogram ``S = H + P + R``
    where ``R`` contains residual components not included in ``H`` or ``P``.

    This implementation is based upon the algorithm described by [#]_ and [#]_.

    .. [#] Fitzgerald, Derry.
        "Harmonic/percussive separation using median filtering."
        13th International Conference on Digital Audio Effects (DAFX10),
        Graz, Austria, 2010.

    .. [#] Driedger, Müller, Disch.
        "Extending harmonic-percussive separation of audio."
        15th International Society for Music Information Retrieval Conference (ISMIR 2014),
        Taipei, Taiwan, 2014.

    Parameters
    ----------
    S : np.ndarray [shape=(..., d, n)]
        input spectrogram. May be real (magnitude) or complex.
        Multi-channel is supported.

    kernel_size : int or tuple (kernel_harmonic, kernel_percussive)
        kernel size(s) for the median filters.

        - If scalar, the same size is used for both harmonic and percussive.
        - If tuple, the first value specifies the width of the
          harmonic filter, and the second value specifies the width
          of the percussive filter.

    power : float > 0 [scalar]
        Exponent for the Wiener filter when constructing soft mask matrices.

    mask : bool
        Return the masking matrices instead of components.

        Masking matrices contain non-negative real values that
        can be used to measure the assignment of energy from ``S``
        into harmonic or percussive components.

        Components can be recovered by multiplying ``S * mask_H``
        or ``S * mask_P``.

    margin : float or tuple (margin_harmonic, margin_percussive)
        margin size(s) for the masks (as described in [2]_)

        - If scalar, the same size is used for both harmonic and percussive.
        - If tuple, the first value specifies the margin of the
          harmonic mask, and the second value specifies the margin
          of the percussive mask.

    Returns
    -------
    harmonic : np.ndarray [shape=(..., d, n)]
        harmonic component (or mask)
    percussive : np.ndarray [shape=(..., d, n)]
        percussive component (or mask)

    See Also
    --------
    librosa.util.softmask

    Notes
    -----
    This function caches at level 30.

    Examples
    --------
    Separate into harmonic and percussive

    >>> y, sr = librosa.load(librosa.ex('choice'), duration=5)
    >>> D = librosa.stft(y)
    >>> H, P = librosa.decompose.hpss(D)

    >>> import matplotlib.pyplot as plt
    >>> fig, ax = plt.subplots(nrows=3, sharex=True, sharey=True)
    >>> img = librosa.display.specshow(librosa.amplitude_to_db(np.abs(D),
    ...                                                        ref=np.max),
    ...                          y_axis='log', x_axis='time', ax=ax[0])
    >>> ax[0].set(title='Full power spectrogram')
    >>> ax[0].label_outer()
    >>> librosa.display.specshow(librosa.amplitude_to_db(np.abs(H),
    ...                                                  ref=np.max(np.abs(D))),
    ...                          y_axis='log', x_axis='time', ax=ax[1])
    >>> ax[1].set(title='Harmonic power spectrogram')
    >>> ax[1].label_outer()
    >>> librosa.display.specshow(librosa.amplitude_to_db(np.abs(P),
    ...                                                  ref=np.max(np.abs(D))),
    ...                          y_axis='log', x_axis='time', ax=ax[2])
    >>> ax[2].set(title='Percussive power spectrogram')
    >>> fig.colorbar(img, ax=ax, format='%+2.0f dB')

    Or with a narrower horizontal filter

    >>> H, P = librosa.decompose.hpss(D, kernel_size=(13, 31))

    Just get harmonic/percussive masks, not the spectra

    >>> mask_H, mask_P = librosa.decompose.hpss(D, mask=True)
    >>> mask_H
    array([[1.853e-03, 1.701e-04, ..., 9.922e-01, 1.000e+00],
           [2.316e-03, 2.127e-04, ..., 9.989e-01, 1.000e+00],
           ...,
           [8.195e-05, 6.939e-05, ..., 3.105e-04, 4.231e-04],
           [3.159e-05, 4.156e-05, ..., 6.216e-04, 6.188e-04]],
          dtype=float32)
    >>> mask_P
    array([[9.981e-01, 9.998e-01, ..., 7.759e-03, 3.201e-05],
           [9.977e-01, 9.998e-01, ..., 1.122e-03, 4.451e-06],
           ...,
           [9.999e-01, 9.999e-01, ..., 9.997e-01, 9.996e-01],
           [1.000e+00, 1.000e+00, ..., 9.994e-01, 9.994e-01]],
          dtype=float32)

    Separate into harmonic/percussive/residual components by using a margin > 1.0

    >>> H, P = librosa.decompose.hpss(D, margin=3.0)
    >>> R = D - (H+P)
    >>> y_harm = librosa.istft(H)
    >>> y_perc = librosa.istft(P)
    >>> y_resi = librosa.istft(R)

    Get a more isolated percussive component by widening its margin

    >>> H, P = librosa.decompose.hpss(D, margin=(1.0,5.0))
    """
    phase: Union[float, np.ndarray]

    if np.iscomplexobj(S):
        S, phase = core.magphase(S)
    else:
        phase = 1

    if isinstance(kernel_size, tuple) or isinstance(kernel_size, list):
        win_harm = kernel_size[0]
        win_perc = kernel_size[1]
    else:
        win_harm = kernel_size
        win_perc = kernel_size

    if isinstance(margin, tuple) or isinstance(margin, list):
        margin_harm = margin[0]
        margin_perc = margin[1]
    else:
        margin_harm = margin
        margin_perc = margin

    # margin minimum is 1.0
    if margin_harm < 1 or margin_perc < 1:
        raise ParameterError(
            "Margins must be >= 1.0. " "A typical range is between 1 and 10."
        )

    # shape for kernels
    harm_shape: List[_IntLike_co] = [1] * S.ndim
    harm_shape[-1] = win_harm

    perc_shape: List[_IntLike_co] = [1] * S.ndim
    perc_shape[-2] = win_perc

    # Compute median filters. Pre-allocation here preserves memory layout.
    harm = np.empty_like(S)
    harm[:] = median_filter(S, size=harm_shape, mode="reflect")

    perc = np.empty_like(S)
    perc[:] = median_filter(S, size=perc_shape, mode="reflect")

    split_zeros = margin_harm == 1 and margin_perc == 1

    mask_harm = util.softmask(
        harm, perc * margin_harm, power=power, split_zeros=split_zeros
    )

    mask_perc = util.softmask(
        perc, harm * margin_perc, power=power, split_zeros=split_zeros
    )

    if mask:
        return mask_harm, mask_perc

    return ((S * mask_harm) * phase, (S * mask_perc) * phase)


@cache(level=30)
def nn_filter(
    S: np.ndarray,
    *,
    rec: Optional[Union[scipy.sparse.spmatrix, np.ndarray]] = None,
    aggregate: Optional[Callable] = None,
    axis: int = -1,
    **kwargs: Any,
) -> np.ndarray:
    """Filter by nearest-neighbor aggregation.

    Each data point (e.g, spectrogram column) is replaced
    by aggregating its nearest neighbors in feature space.

    This can be useful for de-noising a spectrogram or feature matrix.

    The non-local means method [#]_ can be recovered by providing a
    weighted recurrence matrix as input and specifying ``aggregate=np.average``.

    Similarly, setting ``aggregate=np.median`` produces sparse de-noising
    as in REPET-SIM [#]_.

    .. [#] Buades, A., Coll, B., & Morel, J. M.
        (2005, June). A non-local algorithm for image denoising.
        In Computer Vision and Pattern Recognition, 2005.
        CVPR 2005. IEEE Computer Society Conference on (Vol. 2, pp. 60-65). IEEE.

    .. [#] Rafii, Z., & Pardo, B.
        (2012, October).  "Music/Voice Separation Using the Similarity Matrix."
        International Society for Music Information Retrieval Conference, 2012.

    Parameters
    ----------
    S : np.ndarray
        The input data (spectrogram) to filter. Multi-channel is supported.

    rec : (optional) scipy.sparse.spmatrix or np.ndarray
        Optionally, a pre-computed nearest-neighbor matrix
        as provided by `librosa.segment.recurrence_matrix`

    aggregate : function
        aggregation function (default: `np.mean`)

        If ``aggregate=np.average``, then a weighted average is
        computed according to the (per-row) weights in ``rec``.

        For all other aggregation functions, all neighbors
        are treated equally.

    axis : int
        The axis along which to filter (by default, columns)

    **kwargs
        Additional keyword arguments provided to
        `librosa.segment.recurrence_matrix` if ``rec`` is not provided

    Returns
    -------
    S_filtered : np.ndarray
        The filtered data, with shape equivalent to the input ``S``.

    Raises
    ------
    ParameterError
        if ``rec`` is provided and its shape is incompatible with ``S``.

    See Also
    --------
    decompose
    hpss
    librosa.segment.recurrence_matrix

    Notes
    -----
    This function caches at level 30.

    Examples
    --------
    De-noise a chromagram by non-local median filtering.
    By default this would use euclidean distance to select neighbors,
    but this can be overridden directly by setting the ``metric`` parameter.

    >>> y, sr = librosa.load(librosa.ex('brahms'),
    ...                      offset=30, duration=10)
    >>> chroma = librosa.feature.chroma_cqt(y=y, sr=sr)
    >>> chroma_med = librosa.decompose.nn_filter(chroma,
    ...                                          aggregate=np.median,
    ...                                          metric='cosine')

    To use non-local means, provide an affinity matrix and ``aggregate=np.average``.

    >>> rec = librosa.segment.recurrence_matrix(chroma, mode='affinity',
    ...                                         metric='cosine', sparse=True)
    >>> chroma_nlm = librosa.decompose.nn_filter(chroma, rec=rec,
    ...                                          aggregate=np.average)

    >>> import matplotlib.pyplot as plt
    >>> fig, ax = plt.subplots(nrows=5, sharex=True, sharey=True, figsize=(10, 10))
    >>> librosa.display.specshow(chroma, y_axis='chroma', x_axis='time', ax=ax[0])
    >>> ax[0].set(title='Unfiltered')
    >>> ax[0].label_outer()
    >>> librosa.display.specshow(chroma_med, y_axis='chroma', x_axis='time', ax=ax[1])
    >>> ax[1].set(title='Median-filtered')
    >>> ax[1].label_outer()
    >>> imgc = librosa.display.specshow(chroma_nlm, y_axis='chroma', x_axis='time', ax=ax[2])
    >>> ax[2].set(title='Non-local means')
    >>> ax[2].label_outer()
    >>> imgr1 = librosa.display.specshow(chroma - chroma_med,
    ...                          y_axis='chroma', x_axis='time', ax=ax[3])
    >>> ax[3].set(title='Original - median')
    >>> ax[3].label_outer()
    >>> imgr2 = librosa.display.specshow(chroma - chroma_nlm,
    ...                          y_axis='chroma', x_axis='time', ax=ax[4])
    >>> ax[4].label_outer()
    >>> ax[4].set(title='Original - NLM')
    >>> fig.colorbar(imgc, ax=ax[:3])
    >>> fig.colorbar(imgr1, ax=[ax[3]])
    >>> fig.colorbar(imgr2, ax=[ax[4]])
    """
    if aggregate is None:
        aggregate = np.mean

    rec_s: scipy.sparse.spmatrix

    if rec is None:
        kwargs = dict(kwargs)
        kwargs["sparse"] = True
        rec_s = segment.recurrence_matrix(S, axis=axis, **kwargs)
    elif not scipy.sparse.issparse(rec):
        rec_s = scipy.sparse.csc_matrix(rec)
    else:
        rec_s = rec

    if rec_s.shape[0] != S.shape[axis] or rec_s.shape[0] != rec_s.shape[1]:
        raise ParameterError(
            "Invalid self-similarity matrix shape "
            f"rec.shape={rec_s.shape} for S.shape={S.shape}"
        )

    return __nn_filter_helper(
        rec_s.data, rec_s.indices, rec_s.indptr, S.swapaxes(0, axis), aggregate
    ).swapaxes(0, axis)


def __nn_filter_helper(
    R_data, R_indices, R_ptr, S: np.ndarray, aggregate: Callable
) -> np.ndarray:
    """Nearest-neighbor filter helper function.

    This is an internal function, not for use outside of the decompose module.

    It applies the nearest-neighbor filter to S, assuming that the first index
    corresponds to observations.

    Parameters
    ----------
    R_data, R_indices, R_ptr : np.ndarrays
        The ``data``, ``indices``, and ``indptr`` of a scipy.sparse matrix
    S : np.ndarray
        The observation data to filter
    aggregate : callable
        The aggregation operator

    Returns
    -------
    S_out : np.ndarray like S
        The filtered data array
    """
    s_out = np.empty_like(S)

    for i in range(len(R_ptr) - 1):
        # Get the non-zeros out of the recurrence matrix
        targets = R_indices[R_ptr[i] : R_ptr[i + 1]]

        if not len(targets):
            s_out[i] = S[i]
            continue

        neighbors = np.take(S, targets, axis=0)

        if aggregate is np.average:
            weights = R_data[R_ptr[i] : R_ptr[i + 1]]
            s_out[i] = aggregate(neighbors, axis=0, weights=weights)
        else:
            s_out[i] = aggregate(neighbors, axis=0)

    return s_out