Why do I get an error using butterworth filter for low frequency filtering?

Question

I'm trying to analyze the amplitude of each 1/3 octave band frequency, so I'm using many bandpass butterworth filters. However, they only work for 50 Hz when it is a 3rd order. I would like to use a 6th order, but for some reason I got no results below 1 kHz.

[fs, x_raw] = wavfile.read('ruido_rosa.wav')
x_max=np.amax(np.abs(x_raw))
x=x_raw/x_max
L=len(x)

# Creates the vector with all frequencies

f_center=np.array([50.12, 63.10, 79.43, 100, 125.89, 158.49, 199.53, 251.19, 316.23, 398.11, 501.19, 630.96, 794.33, 1000, 1258.9, 1584.9, 1995.3, 2511.9, 3162.3, 3981.1, 5011.9, 6309.6, 7943.3, 10000, 12589.3, 15848.9])
f_low=np.array([44.7, 56.2, 70.8, 89.1, 112, 141, 178, 224, 282, 355, 447, 562, 708, 891, 1120, 1410, 1780, 2240, 2820, 3550, 4470, 5620, 7080, 8910, 11200, 14100])
f_high=np.array([56.2, 70.8, 89.1, 112, 141, 178, 224, 282, 355, 447, 562, 708, 891, 1120, 1410, 1780, 2240, 2820, 3550, 4470, 5620, 7080, 8910, 11200, 14100, 17800])

L2=len(f_center)

x_filtered=np.zeros((L,L2))
for n in range (L2):    
    order=6
    nyq = 0.5*fs
    low = f_low[n]/nyq
    high = f_high[n]/nyq
    b,a = butter(order,[low,high],btype='band')
    x_filtered[:,n] = lfilter(b,a,x)

x_filtered_squared=np.power(x_filtered,2)

x_filtered_sum=np.sqrt(np.sum(x_filtered_squared,axis=0)/L)

pyplot.figure(2)
pyplot.semilogx(f_center,20*np.log10(np.abs(x_filtered_sum)))
pyplot.xlim((50,16000))
pyplot.xlabel('Frequência (Hz)')
pyplot.ylabel('Amplitude (dB)')

How can I properly filter a 50 Hz bandpass with a 6th order butterworth filter?

SleuthEye SleuthEye · Accepted Answer · 2019-05-01T02:08:57

The polynomial "ba" representation of IIR filters can be quite susceptible to filter coefficient quantization errors which can move the poles outside the unit circle and correspondingly result in an unstable filter. This can be particularly problematic for filters with narrow bandwidths.

To better understand what's going on, we can compare the intended pole locations obtained using scipy.butter(..., output='zpk') with the effective pole locations obtained by computing the roots of the feedback polynomial (the a coefficients). This can be done with the following code:

fs = 44100
order = 6
bands = np.array([np.array([2240, 2820]), 
                  np.array([1650,1.2589*1650]), 
                  np.array([1300,1.2589*1300]), 
                  np.array([1120,1410])])
plt.figure()
index = 1
for band in bands:
  low = band[0]/(fs/2)
  high = band[1]/(fs/2)
  b,a=butter(order,[low,high],btype='band',output='ba')
  z,p,k=butter(order,[low,high],btype='band',output='zpk')
  r = np.roots(a)

  t = np.linspace(-np.pi,np.pi,200)
  plt.subplot(2,2,index)
  index += 1
  plt.title("low={:.0f}, high={:.0f} (BW={:.0f}Hz)".format(band[0], band[1], band[1]-band[0]))
  plt.plot(np.cos(t),np.sin(t),'k:')
  plt.plot(np.real(p),np.imag(p),'bx',label='poles')
  plt.plot(np.real(r),np.imag(r),'rs',fillstyle='none',label='roots')

  # Focus on the poles/roots at positive frequencies
  x0 = np.cos(0.5*np.pi*(low+high))
  y0 = np.sin(0.5*np.pi*(low+high))
  delta = 0.05
  plt.ylim(y0-delta,y0+delta)
  plt.xlim(x0-delta,x0+delta)
  plt.grid(True)
  plt.legend()

which shows that the two sets are collocated for filters with the larger bandwidths, but the locations start to differ as the bandwidth decreases until the error is sufficiently large that the roots are pushed outside the unit circle:

To avoid this problem, it is possible to use a cascaded second-order sections filter implementation:

sos = butter(order,[low,high],btype='band',output='sos')
x_filtered[:,n] = sosfilt(sos,x)

This should extend the output to the lower frequencies, as shown in the following graph:

Why do I get an error using butterworth filter for low frequency filtering?

1 Answers