librosa/core/spectrum.py
def pcen(S, sr=22050, hop_length=512, gain=0.98, bias=2, power=0.5,
         time_constant=0.400, eps=1e-6, b=None, max_size=1, ref=None,
         axis=-1, max_axis=None, zi=None, return_zf=False):
    '''Per-channel energy normalization (PCEN) [1]_

    This function normalizes a time-frequency representation `S` by
    performing automatic gain control, followed by nonlinear compression:

        P[f, t] = (S / (eps + M[f, t])**gain + bias)**power - bias**power

    where `M` is the result of applying a low-pass, temporal IIR filter
    to `S`:

        M[f, t] = (1 - b) * M[f, t - 1] + b * S[f, t]

    If `b` is not provided, it is calculated as:

        b = (sqrt(1 + 4* T**2) - 1) / (2 * T**2)

    where `T = time_constant * sr / hop_length`.

    This normalization is designed to suppress background noise and
    emphasize foreground signals, and can be used as an alternative to
    decibel scaling (`amplitude_to_db`).

    This implementation also supports smoothing across frequency bins
    by specifying `max_size > 1`.  If this option is used, the filtered
    spectrogram `M` is computed as

        M[f, t] = (1 - b) * M[f, t - 1] + b * R[f, t]

    where `R` has been max-filtered along the frequency axis, similar to
    the SuperFlux algorithm implemented in `onset.onset_strength`:

        R[f, t] = max(S[f - max_size//2: f + max_size//2, t])

    This can be used to perform automatic gain control on signals that cross
    or span multiple frequency bans, which may be desirable for spectrograms
    with high frequency resolution.

    .. [1] Wang, Y., Getreuer, P., Hughes, T., Lyon, R. F., & Saurous, R. A.
       (2017, March). Trainable frontend for robust and far-field keyword spotting.
       In Acoustics, Speech and Signal Processing (ICASSP), 2017
       IEEE International Conference on (pp. 5670-5674). IEEE.

    Parameters
    ----------
    S : np.ndarray (non-negative)
        The input (magnitude) spectrogram

    sr : number > 0 [scalar]
        The audio sampling rate

    hop_length : int > 0 [scalar]
        The hop length of `S`, expressed in samples

    gain : number >= 0 [scalar]
        The gain factor.  Typical values should be slightly less than 1.

    bias : number >= 0 [scalar]
        The bias point of the nonlinear compression (default: 2)

    power : number > 0 [scalar]
        The compression exponent.  Typical values should be between 0 and 1.
        Smaller values of `power` result in stronger compression.

    time_constant : number > 0 [scalar]
        The time constant for IIR filtering, measured in seconds.

    eps : number > 0 [scalar]
        A small constant used to ensure numerical stability of the filter.

    b : number in [0, 1]  [scalar]
        The filter coefficient for the low-pass filter.
        If not provided, it will be inferred from `time_constant`.

    max_size : int > 0 [scalar]
        The width of the max filter applied to the frequency axis.
        If left as `1`, no filtering is performed.

    ref : None or np.ndarray (shape=S.shape)
        An optional pre-computed reference spectrum (`R` in the above).
        If not provided it will be computed from `S`.

    axis : int [scalar]
        The (time) axis of the input spectrogram.

    max_axis : None or int [scalar]
        The frequency axis of the input spectrogram.
        If `None`, and `S` is two-dimensional, it will be inferred
        as the opposite from `axis`.
        If `S` is not two-dimensional, and `max_size > 1`, an error
        will be raised.

    zi : np.ndarray
        The initial filter delay values.

        This may be the `zf` (final delay values) of a previous call to `pcen`, or
        computed by `scipy.signal.lfilter_zi`.

    return_zf : bool
        If `True`, return the final filter delay values along with the PCEN output `P`.
        This is primarily useful in streaming contexts, where the final state of one
        block of processing should be used to initialize the next block.

        If `False` (default) only the PCEN values `P` are returned.


    Returns
    -------
    P : np.ndarray, non-negative [shape=(n, m)]
        The per-channel energy normalized version of `S`.

    zf : np.ndarray (optional)
        The final filter delay values.  Only returned if `return_zf=True`.

    See Also
    --------
    amplitude_to_db
    librosa.onset.onset_strength

    Examples
    --------

    Compare PCEN to log amplitude (dB) scaling on Mel spectra

    >>> import matplotlib.pyplot as plt
    >>> y, sr = librosa.load(librosa.util.example_audio_file(),
    ...                      offset=10, duration=10)

    >>> # We'll use power=1 to get a magnitude spectrum
    >>> # instead of a power spectrum
    >>> S = librosa.feature.melspectrogram(y, sr=sr, power=1)
    >>> log_S = librosa.amplitude_to_db(S, ref=np.max)
    >>> pcen_S = librosa.pcen(S)
    >>> plt.figure()
    >>> plt.subplot(2,1,1)
    >>> librosa.display.specshow(log_S, x_axis='time', y_axis='mel')
    >>> plt.title('log amplitude (dB)')
    >>> plt.colorbar()
    >>> plt.subplot(2,1,2)
    >>> librosa.display.specshow(pcen_S, x_axis='time', y_axis='mel')
    >>> plt.title('Per-channel energy normalization')
    >>> plt.colorbar()
    >>> plt.tight_layout()

    Compare PCEN with and without max-filtering

    >>> pcen_max = librosa.pcen(S, max_size=3)
    >>> plt.figure()
    >>> plt.subplot(2,1,1)
    >>> librosa.display.specshow(pcen_S, x_axis='time', y_axis='mel')
    >>> plt.title('Per-channel energy normalization (no max-filter)')
    >>> plt.colorbar()
    >>> plt.subplot(2,1,2)
    >>> librosa.display.specshow(pcen_max, x_axis='time', y_axis='mel')
    >>> plt.title('Per-channel energy normalization (max_size=3)')
    >>> plt.colorbar()
    >>> plt.tight_layout()

    '''

    if power <= 0:
        raise ParameterError('power={} must be strictly positive'.format(power))

    if gain < 0:
        raise ParameterError('gain={} must be non-negative'.format(gain))

    if bias < 0:
        raise ParameterError('bias={} must be non-negative'.format(bias))

    if eps <= 0:
        raise ParameterError('eps={} must be strictly positive'.format(eps))

    if time_constant <= 0:
        raise ParameterError('time_constant={} must be strictly positive'.format(time_constant))

    if max_size < 1 or not isinstance(max_size, int):
        raise ParameterError('max_size={} must be a positive integer'.format(max_size))

    if b is None:
        t_frames = time_constant * sr / float(hop_length)
        # By default, this solves the equation for b:
        #   b**2  + (1 - b) / t_frames  - 2 = 0
        # which approximates the full-width half-max of the
        # squared frequency response of the IIR low-pass filter

        b = (np.sqrt(1 + 4 * t_frames**2) - 1) / (2 * t_frames**2)

    if not 0 <= b <= 1:
        raise ParameterError('b={} must be between 0 and 1'.format(b))

    if np.issubdtype(S.dtype, np.complexfloating):
        warnings.warn('pcen was called on complex input so phase '
                      'information will be discarded. To suppress this warning, '
                      'call pcen(np.abs(D)) instead.')
        S = np.abs(S)

    if ref is None:
        if max_size == 1:
            ref = S
        elif S.ndim == 1:
            raise ParameterError('Max-filtering cannot be applied to 1-dimensional input')
        else:
            if max_axis is None:
                if S.ndim != 2:
                    raise ParameterError('Max-filtering a {:d}-dimensional spectrogram '
                                         'requires you to specify max_axis'.format(S.ndim))
                # if axis = 0, max_axis=1
                # if axis = +- 1, max_axis = 0
                max_axis = np.mod(1 - axis, 2)

            ref = scipy.ndimage.maximum_filter1d(S, max_size, axis=max_axis)

    if zi is None:
        # Make sure zi matches dimension to input
        shape = tuple([1] * ref.ndim)
        zi = np.empty(shape)
        zi[:] = scipy.signal.lfilter_zi([b], [1, b - 1])[:]

    S_smooth, zf = scipy.signal.lfilter([b], [1, b - 1], ref, zi=zi,
                                        axis=axis)

    # Working in log-space gives us some stability, and a slight speedup
    smooth = np.exp(-gain * (np.log(eps) + np.log1p(S_smooth / eps)))
    S_out = (S * smooth + bias)**power - bias**power

    if return_zf:
        return S_out, zf
    else:
        return S_out
def pcen(S, sr=22050, hop_length=512, gain=0.98, bias=2, power=0.5,
         time_constant=0.400, eps=1e-6, b=None, max_size=1, ref=None,
         axis=-1, max_axis=None, zi=None, return_zf=False):
    '''Per-channel energy normalization (PCEN) [1]_

    This function normalizes a time-frequency representation `S` by
    performing automatic gain control, followed by nonlinear compression:

        P[f, t] = (S / (eps + M[f, t])**gain + bias)**power - bias**power

    where `M` is the result of applying a low-pass, temporal IIR filter
    to `S`:

        M[f, t] = (1 - b) * M[f, t - 1] + b * S[f, t]

    If `b` is not provided, it is calculated as:

        b = (sqrt(1 + 4* T**2) - 1) / (2 * T**2)

    where `T = time_constant * sr / hop_length`.

    This normalization is designed to suppress background noise and
    emphasize foreground signals, and can be used as an alternative to
    decibel scaling (`amplitude_to_db`).

    This implementation also supports smoothing across frequency bins
    by specifying `max_size > 1`.  If this option is used, the filtered
    spectrogram `M` is computed as

        M[f, t] = (1 - b) * M[f, t - 1] + b * R[f, t]

    where `R` has been max-filtered along the frequency axis, similar to
    the SuperFlux algorithm implemented in `onset.onset_strength`:

        R[f, t] = max(S[f - max_size//2: f + max_size//2, t])

    This can be used to perform automatic gain control on signals that cross
    or span multiple frequency bans, which may be desirable for spectrograms
    with high frequency resolution.

    .. [1] Wang, Y., Getreuer, P., Hughes, T., Lyon, R. F., & Saurous, R. A.
       (2017, March). Trainable frontend for robust and far-field keyword spotting.
       In Acoustics, Speech and Signal Processing (ICASSP), 2017
       IEEE International Conference on (pp. 5670-5674). IEEE.

    Parameters
    ----------
    S : np.ndarray (non-negative)
        The input (magnitude) spectrogram

    sr : number > 0 [scalar]
        The audio sampling rate

    hop_length : int > 0 [scalar]
        The hop length of `S`, expressed in samples

    gain : number >= 0 [scalar]
        The gain factor.  Typical values should be slightly less than 1.

    bias : number >= 0 [scalar]
        The bias point of the nonlinear compression (default: 2)

    power : number > 0 [scalar]
        The compression exponent.  Typical values should be between 0 and 1.
        Smaller values of `power` result in stronger compression.

    time_constant : number > 0 [scalar]
        The time constant for IIR filtering, measured in seconds.

    eps : number > 0 [scalar]
        A small constant used to ensure numerical stability of the filter.

    b : number in [0, 1]  [scalar]
        The filter coefficient for the low-pass filter.
        If not provided, it will be inferred from `time_constant`.

    max_size : int > 0 [scalar]
        The width of the max filter applied to the frequency axis.
        If left as `1`, no filtering is performed.

    ref : None or np.ndarray (shape=S.shape)
        An optional pre-computed reference spectrum (`R` in the above).
        If not provided it will be computed from `S`.

    axis : int [scalar]
        The (time) axis of the input spectrogram.

    max_axis : None or int [scalar]
        The frequency axis of the input spectrogram.
        If `None`, and `S` is two-dimensional, it will be inferred
        as the opposite from `axis`.
        If `S` is not two-dimensional, and `max_size > 1`, an error
        will be raised.

    zi : np.ndarray
        The initial filter delay values.

        This may be the `zf` (final delay values) of a previous call to `pcen`, or
        computed by `scipy.signal.lfilter_zi`.

    return_zf : bool
        If `True`, return the final filter delay values along with the PCEN output `P`.
        This is primarily useful in streaming contexts, where the final state of one
        block of processing should be used to initialize the next block.

        If `False` (default) only the PCEN values `P` are returned.


    Returns
    -------
    P : np.ndarray, non-negative [shape=(n, m)]
        The per-channel energy normalized version of `S`.

    zf : np.ndarray (optional)
        The final filter delay values.  Only returned if `return_zf=True`.

    See Also
    --------
    amplitude_to_db
    librosa.onset.onset_strength

    Examples
    --------

    Compare PCEN to log amplitude (dB) scaling on Mel spectra

    >>> import matplotlib.pyplot as plt
    >>> y, sr = librosa.load(librosa.util.example_audio_file(),
    ...                      offset=10, duration=10)

    >>> # We'll use power=1 to get a magnitude spectrum
    >>> # instead of a power spectrum
    >>> S = librosa.feature.melspectrogram(y, sr=sr, power=1)
    >>> log_S = librosa.amplitude_to_db(S, ref=np.max)
    >>> pcen_S = librosa.pcen(S)
    >>> plt.figure()
    >>> plt.subplot(2,1,1)
    >>> librosa.display.specshow(log_S, x_axis='time', y_axis='mel')
    >>> plt.title('log amplitude (dB)')
    >>> plt.colorbar()
    >>> plt.subplot(2,1,2)
    >>> librosa.display.specshow(pcen_S, x_axis='time', y_axis='mel')
    >>> plt.title('Per-channel energy normalization')
    >>> plt.colorbar()
    >>> plt.tight_layout()

    Compare PCEN with and without max-filtering

    >>> pcen_max = librosa.pcen(S, max_size=3)
    >>> plt.figure()
    >>> plt.subplot(2,1,1)
    >>> librosa.display.specshow(pcen_S, x_axis='time', y_axis='mel')
    >>> plt.title('Per-channel energy normalization (no max-filter)')
    >>> plt.colorbar()
    >>> plt.subplot(2,1,2)
    >>> librosa.display.specshow(pcen_max, x_axis='time', y_axis='mel')
    >>> plt.title('Per-channel energy normalization (max_size=3)')
    >>> plt.colorbar()
    >>> plt.tight_layout()

    '''

    if power <= 0:
        raise ParameterError('power={} must be strictly positive'.format(power))

    if gain < 0:
        raise ParameterError('gain={} must be non-negative'.format(gain))

    if bias < 0:
        raise ParameterError('bias={} must be non-negative'.format(bias))

    if eps <= 0:
        raise ParameterError('eps={} must be strictly positive'.format(eps))

    if time_constant <= 0:
        raise ParameterError('time_constant={} must be strictly positive'.format(time_constant))

    if max_size < 1 or not isinstance(max_size, int):
        raise ParameterError('max_size={} must be a positive integer'.format(max_size))

    if b is None:
        t_frames = time_constant * sr / float(hop_length)
        # By default, this solves the equation for b:
        #   b**2  + (1 - b) / t_frames  - 2 = 0
        # which approximates the full-width half-max of the
        # squared frequency response of the IIR low-pass filter

        b = (np.sqrt(1 + 4 * t_frames**2) - 1) / (2 * t_frames**2)

    if not 0 <= b <= 1:
        raise ParameterError('b={} must be between 0 and 1'.format(b))

    if np.issubdtype(S.dtype, np.complexfloating):
        warnings.warn('pcen was called on complex input so phase '
                      'information will be discarded. To suppress this warning, '
                      'call pcen(np.abs(D)) instead.')
        S = np.abs(S)

    if ref is None:
        if max_size == 1:
            ref = S
        elif S.ndim == 1:
            raise ParameterError('Max-filtering cannot be applied to 1-dimensional input')
        else:
            if max_axis is None:
                if S.ndim != 2:
                    raise ParameterError('Max-filtering a {:d}-dimensional spectrogram '
                                         'requires you to specify max_axis'.format(S.ndim))
                # if axis = 0, max_axis=1
                # if axis = +- 1, max_axis = 0
                max_axis = np.mod(1 - axis, 2)

            ref = scipy.ndimage.maximum_filter1d(S, max_size, axis=max_axis)

    if zi is None:
        # Make sure zi matches dimension to input
        shape = tuple([1] * ref.ndim)
        zi = np.empty(shape)
        zi[:] = scipy.signal.lfilter_zi([b], [1, b - 1])[:]

    S_smooth, zf = scipy.signal.lfilter([b], [1, b - 1], ref, zi=zi,
                                        axis=axis)

    # Working in log-space gives us some stability, and a slight speedup
    smooth = np.exp(-gain * (np.log(eps) + np.log1p(S_smooth / eps)))
    S_out = (S * smooth + bias)**power - bias**power

    if return_zf:
        return S_out, zf
    else:
        return S_out
Per-channel energy normalization (PCEN) [1]_

    This function normalizes a time-frequency representation `S` by
    performing automatic gain control, followed by nonlinear compression:

        P[f, t] = (S / (eps + M[f, t])**gain + bias)**power - bias**power

    where `M` is the result of applying a low-pass, temporal IIR filter
    to `S`:

        M[f, t] = (1 - b) * M[f, t - 1] + b * S[f, t]

    If `b` is not provided, it is calculated as:

        b = (sqrt(1 + 4* T**2) - 1) / (2 * T**2)

    where `T = time_constant * sr / hop_length`.

    This normalization is designed to suppress background noise and
    emphasize foreground signals, and can be used as an alternative to
    decibel scaling (`amplitude_to_db`).

    This implementation also supports smoothing across frequency bins
    by specifying `max_size > 1`.  If this option is used, the filtered
    spectrogram `M` is computed as

        M[f, t] = (1 - b) * M[f, t - 1] + b * R[f, t]

    where `R` has been max-filtered along the frequency axis, similar to
    the SuperFlux algorithm implemented in `onset.onset_strength`:

        R[f, t] = max(S[f - max_size//2: f + max_size//2, t])

    This can be used to perform automatic gain control on signals that cross
    or span multiple frequency bans, which may be desirable for spectrograms
    with high frequency resolution.

    .. [1] Wang, Y., Getreuer, P., Hughes, T., Lyon, R. F., & Saurous, R. A.
       (2017, March). Trainable frontend for robust and far-field keyword spotting.
       In Acoustics, Speech and Signal Processing (ICASSP), 2017
       IEEE International Conference on (pp. 5670-5674). IEEE.

    Parameters
    ----------
    S : np.ndarray (non-negative)
        The input (magnitude) spectrogram

    sr : number > 0 [scalar]
        The audio sampling rate

    hop_length : int > 0 [scalar]
        The hop length of `S`, expressed in samples

    gain : number >= 0 [scalar]
        The gain factor.  Typical values should be slightly less than 1.

    bias : number >= 0 [scalar]
        The bias point of the nonlinear compression (default: 2)

    power : number > 0 [scalar]
        The compression exponent.  Typical values should be between 0 and 1.
        Smaller values of `power` result in stronger compression.

    time_constant : number > 0 [scalar]
        The time constant for IIR filtering, measured in seconds.

    eps : number > 0 [scalar]
        A small constant used to ensure numerical stability of the filter.

    b : number in [0, 1]  [scalar]
        The filter coefficient for the low-pass filter.
        If not provided, it will be inferred from `time_constant`.

    max_size : int > 0 [scalar]
        The width of the max filter applied to the frequency axis.
        If left as `1`, no filtering is performed.

    ref : None or np.ndarray (shape=S.shape)
        An optional pre-computed reference spectrum (`R` in the above).
        If not provided it will be computed from `S`.

    axis : int [scalar]
        The (time) axis of the input spectrogram.

    max_axis : None or int [scalar]
        The frequency axis of the input spectrogram.
        If `None`, and `S` is two-dimensional, it will be inferred
        as the opposite from `axis`.
        If `S` is not two-dimensional, and `max_size > 1`, an error
        will be raised.

    zi : np.ndarray
        The initial filter delay values.

        This may be the `zf` (final delay values) of a previous call to `pcen`, or
        computed by `scipy.signal.lfilter_zi`.

    return_zf : bool
        If `True`, return the final filter delay values along with the PCEN output `P`.
        This is primarily useful in streaming contexts, where the final state of one
        block of processing should be used to initialize the next block.

        If `False` (default) only the PCEN values `P` are returned.


    Returns
    -------
    P : np.ndarray, non-negative [shape=(n, m)]
        The per-channel energy normalized version of `S`.

    zf : np.ndarray (optional)
        The final filter delay values.  Only returned if `return_zf=True`.

    See Also
    --------
    amplitude_to_db
    librosa.onset.onset_strength

    Examples
    --------

    Compare PCEN to log amplitude (dB) scaling on Mel spectra

    >>> import matplotlib.pyplot as plt
    >>> y, sr = librosa.load(librosa.util.example_audio_file(),
    ...                      offset=10, duration=10)

    >>> # We'll use power=1 to get a magnitude spectrum
    >>> # instead of a power spectrum
    >>> S = librosa.feature.melspectrogram(y, sr=sr, power=1)
    >>> log_S = librosa.amplitude_to_db(S, ref=np.max)
    >>> pcen_S = librosa.pcen(S)
    >>> plt.figure()
    >>> plt.subplot(2,1,1)
    >>> librosa.display.specshow(log_S, x_axis='time', y_axis='mel')
    >>> plt.title('log amplitude (dB)')
    >>> plt.colorbar()
    >>> plt.subplot(2,1,2)
    >>> librosa.display.specshow(pcen_S, x_axis='time', y_axis='mel')
    >>> plt.title('Per-channel energy normalization')
    >>> plt.colorbar()
    >>> plt.tight_layout()

    Compare PCEN with and without max-filtering

    >>> pcen_max = librosa.pcen(S, max_size=3)
    >>> plt.figure()
    >>> plt.subplot(2,1,1)
    >>> librosa.display.specshow(pcen_S, x_axis='time', y_axis='mel')
    >>> plt.title('Per-channel energy normalization (no max-filter)')
    >>> plt.colorbar()
    >>> plt.subplot(2,1,2)
    >>> librosa.display.specshow(pcen_max, x_axis='time', y_axis='mel')
    >>> plt.title('Per-channel energy normalization (max_size=3)')
    >>> plt.colorbar()
    >>> plt.tight_layout()