0
$\begingroup$

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: ununderstandable part
  • 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.

Thank you very much in advance for all your kind guidance

GH

Python output

====

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.xlabel('f [Hz]', labelpad = -30)
   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.xlabel('f [Hz]', labelpad = -30)
   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.xlabel('f [Hz]', labelpad = -30)
   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.xlabel('f [Hz]', labelpad = -30)
   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.xlabel('t [s]', labelpad = -30)
   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.xlabel('t [s]', labelpad = -30)
   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)

plt.subplots_adjust(hspace=0.25)
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()
$\endgroup$
1
$\begingroup$

In general, what you are seeing is normal. Butterworth filters are designed in the analog domain where the frequency axis is unlimited. When discretizing the the freqyency axis stops at Nyquist.

Both Matlab and Python use the bilinear transform to map analog frequency to the digital frequency. That means that the infinite analog frequency maps to the Nqyusit frequency and the upper part of the frequency axis is highly compressed (in frequency). That's why your filters look steeper and steeper as the cutoff frequency increases.

Whether that's what you want or not, depends on your application. If you don't like this approach you could try the impulse invariant method by just sampling the impulse response an analog filter but that has it's own set of pros and cons.

$\endgroup$
  • $\begingroup$ It's also worth noting that going from the analog (prototype) filter to the digital filter implementation as a design process can be kind of backwards: If OP wants a digital filter with specific characteristics, it would be wise to design it as such, and not take the detour through e.g. Butterworth, which is simply a design for analog, not digital filters. $\endgroup$ – Marcus Müller Oct 3 at 11:21
  • $\begingroup$ Thanks! 2 questions then, if I may: - Where could I find more info on where this compression (or lower resolution) in the upper-frequency range comes from? - Is there a rule of thumb on what value the normalized parameter is supposed to take? Assume please an actual IIR implementation (for 2nd order LP and HP Butterworth). I mean something that says "for most intents and purposes a Wn <= 0.1 should suffice", where "suffice" would mean something like the lack of precision would be in par with that caused by the numerical calculations involved (e.g. over the common IEEE float) $\endgroup$ – user45441 Oct 3 at 11:30
  • $\begingroup$ First question: that's the relationship between discrete-time signals and analog signals. So, your average EE "signals and signals" or "digital signal processing" textbook. $\endgroup$ – Marcus Müller Oct 3 at 11:56
  • $\begingroup$ @user45441 second question: clear indication you really shouldn't be doing a butterworth design and transform that to digital, but directly go for a digital filter design. $\endgroup$ – Marcus Müller Oct 3 at 11:56
  • 1
    $\begingroup$ @user45441: to emphasize Marcus' point: before we can give you any design guidelines, you need to describe what exactly you want your filter to do and not to do. $\endgroup$ – Hilmar Oct 3 at 13:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.