Group delay calculation
## Description

pyfar.dps.group_delay() is buggy

## What I Did
The example below should yield a group delay of 1000 samples, but yields a group delay of 0 samples. The scipy implementation behaves correct and should be used instead

```
import numpy as np
import pyfar
import scipy as sp
import matplotlib.pyplot as plt

fs = 40e3
x = np.zeros(2000)
x[1000] = 1

s = pyfar.Signal(x, fs)
grp_pyfar = pyfar.dsp.group_delay(s)
grp_scipy = sp.signal.group_delay((x, 1), s.frequencies, fs=fs)[1]

plt.plot(s.frequencies, grp_pyfar, label='pyfar')
plt.plot(s.frequencies, grp_scipy, label='scipy')
ax = plt.gca()
ax.legend()
```
diff --git a/pyfar/dsp/dsp.py b/pyfar/dsp/dsp.py
index 0e5f4924..ad66ee14 100644
--- a/pyfar/dsp/dsp.py
+++ b/pyfar/dsp/dsp.py
@@ -40,28 +40,48 @@ def phase(signal, deg=False, unwrap=False):
     return phase
 
 
-def group_delay(signal):
-    """Returns the group delay for a given signal object.
+def group_delay(signal, frequencies=None):
+    """Returns the group delay of a signal in samples.
 
     Parameters
     ----------
     signal : Signal object
         An audio signal object from the pyfar signal class
+    frequencies : number array like
+        Frequency or frequencies in Hz at which the group delay is calculated.
+        The default is None, in which case signal.frequencies is used.
 
     Returns
     -------
-    group_delay : np.array()
-        Group delay.
+    group_delay : numpy array
+        Frequency dependent group delay in samples. The array is flattened if
+        a single channel signal was passed to the function.
     """
 
+    # check input and default values
     if not isinstance(signal, Signal):
         raise TypeError('Input data has to be of type: Signal.')
 
-    phase_vec = phase(signal, deg=False, unwrap=True)
-    bin_dist = signal.sampling_rate / signal.n_samples
-    group_delay = (- np.diff(phase_vec, 1, -1, prepend=0) /
-                   (bin_dist * 2 * np.pi))
-    return np.squeeze(group_delay)
+    frequencies = signal.frequencies if frequencies is None \\
+        else np.asarray(frequencies)
+
+    # get time signal and reshape for easy looping
+    time = signal.time
+    time = time.reshape((-1, signal.n_samples))
+    # initialize group delay
+    group_delay = np.zeros((np.prod(signal.cshape), frequencies.size))
+    # calculate the group delay
+    for cc in range(time.shape[0]):
+        group_delay[cc] = sgn.group_delay(
+            (time[cc], 1), frequencies, fs=signal.sampling_rate)[1]
+    # reshape to match signal
+    group_delay = group_delay.reshape(signal.cshape + (-1, ))
+
+    # flatten in numpy fashion if a single channel is returned
+    if signal.cshape == (1, ):
+        group_delay = np.squeeze(group_delay)
+
+    return group_delay
 
 
 def spectrogram(signal,
diff --git a/tests/test_dsp.py b/tests/test_dsp.py
index 6a4ec7a4..775c5177 100644
--- a/tests/test_dsp.py
+++ b/tests/test_dsp.py
@@ -1,47 +1,141 @@
+from pyfar.dsp.dsp import group_delay
 import pytest
 import numpy as np
 import numpy.testing as npt
-import matplotlib.pyplot as plt
 from unittest import mock
+import copy
 from pyfar import Signal
 from pyfar import dsp
 
+
 def test_phase_rad(sine_plus_impulse_mock):
     """Test the function returning the phase of a signal in radians."""
     phase = dsp.phase(sine_plus_impulse_mock, deg=False, unwrap=False)
     truth = np.angle(sine_plus_impulse_mock.freq)
     npt.assert_allclose(phase, truth, rtol=1e-7)
 
+
 def test_phase_deg(sine_plus_impulse_mock):
     """Test the function returning the phase of a signal in degrees."""
     phase = dsp.phase(sine_plus_impulse_mock, deg=True, unwrap=False)
     truth = np.degrees(np.angle(sine_plus_impulse_mock.freq))
     npt.assert_allclose(phase, truth, rtol=1e-7)
 
+
 def test_phase_unwrap(sine_plus_impulse_mock):
     """Test the function returning the unwrapped phase of a signal."""
     phase = dsp.phase(sine_plus_impulse_mock, deg=False, unwrap=True)
     truth = np.unwrap(np.angle(sine_plus_impulse_mock.freq))
     npt.assert_allclose(phase, truth, rtol=1e-7)
 
+
 def test_phase_deg_unwrap(sine_plus_impulse_mock):
     """Test the function returning the unwrapped phase of a signal in deg."""
     phase = dsp.phase(sine_plus_impulse_mock, deg=True, unwrap=True)
     truth = np.degrees(np.unwrap(np.angle(sine_plus_impulse_mock.freq)))
     npt.assert_allclose(phase, truth, rtol=1e-7)
 
-def test_group_delay(sine_plus_impulse_mock):
+
+def test_group_delay_single_channel(impulse_mock):
+    """Test the function returning the group delay of a signal."""
+    # test single channel signal
+    signal = impulse_mock[0]
+    grp = dsp.group_delay(signal)
+    assert grp.shape == (signal.n_bins, )
+    npt.assert_allclose(grp, impulse_mock[3][0] * np.ones(signal.n_bins))
+
+
+def test_group_delay_two_channel(impulse_mock):
+    """Test the function returning the group delay of a signal."""
+    # test two channel signal
+    signal = impulse_mock[1]
+    grp = dsp.group_delay(signal)
+    assert grp.shape == (2, signal.n_bins)
+    npt.assert_allclose(grp[0], impulse_mock[3][0] * np.ones(signal.n_bins))
+    npt.assert_allclose(grp[1], impulse_mock[3][1] * np.ones(signal.n_bins))
+
+
+def test_group_delay_two_by_two_channel(impulse_mock):
+    """Test the function returning the group delay of a signal."""
+    # test two by two channel signal
+    signal = impulse_mock[2]
+    grp = dsp.group_delay(signal)
+    assert grp.shape == (2, 2, signal.n_bins)
+    npt.assert_allclose(grp[0, 0], impulse_mock[3][0] * np.ones(signal.n_bins))
+    npt.assert_allclose(grp[0, 1], impulse_mock[3][1] * np.ones(signal.n_bins))
+    npt.assert_allclose(grp[1, 0], impulse_mock[3][2] * np.ones(signal.n_bins))
+    npt.assert_allclose(grp[1, 1], impulse_mock[3][3] * np.ones(signal.n_bins))
+
+
+def test_group_delay_custom_frequencies(impulse_mock):
     """Test the function returning the group delay of a signal."""
-    group_delay = dsp.group_delay(sine_plus_impulse_mock)
-    phase = np.unwrap(np.angle(sine_plus_impulse_mock.freq))
-    bin_dist = (sine_plus_impulse_mock.sampling_rate /
-                    sine_plus_impulse_mock.n_samples)
-    truth = np.squeeze((- np.diff(phase,1,-1, prepend=0) /
-                   (bin_dist * 2 * np.pi)))
-    npt.assert_allclose(group_delay, truth, rtol=1e-7)
+    # test single frequency
+    signal = impulse_mock[0]
+    grp = dsp.group_delay(signal, 1e3)
+    assert grp.shape == ()
+    npt.assert_allclose(grp, 1000)
+
+    # test multiple frequencies
+    signal = impulse_mock[0]
+    grp = dsp.group_delay(signal, [1e3, 2e3])
+    assert grp.shape == (2, )
+    npt.assert_allclose(grp, np.array([1e3, 1e3]))
+
+# def test_wrap_to_2pi():
+# def test_nextpow2():
+
+
+@pytest.fixture
+def impulse_mock():
+    """ Generate impulse signals, in order to test independently of the Signal
+    object.
+
+    Returns
+    -------
+    signal_1 : Signal
+        single channel dirac signal, delay 1000 samples
+    signal_2 : Signal
+        two channel dirac signal, delays 100 and 750 samples
+    signal_3 : Signal
+        two by two channel dirac signal, delays [[1000, 750], [500, 250]]
+
+    """
+    n_samples = 2000
+    n_bins = int(n_samples / 2 + 1)
+    sampling_rate = 4000
+    times = np.arange(0, n_samples) / sampling_rate
+    frequencies = np.arange(n_bins) * sampling_rate / n_bins
+
+    group_delays = [1000, 750, 500, 250]
+
+    # time signal:
+    time = np.zeros((2, 2, n_samples))
+    time[0, 0, group_delays[0]] = 1
+    time[0, 1, group_delays[1]] = 1
+    time[1, 0, group_delays[2]] = 1
+    time[1, 1, group_delays[3]] = 1
+
+    # create a mock object of Signal class to test the plot independently
+    signal_1 = mock.Mock(spec_set=Signal(time, sampling_rate))
+    signal_1.time = np.squeeze(time[0, 0])
+    signal_1.cshape = (1, )
+    signal_1.sampling_rate = sampling_rate
+    signal_1.times = times
+    signal_1.n_samples = n_samples
+    signal_1.frequencies = frequencies
+    signal_1.n_bins = n_bins
+    signal_1.signal_type = 'energy'
+
+    signal_2 = copy.deepcopy(signal_1)
+    signal_2.time = np.squeeze(time[0])
+    signal_2.cshape = (2, )
+
+    signal_3 = copy.deepcopy(signal_1)
+    signal_3.time = time
+    signal_3.cshape = (2, 2)
+
+    return signal_1, signal_2, signal_3, group_delays
 
-#def test_wrap_to_2pi():
-#def test_nextpow2():
 
 @pytest.fixture
 def sine_plus_impulse_mock():
@@ -57,7 +151,6 @@ def sine_plus_impulse_mock():
     """
     n_samples = 2000
     sampling_rate = 4000
-    amplitude_sine = 1
     amplitude_impulse = 1
     idx_impulse = 0
     frequency = 200
diff --git a/tests/test_plot.py b/tests/test_plot.py
index f5c3332c..664e9997 100644
--- a/tests/test_plot.py
+++ b/tests/test_plot.py
@@ -137,6 +137,7 @@ def sine_plus_impulse_mock():
     signal_object.frequencies = frequencies
     signal_object.freq = freq
     signal_object.n_samples = n_samples
+    signal_object.cshape = (1, )
     signal_object.signal_type = 'power'
 
     return signal_object
['tests/test_plot.py', 'tests/test_dsp.py', 'pyfar/dsp/dsp.py']