FIR滤波器设计与实现实验报告

实验二：FIR滤波器设计与实现专业班级：12电子信息工程团队成员： 顾鹏伟 陆遥 张春辉一、【实验目的】1通过实验巩固FIR滤波器的认识和理解。2熟练掌握FIR低通滤波器的窗函数设计方法。3理解FIR的具体应用。二、【实验内容】 在通信、信息处理以及信号检测等应用领域广泛使用滤波器进行去噪和信号的增强。FIR滤波器由于可实现线性相位特性以及固有的稳定特征而等到广泛应用，其典型的设计方法是窗函数设计法。设计流程如下：（1）设定指标：截止频率fc，过渡带宽度f，阻带衰减A。（2）求理想低通滤波器（LPF）的时域响应hd（n）。（3）选择窗函数w（n），确定窗长N。（4）将hd（n）右移（N-1）/2点并加窗获取线性相位FIR滤波器的单位脉冲响应h（n）。（5）求FIR的频域响应H（e ），分析是否满足指标。如不满足，转（3）重新选择，否则继续。（6）求FIR的系统函数H（z）。（7）依据差分方程由软件实现FIR滤波器或依据系统函数由硬件实现。 实验要求采用哈明窗（Hamming）设计一个FIR低通滤波器并由软件实现。哈明窗函数如下：w（n）=0.54-0.46cos（ ），0nN-1；设采样频率为fs=10kHz。实验中，窗长度N和截止频率fc应该都能调节。具体实验内容如下：（1）设计FIR低通滤波器（FIR_LPF）（书面进行）。（2）依据差分方程编程实现FIR低通滤波器。（3）输入信号x（n）=3.0sin（0.16n ）+cos（0.8n ）到fc=2000Hz，N=65的FIR_LPF，求输出信号y（n），理论计算并画出0ffs范围输入信号x（n）和输出信号y（n）的幅度谱，标出峰值频率，观察滤波器的实际输出结果，分析其正确性。（4）输入信号x（n）=1.5sin（0.2n ）-cos（0.4n ）+1.2sin（0.9n）到fc=1100Hz，N=65的FIR_LPF，求输出信号y（n），理论计算并画出0ffs范围输入信号x（n）和输出信号y（n）的幅度谱，标出峰值频率，观察滤波器的实际输出结果，分析其正确性。（5）输入信号x（n）=1.5sin（0.2n ）-cos（0.4n ）+1.2sin（0.9n ）到fc=2100Hz，N=65的FIR_LPF，求输出信号y（n），理论计算并画出0ffs范围输入信号x（n）和输出信号y（n）的幅度谱，标出峰值频率，观察滤波器的实际输出结果，分析其正确性。（6）输入信号x（n）=1.5sin（0.2n ）-cos（0.4n ）+5.0sin（0.9n ）到fc=1100Hz，N=65的FIR_LPF，求输出信号y（n），理论计算并画出0ffs范围输入信号x（n）和输出信号y（n）的幅度谱，标出峰值频率，观察滤波器的实际输出结果，分析其正确性。（7）输入信号x（n）=1.5sin（0.2n ）-cos（0.4n ）+1.2sin（0.9n）到fc=1990Hz，N=65的FIR_LPF，求输出信号y（n），理论计算并画出0ffs范围输入信号x（n）和输出信号y（n）的幅度谱，标出峰值频率，观察滤波器的实际输出结果，分析其正确性。三、【实验原理】本实验采用哈明窗函数设计因果线性相位低通滤波器FIR，所设计的FIR低通滤波器单位脉冲响应h（n）如式。h（n）=hd（n）w（n）=，0nN-1如果输入函数为x（n），则y（n）=x（n）*h(n)。这里要注意N的取值，若取值为奇数，则会导致分母为零，频谱出错，故若N为奇数，则取N+1。因为，所以有。四、【实验分析】1 输入信号x（n）=3.0sin（0.16n ）+cos（0.8n ）到fc=2000Hz，N=65的FIR_LPF时理论分析:幅度谱的峰值频率由正弦分量决定，即=800Hz，输入信号x(n)的频率分量有：f1=800Hz，f2=4000Hz，经过低通滤波器后变成输出信号y(n)的频率分量只剩下800Hz，所以峰值频率应该为fm=800Hz；实际处理：在MATLAB中输入FIR（3,1.6,1,0.8,0,0,2000,65,300）得到如下图，由实际输出结果可知经过滤波器，只有一个频率分量，峰值频率为800Hz，与理论值相同，系统满足要求。2 输入信号x（n）=1.5sin（0.2）-cos（0.4）+1.2sin（0.9）到fc=1100Hz，N=65的FIR_LPF时理论分析：幅度谱的峰值频率由正弦分量决定，即=1000Hz,输入信号x(n)的频率分量有：f1=1000Hz，f2=2000Hz，f3=4500Hz，经过低通滤波器后变成输出信号y(n)的频率分量只剩下1000Hz，所以峰值频率应该为fm=1000Hz；实际处理：在MATLAB中输入FIR(1.5，0.2，-1，0.4，1.2，0.9，1100，65，130)得到如下图，由实际输出结果可知经过滤波器，输出只有一个频率分量，峰值频率为1000Hz，与理论值相同，系统满足要求。 3 输入信号x（n）=1.5sin（0.2）-cos（0.4）+1.2sin（0.9）到fc=2100Hz，N=65的FIR_LPF时理论分析：幅度谱的峰值频率由正弦分量决定，即=1000Hz,输入信号x(n)的频率分量有：f1=1000Hz，f2=2000Hz，f3=4500Hz，经过低通滤波器后变成输出信号y(n)的频率分量只剩下1000Hz，2000Hz，所以峰值频率应该为fm=1000Hz；实际处理：在MATLAB中输入FIR（1.5，0.2，-1，0.4，1.2，0.9，2100，65，130）得到如下图，由实际输出结果可知经过滤波器，输出有两个频率分量，峰值频率为1000Hz，与理论值相同，系统满足要求。4 输入信号x（n）=1.5sin（0.2）-cos（0.4）+5.0sin（0.9）到fc=1100Hz，N=65的FIR_LPF时理论分析：幅度谱的峰值频率由正弦分量决定，即=4500Hz,输入信号x(n)的频率分量有：f1=1000Hz，f2=2000Hz，f3=4500Hz，经过低通滤波器后变成输出信号y(n)的频率分量只剩下1000Hz，所以峰值频率应该仍为fm=1000Hz；实际处理：在MATLAB中输入FIR（1.5,0.2,-1,0.4,5,0.9,1100,65,300）得到如下图，由实际输出结果可知经过滤波器，输出只有一个频率分量，峰值频率为1000Hz，与理论值相同，系统满足要求。5 输入信号x（n）=1.5sin（0.2）-cos（0.4）+1.2sin（0.9）到fc=1990Hz，N=65的FIR_LPF时理论分析：幅度谱的峰值频率由正弦分量决定，即=1000Hz,输入信号x(n)的频率分量有：f1=1000Hz，f2=2000Hz，f3=4500Hz，经过低通滤波器后变成输出信号y(n)的频率分量只剩下1000Hz，所以峰值频率应该仍为fm=1000Hz；实际处理：在MATLAB中输入FIR（1.5,0.2，-1,0.4,1.2,0.9,1990,65,300）得到如下图，通过实验结果图可知输出信号的幅度谱有两个频率分量，而理论分析应只有一个频率f=1000Hz的分量。从滤波器的幅频响应分析可知，2000Hz的频率处在过渡带中，所以该频率分量被衰减而未完全滤除。如果要去除该频率分量，滤波器的截止频率应取1500Hz左右。 五、【思考题】1、当哈明窗长度N比65小（32）或大（129）的话，实验结果如何变化？答：当N=32时，实验结果图如下图所示实验结果变化：减小N值时，输出信号的频率分量与理论值偏差更大。原因分析：当哈明窗长度减小时，窗函数的主瓣变宽，过渡带变宽，但阻带衰减基本没有变化。当N=129时，实验结果图如下图。实验结果变化：增大N值，输出信号的频率频率分量和理论值基本一样。原因分析：由以上图可以看到，当窗长变长，滤波器的过渡带变窄，阻带衰减略微变大。所以其滤波时对频率的选择更加精确，效果更好。2、当采用矩形窗的话，实验（3）（4）的结果是怎样的？答：矩形窗设计的FIR低通滤波器的单位脉冲响应为：h（n）=hd（n）w（n）= ，0nN-1实验结果如下图实验（3）实验（4）原因分析：不同的窗函数设计的滤波器，最小阻带衰减和过渡带的宽度有所不同。矩形窗与哈明窗相比，幅频特性其过渡带变窄，最小阻带衰减仅有21dB左右，为哈明窗的一半左右，所以哈明窗效果更好。矩形窗函数设计的滤波器在通带内有波动，所以上面第二图中信号经滤波后1000Hz对应的幅值变大。虽然两个滤波器都有滤波作用，可以调节滤波器的截止频率或窗函数的长度获得更好效果，但可以通过调节窗函数长度从而调节过渡带带宽，但不能增大最小阻带衰减。3、 实验（6）的实验结果说明什么？答：实验结果表明：该低通滤波器系统，能够有效地滤除高于截止频率（fc=1100Hz）的分量，及时高频分量比较大，也能有效地滤除。 六、【实验总结】通过本次实验，深刻理解窗函数对滤波器最小阻带衰减和过渡带带宽的影响，加深了对窗函数设计滤波器的理解。变换窗长N，能够影响滤波器输出结果的正确性和精确性，窗的长度可以调节过渡带的宽度。同时窗的形状，也会影响所设计的滤波器的最小阻带衰减和过渡带的宽度，如本实验中的哈明窗和矩形窗，显然哈明窗的输出结果更精确。同时在设计滤波器时，对所用的编程工具MATLAB也更为熟悉，将平时学习的变成知识用到了具体事例中，加深了理解。参考资料：1 数字信号处理理论与应用第2版.俞一彪、孙兵.东南大学出版社.20112 MATLAB教程及实训.曹弋.机械工业出版社.2008七、【附录】 源程序1 哈明窗function FIR(A,a,B,b,C,c,fc,N,S)fs=10000;wc=fc*2*pi/fs;%输入信号x，取2N长度以保证输入y有N长度为稳定输出n=0:2*N-1;x=A*sin(a*pi*n)+B*cos(b*pi*n)+C*sin(c*pi*n);%系统时域波形hn=0:(N-1);w=0.54-0.46*cos(2*pi*n/(N-1); %窗函数m=n-(N-1)/2+eps; %加一个小数以免分母为“0”hd=sin(wc*m)./(pi*m);h=hd.*w; %系统函数 %线性卷积得到y，取x从第N个点开始，作为稳定的信号输入for m=0:N-1 y(m+1)=0; %给累加初值“0” for k=0:N-1 y(m+1)=y(m+1)+h(k+1)*x(m-k+1+N); %线性卷积，取稳定输出部分的y endendx=x(n+1+N); %x取第N+1个到2N个信号作为稳定的输入信号%DFTs=0:S-1;f=fs*s/S; %横坐标点X=abs(fft(x,S)/(N/2);H=abs(fft(h,S);%Y=abs(fft(y,S)/(N/2); % y截短后做DFT，与实际输出信号DFT有差异Y=X.*H;maxim=0;maxi=0;maxk=0;max1=0;for k=0:S-1 if Y(k+1)maxim maxim=Y(k+1); maxk=k; endendfor l=0:S-1 if X(l+1)maxi maxi=X(l+1); maxl=l; endendfa=maxl/S*fs;fb=maxk/S*fs;Ma=max(X);Mb=max(Y);figure(1);subplot(3,1,1);plot(x,.-)title(x(n);subplot(3,1,2);plot(h,.-)title(h(n);subplot(3,1,3);plot(y,.-)title(y(n);figure(2);subplot(3,1,1);stem(f,X,.,fill);xlabel(频率/Hz);ylabel(输入信号);X_str=sprintf(nnnk=%dn频率=%.2fHzn幅值=%.2f,maxl,fa,Ma);text(fa+fs/50,Ma,X_str);subplot(3,1,2);plot(f,H,.-)xlabel(频率/Hz);ylabel(单位脉冲响应);H_str=sprintf(nfc=%dHz,fc);text(fc+fs/50,1,H_str);subplot(3,1,3);stem(f,Y,.b,fill);xlabel(频率/Hz);ylabel(输出信号);Y_str=sprintf(nnnk=%dn频率=%.2fHzn幅值=%.2f,maxk,fb,Mb);text(fb+fs/50,Mb,Y_str);figure(3)k = 0 :2*pi/S:2*pi*(1-1/S);plot(k,20*log10(abs(H),r); axis(0 3.14 -110 5); text(3.14,-110,); xlabel(w); ylabel(20log10|H(exp(jw)|); title(FIR滤波器幅频响应);2 矩形窗function FIR1(A,a,B,b,C,c,fc,N,S)fs=10000;wc=fc*2*pi/fs;%输入信号x，取2N长度以保证输入y有N长度为稳定输出n=0:2*N-1;x=A*sin(a*pi*n)+B*cos(b*pi*n)+C*sin(c*pi*n);%系统时域波形hn=0:(N-1);w=1; %窗函数m=n-(N-1)/2+eps; %加一个小数以免分母为“0”hd=sin(wc*m)./(pi*m);h=hd.*w; %系统函数 %线性卷积得到y，取x从第N个点开始，作为稳定的信号输入for m=0:N-1 y(m+1)=0; %给累加初值“0” for k=0:N-1 y(m+1)=y(m+1)+h(k+1)*x(m-k+1+N); %线性卷积，取稳定输出部分的y endendx=x(n+1+N); %x取第N+1个到2N个信号作为稳定的输入信号%DFTs=0:S-1;f=fs*s/S; %横坐标点X=abs(fft(x,S)/(N/2);H=abs(fft(h,S);%Y=abs(fft(y,S)/(N/2); % y截短后做DFT，与实际输出信号DFT有差异 Y=X.*H;maxim=0;maxi=0;maxk=0;max1=0; for k=0:S-1 if Y(k+1)maxim maxim=Y(k+1); maxk=k; endendfor l=0:S-1 if X(l+1)maxi maxi=X(l+1); maxl=l; endendfa=maxl/S*fs;fb=maxk/S*fs;Ma=max(X);Mb=max(Y); figure(1);subplot(3,1,1);plot(x,.-)title(x(n);subplot(3,1,2);plot(h,.-)title(h(n);subplot(3,1,3);plot(y,.-)title(y(n); figure(2);subplot(3,1,1);stem(f,X,.,fill);xlabel(频率/Hz);ylabel(输入信号);X_str=sprintf(nnnk=%dn频率=%.2fHzn幅值=%.2f,maxl,fa,Ma);text(fa+fs/50,Ma,X_str); subplot(3,1,2);plot(f,H,.-)xlabel(频率/Hz);ylabel(单位脉冲响应);H_str=sprintf(nfc=%dHz,fc);text(fc+fs/50,1,H_str); subplot(3,1,3);stem(f,Y,.b,fill);xlabel(频率/Hz);ylabel(输出信号);Y_str=sprintf(nnnk=%dn频率=%.2fHzn幅值=%.2f,maxk,fb,Mb);text(fb+fs/50,Mb,Y_str); figure(3)k = 0 :2*pi/S:2*pi*(1-1/S);plot(k,20*log10(abs(H),r); axis(0 3.14 -110 5); text(3.14,-110,); xlabel(w); ylabel(20log10|H(exp(jw)|); title(FIR滤波器幅频响应);