Bilinear transform on Bessel filter ruins the flat group delay property above fs/4

bessel_group_delay.py

"""
Created on Fri Feb 15 14:17:33 2013
"""

from numpy import pi, log10, abs, logspace, diff, unwrap, angle
from scipy import signal
import matplotlib.pyplot as plt

fc    = 300   # Hz
N     = 6
fs    = 2000  # Hz
f_min = 10    # Hz
f_max = 1000  # Hz

f = logspace(log10(f_min), log10(f_max), 200)

# Analog
plt.subplot(2, 2, 1)

# Bessel defined to have same asymptotes as Butterworth of the same order
b, a = signal.butter(N, 2*pi*fc, 'low', analog=True)
w, h = signal.freqs(b, a, 2*pi*f)
plt.plot(f, 20 * log10(abs(h)), color='silver', ls='dashed')

b, a = signal.bessel(N, 2*pi*fc, 'low', analog=True)
w, h = signal.freqs(b, a, 2*pi*f)
plt.plot(f, 20 * log10(abs(h)))
plt.xscale('log')
plt.title('Analog Bessel frequency response')
plt.ylabel('Amplitude [dB]')
plt.ylim(-40, 10)
plt.grid(True, which='both', axis='both')
plt.axvline(fc, color='green')  # cutoff frequency

ax = plt.subplot(2, 2, 3)
group_delay = -diff(unwrap(angle(h)))/diff(2*pi*f)  # rad/(rad/s) = s
plt.plot(f[1:], group_delay)
plt.xscale('log')
plt.title('Analog Bessel group delay')
plt.xlabel('Frequency [radians / second]')
plt.ylabel('Group delay [seconds]')
plt.margins(0, 0.2)
plt.grid(True, which='both', axis='both')
plt.axvline(fc, color='green')  # cutoff frequency


# Digital
plt.subplot(2, 2, 2)

# Bessel defined to have same asymptotes as Butterworth of the same order
b, a = signal.butter(N, fc/(fs/2), 'low')
w, h = signal.freqz(b, a, 2*pi*f/fs)
plt.plot(f, 20 * log10(abs(h)), color='silver', ls='dashed')

b, a = signal.bessel(N, fc/(fs/2), 'low')
w, h = signal.freqz(b, a, 2*pi*f/fs)
plt.plot(f, 20 * log10(abs(h)))
plt.xscale('log')
plt.title('Digital Bessel frequency response')
plt.ylabel('Amplitude [dB]')
plt.ylim(-40, 10)
plt.grid(True, which='both', axis='both')
plt.axvline(fs/4, color='red')  # not supposed to let fc higher than this
plt.axvline(fc, color='green')  # cutoff frequency
if f_max > fs/2:
    plt.axvline(fs/2, color='orange')  # Nyquist frequency

plt.subplot(2, 2, 4)
group_delay = -diff(unwrap(angle(h)))/diff(2*pi*f)
plt.plot(f[1:], group_delay)
plt.xscale('log')
plt.title('Digital Bessel group delay')
plt.xlabel('Frequency [radians / second]')
plt.ylabel('Group delay [seconds]')
plt.ylim(*ax.axes.get_ylim())
plt.grid(True, which='both', axis='both')
plt.axvline(fs/4, color='red')  # not supposed to let fc higher than this
plt.axvline(fc, color='green')  # cutoff frequency
if f_max > fs/2:
    plt.axvline(fs/2, color='orange')  # Nyquist frequency

# impulse invariant