FIR滤波器窗口设计法的问题及其MATLAB实现

窗口设计法存在的问题

问题一

窗口设计法的工作是计算 $h_d(n)$ 和 $w(n)$，当 $H_d(e^{jw})$ 较为复杂时，$h_d(n)$ 不容易由反傅里叶变换求得

这时一般可用离散傅里叶逆变换代替连续傅里叶逆变换，求近似值
$$h_M(n)=\frac{1}{M}\sum _{k=0}^{M-1}{H}_d(k)e^{j2\pi kn/M}$$
$$H_d(k)={H_d\left(e^{j\omega }\right)|}_{\omega =\frac{2\pi k}{M}}$$
当 $M>>N$ 时，$h_M(n)\approx h_d(n)$

问题二

零阶贝塞尔函数无初等函数的解析表达式。可展开成无穷级数
$$I_0(x)=1+\sum _{k=1}^{\infty }{\left[\frac{(x/2{)}^k}{k!}\right]}^2$$

问题三

经验公式
$$\beta =\{\begin{array}{cc}0.1102(At-8.7),& At\ge 50dB\\ 0.5842(At-21{)}^{0.4}+0.07886(At-21),& 21dB<At<50dB\\ 0,& At\le 21dB\end{array}$$
$$N\approx \frac{At-8}{2.286\Delta \omega }$$
其中 $\Delta \omega$ 是过渡带宽，$At$ 是阻带最小衰减

FIR滤波器窗口设计MATLAB实现

在窗口函数方面，MATLAB可计算以下窗口函数

W=hanning(N);
W=hamming(N);
W=blackman(N);
W=kaiser(N, beta);

N是窗函数长度，beta是凯塞窗的参数 $\beta$；W是返回的窗口函数，长度为N。
对于凯塞窗，MATLAB还提供了一个根据滤波器设计指标计算窗口函数参数的函数文件：

[M, Wc, beta, ftype]=kaiserord(f, a, dev, fs)

f为通带和阻带边界频率，最高为$\frac{f_s}{2}$；a为相应频带的幅度值；dev为波动值；fs缺省为2。M是滤波器阶数（M=N-1）；$W_c$ 边界频率；beta即$\beta$；ftype是滤波器类型，低通为low

B = fir1(N,Wn,'type',win)

设计N阶(长度为N+1)、系数向量为B的线性相位FIR滤波器。
Wn为给定滤波器的边界频率，取值范围为0至1.0。type为滤波器的类型，缺省为低通，‘type’=‘low’(低通)，‘high’(高通)，‘bandpass’(带通)或‘stop’(带阻)。win为窗口类型，缺省为Hamming窗

B= fir2(N, f, a, win)

利用给定的频率向量f与幅度向量a，通过线性内插得理想滤波器的频率响应。可用于设计幅度响应为任意形式的FIR线性相位滤波器。win为窗口类型，缺省为Hamming窗

Ex1:
用凯塞窗设计FIR低通，通带边界频率 ${\omega }_c=0.3\pi$，阻带边界频率 ${\omega }_c=0.5\pi$，阻带衰减不小于50dB。

理论方案

首先求解 $h_d(n)$，根据指标求其边界频率
$${\omega }_c^{\prime }=\frac{{\omega }_c+{\omega }_r}{2}=\frac{0.3\pi +0.5\pi }{2}=0.4\pi$$
$$\begin{array}{rl}{h}_d(n)& =\frac{1}{2\pi }{\int }_{-{\omega }_c^{{}^{\prime }}}^{{\omega }_c^{{}^{\prime }}}{e}^{-j\omega \alpha }{e}^{j\omega n}d\omega \\ & =\{\begin{array}{cc}\frac{\sin \left[{\omega }_c^{\prime }(n-\alpha )\right]}{\pi (n-\alpha )}& n\ne \alpha \\ {\omega }_c^{\prime }/\pi & n=\alpha \end{array}\end{array}$$
$$\Delta \omega ={\omega }_r-{\omega }_c=0.2\pi ,At=50\mathrm{d}\mathrm{B}$$
$$N=\frac{50-8}{2.286\times 0.2\pi }\approx 30\beta =0.1102(50-8.7)=4.55$$

wn=kaiser(30,4.55);
nn=[0:1:29];
alfa=(30-1)/2;
hd=sin(0.4*pi*(nn-alfa))./(pi*(nn-alfa));
h=hd.*wn;
[h1,w1]=freqz(h,1);
plot(w1/pi,20*log10(abs(h1)));
axis([0,1,-80,10]);grid;
xlabel('归一化频率/\pi')
ylabel('幅度/dB')

MATLAB方案

阻带衰减 $At=50dB$，则阻带波动值小于
$$\delta =10^{\frac{At}{20}}=10^{\frac{-50}{20}}=0.0032\approx 0.003$$
代入相关参数

[M, Wc, beta, ftype]=kaiserord([0.3 0.5], [1 0], [0.01 0.003])

得M= 30，Wc = 0.4000，beta = 4.6017，ftype =low，与理论方案结果相近

用窗函数法设计一线性相位高通滤波器。该滤波器在通带和阻带具有相同的波动幅度 $\delta \le 0.01$，通带范围$0.35\pi \sim \pi$，阻带范围 $0\sim 0.2\pi$。

采用fir1的函数，凯塞窗，MATLAB程序如下：

f=[0.2 0.35]; 
a=[0 1]; 
dev=[0.01 0.01]; 
[M, Wc, beta, ftype]=kaiserord(f, a, dev); 
h=fir1(M,Wc,ftype, kaiser(M+1,beta)); 
[h1,w1]=freqz(h,1); 
figure(1);
stem([0:M], h); 
xlabel('n'); ylabel('h(n)'); 
figure(2);
plot(w1/pi,20*log10(abs(h1))); 
axis([0,1,-80,10]);grid; 
xlabel(‘归一化频率/\pi’); ylabel('幅度/dB')

用Kaiser窗设计满足下列指标的具有2个通带FIR滤波器
$$\begin{array}{ll}{w}_{\mathrm{s}1}=0.1\mathrm{p}\mathrm{i},& \\ w_{\mathrm{p}\mathrm{l}}=0.2\mathrm{p}\mathrm{i},& w_{\mathrm{p}2}=0.4\mathrm{p}\mathrm{i}\\ w_{\mathrm{s}2}=0.5\mathrm{p}\mathrm{i},& w_{\mathrm{s}3}=0.6\mathrm{p}\mathrm{i}\\ w_{\mathrm{p}3}=0.7\mathrm{p}\mathrm{i},& w_{\mathrm{p}4}=0.8\mathrm{p}\mathrm{i}\\ w_{\mathrm{s}4}=0.9\mathrm{p}\mathrm{i},& d_{\mathrm{s}}=0.01\end{array}$$

f=[0.1 0.2 0.4 0.5 0.6 0.7 0.8 0.9]; 
a=[0,1,0,1,0];
Rs=0.01; 
dev=Rs*ones(1,length(a)); 
[N,Wc,beta,ftype] = kaiserord(f,a,dev); 
h = fir1(N,Wc,ftype,kaiser(N+1,beta)); 
omega=linspace(0,pi,512); 
mag=freqz(h,1,omega); 
plot(omega/pi,20*log10(abs(mag))); 
xlabel('Normalized frequency'); 
ylabel('Gain, db');
grid; axis([0 1 -80 5]);