Repository: pyfar/pyfar
Pull request #598 (2024-04-29 06:42:20)

Test find impulse response start
### Which issue(s) are closed by this pull request?

Closes #

### Changes proposed in this pull request:

-
-
-



Patch:
diff --git a/pyfar/dsp/dsp.py b/pyfar/dsp/dsp.py
index 07bcfdca..0771cf67 100644
--- a/pyfar/dsp/dsp.py
+++ b/pyfar/dsp/dsp.py
@@ -1342,7 +1342,8 @@ def find_impulse_response_start(
     ``threshold`` level relative to the maximum level of the impulse response.
     For room impulse responses, ISO 3382 [#]_ specifies a threshold of 20 dB.
     This function is primary intended to be used when processing room impulse
-    responses.
+    responses. For complex-valued time signals the onset is computed separately
+    for the real and imaginary part.
     Alternatively see :py:func:`pyfar.dsp.find_impulse_response_delay`.
 
 
@@ -1422,43 +1423,50 @@ def find_impulse_response_start(
         >>> ax.legend()
 
     """
-    ir_squared = np.abs(impulse_response.time)**2
-
-    mask_start = int(0.9*impulse_response.n_samples)
-
-    mask = np.arange(mask_start, ir_squared.shape[-1])
-    noise = np.mean(np.take(ir_squared, mask, axis=-1), axis=-1)
-
-    max_sample = np.argmax(ir_squared, axis=-1)
-    max_value = np.max(ir_squared, axis=-1)
-
-    if np.any(max_value < 10**(threshold/10) * noise) or \
-            np.any(max_sample > mask_start):
-        warnings.warn(
-            "The SNR seems lower than the specified threshold value. Check "
-            "if this is a valid impulse response with sufficient SNR.")
-
-    start_sample = max_sample.copy()
-
-    for ch in np.ndindex(impulse_response.cshape):
-        # Only look for the start sample if the maximum index is bigger than 0
-        if start_sample[ch] > 0:
-            # Check samples before maximum
-            ir_before_max = np.squeeze(
-                ir_squared[ch][:max_sample[ch]+1] / max_value[ch])
-            # First sample above or at the threshold level
-            idx_first_above_thresh = np.where(
-                ir_before_max >= 10**(-threshold/10))[0]
-            if idx_first_above_thresh.size > 0:
-                # The start sample is the last sample below the threshold
-                start_sample[ch] = np.min(idx_first_above_thresh) - 1
-            else:
-                start_sample[ch] = 0
-                warnings.warn(
-                    f'No values below threshold found found for channel {ch}',
-                    'defaulting to 0')
+    modes = ['real', 'complex'] if impulse_response.complex else ['real']
+    ir_start = np.zeros(((len(modes), ) + impulse_response.cshape))
+    for idx, mode in enumerate(modes):
+        ir = impulse_response.time
+        ir_squared = np.real(ir)**2 if mode == 'real' else np.imag(ir)**2
+
+        mask_start = int(0.9*impulse_response.n_samples)
+
+        mask = np.arange(mask_start, ir_squared.shape[-1])
+        noise = np.mean(np.take(ir_squared, mask, axis=-1), axis=-1)
+
+        max_sample = np.argmax(ir_squared, axis=-1)
+        max_value = np.max(ir_squared, axis=-1)
+
+        if np.any(max_value < 10**(threshold/10) * noise) or \
+                np.any(max_sample > mask_start):
+            warnings.warn(
+                "The SNR seems lower than the specified threshold value. "
+                "Check if this is a valid impulse response with sufficient "
+                "SNR.")
+
+        start_sample = max_sample.copy()
+
+        for ch in np.ndindex(impulse_response.cshape):
+            # Only look for the start sample if the maximum index is bigger
+            # than 0
+            if start_sample[ch] > 0:
+                # Check samples before maximum
+                ir_before_max = np.squeeze(
+                    ir_squared[ch][:max_sample[ch]+1] / max_value[ch])
+                # First sample above or at the threshold level
+                idx_first_above_thresh = np.where(
+                    ir_before_max >= 10**(-threshold/10))[0]
+                if idx_first_above_thresh.size > 0:
+                    # The start sample is the last sample below the threshold
+                    start_sample[ch] = np.min(idx_first_above_thresh) - 1
+                else:
+                    start_sample[ch] = 0
+                    warnings.warn(
+                        'No values below threshold found found for channel '
+                        f'{ch} defaulting to 0')
 
-    return np.squeeze(start_sample)
+        ir_start[idx, :] = np.squeeze(start_sample)
+    return np.squeeze(ir_start)
 
 
 @rename_arg({"freq_range": "frequency_range"},


Test patch:
diff --git a/tests/test_dsp.py b/tests/test_dsp.py
index 8a75b8d3..7deb2135 100644
--- a/tests/test_dsp.py
+++ b/tests/test_dsp.py
@@ -875,6 +875,57 @@ def test_impulse_response_start_multidim():
     npt.assert_allclose(start_sample_est, start_samples - 1)
 
 
+def test_impulse_response_start_multidim_complex():
+    n_samples = 2**10
+    n_channels = 3
+    ir = np.zeros((n_channels, n_samples), dtype=complex)
+
+    snr = 60
+
+    noise = pf.Signal(
+        np.random.randn(n_channels, n_samples) * 10**(-snr/20), 44100)
+
+    start_sample_real = [24, 5, 43]
+    ir[[0, 1, 2], start_sample_real] = 1
+
+    start_sample_imag = [14, 3, 20]
+    ir[[0, 1, 2], start_sample_imag] = 1j
+
+    ir = pf.Signal(ir, 44100, is_complex=True)
+
+    ir_awgn = ir + noise
+    start_sample_est = dsp.find_impulse_response_start(ir_awgn)
+
+    npt.assert_allclose(start_sample_est[0, :],
+                        np.array(start_sample_real) - 1)
+    npt.assert_allclose(start_sample_est[1, :],
+                        np.array(start_sample_imag) - 1)
+
+    ir = np.zeros((2, n_channels, n_samples), dtype=complex)
+    noise = pf.Signal(
+        np.random.randn(2, n_channels, n_samples) * 10**(-snr/20), 44100)
+
+    start_sample_real_1 = [24, 5, 43]
+    start_sample_imag_1 = [14, 3, 20]
+    ir[0, [0, 1, 2], start_sample_real_1] = 1
+    ir[0, [0, 1, 2], start_sample_imag_1] = 1j
+    start_sample_real_2 = [14, 12, 16]
+    start_sample_imag_2 = [35, 7, 3]
+    ir[1, [0, 1, 2], start_sample_real_2] = 1
+    ir[1, [0, 1, 2], start_sample_imag_2] = 1j
+
+    ir = pf.Signal(ir, 44100, is_complex=True)
+
+    start_samples_real = np.vstack((start_sample_real_1, start_sample_real_2))
+    start_samples_imag = np.vstack((start_sample_imag_1, start_sample_imag_2))
+
+    ir_awgn = ir + noise
+    start_sample_est = dsp.find_impulse_response_start(ir_awgn)
+
+    npt.assert_allclose(start_sample_est[0, :, :], start_samples_real - 1)
+    npt.assert_allclose(start_sample_est[1, :, :], start_samples_imag - 1)
+
+
 def test_convolve_default():
     x = pf.Signal([1, 0.5, 0.25, 0], 44100)
     y = pf.Signal([1, -1, 0], 44100)