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

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. • What do you mean by "total group delay"? Group delay is a function of frequency, and in general for IIR filters, it is not a constant (since they don't generally have linear phase). – Jason R Apr 4 '13 at 13:12
• @JasonR, Yes, generally speaking the group delay of a filter is a function of frequency. But over all, when I send a signal through the filter, it is delayed by a certain amount, as you can see in the plot, the blue is delayed vs the red. And that shifts my optimum decision point for demodulation. Spectrum always means something in time domain, right? like integral of phase noise spectrum gives timing jitter. The group delay spectrum should relate to time domain delay that I saw. – LWZ Apr 4 '13 at 15:18
• @endolith, Does the group delay relates to rise/fall time in any way? – LWZ Apr 13 '13 at 19:50
• @LWZ: I don't know this stuff that well. I know that Gaussian filters have minimal overshoot, minimal rise/fall time, and minimal group delay. I know that Bessel filters have maximally-flat group delay, which also means their overshoot is small, but not as small as Gaussian filters. I think Bessel filters become more and more like Gaussian filters the more poles you add? I guess I don't really know what a Gaussian filter is, though. – endolith Jun 21 '13 at 1:16
• Why don’t you just make an eye diagram and adjust your sample location to the minimum ISI point (once equalized if that is needed)? – Dan Boschen Sep 26 '18 at 20:25

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?

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).

• Here is your code translated in Python: dsp.stackexchange.com/a/52265/5648 – Basj Sep 28 '18 at 12:04
• Thank you for your answer, I'm just not sure to understand one thing: what makes the question in the bounty (=> group delay) different to OP's question (=> rise time)? Is the former about delay of one specific frequency sinusoid, and the latter about delay of a heavyside step function? – Basj Sep 28 '18 at 12:10
• @Basj: Yes, group delay is a function of frequency, and its interpretation as the delay of the envelope makes only sense for narrow band input signals. The original question is about the rise time of the step response, which cannot be translated to group delay in any obvious way. – Matt L. Sep 28 '18 at 13:40

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() • It's interesting to know that scipy's group_delay function has a bug. Where has this been confirmed? Do you know what the bug is? – Matt L. Sep 29 '18 at 6:03
• @MattL. Here are some more informations. – Basj Sep 29 '18 at 7:31