数字滤波器和模拟滤波器（一）

技术分享 2年前 (2024-06-10) 0 999+

模拟滤波器和数字滤波器（一）

下面介绍模拟滤波器和数字滤波器的频率响应的异同，以及如何使用python地scipy.signal来绘制其频谱响应和冲激阶跃响应。在第二期将谈到如何设计模拟滤波器和数字滤波器。

在正文之间，应该介绍连续时间傅立叶变换（CTFT）和离散时间傅立叶变换（DTFT）。

CTFT 连续时间信号的傅立叶变换

时域连续，且具有非周期性的函数，可以进行傅里叶变换，求出连续的非周期的频谱。

[Large begin{aligned}X(omega) &= int_{-infty}^infty x(t)e^{-j omega t}dt \ x(t) &= frac{1}{2pi}int_{-infty}^infty X(omega)e^{j omega t}domega end{aligned} ]

DTFT 离散时间信号的傅立叶变换

时域离散，且具有非周期性的函数，可以求出连续的周期的频谱。周期为(2pi)

[Large begin{aligned}X(omega) &= sum_{-infty}^infty x[n]e^{-j omega n} \ x[n] &= frac{1}{2pi}int_{-pi}^pi X(omega)e^{j omega n}domega end{aligned} ]

最大的区别是，连续时间信号的频谱从0到无穷大，离散时间信号的频谱从0到(2pi)

下面将介绍python当中的模拟和数字滤波器。

1、模拟滤波器

比如一个二阶系统，其传递函数为：

[H(s) = frac{udnf^2}{s^2+2*udnf*dr*s+udnf^2} = frac{0s^2+0s+1}{s^2+1s+1} ]

该传递函数的时域微分形式为：

[frac{d^2y(t) }{dt^2} + 2zeta w_n frac{dy(t)}{dt} + w_n^2y(t) = w_n^2x(t) ]

import numpy as np from scipy.signal import freqs_zpk,freqs,tf2zpk import matplotlib.pyplot as plt dr = 1/2          # damping ratio udnf = 1          # undamped natural frequency b = [0,0,udnf**2] a = [1,2*udnf*dr,udnf**2] z,p,k = tf2zpk(b,a) w, h = freqs_zpk(z, p, k, worN=np.logspace(-3, 5, 1000))  fig = plt.figure(figsize=(14,7)) ax1 = fig.add_subplot(1, 1, 1) ax1.set_title('Analog filter frequency response') ax1.semilogx(w, 20 * np.log10(abs(h)), 'b') ax1.set_ylabel('Amplitude [dB]', color='b') ax1.set_xlabel('Frequency [Hz]') ax1.grid(True)  ax2 = ax1.twinx() angles = np.unwrap(np.angle(h,deg=True),period=360) ax2.semilogx(w, angles, 'g') ax2.set_ylabel('Angle [degree]', color='g')  plt.axis('tight') plt.show()

from scipy.signal import impulse，step print(z,p,k) t, y = impulse((z,p,k)) t1, y1 = step((z,p,k)) plt.plot(t,y) plt.plot(t1,y1) plt.legend(["impulse response","step response"]) plt.show()

上面用到scipy.signal三个函数：

freqs_zpk：基于零极点的模拟频率响应。
1. worN：频率轴范围。
2. np.logspace：生成对数序列

freqs：基于有理传递函数的模拟频率响应。在本例中没有用到。尤其注意b、a对应传递函数是正幂。

        b[0]*(jw)**M + b[1]*(jw)**(M-1) + ... + b[M] H(w) = ----------------------------------------------         a[0]*(jw)**N + a[1]*(jw)**(N-1) + ... + a[N]

tf2zpk：传递函数转零极点表示。

2、数字滤波器

比如一个二阶系统：

[H(z) = frac{1}{1-(2rcos(theta)z^{-1}+r^2z^{-2}} = frac{z^2}{z^2-(2rcos(theta)z+r^2} ]

其单位脉冲响应为：

[h[n] = r^nfrac{sin(n+1)theta}{sintheta}u[n] ]

差分方程表示为：

[y[n]-2rcos(theta)y[n-1]+r^2y[n-2] = x[n] ]

import numpy as np from scipy.signal import freqz_zpk,freqz,tf2zpk import matplotlib.pyplot as plt  fs = 2*np.pi  r = 3/4             theta = 45/180*np.pi           b = [1,0,0] a = [1,-2*r*np.cos(theta),r**2] z,p,k = tf2zpk(b,a) w, h = freqz_zpk(z, p, k, worN=np.linspace(-2.5*np.pi,2.5*np.pi,1000),fs=fs)  fig = plt.figure(figsize=(14,7)) ax1 = fig.add_subplot(1, 1, 1) ax1.set_title('Digital filter frequency response') ax1.plot(w, 20 * np.log10(abs(h)), 'b') ax1.set_ylabel('Amplitude [dB]', color='b') ax1.set_xlabel('w(radians)') ax1.set_xticks([-3*np.pi,-2*np.pi,-1*np.pi,0,1*np.pi,2*np.pi,3*np.pi],                [r"$-3pi$",r"$-2pi$",r"$-pi$","0",r"$pi$",r"$2pi$",r"$3pi$"]) ax1.grid(True)  ax2 = ax1.twinx() angles = np.unwrap(np.angle(h,deg=True),period=360) ax2.plot(w, angles, 'g') ax2.set_ylabel('Angle [degree]', color='g')  plt.axis('tight') plt.show()

该仿真波形和奥本海姆的教材上面的波形一致。

print(z,p,k) from scipy.signal import dimpulse, dstep dt = 0.1 t, y = dimpulse((z,p,k,dt), n=50) t1, y1 = dstep((z,p,k,dt), n=50) plt.stem(t,np.squeeze(y),'r') plt.plot(t1,np.squeeze(y1),'bo-') plt.legend(["impulse response","step response"]) plt.show()

需要注意：

freqs_zpk：没有采样率这个概念，worN的单位就是Hz
freqz_zpk：有采样率这个概念，fs的默认值为(2pi)，此时横坐标的单位为弧度。

freqz：使用传递函数绘制频谱响应。在scipy.signal的定义里面，此函数为负幂。

            jw                 -jw              -jwM    jw    B(e  )    b[0] + b[1]e    + ... + b[M]e H(e  ) = ------ = -----------------------------------             jw                 -jw              -jwN          A(e  )    a[0] + a[1]e    + ... + a[N]e