Possible error when finding group/phase delay in scipy

Question

After designing a "least square optimal" FIR filter I wanted to find the group/phase delay, defined as: $P(\omega) \triangleq - \frac{\Theta(\omega)}{\omega}$ (Depending on the literature group and phase are used inconsistently). In scipy I've done this as follows:

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt


# Filter design parameters 
fs = 25.0                   # Hz
desired = (1, 1, 0, 0)      # Ideal filter "bands"
bands = (0, 1, 5, 12.5)     # Frequency at band edges
taps = 15                  # Filter order.

# Create LS optimal filter
fir_coeffs = signal.firls(taps, bands, desired, fs=fs)   
normfreq, response = signal.freqz(fir_coeffs)
freq = 0.5*fs*normfreq/np.pi #Convert form rad/sample to Hz

# Find grp delay
normgrpfreq, grpdly = signal.group_delay((fir_coeffs,1), w=512, whole=False)
grpfreq = 0.5*fs*normgrpfreq/np.pi #Convert form rad/sample to Hz
zero, pole, gain = signal.tf2zpk(fir_coeffs, 1)

fig, axs = plt.subplots(3)
axs[0].set(title='Amplitude Response')
axs[0].set_xlabel('Frequency (Hz)')
axs[0].set_yscale("log")
axs[0].set_ylabel('Magnitude (DB)')
axs[0].grid(True)
axs[0].plot(freq, np.abs(response))     # Plot with log y axis


axs[1].set(title='Phase Response')
axs[1].set_xlabel('Frequency (Hz)')
axs[1].set_ylabel('Phase (°)')
axs[1].grid(True)
axs[1].plot(freq, np.angle(response)*180/np.pi)     # Plot with log y axis

axs[2].set(title='Phase Delay')
axs[2].set_xlabel('Frequency (Hz)')
axs[2].set_ylabel('Delay (s)')
axs[2].grid(True)
axs[2].plot(grpfreq, grpdly/fs)     # Plot with log y axis 

fig.tight_layout()

resulting in:

The calculated delay seems wrong. If we limit ourselves to looking at the approximately linear part of the phase plot (0 to ~2 Hz) I would expect the delay in this period to be close to a constant $(120/2)* (\pi/ 180)\approx 1\ \mathrm{samples}$.

Am I misunderstanding something or doing something wrong?

EDIT: Further more when I convert the system to a minimum phase system, thereby ensuring minimum group delay the group delay seem to increase to roughly two samples, which doesn’t make sense either..

fir_coeffs =signal.minimum_phase(fir_coeffs, method='homomorphic', n_fft=None)

Answer 1

The function firls designs a linear phase filter, so there's no doubt about the phase delay or group delay. Both delays are equal in that case, and their value must be $(N-1)/2$, where $N$ is the filter length (number of taps). Note that the input parameter taps is not the filter order (as noted in your code), but the number of taps. With $N=15$ we end up with a phase delay and group delay of exactly $7$ samples. And this is also what I see when I look at your phase plot: at a frequency of about $1.8$ the phase equals $-180$ degrees ($-\pi$ radians), which means that the delay is approximately

$$\tau_g\approx-\frac{\pi}{2\pi \cdot 1.8/f_s}=\frac{25}{3.6}=6.94\textrm{ samples}$$

Note that the phase is indeed linear. There are two kinds of jumps in the plot: the artificial ones due the principal value of the phase (jumps by $2\pi$ are no jumps at all), and the real ones due to the zeros of the frequency response in the stopband (jumps by $\pi$). Note that these jumps occur exactly where the frequency response has stopband zeros. Any linear phase filter with stopband zeros has those phase jumps (where the group delay is undefined), but that doesn't change the fact that the phase is linear.

The delay plot should show a value of $7$ samples (and I'm not sure if it does or does not), except at the frequencies where the frequency response has zeros. At those frequencies the delay is not defined (one could argue that the delay is infinite because a sinusoid at such a frequency never reaches the output). The actual delay values at the zeros of the frequency response returned by the software are just a matter of definition (usually zero).

I would say that the phase and delay plots of the minimum phase filter look reasonable.

Answer 2

A few things are happening here (more of a comment but too long for squeeze it in there).

$P(\omega) \triangleq - \frac{\Theta(\omega)}{\omega}$ (Depending on the literature group and phase are used none consistently)

No. The literature about this is quite consistent. Phase delay is phase divided by frequency, group delay is the first derivative of the phase vs frequency, i.e.

$$ p_D(\omega) = -\frac{\varphi(\omega)}{\omega} \\ g_D(\omega) =-\frac{\partial \varphi(\omega)}{\partial \omega} $$

Since you are designing a linear phase filter both are the same, but that's ONLY true for linear phase filters.

axs[1].plot(freq, np.angle(response)*180/np.pi)

Always unwrap your phase before plotting it. It definitely looks wrong but without unwrapping it's hard to see.

I would expect the delay in this period to be close to a constant (120/2)∗(π/180)≈1 samples.

No. This is a linear phase FIR filter and both group and phase delay should be half half the filter order and constant at all frequencies. In your case you should expect a delay of 7 samples.