將二階巴特沃斯轉換爲一階 - 第二部分 -

Question

我試圖將我工作的二階巴特沃斯低通濾波器轉換爲python中的一階，但它給了我很大的數字，就像flt_y_1st [299]： 26198491071387576370322954146679741443295686950912.0。這是我第二次和1階Butterworth：

import math 
import numpy as np 
import matplotlib.pyplot as plt 

from scipy.signal import lfilter 
from scipy.signal import butter 

def butter_lowpass(cutoff, fs, order=1): 
    nyq = 0.5 * fs 
    normal_cutoff = cutoff/nyq 
    b, a = butter(order, normal_cutoff, btype='low', analog=False) 
    return b, a 

def butter_lowpass_filter(data, cutoff, fs, order=1): 
    b, a = butter_lowpass(cutoff, fs, order=order) 
    y = lfilter(b, a, data) 
    return y 

def bw_2nd(y, fc, fs): 
    filtered_y = np.zeros(len(y)) 

    omega_c = math.tan(np.pi*fc/fs) 
    k1 = np.sqrt(2)*omega_c 
    k2 = omega_c**2 
    a0 = k2/(1+k1+k2) 
    a1 = 2*a0 
    a2 = a0 
    k3 = 2*a0/k2 
    b1 = -(-2*a0+k3) 
    b2 = -(1-2*a0-k3) 

    filtered_y[0] = y[0] 
    filtered_y[1] = y[1] 

    for i in range(2, len(y)): 
     filtered_y[i] = a0*y[i]+a1*y[i-1]+a2*y[i-2]-(b1*filtered_y[i-1]+b2*filtered_y[i-2]) 

    return filtered_y 


def bw_1st(y, fc, fs): 
    filtered_y = np.zeros(len(y)) 

    omega_c = math.tan(np.pi*fc/fs) 
    k1 = np.sqrt(2)*omega_c 
    k2 = omega_c**2 
    a0 = k2/(1+k1+k2) 
    a1 = 2*a0 
    k3 = 2*a0/k2 
    b1 = -(-2*a0+k3) 
# b1 = -(-2*a0) # <= Removing k3 makes better, but still not perfect 

    filtered_y[0] = y[0] 

    for i in range(1, len(y)): 
     filtered_y[i] = a0*y[i]+a1*y[i-1]-(b1*filtered_y[i-1]) 

    return filtered_y 

f = 100 
fs = 2000 
x = np.arange(300) 
y = np.sin(2 * np.pi * f * x/fs) 

flt_y_2nd = bw_2nd(y, 120, 2000) 
flt_y_scipy = butter_lowpass_filter(y, 120, 2000, 1) 
flt_y_1st = bw_1st(y, 120, 2000) 

for i in x: 
    print('y[%d]: %6.3f flt_y_2nd[%d]: %6.3f flt_y_scipy[%d]: %6.3f flt_y_1st[%d]: %8.5f' % (i, y[i], i, flt_y_2nd[i], i, flt_y_scipy[i], i, flt_y_1st[i])) 

plt.subplot(1, 1, 1) 
plt.xlabel('Time [ms]') 
plt.ylabel('Acceleration [g]') 
lines = plt.plot(x, y, x, flt_y_2nd, x, flt_y_scipy, x, flt_y_1st) 
l1, l2, l3, l4 = lines 
plt.setp(l1, linewidth=1, color='g', linestyle='-') 
plt.setp(l2, linewidth=1, color='b', linestyle='-') 
plt.setp(l3, linewidth=1, color='y', linestyle='-') 
plt.setp(l4, linewidth=1, color='r', linestyle='-') 
plt.legend(["y", "flt_y_2nd", "flt_y_scipy", "flt_y_1st"]) 
plt.grid(True) 
plt.xlim(0, 150) 
plt.ylim(-1.5, 1.5) 
plt.title('flt_y_2nd vs. flt_y_scipy vs. flt_y_1st') 
plt.show() 

...我刪除了所有[I-2] s，這是一個前饋和反饋。

但是，似乎這還不夠。我想我需要在a0，b1等中改變一些方程。例如，當我從b1中刪除'+ k3'時，我得到這樣一個情節（看起來更好，不是嗎？）：

我不是專業的過濾器，但至少我知道這個第一順序不同於scipy.butter。所以，請幫我找到正確的係數。先謝謝你。

這裏是我的參考：filtering_considerations.pdf

Answer 1

我來回答自己。

最終的係數爲：

omega_c = math.tan(np.pi*fc/fs) 
k1 = np.sqrt(2)*omega_c 
a0 = k1/(math.sqrt(2)+k1) 
a1 = a0 
b1 = -(1-2*a0) 

這裏是如何。 我反編譯他們從scipy.butter，作爲@sizzzzlerz建議（謝謝）。 scipy.butter吐出這些係數：

b: [ 0.16020035 0.16020035] 
a: [ 1.  -0.6795993] 

注意b和一個從我的參考反轉。他們會是：

a0 = 0.16020035 
a1 = 0.16020035 
b0 = 1 
b1 = -0.6795993 

然後，應用這些係數不完整我的公式：

a1 = a0 = 0.16020035 
b1 = -(1-2*a0) = -{1-2*(0.16020035)} = -(0.6795993) 

到目前爲止，一切都很好。 順便說一句：

k1 = 0.2698 
k2 = 0.0364 

所以：

a0 = k2/(1+k1+k2) = 0.0364/(1+0.2698+0.0364) = 0.0279 

...這是遠遠0.16020035。 在這一點上，我消除K2，並把這樣的：

a0 = k1/(1+k1+x) 

當x = 0.4142，我0.16020164。足夠近。

a0 = k1/(1+k1+0.4142) = k1/(1.4142+k1) 

... 1.4142 ...！？我以前見過這個號碼......：

= k1/(math.sqrt(2)+k1) 

現在的情節是這樣的（flt_y_scipy完全被flt_y_1st覆蓋）：

您可以搜索關鍵字「一階」巴特沃斯「低通濾波器」，「開方（2）「等

這是我週日DIY的結束。 ;-)