How to find the group delay of IIR filters in Python/SciPy

Question

I sent a square wave signal (red) through a IIR filter (butterworth in this case), to make it more "realistic". Then there's a slight group delay in my output signal (blue). Is there a way to find out what the group delay is? I only need the total group delay, not spectrum of group delay. I need the number so I can shift my decision point accordingly, so I am still sampling at the middle of the eye diagram to achieve optimum result.

Answer 1

I only need the total group delay, not spectrum of group delay.

Group delay is a spectrum, so this doesn't make sense. The group delay is the derivative of the phase response of the filter, so in Python it can be calculated as

from scipy import signal
from numpy import pi, diff, unwrap, angle

w, h = signal.freqs(b, a)
group_delay = -diff(unwrap(angle(h))) / diff(w)

The units of phase are radians, and the units of frequency are radians/sec, so the units of group delay are in seconds. Examples at 1 and 2

The group delay of a Bessel filter is maximally flat, so when you say "the" group delay, you probably mean the maximum passband group delay. Like in this filter:

The passband has about 0.032 seconds of group delay. Your Butterworth filter does not attempt to have a flat group delay in the passband, however, so it's up to you to decide what "the" group delay is:

Do you want the group delay at 0 Hz, or the maximum group delay at the cutoff frequency? Or maybe at a specific frequency like the fundamental clock frequency of your digital signal?

Answer 2

For future reference, here is same code as @MattL's answer, but in Python:

import numpy as np, matplotlib.pyplot as plt
from scipy import signal
fc = 0.2
b, a = signal.butter(15, fc)
n = np.arange(500)
c = 0.0005
x = np.exp(-c(n-250) * 2)
H, w = signal.freqz(b, a, 4096)
W, gd = signal.group_delay((b, a), 4096)
w0 = .92 * np.pi * fc  # carrier frequency
y = x * np.cos(w0*n)   # modulated signal
z = signal.lfilter(b,a,y)
I = np.argmin([abs(ww-w0) for ww in W])
tau = int(gd[I])       # tau
plt.plot(W, gd)
plt.show()
plt.subplot(2,1,1)
plt.plot(n,y)
plt.subplot(2,1,2)
plt.plot(n,z)
plt.plot(n[:-tau],z[tau:], '--')
plt.show()

Edit: I tried with a bandpass, and I have two cases:

Working:

import numpy as np, matplotlib.pyplot as plt
from scipy import signal
fc = 0.02
fc1, fc2 = 0.019, 0.021
b, a = signal.butter(2, [fc1, fc2], btype='bandpass')
n = np.arange(5000)
x = (abs(n-2000) < 100)
H, w = signal.freqz(b, a, 50000)
W, gd = signal.group_delay((b, a), 50000)
w0 = 1.0 * np.pi * fc

y = x * np.cos(w0*n)

z = signal.lfilter(b,a,y)
I = np.abs(W-w0).argmin()
tau = int(gd[I]) 
if tau == 0:            # this happens when scipy has problem to compute group delay and puts 0
    tau = int(max(gd))  # then we use the max of group delay, so the shifting 
                        # should sometimes be too much, but never not enough
plt.subplot(2,1,1)
plt.plot(n,y)
plt.subplot(2,1,2)
plt.plot(n,z)
plt.plot(n[:-tau],z[tau:], '--')
plt.show()

Non-workingWorking too:

Same code with a narrower bandwidth:

fc = 0.0200
fc1, fc2 = 0.0197, 0.0203

As you see, here it failed!

The reason is simply that scipy's group_delay function is buggy here (confirmed issue).

Let's compute the group delay by hand and it will work:

import numpy as np, matplotlib.pyplot as plt
from scipy import signal
fc = 0.0200
fc1, fc2 = 0.0197, 0.0203
b, a = signal.butter(2, [fc1, fc2], btype='bandpass')
n = np.arange(10000)
x = (abs(n-2000) < 100)
w, H = signal.freqz(b, a, 50000)
gd = -np.diff(np.angle(H))/np.diff(w)
w0 = 1.0 * np.pi * fc

y = x * np.cos(w0*n)

z = signal.lfilter(b,a,y)
i = np.abs(w-w0).argmin()
tau = int(gd[i])
plt.subplot(2,1,1)
plt.plot(n,y)
plt.subplot(2,1,2)
plt.plot(n,z)
plt.plot(n[:-tau],z[tau:], '--')
plt.show()

Answer 3

The original question is basically about rise time; it is only indirectly related to group delay. There are many definitions of rise time, but they all depend directly on the system's step response, or, equivalently, the system's impulse response. The most straightforward way to solve the original problem is to directly determine the rise time from the filter's step response and use it to delay the decision, if necessary. As far as I know, there is no simple way to determine the rise time directly from the filter coefficients (without computing the step or impulse response). An approximation of the rise time $t_r$ is possible by relating it to the inverse of the (low pass) filter's cut-off frequency $f_c$: $$t_r\approx\frac{c}{f_c}$$ where the constant $c$ needs to be determined experimentally.

The question asked by user Basj is different in the sense that it is not related to rise time but directly to group delay. The envelope of a modulated low pass signal is delayed by the filter's group delay evaluated at the carrier frequency. So if the input signal is given by

$$x[n]=s[n]\cos(\omega_0n)\tag{1}$$

where $s[n]$ is a low pass signal, then the output of the filter is approximately given by

$$y[n]\approx |H(e^{j\omega_0})|\cdot s[n-\lfloor\tau_g(\omega_0)\rfloor]\cdot\cos(\omega_0n+\phi(\omega_0))\tag{2}$$

where $H(e^{j\omega})=|H(e^{j\omega})|e^{j\phi(\omega)}$ is the filter's frequency response, $\phi(\omega)$ is the phase, and $\tau_g(\omega)$ is the group delay (which is the negative derivative of the phase with respect to frequency).

Equation $(2)$ is easily verified by the following simple Octave/Matlab example:

fc = .2;
[b,a] = butter(15,fc);
n = 0:500;
c = .0005;
x = exp(-c*(n-250).^2);

[H,w] = freqz(b,a,4096);
gd = grpdelay(b,a,4096);

w0 = .92 * pi*fc;       % carrier frequency
y = x .* cos(w0*n);     % modulated signal
z = filter(b,a,y);

% find filter's magnitude and group delay at the carrier frequency
[dum,I] = min(abs(w - w0));
H = abs(H(I));
gd = round(gd(I));      % 29 samples

xshift = H*[zeros(1,gd),x(1:end-gd)];       % scaled and shifted envelope

subplot(2,1,1), plot(n,y,n,x,'g',n,-x,'g'), grid on
subplot(2,1,2), plot(n,z,n,xshift,'g',n,-xshift,'g'), grid on

Clearly, the envelope of the output signal (bottom figure) is delayed by the group delay at the carrier frequency, which equals $29$ samples. Note that the envelope of the output signal was not computed from the output signal, but it was predicted from Eq. $(2)$ (as shown in the code).