python 滤波器设计
21 April 2018
FIR和IIR滤波器
在介绍滤波器之前，先介绍一下几个采样相关的概念：
- 采样频率：每秒钟采样多少个点，需大于奈奎斯特频率，符号 $f_s$
- 采样宽度：量化目标的位数，音频处理中一般为16位（pyaudio中用几个byte表示）
- 频率分辨率：DFT转换时每个频点间的频率差值，等于采样频率除以参数DFT运算的个数N，即 $f_s / N$
DFT特性：
- 输出相位是奇对称的，输出幅度是偶对称的（对于实数序列来说）
- 输入序列移位，对应频率相移
- 输入是偶函数，则输出相位为0
- 输入是奇函数，则输出振幅为0
- DFT频率相位是相对与该频率点余弦波的相位
- 当输入含有的频率不是DFT频率分辨率倍数时，将发生频谱泄漏
- 对输入序列补零以增加DFT点数，则频域数据插值（没有提高频率分辨率）
- 采样率不变时，增加输入序列的DFT点数，则提高频率分辨率
- 矩形窗的基频等于采样率除以矩形窗的点数（即 $f_s / N_w$）,画图时主瓣和旁瓣个数之和为 $N_w - 1$
DFT频率幅度响应公式(sinc函数)：
其中 $A_0$ 时振幅，N 是DFT点数, k 是实际频率对应频率分辨率的倍数，为整数时泄漏为0，非整数时有泄漏。
FIR滤波器公式为：
IIR 滤波器公式为：
增益（dB）与振幅（A）公式：
经验分享：
- FIR 低通滤波器设计方法：用理想低通滤波器系数 sin(x)/x 与所选窗函数相乘得到最终 h(k)
- FIR 带通滤波器设计方法：低通滤波器系数乘以正弦信号，以平移低通滤波器达到带通滤波器对目的
- FIR 高通滤波器设计方法：低通滤波器系数乘以 fs/2 对正弦信号（1，-1 交替对序列）
- 时域相乘频域卷积（时域乘以正弦信号后相当于频域移动 — 调频）
- 增加FIR滤波器系数，可以减小过渡带宽
- 对FIR系数加窗，可以减小通带抖动，代价是过渡带宽有所增加
- 频谱响应功率下降一半时对频点称之为截止频率
- $0.7 = \sqrt{0.5}$ 由于功率与振幅对平方成正比，功率下降一半时，振幅下降 0.7
- -3dB = 20 * log10(0.7) 也就是说振幅下降至0.7时，增益下降3dB到达截止频率
- -6dB = 20 * log10(0.5) 也就是说振幅下降一半时，增益下降6dB
最优滤波器设计法
最优滤波器设计法又称为 Parks-McClellan FIR 方法，是当前现实应用中最广为接受对滤波器设计法。 python 中使用 scipy.signal.remez 函数利用最优滤波器设计方法生滤波器系数。
Python中的滤波器函数为：signal.lfilter。
低通滤波器
滤波器系统可以通过手动计算或者使用现成的函数来生成：
- signal.firwin : 生成低通滤波器的 FIR 系数
- signal.remez : 生成带通滤波器的 FIR 系数
滤波器采样点跟频率的关系：
%matplotlib notebook
import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt
def lowpass_filter(n, fc):
return 2 * fc * np.sinc(2 * fc * np.arange(-n, n, 1.0))
# 截止频率 0.1 * fs
b1 = lowpass_filter(30, 0.1)
b2 = signal.firwin(len(b1), 0.2)
# 频率响应
w1, h1 = signal.freqz(b1)
w2, h2 = signal.freqz(b2)
# 画图
plt.figure(figsize=(12,9))
plt.subplot(2,1,1)
plt.plot(w1/2/np.pi, 20*np.log10(np.abs(h1)+0.01), label = 'ideal')
plt.plot(w2/2/np.pi, 20*np.log10(np.abs(h2)+0.01), label = 'firwin')
plt.legend()
plt.ylabel(u"幅值(dB)")
plt.title('低通滤波器频响')
plt.subplot(2,1,2)
plt.title('滤波器系数')
plt.plot(b1, label='ideal coef')
plt.plot(b2, label='firwin coef')
plt.legend()
plt.show()
均值低通滤波器
均值低通滤波器在现实生活中还是比较常用的，我们一般仅仅是根据常识来使用均值滤波器，却并不清楚其频响特性。 这里先总结一下关于均值滤波器的几个结论：
- 均值滤波器是低通滤波器，与矩形窗类似，其旁瓣增益很高
- 均值滤波器的个数越多，其截止频率越低
下面我们用python代码画出均值滤波器的频响曲线及其滤波效果。
%matplotlib notebook
import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt
sample_rate = 8000
time = 0.5
# 均值滤波器系统 FIR
b = np.array([0.25, 0.25, 0.25, 0.25])
a = np.array([1.0])
# 生成 0.5s 的音频数据，采样率 8k，频率从 50 ~ 4000
t = np.arange(0, time, 1.0/sample_rate)
x = signal.chirp(t, 50, time, 4000)
# 滤波器处理
y = signal.lfilter(b, a, x, zi=None)
# 计算滤波器的脉冲响应
impulse = np.zeros(7)
impulse[0] = 1
imp_y = signal.lfilter(b, a, impulse)
# 计算频率响应
w, h = signal.freqz(b,a)
# 画图
plt.figure(figsize=(12,6))
plt.subplot(2,2,1)
plt.title('脉冲响应')
plt.plot(imp_y)
plt.subplot(2,2,2)
plt.title('50～4000Hz数据')
plt.plot(t, x)
plt.subplot(2,2,3)
plt.title('频率响应')
plt.plot(w/2/np.pi, 20*np.log10(np.abs(h)+0.01))
plt.subplot(2,2,4)
plt.title('滤波后数据')
plt.plot(t, y)
plt.show()
半带滤波器
半带滤波器因为其稀疏的非零系数使其成为重采样因子为2的整数次幂（2，4，8等）的采样率转换应用对理想滤波器。 半带 FIR 滤波器是频率响应关于点 $f_s/4$ 对称的一种滤波器, 它有以下特点：
- $f_{pass} + f_{stop} = f_s / 2$
- $\Delta f = f_s / 4 - f_{pass} = f_{stop} - f_s / 4$
- 通带和组带抖动相同
- 其 h(k) 的系数交替为零
- 抽头系统 N + 1 必须是4对整数倍
- 所需的乘法器数量大概为 N / 4，因为其中接近一般为0的系数以及系数的对称性
如果使用最优滤波器设计法（Parks-McClellan FIR），则 $\Delta f = f_s /(N + 1)$
下面我们设计一个 11 个抽头的半带滤波器，因为 $\Delta f = \frac{1}{11 + 1} = 0.08$， 所以归一化的半带滤波器的两个频带：
$f_{pass} = 0.25 - 0.08 = 0.17$
$f_{stop} = 0.25 + 0.08 = 0.33$
%matplotlib notebook
import matplotlib.pyplot as plt
import scipy.signal as signal
import numpy as np
def amp2db(amp):
return 20 * np.log10(np.clip(np.abs(amp/abs(amp).max()), 1e-5, 1e100))
## 使用最有滤波器设计法得到半带滤波器系数
N = 11
taps_origin = signal.remez(N, [0.0, 0.17, 0.33, 0.5], [1,0])
taps = taps_origin.copy()
taps[abs(taps) <= 1e-4] = 0.0
## 打印半带滤波器系数
print(' origin taps ')
print('------------------------------------')
for i in range(N):
print(' index %2d %10.6f %10.6f' % (i+1, taps_origin[i], taps[i]))
## 生成频率响应
W, H = signal.freqz(taps, whole=True)
plt.figure(figsize=(12,9))
plt.subplot(211)
plt.title('11抽头半带滤波器系数')
plt.stem(np.arange(len(taps)), taps)
plt.grid()
plt.subplot(212)
plt.title('11抽头半带滤波器频响')
plt.plot(W/2/np.pi, amp2db(H))
plt.axvspan(0.17, 0.33, facecolor='0.5',alpha=0.7)
plt.axvline(0.25, color='g', linestyle='--')
plt.grid('on')
plt.show()
origin taps
------------------------------------
index 1 0.027289 0.027289
index 2 0.000065 0.000000
index 3 -0.078655 -0.078655
index 4 -0.000009 0.000000
index 5 0.308248 0.308248
index 6 0.500045 0.500045
index 7 0.308248 0.308248
index 8 -0.000009 0.000000
index 9 -0.078655 -0.078655
index 10 0.000065 0.000000
index 11 0.027289 0.027289
从上图可以看出，11个系数的半带滤波器有4个为0，剩余的7个系数关于第6个系数(0.500045 约等于 0.5)对称，又由于 0.5 可以使用移位实现，因此只需3个乘法器就可以实现11抽头的半带滤波器。实践证明半带滤波器在降2和升2采样率转换应用中比较好用。