Caution
You're reading an old version of this documentation. If you want up-to-date information, please have a look at 0.10.2.
Source code for librosa.filters
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Filters
=======
Filter bank construction
------------------------
.. autosummary::
:toctree: generated/
mel
chroma
constant_q
semitone_filterbank
Window functions
----------------
.. autosummary::
:toctree: generated/
window_bandwidth
get_window
Miscellaneous
-------------
.. autosummary::
:toctree: generated/
constant_q_lengths
cq_to_chroma
mr_frequencies
window_sumsquare
diagonal_filter
"""
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 .core.convert import note_to_hz, hz_to_midi, midi_to_hz, hz_to_octs
from .core.convert import fft_frequencies, mel_frequencies
__all__ = [
"mel",
"chroma",
"constant_q",
"constant_q_lengths",
"cq_to_chroma",
"window_bandwidth",
"get_window",
"mr_frequencies",
"semitone_filterbank",
"window_sumsquare",
"diagonal_filter",
]
# 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,
"hanning": 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,
n_fft,
n_mels=128,
fmin=0.0,
fmax=None,
htk=False,
norm="slaney",
dtype=np.float32,
):
"""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(22050, 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(22050, 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 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:
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."
)
return weights
[docs]@cache(level=10)
def chroma(
sr,
n_fft,
n_chroma=12,
tuning=0.0,
ctroct=5.0,
octwidth=2,
norm=2,
base_c=True,
dtype=np.float32,
):
"""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(22050, 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(22050, 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(22050, 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):
"""Decorator function for windows with fractional input.
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):
"""The wrapped 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]@cache(level=10)
def constant_q(
sr,
fmin=None,
n_bins=84,
bins_per_octave=12,
window="hann",
filter_scale=1,
pad_fft=True,
norm=1,
dtype=np.complex64,
gamma=0,
**kwargs,
):
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.
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
--------
constant_q_lengths
librosa.cqt
librosa.vqt
librosa.util.normalize
Examples
--------
Use a shorter window for each filter
>>> basis, lengths = librosa.filters.constant_q(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(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,
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 = np.exp(
np.arange(-ilen // 2, ilen // 2, dtype=float) * 1j * 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, max_len, **kwargs) for filt in filters], dtype=dtype
)
return filters, np.asarray(lengths)
[docs]@cache(level=10)
def constant_q_lengths(
sr, fmin, n_bins=84, bins_per_octave=12, window="hann", filter_scale=1, gamma=0
):
r"""Return length of each filter in a constant-Q basis.
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.
Returns
-------
lengths : np.ndarray
The length of each filter.
Notes
-----
This function caches at level 10.
See Also
--------
constant_q
librosa.cqt
"""
if fmin <= 0:
raise ParameterError("fmin must be 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")
# Q should be capitalized here, so we suppress the name warning
# pylint: disable=invalid-name
alpha = 2.0 ** (1.0 / bins_per_octave) - 1.0
Q = float(filter_scale) / alpha
# Compute the frequencies
freq = fmin * (2.0 ** (np.arange(n_bins, dtype=float) / bins_per_octave))
if freq[-1] * (1 + 0.5 * window_bandwidth(window) / Q) > sr / 2.0:
raise ParameterError("Filter pass-band lies beyond Nyquist")
# Convert frequencies to filter lengths
lengths = Q * sr / (freq + gamma / alpha)
return lengths
[docs]@cache(level=10)
def cq_to_chroma(
n_input,
bins_per_octave=12,
n_chroma=12,
fmin=None,
window=None,
base_c=True,
dtype=np.float32,
):
"""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
if fmin is None:
fmin = note_to_hz("C1")
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), 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(np.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, n=1000):
"""Get the equivalent noise bandwidth of a window function.
Parameters
----------
window : callable or string
A window function, or the name of a window function.
Examples:
- scipy.signal.hann
- '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(np.abs(win)) ** 2
return WINDOW_BANDWIDTHS[key]
[docs]@cache(level=10)
def get_window(window, Nx, fftbins=True):
"""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
return scipy.signal.get_window(window, Nx, fftbins=fftbins)
elif isinstance(window, (np.ndarray, list)):
if len(window) == Nx:
return np.asarray(window)
raise ParameterError(
"Window size mismatch: " "{:d} != {:d}".format(len(window), Nx)
)
else:
raise ParameterError("Invalid window specification: {}".format(window))
@cache(level=10)
def _multirate_fb(
center_freqs=None,
sample_rates=None,
Q=25.0,
passband_ripple=1,
stopband_attenuation=50,
ftype="ellip",
flayout="sos",
):
r"""Helper function to 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 bandwith).
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):
r"""Helper function for generating center frequency 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, ]
+ len(np.arange(36, 70)) * [4410, ]
+ len(np.arange(70, 85)) * [22050, ]
)
return center_freqs, sample_rates
[docs]def semitone_filterbank(
center_freqs=None, tuning=0.0, sample_rates=None, flayout="ba", **kwargs
):
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()
>>> fig, ax = plt.subplots()
>>> for cur_sr, cur_filter in zip(sample_rates, semitone_filterbank):
... w, h = scipy.signal.freqz(cur_filter[0], cur_filter[1], worN=2000)
... ax.semilogx((cur_sr / (2 * np.pi)) * w, 20 * np.log10(abs(h)))
>>> ax.set(xlim=[20, 10e3], ylim=[-60, 3], title='Magnitude Responses of the Pitch Filterbank',
... xlabel='Log-Frequency (Hz)', ylabel='Magnitude (dB)')
"""
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
"""Helper function for window sum-square calculation."""
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,
n_frames,
hop_length=512,
win_length=None,
n_fft=2048,
dtype=np.float32,
norm=None,
):
"""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
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('hann', n_frames, hop_length=256)
>>> wss_512 = librosa.filters.window_sumsquare('hann', n_frames, hop_length=512)
>>> wss_1024 = librosa.filters.window_sumsquare('hann', 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, 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, n, slope=1.0, angle=None, zero_mean=False):
"""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