Repository: pyfar/pyfar
Pull request #325 (2022-06-28 20:06:10)

Soft limiting
### Changes proposed in this pull request:

Add a function for soft-limiting Signal and FrequencyData Objects. Draft until discussed and tests are finished.


Patch:
diff --git a/pyfar/dsp/__init__.py b/pyfar/dsp/__init__.py
index 06832872..ea91af01 100644
--- a/pyfar/dsp/__init__.py
+++ b/pyfar/dsp/__init__.py
@@ -16,6 +16,7 @@ from .dsp import (
     deconvolve,
     convolve,
     decibel,
+    soft_limit_spectrum,
     energy,
     power,
     rms,
@@ -54,6 +55,7 @@ __all__ = [
     'deconvolve',
     'convolve',
     'decibel',
+    'soft_limit_spectrum',
     'energy',
     'power',
     'rms',
diff --git a/pyfar/dsp/dsp.py b/pyfar/dsp/dsp.py
index a586a5a7..1881b4a6 100644
--- a/pyfar/dsp/dsp.py
+++ b/pyfar/dsp/dsp.py
@@ -1797,12 +1797,262 @@ def decibel(signal, domain='freq', log_prefix=None, log_reference=1,
             f"Domain is '{domain}', but has to be 'time', 'freq',"
             " or 'freq_raw'.")
     data[data == 0] = np.finfo(float).eps
+
     if return_prefix is True:
         return log_prefix * np.log10(np.abs(data) / log_reference), log_prefix
     else:
         return log_prefix * np.log10(np.abs(data) / log_reference)
 
 
+def soft_limit_spectrum(signal, limit, knee, frequency_range=None,
+                        direction='upper', log_prefix=None):
+    """
+    Soft limiting the magniude spectrum.
+
+    Soft limiting gradually increases the gain reduction to avoid
+    discontinuities in the data that would appear in hard limiting. The
+    transition between the magnitude where no limiting is applied to the
+    magnitude where the limiting reaches its full effect is termed `knee`
+    (see examples below).
+
+    Note that the limiting is applied on `signal.freq`, i.e. the data after
+    the FFT normalization.
+
+    Parameters
+    ----------
+    signal : Signal, FrequencyData
+        The input data
+    limit : number, array like
+        The gain in dB at which the limiting reaches its full effect. If this
+        is a number, the same limit is applied to all frequencies. If this an
+        array like, it must be broadcastable to ``signal.freq`` to apply
+        frequency-dependent limits.
+    knee : number, string
+        If this is a number, a knee with a width of `number` dB according to
+        [#]_ Eq. (4) is applied. This definition of the knee originates from
+        the classic limiting audio effect. If this is ``'arctan'`` an arcus
+        tangens knee according to [#]_ Section 3.6.4 is applied. This knee
+        definition originates from microphone array signal processing.
+    frequency_range : array like, optional
+        Frequency range in which the limiting is applied. This must be an array
+        like containing the lower and upper limit in Hz. The default ``None``
+        applies the limiting to all frequencies.
+    direction : str, optional
+        Define how the limiting works
+
+        ``'upper'`` (default)
+            Soft limiting `signal` to enforce an aboslute maximum value of
+            `limit`.
+        ``'lower'``
+            Soft limiting `signal` to enforce an absolute minimum value of
+            `limit`
+    log_prefix: float, int
+        The log prefix is used to linearize the limit and knee, e.g.,
+        ``limit_linear = 10**(limit / log_prefix)``.The default ``None``, uses
+        ``10`` for signals with ``'psd'`` and ``'power'`` FFT normalization and
+        ``20`` otherwise.
+
+    Returns
+    -------
+    limited : Signal, FrequencyData
+        The limited copy of the input data.
+
+    Examples
+    --------
+
+    Illustrate effect of limit and knee
+
+    .. plot ::
+
+        >>> import pyfar as pf
+        >>> import numpy as np
+        >>>
+        >>> signal = pf.FrequencyData(
+        ...     10**(np.arange(-20, 21)/20), np.arange(0, 41))
+        >>> for knee in ['arctan', 20, 0]:
+        >>>     limited = pf.dsp.soft_limit_spectrum(signal, 0, knee)
+        >>>     pf.plot.freq(limited, freq_scale='linear', label=f'{knee=}')
+        >>> ax = pf.plot.freq(signal, freq_scale='linear',
+        ...                   linestyle='--', label='input')
+        >>>
+        >>> ax.legend(loc='upper left')
+        >>> ax.set_xlim(0, 40)
+        >>> ax.set_ylim(-20, 20)
+        >>> bbox=dict(facecolor='white', edgecolor='none', alpha=0.5)
+        >>> arrowprops= dict(
+        ...     arrowstyle="<->", shrinkA=0, shrinkB=0,
+        ...     color=pf.plot.color('r'))
+        >>> ax.annotate("limit=0 dB", (20, 0), (10, 0), bbox=bbox,
+        ...             arrowprops=dict(arrowstyle="->"), va='center')
+        >>> ax.text(20, -14, 'knee=20 dB', ha='center',
+        ...         color=pf.plot.color('r'), bbox=bbox)
+        >>> ax.annotate("", (10, -15), (30, -15), arrowprops=arrowprops)
+
+    Apply limiting with knee (soft limiting) and without knee (hard limiting)
+
+    .. plot ::
+
+        >>> import pyfar as pf
+        >>>
+        >>> signal = pf.signals.impulse(4096)
+        >>> signal = pf.dsp.filter.bell(signal, 1e3, 20, 1)
+        >>>
+        >>> soft = pf.dsp.soft_limit_spectrum(signal, limit=10, knee=12)
+        >>> hard = pf.dsp.soft_limit_spectrum(signal, limit=10, knee=0)
+        >>>
+        >>> pf.plot.freq(soft, label='soft (knee > 0)')
+        >>> pf.plot.freq(hard, label='hard (knee = 0)')
+        >>> ax = pf.plot.freq(signal, ls='--', label='original')
+        >>>
+        >>> ax.set_ylim(-5, 25)
+        >>> ax.legend()
+
+    Apply soft limiting above and below +/-10 dB
+
+    .. plot ::
+
+        >>> import pyfar as pf
+        >>>
+        >>> signal = pf.signals.impulse(4096)
+        >>> signal = pf.dsp.filter.bell(signal, 2e2, -20, 3)
+        >>> signal = pf.dsp.filter.bell(signal, 6e3, 20, 3)
+        >>>
+        >>> upper = pf.dsp.soft_limit_spectrum(
+        ...     signal, limit=10, knee=6, direction="upper")
+        >>> lower = pf.dsp.soft_limit_spectrum(
+        ...     signal, limit=-10, knee=6, direction="lower")
+        >>>
+        >>> pf.plot.freq(upper, label='upper limit')
+        >>> pf.plot.freq(lower, label='lower limit')
+        >>> ax = pf.plot.freq(signal, ls='--', label='original')
+        >>>
+        >>> ax.set_ylim(-25, 25)
+        >>> ax.legend()
+
+    Apply frequency dependent soft limiting
+
+    .. plot ::
+
+        >>> import pyfar as pf
+        >>>
+        >>> # input data
+        >>> signal = pf.signals.impulse(2**13)
+        >>> signal = pf.dsp.filter.bell(signal, 2e2, 20, 4)
+        >>> signal = pf.dsp.filter.bell(signal, 2e3, 20, 4)
+        >>>
+        >>> # frequency dependent limit
+        >>> limit = pf.dsp.filter.low_shelf(
+        ...     pf.signals.impulse(2**13), 500, -10, 2) * 10
+        >>> limit = pf.dsp.decibel(limit).flatten()
+        >>>
+        >>> # soft limiting input signal
+        >>> limited = pf.dsp.soft_limit_spectrum(signal, limit, 10)
+        >>>
+        >>> # plot
+        >>> ax = pf.plot.freq(signal, label='original')
+        >>> pf.plot.freq(limited, label='limited')
+        >>> ax.plot(signal.frequencies, limit, label='limit', c='k', ls='--')
+        >>> ax.set_ylim(-5, 25)
+        >>> ax.legend(loc='upper left')
+
+    References
+    -----------
+    .. [#] B. Bernschütz, Microphone arrays and sound field decomposition for
+           dynamic binaural synthesis, Ph.D Thesis, (Berlin, Germany,
+           TU Berlin, 2016).
+
+    .. [#]  D. Giannoulis, M. Massberg, and J. D. Reiss, Digital Dynamic Range
+            Compressor Design—A Tutorial and Analysis, J. Audio Eng. Soc. 60,
+            399-408 (2012).
+    """
+
+    # check input
+    if not isinstance(signal, (pyfar.Signal, pyfar.FrequencyData)):
+        raise TypeError(
+            "input signal must be a pyfar Signal or FrequencyData object")
+
+    if direction not in ['upper', 'lower']:
+        raise ValueError((f"direction is '{direction}' but must be 'upper', "
+                          "or 'lower'"))
+
+    if isinstance(knee, str) and knee != "arctan":
+        raise ValueError(f"knee is '{knee}' but 'arctan'")
+    elif isinstance(knee, (int, float)) and knee < 0:
+        raise ValueError(f"knee is {knee} but must be >= 0")
+    elif not isinstance(knee, (str, int, float)):
+        raise TypeError("knee must be a string or number")
+
+    limit = np.broadcast_to(limit, signal.cshape + (signal.n_bins, ))
+
+    # define frequency range
+    if frequency_range is None:
+        freq_mask = np.full((signal.n_bins, ), True)
+    else:
+        freq_mask = np.full(signal.n_bins, False)
+        freq_mask[signal.find_nearest_frequency(np.min(frequency_range)):
+                  signal.find_nearest_frequency(np.max(frequency_range)) + 1] \
+            = True
+
+    # get spectral data
+    signal_limited = signal.copy()
+    freq = signal_limited.freq
+
+    # handle 'inverse' limiting
+    if direction == 'lower':
+        freq[freq == 0] = np.finfo(float).eps
+        freq = 1 / freq
+
+        limit = -limit
+
+    # de-logarithmize the limit
+    if log_prefix is None:
+        if (type(signal) is pyfar.Signal and
+                signal.fft_norm not in ('power', 'psd')) or \
+                type(signal) is pyfar.FrequencyData:
+            log_prefix = 20
+        else:
+            log_prefix = 10
+
+    limit_lin = 10**(limit/log_prefix)
+
+    # absolute spectrum
+    freq_abs = np.abs(freq)
+
+    # arcus tangens limiting
+    if knee == "arctan":
+        alpha = freq_abs[..., freq_mask] / limit_lin[..., freq_mask]
+        freq[..., freq_mask] *= \
+            2/np.pi / alpha * np.arctan(np.pi/2 * alpha)
+    # classic audio effect limiting
+    else:
+
+        # log spectrum
+        freq_db = log_prefix * np.log10(freq_abs)
+
+        # hard limiting outside the knee (Giannoulis et al. Eq. (4)) ----------
+        hard = np.logical_and(2 * (freq_db - limit) > knee, freq_mask)
+        gain_hard = limit_lin[..., hard]
+        freq[..., hard] *= gain_hard / freq_abs[..., hard]
+
+        # soft limiting inside the knee (Giannoulis et al. Eq. (4)) ----------
+        if knee != 0:
+
+            # frequencies, where limiting is applied
+            soft = np.logical_and(
+                2 * np.abs(freq_db - limit) <= knee, freq_mask)
+            # gain factor
+            gain_soft = -(freq_db[..., soft] - limit[..., soft] + knee / 2)**2\
+                / (2 * knee)
+            # apply limiting
+            freq[..., soft] *= 10**(gain_soft / log_prefix)
+
+    if direction == "lower":
+        freq = 1 / freq
+
+    signal_limited.freq = freq
+    return signal_limited
+
+
 def energy(signal):
     r"""
     Computes the channel wise energy in the time domain


Test patch:
diff --git a/tests/test_dsp_soft_limit_spectrum.py b/tests/test_dsp_soft_limit_spectrum.py
new file mode 100644
index 00000000..a52bd64f
--- /dev/null
+++ b/tests/test_dsp_soft_limit_spectrum.py
@@ -0,0 +1,156 @@
+from pyfar.dsp import soft_limit_spectrum
+import pyfar as pf
+import numpy as np
+import numpy.testing as npt
+import pytest
+
+
+def test_assertions_signal():
+    """Test assertion for passing wrong type of audio object."""
+
+    with pytest.raises(TypeError, match='input signal must be a pyfar'):
+        soft_limit_spectrum(pf.TimeData([1, 2, 3], [0, 1, 3]), 0, 0)
+
+
+def test_assertion_direction():
+    """Test assertion for passing invalid value for `direction` parameter."""
+
+    with pytest.raises(ValueError, match="direction is 'both'"):
+        soft_limit_spectrum(pf.Signal([1, 2, 3], 1), 0, 0, direction='both')
+
+
+def test_assertion_knee():
+    """Test assertion for passing invalid value for `knee` parameter."""
+
+    # wrong string
+    with pytest.raises(ValueError, match="knee is 'tangens'"):
+        soft_limit_spectrum(pf.Signal([1, 2, 3], 1), 0, knee='tangens')
+
+    # wrong number
+    with pytest.raises(ValueError, match="knee is -1"):
+        soft_limit_spectrum(pf.Signal([1, 2, 3], 1), 0, knee=-1)
+
+    # wrong type
+    with pytest.raises(TypeError, match="knee must be"):
+        soft_limit_spectrum(pf.Signal([1, 2, 3], 1), 0, knee=(1, 1))
+
+
+@pytest.mark.parametrize('data_in', [
+    pf.Signal([1, 0, 0], 1), pf.FrequencyData([1, 1, 1], [0, 1, 3])])
+def test_signal_and_frequency_data_input(data_in):
+    """Test with all possible audio objects as input."""
+
+    data_out = soft_limit_spectrum(data_in, limit=0, knee=0)
+    assert type(data_out) is type(data_in)
+
+
+@pytest.mark.parametrize('limit', [-10, 0, 10])
+@pytest.mark.parametrize('direction', ['upper', 'lower'])
+def test_limit_and_direction(limit, direction):
+    """Test if all values are correctly limited"""
+
+    data_in = pf.FrequencyData(10**(np.arange(-100, 101)/20), np.arange(201))
+    data_out = soft_limit_spectrum(data_in, limit, knee=0, direction=direction)
+
+    if direction == 'upper':
+        assert np.all(20*np.log10(np.abs(data_out.freq)) <= limit + 1e-14)
+    else:
+        assert np.all(20*np.log10(np.abs(data_out.freq)) >= limit - 1e-14)
+
+
+def test_frequency_dependend_limit():
+    """Test frequency-dependent limit passed as audio object"""
+
+    limit = np.atleast_2d([0, -6, -12])
+    data_in = pf.FrequencyData(10**(np.array([-3, -3, -3])/20), [0, 1, 3])
+    data_out = soft_limit_spectrum(data_in, limit, knee=0)
+
+    limit_applied = 20*np.log10(np.abs(data_in.freq)) > limit
+
+    # test limited values
+    npt.assert_almost_equal(20*np.log10(np.abs(data_out.freq[limit_applied])),
+                            limit[limit_applied])
+    # test non-limited values
+    npt.assert_equal(data_out.freq[~limit_applied],
+                     data_in.freq[~limit_applied])
+
+
+@pytest.mark.parametrize('limit', [-10, 0, 10])
+@pytest.mark.parametrize('knee', [0, 5, 10])
+def test_knee_width_in_db(limit, knee):
+    "Test the knee width given in decibel."""
+
+    data_in = pf.FrequencyData(10**(np.arange(-100, 101)/20), np.arange(201))
+    data_out = soft_limit_spectrum(data_in, limit, knee)
+
+    freq_in = 20*np.log10(np.abs(data_in.freq))
+    freq_out = 20*np.log10(np.abs(data_out.freq))
+
+    # categorize data
+    lower_knee = freq_in <= (limit - knee / 2)
+    upper_knee = freq_in >= (limit + knee / 2)
+    within_knee = np.logical_not(np.logical_or(lower_knee, upper_knee))
+
+    npt.assert_almost_equal(freq_in[lower_knee], freq_out[lower_knee], 14)
+    assert np.all(freq_in[within_knee] > freq_out[within_knee])
+    assert np.all(freq_out[upper_knee] <= limit + 1e-14)
+
+
+@pytest.mark.parametrize('limit', [-10, 0, 10])
+def test_arctan_knee(limit):
+    """Test the arcus tangens knee."""
+
+    data_in = pf.FrequencyData(10**(np.arange(-100, 101)/20), np.arange(201))
+    data_out = soft_limit_spectrum(data_in, limit, 'arctan')
+
+    # arcus tangens knee does at least a minimal gain reduction everywhere
+    assert np.all(data_out.freq < data_in.freq)
+    assert np.all(20*np.log10(np.abs(data_out.freq)) <= limit + 1e-14)
+
+
+@pytest.mark.parametrize('fft_norm,log_prefix', [
+    ('none', 20), ('unitary', 20), ('amplitude', 20), ('rms', 20),
+    ('power', 10), ('psd', 10)])
+def test_automatic_setting_of_log_prefix(fft_norm, log_prefix):
+    """
+    Test if the `log_prefix` is correctly set depending on the FFT
+    normalization.
+    """
+
+    # results in linearized limits of
+    # - 100 if log_prefix = 10
+    # -  10 if log_prefix = 20
+    limit = 20
+
+    data_in = pf.Signal([1, 0, 0], 44100)
+    data_in.fft_norm = fft_norm
+
+    # directly set data.freq which is used inside soft_limit_spectrum
+    # first value only limited if log_prefix = 20
+    # second value always limited
+    data_in.freq = [20, 200]
+
+    data_out = soft_limit_spectrum(data_in, limit, knee=0)
+
+    if log_prefix == 20:
+        npt.assert_almost_equal(data_out.freq.flatten(), [10, 10])
+    else:
+        npt.assert_almost_equal(data_out.freq.flatten(), [20, 100])
+
+
+def test_frequency_range():
+    """Test limiting inside and outside user specified frequency range."""
+
+    data_in = pf.FrequencyData([2, 2, 2, 2], [0, 1, 2, 3])
+    data_out = soft_limit_spectrum(data_in, 0, knee=0, frequency_range=[1, 2])
+    npt.assert_almost_equal(data_out.freq.flatten(), [2, 1, 1, 2])
+
+
+def test_custom_log_prefix():
+    """Test user specified `log_prefix`."""
+
+    data_in = pf.FrequencyData([5, 15], [0, 1])
+    # results in linearized limit of 10
+    data_out = soft_limit_spectrum(data_in, 1, knee=0, log_prefix=1)
+
+    npt.assert_almost_equal(data_out.freq.flatten(), [5, 10])