数字信号处理第三次实验FIR

FIR数字滤波器的设计
一、实验目的
（1）掌握用窗函数法、频率采样法及优化设计法设计FIR数字滤波器的原理及方法，熟悉相应的MATLAB编程。

（2）熟悉线性相位FIR滤波器的幅频特性和相频特性。

（3）了解各种不同窗函数对滤波器性能的影响。

二、实验内容
1.窗函数归一化的幅度谱
clear;
N=45;
W1=boxcar(N);
W2=hamming(N);
W3=blackman(N);
[h1,p1]=freqz(W1,1);
subplot(3,1,1);
plot(p1/pi,20*log(abs(h1)));title('矩形窗');
xlabel('归一化频率/\pi');ylabel('幅度/dB');
[h2,p2]=freqz(W2,1);
subplot(3,1,2);
plot(p2/pi,20*log(abs(h2)));title('汉明窗');
xlabel('归一化频率/\pi');ylabel('幅度/dB');
[h3,p3]=freqz(W3,1);
subplot(3,1,3);
plot(p3/pi,20*log10(abs(h3)));title('布莱克曼窗');
xlabel('归一化频率/\pi');ylabel('幅度/dB');
00.10.20.3
0.40.50.60.70.80.91
-100
100
矩形窗
归一化频率/π幅度/d B
0.1
0.2
0.3
0.40.50.60.7
0.8
0.9
1
-200
200
汉明窗
归一化频率/π幅度/d B
00.10.20.3
0.40.50.60.70.80.91
-200
200
布莱克曼窗
归一化频率/π
幅度/d B
分析:
特点：
矩形窗:具有最窄的主瓣宽度，最大的旁瓣峰值； 汉明窗：主瓣稍宽，但有较小的旁瓣； 布莱克曼窗：主瓣宽度最宽，旁瓣最小。

2.带通滤波器
clear; N=15; M=N-1;
wc=[0.3 0.5];
a1=fir1(M,wc,'bandpass',hanning(N)); [h1,p1]=freqz(a1,1); subplot(2,2,1);
plot(p1/pi,20*log(abs(h1)));title('N=15;');grid;
axis([0,1,-300,10]);
xlabel('归一化频率/\pi');ylabel('幅度/dB'); subplot(2,2,2);plot(p1/pi,angle(h1));
xlabel('归一化频率/\pi');ylabel('相位/dB');title('N=15;');
N=45; M=N-1;
wc=[0.3 0.5];
a1=fir1(M,wc,'bandpass',hanning(N)); [h1,p1]=freqz(a1,1); subplot(2,2,3);
plot(p1/pi,20*log(abs(h1)));title('N=45;');grid;
axis([0,1,-300,10]);
xlabel('归一化频率/\pi');ylabel('幅度/dB'); subplot(2,2,4);plot(p1/pi,angle(h1));
xlabel('归一化频率/\pi');ylabel('相位/dB');title('N=45;');
0.5
1
-300
-200-100
0N=15;
归一化频率/π
幅度/d B
0.5
1
-4-202
4归一化频率/π
相位/d B
N=15;
-300
-200-1000N=45;
归一化频率/π
幅度/d B
-4-202
4归一化频率/π
相位/d B
N=45;
分析：
N =15时，3dB 带宽约为0.16pi 20dB 带宽约为0.3pi N=45时，3dB 带宽约为0.12pi 20dB 带宽约为0.24pi
N 越大，3dB 带宽、20dB 带宽越窄，阻带衰减越大，越逼近设计要求。

3. N=45;
M=N-1;
wc=[0.3 0.5];
a1=fir1(M,wc,'bandpass',hanning(N)); [h1,p1]=freqz(a1,1); subplot(3,1,1);
plot(p1/pi,20*log10(abs(h1)));title('汉宁窗');grid;
axis([0,1,-150,10]);
xlabel('归一化频率/\pi');ylabel('幅度/dB');
N=45; M=N-1;
wc=[0.3 0.5];
a1=fir1(M,wc,'bandpass',boxcar(N)); [h1,p1]=freqz(a1,1); subplot(3,1,2);
plot(p1/pi,20*log10(abs(h1)));title('矩形窗');grid; axis([0,1,-100,10]);
xlabel('归一化频率/\pi');ylabel('幅度/dB'); N=45; M=N-1;
wc=[0.3 0.5];
a1=fir1(M,wc,'bandpass',blackman(N)); [h1,p1]=freqz(a1,1); subplot(3,1,3);
plot(p1/pi,20*log10(abs(h1)));title('布莱克曼窗');grid; axis([0,1,-150,10]);
xlabel('归一化频率/\pi');ylabel('幅度/dB');
0.1
0.2
0.3
0.40.50.60.7
0.8
0.9
1
-150
-100-50
汉宁窗
归一化频率/π幅度/d B
0.1
0.2
0.3
0.40.50.60.7
0.8
0.9
1
-100
-50
0矩形窗
归一化频率/π幅度/d B
0.1
0.2
0.3
0.40.50.60.7
0.8
0.9
1
-150
-100-50
布莱克曼窗
归一化频率/π
幅度/d B
分析：
矩形窗设计的滤波器过渡带最窄，但阻带最小衰减也最差，约为21dB；
汉宁窗设计的滤波器过渡带稍宽，但有较好的阻带衰减，约为44dB；
布莱克曼窗设计的滤波器阻带衰减最好，约为74dB，但过渡带最宽。

4．clear;
beta=4;
N=40;
M=N-1;
f=[0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 ];
a=[0 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0 0];
h=fir2(M,f,a,kaiser(N,beta));
[h1,w1]=freqz(h,1);
subplot(3,2,1);
plot(w1/pi,20*log(abs(h1)));title('beta=4');grid;
axis([0,1,-200,10]);
xlabel('归一化频率/\pi');ylabel('幅度/dB');
subplot(3,2,2);plot(w1/pi,angle(h1));
xlabel('归一化频率/\pi');ylabel('相位/dB');
beta=6;
N=40;
M=N-1;
f=[0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 ];
a=[0 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0 0];
h=fir2(M,f,a,kaiser(N,beta));
[h1,w1]=freqz(h,1);
subplot(3,2,3);
plot(w1/pi,20*log(abs(h1)));title('beta=6');grid;
axis([0,1,-200,10]);
xlabel('归一化频率/\pi');ylabel('幅度/dB');
subplot(3,2,4);plot(w1/pi,angle(h1));
xlabel('归一化频率/\pi');ylabel('相位/dB');
beta=10;
N=40;
M=N-1;
f=[0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 ];
a=[0 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0 0];
h=fir2(M,f,a,kaiser(N,beta));
[h1,w1]=freqz(h,1);
subplot(3,2,5);
plot(w1/pi,20*log(abs(h1)));title('beta=10');grid;
axis([0,1,-200,10]);
xlabel('归一化频率/\pi');ylabel('幅度/dB');
subplot(3,2,6);plot(w1/pi,angle(h1)); xlabel('归一化频率/\pi');ylabel('相位/dB');
-200
-1000beta=4
归一化频率/π
幅度/d B
0.5
1
-505
归一化频率/π
相位/d B
-200
-1000beta=6
归一化频率/π
幅度/d B
-505归一化频率/π
相位/d B
-200
-1000
beta=10
归一化频率/π
幅度/d B
0.5
1
-5
05
归一化频率/π
相位/d B
分析：
Beta 越小，通带波纹越小，阻带衰减越大，越逼近设计要求。

5. clear;
N=40;
Hk=[zeros(1,4) 0.5 ones(1,3) 0.5 zeros(1,3) 0.5 ones(1,3) 0.5 zeros(1,7) 0.5 ones(1,3) 0.5 zeros(1,3) 0.5 ones(1,3) 0.5 zeros(1,3)];
k=0:N-1;
hn=real(ifft(Hk.*exp(-j*pi*((N-1)*k/N)))); [h1,w1]=freqz(hn,1); subplot(2,1,1);
plot(w1/pi,20*log10(abs(h1)));title('频率采样法');grid; axis([0 1 -50 10]);
xlabel('归一化频率/\pi');ylabel('幅度/dB'); subplot(2,1,2);plot(w1/pi,unwrap(angle(h1))); xlabel('归一化频率/\pi');ylabel('相位/dB');
-40
-200
频率采样法
归一化频率/π
幅度/d B
-4
-202
4归一化频率/π
相位/d B
分析：
过渡带取样的优化设计 ：
为了提高逼近质量，减小在通带边缘由于抽样点的陡然变化而引起的起伏
振荡，在阻带、阻带的交界处认为地安排1到几个过渡点，可以减小样点间幅度值的落差，使过渡平缓，反冲减小，阻带最小衰减增大。

6. clear;
N=40; M=N-1;
f=[0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 ]; a=[0 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0 ]; w=[1 1 1 1 1 1 1 1 1 1]; h=remez(M,f,a,w); [h1,w1]=freqz(h,1); subplot(2,1,1);
plot(w1/pi,20*log10(abs(h1)));title('雷米兹交替法');grid; %axis([0,1,-200,10]);
xlabel('归一化频率/\pi');ylabel('幅度/dB'); subplot(2,1,2);plot(w1/pi,angle(h1));
xlabel('归一化频率/\pi');ylabel('相位/dB');
00.10.20.3
0.40.50.60.70.80.91
-100
-50
50雷米兹交替法
归一化频率/π
幅度/d B
00.10.20.3
0.40.50.60.70.80.91
-4
-202
4归一化频率/π
相位/d B
分析：
用凯塞窗设计滤波器，需要额外的参数beta,beta 取值不同设计出来的滤波器频谱也很不同。

用频率采样法设计滤波器，需要设计合理的过渡点的H(k),可以减小在通带边缘由于采样点的突变而引起的起伏震荡。

用雷米兹交替算法可以精确地控制通带和阻带的边界，但是会出现波纹。

7.高通滤波器 clear; N=40; M=N-1; f=[500 800]; a=[0 1]; o=1; At=40;
d1=1-10^(-o/20); d2=10^(-At/20); dev=[d1 d2]; fs=5000;
[M,fpts,mag,wt]=remezord(f,a,dev,fs); h=remez(M,fpts,mag,wt);
[h1,w1]=freqz(h,1); %subplot(2,1,1);
plot(w1/pi,20*log10(abs(h1)));title('雷米兹交替法');grid; %axis([0,1,-200,10]);
xlabel('归一化频率/\pi');ylabel('幅度/dB'); %subplot(2,1,2);plot(w1/pi,angle(h1));
%xlabel('归一化频率/\pi');ylabel('相位/dB');
00.10.20.3
0.40.50.60.70.80.91
-60
-50
-40
-30
-20
-10
10
雷米兹交替法
归一化频率/π
幅度/d B
三、 实验思考
1. 定性地说明用本实验程序设计的FIR 滤波器的3dB 截止频率在什么位置？它等于理想频率响应H d （e jw ）的截止频率吗？
答：3dB 截止频率在最高幅度下降一半是所对应的频率。

因为设计的时候总会有误差，所以并不等于理想Hd(ejw)的截止频率。

2.如果没有给定h(n)的长度N ，而是给定了通带边缘截止频率 c ω和阻带临界频率
p ω，以及相应的衰减，能根据这些条件用窗函数法设计线性相位FIR 低通滤波器
吗？
答：可以，可以用雷米兹交替算法来设计。

四、实验结论
总结三种不同设计方法的特点，归纳设计中的主要公式。

窗口设计法
优点：
⏹ 简单、方便、实用 缺点：
⏹ 若H d （e jw ） 不能用简单的涵数表达，则很难求出h d (n) ⏹ 采用加窗法很难精确控制滤波器的通、阻带边界频率。

主要公式
∑⎰-=--
=→=→=
→1
)()()()()()(21
)()(N n n
j j d n
j j d d j d e n h e H n h n w n h d e
e H n h e H ωω
ωπ
π
ωω
ωπ
频率采样法
优点：
⏹ 可以在频域直接诶设计，并且适合于优化设计。

⏹ 特别适用于设计窄带选频滤波器。

只是只有少数几个非零值的H(k)，因而
设计计算量小。

缺点：
⏹ 抽样频率只能等于
N π2的整数倍或等于N π2的整数倍加上N
π，因而不能确保截止频率的自由取值。

要想实现自由地选择截止频率，必须增加采样点
数N ，但这又使计算量加大。

主要公式
()()[]k H IDFT n h =, 1,,1,0-=N n Λ
()()∑-=-----=
1
01
11N k k N N
z
W k H N
z z H ()()
k N
j d d e H k H πωω2=
=, 1,,1,0-=N k Λ
N 的估算公式为
()12+∆=
m N ω
π
最优化设计
数字滤波器的两种逼近设计方法，即窗口法（时域逼近法）和频率采样法（频域逼近法），
用这两种方法设计出的滤波器的频率特性都是在不同意义上对给定理想频率特性H d (e j ω
)的逼近。

说到逼近，就有一个逼近得好坏的问题，对“好”“坏”的恒量标准不同，也会得出不同的结论，窗口法和频率采样法都是先给出逼近方法，所需变量，然后再讨论其逼近特性，如果反过来要求在某种准则下设计滤波器各参数，以获取最优的结果，这就引出了最优化设计的概念，最优化设计就是将所有采样值皆为变量，以获得最优结果。

最优化设计一般需要大量的计算，所以一般需要依靠计算机进行辅助设计。

主要公式
数字信号处理第三次实验
11。