模擬濾波器和數字濾波器（一）
下面介紹模擬濾波器和數字濾波器的頻率響應的異同，以及如何使用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}{2\pi}\int_{-\infty}^\infty X(\omega)e^{j \omega t}d\omega \end{aligned} \]
- DTFT 離散時間訊號的傅立葉變換
時域離散，且具有非週期性的函式，可以求出連續的週期的頻譜。週期為\(2\pi\)
\[\Large \begin{aligned}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{aligned} \]最大的區別是，連續時間訊號的頻譜從0到無窮大，離散時間訊號的頻譜從0到\(2\pi\)
下面將介紹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} + 2\zeta 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：基於零極點的模擬頻率響應。
- worN：頻率軸範圍。
- 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-(2r\cos(\theta)z^{-1}+r^2z^{-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^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"$-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()
需要注意：
-
freqs_zpk：沒有采樣率這個概念，worN的單位就是Hz
-
freqz_zpk：有采樣率這個概念，fs的預設值為\(2\pi\)，此時橫座標的單位為弧度。
-
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 -
弧度和頻率換算舉例：設定\(worN=[-2\pi,2\pi]\)，如果fs使用預設值\(2\pi Hz\)，那麼實際橫座標的範圍為\([-2\pi,2\pi]\)，即兩個週期；如果fs使用\(\pi Hz\)，那麼實際的橫座標範圍為\([-4\pi,4\pi]\)。其中\(\pi\)弧度對應\(fs/2\) Hz.