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

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

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

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

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

时域连续，且具有非周期性的函数，可以进行傅里叶变换，求出连续的非周期的频谱。
$$\begin{array}{rl}X(\omega )& ={\int }_{-\infty }^{\infty }x(t)e^{-j\omega t}dt\\ x(t)& =\frac{1}{2\pi }{\int }_{-\infty }^{\infty }X(\omega )e^{j\omega t}d\omega \end{array}$$

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

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

$$\begin{array}{rl}X(\omega )& =\sum _{-\infty }^{\infty }x[n]e^{-j\omega n}\\ x[n]& =\frac{1}{2\pi }{\int }_{-\pi }^{\pi }X(\omega )e^{j\omega n}d\omega \end{array}$$

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

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

1、模拟滤波器

比如一个二阶系统，其传递函数为：
$$H(s)=\frac{udn{f}^2}{s^2+2\ast udnf\ast dr\ast s+udn{f}^2}=\frac{0s^2+0s+1}{s^2+1s+1}$$

该传递函数的时域微分形式为：
$$\frac{d^{2}y(t)}{d{t}^2}+2\zeta w_n\frac{dy(t)}{dt}+w_n^{2}y(t)=w_n^{2}x(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三个函数：

1. freqs_zpk：基于零极点的模拟频率响应。

   1. worN：频率轴范围。
   2. np.logspace：生成对数序列

2. 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]
   

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

2、数字滤波器

比如一个二阶系统：
$$H(z)=\frac{1}{1-(2r\cos (\theta )z^{-1}+r^2{z}^{-2}}=\frac{z^2}{z^2-(2r\cos (\theta )z+r^2}$$
其单位脉冲响应为：
$$h[n]=r^n\frac{\sin (n+1)\theta }{\sin \theta }u[n]$$
差分方程表示为：
$$y[n]-2r\cos (\theta )y[n-1]+r^{2}y[n-2]=x[n]$$

import numpy as np
from scipy.signal import freqz_zpk,freqz,tf2zpk
import matplotlib.pyplot as pltfs = 2*np.pir = 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"$-3\pi$",r"$-2\pi$",r"$-\pi$","0",r"$\pi$",r"$2\pi$",r"$3\pi$"])
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()

需要注意：

1. freqs_zpk：没有采样率这个概念，worN的单位就是Hz

2. freqz_zpk：有采样率这个概念，fs的默认值为$2\pi$​，此时横坐标的单位为弧度。

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

               jw                 -jw              -jwMjw    B(e  )    b[0] + b[1]e    + ... + b[M]e
   H(e  ) = ------ = -----------------------------------jw                 -jw              -jwNA(e  )    a[0] + a[1]e    + ... + a[N]e
   

4. 弧度和频率换算举例：设置$worN=[-2\pi ,2\pi ]$，如果fs使用默认值$2\pi Hz$，那么实际横坐标的范围为$[-2\pi ,2\pi ]$，即两个周期；如果fs使用$\pi Hz$，那么实际的横坐标范围为$[-4\pi ,4\pi ]$。其中$\pi$弧度对应$fs/2$ Hz.