浙大本科毕设-FIR滤波器设计分析

一 、实验目的

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

（2）熟悉线性相位FIR滤波器的幅频特性和相频特性。 （3）了解各种不同窗函数对滤波器性能的影响。

二、实验内容

（1）N=45，，计算并画出矩形窗、汉明窗、布莱克曼的归一化额副频谱，并比较各自的主要特点。 程序如下： N=45;

b1=boxcar(N); b2=hamming(N); b3=blackman(N); [h1,w1]=freqz(b1,1); [h2,w2]=freqz(b2,1); [h3,w3]=freqz(b3,1);

plot(w1/pi,20*log10(abs(h1)),'r',w2/pi,20*log10(abs(h2)),'b',w3/pi,20*log10(abs(h3)),'g'); axis([0,1,-100,50]);grid;

xlabel('归一化频率');ylabel('幅度/dB'); 图形如下：

各自特点：

矩形窗：过渡带较窄，主瓣也比较窄，约为汉明窗的一半，旁瓣也幅度较大。 汉明窗：比起矩形窗和布莱克曼窗过渡带，主瓣，旁瓣幅度都居两者之间。 布莱克曼窗：主瓣较宽，旁瓣幅度小，但过渡带宽。

（2）N=15，带通滤波器的两个通带边界分别是w1=0.3,w2=0.5.用汗宁窗设计此线性相位带通滤波器，观察他的实际3dB和20dB带宽。N=45，重复这一设计，观察幅频和相频特性的变化，注意长度N变化的影响。 程序如下： N1=15;N2=45; wc1=0.3;wc2=0.5; Wc=[wc1,wc2];

b1=fir1(N1,Wc,hanning(N1+1)); b2=fir1(N2,Wc,hanning(N2+1)); [h1,w1]=freqz(b1,1); [h2,w2]=freqz(b2,1); figure(1);title('hanning');

subplot(2,2,1);plot(w1/pi,20*log(abs(h1)));

grid;xlabel('N1=15：归一化频率');ylabel('幅度/dB'); subplot(2,2,2);plot(w2/pi,20*log(abs(h2)));

grid;xlabel('N2=45：归一化频率');ylabel('幅度/dB'); subplot(2,2,3);plot(w1/pi,angle(h1));

grid;xlabel('N1=15：归一化频率');ylabel('相位'); subplot(2,2,4);plot(w2/pi,angle(h2));

grid;xlabel('N2=45：归一化频率');ylabel('相位'); 图形如下：

比较图形可知：

N增大时，主瓣变窄，因为主瓣宽为8/N，与N成反比。 旁瓣幅度不变，过渡带宽变窄。相位变化更频繁。

N=15时，3dB带宽0.2，20dB带宽约为0.3 N=45时，3dB带宽0.2，20dB带宽略大于0.2

（3）分别改用矩形窗和布莱克曼窗，设计（2）中的带通滤波器，观察并记录窗函数对滤波器幅频特性的影响，比较三种窗的特点。 程序如下： N1=15;N2=45;wc1=0.3;wc2=0.5; Wc=[wc1,wc2]; b1=fir1(N1,Wc,boxcar(N1+1)); b2=fir1(N2,Wc,boxcar(N2+1)); [h1,w1]=freqz(b1,1);[h2,w2]=freqz(b2,1); figure(2); subplot(2,1,1);plot(w1/pi,20*log(abs(h1))); grid;xabel('N1=15：归一化频率'); ylabel('幅度/dB');title('Boxcar'); subplot(2,1,2);plot(w2/pi,20*log(abs(h2))); grid;xlabel('N2=45：归一化频率'); ylabel('幅度/dB');

图形如下：

b1=fir1(N1,Wc,blackman(N1+1)); b2=fir1(N2,Wc,blackman(N2+1)); [h1,w1]=freqz(b1,1); [h2,w2]=freqz(b2,1); figure(3); subplot(2,1,1);pot(w1/pi,20*log(abs(h1))); grid;xlabel('N1=15：归一化频率'); ylabel('幅度/dB');title('Blackman'); subplot(2,1,2);plot(w2/pi,20*log(abs(h2))); grid;xlabel('N2=45：归一化频率'); ylabel('幅度/dB');

比较三种窗，各自特点如下：

汗宁窗：主瓣较宽，过渡带宽带，旁瓣幅度居中 矩形窗：主瓣最窄，过渡带最窄，旁瓣也幅度最大 布莱克曼窗：主瓣宽，过渡带最宽，旁瓣也幅度最小

（4）用凯塞窗设计一专用线性相位滤波器，N=40，|Hd（ej）|如实验四图，当β=4、6、

10时，分别设计、比较他们的幅频特性的影响，比较三种窗的特点。 程序如下：

N=40;bt1=4;bt2=6;bt3=10;n=0:1:39;af=(N-1)/2; wn1=kaiser(N,bt1);wn2=kaiser(N,bt2);wn3=kaiser(N,bt3); hd=(sin(0.4*pi*(n-af))-sin(0.2*pi*(n-af))+sin(0.8*pi*(n-af))-sin(0.6*pi*(n-af)))./(pi*(n-af)); b1=wn1'.*hd;b2=wn2'.*hd;b3=wn3'.*hd; [h1,w1]=freqz(b1);[h2,w2]=freqz(b2);[h3,w3]=freqz(b3); figure(1); subplot(2,1,1); plot(w1/pi,20*log10(abs(h1))); axis([0,1,-80,10]);grid; xlabel('归一化频率'); ylabel('幅度/dB');title('β=4'); subplot(2,1,2); plot(w1/pi,angle(h1));grid; xlabel('归一化频率'); ylabel('相位'); 图形如下：

figure(2); subplot(2,1,1); plot(w2/pi,20*log10(abs(h2))); axis([0,1,-80,10]);grid; xlabel('归一化频率'); ylabel('幅度/dB');title('β=6'); subplot(2,1,2); plot(w2/pi,angle(h2));grid; xlabel('归一化频率'); ylabel('相位'); figure(3); subplot(2,1,1); plot(w3/pi,20*log10(abs(h3))); axis([0,1,-80,10]);grid; xlabel('归一化频率'); ylabel('幅度/dB');title('β=10'); subplot(2,1,2); plot(w3/pi,angle(h3));grid; xlabel('归一化频率'); ylabel('相位');

比较三幅图形得：

随着β的增大，主瓣宽度变小，过渡带也逐步变宽，但旁瓣幅度也逐渐变小。

（5）用频率采样法设计（4）中的滤波器，过渡带分别设一个过渡点，令H（k）=0.5.比较

两种方法的结果。 程序如下： N=40;

alfa=(40-1)/2; k=0:N-1;

w=(2*pi/N)*k;

hrs=[zeros(1,2),0.5,ones(1,5),0.5,0,0.5,ones(1,5),0.5,zeros(1,5),0.5,ones(1,5),0.5,0,0.5,ones(1,5),0.5,zeros(1,3)];

k1=0:floor((N-1)/2);k2=floor((N-1)/2)+1:N-1; angH=[-alfa*(2*pi)/N*k1,alfa*(2*pi/N*(N-k2))]; H=hrs.*exp(j*angH); b=real(ifft(H)); [h,w]=freqz(b,1); figure (2);

subplot(2,1,1);plot(w/pi,20*log10(abs(h))); axis([0,1,-80,10]);grid;

xlabel('归一化频率');ylabel('幅度/dB'); subplot(2,1,2);plot(w/pi,angle(h));grid; xlabel('归一化频率');ylabel('相位'); 图形如下：

用这种方法设计，相对于(4)中主瓣变宽，过渡带也变宽，旁瓣幅度变大，

（6）用雷米兹交替算法，设计一个线性相位高通FIR数字滤波器，并比较（4）、（5）、（6）

三种方法的结果。 程序如下： N=40; M=N-1;

f=[0 0.15 0.2 0.4 0.45 0.55 0.6 0.8 0.85 1]; a=[0 0 1 1 0 0 1 1 0 0]; b=remez(M,f,a); [h,w]=freqz(b,1); figure (5);

subplot(2,1,1);plot(w/pi,20*log10(abs(h))); axis([0,1,-80,10]);grid;

xlabel('频率/Hz');ylabel('幅度/dB'); subplot(2,1,2);plot(w/pi,angle(h));grid; xlabel('频率/Hz');ylabel('相位'); 图形如下：

比较三种方法：

（4）中旁瓣幅度最小，主瓣较窄

（5）中过渡带最宽，主瓣和（6）中差不多。 （6）中过渡带最窄，旁瓣幅度最大

（7）用雷米兹交替算法，设计一个线性相位高通FIR滤波器，其指标为：通带边界频率为fc=800Hz，阻带边界fr=500Hz，通带波动δ=1dB阻带最小衰减At=40dB，采样频率fs=5000Hz。 程序如下：

δ1=1-10^(-δ/20)=0.109；δ2=10^(-At/20)=0.01；

fedge=[500,800]; maval=[0,1]; dev=[δ1,δ2]; fs=5000;

[N,fpts,mag,wt]=remezord(fedge,maval,dev,fs); b=remez(N,fpts,mag,wt); [h,w]=freqz(b,1,512);

plot(w*2500/pi,20*log10(abs(h))); grid;

xlabel('频率/Hz');ylabel('幅度/dB'); 图形如下：

三、思考题

（1）定性的说明本实验程序程序设计的FIR滤波器的3dB截止频率在什么位置？它等于理想频率响应|Hd（ej）|的截止频率吗？

答：不等于。

（2）如果没有给定的h（n）的长度N，而是给定了通带边缘截止频率wc和阻带临界频率wp，以及相应的衰减，能根据这些条件用窗函数设计线性相位FIR低通滤波器吗？

答：可以，用凯塞窗或雷米兹交替算法设计，可以估算出相应的N，画出幅频特性曲线，调整N可得到所需结果。

网上找了半天fir滤波器应用的例子竟然没有，自己写一个共同交流。

设采样频率为fs=1000Hz, 已知原始信号为x=sin（2*pi*80*t）+2*sin(2*pi*140*t),加噪声后的信号为xn=x+randn(size(t)), 下面通过设计数字FIR滤波器恢复出信号。

设计的频带 [0, 65/500]=[0, 0.13] [75/500, 85/500]=[0.15, 0.17] [95/500, 125/500]=[0.19, 0.25] [135/500, 145/500]=[0.27, 0.29] [155/500, 1]=[0.31, 1] “通/阻”标记 用“0”表示阻带 用“1”表示阻带 用“0”表示阻带 用“1”表示阻带 用“0”表示阻带

程序如下： close all ; % 采样频率

Fs = 1000 ; % Hz

% 根据采样频率确定采样时间点 t = 0:1/Fs:5 ; % 产生信号

x = sin(2*pi*80*t) + 2*sin(2*pi*140*t) ; % 对信号加噪声

xn = x + 3*randn( size(t) ) ; % 滤波器参数

n = 100 ; % 滤波器的阶数

f = [0,0.13,0.15,0.17,0.19,0.25,0.27,0.29,0.31,1] ; % 通、阻带的频点 m = [0,0,1,1,0,0,1,1,0,0] ; % 通、阻带标示 % 构造滤波器

b = firls(n, f, m) ; % 这个函数通过最小方差生成滤波器，具有线性相位。 % b为生成的滤波器抽头系数

% b = fir2(n,f,m); % 若将firls换成fir2，同样可达到目的

% 求出滤波器的频率响应

[H,W] = freqz(b,1,512,Fs) ; % 求出b/1的频率响应，设置长度为512点，采样率Fs为1000

% 画出频谱图

plot(W, abs(H) ) ; grid

xlabel('频率（采样频率=1000）') ; ylabel('幅度') ;

% 进行滤波

xo = filter(b, 1, xn) ; % 画出滤波前后对比图 nn = 500:800 ; tt = nn/Fs ; figure ;

subplot(311) ; plot(tt,x(nn))

ylabel('原始信号'), grid ; subplot(312) ; plot(tt,xn(nn)) ;

ylabel('污染信号'), grid ; subplot(313) ; plot(tt,xo(nn)) ;

ylabel('滤波信号'), grid ; xlabel('时间（秒）');