Repository: pyfar/pyrato

Issue:
[BUG] EDC-Functions can't handle Signals with cdim > 1
### Is there an existing issue for this?

- [x] I have searched the existing issues

### Problem description

The energy decay curve functions can't seem to handle signals with `cdim > 1`. This becomes relevant when for example using the output of `fractional_octave_bands` where the returned signal is of cshape `(n_bands, *cshape)`. The example below works when squeezing the fractional octave signal to cshape `(n_bands,)`

### Code sample

Code run:

```
rir = pf.signals.files.room_impulse_response(sampling_rate=44100)
frac = pf.dsp.filter.fractional_octave_bands(rir, 1, frequency_range=(100, 16e3))
edc = ra.edc.energy_decay_curve_lundeby(frac)
```

Traceback:

```
---------------------------------------------------------------------------
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[105], line 4
      2 rir = pf.signals.files.room_impulse_response(sampling_rate=44100)
      3 frac = pf.dsp.filter.fractional_octave_bands(rir, 1, frequency_range=(100, 16e3))
----> 4 edc = ra.edc.energy_decay_curve_lundeby(frac)
      5 pf.plot.time_freq(frac)
      6 plt.show()

File ~/Documents/git-repositories/pyrato/pyrato/edc.py:414, in energy_decay_curve_lundeby(data, freq, noise_level, is_energy, time_shift, channel_independent, normalize, plot)
    410 n_samples = energy_data.n_samples
    411 sampling_rate = data.sampling_rate
    413 intersection_time, late_reverberation_time, noise_estimation = \
--> 414     intersection_time_lundeby(
    415         energy_data,
    416         freq=freq,
    417         initial_noise_power=noise_level,
    418         is_energy=True,
    419         time_shift=False,
    420         channel_independent=False,
    421         plot=False)
    422 time_vector = data.times
    424 energy_decay_curve = np.zeros([*data.cshape, n_samples])

File ~/Documents/git-repositories/pyrato/pyrato/edc.py:872, in intersection_time_lundeby(data, freq, initial_noise_power, is_energy, time_shift, channel_independent, plot)
    869     freq_dependent_window_time = (800/freq+10) / 1000
    871 # (1) SMOOTH
--> 872 time_window_data, time_vector_window, time_vector = dsp._smooth_rir(
    873     energy_data, sampling_rate, freq_dependent_window_time)
    875 # (2) ESTIMATE NOISE
    876 if initial_noise_power == 'auto':

File ~/Documents/git-repositories/pyrato/pyrato/dsp.py:444, in _smooth_rir(data, sampling_rate, smooth_block_length)
    442 n_blocks_min = int(np.min(n_blocks))
    443 n_samples_actual = int(n_blocks_min*n_samples_per_block)
--> 444 reshaped_array = np.reshape(
    445     data[..., :n_samples_actual],
    446     (-1, n_blocks_min, n_samples_per_block))
    447 time_window_data = np.mean(reshaped_array, axis=-1)
    449 # Use average time instances corresponding to the average energy level
    450 # instead of time for the first sample of the block

File /opt/anaconda3/envs/pyrato/lib/python3.13/site-packages/numpy/_core/fromnumeric.py:324, in reshape(a, shape, order, newshape, copy)
    322 if copy is not None:
    323     return _wrapfunc(a, 'reshape', shape, order=order, copy=copy)
--> 324 return _wrapfunc(a, 'reshape', shape, order=order)

File /opt/anaconda3/envs/pyrato/lib/python3.13/site-packages/numpy/_core/fromnumeric.py:57, in _wrapfunc(obj, method, *args, **kwds)
     54     return _wrapit(obj, method, *args, **kwds)
     56 try:
---> 57     return bound(*args, **kwds)
     58 except TypeError:
     59     # A TypeError occurs if the object does have such a method in its
     60     # class, but its signature is not identical to that of NumPy's. This
   (...)     64     # Call _wrapit from within the except clause to ensure a potential
     65     # exception has a traceback chain.
     66     return _wrapit(obj, method, *args, **kwds)

ValueError: cannot reshape array of size 0 into shape (0,1323)

```

</details>


### Anything else?

_No response_

Discussion:

Fix (patch):
diff --git a/pyrato/edc.py b/pyrato/edc.py
index 9bc300d..6fd5310 100644
--- a/pyrato/edc.py
+++ b/pyrato/edc.py
@@ -263,6 +263,10 @@ def energy_decay_curve_truncation(
         >>> ax.legend()
 
     """
+    # flatten to allow signals with cdim > 1
+    shape = data.time.shape
+    data = data.flatten()
+
     energy_data = dsp.preprocess_rir(
         data,
         is_energy=is_energy,
@@ -311,7 +315,7 @@ def energy_decay_curve_truncation(
                     / np.nanmax(energy_decay_curve, axis=-1, keepdims=True)
 
     edc = pf.TimeData(
-        energy_decay_curve, data.times, comment=data.comment)
+        energy_decay_curve.reshape(shape), data.times, comment=data.comment)
 
     if plot:
         ax = pf.plot.time(data, dB=True, label='RIR')
@@ -403,6 +407,9 @@ def energy_decay_curve_lundeby(
         >>> ax.legend()
 
     """
+    # flatten to allow signals with cdim > 1
+    shape = data.time.shape
+    data = data.flatten()
 
     energy_data = dsp.preprocess_rir(
         data,
@@ -458,7 +465,7 @@ def energy_decay_curve_lundeby(
                     / np.nanmax(energy_decay_curve, axis=-1, keepdims=True)
 
     edc = pf.TimeData(
-        energy_decay_curve, data.times, comment=data.comment)
+        energy_decay_curve.reshape(shape), data.times, comment=data.comment)
 
     if plot:
         ax = pf.plot.time(data, dB=True, label='RIR')
@@ -551,6 +558,10 @@ def energy_decay_curve_chu(
         >>> ax.legend()
 
     """
+    # flatten to allow signals with cdim > 1
+    shape = data.cshape
+    data = data.flatten()
+
     energy_data = dsp.preprocess_rir(
         data,
         is_energy=is_energy,
@@ -588,6 +599,8 @@ def energy_decay_curve_chu(
         trunc_levels = 10*np.log10((psnr)) - threshold
         edc = truncate_energy_decay_curve(edc, trunc_levels)
 
+    edc = edc.reshape(shape)
+
     if plot:
         plt.figure(figsize=(15, 3))
         pf.plot.use('light')
@@ -693,6 +706,9 @@ def energy_decay_curve_chu_lundeby(
         >>> ax.legend()
 
     """
+    # flatten to allow signals with cdim > 1
+    shape = data.time.shape
+    data = data.flatten()
 
     energy_data = dsp.preprocess_rir(
         data,
@@ -754,7 +770,7 @@ def energy_decay_curve_chu_lundeby(
                     / np.nanmax(energy_decay_curve, axis=-1, keepdims=True)
 
     edc = pf.TimeData(
-        energy_decay_curve, data.times, comment=data.comment)
+        energy_decay_curve.reshape(shape), data.times, comment=data.comment)
 
     if plot:
         ax = pf.plot.time(data, dB=True, label='RIR')


Test changes:
diff --git a/tests/test_edc.py b/tests/test_edc.py
index 17f2cf1..7770ae9 100644
--- a/tests/test_edc.py
+++ b/tests/test_edc.py
@@ -96,3 +96,30 @@ def test_edc_sabine():
         6.37107964e-04, 5.46555336e-04,
         ])
     npt.assert_almost_equal(edc, truth)
+
+
+@pytest.mark.parametrize(
+        "edc_function",
+        [ra.edc.energy_decay_curve_chu,
+         ra.edc.energy_decay_curve_lundeby,
+         ra.edc.energy_decay_curve_chu_lundeby,
+         ra.energy_decay_curve_lundeby],
+)
+def test_multidim_edc(edc_function):
+    """
+    Test if edcs from multichannel signal are equal to corresponding single
+    channel edcs.
+    """
+    rir = pf.signals.files.room_impulse_response()
+    rir_oct = pf.dsp.filter.fractional_octave_bands(rir, 1)
+    shape = rir_oct.time.shape
+    edc = edc_function(rir_oct, channel_independent=True)
+
+    assert shape == edc.time.shape
+
+    edc = edc.flatten()
+    rir_oct = rir_oct.flatten()
+
+    for i in range(edc.cshape[0]):
+        baseline = edc_function(rir_oct[i])
+        npt.assert_array_equal(edc[i].time, baseline.time)