FIR滤波器窗口设计法的问题及其MATLAB实现
窗口设计法存在的问题
问题一
窗口设计法的工作是计算 h d ( n ) h_d(n) hd(n) 和 w ( n ) w(n) w(n),当 H d ( e j w ) H_d(e^{jw}) Hd(ejw) 较为复杂时, h d ( n ) h_d(n) hd(n) 不容易由反傅里叶变换求得
这时一般可用离散傅里叶逆变换代替连续傅里叶逆变换,求近似值
h
M
(
n
)
=
1
M
∑
k
=
0
M
−
1
H
d
(
k
)
e
j
2
π
k
n
/
M
h_{M}(n)=\frac{1}{M} \sum_{k=0}^{M-1} H_{d}(k) e^{j 2 \pi k n / M}
hM(n)=M1k=0∑M−1Hd(k)ej2πkn/M
H
d
(
k
)
=
H
d
(
e
j
ω
)
∣
ω
=
2
π
k
M
H_{d}(k)=\left.H_{d}\left(e^{j \omega}\right)\right|_{\omega=\frac{2 \pi k}{M}}
Hd(k)=Hd(ejω)∣∣ω=M2πk
当
M
>
>
N
M>>N
M>>N 时,
h
M
(
n
)
≈
h
d
(
n
)
h_M(n) \approx h_d(n)
hM(n)≈hd(n)
问题二
零阶贝塞尔函数无初等函数的解析表达式。可展开成无穷级数
I
0
(
x
)
=
1
+
∑
k
=
1
∞
[
(
x
/
2
)
k
k
!
]
2
I_{0}(x)=1+\sum_{k=1}^{\infty}\left[\frac{(x / 2)^{k}}{k !}\right]^{2}
I0(x)=1+k=1∑∞[k!(x/2)k]2
问题三
经验公式
β
=
{
0.1102
(
A
t
−
8.7
)
,
A
t
≥
50
d
B
0.5842
(
A
t
−
21
)
0.4
+
0.07886
(
A
t
−
21
)
,
21
d
B
<
A
t
<
50
d
B
0
,
A
t
≤
21
d
B
\beta=\left\{\begin{array}{cc} 0.1102(A t-8.7), & A t \geq 50 d B \\ 0.5842(A t-21)^{0.4}+0.07886(A t-21), & 21 d B<A t<50 d B \\ 0, & A t \leq 21 d B \end{array}\right.
β=⎩⎨⎧0.1102(At−8.7),0.5842(At−21)0.4+0.07886(At−21),0,At≥50dB21dB<At<50dBAt≤21dB
N
≈
A
t
−
8
2.286
Δ
ω
N \approx \frac{A t-8}{2.286 \Delta \omega}
N≈2.286ΔωAt−8
其中
Δ
ω
\Delta \omega
Δω 是过渡带宽,
A
t
At
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为通带和阻带边界频率,最高为 f s 2 \frac{f_s}{2} 2fs;a为相应频带的幅度值;dev为波动值;fs缺省为2。M是滤波器阶数(M=N-1); W c W_c Wc 边界频率;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低通,通带边界频率 ω c = 0.3 π \omega_c=0.3\pi ωc=0.3π,阻带边界频率 ω c = 0.5 π \omega_c=0.5\pi ωc=0.5π,阻带衰减不小于50dB。
理论方案
首先求解
h
d
(
n
)
h_d(n)
hd(n),根据指标求其边界频率
ω
c
′
=
ω
c
+
ω
r
2
=
0.3
π
+
0.5
π
2
=
0.4
π
\omega_{c}^{\prime}=\frac{\omega_{c}+\omega_{r}}{2}=\frac{0.3 \pi+0.5 \pi}{2}=0.4 \pi
ωc′=2ωc+ωr=20.3π+0.5π=0.4π
h
d
(
n
)
=
1
2
π
∫
−
ω
c
′
ω
c
′
e
−
j
ω
α
e
j
ω
n
d
ω
=
{
sin
[
ω
c
′
(
n
−
α
)
]
π
(
n
−
α
)
n
≠
α
ω
c
′
/
π
n
=
α
\begin{aligned} h_{d}(n) &=\frac{1}{2 \pi} \int_{-\omega_{c}^{'}}^{\omega_{c}^{'}} e^{-j \omega \alpha} e^{j \omega n} d \omega \\ &=\left\{\begin{array}{cc} \frac{\sin \left[\omega_{c}^{\prime}(n-\alpha)\right]}{\pi(n-\alpha)} & n \neq \alpha \\ \omega_{c}^{\prime} / \pi & n=\alpha \end{array}\right. \end{aligned}
hd(n)=2π1∫−ωc′ωc′e−jωαejωndω={π(n−α)sin[ωc′(n−α)]ωc′/πn=αn=α
Δ
ω
=
ω
r
−
ω
c
=
0.2
π
,
A
t
=
50
d
B
\Delta \omega=\omega_{r}-\omega_{c}=0.2 \pi, \quad A t=50 \mathrm{dB}
Δω=ωr−ωc=0.2π,At=50dB
N
=
50
−
8
2.286
×
0.2
π
≈
30
β
=
0.1102
(
50
−
8.7
)
=
4.55
N=\frac{50-8}{2.286 \times 0.2 \pi} \approx 30 \quad \beta=0.1102(50-8.7)=4.55
N=2.286×0.2π50−8≈30β=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方案
阻带衰减
A
t
=
50
d
B
At=50dB
At=50dB,则阻带波动值小于
δ
=
1
0
A
t
20
=
1
0
−
50
20
=
0.0032
≈
0.003
\delta=10^{\frac{A t}{20}}=10^{\frac{-50}{20}}=0.0032 \approx 0.003
δ=1020At=1020−50=0.0032≈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,与理论方案结果相近
用窗函数法设计一线性相位高通滤波器。该滤波器在通带和阻带具有相同的波动幅度 δ ≤ 0.01 \delta \leq 0.01 δ≤0.01,通带范围 0.35 π ∼ π 0.35\pi\sim\pi 0.35π∼π,阻带范围 0 ∼ 0.2 π 0\sim0.2\pi 0∼0.2π。
采用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滤波器
w s 1 = 0.1 p i , w p l = 0.2 p i , w p 2 = 0.4 p i w s 2 = 0.5 p i , w s 3 = 0.6 p i w p 3 = 0.7 p i , w p 4 = 0.8 p i w s 4 = 0.9 p i , d s = 0.01 \begin{array}{ll} w_{\mathrm{s}1} =0.1 \mathrm{pi}, & \\ w_{\mathrm{pl}}=0.2 \mathrm{pi}, & w_{\mathrm{p} 2}=0.4 \mathrm{pi} \\ w_{\mathrm{s} 2}=0.5 \mathrm{pi}, & w_{\mathrm{s} 3} =0.6 \mathrm{pi} \\ w_{\mathrm{p} 3}=0.7 \mathrm{pi}, & w_{\mathrm{p} 4}=0.8 \mathrm{pi} \\ w_{\mathrm{s} 4}=0.9 \mathrm{pi}, & d_{\mathrm{s}}=0.01 \end{array} ws1=0.1pi,wpl=0.2pi,ws2=0.5pi,wp3=0.7pi,ws4=0.9pi,wp2=0.4piws3=0.6piwp4=0.8pids=0.01
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]);