窗函数法FIR滤波器设计

FIR滤波器设计通常可以分为窗函数法和频率采样法两类，这里先介绍窗函数法

1. FIR 滤波器简介

考虑一个$N-1$阶FIR滤波器
$z变换$：
$$H(z)=\sum _{n=0}^{N-1}h(n)z^{-n}$$
差分方程：
$$y(n)=b_{0}x(n)+b_{1}x(n-1)+.....b_{N-1}x(n-N+1)$$
FIR滤波器只有零点，也叫全零点滤波器(all-zero filter)，也有其它的叫法，比如feedforward、non-recursive或者transversal filter等等
FIR滤波器因为只有零点，没有反馈环节，因此是稳定的，另外一个最大的特点就是容易实现严格的线性相位，当FIR的系数满足对称结构，无论是奇对称还是偶对称
$$奇对称：h(n)=h(N-1-n)$$
$$偶对称：h(n)=-h(N-1-n)$$
按照冲击响应的奇偶和对称性的奇偶，FIR可以分为以下四种类型

• Type 1:$h(n)$偶对称，系数为奇数(偶数阶)，低通、高通、带通
• Type 2:$h(n)$偶对称，系数为偶数(奇数阶)，低通、带通
• Type 3:$h(n)$奇对称，系数为奇数(偶数阶)，带通
• Type 4:$h(n)$奇对称，系数为偶数(奇数阶)，高通、带通
  　　每种Type适合什么类型的滤波器是根据$H(w)$在$0、\pi 、2\pi$等几个关键地方是否为零或对称性来判断的，在任意一本数字信号处理教材上都会有介绍
  以上几种类型主要是概念上的区分，matlab在设计出一个滤波器后也会给出这些信息，如随便用FDA工具设计一个高通后，给出的信息里可以看到奇数个系数（56，偶数阶），线性相位，Type 1型

2.窗函数法设计FIR滤波器

理想滤波器的冲击响应是非因果、无限长的，窗函数法是从时域出发，用因果、有限长的冲击响应去逼近非因果无限长的理想滤波器。
　　设有一理想低通滤波器，$H_d(e^{j\omega })、h_d(n)$,
　　$$h_d(n)=\frac{1}{2\pi }{\int }_{-\pi }^{\pi }{H}_d(e^{j\omega })e^{j\omega }d\omega$$
　　因为$H_d(e^{j\omega })$是矩形，因此$h_d(n)$是非因果无限长的，现在我们用$h(n)$来逼近，最直接的方法就是用有限长的窗来截断
　　$$h(n)=h_d(n)w(n)$$
　　这里再把$H_d(e^{j\omega })$定义的具体些，如一个截止频率为$w_c$的低通滤波器，通带增益为0dB、具有线性相位且群延迟为$\alpha$，则$H_d(e^{j\omega })$表示为以下形式
　　
$$H_d(e^{j\omega })=\{\begin{array}{c}{e}^{-j\omega \alpha },-{\omega }_c\le \omega \le {\omega }_c\\ 0,{\omega }_c<|\omega |<{\omega }_c\end{array}$$
则相应的时域$h_d(n)$为
$$h_d(n)=\frac{1}{2\pi }{\int }_{-\pi }^{\pi }{e}^{-j\omega \alpha }{e}^{j\omega }d\omega =\frac{sin[{\omega }_c(n-\alpha )]}{\pi (n-\alpha )}$$
　　这个$h_d(n)$还是非因果无限长的，最直接的方法，用一个因果的矩形窗$R_N(n)$去截断，同时，我们还要满足FIR的对称结构，因为因果的矩形窗对称中心为$(N-1)/2$，先将理想滤波器也延迟$(N-1)/2$，即
$$\alpha =(N-1)/2$$
因此，延迟后$h_d(n)$变为下式
$$h_d(n)=\frac{sin[{\omega }_c(n-\frac{(N-1)}{2})]}{\pi (n-\frac{(N-1)}{2}))}，-\infty <n<+\infty$$
再截断变为有限长，因为矩形窗非零部分等于1，逼近得到的因果有限长$h(n)$直接
$$h(n)=\frac{sin[{\omega }_c(n-\frac{(N-1)}{2})]}{\pi (n-\frac{(N-1)}{2}))}，0\le n\le N-1$$

２.1 talk is cheap,show me code

以上步骤其实就基本完成了一个窗函数法设计低通滤波器的过程，看起来挺简单，用代码验证下

close all
Nd = 16385;             % 很长很长，代表理想滤波器的无限长
fs = 16000;
fc = 4000.0;               %cut-off frequency
wc = 2*pi*fc/fs;
n = -(Nd-1)/2:1:(Nd-1)/2;
% n = 
hd = sin(n*wc)./(pi*n);
center = (Nd+1)/2;
hd(center)=wc/pi;
figure,plot(n,hd),title('理想低通')
N_FFT_hd = 16384*2;
Hd = fft(hd,N_FFT_hd);
Hd_freResp = abs(Hd(1:N_FFT_hd/2+1));
Hd_freResp = 2*pow2db(Hd_freResp/max(Hd_freResp));
figure,plot((1:N_FFT_hd/2+1)*fs/N_FFT_hd,Hd_freResp),title('理想低通频率响应'),grid on

N = 257;  %截断长度
win = hamming(N)';
cut_index = center-(257-1)/2:center+(257-1)/2;
h = hd(cut_index).*win;
figure,plot(h),title('截断h（n)')
N_FFT_h = 16384;
H = fft(h,N_FFT_h);
H_freResp = (abs(H(1:N_FFT_h/2+1)));
H_freResp = 2*pow2db(H_freResp/max(H_freResp));
figure,plot(H_freResp(1:N_FFT_h/2+1)),title('逼近低通频率响应'),grid on

可以对比看到截断后的逼近得到的滤波器阻带衰减更少，并且过渡带纹波更大

3.窗函数选择

在第二部分中实现利用窗函数设计低通滤波器的过程，截断过程是窗与理想滤波器相乘，由复卷积定理可知时域相乘等于频域卷积
　　$$H(e^{j\omega })=\frac{1}{2\pi }{\int }_{-\pi }^{\pi }{H}_d(\theta )W_R(\omega -\theta )d\theta$$
　　贴一张图说明这个卷积过程
　　
　　从这个卷积过程可以看到，窗的形状会对过渡带以及肩峰（吉布斯效应）产生影响，还是上面的代码，改变截断窗为汉明窗，结果得到的低通滤波器如下
　　
　　可以看到，最终得到的结果与使用矩形窗相比有了较大的改善，过渡带更平稳、阻带衰减更大，更接近理想滤波器的性能。这个结果也可以与直接使用FDA工具设计出来的做对比，结果是一样的。
　　在窗的选择上，只有当窗函数的能量集中在主瓣，旁瓣能量越小，即窗的形状越接近冲击函数时，$H(\omega )$才能逼近$H_d(\omega )$。
　　因此，在实际设计中，需要根据不同窗的主瓣宽带、旁瓣峰值、过渡带宽、阻带衰减等指标选择合适的窗，各个窗的一些指标在各书籍中都能查到。矩形窗的主瓣宽度最小，但是旁瓣衰减也最小，hamming窗是折中主瓣宽度与旁瓣衰减常用的一种窗。

Design of FIR Filters