Frequency response plot and actual behavior of a digital Butterworth LPF

I have been struggling with a not-too-high (2nd) order digital Butterworth low-pass filter both in SciPy and MatLab (latest of both).

The intention is to apply this filter on signals with different sampling rate, some with a Nyquist frequency very close to the filter's critical frequency, some with that farther away.

When I plotted the frequency response of these filters, using freqz(), I got something else than I expected: - The -3dB point was at the correct location in all cases -> good - However the actual response curve looked different: the closer I got to 1 with the normalized critical frequency Wn, which is in [0, 1), the "sharper" the cuttoff looked like, as if the filter had been much higher (with Wn around 0.99 I got a perfect-looking filter) -> something is obviously wrong here

The problem seems to be two-fold:

• Frequency response plot: with Wn above approx. 0.1 the curve do seem to diverge from the ideal (being 6 dB drop in gain / octave / order, or 20 dB drop in gain / decade / order). Part of it is understandable:
• frequencies above the critical are expected to get attenuated faster, simply because those are beyond the Nyquist frequency
• however I expected the frequencies below the critical to get attenuated (near-)the same way in all cases, but closer Wn gets to 1, the sharper this gets, see the area/curve shown by the red arrow:
• Actual filter behavior: with Wn above approx. 0.8 the filter behavior is incorrect (attenuation is too small or none)

Allow me please to have multiple questions with regards to this:

• Did I make a mistake somewhere? I obviously did (= some fundamental piece of the puzzle is missing), otherwise this would be the best solution for implementing the perfect filter (no attenuation fc <= NF, infinite attenuation for fc > NF)
• If not, where the explanation for these phenomena lie (in the theory)?
• Is the discrepancy in the frequency response real (= communication bus filters indeed works that way) or just the outcome of some underlying mechanism (e.g. numeric precision)?
• Is the discrepancy in the filter behavior real?
• What is the "safe range" (or limit) for Wn, to get a filter in real life that performs as expected for simulations and that for implementations (of IIRs)?

Complete Python code below, if one wishes to run it.

GH

====

import math as m
import numpy as np
import scipy.signal as sps
import matplotlib.pyplot as plt

# fc_all: the -3dB point for all filters [Hz]
fc_all = 1

# srs: sampling rate for filter the filters [Hz]
srs = [2.1, 3, 4, 5, 7, 8, 10, 15, 20, 50, 100]

def do_plot(i, sr):
global sp_411
global sp_412
global sp_425
global sp_426
global sp_427
global sp_428

b_last = (i == (len(srs) - 1))

Wn = (fc_all / (sr / 2))
s_Wn = str(round(Wn, 2))

btl, atl = sps.butter(2, Wn, 'low', analog = False)
wtl, htl = sps.freqz(btl, atl, worN = 65536)
wtx = ((sr * wtl) / (2 * m.pi))

plt.sca(sp_411)
plt.title('Frequency-domain response: log ferquency vs. linear gain')
plt.grid(which='both', axis='both')
plt.xscale('log')
plt.ylabel('gain')
plt.xlim(0.35, 10)
plt.plot(wtx, abs(htl), label='Wn=' + s_Wn)
if b_last:
plt.plot(wtx, np.full(len(wtx), (1 / m.sqrt(2))), label = '-3dB', color = 'red')
plt.legend(loc = 'upper right', ncol = 2)

plt.sca(sp_412)
plt.title('Frequency-domain response: log ferquency vs. log gain')
plt.grid(which='both', axis='both')
plt.xscale('log')
plt.ylabel('gain [dB]')
plt.xlim(0.9, 60)
plt.plot(wtx, (20 * np.log10(abs(htl))))

plt.sca(sp_425)
plt.title('Frequency-domain response: log ferquency vs. log gain (vertically zoomed)')
plt.grid(which='both', axis='both')
plt.xscale('log')
plt.ylabel('gain [dB]')
plt.ylim(-21, 1)
plt.plot(wtx, (20 * np.log10(abs(htl))))
if b_last:
plt.plot(wtx, np.full(len(wtx), -3), color = 'red')

plt.sca(sp_426)
plt.title('Frequency-domain response: log ferquency vs. log gain (horizontally and vertically zoomed)')
plt.grid(which='both', axis='both')
plt.xscale('log')
plt.ylabel('gain [dB]')
plt.ylim(-21, 1)
plt.xlim(0.5, 4)
plt.plot(wtx, (20 * np.log10(abs(htl))))
if b_last:
plt.plot(wtx, np.full(len(wtx), -3), color = 'red')

plt.sca(sp_427)
plt.title('Time-domain response')
plt.ylabel('signal [V]')
x = np.linspace(0, 10, (((10 * sr) / fc_all) + 1), endpoint = False)
y = np.sin(x * (2 * m.pi))
plt.plot(x, sps.lfilter(btl, atl, y), zorder = 1)
if b_last:
plt.plot(x, y, zorder = 0, linestyle='dashed', label = '1V input', color = 'grey', linewidth = 1.5)
plt.plot(x, np.full(len(x), (1 / m.sqrt(2))), label = '+/-1/sqrt(2)', color = 'red')
plt.plot(x, np.full(len(x), (-1 / m.sqrt(2))), color = 'red')
plt.legend(loc = 'upper right', ncol = 2)

plt.sca(sp_428)
plt.title('Time-domain response (horizontally zoomed)')
plt.ylabel('signal [V]')
plt.xlim(0.25, 1.75)
plt.plot(x, sps.lfilter(btl, atl, y), zorder = 1)
if b_last:
plt.plot(x, np.full(len(x), (1 / m.sqrt(2))), color = 'red')
plt.plot(x, np.full(len(x), (-1 / m.sqrt(2))), color = 'red')
plt.plot(x, y, zorder = 0, linestyle='dashed', color = 'grey', linewidth = 1.5)

sp_411 = plt.subplot(4, 1, 1)
sp_412 = plt.subplot(4, 1, 2)
sp_425 = plt.subplot(4, 2, 5)
sp_426 = plt.subplot(4, 2, 6)
sp_427 = plt.subplot(4, 2, 7)
sp_428 = plt.subplot(4, 2, 8)

for i, v in enumerate(srs):
do_plot(i, v)

plt.get_current_fig_manager().window.state('zoomed')
plt.show()