Repository: pyfar/pyfar
Pull request #99 (2021-02-12 19:59:18)

Regularized spectrum inversion
### Changes proposed in this pull request:

- Add a function for the regularized inversion of signal spectra 



Patch:
diff --git a/pyfar/dsp/__init__.py b/pyfar/dsp/__init__.py
index ffb5a3c5..6f6999c3 100644
--- a/pyfar/dsp/__init__.py
+++ b/pyfar/dsp/__init__.py
@@ -1,9 +1,11 @@
 from .classes import (Filter, FilterFIR, FilterIIR, FilterSOS)
 from .dsp import (
-    phase, group_delay, wrap_to_2pi, nextpow2, spectrogram)
+    phase, group_delay, wrap_to_2pi, nextpow2, spectrogram,
+    regularized_spectrum_inversion)
 
 
 __all__ = [
-    Filter, FilterFIR, FilterIIR, FilterSOS,
-    phase, group_delay, wrap_to_2pi, nextpow2, spectrogram
+    'Filter', 'FilterFIR', 'FilterIIR', 'FilterSOS',
+    'phase', 'group_delay', 'wrap_to_2pi', 'nextpow2', 'spectrogram',
+    'regularized_spectrum_inversion'
 ]
diff --git a/pyfar/dsp/dsp.py b/pyfar/dsp/dsp.py
index 78417233..420cc15f 100644
--- a/pyfar/dsp/dsp.py
+++ b/pyfar/dsp/dsp.py
@@ -1,5 +1,6 @@
 import numpy as np
 from scipy import signal as sgn
+import pyfar
 from pyfar import Signal, FrequencyData
 import pyfar.fft as fft
 
@@ -225,3 +226,133 @@ def spectrogram(signal, dB=True, log_prefix=20, log_reference=1,
     times -= times[0]
 
     return frequencies, times, spectrogram
+
+
+def regularized_spectrum_inversion(
+        signal, freq_range,
+        regu_outside=1., regu_inside=10**(-200/20), regu_final=None):
+    r"""Invert the spectrum of a signal applying frequency dependent
+    regularization. Regularization can either be specified within a given
+    frequency range using two different regularization factors, or for each
+    frequency individually using the parameter `regu_final`. In the first case
+    the regularization factors for the frequency regions are cross-faded using
+    a raised cosine window function with a width of `math:f*\sqrt(2)` above and
+    below the given frequency range. Note that the resulting regularization
+    function is adjusted to the quadratic maximum of the given signal.
+    In case the `regu_final` parameter is used, all remaining options are
+    ignored and an array matching the number of frequency bins of the signal
+    needs to be given. In this case, no normalization of the regularization
+    function is applied. Finally, the inverse spectrum is calculated as
+    [1]_, [2]_,
+
+    .. math::
+
+        S^{-1}(f) = \frac{S^*(f)}{S^*(f)S(f) + \epsilon(f)}
+
+
+    Parameters
+    ----------
+    signal : pyfar.Signal
+        The signals which spectra are to be inverted.
+    freq_range : tuple, array_like, double
+        The upper and lower frequency limits outside of which the
+        regularization factor is to be applied.
+    regu_outside : float, optional
+        The normalized regularization factor outside the frequency range.
+        The default is 1.
+    regu_inside : float, optional
+        The normalized regularization factor inside the frequency range.
+        The default is 10**(-200/20).
+    regu_final : float, array_like, optional
+        The final regularization factor for each frequency, by default None.
+        If this parameter is set, the remaining regularization factors are
+        ignored.
+
+    Returns
+    -------
+    pyfar.Signal
+        The resulting signal after inversion.
+
+    References
+    ----------
+    .. [1]  O. Kirkeby and P. A. Nelson, “Digital Filter Designfor Inversion
+            Problems in Sound Reproduction,” J. Audio Eng. Soc., vol. 47,
+            no. 7, p. 13, 1999.
+
+    .. [2]  P. C. Hansen, Rank-deficient and discrete ill-posed problems:
+            numerical aspects of linear inversion. Philadelphia: SIAM, 1998.
+
+    """
+    if not isinstance(signal, pyfar.Signal):
+        raise ValueError("The input signal needs to be of type pyfar.Signal.")
+
+    data = signal.freq
+    freq_range = np.asarray(freq_range)
+
+    if freq_range.size < 2:
+        raise ValueError(
+            "The frequency range needs to specify lower and upper limits.")
+
+    if regu_final is None:
+        regu_inside = np.ones(signal.n_bins, dtype=np.double) * regu_inside
+        regu_outside = np.ones(signal.n_bins, dtype=np.double) * regu_outside
+
+        idx_xfade_lower = signal.find_nearest_frequency(
+                [freq_range[0]/np.sqrt(2), freq_range[0]])
+
+        regu_final = _cross_fade(regu_outside, regu_inside, idx_xfade_lower)
+
+        if freq_range[1] < signal.sampling_rate/2:
+            idx_xfade_upper = signal.find_nearest_frequency([
+                freq_range[1],
+                np.min([freq_range[1]*np.sqrt(2), signal.sampling_rate/2])])
+
+            regu_final = _cross_fade(regu_final, regu_outside, idx_xfade_upper)
+
+        regu_final *= np.max(np.abs(data)**2)
+
+    inverse = signal.copy()
+    inverse.freq = np.conj(data) / (np.conj(data)*data + regu_final)
+
+    return inverse
+
+
+def _cross_fade(first, second, indices):
+    """Cross-fade two numpy arrays by multiplication with a raised cosine
+    window inside the range specified by the indices. Outside the range, the
+    result will be the respective first or second array, without distortions.
+
+    Parameters
+    ----------
+    first : array, double
+        The first array.
+    second : array, double
+        The second array.
+    indices : array-like, tuple, int
+        The lower and upper cross-fade indices.
+
+    Returns
+    -------
+    result : array, double
+        The resulting array after cross-fading.
+    """
+    indices = np.asarray(indices)
+    if np.shape(first)[-1] != np.shape(second)[-1]:
+        raise ValueError("Both arrays need to be of same length.")
+    len_arrays = np.shape(first)[-1]
+    if np.any(indices > np.shape(first)[-1]):
+        raise IndexError("Index is out of range.")
+
+    len_xfade = np.squeeze(np.abs(np.diff(indices)))
+    window = sgn.windows.windows.hann(len_xfade*2 + 1, sym=True)
+    window_rising = window[:len_xfade]
+    window_falling = window[len_xfade+1:]
+
+    window_first = np.concatenate(
+        (np.ones(indices[0]), window_falling, np.zeros(len_arrays-indices[1])))
+    window_second = np.concatenate(
+        (np.zeros(indices[0]), window_rising, np.ones(len_arrays-indices[1])))
+
+    result = first * window_first + second * window_second
+
+    return result


Test patch:
diff --git a/pyfar/testing/stub_utils.py b/pyfar/testing/stub_utils.py
index e2f10f57..67ecb8b5 100644
--- a/pyfar/testing/stub_utils.py
+++ b/pyfar/testing/stub_utils.py
@@ -47,6 +47,18 @@ def signal_stub(time, freq, sampling_rate, fft_norm):
             signal.fft_norm)
         return item
 
+    def find_nearest_time(times):
+        samples = np.zeros(len(times), dtype=int)
+        for idx, time in enumerate(times):
+            samples[idx] = np.argmin(np.abs(signal.times-time))
+        return samples
+
+    def find_nearest_frequency(freqs):
+        bin = np.zeros(len(freqs), dtype=int)
+        for idx, freq in enumerate(freqs):
+            bin[idx] = np.argmin(np.abs(signal.frequencies-freq))
+        return bin
+
     signal = mock.MagicMock(
         spec_set=Signal(time, sampling_rate, domain='time'))
     signal.time = np.atleast_2d(time)
@@ -61,6 +73,8 @@ def signal_stub(time, freq, sampling_rate, fft_norm):
     signal.frequencies = np.atleast_1d(
         np.fft.rfftfreq(signal.n_samples, 1 / sampling_rate))
     signal.__getitem__.side_effect = getitem
+    signal.find_nearest_time = find_nearest_time
+    signal.find_nearest_frequency = find_nearest_frequency
 
     return signal
 
diff --git a/tests/test_dsp.py b/tests/test_dsp.py
index 69a7c8ce..b9336c02 100644
--- a/tests/test_dsp.py
+++ b/tests/test_dsp.py
@@ -100,3 +100,44 @@ def test_group_delay_custom_frequencies(impulse_group_delay):
     grp = dsp.group_delay(signal, frequency, method='scipy')
     assert grp.shape == (2,)
     npt.assert_allclose(grp, impulse_group_delay[1][0, frequency_idx])
+
+
+def test_xfade(impulse):
+    first = np.ones(5001)
+    idx_1 = 500
+    second = np.ones(5001)*2
+    idx_2 = 1000
+
+    res = dsp.dsp._cross_fade(first, second, [idx_1, idx_2])
+    np.testing.assert_array_almost_equal(first[:idx_1], res[:idx_1])
+    np.testing.assert_array_almost_equal(second[idx_2:], res[idx_2:])
+
+    idx_1 = 501
+    idx_2 = 1000
+    res = dsp.dsp._cross_fade(first, second, [idx_1, idx_2])
+    np.testing.assert_array_almost_equal(first[:idx_1], res[:idx_1])
+    np.testing.assert_array_almost_equal(second[idx_2:], res[idx_2:])
+
+
+def test_regu_inversion(impulse):
+
+    with pytest.raises(
+            ValueError, match='needs to be of type pyfar.Signal'):
+        dsp.regularized_spectrum_inversion('error', (1, 2))
+
+    with pytest.raises(
+            ValueError, match='lower and upper limits'):
+        dsp.regularized_spectrum_inversion(impulse, (2))
+
+    impulse.freq = impulse.freq*2
+    impulse.time = impulse.time*2
+
+    res = dsp.regularized_spectrum_inversion(impulse, [200, 10e3])
+
+    ind = impulse.find_nearest_frequency([200, 10e3])
+    npt.assert_allclose(
+        res.freq[:, ind[0]:ind[1]],
+        np.ones((1, ind[1]-ind[0]), dtype=complex)*0.5)
+
+    npt.assert_allclose(res.freq[:, 0], [0.25])
+    npt.assert_allclose(res.freq[:, -1], [0.25])