数字滤波器设计是信号处理原理的必考考题,但此题的过程相对来说非常机械,只需要按照步骤照搬即可,这里首先记录滤波器设计的步骤,然后对原理进行简单分析

滤波器设计的问题形状

滤波器设计的问题往往会给出以下参数:

  • 滤波器类型,例如高通、低通、带通等
  • 通带边缘频率
  • 阻带边缘频率
  • 采样频率
  • 阻带衰减,单位为分贝(dB)
  • 窗函数表,给出了可以选择的窗函数的阻带衰减窗内项数等参数

FIR滤波器设计

低通FIR滤波器

步骤如下:

  1. 计算滤波器的截止频率fcf_c,使用通带边缘频率阻带边缘频率的平均值
  2. 计算截止频率对应的数字频率,注意采样频率fsf_s对应的数字频率(角频率)为2π2\pi,因此fcf_c对应的数字频率为ωc=2πfc/fs\omega_c=2\pi f_c/f_s,计算出h(n)=sin(nωc)nπh(n)=\frac{\sin(n\omega_c)}{n\pi}
  3. 根据窗函数表选择合适的窗函数,选择阻带衰减不低于于题目要求的最简单的窗函数即可
  4. 根据窗函数的窗内项数计算出NN,得到窗函数的表达式
  5. 将窗函数与h(n)h(n)相乘,得到有限长脉冲响应
  6. 将脉冲响应右移,使得第一个非零值在n=0n=0处(这相当于对信号进行时移,根据DTFT时移性质,这不会改变幅频响应,只会改变相频响应,但我们对后者并不关心,因此时移不会带来滤波效果的改变)

下面对上述步骤进行实例演示

例: 对以下要求设计低通FIR滤波器

要求: 通带边缘频率10kHz\rm 10kHz,阻带边缘频率22kHz\rm 22kHz,采样频率50kHz\rm 50kHz,阻带衰减75dB\rm 75dB,可选窗函数如下表

image-20250103130648654

Step1: 截止频率为

fc=10+222=16kHzf_c=\frac{10+22}{2}=16\rm kHz

Step2: 截止频率对应的数字频率

ωc=2πfcfs=1625π\omega_c=\frac{2\pi f_c}{f_s}=\frac{16}{25}\pi

于是

h(n)=sin(16nπ/25)nπh(n)=\frac{\sin(16n\pi/25)}{n\pi}

Step3: 阻带衰减为75dB\rm 75dB,选择布莱克曼窗

Step4: 窗内项数

N=5.98fsT.W.=5.985012=24.92N=5.98\frac{f_s}{T.W.}=5.98\frac{50}{12}=24.92

取不小于NN的最小奇数,有

N=25N=25

于是窗函数为

0.42+0.5cosπn12+0.08cosπn60.42+0.5\cos\frac{\pi n}{12}+0.08\cos\frac{\pi n}{6}

Step5: 有限长脉冲响应为

h(n)=h(n)w(n)=sin(16nπ/25)nπ(0.42+0.5cosπn12+0.08cosπn6)h'(n)=h(n)w(n)=\frac{\sin(16n\pi/25)}{n\pi}(0.42+0.5\cos\frac{\pi n}{12}+0.08\cos\frac{\pi n}{6})

Step6: 脉冲响应第一个非零值为n=12n=-12,于是右移1212单位得到

h(n)=h(n12)w(n12),  0n24h'(n)=h(n-12)w(n-12),\ \ 0\leq n \leq 24

于是完成了滤波器设计

带通、高通FIR滤波器

设计带通、高通FIR滤波器只需要对低通FIR滤波器进行频移,由于DTFT的频移特性:

x(n)ejω0nX(ωω0)x(n)e^{j\omega_0 n}\Leftrightarrow X(\omega-\omega_0)

又由于数字滤波器的冲激响应需要是实偶对称序列,因此这里乘的频移变换核需要是实偶序列,取出其实部,则使用

x(n)cos(ω0n)x(n)\cos(\omega_0n)

作为频移后的结果

实际计算中应在Step5中进行频移,然后再计算Step6的时移


对于带通滤波器,设其通带中心频率为f0f_0,则

ω0=2πf0fs\omega_0=\frac{2\pi f_0}{f_s}

对于高通滤波器,则为

ω0=2πfs2fs=π\omega_0=\frac{2\pi\frac{f_s}{2}}{f_s}=\pi

滤波器设计的理论分析

阻带衰减的含义

分贝不止可以用来表示声音大小,还可以表示功率的比值

事实上,声音的分贝数表示的也是声音功率和最低可感知功率的比值

其定义为,功率P1P_1P2P_2高的分贝数为

10log10P1P210\log_{10}\frac{P_1}{P_2}

即:10dB=1B\rm 10dB=1B代表功率比值为1010,于是每1dB\rm 1dB代表101/1010^{1/10}的功率比值

滤波器的阻带衰减就是此意,例如阻带衰减40dB\rm 40dB表示信号通过滤波器后,阻带部分至少要衰减到原来的10410^{-4}

h(n)的含义

这里的h(n)h(n)实际上就是理想低通滤波器的单位冲激响应。考虑理想低通滤波器:

H(ω)={1,ωωc0,ωc<ωπH(\omega)=\begin{cases}1,&|\omega|\leq \omega_c\\0,&\omega_c<|\omega|\leq \pi\end{cases}

对其应用IDTFT,得到

h(n)=12πωcωcejωndω=12πjnejωnωcωc=2jsin(ωcn)2πjn=sin(ωcn)nπh(n)=\frac{1}{2\pi}\int_{-\omega_c}^{\omega_c}e^{j\omega n}{\rm d}\omega=\frac{1}{2\pi jn}e^{j\omega n}\Bigg|_{-\omega_c}^{\omega_c}=\frac{2j\sin(\omega_c n)}{2\pi jn}=\frac{\sin(\omega_c n)}{n\pi}

窗函数

按理说,可以直接使用h(n)h(n)作为数字滤波器,但其序列有无限项,在实践中无法实现,窗函数的目的就是通过对无限项的序列加窗,使得其在尽量不改变原有滤波能力的前提下转变为有限项序列。下面以矩形窗函数为例进行分析


对矩形窗函数而言,其窗函数形式为

w(n)={1,nN120,elsew(n)=\begin{cases}1,&|n|\leq \frac{N-1}{2}\\0,&\rm else\end{cases}

其频率响应为

W(ω)=n=(N1)/2(N1)/2ejωn=ejω(N1)/2(1ejωN)1ejωW(\omega)=\sum\limits_{n=-(N-1)/2}^{(N-1)/2}e^{-j\omega n}=\frac{e^{j\omega(N-1)/2}(1-e^{-j\omega N})}{1-e^{-j\omega}}

其幅频响应为

W(ω)=1ejωN1ejω|W(\omega)|=\frac{|1-e^{-j\omega N}|}{|1-e^{-j\omega}|}

其中

1ejω=(1cosω)2+sin2ω=22cosω=2sinω2|1-e^{-j\omega}|=\sqrt{(1-\cos\omega)^2+\sin^2\omega}=\sqrt{2-2\cos\omega}=2|\sin\frac{\omega}{2}|

于是

W(ω)=sin(Nω/2)sin(ω/2)|W(\omega)|=|\frac{\sin(N\omega/2)}{\sin(\omega/2)}|

加窗得到的滤波器的幅频响应为

H(ω)=H(ω)W(ω)=ππH(λ)W(ωλ)dλ|H'(\omega)|=|H(\omega)\otimes W(\omega)|=\int_{-\pi}^{\pi}H(\lambda)|W(\omega-\lambda)|{\rm d}\lambda

=ωcωcsin[N(ωλ)/2]sin[(ωλ)/2]dλ=\int_{-\omega_c}^{\omega_c}\bigg|\frac{\sin[N(\omega-\lambda)/2]}{\sin[(\omega-\lambda)/2]}\bigg|{\rm d}\lambda

这相当于将窗函数平移叠加,于是得到的新的幅频相应的过渡带近似为端点位置窗函数的主瓣宽度

即:过渡带宽度ftwf_{tw}满足

2πftwfs=2πN\frac{2\pi f_{tw}}{f_s}=\frac{2\pi}{N}

于是

N=fsftwN=\frac{f_s}{f_{tw}}

实际的过渡带宽度需要根据通带、阻带的定义对上述积分进行数值计算,比较复杂,但上面的估计结果与表格值

N=0.91fsftwN=0.91\frac{f_s}{f_{tw}}

比较接近

阻带衰减可以近似为窗函数幅频响应主瓣高度和最大旁瓣高度的比值,约为

W(0)W(3π/N)=N2N/3π=3π2=4.712\frac{|W(0)|}{|W(3\pi/N)|}=\frac{N}{2N/3\pi}=\frac{3\pi}{2}=4.712

功率为幅度的平方,因此阻带衰减估算为

4.712213.46dB4.712^2\approx \rm 13.46dB

这与表中的数据有较大出入,但这是因为我们在计算过程中使用了较大幅度的近似