我正在编写一个 python 脚本,以使用非线性最小二乘优化将级联双二阶滤波器拟合到一些所需的滤波器形状。滤波器 S 域传递函数如下所示:
4),则可以分解为 N-1 个二阶低通和一个二阶带通部分的级联,例如:
此滤波器专门用于模拟听觉滤波器形状。所有部分都使用一个相同的 Q,以保持参数的数量尽可能少并避免过度拟合(这样过滤器可以很容易地泛化而无需重新拟合)。我的无约束优化例程(Levenberg-Marquardt 方法)最小化了几个成本函数以找到最佳滤波器参数:Q(品质因数)、N(滤波器阶数)、a(低频斜率)和滤波器增益。它在拉普拉斯域中进行了优化,使参数变化(不同配件之间)的建模更加直接和有意义。
您可以在下面看到该过滤器的实现。其产生的频率响应适合所需的滤波器形状:
def get_filter_frequency_response_s(worN,Q,a,N,gain,fc,fs):
wc = np.pi*fc / (fs/2)
w0=wc/(np.sqrt(((N-1)/(2*N-1))*(1-1/(2*Q**2)))* \
np.sqrt(1+np.sqrt(1+(1/((N-1)**2/(2*N-1)*(1-1/(2*Q**2))**2)))))
wz =a*w0
# LP transfer function
num_LP=[0, 0, w0**2]
den_LP=[1, w0/Q, w0**2]
from scipy import signal
w_LP, h_LP=signal.freqs(num_LP,den_LP,worN=worN)
# BP transfer function
num_BP=[0, w0, wz*w0]
den_BP=[1, w0/Q, w0**2]
w_BP, h_BP =signal.freqs(num_BP,den_BP,worN=worN)
# Cascad 3 LP biquads with 1 BP biquad filter
h_casc=(h_LP**(N-1)*h_BP)*gain
return h_casc
一旦找到最佳拟合的参数,将 S 域 TF 转换为 Z 域二阶部分,然后实现实际滤波器:
wc = np.pi*fc / (fs/2)
w0=wc/(np.sqrt(((N-1)/(2*N-1))*(1-1/(2*Q**2)))* \
np.sqrt(1+np.sqrt(1+(1/((N-1)**2/(2*N-1)*(1-1/(2*Q**2))**2)))))
wz =a*w0
# LP transfer function
num_LP=[0, 0, w0**2]
den_LP=[1, w0/Q, w0**2]
filtz_LP = signal.lti(*signal.bilinear(num_LP, den_LP)) # convert to Z-domain
# BP transfer function
num_BP=[0, w0, wz*w0]
den_BP=[1, w0/Q, w0**2]
filtz_BP = signal.lti(*signal.bilinear(num_BP, den_BP)) # convert to Z-domain
# Convert to SOS form
sos_LP=signal.tf2sos(filtz_LP.num, filtz_LP.den)
sos_BP=signal.tf2sos(filtz_BP.num, filtz_BP.den)
sos_casc=sos_LP # add a first LP filter
for i in range(N-2): # add N-2 LP filters
sos_casc=np.vstack((sos_casc, sos_LP))
sos_casc=np.vstack((sos_casc, sos_BP)) # add the BP filter
y=signal.sosfiltfilt(sos_casc,x)*gain
拟合工作正常,但它确实返回一个非整数值作为过滤顺序。因此,我的问题是:有没有办法将过滤器顺序限制为仅整数?
非常感谢您的时间和帮助!