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

1、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

2、),'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.用汗宁窗设计

3、此线性相位带通滤波器，观察他的实际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：归

4、一化频率');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：归一化频率');yl

5、abel('相位');图形如下：比较图形可知：N增大时，主瓣变窄，因为主瓣宽为8/N，与N成反比。旁瓣幅度不变，过渡带宽变窄。相位变化更频繁。N=15时，3dB带宽0.2，20dB带宽约为0.3N=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

6、(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

7、,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');图形如下：比较三种窗，各自

8、特点如下：汗宁窗：主瓣较宽，过渡带宽带，旁瓣幅度居中矩形窗：主瓣最窄，过渡带最窄，旁瓣也幅度最大布莱克曼窗：主瓣宽，过渡带最宽，旁瓣也幅度最小（4）用凯塞窗设计一专用线性相位滤波器，N=40，|Hd（）|如实验四图，当=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*(

9、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

10、,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);

11、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

12、（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=f

13、reqz(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）

14、三种方法的结果。程序如下：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'

15、;);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,ma

16、val,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（）|的截止频率吗？答：不等于。（2）如果没有给定的h（n）的长度N，而是给定了通带边缘截止频率wc和阻带临界频率wp，以及相应的衰减，能根据这些条件用窗函数设计线性相位FIR低通滤波器吗？答：可以，用凯

17、塞窗或雷米兹交替算法设计，可以估算出相应的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用“0”表示阻带75/500, 85/500=0.15, 0.17用“1”表示阻带95/500, 125/500=0.19, 0.25用“0”表示阻带135/500, 145/5

18、00=0.27, 0.29用“1”表示阻带155/500, 1=0.31, 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) ) ;gridxlabel('频率（采样频率=1000）') ;