Repository: pyfar/pyfar (python)

Issue:
Restructure `pf.dsp.regularized_spectrum_inversion`
## General

* pyfar version: 0.5.2

## Description

We discussed that function in its current form is a bit too complicated. It has a special case of a regularization function implement through `freq_range`, etc. and a general case through `regu_final`. We decided to refactor the function on the long run:

- move the regularization function generated by `freq_range` et al. to a separate function/class
- always require to pass a regularization function to `pf.dsp.regularized_spectrum_inversion`
- regularization function must be passed as `FrequencyData`object

Pull request description:
Regularization class
New class for regularization and spectrum inversion to replace `dsp.regularized_spectrum_inversion()`.
Each method of defining a regularization function has it's own class method, so more regularization types can easily be added.

Closes #440 #635 

### Changes proposed in this pull request:

- move spectrum inversion to `Regularization.invert()` and deprecate dsp function
- add possibility to pass target function for inversion. If provided the convolution of the inverse filter with the system approaches the target function.
- simplified regularization from frequency range by replacing regu_inside and regu_outside by a single regularization factor. To this end regu_inside is set to 0, and the `normalize` parameter is removed. Results using the new parameters are close to identical but the interface is cleaner (see example below).
- regularization from a passed signal, containing the regularization factors as `FrequencyData` (_regu_final_ in dsp-function)

### Previous vs. new implementation

The regularization function is now squared, which was not done previously but seems to be more common. This causes slight differences, with respect to the old method:

```python
import pyfar as pf
%matplotlib inline

sweep = pf.signals.exponential_sweep_freq(2**16, [50, 16e3], 1000, 10e3)

inversion_1 = pf.dsp.regularized_spectrum_inversion(sweep, [50, 16e3])

Inversion = pf.dsp.RegularizedSpectrumInversion.from_frequency_range(
    [50, 16e3])
inversion_2 = Inversion.invert(sweep, normalize_regularization='max')

ax = pf.plot.freq(inversion_1, label='previous')
pf.plot.freq(inversion_2, linestyle=':', label='new')
ax.legend()
```

![output](https://github.com/user-attachments/assets/7ca66562-853e-466d-b8fb-4b660d1b346d)

### tbc
- Add more regularization types (e.g. maximum gain, automatic regularization parameter)

### For discussion:
- [x] Rename class to `RegularizedSpectrumInversion`, `from_signal` to `from_magnitude_spectrum`, and `get_regularization` to `regularization`
- [x] Use magnitude and phase of target function
- [x] Make the normalization of the regularization function a parameter
- [x] Check warnings for using `len(Signal)` when using `if signal`
- [ ] Examples Notebook

Patch:
diff --git a/pyfar/dsp/__init__.py b/pyfar/dsp/__init__.py
index 2acbf014..de01d971 100644
--- a/pyfar/dsp/__init__.py
+++ b/pyfar/dsp/__init__.py
@@ -37,6 +37,8 @@ from .interpolation import (
 from . import filter
 from . import fft
 
+from .regularization import RegularizedSpectrumInversion
+
 
 __all__ = [
     'fft',
@@ -69,4 +71,5 @@ __all__ = [
     'normalize',
     'fractional_time_shift',
     'correlate',
+    'RegularizedSpectrumInversion',
 ]
diff --git a/pyfar/dsp/regularization.py b/pyfar/dsp/regularization.py
new file mode 100644
index 00000000..95bf3255
--- /dev/null
+++ b/pyfar/dsp/regularization.py
@@ -0,0 +1,428 @@
+"""
+Documentation for regularization class.
+"""
+
+import pyfar as pf
+from pyfar.dsp.dsp import _cross_fade
+from pyfar.dsp.fft import _n_samples_from_n_bins
+import numpy as np
+import numbers
+
+class RegularizedSpectrumInversion():
+    r"""Class for frequency-dependent regularized inversion.
+
+    Regularized inversion limits the gain applied by an inverse filter.
+    This can be useful to avoid extreme amplification when inverting a signal
+    of limited bandwidth or containing notches in the spectrum.
+
+    The inverse is computed as [#]_:
+
+    .. math::
+        :label: regularized_inversion
+
+        S^{-1}(f) = \frac{S^*(f)}{S^*(f)S(f) + \beta |\epsilon(f)|^2} D(f)
+
+    with :math:`f` being the frequency, :math:`S(f)` the spectrum of the signal
+    to be inverted, :math:`^*` the complex conjugate,
+    :math:`\epsilon(f)` the regularization, and :math:`\beta` a scalar to
+    control the amount of regularization. :math:`D(f)` denotes an optional
+    target function.
+
+    The compensated system :math:`C = S(f)S^{-1}(f)` approaches the target
+    function in magnitude and phase. In many applications, the target function
+    should contain a delay to make sure that :math:`S^{-1}(f)` is causal. The
+    larger :math:`\beta` and :math:`\epsilon(f)` are, the larger the deviation
+    of :math:`C(f)` from the target :math:`D(f)`.
+
+    The inversion is done in two steps:
+
+    1. Define :math:`S(f)`, :math:`\epsilon(f)`, :math:`\beta`, and
+       :math:`D(f)` using one of the ``from_...()`` methods listed below.
+    2. Compute the inverse :math:`S(f)^{-1}` using :py:func:`~invert`
+
+    The parameters that are defined in the first step are often iteratively
+    adjusted. Examples are given in the documentation of the ``from_...()``
+    methods.
+
+    References
+    ----------
+    .. [#]  S. G. Norcross, M. Bouchard, and G. A. Soulodre, “Inverse Filtering
+            Design Using a Minimal-Phase Target Function from Regularization,”
+            Convention Paper 6929 Presented at the 121st Convention, 2006
+            October 5–8, San Francisco, CA, USA.
+    """
+
+    def __init__(self) -> None:
+
+        # throw error if object is instanced without classmethod
+        raise RuntimeError("Regularization objects must be created "
+                           "using one of the 'from_()' classmethods.")
+
+    def __repr__(self):
+        """String representation of RegularizedSpectrumInversion class."""
+        return f"Regularization object with regularization "\
+            f"of type '{self._regularization_type}'."
+
+    @classmethod
+    def from_frequency_range(cls, signal, frequency_range,
+                             regularization_within=0, beta=1, target=None):
+        r"""
+        Regularization from a given frequency range.
+
+        Defines a frequency range within which the regularization
+        :math:`\epsilon(f)` is set to `regularization_within`.
+        Outside the frequency range the regularization is
+        :math:`\epsilon(f)=1` and can be controlled using the `beta`
+        parameter.
+        The regularization factors are cross-faded using a raised cosine window
+        function with a width of :math:`\sqrt{2}f` above and below the given
+        frequency range.
+
+        Parameters
+        ----------
+        signal : Signal
+            Signal to be inverted.
+        frequency_range : array like
+            Array like containing the lower and upper frequency limit in Hz.
+        regularization_within: float, optional
+            Set regularization inside frequency range. The default is `0`.
+        beta : float, string, optional
+            Beta parameter to control the scaling of the regularization as in
+            equation :eq:`regularized_inversion`. Can be a
+
+            ``numerical value``
+                Usually between ``0`` and ``1``, with ``0`` being no
+                regularization applied.
+            ``'energy'``
+                Normalize the regularization to match the signal's energy.
+
+                :math:`\beta = \frac{\frac{1}{N}\sum_{k=0}^{N-1}|S[k]|^2}{\frac{1}{N}\sum_{k=0}^{N-1}|\epsilon[k]|^2}`
+            ``'max'``
+                Normalize the regularization :math:`\epsilon(f)` to the maximum
+                magnitude of a given signal :math:`S(f)`.
+
+                :math:`\beta = \frac{\max(|S(f)|)}{\max(|\epsilon(f)|)}`
+            ``'mean'``
+                Normalize the regularization :math:`\epsilon(f)` to the mean
+                magnitude of a given signal :math:`S(f).`
+
+                :math:`\beta = \frac{\text{mean}(|S(f)|)}{\text{mean}(|\epsilon(f)|)}`
+
+            The default is ``1``.
+        target : Signal, optional
+            Target function for the regularization. The default ``None`` uses a
+            zero-phase spectrum with an amplitude of 1 as target, equal to an
+            impulse in the time domain.
+
+
+        Examples
+        --------
+        Invert a sine sweep with limited bandwidth and apply maximum
+        normalization to the regularization function.
+
+        .. plot::
+
+            >>> import pyfar as pf
+            >>> import matplotlib.pyplot as plt
+            ...
+            >>> sweep = pf.signals.exponential_sweep_freq(
+            ...     2**16, [50, 16e3], 1000, 10e3)
+            ...
+            >>> # Inversion
+            >>> Inversion = pf.dsp.RegularizedSpectrumInversion.from_frequency_range(
+            >>>     sweep, [50, 16e3], beta='max')
+            >>> inverted = Inversion.invert
+            ...
+            >>> # Obtain the scaled regularization function
+            >>> regularization = Inversion.regularization * Inversion.beta_value
+            ...
+            >>> pf.plot.use()
+            >>> fig, axes = plt.subplots(2,1)
+            >>> pf.plot.freq(sweep, ax=axes[0], label='sweep')
+            >>> pf.plot.freq(regularization, ax=axes[0], label='regularization')
+            >>> axes[0].axvline(50, color='k', linestyle='--',
+            ...     label='frequency range')
+            >>> axes[0].axvline(16e3, color='k', linestyle='--')
+            >>> axes[0].legend(loc='lower center')
+            ...
+            >>> pf.plot.freq(inverted, ax=axes[1], color='p',
+            ...     label='inverted with regulariztion')
+            >>> pf.plot.freq(1 / (sweep+1e-10), ax=axes[1], color='y',
+            ...     linestyle=':', label='inverted without regulariztion')
+            >>> axes[1].axvline(50, color='k', linestyle='--',
+            ...     label='frequency range')
+            >>> axes[1].axvline(16e3, color='k', linestyle='--')
+            >>> axes[1].set_ylim(-120, -20)
+            >>> axes[1].legend(loc='lower center')
+        """# noqa: E501
+        instance = cls.__new__(cls)
+
+        instance._frequency_range = np.asarray(frequency_range).flatten()
+
+        if instance._frequency_range.shape != (2, ):
+            raise ValueError(
+                "The frequency range needs to specify lower and upper limits.")
+
+        instance._regularization_type = "frequency range"
+        instance._regularization_within = regularization_within
+        instance.signal = signal
+        instance.beta = beta
+        instance.target = target
+
+        return instance
+
+    @classmethod
+    def from_magnitude_spectrum(cls, signal, regularization, beta=1,
+                                target=None):
+        r"""
+        Regularization from a given magnitude spectrum.
+
+        Regularization passed as :py:class:`~pyfar.Signal`. The length and
+        sampling rate of `regularization` must match the length and sampling
+        rate of the signal to be inverted.
+
+        Parameters
+        ----------
+        signal : Signal
+            Signal to be inverted.
+        regularization : Signal
+            Signal defining the regularization :math:`\epsilon(f)` in
+            ``regularization.freq``.
+        beta : float, string, optional
+            Beta parameter to control the scaling of the regularization as in
+            :eq:`regularized_inversion`. Can be a
+
+            ``numerical value``
+                Usually between ``0`` and ``1``, with ``0`` being no
+                regularization applied.
+            ``'energy'``
+                Normalize the regularization to match the signal's energy.
+
+                :math:`\beta = \frac{\frac{1}{N}\sum_{k=0}^{N-1}|S[k]|^2}{\frac{1}{N}\sum_{k=0}^{N-1}|\epsilon[k]|^2}`
+            ``'max'``
+                Normalize the regularization :math:`\epsilon(f)` to the maximum
+                magnitude of a given signal :math:`S(f)`.
+
+                :math:`\beta = \frac{\max(|S(f)|)}{\max(|\epsilon(f)|)}`
+            ``'mean'``
+                Normalize the regularization :math:`\epsilon(f)` to the mean
+                magnitude of a given signal :math:`S(f).`
+
+                :math:`\beta = \frac{\text{mean}(|S(f)|)}{\text{mean}(|\epsilon(f)|)}`
+
+            The default is ``1``.
+        target : Signal, optional
+            Target function for the regularization. The default ``None`` uses a
+            zero-phase spectrum with an amplitude of 1 as target, equal to an
+            impulse in the time domain.
+
+        Examples
+        --------
+        Invert a headphone transfer function (HpTF), regularize the inversion at
+        high frequencies, and use a band-pass as target function. Note that the
+        equalized HpTF, which is obtained from a convolution of the HpTF with its
+        inverse filter, approaches the target function in both, time and frequency.
+
+        .. plot::
+
+            >>> import pyfar as pf
+            >>> import numpy as np
+            >>>
+            >>> # Headphone transfer function to be inverted
+            >>> hptf = pf.signals.files.headphone_impulse_responses()[0, 0].flatten()
+            >>> # Regularize inversion at high frequencies
+            >>> regularization = pf.dsp.filter.low_shelf(
+            ...     pf.signals.impulse(hptf.n_samples), 4e3, -20, 2, 'II')
+            >>> # Use band-pass as target function
+            >>> target = pf.dsp.filter.butterworth(
+            ...     pf.signals.impulse(hptf.n_samples), 4, [100, 16e3], 'bandpass')
+            >>> target = pf.dsp.time_shift(target, 125, 'cyclic')
+            >>>
+            >>> # Regulated inversion
+            >>> Inversion = pf.dsp.RegularizedSpectrumInversion.from_magnitude_spectrum(
+            ...     hptf, regularization, beta=.1, target=target)
+            >>> inverted = Inversion.invert
+            >>>
+            >>> # Plot the result
+            >>> ax = pf.plot.time_freq(hptf, label='HpTF')
+            >>> pf.plot.time_freq(inverted, label='HpTF inverted')
+            >>> pf.plot.time_freq(hptf * inverted, label='Equalized HpTF')
+            >>> pf.plot.time_freq(target, linestyle='--', label='Target')
+            >>> pf.plot.freq(regularization, ax=ax[1], label='Regularization')
+            >>> ax[0].set_xlim(0 , .005)
+            >>> ax[1].set_ylim(-40 , 15)
+            >>> ax[1].legend(loc='lower center', ncols=3)
+        """# noqa: E501
+        if not isinstance(regularization, pf.Signal):
+            raise ValueError(
+                "Regularization must be a pyfar.Signal object.")
+        if signal.n_samples != regularization.n_samples:
+            raise ValueError(
+                "The number of samples in the signal and the regularization "
+                "function differs but must be equal.")
+        if signal.sampling_rate != regularization.sampling_rate:
+            raise ValueError(
+                "The sampling rate of the signal and the regularization "
+                "function differs but must be equal.")
+
+        instance = cls.__new__(cls)
+        instance._regularization_signal = regularization
+        instance._regularization_type = "magnitude spectrum"
+        instance.signal = signal
+        instance.beta = beta
+        instance.target = target
+
+        return instance
+
+    @property
+    def beta(self):
+        r"""
+        Scaling :math:`\beta` of regularization function :math:`\epsilon(f)`.
+
+        This can be ``'energy'``, ``'max'``, ``'mean'``, or a number.
+
+        To return the numeric :math:`\beta` value used in the inversion as in
+        equation :eq:`regularized_inversion`, use the property
+        :py:attr:`beta_value`.
+        """
+        return self._beta
+
+    @beta.setter
+    def beta(self, beta):
+        """Set beta parameter for scaling of regularization function."""
+        if not isinstance(beta, numbers.Number) \
+            and beta not in ['energy', 'mean', 'max']:
+            raise ValueError("Beta must be a scalar or 'energy', 'mean' or "
+                             "'max'.")
+        self._beta = beta
+        return self._beta
+
+    @property
+    def beta_value(self):
+        r"""Numeric :math:`\beta` value."""
+        if self.beta == 'mean':
+            self._beta_value = np.mean(np.abs(self.signal.freq)) / np.mean(
+                np.abs(self.regularization.freq))
+        elif self.beta == 'max':
+            self._beta_value = np.max(np.abs(self.signal.freq)) / np.max(
+                np.abs(self.regularization.freq))
+        elif self.beta == 'energy':
+            self._beta_value = \
+                pf.dsp.energy(self.signal) / pf.dsp.energy(self.regularization)
+        else:
+            self._beta_value = self.beta
+        return self._beta_value
+
+    @property
+    def regularization(self):
+        r"""
+        Regularization :math:`\epsilon(f)` without scaling by :math:`\beta`.
+        """
+        # Call private method to get regularization factors
+        if self._regularization_type == "frequency range":
+            self._regularization = \
+                self._get_regularization_from_frequency_range(self.signal)
+        elif self._regularization_type == "magnitude spectrum":
+            self._regularization =  \
+                self._get_regularization_from_regularization_signal(
+                    self._regularization_signal)
+        return self._regularization
+
+    @property
+    def signal(self):
+        r"""Signal :math:`S(f)` to be inverted."""
+        return self._signal
+
+    @signal.setter
+    def signal(self, signal):
+        """Set signal to be inverted."""
+        if not isinstance(signal, pf.Signal):
+            raise TypeError("Regularization of type "
+                            f"'{self._regularization_type}' requires an input"
+                            " signal of type pyfar.Signal.")
+        self._signal = signal
+
+    @property
+    def target(self):
+        """Target function :math:`D(f)`."""
+        return self._target
+
+    @target.setter
+    def target(self, target):
+        """Set target function."""
+        if target is not None and not isinstance(target, pf.Signal):
+            raise ValueError(
+                "Target function must be a pyfar.Signal object.")
+        if target is not None and target.n_samples != self.signal.n_samples:
+            raise ValueError(
+                "The number of samples in the signal and the target function"
+                " differs but must be equal.")
+        if (target is not None and
+            target.sampling_rate != self.signal.sampling_rate):
+            raise ValueError(
+                "The sampling rate of the signal and the target function"
+                " differs but must be equal.")
+        self._target = target
+
+    @property
+    def invert(self):
+        """
+        Invert signal using frequency-dependent regularized inversion.
+
+        Returns
+        -------
+        inverse : Signal
+            The inverted signal.
+        """# noqa: E501
+        data = self.signal.freq
+
+        # calculate inverse filter
+        inverse = self.signal.copy()
+        inverse.freq = \
+            np.conj(data) / (np.abs(data)**2 + self.beta_value *
+                             np.abs(self.regularization.freq)**2)
+
+        # Apply target function
+        if self.target is not None:
+            inverse.freq *= self.target.freq
+
+        return inverse
+
+
+    def _get_regularization_from_frequency_range(self, signal):
+        """Get regularization factors from frequency range."""
+        # arrange regularization factors for inside and outside frequency range
+        regu_inside = \
+            np.ones(signal.n_bins, dtype=np.double)*self._regularization_within
+        regu_outside = np.ones(signal.n_bins, dtype=np.double)
+
+        # index for crossfade at lower frequency limit
+        idx_xfade_lower = signal.find_nearest_frequency(
+            [self._frequency_range[0]/np.sqrt(2), self._frequency_range[0]])
+
+        regu_final = _cross_fade(regu_outside, regu_inside, idx_xfade_lower)
+
+        # index for crossfade at upper frequency limit
+        if self._frequency_range[1] < signal.sampling_rate/2:
+            idx_xfade_upper = signal.find_nearest_frequency([
+                self._frequency_range[1],
+                np.min([self._frequency_range[1]*np.sqrt(2),
+                        signal.sampling_rate/2])])
+
+            # crossfade regularization factors at frequency limits
+            regu_final = _cross_fade(regu_final, regu_outside, idx_xfade_upper)
+
+        # get number of samples according to the number of bins
+        n_samples = _n_samples_from_n_bins(signal.n_bins, signal.complex)
+
+        return pf.Signal(regu_final, signal.sampling_rate,
+                         n_samples, domain='freq')
+
+
+    def _get_regularization_from_regularization_signal(self, signal):
+        """Get regularization from signal."""
+        # get number of samples according to the number of bins
+        n_samples = _n_samples_from_n_bins(signal.n_bins, signal.complex)
+
+        return pf.Signal(np.abs(self._regularization_signal.freq),
+                         signal.sampling_rate, n_samples, domain='freq')


Test patch:
diff --git a/tests/test_dsp_regularization.py b/tests/test_dsp_regularization.py
new file mode 100644
index 00000000..ed96bf59
--- /dev/null
+++ b/tests/test_dsp_regularization.py
@@ -0,0 +1,152 @@
+import pyfar as pf
+import pytest
+import numpy as np
+import numpy.testing as npt
+import re
+
+
+def test_regularization_errors(impulse):
+    """Test errors."""
+    with pytest.raises(RuntimeError,
+                       match=re.escape("Regularization objects must be created"
+                                       " using one of the 'from_()' "
+                                       "classmethods.")):
+        assert pf.dsp.RegularizedSpectrumInversion()
+
+    with pytest.raises(ValueError,
+                       match="The frequency range needs to specify lower and"
+                        " upper limits."):
+        assert pf.dsp.RegularizedSpectrumInversion.from_frequency_range(
+            impulse, (200,))
+
+    with pytest.raises(ValueError,
+                       match="Beta must be a scalar or 'energy', 'mean' or "
+                             "'max'."):
+        assert pf.dsp.RegularizedSpectrumInversion.from_frequency_range(
+            impulse, (200, 15e3), beta='bla')
+
+    with pytest.raises(ValueError,
+                       match="Target function must be a pyfar.Signal object."):
+        assert pf.dsp.RegularizedSpectrumInversion.from_frequency_range(
+            impulse, (200, 20e3), target=1)
+
+    with pytest.raises(ValueError,
+                       match="The number of samples in the signal and the "
+                       "target function differs but must be equal."):
+        assert pf.dsp.RegularizedSpectrumInversion.from_frequency_range(
+            impulse, (200, 20e3), target=pf.signals.impulse(12))
+
+    with pytest.raises(ValueError,
+                       match="The sampling rate of the signal and the "
+                       "target function differs but must be equal."):
+        assert pf.dsp.RegularizedSpectrumInversion.from_frequency_range(
+            impulse, (200, 20e3),
+            target=pf.signals.impulse(impulse.n_samples, sampling_rate=48e3))
+
+    with pytest.raises(ValueError,
+                       match="Regularization must be a pyfar.Signal"
+                       " object."):
+        assert pf.dsp.RegularizedSpectrumInversion.from_magnitude_spectrum(
+            impulse, 1)
+
+    with pytest.raises(ValueError,
+                       match="The number of samples in the signal and the "
+                       "regularization function differs but must be equal."):
+        assert pf.dsp.RegularizedSpectrumInversion.from_magnitude_spectrum(
+            impulse, pf.signals.noise(12))
+
+    with pytest.raises(ValueError,
+                       match="The sampling rate of the signal and the "
+                       "regularization function differs but must be equal."):
+        assert pf.dsp.RegularizedSpectrumInversion.from_magnitude_spectrum(
+            impulse, pf.signals.noise(impulse.n_samples, sampling_rate=48e3))
+
+
+@pytest.mark.parametrize(("beta", "expected"), [(1, 0.5), (0.5, 2/3), (0, 1)])
+def test_regularization_frequency_range(impulse, beta, expected):
+    """Test Regularization from a frequency range using different beta
+    values.
+    """
+    Regu = pf.dsp.RegularizedSpectrumInversion.from_frequency_range(impulse,
+        (200, 10e3), beta=beta)
+    inv = Regu.invert
+
+    idx = inv.find_nearest_frequency([200, 10e3])
+
+    npt.assert_allclose(inv.freq[:, idx[0]:idx[1]], 1)
+    npt.assert_allclose(inv.freq[:, 0], expected)
+    npt.assert_allclose(inv.freq[:, -1], expected)
+
+
+@pytest.mark.parametrize(("norm", "expected"), [(1, 0.4),
+                                                ('max', 1/3)])
+def test_regularization_normalization(impulse, norm, expected):
+    """Test normalization parameter of Regularization."""
+    Regu = pf.dsp.RegularizedSpectrumInversion.from_frequency_range(impulse*2,
+        (200, 10e3), beta=norm)
+    inv = Regu.invert
+
+    idx = inv.find_nearest_frequency([200, 10e3])
+
+    npt.assert_allclose(inv.freq[:, idx[0]:idx[1]], 0.5)
+    npt.assert_allclose(inv.freq[:, 0], expected)
+    npt.assert_allclose(inv.freq[:, -1], expected)
+
+
+def test_regularization_within(impulse):
+    """Test regularization_within parameter."""
+    Regu = pf.dsp.RegularizedSpectrumInversion.from_frequency_range(impulse,
+        (200, 10e3), regularization_within=.5)
+    inv = Regu.invert
+
+    idx = inv.find_nearest_frequency([200, 10e3])
+
+    npt.assert_allclose(inv.freq[:, idx[0]:idx[1]], 0.8)
+
+
+def test_regularization_beta_mean_max(impulse):
+    """Test normalization factor of caclulated during inversion."""
+    impulse = impulse*2
+    Regu = pf.dsp.RegularizedSpectrumInversion.from_frequency_range(impulse,
+        (200, 10e3))
+
+    Regu.beta = 'energy'
+    assert Regu.beta_value == \
+        pf.dsp.energy(impulse) / pf.dsp.energy(Regu.regularization)
+    assert Regu.beta == 'energy'
+
+    Regu.beta = 'max'
+    assert Regu.beta_value == \
+        np.max(np.abs(impulse.freq)) / np.max(np.abs(Regu.regularization.freq))
+    assert Regu.beta == 'max'
+
+    Regu.beta = 'mean'
+    assert Regu.beta_value == \
+        np.mean(np.abs(impulse.freq)) / \
+            np.mean(np.abs(Regu.regularization.freq))
+    assert Regu.beta == 'mean'
+
+
+def test_regularization_target(impulse):
+    """Test Regularization from a frequency range using a bandpass filter
+    target function.
+    """
+    bp = pf.dsp.filter.butterworth(None, 4, (20, 15.5e3), 'bandpass',
+                                   impulse.sampling_rate)
+    target = bp.process(impulse)
+
+    ReguTarget = pf.dsp.RegularizedSpectrumInversion.from_frequency_range(
+        impulse, (20, 15.5e3), target=target)
+
+    inv = ReguTarget.invert
+    idx = inv.find_nearest_frequency([20, 15.5e3])
+
+    npt.assert_allclose(np.abs(inv.freq[0, idx]), np.abs(target.freq[0, idx]))
+
+
+def test_regularization_from_magnitude_spectrum(noise):
+    """Test Regularization from an arbitrary signal."""
+    Regu = pf.dsp.RegularizedSpectrumInversion.from_magnitude_spectrum(
+        pf.signals.impulse(noise.n_samples), noise)
+
+    npt.assert_allclose(Regu.regularization.freq, np.abs(noise.freq))