Почему я получаю ошибку, используя фильтр Баттерворта для низкочастотной фильтрации?

Question

Я пытаюсь проанализировать амплитуду каждой частоты в 1/3 октавы, поэтому я использую много полосовых фильтров Баттерворта. Тем не менее, они работают только на 50 Гц, когда это 3-й порядок. Я хотел бы использовать 6-й порядок, но по какой-то причине я не получил результатов ниже 1 кГц.

[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)') 

Как правильно отфильтровать полосу пропускания 50 Гц с помощью фильтра Баттерворта 6-го порядка?

Answer 1

Полиномиальное представление «ba» БИХ-фильтров может быть весьма чувствительным к ошибкам квантования коэффициентов фильтра, которые могут перемещать полюса за пределы единичного круга и, соответственно, приводить к нестабильному фильтру.Это может быть особенно проблематично для фильтров с узкой полосой пропускания.

Чтобы лучше понять, что происходит, мы можем сравнить предполагаемые местоположения полюсов, полученные с помощью scipy.butter(..., output='zpk'), с эффективными местоположениями полюсов, полученными путем вычисления корней полинома обратной связи.(a коэффициенты).Это можно сделать с помощью следующего кода:

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

, который показывает, что два набора расположены совместно для фильтров с большей пропускной способностью, но местоположения начинают различаться по мере уменьшения пропускной способности, пока ошибка не станет достаточно большой, чтобыкорни выдвигаются за пределы окружности:

Чтобы избежать этой проблемы, можно использовать каскадную реализацию фильтра секций второго порядка:

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

Это должно расширить выход на более низкие частоты, как показано на следующем графике: