FIR滤波器

FIR滤波器定义

FIR(Finite Impulse Response)滤波器：有限长单位冲激响应滤波器，一个$M$阶FIR滤波器的定义如下：
$$y(n)=\sum _{k=0}^{M-1}h(k)x(n-k)$$
其系统函数为:
$$H(w)=\sum _{k=0}^{M-1}{h}_k{e}^{-jw}$$
FIR滤波器的特点:

1）只存在零点，不存在极点，系统比较稳定。

2）有限长的序列容易设计为对称结构（偶对称和奇对称），对称结构的滤波器具有线性相位的特点。

3）与IIR滤波器相比，滤波器阶数比较高。

窗函数设计FIR滤波器

窗函数法

从理想滤波器的频域响应出发，利用逆傅里叶变换求出理想滤波器的时域响应，由于理想滤波器是非因果的，所以理想滤波器时域响应是无限长的，无法实现的，窗函数方法就是用有限长的窗函数对理想滤波器时域响应进行截取。

理想低通滤波器响应

理想低通滤波器的频域响应为，
$$H(w)=\{\begin{array}{ll}& 1,|w|<w_c\\ & 0,otherwise\end{array}$$

其中$w_c$为低通滤波器截止频率，进行傅里叶逆变换，

$$\begin{array}{rl}h(n)& =\frac{1}{2\pi }{\int }_{-w_c}^{w_c}{e}^{jwn}dw\\ & =\frac{1}{2jn\pi }{e}^{jwn}{|}_{-w_c}^{w_c}\\ & =\frac{1}{n\pi }sin(w_{c}n)\\ & =\frac{w_c}{\pi }sinc(\frac{w_{c}n}{\pi })\end{array}$$

其中，$n\in (-\infty ,\infty )$,$sinc(x)=\frac{sin(\pi x)}{\pi x}$。

窗函数截取

利用窗函数$w(n)$对理想低通滤波器时域响应进行截取，只取其中的主要部分。最简单的窗函数就是矩形窗，选取的区间保持原样，其余区间丢掉，但矩形窗的旁瓣比较高。不适合利用。
除了矩形窗，常见的窗函数有hanning、hamming、Blcakman等窗函数。
1)hanning窗：
$w(n)=0.5+0.5cos(2\pi n/M),0\le n\le M-1$
2)hamming窗：
$w(n)=0.54+0.46cos(2\pi n/M),0\le n\le M-1$
3)blackman窗：
$w(n) = 0.42 - 0.5cos(2\pi n/M) + 0.08cos(4\pi n/M) ,\ 0 \leq n \leq M-1 $
下图例举了常用窗函数的特点。

频率采样法

首先给出期望的频率响应$H_d(w)$，再利用fir滤波器的系统函数$H(w)$去逼近期望的频率响应，利用优化的方法求解fir滤波器的系数，存在多种不同的优化准则，这里使用的优化其实和波束优化是一样的，感兴趣的可以看我之前写的波束优化相关的程序。

FIR滤波器系统函数为:

$$\begin{array}{rl}H(w)& =\sum _{k=0}^{M-1}{h}_k{e}^{-jw}\\ & =\sum _{k=0}^{M-1}{h}_{k}cos(kw)+j\sum _{k=0}^{M-1}{h}_{k}sin(kw)\end{array}$$
可以看出系统函数是周期共轭对称函数，所以优化时只需要设计优化$0\le w\le \pi$区间就可以了。

最大最小准则

最大最小准则是指通带内与期望响应的纹波的最大值${\delta }_1$最小化，同时限制频率响应小于${\delta }_2$
频率区间$[0,\pi ]$频率采样为$N$个点，$N$要远大于滤波器长度$M$。根据通带截止频率和阻带开始频率对频率区间进行划分为通带区间$[0,w_p]$和阻带区间$[w_s,\pi ]$。
$minimize({\delta }_1)$
subject to:
$$\begin{gathered} max(abs(H(w_k)-H_d(w_k))),0\le w_k\le w_p \\ abs(H(w_k)\le {\delta }_2,0\le w_k\le w_p \end{gathered}$$

更多的准则可以参考波束优化的书籍《传感器阵列波束优化设计与应用》，原理都是一致的。

程序仿真

窗函数设计法

%% window function approximation 
fc = 1/2;                     %% w_c/pi
n = floor(6.2/(0.15 * fc));   
n = floor(n/2) * 2 + 1;       % set window len to old
w = hamming(n);
idx = 0:(n - 1);
idx = idx - (n - 1)/2 ;
idx = idx(:);
h = fc * sinc(fc * idx).* w;
figure
freqz(h)

频率采样法

fc = 1/2;                     %% w_c/pi
n = floor(6.2/(0.15 * fc));   
N = floor(n/2) * 2 + 1;       % set window len to old
n = (N - 1)/2;

w = linspace(0,pi,360);
A = [ones(360,1) 2*cos(kron(w',[1:n]))];
w_pass_end = 0.5 * pi;
w_stop_start = 0.7 * pi;
pass_dB = 0;
stop_dB = -30;

idx = find((w >=0) & (w<= w_pass_end));
A_pass = A(idx,:)
H_pass = 10^(stop_dB/20) * ones(length(idx),1);

idx = find((w >=w_stop_start) & (w<= pi));
A_stop = A(idx,:)
H_stop = 10^(stop_dB/20) * ones(length(idx),1);

cvx_begin
    variable d
    variable h(n+1,1);
    
    minimize (d)
    subject to
        max(abs(A_pass*h - H_pass)) <=d
        abs(A_stop*h) <= H_stop;
cvx_end

h = [flipud(h(2:end));h];
figure
freqz(h)

代码链接
github仓库

参考

https://ccrma.stanford.edu/~jos/sasp/FIR_Digital_Filter_Design.html
https://blog.csdn.net/suijue9389/article/details/87900091
鄢社锋, 马远良. 传感器阵列波束优化设计及应用[M]. Ke xue chu ban she, 2009.