在Python中使用双线性插值从连续到离散

Question

我在频域中有一个一阶低通滤波器(LPF(，我想将其数字化。我正在比较频率响应图进行测试，但我得到了奇怪的结果。。。

虽然很基本，但我无法通过阅读scipy.signal.bilinear帮助页面或网络来正确理解。

• num，den：S平面上的分子和分母
• b，a：我希望它是数字差分方程滤波器(IIR(的b，a系数，形状为：Y[n]=a0*X[n]+a1*X[n-1]+…-b1*Y[n-1]

下面是一个代码示例：

Fs = 48000.0
f = 2 * np.logspace(1,4,1024)
num =  [0 , 1]
den = [0.001 , 1]
tmp, H = sig.freqs(num, den, worN=1024)
b, a = sig.bilinear(num, den, 1.0)
tmp, Hd = sig.freqz(b,a, worN=1024)
plt.semilogx(f, 20*np.log10(np.abs(H)))
plt.semilogx(f, 20*np.log10(np.abs(Hd)))

我做错了什么？

Answer 1

问题是在绘制时没有在x轴上使用tmp，而且freqz给出了以弧度/sample:为单位的归一化tmp矢量

import numpy as np
import scipy.signal as sig
import matplotlib.pyplot as plt
Fs = 48000
num =  [0 , 1000]
den = [1 , 1000]
w1, H = sig.freqs(num, den, worN=1024)
b, a = sig.bilinear(num, den, Fs)
w2, Hd = sig.freqz(b, a, worN=1024)
fig = plt.figure()
plt.title('Filter frequency response')
plt.semilogx(w1, 20*np.log10(np.abs(H)),'b')
plt.semilogx(w2*Fs, 20*np.log10(np.abs(Hd)),'k')
plt.ylabel('magnitude [dB]')
plt.xlabel('frequency [Hz]')
plt.grid()
plt.axis('tight')
plt.xlim([0.001, Fs/2])
plt.show()

这个代码运行得非常好。希望能有所帮助。

Answer 2

是的，这非常有帮助。然后我进入下一步，使两个"w"向量的长度相同，并覆盖所有有趣的采样点，如下所示：

Fs = 48000.0
f = 2 * np.logspace(1,4,1024)
w = 2 * np.pi * f
Snum =  [0 , 1]
Sden = [0.001 , 1]
w1, H = sig.freqs(Snum, Sden, worN=w)
b, a = sig.bilinear(Snum, Sden, Fs)
w2, Hd = sig.freqz(b,a, worN=w1/Fs)

此外，如果所示图形以Hz为单位(而不是Rad/Sec0(，则w应除以2*PI:

plt.semilogx(w1/np.pi/2, 20*np.log10(np.abs(H)), 'r')