#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""Pitch-tracking and tuning estimation"""
import warnings
import numpy as np
import scipy
import numba
from .spectrum import _spectrogram
from . import convert
from .._cache import cache
from .. import util
from .. import sequence
from ..util.exceptions import ParameterError
from numpy.typing import ArrayLike
from typing import Any, Callable, Optional, Tuple, Union
from .._typing import _WindowSpec, _PadMode, _PadModeSTFT
__all__ = ["estimate_tuning", "pitch_tuning", "piptrack", "yin", "pyin"]
[docs]def estimate_tuning(
*,
y: Optional[np.ndarray] = None,
sr: float = 22050,
S: Optional[np.ndarray] = None,
n_fft: Optional[int] = 2048,
resolution: float = 0.01,
bins_per_octave: int = 12,
**kwargs: Any,
) -> float:
"""Estimate the tuning of an audio time series or spectrogram input.
Parameters
----------
y : np.ndarray [shape=(..., n)] or None
audio signal. Multi-channel is supported..
sr : number > 0 [scalar]
audio sampling rate of ``y``
S : np.ndarray [shape=(..., d, t)] or None
magnitude or power spectrogram
n_fft : int > 0 [scalar] or None
number of FFT bins to use, if ``y`` is provided.
resolution : float in `(0, 1)`
Resolution of the tuning as a fraction of a bin.
0.01 corresponds to measurements in cents.
bins_per_octave : int > 0 [scalar]
How many frequency bins per octave
**kwargs : additional keyword arguments
Additional arguments passed to `piptrack`
Returns
-------
tuning: float in `[-0.5, 0.5)`
estimated tuning deviation (fractions of a bin).
Note that if multichannel input is provided, a single tuning estimate is provided spanning all
channels.
See Also
--------
piptrack : Pitch tracking by parabolic interpolation
Examples
--------
With time-series input
>>> y, sr = librosa.load(librosa.ex('trumpet'))
>>> librosa.estimate_tuning(y=y, sr=sr)
-0.08000000000000002
In tenths of a cent
>>> librosa.estimate_tuning(y=y, sr=sr, resolution=1e-3)
-0.016000000000000014
Using spectrogram input
>>> S = np.abs(librosa.stft(y))
>>> librosa.estimate_tuning(S=S, sr=sr)
-0.08000000000000002
Using pass-through arguments to `librosa.piptrack`
>>> librosa.estimate_tuning(y=y, sr=sr, n_fft=8192,
... fmax=librosa.note_to_hz('G#9'))
-0.08000000000000002
"""
pitch, mag = piptrack(y=y, sr=sr, S=S, n_fft=n_fft, **kwargs)
# Only count magnitude where frequency is > 0
pitch_mask = pitch > 0
if pitch_mask.any():
threshold = np.median(mag[pitch_mask])
else:
threshold = 0.0
return pitch_tuning(
pitch[(mag >= threshold) & pitch_mask],
resolution=resolution,
bins_per_octave=bins_per_octave,
)
[docs]def pitch_tuning(
frequencies: ArrayLike, *, resolution: float = 0.01, bins_per_octave: int = 12
) -> float:
"""Given a collection of pitches, estimate its tuning offset
(in fractions of a bin) relative to A440=440.0Hz.
Parameters
----------
frequencies : array-like, float
A collection of frequencies detected in the signal.
See `piptrack`
resolution : float in `(0, 1)`
Resolution of the tuning as a fraction of a bin.
0.01 corresponds to cents.
bins_per_octave : int > 0 [scalar]
How many frequency bins per octave
Returns
-------
tuning: float in `[-0.5, 0.5)`
estimated tuning deviation (fractions of a bin)
See Also
--------
estimate_tuning : Estimating tuning from time-series or spectrogram input
Examples
--------
>>> # Generate notes at +25 cents
>>> freqs = librosa.cqt_frequencies(n_bins=24, fmin=55, tuning=0.25)
>>> librosa.pitch_tuning(freqs)
0.25
>>> # Track frequencies from a real spectrogram
>>> y, sr = librosa.load(librosa.ex('trumpet'))
>>> freqs, times, mags = librosa.reassigned_spectrogram(y, sr=sr,
... fill_nan=True)
>>> # Select out pitches with high energy
>>> freqs = freqs[mags > np.median(mags)]
>>> librosa.pitch_tuning(freqs)
-0.07
"""
frequencies = np.atleast_1d(frequencies)
# Trim out any DC components
frequencies = frequencies[frequencies > 0]
if not np.any(frequencies):
warnings.warn(
"Trying to estimate tuning from empty frequency set.", stacklevel=2
)
return 0.0
# Compute the residual relative to the number of bins
residual = np.mod(bins_per_octave * convert.hz_to_octs(frequencies), 1.0)
# Are we on the wrong side of the semitone?
# A residual of 0.95 is more likely to be a deviation of -0.05
# from the next tone up.
residual[residual >= 0.5] -= 1.0
bins = np.linspace(-0.5, 0.5, int(np.ceil(1.0 / resolution)) + 1)
counts, tuning = np.histogram(residual, bins)
# return the histogram peak
tuning_est: float = tuning[np.argmax(counts)]
return tuning_est
[docs]@cache(level=30)
def piptrack(
*,
y: Optional[np.ndarray] = None,
sr: float = 22050,
S: Optional[np.ndarray] = None,
n_fft: Optional[int] = 2048,
hop_length: Optional[int] = None,
fmin: float = 150.0,
fmax: float = 4000.0,
threshold: float = 0.1,
win_length: Optional[int] = None,
window: _WindowSpec = "hann",
center: bool = True,
pad_mode: _PadModeSTFT = "constant",
ref: Optional[Union[float, Callable]] = None,
) -> Tuple[np.ndarray, np.ndarray]:
"""Pitch tracking on thresholded parabolically-interpolated STFT.
This implementation uses the parabolic interpolation method described by [#]_.
.. [#] https://ccrma.stanford.edu/~jos/sasp/Sinusoidal_Peak_Interpolation.html
Parameters
----------
y : np.ndarray [shape=(..., n)] or None
audio signal. Multi-channel is supported..
sr : number > 0 [scalar]
audio sampling rate of ``y``
S : np.ndarray [shape=(..., d, t)] or None
magnitude or power spectrogram
n_fft : int > 0 [scalar] or None
number of FFT bins to use, if ``y`` is provided.
hop_length : int > 0 [scalar] or None
number of samples to hop
threshold : float in `(0, 1)`
A bin in spectrum ``S`` is considered a pitch when it is greater than
``threshold * ref(S)``.
By default, ``ref(S)`` is taken to be ``max(S, axis=0)`` (the maximum value in
each column).
fmin : float > 0 [scalar]
lower frequency cutoff.
fmax : float > 0 [scalar]
upper frequency cutoff.
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.
See also: `np.pad`.
ref : scalar or callable [default=np.max]
If scalar, the reference value against which ``S`` is compared for determining
pitches.
If callable, the reference value is computed as ``ref(S, axis=0)``.
Returns
-------
pitches, magnitudes : np.ndarray [shape=(..., d, t)]
Where ``d`` is the subset of FFT bins within ``fmin`` and ``fmax``.
``pitches[..., f, t]`` contains instantaneous frequency at bin
``f``, time ``t``
``magnitudes[..., f, t]`` contains the corresponding magnitudes.
Both ``pitches`` and ``magnitudes`` take value 0 at bins
of non-maximal magnitude.
Notes
-----
This function caches at level 30.
One of ``S`` or ``y`` must be provided.
If ``S`` is not given, it is computed from ``y`` using
the default parameters of `librosa.stft`.
Examples
--------
Computing pitches from a waveform input
>>> y, sr = librosa.load(librosa.ex('trumpet'))
>>> pitches, magnitudes = librosa.piptrack(y=y, sr=sr)
Or from a spectrogram input
>>> S = np.abs(librosa.stft(y))
>>> pitches, magnitudes = librosa.piptrack(S=S, sr=sr)
Or with an alternate reference value for pitch detection, where
values above the mean spectral energy in each frame are counted as pitches
>>> pitches, magnitudes = librosa.piptrack(S=S, sr=sr, threshold=1,
... ref=np.mean)
"""
# Check that we received an audio time series or STFT
S, n_fft = _spectrogram(
y=y,
S=S,
n_fft=n_fft,
hop_length=hop_length,
win_length=win_length,
window=window,
center=center,
pad_mode=pad_mode,
)
# Make sure we're dealing with magnitudes
S = np.abs(S)
# Truncate to feasible region
fmin = np.maximum(fmin, 0)
fmax = np.minimum(fmax, float(sr) / 2)
fft_freqs = convert.fft_frequencies(sr=sr, n_fft=n_fft)
# Do the parabolic interpolation everywhere,
# then figure out where the peaks are
# then restrict to the feasible range (fmin:fmax)
avg = np.gradient(S, axis=-2)
shift = _parabolic_interpolation(S, axis=-2)
# this will get us the interpolated peak value
dskew = 0.5 * avg * shift
# Pre-allocate output
pitches = np.zeros_like(S)
mags = np.zeros_like(S)
# Clip to the viable frequency range
freq_mask = (fmin <= fft_freqs) & (fft_freqs < fmax)
freq_mask = util.expand_to(freq_mask, ndim=S.ndim, axes=-2)
# Compute the column-wise local max of S after thresholding
# Find the argmax coordinates
if ref is None:
ref = np.max
if callable(ref):
ref_value = threshold * ref(S, axis=-2)
# Reinsert the frequency axis here, in case the callable doesn't
# support keepdims=True
ref_value = np.expand_dims(ref_value, -2)
else:
ref_value = np.abs(ref)
# Store pitch and magnitude
idx = np.nonzero(freq_mask & util.localmax(S * (S > ref_value), axis=-2))
pitches[idx] = (idx[-2] + shift[idx]) * float(sr) / n_fft
mags[idx] = S[idx] + dskew[idx]
return pitches, mags
def _cumulative_mean_normalized_difference(
y_frames: np.ndarray,
frame_length: int,
win_length: int,
min_period: int,
max_period: int,
) -> np.ndarray:
"""Cumulative mean normalized difference function (equation 8 in [#]_)
.. [#] De Cheveigné, Alain, and Hideki Kawahara.
"YIN, a fundamental frequency estimator for speech and music."
The Journal of the Acoustical Society of America 111.4 (2002): 1917-1930.
Parameters
----------
y_frames : np.ndarray [shape=(frame_length, n_frames)]
framed audio time series.
frame_length : int > 0 [scalar]
length of the frames in samples.
win_length : int > 0 [scalar]
length of the window for calculating autocorrelation in samples.
min_period : int > 0 [scalar]
minimum period.
max_period : int > 0 [scalar]
maximum period.
Returns
-------
yin_frames : np.ndarray [shape=(max_period-min_period+1,n_frames)]
Cumulative mean normalized difference function for each frame.
"""
# Autocorrelation.
a = np.fft.rfft(y_frames, frame_length, axis=-2)
b = np.fft.rfft(y_frames[..., win_length:0:-1, :], frame_length, axis=-2)
acf_frames = np.fft.irfft(a * b, frame_length, axis=-2)[..., win_length:, :]
acf_frames[np.abs(acf_frames) < 1e-6] = 0
# Energy terms.
energy_frames = np.cumsum(y_frames**2, axis=-2)
energy_frames = (
energy_frames[..., win_length:, :] - energy_frames[..., :-win_length, :]
)
energy_frames[np.abs(energy_frames) < 1e-6] = 0
# Difference function.
yin_frames = energy_frames[..., :1, :] + energy_frames - 2 * acf_frames
# Cumulative mean normalized difference function.
yin_numerator = yin_frames[..., min_period : max_period + 1, :]
# broadcast this shape to have leading ones
tau_range = util.expand_to(
np.arange(1, max_period + 1), ndim=yin_frames.ndim, axes=-2
)
cumulative_mean = (
np.cumsum(yin_frames[..., 1 : max_period + 1, :], axis=-2) / tau_range
)
yin_denominator = cumulative_mean[..., min_period - 1 : max_period, :]
yin_frames: np.ndarray = yin_numerator / (
yin_denominator + util.tiny(yin_denominator)
)
return yin_frames
@numba.stencil # type: ignore
def _pi_stencil(x: np.ndarray) -> np.ndarray:
"""Stencil to compute local parabolic interpolation"""
a = x[1] + x[-1] - 2 * x[0]
b = (x[1] - x[-1]) / 2
if np.abs(b) >= np.abs(a):
# If this happens, we'll shift by more than 1 bin
# Suppressing types because mypy has no idea about stencils
return 0 # type: ignore
return -b / a # type: ignore
@numba.guvectorize(
["void(float32[:], float32[:])", "void(float64[:], float64[:])"],
"(n)->(n)",
cache=True,
nopython=True,
) # type: ignore
def _pi_wrapper(x: np.ndarray, y: np.ndarray) -> None: # pragma: no cover
"""Vectorized wrapper for the parabolic interpolation stencil"""
y[:] = _pi_stencil(x)
def _parabolic_interpolation(x: np.ndarray, *, axis: int = -2) -> np.ndarray:
"""Piecewise parabolic interpolation for yin and pyin.
Parameters
----------
x : np.ndarray
array to interpolate
axis : int
axis along which to interpolate
Returns
-------
parabolic_shifts : np.ndarray [shape=x.shape]
position of the parabola optima (relative to bin indices)
Note: the shift at bin `n` is determined as 0 if the estimated
optimum is outside the range `[n-1, n+1]`.
"""
# Rotate the target axis to the end
xi = x.swapaxes(-1, axis)
# Allocate the output array and rotate target axis
shifts = np.empty_like(x)
shiftsi = shifts.swapaxes(-1, axis)
# Call the vectorized stencil
_pi_wrapper(xi, shiftsi)
# Handle the edge condition not covered by the stencil
shiftsi[..., -1] = 0
shiftsi[..., 0] = 0
return shifts
[docs]def yin(
y: np.ndarray,
*,
fmin: float,
fmax: float,
sr: float = 22050,
frame_length: int = 2048,
win_length: Optional[int] = None,
hop_length: Optional[int] = None,
trough_threshold: float = 0.1,
center: bool = True,
pad_mode: _PadMode = "constant",
) -> np.ndarray:
"""Fundamental frequency (F0) estimation using the YIN algorithm.
YIN is an autocorrelation based method for fundamental frequency estimation [#]_.
First, a normalized difference function is computed over short (overlapping) frames of audio.
Next, the first minimum in the difference function below ``trough_threshold`` is selected as
an estimate of the signal's period.
Finally, the estimated period is refined using parabolic interpolation before converting
into the corresponding frequency.
.. [#] De Cheveigné, Alain, and Hideki Kawahara.
"YIN, a fundamental frequency estimator for speech and music."
The Journal of the Acoustical Society of America 111.4 (2002): 1917-1930.
Parameters
----------
y : np.ndarray [shape=(..., n)]
audio time series. Multi-channel is supported..
fmin : number > 0 [scalar]
minimum frequency in Hertz.
The recommended minimum is ``librosa.note_to_hz('C2')`` (~65 Hz)
though lower values may be feasible.
fmax : number > fmin, <= sr/2 [scalar]
maximum frequency in Hertz.
The recommended maximum is ``librosa.note_to_hz('C7')`` (~2093 Hz)
though higher values may be feasible.
sr : number > 0 [scalar]
sampling rate of ``y`` in Hertz.
frame_length : int > 0 [scalar]
length of the frames in samples.
By default, ``frame_length=2048`` corresponds to a time scale of about 93 ms at
a sampling rate of 22050 Hz.
win_length : None or int > 0 [scalar]
length of the window for calculating autocorrelation in samples.
If ``None``, defaults to ``frame_length // 2``
hop_length : None or int > 0 [scalar]
number of audio samples between adjacent YIN predictions.
If ``None``, defaults to ``frame_length // 4``.
trough_threshold : number > 0 [scalar]
absolute threshold for peak estimation.
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.core.frames_to_samples``.
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.
If ``center=False``, this argument is ignored.
.. see also:: `np.pad`
Returns
-------
f0: np.ndarray [shape=(..., n_frames)]
time series of fundamental frequencies in Hertz.
If multi-channel input is provided, f0 curves are estimated separately for each channel.
See Also
--------
librosa.pyin :
Fundamental frequency (F0) estimation using probabilistic YIN (pYIN).
Examples
--------
Computing a fundamental frequency (F0) curve from an audio input
>>> y = librosa.chirp(fmin=440, fmax=880, duration=5.0, sr=22050)
>>> librosa.yin(y, fmin=440, fmax=880, sr=22050)
array([442.66354675, 441.95299983, 441.58010963, ...,
871.161732 , 873.99001454, 877.04297681])
"""
if fmin is None or fmax is None:
raise ParameterError('both "fmin" and "fmax" must be provided')
# Set the default window length if it is not already specified.
if win_length is None:
win_length = frame_length // 2
__check_yin_params(
sr=sr, fmax=fmax, fmin=fmin, frame_length=frame_length, win_length=win_length
)
# Set the default hop if it is not already specified.
if hop_length is None:
hop_length = frame_length // 4
# Check that audio is valid.
util.valid_audio(y, mono=False)
# Pad the time series so that frames are centered
if center:
padding = [(0, 0)] * y.ndim
padding[-1] = (frame_length // 2, frame_length // 2)
y = np.pad(y, padding, mode=pad_mode)
# Frame audio.
y_frames = util.frame(y, frame_length=frame_length, hop_length=hop_length)
# Calculate minimum and maximum periods
min_period = int(np.floor(sr / fmax))
max_period = min(int(np.ceil(sr / fmin)), frame_length - win_length - 1)
# Calculate cumulative mean normalized difference function.
yin_frames = _cumulative_mean_normalized_difference(
y_frames, frame_length, win_length, min_period, max_period
)
# Parabolic interpolation.
parabolic_shifts = _parabolic_interpolation(yin_frames)
# Find local minima.
is_trough = util.localmin(yin_frames, axis=-2)
is_trough[..., 0, :] = yin_frames[..., 0, :] < yin_frames[..., 1, :]
# Find minima below peak threshold.
is_threshold_trough = np.logical_and(is_trough, yin_frames < trough_threshold)
# Absolute threshold.
# "The solution we propose is to set an absolute threshold and choose the
# smallest value of tau that gives a minimum of d' deeper than
# this threshold. If none is found, the global minimum is chosen instead."
target_shape = list(yin_frames.shape)
target_shape[-2] = 1
global_min = np.argmin(yin_frames, axis=-2)
yin_period = np.argmax(is_threshold_trough, axis=-2)
global_min = global_min.reshape(target_shape)
yin_period = yin_period.reshape(target_shape)
no_trough_below_threshold = np.all(~is_threshold_trough, axis=-2, keepdims=True)
yin_period[no_trough_below_threshold] = global_min[no_trough_below_threshold]
# Refine peak by parabolic interpolation.
yin_period = (
min_period
+ yin_period
+ np.take_along_axis(parabolic_shifts, yin_period, axis=-2)
)[..., 0, :]
# Convert period to fundamental frequency.
f0: np.ndarray = sr / yin_period
return f0
[docs]def pyin(
y: np.ndarray,
*,
fmin: float,
fmax: float,
sr: float = 22050,
frame_length: int = 2048,
win_length: Optional[int] = None,
hop_length: Optional[int] = None,
n_thresholds: int = 100,
beta_parameters: Tuple[float, float] = (2, 18),
boltzmann_parameter: float = 2,
resolution: float = 0.1,
max_transition_rate: float = 35.92,
switch_prob: float = 0.01,
no_trough_prob: float = 0.01,
fill_na: Optional[float] = np.nan,
center: bool = True,
pad_mode: _PadMode = "constant",
) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
"""Fundamental frequency (F0) estimation using probabilistic YIN (pYIN).
pYIN [#]_ is a modificatin of the YIN algorithm [#]_ for fundamental frequency (F0) estimation.
In the first step of pYIN, F0 candidates and their probabilities are computed using the YIN algorithm.
In the second step, Viterbi decoding is used to estimate the most likely F0 sequence and voicing flags.
.. [#] Mauch, Matthias, and Simon Dixon.
"pYIN: A fundamental frequency estimator using probabilistic threshold distributions."
2014 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2014.
.. [#] De Cheveigné, Alain, and Hideki Kawahara.
"YIN, a fundamental frequency estimator for speech and music."
The Journal of the Acoustical Society of America 111.4 (2002): 1917-1930.
Parameters
----------
y : np.ndarray [shape=(..., n)]
audio time series. Multi-channel is supported.
fmin : number > 0 [scalar]
minimum frequency in Hertz.
The recommended minimum is ``librosa.note_to_hz('C2')`` (~65 Hz)
though lower values may be feasible.
fmax : number > fmin, <= sr/2 [scalar]
maximum frequency in Hertz.
The recommended maximum is ``librosa.note_to_hz('C7')`` (~2093 Hz)
though higher values may be feasible.
sr : number > 0 [scalar]
sampling rate of ``y`` in Hertz.
frame_length : int > 0 [scalar]
length of the frames in samples.
By default, ``frame_length=2048`` corresponds to a time scale of about 93 ms at
a sampling rate of 22050 Hz.
win_length : None or int > 0 [scalar]
length of the window for calculating autocorrelation in samples.
If ``None``, defaults to ``frame_length // 2``
hop_length : None or int > 0 [scalar]
number of audio samples between adjacent pYIN predictions.
If ``None``, defaults to ``frame_length // 4``.
n_thresholds : int > 0 [scalar]
number of thresholds for peak estimation.
beta_parameters : tuple
shape parameters for the beta distribution prior over thresholds.
boltzmann_parameter : number > 0 [scalar]
shape parameter for the Boltzmann distribution prior over troughs.
Larger values will assign more mass to smaller periods.
resolution : float in `(0, 1)`
Resolution of the pitch bins.
0.01 corresponds to cents.
max_transition_rate : float > 0
maximum pitch transition rate in octaves per second.
switch_prob : float in ``(0, 1)``
probability of switching from voiced to unvoiced or vice versa.
no_trough_prob : float in ``(0, 1)``
maximum probability to add to global minimum if no trough is below threshold.
fill_na : None, float, or ``np.nan``
default value for unvoiced frames of ``f0``.
If ``None``, the unvoiced frames will contain a best guess value.
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.core.frames_to_samples``.
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.
If ``center=False``, this argument is ignored.
.. see also:: `np.pad`
Returns
-------
f0: np.ndarray [shape=(..., n_frames)]
time series of fundamental frequencies in Hertz.
voiced_flag: np.ndarray [shape=(..., n_frames)]
time series containing boolean flags indicating whether a frame is voiced or not.
voiced_prob: np.ndarray [shape=(..., n_frames)]
time series containing the probability that a frame is voiced.
.. note:: If multi-channel input is provided, f0 and voicing are estimated separately for each channel.
See Also
--------
librosa.yin :
Fundamental frequency (F0) estimation using the YIN algorithm.
Examples
--------
Computing a fundamental frequency (F0) curve from an audio input
>>> y, sr = librosa.load(librosa.ex('trumpet'))
>>> f0, voiced_flag, voiced_probs = librosa.pyin(y,
... sr=sr,
... fmin=librosa.note_to_hz('C2'),
... fmax=librosa.note_to_hz('C7'))
>>> times = librosa.times_like(f0, sr=sr)
Overlay F0 over a spectrogram
>>> import matplotlib.pyplot as plt
>>> D = librosa.amplitude_to_db(np.abs(librosa.stft(y)), ref=np.max)
>>> fig, ax = plt.subplots()
>>> img = librosa.display.specshow(D, x_axis='time', y_axis='log', ax=ax)
>>> ax.set(title='pYIN fundamental frequency estimation')
>>> fig.colorbar(img, ax=ax, format="%+2.f dB")
>>> ax.plot(times, f0, label='f0', color='cyan', linewidth=3)
>>> ax.legend(loc='upper right')
"""
if fmin is None or fmax is None:
raise ParameterError('both "fmin" and "fmax" must be provided')
# Set the default window length if it is not already specified.
if win_length is None:
win_length = frame_length // 2
__check_yin_params(
sr=sr, fmax=fmax, fmin=fmin, frame_length=frame_length, win_length=win_length
)
# Set the default hop if it is not already specified.
if hop_length is None:
hop_length = frame_length // 4
# Check that audio is valid.
util.valid_audio(y, mono=False)
# Pad the time series so that frames are centered
if center:
padding = [(0, 0) for _ in y.shape]
padding[-1] = (frame_length // 2, frame_length // 2)
y = np.pad(y, padding, mode=pad_mode)
# Frame audio.
y_frames = util.frame(y, frame_length=frame_length, hop_length=hop_length)
# Calculate minimum and maximum periods
min_period = int(np.floor(sr / fmax))
max_period = min(int(np.ceil(sr / fmin)), frame_length - win_length - 1)
# Calculate cumulative mean normalized difference function.
yin_frames = _cumulative_mean_normalized_difference(
y_frames, frame_length, win_length, min_period, max_period
)
# Parabolic interpolation.
parabolic_shifts = _parabolic_interpolation(yin_frames)
# Find Yin candidates and probabilities.
# The implementation here follows the official pYIN software which
# differs from the method described in the paper.
# 1. Define the prior over the thresholds.
thresholds = np.linspace(0, 1, n_thresholds + 1)
beta_cdf = scipy.stats.beta.cdf(thresholds, beta_parameters[0], beta_parameters[1])
beta_probs = np.diff(beta_cdf)
n_bins_per_semitone = int(np.ceil(1.0 / resolution))
n_pitch_bins = int(np.floor(12 * n_bins_per_semitone * np.log2(fmax / fmin))) + 1
def _helper(a, b):
return __pyin_helper(
a,
b,
sr,
thresholds,
boltzmann_parameter,
beta_probs,
no_trough_prob,
min_period,
fmin,
n_pitch_bins,
n_bins_per_semitone,
)
helper = np.vectorize(_helper, signature="(f,t),(k,t)->(1,d,t),(j,t)")
observation_probs, voiced_prob = helper(yin_frames, parabolic_shifts)
# Construct transition matrix.
max_semitones_per_frame = round(max_transition_rate * 12 * hop_length / sr)
transition_width = max_semitones_per_frame * n_bins_per_semitone + 1
# Construct the within voicing transition probabilities
transition = sequence.transition_local(
n_pitch_bins, transition_width, window="triangle", wrap=False
)
# Include across voicing transition probabilities
t_switch = sequence.transition_loop(2, 1 - switch_prob)
transition = np.kron(t_switch, transition)
p_init = np.zeros(2 * n_pitch_bins)
p_init[n_pitch_bins:] = 1 / n_pitch_bins
states = sequence.viterbi(observation_probs, transition, p_init=p_init)
# Find f0 corresponding to each decoded pitch bin.
freqs = fmin * 2 ** (np.arange(n_pitch_bins) / (12 * n_bins_per_semitone))
f0 = freqs[states % n_pitch_bins]
voiced_flag = states < n_pitch_bins
if fill_na is not None:
f0[~voiced_flag] = fill_na
return f0[..., 0, :], voiced_flag[..., 0, :], voiced_prob[..., 0, :]
def __pyin_helper(
yin_frames,
parabolic_shifts,
sr,
thresholds,
boltzmann_parameter,
beta_probs,
no_trough_prob,
min_period,
fmin,
n_pitch_bins,
n_bins_per_semitone,
):
yin_probs = np.zeros_like(yin_frames)
for i, yin_frame in enumerate(yin_frames.T):
# 2. For each frame find the troughs.
is_trough = util.localmin(yin_frame)
is_trough[0] = yin_frame[0] < yin_frame[1]
(trough_index,) = np.nonzero(is_trough)
if len(trough_index) == 0:
continue
# 3. Find the troughs below each threshold.
# these are the local minima of the frame, could get them directly without the trough index
trough_heights = yin_frame[trough_index]
trough_thresholds = np.less.outer(trough_heights, thresholds[1:])
# 4. Define the prior over the troughs.
# Smaller periods are weighted more.
trough_positions = np.cumsum(trough_thresholds, axis=0) - 1
n_troughs = np.count_nonzero(trough_thresholds, axis=0)
trough_prior = scipy.stats.boltzmann.pmf(
trough_positions, boltzmann_parameter, n_troughs
)
trough_prior[~trough_thresholds] = 0
# 5. For each threshold add probability to global minimum if no trough is below threshold,
# else add probability to each trough below threshold biased by prior.
probs = trough_prior.dot(beta_probs)
global_min = np.argmin(trough_heights)
n_thresholds_below_min = np.count_nonzero(~trough_thresholds[global_min, :])
probs[global_min] += no_trough_prob * np.sum(
beta_probs[:n_thresholds_below_min]
)
yin_probs[trough_index, i] = probs
yin_period, frame_index = np.nonzero(yin_probs)
# Refine peak by parabolic interpolation.
period_candidates = min_period + yin_period
period_candidates = period_candidates + parabolic_shifts[yin_period, frame_index]
f0_candidates = sr / period_candidates
# Find pitch bin corresponding to each f0 candidate.
bin_index = 12 * n_bins_per_semitone * np.log2(f0_candidates / fmin)
bin_index = np.clip(np.round(bin_index), 0, n_pitch_bins).astype(int)
# Observation probabilities.
observation_probs = np.zeros((2 * n_pitch_bins, yin_frames.shape[1]))
observation_probs[bin_index, frame_index] = yin_probs[yin_period, frame_index]
voiced_prob = np.clip(
np.sum(observation_probs[:n_pitch_bins, :], axis=0, keepdims=True), 0, 1
)
observation_probs[n_pitch_bins:, :] = (1 - voiced_prob) / n_pitch_bins
return observation_probs[np.newaxis], voiced_prob
def __check_yin_params(
*, sr: float, fmax: float, fmin: float, frame_length: int, win_length: int
):
"""Check the feasibility of yin/pyin parameters against
the following conditions:
1. 0 < fmin < fmax <= sr/2
2. frame_length - win_length - 1 > sr/fmax
"""
if fmax > sr / 2:
raise ParameterError(f"fmax={fmax:.3f} cannot exceed Nyquist frequency {sr/2}")
if fmin >= fmax:
raise ParameterError(f"fmin={fmin:.3f} must be less than fmax={fmax:.3f}")
if fmin <= 0:
raise ParameterError(f"fmin={fmin:.3f} must be strictly positive")
if win_length >= frame_length:
raise ParameterError(
f"win_length={win_length} must be less than frame_length={frame_length}"
)
if frame_length - win_length - 1 <= sr // fmax:
fmax_feasible = sr / (frame_length - win_length - 1)
frame_length_feasible = int(np.ceil(sr/fmax) + win_length + 1)
raise ParameterError(
f"fmax={fmax:.3f} is too small for frame_length={frame_length}, win_length={win_length}, and sr={sr}. "
f"Either increase to fmax={fmax_feasible:.3f} or frame_length={frame_length_feasible}"
)