#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""Utilities for spectral processing"""
from __future__ import annotations
import warnings
import numpy as np
import scipy
import scipy.ndimage
import scipy.signal
import scipy.interpolate
from numba import jit
from . import convert
from .fft import get_fftlib
from .audio import resample
from .._cache import cache
from .. import util
from ..util.exceptions import ParameterError
from ..filters import get_window, semitone_filterbank
from ..filters import window_sumsquare
from numpy.typing import DTypeLike
from typing import Any, Callable, Optional, Tuple, List, Union, overload
from typing_extensions import Literal
from .._typing import (
_WindowSpec,
_PadMode,
_PadModeSTFT,
_SequenceLike,
_ScalarOrSequence,
_ComplexLike_co,
_FloatLike_co
)
__all__ = [
"stft",
"istft",
"magphase",
"iirt",
"reassigned_spectrogram",
"phase_vocoder",
"perceptual_weighting",
"power_to_db",
"db_to_power",
"amplitude_to_db",
"db_to_amplitude",
"fmt",
"pcen",
"griffinlim",
]
[docs]@cache(level=20)
def stft(
y: np.ndarray,
*,
n_fft: int = 2048,
hop_length: Optional[int] = None,
win_length: Optional[int] = None,
window: _WindowSpec = "hann",
center: bool = True,
dtype: Optional[DTypeLike] = None,
pad_mode: _PadModeSTFT = "constant",
out: Optional[np.ndarray] = None,
) -> np.ndarray:
"""Short-time Fourier transform (STFT).
The STFT represents a signal in the time-frequency domain by
computing discrete Fourier transforms (DFT) over short overlapping
windows.
This function returns a complex-valued matrix D such that
- ``np.abs(D[..., f, t])`` is the magnitude of frequency bin ``f``
at frame ``t``, and
- ``np.angle(D[..., f, t])`` is the phase of frequency bin ``f``
at frame ``t``.
The integers ``t`` and ``f`` can be converted to physical units by means
of the utility functions `frames_to_samples` and `fft_frequencies`.
Parameters
----------
y : np.ndarray [shape=(..., n)], real-valued
input signal. Multi-channel is supported.
n_fft : int > 0 [scalar]
length of the windowed signal after padding with zeros.
The number of rows in the STFT matrix ``D`` is ``(1 + n_fft/2)``.
The default value, ``n_fft=2048`` samples, corresponds to a physical
duration of 93 milliseconds at a sample rate of 22050 Hz, i.e. the
default sample rate in librosa. This value is well adapted for music
signals. However, in speech processing, the recommended value is 512,
corresponding to 23 milliseconds at a sample rate of 22050 Hz.
In any case, we recommend setting ``n_fft`` to a power of two for
optimizing the speed of the fast Fourier transform (FFT) algorithm.
hop_length : int > 0 [scalar]
number of audio samples between adjacent STFT columns.
Smaller values increase the number of columns in ``D`` without
affecting the frequency resolution of the STFT.
If unspecified, defaults to ``win_length // 4`` (see below).
win_length : int <= n_fft [scalar]
Each frame of audio is windowed by ``window`` of length ``win_length``
and then padded with zeros to match ``n_fft``. Padding is added on
both the left- and the right-side of the window so that the window
is centered within the frame.
Smaller values improve the temporal resolution of the STFT (i.e. the
ability to discriminate impulses that are closely spaced in time)
at the expense of frequency resolution (i.e. the ability to discriminate
pure tones that are closely spaced in frequency). This effect is known
as the time-frequency localization trade-off and needs to be adjusted
according to the properties of the input signal ``y``.
If unspecified, defaults to ``win_length = n_fft``.
window : string, tuple, number, function, or np.ndarray [shape=(n_fft,)]
Either:
- a window specification (string, tuple, or number);
see `scipy.signal.get_window`
- a window function, such as `scipy.signal.windows.hann`
- a vector or array of length ``n_fft``
Defaults to a raised cosine window (`'hann'`), which is adequate for
most applications in audio signal processing.
.. see also:: `filters.get_window`
center : boolean
If ``True``, the signal ``y`` is padded so that frame
``D[:, t]`` is centered at ``y[t * hop_length]``.
If ``False``, then ``D[:, t]`` begins at ``y[t * hop_length]``.
Defaults to ``True``, which simplifies the alignment of ``D`` onto a
time grid by means of `librosa.frames_to_samples`.
Note, however, that ``center`` must be set to `False` when analyzing
signals with `librosa.stream`.
.. see also:: `librosa.stream`
dtype : np.dtype, optional
Complex numeric type for ``D``. Default is inferred to match the
precision of the input signal.
pad_mode : string or function
If ``center=True``, this argument is passed to `np.pad` for padding
the edges of the signal ``y``. By default (``pad_mode="constant"``),
``y`` is padded on both sides with zeros.
.. note:: Not all padding modes supported by `numpy.pad` are supported here.
`wrap`, `mean`, `maximum`, `median`, and `minimum` are not supported.
Other modes that depend at most on input values at the edges of the
signal (e.g., `constant`, `edge`, `linear_ramp`) are supported.
If ``center=False``, this argument is ignored.
.. see also:: `numpy.pad`
out : np.ndarray or None
A pre-allocated, complex-valued array to store the STFT results.
This must be of compatible shape and dtype for the given input parameters.
If `out` is larger than necessary for the provided input signal, then only
a prefix slice of `out` will be used.
If not provided, a new array is allocated and returned.
Returns
-------
D : np.ndarray [shape=(..., 1 + n_fft/2, n_frames), dtype=dtype]
Complex-valued matrix of short-term Fourier transform
coefficients.
If a pre-allocated `out` array is provided, then `D` will be
a reference to `out`.
If `out` is larger than necessary, then `D` will be a sliced
view: `D = out[..., :n_frames]`.
See Also
--------
istft : Inverse STFT
reassigned_spectrogram : Time-frequency reassigned spectrogram
Notes
-----
This function caches at level 20.
Examples
--------
>>> y, sr = librosa.load(librosa.ex('trumpet'))
>>> S = np.abs(librosa.stft(y))
>>> S
array([[5.395e-03, 3.332e-03, ..., 9.862e-07, 1.201e-05],
[3.244e-03, 2.690e-03, ..., 9.536e-07, 1.201e-05],
...,
[7.523e-05, 3.722e-05, ..., 1.188e-04, 1.031e-03],
[7.640e-05, 3.944e-05, ..., 5.180e-04, 1.346e-03]],
dtype=float32)
Use left-aligned frames, instead of centered frames
>>> S_left = librosa.stft(y, center=False)
Use a shorter hop length
>>> D_short = librosa.stft(y, hop_length=64)
Display a spectrogram
>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots()
>>> img = librosa.display.specshow(librosa.amplitude_to_db(S,
... ref=np.max),
... y_axis='log', x_axis='time', ax=ax)
>>> ax.set_title('Power spectrogram')
>>> fig.colorbar(img, ax=ax, format="%+2.0f dB")
"""
# By default, use the entire frame
if win_length is None:
win_length = n_fft
# Set the default hop, if it's not already specified
if hop_length is None:
hop_length = int(win_length // 4)
elif not util.is_positive_int(hop_length):
raise ParameterError(f"hop_length={hop_length} must be a positive integer")
# Check audio is valid
util.valid_audio(y, mono=False)
fft_window = get_window(window, win_length, fftbins=True)
# Pad the window out to n_fft size
fft_window = util.pad_center(fft_window, size=n_fft)
# Reshape so that the window can be broadcast
fft_window = util.expand_to(fft_window, ndim=1 + y.ndim, axes=-2)
# Pad the time series so that frames are centered
if center:
if pad_mode in ("wrap", "maximum", "mean", "median", "minimum"):
# Note: padding with a user-provided function "works", but
# use at your own risk.
# Since we don't pass-through kwargs here, any arguments
# to a user-provided pad function should be encapsulated
# by using functools.partial:
#
# >>> my_pad_func = functools.partial(pad_func, foo=x, bar=y)
# >>> librosa.stft(..., pad_mode=my_pad_func)
raise ParameterError(
f"pad_mode='{pad_mode}' is not supported by librosa.stft"
)
if n_fft > y.shape[-1]:
warnings.warn(
f"n_fft={n_fft} is too large for input signal of length={y.shape[-1]}"
)
# Set up the padding array to be empty, and we'll fix the target dimension later
padding = [(0, 0) for _ in range(y.ndim)]
# How many frames depend on left padding?
start_k = int(np.ceil(n_fft // 2 / hop_length))
# What's the first frame that depends on extra right-padding?
tail_k = (y.shape[-1] + n_fft // 2 - n_fft) // hop_length + 1
if tail_k <= start_k:
# If tail and head overlap, then just copy-pad the signal and carry on
start = 0
extra = 0
padding[-1] = (n_fft // 2, n_fft // 2)
y = np.pad(y, padding, mode=pad_mode)
else:
# If tail and head do not overlap, then we can implement padding on each part separately
# and avoid a full copy-pad
# "Middle" of the signal starts here, and does not depend on head padding
start = start_k * hop_length - n_fft // 2
padding[-1] = (n_fft // 2, 0)
# +1 here is to ensure enough samples to fill the window
# fixes bug #1567
y_pre = np.pad(
y[..., : (start_k - 1) * hop_length - n_fft // 2 + n_fft + 1],
padding,
mode=pad_mode,
)
y_frames_pre = util.frame(y_pre, frame_length=n_fft, hop_length=hop_length)
# Trim this down to the exact number of frames we should have
y_frames_pre = y_frames_pre[..., :start_k]
# How many extra frames do we have from the head?
extra = y_frames_pre.shape[-1]
# Determine if we have any frames that will fit inside the tail pad
if tail_k * hop_length - n_fft // 2 + n_fft <= y.shape[-1] + n_fft // 2:
padding[-1] = (0, n_fft // 2)
y_post = np.pad(
y[..., (tail_k) * hop_length - n_fft // 2 :], padding, mode=pad_mode
)
y_frames_post = util.frame(
y_post, frame_length=n_fft, hop_length=hop_length
)
# How many extra frames do we have from the tail?
extra += y_frames_post.shape[-1]
else:
# In this event, the first frame that touches tail padding would run off
# the end of the padded array
# We'll circumvent this by allocating an empty frame buffer for the tail
# this keeps the subsequent logic simple
post_shape = list(y_frames_pre.shape)
post_shape[-1] = 0
y_frames_post = np.empty_like(y_frames_pre, shape=post_shape)
else:
if n_fft > y.shape[-1]:
raise ParameterError(
f"n_fft={n_fft} is too large for uncentered analysis of input signal of length={y.shape[-1]}"
)
# "Middle" of the signal starts at sample 0
start = 0
# We have no extra frames
extra = 0
fft = get_fftlib()
if dtype is None:
dtype = util.dtype_r2c(y.dtype)
# Window the time series.
y_frames = util.frame(y[..., start:], frame_length=n_fft, hop_length=hop_length)
# Pre-allocate the STFT matrix
shape = list(y_frames.shape)
# This is our frequency dimension
shape[-2] = 1 + n_fft // 2
# If there's padding, there will be extra head and tail frames
shape[-1] += extra
if out is None:
stft_matrix = np.zeros(shape, dtype=dtype, order="F")
elif not (np.allclose(out.shape[:-1], shape[:-1]) and out.shape[-1] >= shape[-1]):
raise ParameterError(
f"Shape mismatch for provided output array out.shape={out.shape} and target shape={shape}"
)
elif not np.iscomplexobj(out):
raise ParameterError(f"output with dtype={out.dtype} is not of complex type")
else:
if np.allclose(shape, out.shape):
stft_matrix = out
else:
stft_matrix = out[..., : shape[-1]]
# Fill in the warm-up
if center and extra > 0:
off_start = y_frames_pre.shape[-1]
stft_matrix[..., :off_start] = fft.rfft(fft_window * y_frames_pre, axis=-2)
off_end = y_frames_post.shape[-1]
if off_end > 0:
stft_matrix[..., -off_end:] = fft.rfft(fft_window * y_frames_post, axis=-2)
else:
off_start = 0
n_columns = int(
util.MAX_MEM_BLOCK // (np.prod(y_frames.shape[:-1]) * y_frames.itemsize)
)
n_columns = max(n_columns, 1)
for bl_s in range(0, y_frames.shape[-1], n_columns):
bl_t = min(bl_s + n_columns, y_frames.shape[-1])
stft_matrix[..., bl_s + off_start : bl_t + off_start] = fft.rfft(
fft_window * y_frames[..., bl_s:bl_t], axis=-2
)
return stft_matrix
[docs]@cache(level=30)
def istft(
stft_matrix: np.ndarray,
*,
hop_length: Optional[int] = None,
win_length: Optional[int] = None,
n_fft: Optional[int] = None,
window: _WindowSpec = "hann",
center: bool = True,
dtype: Optional[DTypeLike] = None,
length: Optional[int] = None,
out: Optional[np.ndarray] = None,
) -> np.ndarray:
"""
Inverse short-time Fourier transform (ISTFT).
Converts a complex-valued spectrogram ``stft_matrix`` to time-series ``y``
by minimizing the mean squared error between ``stft_matrix`` and STFT of
``y`` as described in [#]_ up to Section 2 (reconstruction from MSTFT).
In general, window function, hop length and other parameters should be same
as in stft, which mostly leads to perfect reconstruction of a signal from
unmodified ``stft_matrix``.
.. [#] D. W. Griffin and J. S. Lim,
"Signal estimation from modified short-time Fourier transform,"
IEEE Trans. ASSP, vol.32, no.2, pp.236–243, Apr. 1984.
Parameters
----------
stft_matrix : np.ndarray [shape=(..., 1 + n_fft//2, t)]
STFT matrix from ``stft``
hop_length : int > 0 [scalar]
Number of frames between STFT columns.
If unspecified, defaults to ``win_length // 4``.
win_length : int <= n_fft = 2 * (stft_matrix.shape[0] - 1)
When reconstructing the time series, each frame is windowed
and each sample is normalized by the sum of squared window
according to the ``window`` function (see below).
If unspecified, defaults to ``n_fft``.
n_fft : int > 0 or None
The number of samples per frame in the input spectrogram.
By default, this will be inferred from the shape of ``stft_matrix``.
However, if an odd frame length was used, you can specify the correct
length by setting ``n_fft``.
window : string, tuple, number, function, np.ndarray [shape=(n_fft,)]
- a window specification (string, tuple, or number);
see `scipy.signal.get_window`
- a window function, such as `scipy.signal.windows.hann`
- a user-specified window vector of length ``n_fft``
.. see also:: `filters.get_window`
center : boolean
- If ``True``, ``D`` is assumed to have centered frames.
- If ``False``, ``D`` is assumed to have left-aligned frames.
dtype : numeric type
Real numeric type for ``y``. Default is to match the numerical
precision of the input spectrogram.
length : int > 0, optional
If provided, the output ``y`` is zero-padded or clipped to exactly
``length`` samples.
out : np.ndarray or None
A pre-allocated, complex-valued array to store the reconstructed signal
``y``. This must be of the correct shape for the given input parameters.
If not provided, a new array is allocated and returned.
Returns
-------
y : np.ndarray [shape=(..., n)]
time domain signal reconstructed from ``stft_matrix``.
If ``stft_matrix`` contains more than two axes
(e.g., from a stereo input signal), then ``y`` will match shape on the leading dimensions.
See Also
--------
stft : Short-time Fourier Transform
Notes
-----
This function caches at level 30.
Examples
--------
>>> y, sr = librosa.load(librosa.ex('trumpet'))
>>> D = librosa.stft(y)
>>> y_hat = librosa.istft(D)
>>> y_hat
array([-1.407e-03, -4.461e-04, ..., 5.131e-06, -1.417e-05],
dtype=float32)
Exactly preserving length of the input signal requires explicit padding.
Otherwise, a partial frame at the end of ``y`` will not be represented.
>>> n = len(y)
>>> n_fft = 2048
>>> y_pad = librosa.util.fix_length(y, size=n + n_fft // 2)
>>> D = librosa.stft(y_pad, n_fft=n_fft)
>>> y_out = librosa.istft(D, length=n)
>>> np.max(np.abs(y - y_out))
8.940697e-08
"""
if n_fft is None:
n_fft = 2 * (stft_matrix.shape[-2] - 1)
# By default, use the entire frame
if win_length is None:
win_length = n_fft
# Set the default hop, if it's not already specified
if hop_length is None:
hop_length = int(win_length // 4)
ifft_window = get_window(window, win_length, fftbins=True)
# Pad out to match n_fft, and add broadcasting axes
ifft_window = util.pad_center(ifft_window, size=n_fft)
ifft_window = util.expand_to(ifft_window, ndim=stft_matrix.ndim, axes=-2)
# For efficiency, trim STFT frames according to signal length if available
if length:
if center:
padded_length = length + 2 * (n_fft // 2)
else:
padded_length = length
n_frames = min(stft_matrix.shape[-1], int(np.ceil(padded_length / hop_length)))
else:
n_frames = stft_matrix.shape[-1]
if dtype is None:
dtype = util.dtype_c2r(stft_matrix.dtype)
shape = list(stft_matrix.shape[:-2])
expected_signal_len = n_fft + hop_length * (n_frames - 1)
if length:
expected_signal_len = length
elif center:
expected_signal_len -= 2 * (n_fft // 2)
shape.append(expected_signal_len)
if out is None:
y = np.zeros(shape, dtype=dtype)
elif not np.allclose(out.shape, shape):
raise ParameterError(
f"Shape mismatch for provided output array out.shape={out.shape} != {shape}"
)
else:
y = out
# Since we'll be doing overlap-add here, this needs to be initialized to zero.
y.fill(0.0)
fft = get_fftlib()
if center:
# First frame that does not depend on padding
# k * hop_length - n_fft//2 >= 0
# k * hop_length >= n_fft // 2
# k >= (n_fft//2 / hop_length)
start_frame = int(np.ceil((n_fft // 2) / hop_length))
# Do overlap-add on the head block
ytmp = ifft_window * fft.irfft(stft_matrix[..., :start_frame], n=n_fft, axis=-2)
shape[-1] = n_fft + hop_length * (start_frame - 1)
head_buffer = np.zeros(shape, dtype=dtype)
__overlap_add(head_buffer, ytmp, hop_length)
# If y is smaller than the head buffer, take everything
if y.shape[-1] < shape[-1] - n_fft // 2:
y[..., :] = head_buffer[..., n_fft // 2 : y.shape[-1] + n_fft // 2]
else:
# Trim off the first n_fft//2 samples from the head and copy into target buffer
y[..., : shape[-1] - n_fft // 2] = head_buffer[..., n_fft // 2 :]
# This offset compensates for any differences between frame alignment
# and padding truncation
offset = start_frame * hop_length - n_fft // 2
else:
start_frame = 0
offset = 0
n_columns = int(
util.MAX_MEM_BLOCK // (np.prod(stft_matrix.shape[:-1]) * stft_matrix.itemsize)
)
n_columns = max(n_columns, 1)
frame = 0
for bl_s in range(start_frame, n_frames, n_columns):
bl_t = min(bl_s + n_columns, n_frames)
# invert the block and apply the window function
ytmp = ifft_window * fft.irfft(stft_matrix[..., bl_s:bl_t], n=n_fft, axis=-2)
# Overlap-add the istft block starting at the i'th frame
__overlap_add(y[..., frame * hop_length + offset :], ytmp, hop_length)
frame += bl_t - bl_s
# Normalize by sum of squared window
ifft_window_sum = window_sumsquare(
window=window,
n_frames=n_frames,
win_length=win_length,
n_fft=n_fft,
hop_length=hop_length,
dtype=dtype,
)
if center:
start = n_fft // 2
else:
start = 0
ifft_window_sum = util.fix_length(ifft_window_sum[..., start:], size=y.shape[-1])
approx_nonzero_indices = ifft_window_sum > util.tiny(ifft_window_sum)
y[..., approx_nonzero_indices] /= ifft_window_sum[approx_nonzero_indices]
return y
@jit(nopython=True, cache=True)
def __overlap_add(y, ytmp, hop_length):
# numba-accelerated overlap add for inverse stft
# y is the pre-allocated output buffer
# ytmp is the windowed inverse-stft frames
# hop_length is the hop-length of the STFT analysis
n_fft = ytmp.shape[-2]
N = n_fft
for frame in range(ytmp.shape[-1]):
sample = frame * hop_length
if N > y.shape[-1] - sample:
N = y.shape[-1] - sample
y[..., sample : (sample + N)] += ytmp[..., :N, frame]
def __reassign_frequencies(
y: np.ndarray,
sr: float = 22050,
S: Optional[np.ndarray] = None,
n_fft: int = 2048,
hop_length: Optional[int] = None,
win_length: Optional[int] = None,
window: _WindowSpec = "hann",
center: bool = True,
dtype: Optional[DTypeLike] = None,
pad_mode: _PadModeSTFT = "constant",
) -> Tuple[np.ndarray, np.ndarray]:
"""Instantaneous frequencies based on a spectrogram representation.
The reassignment vector is calculated using equation 5.20 in Flandrin,
Auger, & Chassande-Mottin 2002::
omega_reassigned = omega - np.imag(S_dh/S_h)
where ``S_h`` is the complex STFT calculated using the original window, and
``S_dh`` is the complex STFT calculated using the derivative of the original
window.
See `reassigned_spectrogram` for references.
It is recommended to use ``pad_mode="wrap"`` or else ``center=False``, rather
than the defaults. Frequency reassignment assumes that the energy in each
FFT bin is associated with exactly one signal component. Reflection padding
at the edges of the signal may invalidate the reassigned estimates in the
boundary frames.
Parameters
----------
y : np.ndarray [shape=(..., n,)], real-valued
audio time series. Multi-channel is supported.
sr : number > 0 [scalar]
sampling rate of ``y``
S : np.ndarray [shape=(..., d, t)] or None
(optional) complex STFT calculated using the other arguments provided
to `__reassign_frequencies`
n_fft : int > 0 [scalar]
FFT window size. Defaults to 2048.
hop_length : int > 0 [scalar]
hop length, number samples between subsequent frames.
If not supplied, defaults to ``win_length // 4``.
win_length : int > 0, <= n_fft
Window length. Defaults to ``n_fft``.
See ``stft`` for details.
window : string, tuple, number, function, or np.ndarray [shape=(n_fft,)]
- a window specification (string, tuple, number);
see `scipy.signal.get_window`
- a window function, such as `scipy.signal.windows.hann`
- a user-specified window vector of length ``n_fft``
See `stft` for details.
.. see also:: `filters.get_window`
center : boolean
- If ``True``, the signal ``y`` is padded so that frame
``S[:, t]`` is centered at ``y[t * hop_length]``.
- If ``False``, then ``S[:, t]`` begins at ``y[t * hop_length]``.
dtype : numeric type
Complex numeric type for ``S``. Default is inferred to match
the numerical precision of the input signal.
pad_mode : string
If ``center=True``, the padding mode to use at the edges of the signal.
By default, STFT uses zero padding.
Returns
-------
freqs : np.ndarray [shape=(..., 1 + n_fft/2, t), dtype=real]
Instantaneous frequencies:
``freqs[f, t]`` is the frequency for bin ``f``, frame ``t``.
S : np.ndarray [shape=(..., 1 + n_fft/2, t), dtype=complex]
Short-time Fourier transform
Warns
-----
RuntimeWarning
Frequencies with zero support will produce a divide-by-zero warning and
will be returned as `np.nan`.
See Also
--------
stft : Short-time Fourier Transform
reassigned_spectrogram : Time-frequency reassigned spectrogram
Examples
--------
>>> y, sr = librosa.load(librosa.ex('trumpet'))
>>> frequencies, S = librosa.core.spectrum.__reassign_frequencies(y, sr=sr)
>>> frequencies
array([[0.000e+00, 0.000e+00, ..., 0.000e+00, 0.000e+00],
[3.628e+00, 4.698e+00, ..., 1.239e+01, 1.072e+01],
...,
[1.101e+04, 1.102e+04, ..., 1.105e+04, 1.102e+04],
[1.102e+04, 1.102e+04, ..., 1.102e+04, 1.102e+04]])
"""
# retrieve window samples if needed so that the window derivative can be
# calculated
if win_length is None:
win_length = n_fft
window = get_window(window, win_length, fftbins=True)
window = util.pad_center(window, size=n_fft)
if S is None:
if dtype is None:
dtype = util.dtype_r2c(y.dtype)
S_h = stft(
y=y,
n_fft=n_fft,
hop_length=hop_length,
window=window,
center=center,
dtype=dtype,
pad_mode=pad_mode,
)
else:
if dtype is None:
dtype = S.dtype
S_h = S
# cyclic gradient to correctly handle edges of a periodic window
window_derivative = util.cyclic_gradient(window)
S_dh = stft(
y=y,
n_fft=n_fft,
hop_length=hop_length,
window=window_derivative,
center=center,
dtype=dtype,
pad_mode=pad_mode,
)
# equation 5.20 of Flandrin, Auger, & Chassande-Mottin 2002
# the sign of the correction is reversed in some papers - see Plante,
# Meyer, & Ainsworth 1998 pp. 283-284
with np.errstate(invalid="ignore"):
# We can ignore divide-by-zero here because NaN is an acceptable correction value
correction = -np.imag(S_dh / S_h)
freqs = convert.fft_frequencies(sr=sr, n_fft=n_fft)
freqs = util.expand_to(freqs, ndim=correction.ndim, axes=-2) + correction * (
0.5 * sr / np.pi
)
return freqs, S_h
def __reassign_times(
y: np.ndarray,
sr: float = 22050,
S: Optional[np.ndarray] = None,
n_fft: int = 2048,
hop_length: Optional[int] = None,
win_length: Optional[int] = None,
window: _WindowSpec = "hann",
center: bool = True,
dtype: Optional[DTypeLike] = None,
pad_mode: _PadModeSTFT = "constant",
) -> Tuple[np.ndarray, np.ndarray]:
"""Time reassignments based on a spectrogram representation.
The reassignment vector is calculated using equation 5.23 in Flandrin,
Auger, & Chassande-Mottin 2002::
t_reassigned = t + np.real(S_th/S_h)
where ``S_h`` is the complex STFT calculated using the original window, and
``S_th`` is the complex STFT calculated using the original window multiplied
by the time offset from the window center.
See `reassigned_spectrogram` for references.
It is recommended to use ``pad_mode="constant"`` (zero padding) or else
``center=False``, rather than the defaults. Time reassignment assumes that
the energy in each FFT bin is associated with exactly one impulse event.
Reflection padding at the edges of the signal may invalidate the reassigned
estimates in the boundary frames.
Parameters
----------
y : np.ndarray [shape=(..., n,)], real-valued
audio time series. Multi-channel is supported.
sr : number > 0 [scalar]
sampling rate of ``y``
S : np.ndarray [shape=(..., d, t)] or None
(optional) complex STFT calculated using the other arguments provided
to `__reassign_times`
n_fft : int > 0 [scalar]
FFT window size. Defaults to 2048.
hop_length : int > 0 [scalar]
hop length, number samples between subsequent frames.
If not supplied, defaults to ``win_length // 4``.
win_length : int > 0, <= n_fft
Window length. Defaults to ``n_fft``.
See `stft` for details.
window : string, tuple, number, function, or np.ndarray [shape=(n_fft,)]
- a window specification (string, tuple, number);
see `scipy.signal.get_window`
- a window function, such as `scipy.signal.windows.hann`
- a user-specified window vector of length ``n_fft``
See `stft` for details.
.. see also:: `filters.get_window`
center : boolean
- If ``True``, the signal ``y`` is padded so that frame
``S[:, t]`` is centered at ``y[t * hop_length]``.
- If ``False``, then ``S[:, t]`` begins at ``y[t * hop_length]``.
dtype : numeric type
Complex numeric type for ``S``. Default is inferred to match
the precision of the input signal.
pad_mode : string
If ``center=True``, the padding mode to use at the edges of the signal.
By default, STFT uses zero padding.
Returns
-------
times : np.ndarray [shape=(..., 1 + n_fft/2, t), dtype=real]
Reassigned times:
``times[f, t]`` is the time for bin ``f``, frame ``t``.
S : np.ndarray [shape=(..., 1 + n_fft/2, t), dtype=complex]
Short-time Fourier transform
Warns
-----
RuntimeWarning
Time estimates with zero support will produce a divide-by-zero warning
and will be returned as `np.nan`.
See Also
--------
stft : Short-time Fourier Transform
reassigned_spectrogram : Time-frequency reassigned spectrogram
Examples
--------
>>> y, sr = librosa.load(librosa.ex('trumpet'))
>>> times, S = librosa.core.spectrum.__reassign_times(y, sr=sr)
>>> times
array([[ 2.268e-05, 1.144e-02, ..., 5.332e+00, 5.333e+00],
[ 2.268e-05, 1.451e-02, ..., 5.334e+00, 5.333e+00],
...,
[ 2.268e-05, -6.177e-04, ..., 5.368e+00, 5.327e+00],
[ 2.268e-05, 1.420e-03, ..., 5.307e+00, 5.328e+00]])
"""
# retrieve window samples if needed so that the time-weighted window can be
# calculated
if win_length is None:
win_length = n_fft
window = get_window(window, win_length, fftbins=True)
window = util.pad_center(window, size=n_fft)
# retrieve hop length if needed so that the frame times can be calculated
if hop_length is None:
hop_length = int(win_length // 4)
if S is None:
if dtype is None:
dtype = util.dtype_r2c(y.dtype)
S_h = stft(
y=y,
n_fft=n_fft,
hop_length=hop_length,
window=window,
center=center,
dtype=dtype,
pad_mode=pad_mode,
)
else:
if dtype is None:
dtype = S.dtype
S_h = S
# calculate window weighted by time
half_width = n_fft // 2
window_times: np.ndarray
if n_fft % 2:
window_times = np.arange(-half_width, half_width + 1)
else:
window_times = np.arange(0.5 - half_width, half_width)
window_time_weighted = window * window_times
S_th = stft(
y=y,
n_fft=n_fft,
hop_length=hop_length,
window=window_time_weighted,
center=center,
dtype=dtype,
pad_mode=pad_mode,
)
# equation 5.23 of Flandrin, Auger, & Chassande-Mottin 2002
# the sign of the correction is reversed in some papers - see Plante,
# Meyer, & Ainsworth 1998 pp. 283-284
with np.errstate(invalid="ignore"):
# We can ignore divide-by-zero here because NaN is an acceptable correction value
correction = np.real(S_th / S_h)
if center:
pad_length = None
else:
pad_length = n_fft
times = convert.frames_to_time(
np.arange(S_h.shape[-1]), sr=sr, hop_length=hop_length, n_fft=pad_length
)
times = util.expand_to(times, ndim=correction.ndim, axes=-1) + correction / sr
return times, S_h
[docs]def reassigned_spectrogram(
y: np.ndarray,
*,
sr: float = 22050,
S: Optional[np.ndarray] = None,
n_fft: int = 2048,
hop_length: Optional[int] = None,
win_length: Optional[int] = None,
window: _WindowSpec = "hann",
center: bool = True,
reassign_frequencies: bool = True,
reassign_times: bool = True,
ref_power: Union[float, Callable] = 1e-6,
fill_nan: bool = False,
clip: bool = True,
dtype: Optional[DTypeLike] = None,
pad_mode: _PadModeSTFT = "constant",
) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
r"""Time-frequency reassigned spectrogram.
The reassignment vectors are calculated using equations 5.20 and 5.23 in
[#]_::
t_reassigned = t + np.real(S_th/S_h)
omega_reassigned = omega - np.imag(S_dh/S_h)
where ``S_h`` is the complex STFT calculated using the original window,
``S_dh`` is the complex STFT calculated using the derivative of the original
window, and ``S_th`` is the complex STFT calculated using the original window
multiplied by the time offset from the window center. See [#]_ for
additional algorithms, and [#]_ and [#]_ for history and discussion of the
method.
.. [#] Flandrin, P., Auger, F., & Chassande-Mottin, E. (2002).
Time-Frequency reassignment: From principles to algorithms. In
Applications in Time-Frequency Signal Processing (Vol. 10, pp.
179-204). CRC Press.
.. [#] Fulop, S. A., & Fitz, K. (2006). Algorithms for computing the
time-corrected instantaneous frequency (reassigned) spectrogram, with
applications. The Journal of the Acoustical Society of America, 119(1),
360. doi:10.1121/1.2133000
.. [#] Auger, F., Flandrin, P., Lin, Y.-T., McLaughlin, S., Meignen, S.,
Oberlin, T., & Wu, H.-T. (2013). Time-Frequency Reassignment and
Synchrosqueezing: An Overview. IEEE Signal Processing Magazine, 30(6),
32-41. doi:10.1109/MSP.2013.2265316
.. [#] Hainsworth, S., Macleod, M. (2003). Time-frequency reassignment: a
review and analysis. Tech. Rep. CUED/FINFENG/TR.459, Cambridge
University Engineering Department
Parameters
----------
y : np.ndarray [shape=(..., n)], real-valued
audio time series. Multi-channel is supported.
sr : number > 0 [scalar]
sampling rate of ``y``
S : np.ndarray [shape=(..., d, t)] or None
(optional) complex STFT calculated using the other arguments provided
to ``reassigned_spectrogram``
n_fft : int > 0 [scalar]
FFT window size. Defaults to 2048.
hop_length : int > 0 [scalar]
hop length, number samples between subsequent frames.
If not supplied, defaults to ``win_length // 4``.
win_length : int > 0, <= n_fft
Window length. Defaults to ``n_fft``.
See `stft` for details.
window : string, tuple, number, function, or np.ndarray [shape=(n_fft,)]
- a window specification (string, tuple, number);
see `scipy.signal.get_window`
- a window function, such as `scipy.signal.windows.hann`
- a user-specified window vector of length ``n_fft``
See `stft` for details.
.. see also:: `filters.get_window`
center : boolean
- If ``True`` (default), the signal ``y`` is padded so that frame
``S[:, t]`` is centered at ``y[t * hop_length]``. See `Notes` for
recommended usage in this function.
- If ``False``, then ``S[:, t]`` begins at ``y[t * hop_length]``.
reassign_frequencies : boolean
- If ``True`` (default), the returned frequencies will be instantaneous
frequency estimates.
- If ``False``, the returned frequencies will be a read-only view of the
STFT bin frequencies for all frames.
reassign_times : boolean
- If ``True`` (default), the returned times will be corrected
(reassigned) time estimates for each bin.
- If ``False``, the returned times will be a read-only view of the STFT
frame times for all bins.
ref_power : float >= 0 or callable
Minimum power threshold for estimating time-frequency reassignments.
Any bin with ``np.abs(S[f, t])**2 < ref_power`` will be returned as
`np.nan` in both frequency and time, unless ``fill_nan`` is ``True``. If 0
is provided, then only bins with zero power will be returned as
`np.nan` (unless ``fill_nan=True``).
fill_nan : boolean
- If ``False`` (default), the frequency and time reassignments for bins
below the power threshold provided in ``ref_power`` will be returned as
`np.nan`.
- If ``True``, the frequency and time reassignments for these bins will
be returned as the bin center frequencies and frame times.
clip : boolean
- If ``True`` (default), estimated frequencies outside the range
`[0, 0.5 * sr]` or times outside the range `[0, len(y) / sr]` will be
clipped to those ranges.
- If ``False``, estimated frequencies and times beyond the bounds of the
spectrogram may be returned.
dtype : numeric type
Complex numeric type for STFT calculation. Default is inferred to match
the precision of the input signal.
pad_mode : string
If ``center=True``, the padding mode to use at the edges of the signal.
By default, STFT uses zero padding.
Returns
-------
freqs, times, mags : np.ndarray [shape=(..., 1 + n_fft/2, t), dtype=real]
Instantaneous frequencies:
``freqs[..., f, t]`` is the frequency for bin ``f``, frame ``t``.
If ``reassign_frequencies=False``, this will instead be a read-only array
of the same shape containing the bin center frequencies for all frames.
Reassigned times:
``times[..., f, t]`` is the time for bin ``f``, frame ``t``.
If ``reassign_times=False``, this will instead be a read-only array of
the same shape containing the frame times for all bins.
Magnitudes from short-time Fourier transform:
``mags[..., f, t]`` is the magnitude for bin ``f``, frame ``t``.
Warns
-----
RuntimeWarning
Frequency or time estimates with zero support will produce a
divide-by-zero warning, and will be returned as `np.nan` unless
``fill_nan=True``.
See Also
--------
stft : Short-time Fourier Transform
Notes
-----
It is recommended to use ``center=False`` with this function rather than the
librosa default ``True``. Unlike ``stft``, reassigned times are not aligned to
the left or center of each frame, so padding the signal does not affect the
meaning of the reassigned times. However, reassignment assumes that the
energy in each FFT bin is associated with exactly one signal component and
impulse event.
If ``reassign_times`` is ``False``, the frame times that are returned will be
aligned to the left or center of the frame, depending on the value of
``center``. In this case, if ``center`` is ``True``, then ``pad_mode="wrap"`` is
recommended for valid estimation of the instantaneous frequencies in the
boundary frames.
Examples
--------
>>> import matplotlib.pyplot as plt
>>> amin = 1e-10
>>> n_fft = 64
>>> sr = 4000
>>> y = 1e-3 * librosa.clicks(times=[0.3], sr=sr, click_duration=1.0,
... click_freq=1200.0, length=8000) +\
... 1e-3 * librosa.clicks(times=[1.5], sr=sr, click_duration=0.5,
... click_freq=400.0, length=8000) +\
... 1e-3 * librosa.chirp(fmin=200, fmax=1600, sr=sr, duration=2.0) +\
... 1e-6 * np.random.randn(2*sr)
>>> freqs, times, mags = librosa.reassigned_spectrogram(y=y, sr=sr,
... n_fft=n_fft)
>>> mags_db = librosa.amplitude_to_db(mags, ref=np.max)
>>> fig, ax = plt.subplots(nrows=2, sharex=True, sharey=True)
>>> img = librosa.display.specshow(mags_db, x_axis="s", y_axis="linear", sr=sr,
... hop_length=n_fft//4, ax=ax[0])
>>> ax[0].set(title="Spectrogram", xlabel=None)
>>> ax[0].label_outer()
>>> ax[1].scatter(times, freqs, c=mags_db, cmap="magma", alpha=0.1, s=5)
>>> ax[1].set_title("Reassigned spectrogram")
>>> fig.colorbar(img, ax=ax, format="%+2.f dB")
"""
if not callable(ref_power) and ref_power < 0:
raise ParameterError("ref_power must be non-negative or callable.")
if not reassign_frequencies and not reassign_times:
raise ParameterError("reassign_frequencies or reassign_times must be True.")
if win_length is None:
win_length = n_fft
if hop_length is None:
hop_length = int(win_length // 4)
# frequency and time reassignment if requested
if reassign_frequencies:
freqs, S = __reassign_frequencies(
y=y,
sr=sr,
S=S,
n_fft=n_fft,
hop_length=hop_length,
win_length=win_length,
window=window,
center=center,
dtype=dtype,
pad_mode=pad_mode,
)
if reassign_times:
times, S = __reassign_times(
y=y,
sr=sr,
S=S,
n_fft=n_fft,
hop_length=hop_length,
win_length=win_length,
window=window,
center=center,
dtype=dtype,
pad_mode=pad_mode,
)
assert S is not None
mags: np.ndarray = np.abs(S)
# clean up reassignment issues: divide-by-zero, bins with near-zero power,
# and estimates outside the spectrogram bounds
# retrieve bin frequencies and frame times to replace missing estimates
if fill_nan or not reassign_frequencies or not reassign_times:
if center:
pad_length = None
else:
pad_length = n_fft
bin_freqs = convert.fft_frequencies(sr=sr, n_fft=n_fft)
frame_times = convert.frames_to_time(
frames=np.arange(S.shape[-1]),
sr=sr,
hop_length=hop_length,
n_fft=pad_length,
)
# find bins below the power threshold
# reassigned bins with zero power will already be NaN
if callable(ref_power):
ref_p = ref_power(mags**2)
else:
ref_p = ref_power
mags_low = np.less(mags, ref_p**0.5, where=~np.isnan(mags))
# for reassigned estimates, optionally set thresholded bins to NaN, return
# bin frequencies and frame times in place of NaN generated by
# divide-by-zero and power threshold, and clip to spectrogram bounds
if reassign_frequencies:
if ref_p > 0:
freqs[mags_low] = np.nan
if fill_nan:
freqs = np.where(np.isnan(freqs), bin_freqs[:, np.newaxis], freqs)
if clip:
np.clip(freqs, 0, sr / 2.0, out=freqs)
# or if reassignment was not requested, return bin frequencies and frame
# times for every cell is the spectrogram
else:
freqs = np.broadcast_to(bin_freqs[:, np.newaxis], S.shape)
if reassign_times:
if ref_p > 0:
times[mags_low] = np.nan
if fill_nan:
times = np.where(np.isnan(times), frame_times[np.newaxis, :], times)
if clip:
np.clip(times, 0, y.shape[-1] / float(sr), out=times)
else:
times = np.broadcast_to(frame_times[np.newaxis, :], S.shape)
return freqs, times, mags
[docs]def magphase(D: np.ndarray, *, power: float = 1) -> Tuple[np.ndarray, np.ndarray]:
"""Separate a complex-valued spectrogram D into its magnitude (S)
and phase (P) components, so that ``D = S * P``.
Parameters
----------
D : np.ndarray [shape=(..., d, t), dtype=complex]
complex-valued spectrogram
power : float > 0
Exponent for the magnitude spectrogram,
e.g., 1 for energy, 2 for power, etc.
Returns
-------
D_mag : np.ndarray [shape=(..., d, t), dtype=real]
magnitude of ``D``, raised to ``power``
D_phase : np.ndarray [shape=(..., d, t), dtype=complex]
``exp(1.j * phi)`` where ``phi`` is the phase of ``D``
Examples
--------
>>> y, sr = librosa.load(librosa.ex('trumpet'))
>>> D = librosa.stft(y)
>>> magnitude, phase = librosa.magphase(D)
>>> magnitude
array([[5.395e-03, 3.332e-03, ..., 9.862e-07, 1.201e-05],
[3.244e-03, 2.690e-03, ..., 9.536e-07, 1.201e-05],
...,
[7.523e-05, 3.722e-05, ..., 1.188e-04, 1.031e-03],
[7.640e-05, 3.944e-05, ..., 5.180e-04, 1.346e-03]],
dtype=float32)
>>> phase
array([[ 1. +0.000e+00j, 1. +0.000e+00j, ...,
-1. -8.742e-08j, -1. -8.742e-08j],
[-1. -8.742e-08j, -0.775-6.317e-01j, ...,
-0.885-4.648e-01j, 0.472-8.815e-01j],
...,
[ 1. -4.342e-12j, 0.028-9.996e-01j, ...,
-0.222-9.751e-01j, -0.75 -6.610e-01j],
[-1. -8.742e-08j, -1. -8.742e-08j, ...,
1. +0.000e+00j, 1. +0.000e+00j]], dtype=complex64)
Or get the phase angle (in radians)
>>> np.angle(phase)
array([[ 0.000e+00, 0.000e+00, ..., -3.142e+00, -3.142e+00],
[-3.142e+00, -2.458e+00, ..., -2.658e+00, -1.079e+00],
...,
[-4.342e-12, -1.543e+00, ..., -1.794e+00, -2.419e+00],
[-3.142e+00, -3.142e+00, ..., 0.000e+00, 0.000e+00]],
dtype=float32)
"""
mag = np.abs(D)
# Prevent NaNs and return magnitude 0, phase 1+0j for zero
zeros_to_ones = mag == 0
mag_nonzero = mag + zeros_to_ones
# Compute real and imaginary separately, because complex division can
# produce NaNs when denormalized numbers are involved (< ~2e-39 for
# complex64, ~5e-309 for complex128)
phase = np.empty_like(D, dtype=util.dtype_r2c(D.dtype))
phase.real = D.real / mag_nonzero + zeros_to_ones
phase.imag = D.imag / mag_nonzero
mag **= power
return mag, phase
[docs]def phase_vocoder(
D: np.ndarray,
*,
rate: float,
hop_length: Optional[int] = None,
n_fft: Optional[int] = None,
) -> np.ndarray:
"""Phase vocoder. Given an STFT matrix D, speed up by a factor of ``rate``
Based on the implementation provided by [#]_.
This is a simplified implementation, intended primarily for
reference and pedagogical purposes. It makes no attempt to
handle transients, and is likely to produce many audible
artifacts. For a higher quality implementation, we recommend
the RubberBand library [#]_ and its Python wrapper `pyrubberband`.
.. [#] Ellis, D. P. W. "A phase vocoder in Matlab."
Columbia University, 2002.
https://www.ee.columbia.edu/~dpwe/resources/matlab/pvoc/
.. [#] https://breakfastquay.com/rubberband/
Examples
--------
>>> # Play at double speed
>>> y, sr = librosa.load(librosa.ex('trumpet'))
>>> D = librosa.stft(y, n_fft=2048, hop_length=512)
>>> D_fast = librosa.phase_vocoder(D, rate=2.0, hop_length=512)
>>> y_fast = librosa.istft(D_fast, hop_length=512)
>>> # Or play at 1/3 speed
>>> y, sr = librosa.load(librosa.ex('trumpet'))
>>> D = librosa.stft(y, n_fft=2048, hop_length=512)
>>> D_slow = librosa.phase_vocoder(D, rate=1./3, hop_length=512)
>>> y_slow = librosa.istft(D_slow, hop_length=512)
Parameters
----------
D : np.ndarray [shape=(..., d, t), dtype=complex]
STFT matrix
rate : float > 0 [scalar]
Speed-up factor: ``rate > 1`` is faster, ``rate < 1`` is slower.
hop_length : int > 0 [scalar] or None
The number of samples between successive columns of ``D``.
If None, defaults to ``n_fft//4 = (D.shape[0]-1)//2``
n_fft : int > 0 or None
The number of samples per frame in D.
By default (None), this will be inferred from the shape of D.
However, if D was constructed using an odd-length window, the correct
frame length can be specified here.
Returns
-------
D_stretched : np.ndarray [shape=(..., d, t / rate), dtype=complex]
time-stretched STFT
See Also
--------
pyrubberband
"""
if n_fft is None:
n_fft = 2 * (D.shape[-2] - 1)
if hop_length is None:
hop_length = int(n_fft // 4)
time_steps = np.arange(0, D.shape[-1], rate, dtype=np.float64)
# Create an empty output array
shape = list(D.shape)
shape[-1] = len(time_steps)
d_stretch = np.zeros_like(D, shape=shape)
# Expected phase advance in each bin per frame
phi_advance = hop_length * convert.fft_frequencies(sr=2 * np.pi, n_fft=n_fft)
# Phase accumulator; initialize to the first sample
phase_acc = np.angle(D[..., 0])
# Pad 0 columns to simplify boundary logic
padding = [(0, 0) for _ in D.shape]
padding[-1] = (0, 2)
D = np.pad(D, padding, mode="constant")
for t, step in enumerate(time_steps):
columns = D[..., int(step) : int(step + 2)]
# Weighting for linear magnitude interpolation
alpha = np.mod(step, 1.0)
mag = (1.0 - alpha) * np.abs(columns[..., 0]) + alpha * np.abs(columns[..., 1])
# Store to output array
d_stretch[..., t] = util.phasor(phase_acc, mag=mag)
# Compute phase advance
dphase = np.angle(columns[..., 1]) - np.angle(columns[..., 0]) - phi_advance
# Wrap to -pi:pi range
dphase = dphase - 2.0 * np.pi * np.round(dphase / (2.0 * np.pi))
# Accumulate phase
phase_acc += phi_advance + dphase
return d_stretch
[docs]@cache(level=20)
def iirt(
y: np.ndarray,
*,
sr: float = 22050,
win_length: int = 2048,
hop_length: Optional[int] = None,
center: bool = True,
tuning: float = 0.0,
pad_mode: _PadMode = "constant",
flayout: str = "sos",
res_type: str = "soxr_hq",
**kwargs: Any,
) -> np.ndarray:
r"""Time-frequency representation using IIR filters
This function will return a time-frequency representation
using a multirate filter bank consisting of IIR filters. [#]_
First, ``y`` is resampled as needed according to the provided ``sample_rates``.
Then, a filterbank with with ``n`` band-pass filters is designed.
The resampled input signals are processed by the filterbank as a whole.
(`scipy.signal.filtfilt` resp. `sosfiltfilt` is used to make the phase linear.)
The output of the filterbank is cut into frames.
For each band, the short-time mean-square power (STMSP) is calculated by
summing ``win_length`` subsequent filtered time samples.
When called with the default set of parameters, it will generate the TF-representation
(pitch filterbank):
* 85 filters with MIDI pitches [24, 108] as ``center_freqs``.
* each filter having a bandwidth of one semitone.
.. [#] Müller, Meinard.
"Information Retrieval for Music and Motion."
Springer Verlag. 2007.
Parameters
----------
y : np.ndarray [shape=(..., n)]
audio time series. Multi-channel is supported.
sr : number > 0 [scalar]
sampling rate of ``y``
win_length : int > 0, <= n_fft
Window length.
hop_length : int > 0 [scalar]
Hop length, number samples between subsequent frames.
If not supplied, defaults to ``win_length // 4``.
center : boolean
- If ``True``, the signal ``y`` is padded so that frame
``D[..., :, t]`` is centered at ``y[t * hop_length]``.
- If ``False``, then `D[..., :, t]`` begins at ``y[t * hop_length]``
tuning : float [scalar]
Tuning deviation from A440 in fractions of a bin.
pad_mode : string
If ``center=True``, the padding mode to use at the edges of the signal.
By default, this function uses zero padding.
flayout : string
- If `sos` (default), a series of second-order filters is used for filtering with `scipy.signal.sosfiltfilt`.
Minimizes numerical precision errors for high-order filters, but is slower.
- If `ba`, the standard difference equation is used for filtering with `scipy.signal.filtfilt`.
Can be unstable for high-order filters.
res_type : string
The resampling mode. See `librosa.resample` for details.
**kwargs : additional keyword arguments
Additional arguments for `librosa.filters.semitone_filterbank`
(e.g., could be used to provide another set of ``center_freqs`` and ``sample_rates``).
Returns
-------
bands_power : np.ndarray [shape=(..., n, t), dtype=dtype]
Short-time mean-square power for the input signal.
Raises
------
ParameterError
If ``flayout`` is not None, `ba`, or `sos`.
See Also
--------
librosa.filters.semitone_filterbank
librosa.filters.mr_frequencies
librosa.cqt
scipy.signal.filtfilt
scipy.signal.sosfiltfilt
Examples
--------
>>> import matplotlib.pyplot as plt
>>> y, sr = librosa.load(librosa.ex('trumpet'), duration=3)
>>> D = np.abs(librosa.iirt(y, sr=sr))
>>> C = np.abs(librosa.cqt(y=y, sr=sr))
>>> fig, ax = plt.subplots(nrows=2, sharex=True, sharey=True)
>>> img = librosa.display.specshow(librosa.amplitude_to_db(C, ref=np.max),
... y_axis='cqt_hz', x_axis='time', ax=ax[0])
>>> ax[0].set(title='Constant-Q transform')
>>> ax[0].label_outer()
>>> img = librosa.display.specshow(librosa.amplitude_to_db(D, ref=np.max),
... y_axis='cqt_hz', x_axis='time', ax=ax[1])
>>> ax[1].set_title('Semitone spectrogram (iirt)')
>>> fig.colorbar(img, ax=ax, format="%+2.0f dB")
"""
if flayout not in ("ba", "sos"):
raise ParameterError(f"Unsupported flayout={flayout}")
# check audio input
util.valid_audio(y, mono=False)
# Set the default hop, if it's not already specified
if hop_length is None:
hop_length = win_length // 4
# Pad the time series so that frames are centered
if center:
padding = [(0, 0) for _ in y.shape]
padding[-1] = (win_length // 2, win_length // 2)
y = np.pad(y, padding, mode=pad_mode)
# get the semitone filterbank
filterbank_ct, sample_rates = semitone_filterbank(
tuning=tuning, flayout=flayout, **kwargs
)
# create three downsampled versions of the audio signal
y_resampled = []
y_srs = np.unique(sample_rates)
for cur_sr in y_srs:
y_resampled.append(resample(y, orig_sr=sr, target_sr=cur_sr, res_type=res_type))
# Compute the number of frames that will fit. The end may get truncated.
n_frames = int(1 + (y.shape[-1] - win_length) // hop_length)
# Pre-allocate the output array
shape = list(y.shape)
# Time dimension reduces to n_frames
shape[-1] = n_frames
# Insert a new axis at position -2 for filter response
shape.insert(-1, len(filterbank_ct))
bands_power = np.empty_like(y, shape=shape)
slices: List[Union[int, slice]] = [slice(None) for _ in bands_power.shape]
for i, (cur_sr, cur_filter) in enumerate(zip(sample_rates, filterbank_ct)):
slices[-2] = i
# filter the signal
cur_sr_idx = np.flatnonzero(y_srs == cur_sr)[0]
if flayout == "ba":
cur_filter_output = scipy.signal.filtfilt(
cur_filter[0], cur_filter[1], y_resampled[cur_sr_idx], axis=-1
)
elif flayout == "sos":
cur_filter_output = scipy.signal.sosfiltfilt(
cur_filter, y_resampled[cur_sr_idx], axis=-1
)
factor = sr / cur_sr
hop_length_STMSP = hop_length / factor
win_length_STMSP_round = int(round(win_length / factor))
# hop_length_STMSP is used here as a floating-point number.
# The discretization happens at the end to avoid accumulated rounding errors.
start_idx = np.arange(
0, cur_filter_output.shape[-1] - win_length_STMSP_round, hop_length_STMSP
)
if len(start_idx) < n_frames:
min_length = (
int(np.ceil(n_frames * hop_length_STMSP)) + win_length_STMSP_round
)
cur_filter_output = util.fix_length(cur_filter_output, size=min_length)
start_idx = np.arange(
0,
cur_filter_output.shape[-1] - win_length_STMSP_round,
hop_length_STMSP,
)
start_idx = np.round(start_idx).astype(int)[:n_frames]
idx = np.add.outer(start_idx, np.arange(win_length_STMSP_round))
bands_power[tuple(slices)] = factor * np.sum(
cur_filter_output[..., idx] ** 2, axis=-1
)
return bands_power
@overload
def power_to_db(
S: _ComplexLike_co,
*,
ref: Union[float, Callable] = ...,
amin: float = ...,
top_db: Optional[float] = ...,
) -> np.floating[Any]:
...
@overload
def power_to_db(
S: _SequenceLike[_ComplexLike_co],
*,
ref: Union[float, Callable] = ...,
amin: float = ...,
top_db: Optional[float] = ...,
) -> np.ndarray:
...
@overload
def power_to_db(
S: _ScalarOrSequence[_ComplexLike_co],
*,
ref: Union[float, Callable] = ...,
amin: float = ...,
top_db: Optional[float] = ...,
) -> Union[np.floating[Any], np.ndarray]:
...
[docs]@cache(level=30)
def power_to_db(
S: _ScalarOrSequence[_ComplexLike_co],
*,
ref: Union[float, Callable] = 1.0,
amin: float = 1e-10,
top_db: Optional[float] = 80.0,
) -> Union[np.floating[Any], np.ndarray]:
"""Convert a power spectrogram (amplitude squared) to decibel (dB) units
This computes the scaling ``10 * log10(S / ref)`` in a numerically
stable way.
Parameters
----------
S : np.ndarray
input power
ref : scalar or callable
If scalar, the amplitude ``abs(S)`` is scaled relative to ``ref``::
10 * log10(S / ref)
Zeros in the output correspond to positions where ``S == ref``.
If callable, the reference value is computed as ``ref(S)``.
amin : float > 0 [scalar]
minimum threshold for ``abs(S)`` and ``ref``
top_db : float >= 0 [scalar]
threshold the output at ``top_db`` below the peak:
``max(10 * log10(S/ref)) - top_db``
Returns
-------
S_db : np.ndarray
``S_db ~= 10 * log10(S) - 10 * log10(ref)``
See Also
--------
perceptual_weighting
db_to_power
amplitude_to_db
db_to_amplitude
Notes
-----
This function caches at level 30.
Examples
--------
Get a power spectrogram from a waveform ``y``
>>> y, sr = librosa.load(librosa.ex('trumpet'))
>>> S = np.abs(librosa.stft(y))
>>> librosa.power_to_db(S**2)
array([[-41.809, -41.809, ..., -41.809, -41.809],
[-41.809, -41.809, ..., -41.809, -41.809],
...,
[-41.809, -41.809, ..., -41.809, -41.809],
[-41.809, -41.809, ..., -41.809, -41.809]], dtype=float32)
Compute dB relative to peak power
>>> librosa.power_to_db(S**2, ref=np.max)
array([[-80., -80., ..., -80., -80.],
[-80., -80., ..., -80., -80.],
...,
[-80., -80., ..., -80., -80.],
[-80., -80., ..., -80., -80.]], dtype=float32)
Or compare to median power
>>> librosa.power_to_db(S**2, ref=np.median)
array([[16.578, 16.578, ..., 16.578, 16.578],
[16.578, 16.578, ..., 16.578, 16.578],
...,
[16.578, 16.578, ..., 16.578, 16.578],
[16.578, 16.578, ..., 16.578, 16.578]], dtype=float32)
And plot the results
>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots(nrows=2, sharex=True, sharey=True)
>>> imgpow = librosa.display.specshow(S**2, sr=sr, y_axis='log', x_axis='time',
... ax=ax[0])
>>> ax[0].set(title='Power spectrogram')
>>> ax[0].label_outer()
>>> imgdb = librosa.display.specshow(librosa.power_to_db(S**2, ref=np.max),
... sr=sr, y_axis='log', x_axis='time', ax=ax[1])
>>> ax[1].set(title='Log-Power spectrogram')
>>> fig.colorbar(imgpow, ax=ax[0])
>>> fig.colorbar(imgdb, ax=ax[1], format="%+2.0f dB")
"""
S = np.asarray(S)
if amin <= 0:
raise ParameterError("amin must be strictly positive")
if np.issubdtype(S.dtype, np.complexfloating):
warnings.warn(
"power_to_db was called on complex input so phase "
"information will be discarded. To suppress this warning, "
"call power_to_db(np.abs(D)**2) instead.",
stacklevel=2,
)
magnitude = np.abs(S)
else:
magnitude = S
if callable(ref):
# User supplied a function to calculate reference power
ref_value = ref(magnitude)
else:
ref_value = np.abs(ref)
log_spec: np.ndarray = 10.0 * np.log10(np.maximum(amin, magnitude))
log_spec -= 10.0 * np.log10(np.maximum(amin, ref_value))
if top_db is not None:
if top_db < 0:
raise ParameterError("top_db must be non-negative")
log_spec = np.maximum(log_spec, log_spec.max() - top_db)
return log_spec
@overload
def db_to_power(
S_db: _FloatLike_co,
*,
ref: float = ...,
) -> np.floating[Any]:
...
@overload
def db_to_power(
S_db: np.ndarray,
*,
ref: float = ...,
) -> np.ndarray:
...
@overload
def db_to_power(
S_db: Union[_FloatLike_co, np.ndarray],
*,
ref: float = ...,
) -> Union[np.floating[Any], np.ndarray]:
...
[docs]@cache(level=30)
def db_to_power(S_db: Union[_FloatLike_co, np.ndarray], *, ref: float = 1.0) -> Union[np.floating[Any], np.ndarray]:
"""Convert dB-scale values to a power values.
This effectively inverts ``power_to_db``::
db_to_power(S_db) ~= ref * 10.0**(S_db / 10)
Parameters
----------
S_db : np.ndarray
dB-scaled values
ref : number > 0
Reference power: output will be scaled by this value
Returns
-------
S : np.ndarray
Power values
Notes
-----
This function caches at level 30.
"""
return ref * np.power(10.0, S_db * 0.1)
@overload
def amplitude_to_db(
S: _ComplexLike_co,
*,
ref: Union[float, Callable] = ...,
amin: float = ...,
top_db: Optional[float] = ...,
) -> np.floating[Any]:
...
@overload
def amplitude_to_db(
S: _SequenceLike[_ComplexLike_co],
*,
ref: Union[float, Callable] = ...,
amin: float = ...,
top_db: Optional[float] = ...,
) -> np.ndarray:
...
@overload
def amplitude_to_db(
S: _ScalarOrSequence[_ComplexLike_co],
*,
ref: Union[float, Callable] = ...,
amin: float = ...,
top_db: Optional[float] = ...,
) -> Union[np.floating[Any], np.ndarray]:
...
[docs]@cache(level=30)
def amplitude_to_db(
S: _ScalarOrSequence[_ComplexLike_co],
*,
ref: Union[float, Callable] = 1.0,
amin: float = 1e-5,
top_db: Optional[float] = 80.0,
) -> Union[np.floating[Any], np.ndarray]:
"""Convert an amplitude spectrogram to dB-scaled spectrogram.
This is equivalent to ``power_to_db(S**2, ref=ref**2, amin=amin**2, top_db=top_db)``,
but is provided for convenience.
Parameters
----------
S : np.ndarray
input amplitude
ref : scalar or callable
If scalar, the amplitude ``abs(S)`` is scaled relative to ``ref``:
``20 * log10(S / ref)``.
Zeros in the output correspond to positions where ``S == ref``.
If callable, the reference value is computed as ``ref(S)``.
amin : float > 0 [scalar]
minimum threshold for ``S`` and ``ref``
top_db : float >= 0 [scalar]
threshold the output at ``top_db`` below the peak:
``max(20 * log10(S/ref)) - top_db``
Returns
-------
S_db : np.ndarray
``S`` measured in dB
See Also
--------
power_to_db, db_to_amplitude
Notes
-----
This function caches at level 30.
"""
S = np.asarray(S)
if np.issubdtype(S.dtype, np.complexfloating):
warnings.warn(
"amplitude_to_db was called on complex input so phase "
"information will be discarded. To suppress this warning, "
"call amplitude_to_db(np.abs(S)) instead.",
stacklevel=2,
)
magnitude = np.abs(S)
if callable(ref):
# User supplied a function to calculate reference power
ref_value = ref(magnitude)
else:
ref_value = np.abs(ref)
out_array = magnitude if isinstance(magnitude, np.ndarray) else None
power = np.square(magnitude, out=out_array)
db: np.ndarray = power_to_db(power, ref=ref_value**2, amin=amin**2, top_db=top_db)
return db
@overload
def db_to_amplitude(
S_db: _FloatLike_co,
*,
ref: float = ...,
) -> np.floating[Any]:
...
@overload
def db_to_amplitude(
S_db: np.ndarray,
*,
ref: float = ...,
) -> np.ndarray:
...
@overload
def db_to_amplitude(
S_db: Union[_FloatLike_co, np.ndarray],
*,
ref: float = ...,
) -> Union[np.floating[Any], np.ndarray]:
...
[docs]@cache(level=30)
def db_to_amplitude(S_db: Union[_FloatLike_co, np.ndarray], *, ref: float = 1.0) -> Union[np.floating[Any], np.ndarray]:
"""Convert a dB-scaled spectrogram to an amplitude spectrogram.
This effectively inverts `amplitude_to_db`::
db_to_amplitude(S_db) ~= 10.0**(0.5 * S_db/10 + log10(ref))
Parameters
----------
S_db : np.ndarray
dB-scaled values
ref : number > 0
Optional reference amplitude.
Returns
-------
S : np.ndarray
Linear magnitude values
Notes
-----
This function caches at level 30.
"""
return db_to_power(S_db, ref=ref**2) ** 0.5
[docs]@cache(level=30)
def perceptual_weighting(
S: np.ndarray, frequencies: np.ndarray, *, kind: str = "A", **kwargs: Any
) -> np.ndarray:
"""Perceptual weighting of a power spectrogram::
S_p[..., f, :] = frequency_weighting(f, 'A') + 10*log(S[..., f, :] / ref)
Parameters
----------
S : np.ndarray [shape=(..., d, t)]
Power spectrogram
frequencies : np.ndarray [shape=(d,)]
Center frequency for each row of` `S``
kind : str
The frequency weighting curve to use.
e.g. `'A'`, `'B'`, `'C'`, `'D'`, `None or 'Z'`
**kwargs : additional keyword arguments
Additional keyword arguments to `power_to_db`.
Returns
-------
S_p : np.ndarray [shape=(..., d, t)]
perceptually weighted version of ``S``
See Also
--------
power_to_db
Notes
-----
This function caches at level 30.
Examples
--------
Re-weight a CQT power spectrum, using peak power as reference
>>> y, sr = librosa.load(librosa.ex('trumpet'))
>>> C = np.abs(librosa.cqt(y, sr=sr, fmin=librosa.note_to_hz('A1')))
>>> freqs = librosa.cqt_frequencies(C.shape[0],
... fmin=librosa.note_to_hz('A1'))
>>> perceptual_CQT = librosa.perceptual_weighting(C**2,
... freqs,
... ref=np.max)
>>> perceptual_CQT
array([[ -96.528, -97.101, ..., -108.561, -108.561],
[ -95.88 , -96.479, ..., -107.551, -107.551],
...,
[ -65.142, -53.256, ..., -80.098, -80.098],
[ -71.542, -53.197, ..., -80.311, -80.311]])
>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots(nrows=2, sharex=True, sharey=True)
>>> img = librosa.display.specshow(librosa.amplitude_to_db(C,
... ref=np.max),
... fmin=librosa.note_to_hz('A1'),
... y_axis='cqt_hz', x_axis='time',
... ax=ax[0])
>>> ax[0].set(title='Log CQT power')
>>> ax[0].label_outer()
>>> imgp = librosa.display.specshow(perceptual_CQT, y_axis='cqt_hz',
... fmin=librosa.note_to_hz('A1'),
... x_axis='time', ax=ax[1])
>>> ax[1].set(title='Perceptually weighted log CQT')
>>> fig.colorbar(img, ax=ax[0], format="%+2.0f dB")
>>> fig.colorbar(imgp, ax=ax[1], format="%+2.0f dB")
"""
offset = convert.frequency_weighting(frequencies, kind=kind).reshape((-1, 1))
result: np.ndarray = offset + power_to_db(S, **kwargs)
return result
[docs]@cache(level=30)
def fmt(
y: np.ndarray,
*,
t_min: float = 0.5,
n_fmt: Optional[int] = None,
kind: str = "cubic",
beta: float = 0.5,
over_sample: float = 1,
axis: int = -1,
) -> np.ndarray:
"""Fast Mellin transform (FMT)
The Mellin of a signal `y` is performed by interpolating `y` on an exponential time
axis, applying a polynomial window, and then taking the discrete Fourier transform.
When the Mellin parameter (beta) is 1/2, it is also known as the scale transform. [#]_
The scale transform can be useful for audio analysis because its magnitude is invariant
to scaling of the domain (e.g., time stretching or compression). This is analogous
to the magnitude of the Fourier transform being invariant to shifts in the input domain.
.. [#] De Sena, Antonio, and Davide Rocchesso.
"A fast Mellin and scale transform."
EURASIP Journal on Applied Signal Processing 2007.1 (2007): 75-75.
.. [#] Cohen, L.
"The scale representation."
IEEE Transactions on Signal Processing 41, no. 12 (1993): 3275-3292.
Parameters
----------
y : np.ndarray, real-valued
The input signal(s). Can be multidimensional.
The target axis must contain at least 3 samples.
t_min : float > 0
The minimum time spacing (in samples).
This value should generally be less than 1 to preserve as much information as
possible.
n_fmt : int > 2 or None
The number of scale transform bins to use.
If None, then ``n_bins = over_sample * ceil(n * log((n-1)/t_min))`` is taken,
where ``n = y.shape[axis]``
kind : str
The type of interpolation to use when re-sampling the input.
See `scipy.interpolate.interp1d` for possible values.
Note that the default is to use high-precision (cubic) interpolation.
This can be slow in practice; if speed is preferred over accuracy,
then consider using ``kind='linear'``.
beta : float
The Mellin parameter. ``beta=0.5`` provides the scale transform.
over_sample : float >= 1
Over-sampling factor for exponential resampling.
axis : int
The axis along which to transform ``y``
Returns
-------
x_scale : np.ndarray [dtype=complex]
The scale transform of ``y`` along the ``axis`` dimension.
Raises
------
ParameterError
if ``n_fmt < 2`` or ``t_min <= 0``
or if ``y`` is not finite
or if ``y.shape[axis] < 3``.
Notes
-----
This function caches at level 30.
Examples
--------
>>> # Generate a signal and time-stretch it (with energy normalization)
>>> scale = 1.25
>>> freq = 3.0
>>> x1 = np.linspace(0, 1, num=1024, endpoint=False)
>>> x2 = np.linspace(0, 1, num=int(scale * len(x1)), endpoint=False)
>>> y1 = np.sin(2 * np.pi * freq * x1)
>>> y2 = np.sin(2 * np.pi * freq * x2) / np.sqrt(scale)
>>> # Verify that the two signals have the same energy
>>> np.sum(np.abs(y1)**2), np.sum(np.abs(y2)**2)
(255.99999999999997, 255.99999999999969)
>>> scale1 = librosa.fmt(y1, n_fmt=512)
>>> scale2 = librosa.fmt(y2, n_fmt=512)
>>> # And plot the results
>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots(nrows=2)
>>> ax[0].plot(y1, label='Original')
>>> ax[0].plot(y2, linestyle='--', label='Stretched')
>>> ax[0].set(xlabel='time (samples)', title='Input signals')
>>> ax[0].legend()
>>> ax[1].semilogy(np.abs(scale1), label='Original')
>>> ax[1].semilogy(np.abs(scale2), linestyle='--', label='Stretched')
>>> ax[1].set(xlabel='scale coefficients', title='Scale transform magnitude')
>>> ax[1].legend()
>>> # Plot the scale transform of an onset strength autocorrelation
>>> y, sr = librosa.load(librosa.ex('choice'))
>>> odf = librosa.onset.onset_strength(y=y, sr=sr)
>>> # Auto-correlate with up to 10 seconds lag
>>> odf_ac = librosa.autocorrelate(odf, max_size=10 * sr // 512)
>>> # Normalize
>>> odf_ac = librosa.util.normalize(odf_ac, norm=np.inf)
>>> # Compute the scale transform
>>> odf_ac_scale = librosa.fmt(librosa.util.normalize(odf_ac), n_fmt=512)
>>> # Plot the results
>>> fig, ax = plt.subplots(nrows=3)
>>> ax[0].plot(odf, label='Onset strength')
>>> ax[0].set(xlabel='Time (frames)', title='Onset strength')
>>> ax[1].plot(odf_ac, label='Onset autocorrelation')
>>> ax[1].set(xlabel='Lag (frames)', title='Onset autocorrelation')
>>> ax[2].semilogy(np.abs(odf_ac_scale), label='Scale transform magnitude')
>>> ax[2].set(xlabel='scale coefficients')
"""
n = y.shape[axis]
if n < 3:
raise ParameterError(f"y.shape[{axis}]=={n} < 3")
if t_min <= 0:
raise ParameterError(f"t_min={t_min} must be a positive number")
if n_fmt is None:
if over_sample < 1:
raise ParameterError(f"over_sample={over_sample} must be >= 1")
# The base is the maximum ratio between adjacent samples
# Since the sample spacing is increasing, this is simply the
# ratio between the positions of the last two samples: (n-1)/(n-2)
log_base = np.log(n - 1) - np.log(n - 2)
n_fmt = int(np.ceil(over_sample * (np.log(n - 1) - np.log(t_min)) / log_base))
elif n_fmt < 3:
raise ParameterError(f"n_fmt=={n_fmt} < 3")
else:
log_base = (np.log(n_fmt - 1) - np.log(n_fmt - 2)) / over_sample
if not np.all(np.isfinite(y)):
raise ParameterError("y must be finite everywhere")
base = np.exp(log_base)
# original grid: signal covers [0, 1). This range is arbitrary, but convenient.
# The final sample is positioned at (n-1)/n, so we omit the endpoint
x = np.linspace(0, 1, num=n, endpoint=False)
# build the interpolator
f_interp = scipy.interpolate.interp1d(x, y, kind=kind, axis=axis)
# build the new sampling grid
# exponentially spaced between t_min/n and 1 (exclusive)
# we'll go one past where we need, and drop the last sample
# When over-sampling, the last input sample contributions n_over samples.
# To keep the spacing consistent, we over-sample by n_over, and then
# trim the final samples.
n_over = int(np.ceil(over_sample))
x_exp = np.logspace(
(np.log(t_min) - np.log(n)) / log_base,
0,
num=n_fmt + n_over,
endpoint=False,
base=base,
)[:-n_over]
# Clean up any rounding errors at the boundaries of the interpolation
# The interpolator gets angry if we try to extrapolate, so clipping is necessary here.
if x_exp[0] < t_min or x_exp[-1] > float(n - 1.0) / n:
x_exp = np.clip(x_exp, float(t_min) / n, x[-1])
# Make sure that all sample points are unique
# This should never happen!
if len(np.unique(x_exp)) != len(x_exp):
raise ParameterError("Redundant sample positions in Mellin transform")
# Resample the signal
y_res = f_interp(x_exp)
# Broadcast the window correctly
shape = [1] * y_res.ndim
shape[axis] = -1
# Apply the window and fft
# Normalization is absorbed into the window here for expedience
fft = get_fftlib()
result: np.ndarray = fft.rfft(
y_res * ((x_exp**beta).reshape(shape) * np.sqrt(n) / n_fmt), axis=axis
)
return result
@overload
def pcen(
S: np.ndarray,
*,
sr: float = ...,
hop_length: int = ...,
gain: float = ...,
bias: float = ...,
power: float = ...,
time_constant: float = ...,
eps: float = ...,
b: Optional[float] = ...,
max_size: int = ...,
ref: Optional[np.ndarray] = ...,
axis: int = ...,
max_axis: Optional[int] = ...,
zi: Optional[np.ndarray] = ...,
return_zf: Literal[False] = ...,
) -> np.ndarray:
...
@overload
def pcen(
S: np.ndarray,
*,
sr: float = ...,
hop_length: int = ...,
gain: float = ...,
bias: float = ...,
power: float = ...,
time_constant: float = ...,
eps: float = ...,
b: Optional[float] = ...,
max_size: int = ...,
ref: Optional[np.ndarray] = ...,
axis: int = ...,
max_axis: Optional[int] = ...,
zi: Optional[np.ndarray] = ...,
return_zf: Literal[True],
) -> Tuple[np.ndarray, np.ndarray]:
...
@overload
def pcen(
S: np.ndarray,
*,
sr: float = ...,
hop_length: int = ...,
gain: float = ...,
bias: float = ...,
power: float = ...,
time_constant: float = ...,
eps: float = ...,
b: Optional[float] = ...,
max_size: int = ...,
ref: Optional[np.ndarray] = ...,
axis: int = ...,
max_axis: Optional[int] = ...,
zi: Optional[np.ndarray] = ...,
return_zf: bool = ...,
) -> Union[np.ndarray, Tuple[np.ndarray, np.ndarray]]:
...
[docs]@cache(level=30)
def pcen(
S: np.ndarray,
*,
sr: float = 22050,
hop_length: int = 512,
gain: float = 0.98,
bias: float = 2,
power: float = 0.5,
time_constant: float = 0.400,
eps: float = 1e-6,
b: Optional[float] = None,
max_size: int = 1,
ref: Optional[np.ndarray] = None,
axis: int = -1,
max_axis: Optional[int] = None,
zi: Optional[np.ndarray] = None,
return_zf: bool = False,
) -> Union[np.ndarray, Tuple[np.ndarray, np.ndarray]]:
"""Per-channel energy normalization (PCEN)
This function normalizes a time-frequency representation ``S`` by
performing automatic gain control, followed by nonlinear compression [#]_ ::
P[f, t] = (S / (eps + M[f, t])**gain + bias)**power - bias**power
IMPORTANT: the default values of eps, gain, bias, and power match the
original publication, in which ``S`` is a 40-band mel-frequency
spectrogram with 25 ms windowing, 10 ms frame shift, and raw audio values
in the interval [-2**31; 2**31-1[. If you use these default values, we
recommend to make sure that the raw audio is properly scaled to this
interval, and not to [-1, 1[ as is most often the case.
The matrix ``M`` is the result of applying a low-pass, temporal IIR filter
to ``S``::
M[f, t] = (1 - b) * M[f, t - 1] + b * S[f, t]
If ``b`` is not provided, it is calculated as::
b = (sqrt(1 + 4* T**2) - 1) / (2 * T**2)
where ``T = time_constant * sr / hop_length``. [#]_
This normalization is designed to suppress background noise and
emphasize foreground signals, and can be used as an alternative to
decibel scaling (`amplitude_to_db`).
This implementation also supports smoothing across frequency bins
by specifying ``max_size > 1``. If this option is used, the filtered
spectrogram ``M`` is computed as::
M[f, t] = (1 - b) * M[f, t - 1] + b * R[f, t]
where ``R`` has been max-filtered along the frequency axis, similar to
the SuperFlux algorithm implemented in `onset.onset_strength`::
R[f, t] = max(S[f - max_size//2: f + max_size//2, t])
This can be used to perform automatic gain control on signals that cross
or span multiple frequency bans, which may be desirable for spectrograms
with high frequency resolution.
.. [#] Wang, Y., Getreuer, P., Hughes, T., Lyon, R. F., & Saurous, R. A.
(2017, March). Trainable frontend for robust and far-field keyword spotting.
In Acoustics, Speech and Signal Processing (ICASSP), 2017
IEEE International Conference on (pp. 5670-5674). IEEE.
.. [#] Lostanlen, V., Salamon, J., McFee, B., Cartwright, M., Farnsworth, A.,
Kelling, S., and Bello, J. P. Per-Channel Energy Normalization: Why and How.
IEEE Signal Processing Letters, 26(1), 39-43.
Parameters
----------
S : np.ndarray (non-negative)
The input (magnitude) spectrogram
sr : number > 0 [scalar]
The audio sampling rate
hop_length : int > 0 [scalar]
The hop length of ``S``, expressed in samples
gain : number >= 0 [scalar]
The gain factor. Typical values should be slightly less than 1.
bias : number >= 0 [scalar]
The bias point of the nonlinear compression (default: 2)
power : number >= 0 [scalar]
The compression exponent. Typical values should be between 0 and 0.5.
Smaller values of ``power`` result in stronger compression.
At the limit ``power=0``, polynomial compression becomes logarithmic.
time_constant : number > 0 [scalar]
The time constant for IIR filtering, measured in seconds.
eps : number > 0 [scalar]
A small constant used to ensure numerical stability of the filter.
b : number in [0, 1] [scalar]
The filter coefficient for the low-pass filter.
If not provided, it will be inferred from ``time_constant``.
max_size : int > 0 [scalar]
The width of the max filter applied to the frequency axis.
If left as `1`, no filtering is performed.
ref : None or np.ndarray (shape=S.shape)
An optional pre-computed reference spectrum (``R`` in the above).
If not provided it will be computed from ``S``.
axis : int [scalar]
The (time) axis of the input spectrogram.
max_axis : None or int [scalar]
The frequency axis of the input spectrogram.
If `None`, and ``S`` is two-dimensional, it will be inferred
as the opposite from ``axis``.
If ``S`` is not two-dimensional, and ``max_size > 1``, an error
will be raised.
zi : np.ndarray
The initial filter delay values.
This may be the ``zf`` (final delay values) of a previous call to ``pcen``, or
computed by `scipy.signal.lfilter_zi`.
return_zf : bool
If ``True``, return the final filter delay values along with the PCEN output ``P``.
This is primarily useful in streaming contexts, where the final state of one
block of processing should be used to initialize the next block.
If ``False`` (default) only the PCEN values ``P`` are returned.
Returns
-------
P : np.ndarray, non-negative [shape=(n, m)]
The per-channel energy normalized version of ``S``.
zf : np.ndarray (optional)
The final filter delay values. Only returned if ``return_zf=True``.
See Also
--------
amplitude_to_db
librosa.onset.onset_strength
Examples
--------
Compare PCEN to log amplitude (dB) scaling on Mel spectra
>>> import matplotlib.pyplot as plt
>>> y, sr = librosa.load(librosa.ex('robin'))
>>> # We recommend scaling y to the range [-2**31, 2**31[ before applying
>>> # PCEN's default parameters. Furthermore, we use power=1 to get a
>>> # magnitude spectrum instead of a power spectrum.
>>> S = librosa.feature.melspectrogram(y=y, sr=sr, power=1)
>>> log_S = librosa.amplitude_to_db(S, ref=np.max)
>>> pcen_S = librosa.pcen(S * (2**31))
>>> fig, ax = plt.subplots(nrows=2, sharex=True, sharey=True)
>>> img = librosa.display.specshow(log_S, x_axis='time', y_axis='mel', ax=ax[0])
>>> ax[0].set(title='log amplitude (dB)', xlabel=None)
>>> ax[0].label_outer()
>>> imgpcen = librosa.display.specshow(pcen_S, x_axis='time', y_axis='mel', ax=ax[1])
>>> ax[1].set(title='Per-channel energy normalization')
>>> fig.colorbar(img, ax=ax[0], format="%+2.0f dB")
>>> fig.colorbar(imgpcen, ax=ax[1])
Compare PCEN with and without max-filtering
>>> pcen_max = librosa.pcen(S * (2**31), max_size=3)
>>> fig, ax = plt.subplots(nrows=2, sharex=True, sharey=True)
>>> librosa.display.specshow(pcen_S, x_axis='time', y_axis='mel', ax=ax[0])
>>> ax[0].set(title='Per-channel energy normalization (no max-filter)')
>>> ax[0].label_outer()
>>> img = librosa.display.specshow(pcen_max, x_axis='time', y_axis='mel', ax=ax[1])
>>> ax[1].set(title='Per-channel energy normalization (max_size=3)')
>>> fig.colorbar(img, ax=ax)
"""
if power < 0:
raise ParameterError(f"power={power} must be nonnegative")
if gain < 0:
raise ParameterError(f"gain={gain} must be non-negative")
if bias < 0:
raise ParameterError(f"bias={bias} must be non-negative")
if eps <= 0:
raise ParameterError(f"eps={eps} must be strictly positive")
if time_constant <= 0:
raise ParameterError(f"time_constant={time_constant} must be strictly positive")
if not util.is_positive_int(max_size):
raise ParameterError(f"max_size={max_size} must be a positive integer")
if b is None:
t_frames = time_constant * sr / float(hop_length)
# By default, this solves the equation for b:
# b**2 + (1 - b) / t_frames - 2 = 0
# which approximates the full-width half-max of the
# squared frequency response of the IIR low-pass filter
b = (np.sqrt(1 + 4 * t_frames**2) - 1) / (2 * t_frames**2)
if not 0 <= b <= 1:
raise ParameterError(f"b={b} must be between 0 and 1")
if np.issubdtype(S.dtype, np.complexfloating):
warnings.warn(
"pcen was called on complex input so phase "
"information will be discarded. To suppress this warning, "
"call pcen(np.abs(D)) instead.",
stacklevel=2,
)
S = np.abs(S)
if ref is None:
if max_size == 1:
ref = S
elif S.ndim == 1:
raise ParameterError(
"Max-filtering cannot be applied to 1-dimensional input"
)
else:
if max_axis is None:
if S.ndim != 2:
raise ParameterError(
f"Max-filtering a {S.ndim:d}-dimensional spectrogram "
"requires you to specify max_axis"
)
# if axis = 0, max_axis=1
# if axis = +- 1, max_axis = 0
max_axis = np.mod(1 - axis, 2)
ref = scipy.ndimage.maximum_filter1d(S, max_size, axis=max_axis)
if zi is None:
# Make sure zi matches dimension to input
shape = tuple([1] * ref.ndim)
zi = np.empty(shape)
zi[:] = scipy.signal.lfilter_zi([b], [1, b - 1])[:]
# Temporal integration
S_smooth: np.ndarray
zf: np.ndarray
S_smooth, zf = scipy.signal.lfilter([b], [1, b - 1], ref, zi=zi, axis=axis)
# Adaptive gain control
# Working in log-space gives us some stability, and a slight speedup
smooth = np.exp(-gain * (np.log(eps) + np.log1p(S_smooth / eps)))
# Dynamic range compression
S_out: np.ndarray
if power == 0:
S_out = np.log1p(S * smooth)
elif bias == 0:
S_out = np.exp(power * (np.log(S) + np.log(smooth)))
else:
S_out = (bias**power) * np.expm1(power * np.log1p(S * smooth / bias))
if return_zf:
return S_out, zf
else:
return S_out
[docs]def griffinlim(
S: np.ndarray,
*,
n_iter: int = 32,
hop_length: Optional[int] = None,
win_length: Optional[int] = None,
n_fft: Optional[int] = None,
window: _WindowSpec = "hann",
center: bool = True,
dtype: Optional[DTypeLike] = None,
length: Optional[int] = None,
pad_mode: _PadModeSTFT = "constant",
momentum: float = 0.99,
init: Optional[str] = "random",
random_state: Optional[
Union[int, np.random.RandomState, np.random.Generator]
] = None,
) -> np.ndarray:
"""Approximate magnitude spectrogram inversion using the "fast" Griffin-Lim algorithm.
Given a short-time Fourier transform magnitude matrix (``S``), the algorithm randomly
initializes phase estimates, and then alternates forward- and inverse-STFT
operations. [#]_
Note that this assumes reconstruction of a real-valued time-domain signal, and
that ``S`` contains only the non-negative frequencies (as computed by
`stft`).
The "fast" GL method [#]_ uses a momentum parameter to accelerate convergence.
.. [#] D. W. Griffin and J. S. Lim,
"Signal estimation from modified short-time Fourier transform,"
IEEE Trans. ASSP, vol.32, no.2, pp.236–243, Apr. 1984.
.. [#] Perraudin, N., Balazs, P., & Søndergaard, P. L.
"A fast Griffin-Lim algorithm,"
IEEE Workshop on Applications of Signal Processing to Audio and Acoustics (pp. 1-4),
Oct. 2013.
Parameters
----------
S : np.ndarray [shape=(..., n_fft // 2 + 1, t), non-negative]
An array of short-time Fourier transform magnitudes as produced by
`stft`.
n_iter : int > 0
The number of iterations to run
hop_length : None or int > 0
The hop length of the STFT. If not provided, it will default to ``n_fft // 4``
win_length : None or int > 0
The window length of the STFT. By default, it will equal ``n_fft``
n_fft : None or int > 0
The number of samples per frame.
By default, this will be inferred from the shape of ``S`` as an even number.
However, if an odd frame length was used, you can explicitly set ``n_fft``.
window : string, tuple, number, function, or np.ndarray [shape=(n_fft,)]
A window specification as supported by `stft` or `istft`
center : boolean
If ``True``, the STFT is assumed to use centered frames.
If ``False``, the STFT is assumed to use left-aligned frames.
dtype : np.dtype
Real numeric type for the time-domain signal. Default is inferred
to match the precision of the input spectrogram.
length : None or int > 0
If provided, the output ``y`` is zero-padded or clipped to exactly ``length``
samples.
pad_mode : string
If ``center=True``, the padding mode to use at the edges of the signal.
By default, STFT uses zero padding.
momentum : number >= 0
The momentum parameter for fast Griffin-Lim.
Setting this to 0 recovers the original Griffin-Lim method [1]_.
Values near 1 can lead to faster convergence, but above 1 may not converge.
init : None or 'random' [default]
If 'random' (the default), then phase values are initialized randomly
according to ``random_state``. This is recommended when the input ``S`` is
a magnitude spectrogram with no initial phase estimates.
If `None`, then the phase is initialized from ``S``. This is useful when
an initial guess for phase can be provided, or when you want to resume
Griffin-Lim from a previous output.
random_state : None, int, np.random.RandomState, or np.random.Generator
If int, random_state is the seed used by the random number generator
for phase initialization.
If `np.random.RandomState` or `np.random.Generator` instance, the random number
generator itself.
If `None`, defaults to the `np.random.default_rng()` object.
Returns
-------
y : np.ndarray [shape=(..., n)]
time-domain signal reconstructed from ``S``
See Also
--------
stft
istft
magphase
filters.get_window
Examples
--------
A basic STFT inverse example
>>> y, sr = librosa.load(librosa.ex('trumpet'))
>>> # Get the magnitude spectrogram
>>> S = np.abs(librosa.stft(y))
>>> # Invert using Griffin-Lim
>>> y_inv = librosa.griffinlim(S)
>>> # Invert without estimating phase
>>> y_istft = librosa.istft(S)
Wave-plot the results
>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots(nrows=3, sharex=True, sharey=True)
>>> librosa.display.waveshow(y, sr=sr, color='b', ax=ax[0])
>>> ax[0].set(title='Original', xlabel=None)
>>> ax[0].label_outer()
>>> librosa.display.waveshow(y_inv, sr=sr, color='g', ax=ax[1])
>>> ax[1].set(title='Griffin-Lim reconstruction', xlabel=None)
>>> ax[1].label_outer()
>>> librosa.display.waveshow(y_istft, sr=sr, color='r', ax=ax[2])
>>> ax[2].set_title('Magnitude-only istft reconstruction')
"""
if random_state is None:
rng = np.random.default_rng()
elif isinstance(random_state, int):
rng = np.random.RandomState(seed=random_state) # type: ignore
elif isinstance(random_state, (np.random.RandomState, np.random.Generator)):
rng = random_state # type: ignore
else:
raise ParameterError(f"Unsupported random_state={random_state!r}")
if momentum > 1:
warnings.warn(
f"Griffin-Lim with momentum={momentum} > 1 can be unstable. "
"Proceed with caution!",
stacklevel=2,
)
elif momentum < 0:
raise ParameterError(f"griffinlim() called with momentum={momentum} < 0")
# Infer n_fft from the spectrogram shape
if n_fft is None:
n_fft = 2 * (S.shape[-2] - 1)
# Infer the dtype from S
angles = np.empty(S.shape, dtype=util.dtype_r2c(S.dtype))
eps = util.tiny(angles)
if init == "random":
# randomly initialize the phase
angles[:] = util.phasor((2 * np.pi * rng.random(size=S.shape)))
elif init is None:
# Initialize an all ones complex matrix
angles[:] = 1.0
else:
raise ParameterError(f"init={init} must either None or 'random'")
# Place-holders for temporary data and reconstructed buffer
rebuilt = None
tprev = None
inverse = None
# Absorb magnitudes into angles
angles *= S
for _ in range(n_iter):
# Invert with our current estimate of the phases
inverse = istft(
angles,
hop_length=hop_length,
win_length=win_length,
n_fft=n_fft,
window=window,
center=center,
dtype=dtype,
length=length,
out=inverse,
)
# Rebuild the spectrogram
rebuilt = stft(
inverse,
n_fft=n_fft,
hop_length=hop_length,
win_length=win_length,
window=window,
center=center,
pad_mode=pad_mode,
out=rebuilt,
)
# Update our phase estimates
angles[:] = rebuilt
if tprev is not None:
angles -= (momentum / (1 + momentum)) * tprev
angles /= np.abs(angles) + eps
angles *= S
# Store
rebuilt, tprev = tprev, rebuilt
# Return the final phase estimates
return istft(
angles,
hop_length=hop_length,
win_length=win_length,
n_fft=n_fft,
window=window,
center=center,
dtype=dtype,
length=length,
out=inverse,
)
def _spectrogram(
*,
y: Optional[np.ndarray] = None,
S: Optional[np.ndarray] = None,
n_fft: Optional[int] = 2048,
hop_length: Optional[int] = 512,
power: float = 1,
win_length: Optional[int] = None,
window: _WindowSpec = "hann",
center: bool = True,
pad_mode: _PadModeSTFT = "constant",
) -> Tuple[np.ndarray, int]:
"""Retrieve a magnitude spectrogram.
This is primarily used in feature extraction functions that can operate on
either audio time-series or spectrogram input.
Parameters
----------
y : None or np.ndarray
If provided, an audio time series
S : None or np.ndarray
Spectrogram input, optional
n_fft : int > 0
STFT window size
hop_length : int > 0
STFT hop length
power : float > 0
Exponent for the magnitude spectrogram,
e.g., 1 for energy, 2 for power, etc.
win_length : int <= n_fft [scalar]
Each frame of audio is windowed by ``window``.
The window will be of length ``win_length`` and then padded
with zeros to match ``n_fft``.
If unspecified, defaults to ``win_length = n_fft``.
window : string, tuple, number, function, or np.ndarray [shape=(n_fft,)]
- a window specification (string, tuple, or number);
see `scipy.signal.get_window`
- a window function, such as `scipy.signal.windows.hann`
- a vector or array of length ``n_fft``
.. see also:: `filters.get_window`
center : boolean
- If ``True``, the signal ``y`` is padded so that frame
``t`` is centered at ``y[t * hop_length]``.
- If ``False``, then frame ``t`` begins at ``y[t * hop_length]``
pad_mode : string
If ``center=True``, the padding mode to use at the edges of the signal.
By default, STFT uses zero padding.
Returns
-------
S_out : np.ndarray [dtype=np.float]
- If ``S`` is provided as input, then ``S_out == S``
- Else, ``S_out = |stft(y, ...)|**power``
n_fft : int > 0
- If ``S`` is provided, then ``n_fft`` is inferred from ``S``
- Else, copied from input
"""
if S is not None:
# Infer n_fft from spectrogram shape, but only if it mismatches
if n_fft is None or n_fft // 2 + 1 != S.shape[-2]:
n_fft = 2 * (S.shape[-2] - 1)
else:
# Otherwise, compute a magnitude spectrogram from input
if n_fft is None:
raise ParameterError(f"Unable to compute spectrogram with n_fft={n_fft}")
if y is None:
raise ParameterError(
"Input signal must be provided to compute a spectrogram"
)
S = (
np.abs(
stft(
y,
n_fft=n_fft,
hop_length=hop_length,
win_length=win_length,
center=center,
window=window,
pad_mode=pad_mode,
)
)
** power
)
return S, n_fft