将二阶 Butterworth 转换为一阶 - 第二部分 -
Convert 2nd order Butterworth to 1st order - Part II -
我正在尝试将我工作的二阶巴特沃斯低通滤波器转换为 python 中的一阶,但它给了我非常大的数字,例如 flt_y_1st[299]:26198491071387576370322954146679741443295686950912.0 .这是我的二阶和一阶 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],它们是前馈和反馈。
然而,这似乎还不够。我想我需要更改 a0、b1 等中的一些方程式。例如,当我从 b1 中删除“+k3”时,我得到这样的图(看起来更好,不是吗?):
我不是专门研究过滤器的,但至少我知道这个一阶不同于 scipy.butter。所以,请帮我找到正确的系数。提前谢谢你。
这是我的参考:filtering_considerations.pdf
让我自己回答。
最终系数为:
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)
方法如下。
正如@sizzzzlerz 所建议的(谢谢),我从 scipy.butter 对它们进行了逆向工程。
scipy.butter 吐出这些系数:
b: [ 0.16020035 0.16020035]
a: [ 1. -0.6795993]
请注意,b 和 a 与我的参考相反。他们将是:
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 ...!?我以前见过这个号码...:[=20=]
= k1/(math.sqrt(2)+k1)
现在的情节是这样的(flt_y_scipy 完全被 flt_y_1st 覆盖):
您可以搜索关键字"first order" butterworth "low pass filter" "sqrt(2)"等
我的星期天 DIY 到此结束。 ;-)
我正在尝试将我工作的二阶巴特沃斯低通滤波器转换为 python 中的一阶,但它给了我非常大的数字,例如 flt_y_1st[299]:26198491071387576370322954146679741443295686950912.0 .这是我的二阶和一阶 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],它们是前馈和反馈。
然而,这似乎还不够。我想我需要更改 a0、b1 等中的一些方程式。例如,当我从 b1 中删除“+k3”时,我得到这样的图(看起来更好,不是吗?):
我不是专门研究过滤器的,但至少我知道这个一阶不同于 scipy.butter。所以,请帮我找到正确的系数。提前谢谢你。
这是我的参考:filtering_considerations.pdf
让我自己回答。
最终系数为:
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)
方法如下。 正如@sizzzzlerz 所建议的(谢谢),我从 scipy.butter 对它们进行了逆向工程。 scipy.butter 吐出这些系数:
b: [ 0.16020035 0.16020035]
a: [ 1. -0.6795993]
请注意,b 和 a 与我的参考相反。他们将是:
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 ...!?我以前见过这个号码...:[=20=]
= k1/(math.sqrt(2)+k1)
现在的情节是这样的(flt_y_scipy 完全被 flt_y_1st 覆盖):
您可以搜索关键字"first order" butterworth "low pass filter" "sqrt(2)"等
我的星期天 DIY 到此结束。 ;-)