#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Filters
=======
Filter bank construction
------------------------
.. autosummary::
:toctree: generated/
mel
chroma
wavelet
semitone_filterbank
Window functions
----------------
.. autosummary::
:toctree: generated/
window_bandwidth
get_window
Miscellaneous
-------------
.. autosummary::
:toctree: generated/
wavelet_lengths
cq_to_chroma
mr_frequencies
window_sumsquare
diagonal_filter
Deprecated
----------
.. autosummary::
:toctree: generated/
constant_q
constant_q_lengths
"""
import warnings
import numpy as np
import scipy
import scipy.signal
import scipy.ndimage
from numba import jit
from ._cache import cache
from . import util
from .util.exceptions import ParameterError
from .util.decorators import deprecated
from .core.convert import note_to_hz, hz_to_midi, midi_to_hz, hz_to_octs
from .core.convert import fft_frequencies, mel_frequencies
from numpy.typing import ArrayLike, DTypeLike
from typing import Any, List, Optional, Tuple, Union
from typing_extensions import Literal
from ._typing import _WindowSpec, _FloatLike_co
__all__ = [
"mel",
"chroma",
"constant_q",
"constant_q_lengths",
"cq_to_chroma",
"window_bandwidth",
"get_window",
"mr_frequencies",
"semitone_filterbank",
"window_sumsquare",
"diagonal_filter",
"wavelet",
"wavelet_lengths",
]
# Dictionary of window function bandwidths
WINDOW_BANDWIDTHS = {
"bart": 1.3334961334912805,
"barthann": 1.4560255965133932,
"bartlett": 1.3334961334912805,
"bkh": 2.0045975283585014,
"black": 1.7269681554262326,
"blackharr": 2.0045975283585014,
"blackman": 1.7269681554262326,
"blackmanharris": 2.0045975283585014,
"blk": 1.7269681554262326,
"bman": 1.7859588613860062,
"bmn": 1.7859588613860062,
"bohman": 1.7859588613860062,
"box": 1.0,
"boxcar": 1.0,
"brt": 1.3334961334912805,
"brthan": 1.4560255965133932,
"bth": 1.4560255965133932,
"cosine": 1.2337005350199792,
"flat": 2.7762255046484143,
"flattop": 2.7762255046484143,
"flt": 2.7762255046484143,
"halfcosine": 1.2337005350199792,
"ham": 1.3629455320350348,
"hamm": 1.3629455320350348,
"hamming": 1.3629455320350348,
"han": 1.50018310546875,
"hann": 1.50018310546875,
"nut": 1.9763500280946082,
"nutl": 1.9763500280946082,
"nuttall": 1.9763500280946082,
"ones": 1.0,
"par": 1.9174603174603191,
"parz": 1.9174603174603191,
"parzen": 1.9174603174603191,
"rect": 1.0,
"rectangular": 1.0,
"tri": 1.3331706523555851,
"triang": 1.3331706523555851,
"triangle": 1.3331706523555851,
}
[docs]@cache(level=10)
def mel(
*,
sr: float,
n_fft: int,
n_mels: int = 128,
fmin: float = 0.0,
fmax: Optional[float] = None,
htk: bool = False,
norm: Optional[Union[Literal["slaney"], float]] = "slaney",
dtype: DTypeLike = np.float32,
) -> np.ndarray:
"""Create a Mel filter-bank.
This produces a linear transformation matrix to project
FFT bins onto Mel-frequency bins.
Parameters
----------
sr : number > 0 [scalar]
sampling rate of the incoming signal
n_fft : int > 0 [scalar]
number of FFT components
n_mels : int > 0 [scalar]
number of Mel bands to generate
fmin : float >= 0 [scalar]
lowest frequency (in Hz)
fmax : float >= 0 [scalar]
highest frequency (in Hz).
If `None`, use ``fmax = sr / 2.0``
htk : bool [scalar]
use HTK formula instead of Slaney
norm : {None, 'slaney', or number} [scalar]
If 'slaney', divide the triangular mel weights by the width of the mel band
(area normalization).
If numeric, use `librosa.util.normalize` to normalize each filter by to unit l_p norm.
See `librosa.util.normalize` for a full description of supported norm values
(including `+-np.inf`).
Otherwise, leave all the triangles aiming for a peak value of 1.0
dtype : np.dtype
The data type of the output basis.
By default, uses 32-bit (single-precision) floating point.
Returns
-------
M : np.ndarray [shape=(n_mels, 1 + n_fft/2)]
Mel transform matrix
See Also
--------
librosa.util.normalize
Notes
-----
This function caches at level 10.
Examples
--------
>>> melfb = librosa.filters.mel(sr=22050, n_fft=2048)
>>> melfb
array([[ 0. , 0.016, ..., 0. , 0. ],
[ 0. , 0. , ..., 0. , 0. ],
...,
[ 0. , 0. , ..., 0. , 0. ],
[ 0. , 0. , ..., 0. , 0. ]])
Clip the maximum frequency to 8KHz
>>> librosa.filters.mel(sr=22050, n_fft=2048, fmax=8000)
array([[ 0. , 0.02, ..., 0. , 0. ],
[ 0. , 0. , ..., 0. , 0. ],
...,
[ 0. , 0. , ..., 0. , 0. ],
[ 0. , 0. , ..., 0. , 0. ]])
>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots()
>>> img = librosa.display.specshow(melfb, x_axis='linear', ax=ax)
>>> ax.set(ylabel='Mel filter', title='Mel filter bank')
>>> fig.colorbar(img, ax=ax)
"""
if fmax is None:
fmax = float(sr) / 2
# Initialize the weights
n_mels = int(n_mels)
weights = np.zeros((n_mels, int(1 + n_fft // 2)), dtype=dtype)
# Center freqs of each FFT bin
fftfreqs = fft_frequencies(sr=sr, n_fft=n_fft)
# 'Center freqs' of mel bands - uniformly spaced between limits
mel_f = mel_frequencies(n_mels + 2, fmin=fmin, fmax=fmax, htk=htk)
fdiff = np.diff(mel_f)
ramps = np.subtract.outer(mel_f, fftfreqs)
for i in range(n_mels):
# lower and upper slopes for all bins
lower = -ramps[i] / fdiff[i]
upper = ramps[i + 2] / fdiff[i + 1]
# .. then intersect them with each other and zero
weights[i] = np.maximum(0, np.minimum(lower, upper))
if isinstance(norm, str):
if norm == "slaney":
# Slaney-style mel is scaled to be approx constant energy per channel
enorm = 2.0 / (mel_f[2 : n_mels + 2] - mel_f[:n_mels])
weights *= enorm[:, np.newaxis]
else:
raise ParameterError(f"Unsupported norm={norm}")
else:
weights = util.normalize(weights, norm=norm, axis=-1)
# Only check weights if f_mel[0] is positive
if not np.all((mel_f[:-2] == 0) | (weights.max(axis=1) > 0)):
# This means we have an empty channel somewhere
warnings.warn(
"Empty filters detected in mel frequency basis. "
"Some channels will produce empty responses. "
"Try increasing your sampling rate (and fmax) or "
"reducing n_mels.",
stacklevel=2,
)
return weights
[docs]@cache(level=10)
def chroma(
*,
sr: float,
n_fft: int,
n_chroma: int = 12,
tuning: float = 0.0,
ctroct: float = 5.0,
octwidth: Union[float, None] = 2,
norm: Optional[float] = 2,
base_c: bool = True,
dtype: DTypeLike = np.float32,
) -> np.ndarray:
"""Create a chroma filter bank.
This creates a linear transformation matrix to project
FFT bins onto chroma bins (i.e. pitch classes).
Parameters
----------
sr : number > 0 [scalar]
audio sampling rate
n_fft : int > 0 [scalar]
number of FFT bins
n_chroma : int > 0 [scalar]
number of chroma bins
tuning : float
Tuning deviation from A440 in fractions of a chroma bin.
ctroct : float > 0 [scalar]
octwidth : float > 0 or None [scalar]
``ctroct`` and ``octwidth`` specify a dominance window:
a Gaussian weighting centered on ``ctroct`` (in octs, A0 = 27.5Hz)
and with a gaussian half-width of ``octwidth``.
Set ``octwidth`` to `None` to use a flat weighting.
norm : float > 0 or np.inf
Normalization factor for each filter
base_c : bool
If True, the filter bank will start at 'C'.
If False, the filter bank will start at 'A'.
dtype : np.dtype
The data type of the output basis.
By default, uses 32-bit (single-precision) floating point.
Returns
-------
wts : ndarray [shape=(n_chroma, 1 + n_fft / 2)]
Chroma filter matrix
See Also
--------
librosa.util.normalize
librosa.feature.chroma_stft
Notes
-----
This function caches at level 10.
Examples
--------
Build a simple chroma filter bank
>>> chromafb = librosa.filters.chroma(sr=22050, n_fft=4096)
array([[ 1.689e-05, 3.024e-04, ..., 4.639e-17, 5.327e-17],
[ 1.716e-05, 2.652e-04, ..., 2.674e-25, 3.176e-25],
...,
[ 1.578e-05, 3.619e-04, ..., 8.577e-06, 9.205e-06],
[ 1.643e-05, 3.355e-04, ..., 1.474e-10, 1.636e-10]])
Use quarter-tones instead of semitones
>>> librosa.filters.chroma(sr=22050, n_fft=4096, n_chroma=24)
array([[ 1.194e-05, 2.138e-04, ..., 6.297e-64, 1.115e-63],
[ 1.206e-05, 2.009e-04, ..., 1.546e-79, 2.929e-79],
...,
[ 1.162e-05, 2.372e-04, ..., 6.417e-38, 9.923e-38],
[ 1.180e-05, 2.260e-04, ..., 4.697e-50, 7.772e-50]])
Equally weight all octaves
>>> librosa.filters.chroma(sr=22050, n_fft=4096, octwidth=None)
array([[ 3.036e-01, 2.604e-01, ..., 2.445e-16, 2.809e-16],
[ 3.084e-01, 2.283e-01, ..., 1.409e-24, 1.675e-24],
...,
[ 2.836e-01, 3.116e-01, ..., 4.520e-05, 4.854e-05],
[ 2.953e-01, 2.888e-01, ..., 7.768e-10, 8.629e-10]])
>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots()
>>> img = librosa.display.specshow(chromafb, x_axis='linear', ax=ax)
>>> ax.set(ylabel='Chroma filter', title='Chroma filter bank')
>>> fig.colorbar(img, ax=ax)
"""
wts = np.zeros((n_chroma, n_fft))
# Get the FFT bins, not counting the DC component
frequencies = np.linspace(0, sr, n_fft, endpoint=False)[1:]
frqbins = n_chroma * hz_to_octs(
frequencies, tuning=tuning, bins_per_octave=n_chroma
)
# make up a value for the 0 Hz bin = 1.5 octaves below bin 1
# (so chroma is 50% rotated from bin 1, and bin width is broad)
frqbins = np.concatenate(([frqbins[0] - 1.5 * n_chroma], frqbins))
binwidthbins = np.concatenate((np.maximum(frqbins[1:] - frqbins[:-1], 1.0), [1]))
D = np.subtract.outer(frqbins, np.arange(0, n_chroma, dtype="d")).T
n_chroma2 = np.round(float(n_chroma) / 2)
# Project into range -n_chroma/2 .. n_chroma/2
# add on fixed offset of 10*n_chroma to ensure all values passed to
# rem are positive
D = np.remainder(D + n_chroma2 + 10 * n_chroma, n_chroma) - n_chroma2
# Gaussian bumps - 2*D to make them narrower
wts = np.exp(-0.5 * (2 * D / np.tile(binwidthbins, (n_chroma, 1))) ** 2)
# normalize each column
wts = util.normalize(wts, norm=norm, axis=0)
# Maybe apply scaling for fft bins
if octwidth is not None:
wts *= np.tile(
np.exp(-0.5 * (((frqbins / n_chroma - ctroct) / octwidth) ** 2)),
(n_chroma, 1),
)
if base_c:
wts = np.roll(wts, -3 * (n_chroma // 12), axis=0)
# remove aliasing columns, copy to ensure row-contiguity
return np.ascontiguousarray(wts[:, : int(1 + n_fft / 2)], dtype=dtype)
def __float_window(window_spec):
"""Decorate a window function to support fractional input lengths.
This function guarantees that for fractional ``x``, the following hold:
1. ``__float_window(window_function)(x)`` has length ``np.ceil(x)``
2. all values from ``np.floor(x)`` are set to 0.
For integer-valued ``x``, there should be no change in behavior.
"""
def _wrap(n, *args, **kwargs):
"""Wrap the window"""
n_min, n_max = int(np.floor(n)), int(np.ceil(n))
window = get_window(window_spec, n_min)
if len(window) < n_max:
window = np.pad(window, [(0, n_max - len(window))], mode="constant")
window[n_min:] = 0.0
return window
return _wrap
[docs]@deprecated(version="0.9.0", version_removed="1.0")
def constant_q(
*,
sr: float,
fmin: Optional[_FloatLike_co] = None,
n_bins: int = 84,
bins_per_octave: int = 12,
window: _WindowSpec = "hann",
filter_scale: float = 1,
pad_fft: bool = True,
norm: Optional[float] = 1,
dtype: DTypeLike = np.complex64,
gamma: float = 0,
**kwargs: Any,
) -> Tuple[np.ndarray, np.ndarray]:
r"""Construct a constant-Q basis.
This function constructs a filter bank similar to Morlet wavelets,
where complex exponentials are windowed to different lengths
such that the number of cycles remains fixed for all frequencies.
By default, a Hann window (rather than the Gaussian window of Morlet wavelets)
is used, but this can be controlled by the ``window`` parameter.
Frequencies are spaced geometrically, increasing by a factor of
``(2**(1./bins_per_octave))`` at each successive band.
.. warning:: This function is deprecated as of v0.9 and will be removed in 1.0.
See `librosa.filters.wavelet`.
Parameters
----------
sr : number > 0 [scalar]
Audio sampling rate
fmin : float > 0 [scalar]
Minimum frequency bin. Defaults to `C1 ~= 32.70`
n_bins : int > 0 [scalar]
Number of frequencies. Defaults to 7 octaves (84 bins).
bins_per_octave : int > 0 [scalar]
Number of bins per octave
window : string, tuple, number, or function
Windowing function to apply to filters.
filter_scale : float > 0 [scalar]
Scale of filter windows.
Small values (<1) use shorter windows for higher temporal resolution.
pad_fft : boolean
Center-pad all filters up to the nearest integral power of 2.
By default, padding is done with zeros, but this can be overridden
by setting the ``mode=`` field in *kwargs*.
norm : {inf, -inf, 0, float > 0}
Type of norm to use for basis function normalization.
See librosa.util.normalize
gamma : number >= 0
Bandwidth offset for variable-Q transforms.
``gamma=0`` produces a constant-Q filterbank.
dtype : np.dtype
The data type of the output basis.
By default, uses 64-bit (single precision) complex floating point.
**kwargs : additional keyword arguments
Arguments to `np.pad()` when ``pad==True``.
Returns
-------
filters : np.ndarray, ``len(filters) == n_bins``
``filters[i]`` is ``i``\ th time-domain CQT basis filter
lengths : np.ndarray, ``len(lengths) == n_bins``
The (fractional) length of each filter
Notes
-----
This function caches at level 10.
See Also
--------
wavelet
constant_q_lengths
librosa.cqt
librosa.vqt
librosa.util.normalize
Examples
--------
Use a shorter window for each filter
>>> basis, lengths = librosa.filters.constant_q(sr=22050, filter_scale=0.5)
Plot one octave of filters in time and frequency
>>> import matplotlib.pyplot as plt
>>> basis, lengths = librosa.filters.constant_q(sr=22050)
>>> fig, ax = plt.subplots(nrows=2, figsize=(10, 6))
>>> notes = librosa.midi_to_note(np.arange(24, 24 + len(basis)))
>>> for i, (f, n) in enumerate(zip(basis, notes[:12])):
... f_scale = librosa.util.normalize(f) / 2
... ax[0].plot(i + f_scale.real)
... ax[0].plot(i + f_scale.imag, linestyle=':')
>>> ax[0].set(yticks=np.arange(len(notes[:12])), yticklabels=notes[:12],
... ylabel='CQ filters',
... title='CQ filters (one octave, time domain)',
... xlabel='Time (samples at 22050 Hz)')
>>> ax[0].legend(['Real', 'Imaginary'])
>>> F = np.abs(np.fft.fftn(basis, axes=[-1]))
>>> # Keep only the positive frequencies
>>> F = F[:, :(1 + F.shape[1] // 2)]
>>> librosa.display.specshow(F, x_axis='linear', y_axis='cqt_note', ax=ax[1])
>>> ax[1].set(ylabel='CQ filters', title='CQ filter magnitudes (frequency domain)')
"""
if fmin is None:
fmin = note_to_hz("C1")
# Pass-through parameters to get the filter lengths
lengths = constant_q_lengths(
sr=sr,
fmin=fmin,
n_bins=n_bins,
bins_per_octave=bins_per_octave,
window=window,
filter_scale=filter_scale,
gamma=gamma,
)
freqs = fmin * (2.0 ** (np.arange(n_bins, dtype=float) / bins_per_octave))
# Build the filters
filters = []
for ilen, freq in zip(lengths, freqs):
# Build the filter: note, length will be ceil(ilen)
sig = util.phasor(
np.arange(-ilen // 2, ilen // 2, dtype=float) * 2 * np.pi * freq / sr
)
# Apply the windowing function
sig = sig * __float_window(window)(len(sig))
# Normalize
sig = util.normalize(sig, norm=norm)
filters.append(sig)
# Pad and stack
max_len = max(lengths)
if pad_fft:
max_len = int(2.0 ** (np.ceil(np.log2(max_len))))
else:
max_len = int(np.ceil(max_len))
filters = np.asarray(
[util.pad_center(filt, size=max_len, **kwargs) for filt in filters], dtype=dtype
)
return filters, np.asarray(lengths)
[docs]@deprecated(version="0.9.0", version_removed="1.0")
@cache(level=10)
def constant_q_lengths(
*,
sr: float,
fmin: _FloatLike_co,
n_bins: int = 84,
bins_per_octave: int = 12,
window: _WindowSpec = "hann",
filter_scale: float = 1,
gamma: float = 0,
) -> np.ndarray:
r"""Return length of each filter in a constant-Q basis.
.. warning:: This function is deprecated as of v0.9 and will be removed in 1.0.
See `librosa.filters.wavelet_lengths`.
Parameters
----------
sr : number > 0 [scalar]
Audio sampling rate
fmin : float > 0 [scalar]
Minimum frequency bin.
n_bins : int > 0 [scalar]
Number of frequencies. Defaults to 7 octaves (84 bins).
bins_per_octave : int > 0 [scalar]
Number of bins per octave
window : str or callable
Window function to use on filters
filter_scale : float > 0 [scalar]
Resolution of filter windows. Larger values use longer windows.
gamma : number >= 0
Bandwidth offset for variable-Q transforms.
``gamma=0`` produces a constant-Q filterbank.
Returns
-------
lengths : np.ndarray
The length of each filter.
Notes
-----
This function caches at level 10.
See Also
--------
wavelet_lengths
"""
if fmin <= 0:
raise ParameterError("fmin must be strictly positive")
if bins_per_octave <= 0:
raise ParameterError("bins_per_octave must be positive")
if filter_scale <= 0:
raise ParameterError("filter_scale must be positive")
if n_bins <= 0 or not isinstance(n_bins, (int, np.integer)):
raise ParameterError("n_bins must be a positive integer")
# Compute the frequencies
freq = fmin * (2.0 ** (np.arange(n_bins, dtype=float) / bins_per_octave))
# Q should be capitalized here, so we suppress the name warning
# pylint: disable=invalid-name
#
# Balance filter bandwidths
alpha = (2.0 ** (2 / bins_per_octave) - 1) / (2.0 ** (2 / bins_per_octave) + 1)
Q = float(filter_scale) / alpha
if max(freq * (1 + 0.5 * window_bandwidth(window) / Q)) > sr / 2.0:
raise ParameterError(
f"Maximum filter frequency={max(freq):.2f} would exceed Nyquist={sr/2}"
)
# Convert frequencies to filter lengths
lengths: np.ndarray = Q * sr / (freq + gamma / alpha)
return lengths
[docs]@cache(level=10)
def wavelet_lengths(
*,
freqs: ArrayLike,
sr: float = 22050,
window: _WindowSpec = "hann",
filter_scale: float = 1,
gamma: Optional[float] = 0,
alpha: Optional[Union[float, np.ndarray]] = None,
) -> Tuple[np.ndarray, float]:
"""Return length of each filter in a wavelet basis.
Parameters
----------
freqs : np.ndarray (positive)
Center frequencies of the filters (in Hz).
Must be in ascending order.
sr : number > 0 [scalar]
Audio sampling rate
window : str or callable
Window function to use on filters
filter_scale : float > 0 [scalar]
Resolution of filter windows. Larger values use longer windows.
gamma : number >= 0 [scalar, optional]
Bandwidth offset for determining filter lengths, as used in
Variable-Q transforms.
Bandwidth for the k'th filter is determined by::
B[k] = alpha[k] * freqs[k] + gamma
``alpha[k]`` is twice the relative difference between ``freqs[k+1]`` and ``freqs[k-1]``::
alpha[k] = (freqs[k+1]-freqs[k-1]) / (freqs[k+1]+freqs[k-1])
If ``freqs`` follows a geometric progression (as in CQT and VQT), the vector
``alpha`` is constant and such that::
(1 + alpha) * freqs[k-1] = (1 - alpha) * freqs[k+1]
Furthermore, if ``gamma=0`` (default), ``alpha`` is such that even-``k`` and
odd-``k`` filters are interleaved::
freqs[k-1] + B[k-1] = freqs[k+1] - B[k+1]
If ``gamma=None`` is specified, then ``gamma`` is computed such
that each filter has bandwidth proportional to the equivalent
rectangular bandwidth (ERB) at frequency ``freqs[k]``::
gamma[k] = 24.7 * alpha[k] / 0.108
as derived by [#]_.
.. [#] Glasberg, Brian R., and Brian CJ Moore.
"Derivation of auditory filter shapes from notched-noise data."
Hearing research 47.1-2 (1990): 103-138.
alpha : number > 0 [optional]
Optional pre-computed relative bandwidth parameter.
Note that this must be provided if ``len(freqs)==1`` because bandwidth cannot be
inferred from a single frequency.
Otherwise, if left unspecified, it will be automatically derived by the rules
specified above.
Returns
-------
lengths : np.ndarray
The length of each filter.
f_cutoff : float
The lowest frequency at which all filters' main lobes have decayed by
at least 3dB.
This second output serves in cqt and vqt to ensure that all wavelet
bands remain below the Nyquist frequency.
Notes
-----
This function caches at level 10.
Raises
------
ParameterError
- If ``filter_scale`` is not strictly positive
- If ``gamma`` is a negative number
- If any frequencies are <= 0
- If the frequency array is not sorted in ascending order
"""
freqs = np.asarray(freqs)
if filter_scale <= 0:
raise ParameterError(f"filter_scale={filter_scale} must be positive")
if gamma is not None and gamma < 0:
raise ParameterError(f"gamma={gamma} must be non-negative")
if np.any(freqs <= 0):
raise ParameterError("frequencies must be strictly positive")
if len(freqs) > 1 and np.any(freqs[:-1] > freqs[1:]):
raise ParameterError(
f"Frequency array={freqs} must be in strictly ascending order"
)
if alpha is None:
alpha = _relative_bandwidth(freqs=freqs)
else:
alpha = np.asarray(alpha)
gamma_: Union[_FloatLike_co, np.ndarray]
if gamma is None:
gamma_ = alpha * 24.7 / 0.108
else:
gamma_ = gamma
# Q should be capitalized here, so we suppress the name warning
# pylint: disable=invalid-name
Q = float(filter_scale) / alpha
# How far up does our highest frequency reach?
f_cutoff = max(freqs * (1 + 0.5 * window_bandwidth(window) / Q) + 0.5 * gamma_)
# Convert frequencies to filter lengths
lengths = Q * sr / (freqs + gamma_ / alpha)
return lengths, f_cutoff
def _relative_bandwidth(*, freqs: np.ndarray) -> np.ndarray:
"""Compute the relative bandwidth for each of a set of specified frequencies.
This function is used as a helper in wavelet basis construction.
Parameters
----------
freqs : np.ndarray
The array of frequencies
Returns
-------
alpha : np.ndarray
Relative bandwidth
"""
if len(freqs) <= 1:
raise ParameterError(f"2 or more frequencies are required to compute bandwidths. Given freqs={freqs}")
# Approximate the local octave resolution around each frequency
bpo = np.empty_like(freqs)
logf = np.log2(freqs)
# Reflect at the lowest and highest frequencies
bpo[0] = 1 / (logf[1] - logf[0])
bpo[-1] = 1 / (logf[-1] - logf[-2])
# For everything else, do a centered difference
bpo[1:-1] = 2 / (logf[2:] - logf[:-2])
# Compute relative bandwidths
alpha = (2.0 ** (2 / bpo) - 1) / (2.0 ** (2 / bpo) + 1)
return alpha
[docs]@cache(level=10)
def wavelet(
*,
freqs: np.ndarray,
sr: float = 22050,
window: _WindowSpec = "hann",
filter_scale: float = 1,
pad_fft: bool = True,
norm: Optional[float] = 1,
dtype: DTypeLike = np.complex64,
gamma: float = 0,
alpha: Optional[float] = None,
**kwargs: Any,
) -> Tuple[np.ndarray, np.ndarray]:
"""Construct a wavelet basis using windowed complex sinusoids.
This function constructs a wavelet filterbank at a specified set of center
frequencies.
Parameters
----------
freqs : np.ndarray (positive)
Center frequencies of the filters (in Hz).
Must be in ascending order.
sr : number > 0 [scalar]
Audio sampling rate
window : string, tuple, number, or function
Windowing function to apply to filters.
filter_scale : float > 0 [scalar]
Scale of filter windows.
Small values (<1) use shorter windows for higher temporal resolution.
pad_fft : boolean
Center-pad all filters up to the nearest integral power of 2.
By default, padding is done with zeros, but this can be overridden
by setting the ``mode=`` field in *kwargs*.
norm : {inf, -inf, 0, float > 0}
Type of norm to use for basis function normalization.
See librosa.util.normalize
gamma : number >= 0
Bandwidth offset for variable-Q transforms.
dtype : np.dtype
The data type of the output basis.
By default, uses 64-bit (single precision) complex floating point.
alpha : number > 0 [optional]
Optional pre-computed relative bandwidth parameter.
Note that this must be provided if ``len(freqs)==1`` because bandwidth cannot be
inferred from a single frequency.
Otherwise, if left unspecified, it will be automatically derived by the rules
specified above.
**kwargs : additional keyword arguments
Arguments to `np.pad()` when ``pad==True``.
Returns
-------
filters : np.ndarray, ``len(filters) == n_bins``
each ``filters[i]`` is a (complex) time-domain filter
lengths : np.ndarray, ``len(lengths) == n_bins``
The (fractional) length of each filter in samples
Notes
-----
This function caches at level 10.
See Also
--------
wavelet_lengths
librosa.cqt
librosa.vqt
librosa.util.normalize
Examples
--------
Create a constant-Q basis
>>> freqs = librosa.cqt_frequencies(n_bins=84, fmin=librosa.note_to_hz('C1'))
>>> basis, lengths = librosa.filters.wavelet(freqs=freqs, sr=22050)
Plot one octave of filters in time and frequency
>>> import matplotlib.pyplot as plt
>>> basis, lengths = librosa.filters.wavelet(freqs=freqs, sr=22050)
>>> fig, ax = plt.subplots(nrows=2, figsize=(10, 6))
>>> notes = librosa.midi_to_note(np.arange(24, 24 + len(basis)))
>>> for i, (f, n) in enumerate(zip(basis, notes[:12])):
... f_scale = librosa.util.normalize(f) / 2
... ax[0].plot(i + f_scale.real)
... ax[0].plot(i + f_scale.imag, linestyle=':')
>>> ax[0].set(yticks=np.arange(len(notes[:12])), yticklabels=notes[:12],
... ylabel='CQ filters',
... title='CQ filters (one octave, time domain)',
... xlabel='Time (samples at 22050 Hz)')
>>> ax[0].legend(['Real', 'Imaginary'])
>>> F = np.abs(np.fft.fftn(basis, axes=[-1]))
>>> # Keep only the positive frequencies
>>> F = F[:, :(1 + F.shape[1] // 2)]
>>> librosa.display.specshow(F, x_axis='linear', y_axis='cqt_note', ax=ax[1])
>>> ax[1].set(ylabel='CQ filters', title='CQ filter magnitudes (frequency domain)')
"""
# Pass-through parameters to get the filter lengths
lengths, _ = wavelet_lengths(
freqs=freqs,
sr=sr,
window=window,
filter_scale=filter_scale,
gamma=gamma,
alpha=alpha,
)
# Build the filters
filters = []
for ilen, freq in zip(lengths, freqs):
# Build the filter: note, length will be ceil(ilen)
sig = util.phasor(
np.arange(-ilen // 2, ilen // 2, dtype=float) * 2 * np.pi * freq / sr
)
# Apply the windowing function
sig *= __float_window(window)(len(sig))
# Normalize
sig = util.normalize(sig, norm=norm)
filters.append(sig)
# Pad and stack
max_len = max(lengths)
if pad_fft:
max_len = int(2.0 ** (np.ceil(np.log2(max_len))))
else:
max_len = int(np.ceil(max_len))
filters = np.asarray(
[util.pad_center(filt, size=max_len, **kwargs) for filt in filters], dtype=dtype
)
return filters, lengths
[docs]@cache(level=10)
def cq_to_chroma(
n_input: int,
*,
bins_per_octave: int = 12,
n_chroma: int = 12,
fmin: Optional[_FloatLike_co] = None,
window: Optional[np.ndarray] = None,
base_c: bool = True,
dtype: DTypeLike = np.float32,
) -> np.ndarray:
"""Construct a linear transformation matrix to map Constant-Q bins
onto chroma bins (i.e., pitch classes).
Parameters
----------
n_input : int > 0 [scalar]
Number of input components (CQT bins)
bins_per_octave : int > 0 [scalar]
How many bins per octave in the CQT
n_chroma : int > 0 [scalar]
Number of output bins (per octave) in the chroma
fmin : None or float > 0
Center frequency of the first constant-Q channel.
Default: 'C1' ~= 32.7 Hz
window : None or np.ndarray
If provided, the cq_to_chroma filter bank will be
convolved with ``window``.
base_c : bool
If True, the first chroma bin will start at 'C'
If False, the first chroma bin will start at 'A'
dtype : np.dtype
The data type of the output basis.
By default, uses 32-bit (single-precision) floating point.
Returns
-------
cq_to_chroma : np.ndarray [shape=(n_chroma, n_input)]
Transformation matrix: ``Chroma = np.dot(cq_to_chroma, CQT)``
Raises
------
ParameterError
If ``n_input`` is not an integer multiple of ``n_chroma``
Notes
-----
This function caches at level 10.
Examples
--------
Get a CQT, and wrap bins to chroma
>>> y, sr = librosa.load(librosa.ex('trumpet'))
>>> CQT = np.abs(librosa.cqt(y, sr=sr))
>>> chroma_map = librosa.filters.cq_to_chroma(CQT.shape[0])
>>> chromagram = chroma_map.dot(CQT)
>>> # Max-normalize each time step
>>> chromagram = librosa.util.normalize(chromagram, axis=0)
>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots(nrows=3, sharex=True)
>>> imgcq = librosa.display.specshow(librosa.amplitude_to_db(CQT,
... ref=np.max),
... y_axis='cqt_note', x_axis='time',
... ax=ax[0])
>>> ax[0].set(title='CQT Power')
>>> ax[0].label_outer()
>>> librosa.display.specshow(chromagram, y_axis='chroma', x_axis='time',
... ax=ax[1])
>>> ax[1].set(title='Chroma (wrapped CQT)')
>>> ax[1].label_outer()
>>> chroma = librosa.feature.chroma_stft(y=y, sr=sr)
>>> imgchroma = librosa.display.specshow(chroma, y_axis='chroma', x_axis='time', ax=ax[2])
>>> ax[2].set(title='librosa.feature.chroma_stft')
"""
# How many fractional bins are we merging?
n_merge = float(bins_per_octave) / n_chroma
fmin_: _FloatLike_co
if fmin is None:
fmin_ = note_to_hz("C1")
else:
fmin_ = fmin
if np.mod(n_merge, 1) != 0:
raise ParameterError(
"Incompatible CQ merge: "
"input bins must be an "
"integer multiple of output bins."
)
# Tile the identity to merge fractional bins
cq_to_ch = np.repeat(np.eye(n_chroma), int(n_merge), axis=1)
# Roll it left to center on the target bin
cq_to_ch = np.roll(cq_to_ch, -int(n_merge // 2), axis=1)
# How many octaves are we repeating?
n_octaves = np.ceil(float(n_input) / bins_per_octave)
# Repeat and trim
cq_to_ch = np.tile(cq_to_ch, int(n_octaves))[:, :n_input]
# What's the note number of the first bin in the CQT?
# midi uses 12 bins per octave here
midi_0 = np.mod(hz_to_midi(fmin_), 12)
if base_c:
# rotate to C
roll = midi_0
else:
# rotate to A
roll = midi_0 - 9
# Adjust the roll in terms of how many chroma we want out
# We need to be careful with rounding here
roll = int(np.round(roll * (n_chroma / 12.0)))
# Apply the roll
cq_to_ch = np.roll(cq_to_ch, roll, axis=0).astype(dtype)
if window is not None:
cq_to_ch = scipy.signal.convolve(cq_to_ch, np.atleast_2d(window), mode="same")
return cq_to_ch
[docs]@cache(level=10)
def window_bandwidth(window: _WindowSpec, n: int = 1000) -> float:
"""Get the equivalent noise bandwidth (ENBW) of a window function.
The ENBW of a window is defined by [#]_ (equation 11) as the normalized
ratio of the sum of squares to the square of sums::
enbw = n * sum(window**2) / sum(window)**2
.. [#] Harris, F. J.
"On the use of windows for harmonic analysis with the discrete Fourier transform."
Proceedings of the IEEE, 66(1), 51-83. 1978.
Parameters
----------
window : callable or string
A window function, or the name of a window function,
e.g.: `scipy.signal.hann` or `'boxcar'`
n : int > 0
The number of coefficients to use in estimating the
window bandwidth
Returns
-------
bandwidth : float
The equivalent noise bandwidth (in FFT bins) of the
given window function
Notes
-----
This function caches at level 10.
See Also
--------
get_window
"""
if hasattr(window, "__name__"):
key = window.__name__
else:
key = window
if key not in WINDOW_BANDWIDTHS:
win = get_window(window, n)
WINDOW_BANDWIDTHS[key] = (
n * np.sum(win**2) / (np.sum(win) ** 2 + util.tiny(win))
)
return WINDOW_BANDWIDTHS[key]
[docs]@cache(level=10)
def get_window(
window: _WindowSpec,
Nx: int,
*,
fftbins: Optional[bool] = True,
) -> np.ndarray:
"""Compute a window function.
This is a wrapper for `scipy.signal.get_window` that additionally
supports callable or pre-computed windows.
Parameters
----------
window : string, tuple, number, callable, or list-like
The window specification:
- If string, it's the name of the window function (e.g., `'hann'`)
- If tuple, it's the name of the window function and any parameters
(e.g., `('kaiser', 4.0)`)
- If numeric, it is treated as the beta parameter of the `'kaiser'`
window, as in `scipy.signal.get_window`.
- If callable, it's a function that accepts one integer argument
(the window length)
- If list-like, it's a pre-computed window of the correct length `Nx`
Nx : int > 0
The length of the window
fftbins : bool, optional
If True (default), create a periodic window for use with FFT
If False, create a symmetric window for filter design applications.
Returns
-------
get_window : np.ndarray
A window of length `Nx` and type `window`
See Also
--------
scipy.signal.get_window
Notes
-----
This function caches at level 10.
Raises
------
ParameterError
If `window` is supplied as a vector of length != `n_fft`,
or is otherwise mis-specified.
"""
if callable(window):
return window(Nx)
elif isinstance(window, (str, tuple)) or np.isscalar(window):
# TODO: if we add custom window functions in librosa, call them here
win: np.ndarray = scipy.signal.get_window(window, Nx, fftbins=fftbins)
return win
elif isinstance(window, (np.ndarray, list)):
if len(window) == Nx:
return np.asarray(window)
raise ParameterError(f"Window size mismatch: {len(window):d} != {Nx:d}")
else:
raise ParameterError(f"Invalid window specification: {window!r}")
@cache(level=10)
def _multirate_fb(
center_freqs: Optional[np.ndarray] = None,
sample_rates: Optional[np.ndarray] = None,
Q: float = 25.0,
passband_ripple: float = 1,
stopband_attenuation: float = 50,
ftype: str = "ellip",
flayout: str = "sos",
) -> Tuple[List[Any], np.ndarray]:
r"""Construct a multirate filterbank.
A filter bank consists of multiple band-pass filters which divide the input signal
into subbands. In the case of a multirate filter bank, the band-pass filters
operate with resampled versions of the input signal, e.g. to keep the length
of a filter constant while shifting its center frequency.
This implementation uses `scipy.signal.iirdesign` to design the filters.
Parameters
----------
center_freqs : np.ndarray [shape=(n,), dtype=float]
Center frequencies of the filter kernels.
Also defines the number of filters in the filterbank.
sample_rates : np.ndarray [shape=(n,), dtype=float]
Samplerate for each filter (used for multirate filterbank).
Q : float
Q factor (influences the filter bandwidth).
passband_ripple : float
The maximum loss in the passband (dB)
See `scipy.signal.iirdesign` for details.
stopband_attenuation : float
The minimum attenuation in the stopband (dB)
See `scipy.signal.iirdesign` for details.
ftype : str
The type of IIR filter to design
See `scipy.signal.iirdesign` for details.
flayout : string
Valid `output` argument for `scipy.signal.iirdesign`.
- If `ba`, returns numerators/denominators of the transfer functions,
used for filtering with `scipy.signal.filtfilt`.
Can be unstable for high-order filters.
- If `sos`, returns a series of second-order filters,
used for filtering with `scipy.signal.sosfiltfilt`.
Minimizes numerical precision errors for high-order filters, but is slower.
- If `zpk`, returns zeros, poles, and system gains of the transfer functions.
Returns
-------
filterbank : list [shape=(n,), dtype=float]
Each list entry comprises the filter coefficients for a single filter.
sample_rates : np.ndarray [shape=(n,), dtype=float]
Samplerate for each filter.
Notes
-----
This function caches at level 10.
See Also
--------
scipy.signal.iirdesign
Raises
------
ParameterError
If ``center_freqs`` is ``None``.
If ``sample_rates`` is ``None``.
If ``center_freqs.shape`` does not match ``sample_rates.shape``.
"""
if center_freqs is None:
raise ParameterError("center_freqs must be provided.")
if sample_rates is None:
raise ParameterError("sample_rates must be provided.")
if center_freqs.shape != sample_rates.shape:
raise ParameterError(
"Number of provided center_freqs and sample_rates must be equal."
)
nyquist = 0.5 * sample_rates
filter_bandwidths = center_freqs / float(Q)
filterbank = []
for cur_center_freq, cur_nyquist, cur_bw in zip(
center_freqs, nyquist, filter_bandwidths
):
passband_freqs = [
cur_center_freq - 0.5 * cur_bw,
cur_center_freq + 0.5 * cur_bw,
] / cur_nyquist
stopband_freqs = [
cur_center_freq - cur_bw,
cur_center_freq + cur_bw,
] / cur_nyquist
cur_filter = scipy.signal.iirdesign(
passband_freqs,
stopband_freqs,
passband_ripple,
stopband_attenuation,
analog=False,
ftype=ftype,
output=flayout,
)
filterbank.append(cur_filter)
return filterbank, sample_rates
[docs]@cache(level=10)
def mr_frequencies(tuning: float) -> Tuple[np.ndarray, np.ndarray]:
r"""Generate center frequencies and sample rate pairs.
This function will return center frequency and corresponding sample rates
to obtain similar pitch filterbank settings as described in [#]_.
Instead of starting with MIDI pitch `A0`, we start with `C0`.
.. [#] Müller, Meinard.
"Information Retrieval for Music and Motion."
Springer Verlag. 2007.
Parameters
----------
tuning : float [scalar]
Tuning deviation from A440, measure as a fraction of the equally
tempered semitone (1/12 of an octave).
Returns
-------
center_freqs : np.ndarray [shape=(n,), dtype=float]
Center frequencies of the filter kernels.
Also defines the number of filters in the filterbank.
sample_rates : np.ndarray [shape=(n,), dtype=float]
Sample rate for each filter, used for multirate filterbank.
Notes
-----
This function caches at level 10.
See Also
--------
librosa.filters.semitone_filterbank
"""
center_freqs = midi_to_hz(np.arange(24 + tuning, 109 + tuning))
sample_rates = np.asarray(
len(np.arange(0, 36))
* [
882.0,
]
+ len(np.arange(36, 70))
* [
4410.0,
]
+ len(np.arange(70, 85))
* [
22050.0,
]
)
return center_freqs, sample_rates
[docs]def semitone_filterbank(
*,
center_freqs: Optional[np.ndarray] = None,
tuning: float = 0.0,
sample_rates: Optional[np.ndarray] = None,
flayout: str = "ba",
**kwargs: Any,
) -> Tuple[List[Any], np.ndarray]:
r"""Construct a multi-rate bank of infinite-impulse response (IIR)
band-pass filters at user-defined center frequencies and sample rates.
By default, these center frequencies are set equal to the 88 fundamental
frequencies of the grand piano keyboard, according to a pitch tuning standard
of A440, that is, note A above middle C set to 440 Hz. The center frequencies
are tuned to the twelve-tone equal temperament, which means that they grow
exponentially at a rate of 2**(1/12), that is, twelve notes per octave.
The A440 tuning can be changed by the user while keeping twelve-tone equal
temperament. While A440 is currently the international standard in the music
industry (ISO 16), some orchestras tune to A441-A445, whereas baroque musicians
tune to A415.
See [#]_ for details.
.. [#] Müller, Meinard.
"Information Retrieval for Music and Motion."
Springer Verlag. 2007.
Parameters
----------
center_freqs : np.ndarray [shape=(n,), dtype=float]
Center frequencies of the filter kernels.
Also defines the number of filters in the filterbank.
tuning : float [scalar]
Tuning deviation from A440 as a fraction of a semitone (1/12 of an octave
in equal temperament).
sample_rates : np.ndarray [shape=(n,), dtype=float]
Sample rates of each filter in the multirate filterbank.
flayout : string
- If `ba`, the standard difference equation is used for filtering with `scipy.signal.filtfilt`.
Can be unstable for high-order filters.
- If `sos`, 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.
**kwargs : additional keyword arguments
Additional arguments to the private function `_multirate_fb()`.
Returns
-------
filterbank : list [shape=(n,), dtype=float]
Each list entry contains the filter coefficients for a single filter.
fb_sample_rates : np.ndarray [shape=(n,), dtype=float]
Sample rate for each filter.
See Also
--------
librosa.cqt
librosa.iirt
librosa.filters.mr_frequencies
scipy.signal.iirdesign
Examples
--------
>>> import matplotlib.pyplot as plt
>>> import numpy as np
>>> import scipy.signal
>>> semitone_filterbank, sample_rates = librosa.filters.semitone_filterbank(
... center_freqs=librosa.midi_to_hz(np.arange(60, 72)),
... sample_rates=np.repeat(4410.0, 12),
... flayout='sos'
... )
>>> magnitudes = []
>>> for cur_sr, cur_filter in zip(sample_rates, semitone_filterbank):
... w, h = scipy.signal.sosfreqz(cur_filter,fs=cur_sr, worN=1025)
... magnitudes.append(20 * np.log10(np.abs(h)))
>>> fig, ax = plt.subplots(figsize=(12,6))
>>> img = librosa.display.specshow(
... np.array(magnitudes),
... x_axis="hz",
... sr=4410,
... y_coords=librosa.midi_to_hz(np.arange(60, 72)),
... vmin=-60,
... vmax=3,
... ax=ax
... )
>>> fig.colorbar(img, ax=ax, format="%+2.f dB", label="Magnitude (dB)")
>>> ax.set(
... xlim=[200, 600],
... yticks=librosa.midi_to_hz(np.arange(60, 72)),
... title='Magnitude Responses of the Pitch Filterbank',
... xlabel='Frequency (Hz)',
... ylabel='Semitone filter center frequency (Hz)'
... )
"""
if (center_freqs is None) and (sample_rates is None):
center_freqs, sample_rates = mr_frequencies(tuning)
filterbank, fb_sample_rates = _multirate_fb(
center_freqs=center_freqs, sample_rates=sample_rates, flayout=flayout, **kwargs
)
return filterbank, fb_sample_rates
@jit(nopython=True, cache=True)
def __window_ss_fill(x, win_sq, n_frames, hop_length): # pragma: no cover
"""Compute the sum-square envelope of a window."""
n = len(x)
n_fft = len(win_sq)
for i in range(n_frames):
sample = i * hop_length
x[sample : min(n, sample + n_fft)] += win_sq[: max(0, min(n_fft, n - sample))]
[docs]def window_sumsquare(
*,
window: _WindowSpec,
n_frames: int,
hop_length: int = 512,
win_length: Optional[int] = None,
n_fft: int = 2048,
dtype: DTypeLike = np.float32,
norm: Optional[float] = None,
) -> np.ndarray:
"""Compute the sum-square envelope of a window function at a given hop length.
This is used to estimate modulation effects induced by windowing observations
in short-time Fourier transforms.
Parameters
----------
window : string, tuple, number, callable, or list-like
Window specification, as in `get_window`
n_frames : int > 0
The number of analysis frames
hop_length : int > 0
The number of samples to advance between frames
win_length : [optional]
The length of the window function. By default, this matches ``n_fft``.
n_fft : int > 0
The length of each analysis frame.
dtype : np.dtype
The data type of the output
norm : {np.inf, -np.inf, 0, float > 0, None}
Normalization mode used in window construction.
Note that this does not affect the squaring operation.
Returns
-------
wss : np.ndarray, shape=``(n_fft + hop_length * (n_frames - 1))``
The sum-squared envelope of the window function
Examples
--------
For a fixed frame length (2048), compare modulation effects for a Hann window
at different hop lengths:
>>> n_frames = 50
>>> wss_256 = librosa.filters.window_sumsquare(window='hann', n_frames=n_frames, hop_length=256)
>>> wss_512 = librosa.filters.window_sumsquare(window='hann', n_frames=n_frames, hop_length=512)
>>> wss_1024 = librosa.filters.window_sumsquare(window='hann', n_frames=n_frames, hop_length=1024)
>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots(nrows=3, sharey=True)
>>> ax[0].plot(wss_256)
>>> ax[0].set(title='hop_length=256')
>>> ax[1].plot(wss_512)
>>> ax[1].set(title='hop_length=512')
>>> ax[2].plot(wss_1024)
>>> ax[2].set(title='hop_length=1024')
"""
if win_length is None:
win_length = n_fft
n = n_fft + hop_length * (n_frames - 1)
x = np.zeros(n, dtype=dtype)
# Compute the squared window at the desired length
win_sq = get_window(window, win_length)
win_sq = util.normalize(win_sq, norm=norm) ** 2
win_sq = util.pad_center(win_sq, size=n_fft)
# Fill the envelope
__window_ss_fill(x, win_sq, n_frames, hop_length)
return x
[docs]@cache(level=10)
def diagonal_filter(
window: _WindowSpec,
n: int,
*,
slope: float = 1.0,
angle: Optional[float] = None,
zero_mean: bool = False,
) -> np.ndarray:
"""Build a two-dimensional diagonal filter.
This is primarily used for smoothing recurrence or self-similarity matrices.
Parameters
----------
window : string, tuple, number, callable, or list-like
The window function to use for the filter.
See `get_window` for details.
Note that the window used here should be non-negative.
n : int > 0
the length of the filter
slope : float
The slope of the diagonal filter to produce
angle : float or None
If given, the slope parameter is ignored,
and angle directly sets the orientation of the filter (in radians).
Otherwise, angle is inferred as `arctan(slope)`.
zero_mean : bool
If True, a zero-mean filter is used.
Otherwise, a non-negative averaging filter is used.
This should be enabled if you want to enhance paths and suppress
blocks.
Returns
-------
kernel : np.ndarray, shape=[(m, m)]
The 2-dimensional filter kernel
Notes
-----
This function caches at level 10.
"""
if angle is None:
angle = np.arctan(slope)
win = np.diag(get_window(window, n, fftbins=False))
if not np.isclose(angle, np.pi / 4):
win = scipy.ndimage.rotate(
win, 45 - angle * 180 / np.pi, order=5, prefilter=False
)
np.clip(win, 0, None, out=win)
win /= win.sum()
if zero_mean:
win -= win.mean()
return win