Upload 3363 files

4cef980 verified about 2 years ago

19.7 kB

#!/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 .util.decorators import deprecate_positional_args

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


@deprecate_positional_args
def decompose(
    S, *, n_components=None, transformer=None, sort=False, fit=True, **kwargs
):
    """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)

    if fit:
        activations = transformer.fit_transform(S).T
    else:
        activations = transformer.transform(S).T

    components = transformer.components_
    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)
@deprecate_positional_args
def hpss(S, *, kernel_size=31, power=2.0, mask=False, margin=1.0):
    """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))
    """

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

    if np.isscalar(kernel_size):
        win_harm = kernel_size
        win_perc = kernel_size
    else:
        win_harm = kernel_size[0]
        win_perc = kernel_size[1]

    if np.isscalar(margin):
        margin_harm = margin
        margin_perc = margin
    else:
        margin_harm = margin[0]
        margin_perc = margin[1]

    # 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 = [1 for _ in S.shape]
    harm_shape[-1] = win_harm

    perc_shape = [1 for _ in S.shape]
    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)
@deprecate_positional_args
def nn_filter(S, *, rec=None, aggregate=None, axis=-1, **kwargs):
    """Filtering by nearest-neighbors.
    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

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

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

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


def __nn_filter_helper(R_data, R_indices, R_ptr, S, aggregate):
    """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

	#!/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 .util.decorators import deprecate_positional_args

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


	@deprecate_positional_args
	def decompose(
	S, , n_components=None, transformer=None, sort=False, fit=True, *kwargs
	):
	"""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)

	if fit:
	activations = transformer.fit_transform(S).T
	else:
	activations = transformer.transform(S).T

	components = transformer.components_
	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)
	@deprecate_positional_args
	def hpss(S, *, kernel_size=31, power=2.0, mask=False, margin=1.0):
	"""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))

	"""

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

	if np.isscalar(kernel_size):
	win_harm = kernel_size
	win_perc = kernel_size
	else:
	win_harm = kernel_size[0]
	win_perc = kernel_size[1]

	if np.isscalar(margin):
	margin_harm = margin
	margin_perc = margin
	else:
	margin_harm = margin[0]
	margin_perc = margin[1]

	# 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 = [1 for _ in S.shape]
	harm_shape[-1] = win_harm

	perc_shape = [1 for _ in S.shape]
	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)
	@deprecate_positional_args
	def nn_filter(S, , rec=None, aggregate=None, axis=-1, *kwargs):
	"""Filtering by nearest-neighbors.

	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

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

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

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


	def __nn_filter_helper(R_data, R_indices, R_ptr, S, aggregate):
	"""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