#!/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"]
[docs]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
[docs]@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)
[docs]@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