python 5阶切比雪夫滤波器

%matplotlib inline
import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt
from matplotlib import patches

worN = 512

## 获取角度，转换成 360 度形式
def angle(h):
    return np.unwrap(np.angle(h)) * 360 / 2 / np.pi


## 画零极点图（matlab 使用 zplane 函数）
def plot_plane(ax, z, p):
    unit_circle = patches.Circle((0,0), radius=1, fill=False, 
                                 color='black', ls='solid', alpha=0.1)
    ax.add_patch(unit_circle)
    poles = plt.plot(p.real, p.imag, 'x', markersize=9, alpha=0.5)
    zeros = plt.plot(z.real, z.imag, 'o', markersize=9, 
             color='none', alpha=0.5,
             markeredgecolor=poles[0].get_color())
    plt.axis('scaled')
    r = 1.1 * np.amax(np.concatenate((abs(z), abs(p), [1])))
    plt.axis([-r, r, -r, r])
                     
    
## 切比雪夫1型低通滤波器，增益为1.0, 截止频率为 0.4
B, A = signal.cheby1(5, 1.0, 0.4)

print('B = ', B)
print('A = ', A)

## 频率响应和相位响应
w, H = signal.freqz(B, A, worN=worN, whole=False)


## 群延迟
w1, GD =  signal.group_delay((B, A), w=worN, whole=False)


## 获取零极点
z, p, k = signal.tf2zpk(B, A)

## 频率归一化
f = w/np.pi


plt.figure(figsize=(10,6))
ax = plt.subplot(221)
plt.title('零极点')
plot_plane(ax, z, p)
plt.plot()

plt.subplot(222)
plt.title('频率响应')
plt.xlabel('归一化频率(w/pi)')
plt.ylabel('振幅(dB)')
plt.plot(f, np.abs(H))

plt.subplot(223)
plt.title('相位响应')
plt.xlabel('归一化频率(w/pi)')
plt.ylabel('相位(角度)')
plt.plot(f, angle(H))

plt.subplot(224)
plt.title('群延迟')
plt.xlabel('归一化频率(w/pi)')
plt.ylabel('群延迟(采样点)')
plt.plot(f, GD)

plt.subplots_adjust(hspace=0.6)
plt.show()