实验四参考 基于MATLAB的FIR数字滤波器设计

1、实验4 基于MATLAB的FIR数字滤波器设计一、 实验目的（1） 加深对数字滤波器的常用指标理解。（2） 学习数字滤波器的设计方法。二、 实验原理低通滤波器 低通滤波器的常用指标：1、 通带边缘频率，2、 阻带边缘频率 ，3、 通带起伏，4、 通带峰值起伏，5、 阻带起伏，最小阻带衰减。 三、 要求（1）在MATLAB中，熟悉函数fir1、kaiserord 、remezord、remez的使用；B = fir1(n,Wn,high,noscale)设计滤波器；n,Wn,beta,ftype = kaiserord(f,a,dev)估计滤波器阶数；n,fo,ao,w = remezord (

2、f,a,dev,fs)计算等波纹滤波器阶数n和加权函数w()；B=remez(n,f,a)进行等波纹滤波器的设计（2）阅读附录中的实例，学习FIR滤波器的设计方法及其在MATLAB中的实现；（3）给出FIR数字滤波器的冲激响应，绘出它们的幅度和相位频响曲线，讨论它们各自的实现形式和特点。四、 实验内容利用MATLAB编程，分别用窗函数法和等波纹滤波器法设计两种FIR数字滤波器，指标要求如下：通带边缘频率：，通带峰值起伏：。阻带边缘频率：，最小阻带衰减： 。附录：例1 用凯塞窗设计一FIR低通滤波器，通带边界频率 ，阻带边界频率 ，阻带衰减 不小于50dB。 解 首先由过渡带宽和阻带衰减 来决定

3、凯塞窗的N和 ， ， 上图给出了以上设计的频率特性，(a) 为N=30直接截取的频率特性(b)为凯塞窗设计的频率特性。凯塞窗设计对应的MATLAB程序为： 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); 或者：b = fir1(29,0.4,kaiser(30,4.55); h1,w1=freqz(b,1); plot(w1/pi,20*log10(abs(h1); axis(0,1,-80,10); grid;xl

4、abel(归一化频率/p) ;ylabel(幅度/dB) ;此外，还可以使用n,Wn,beta,ftype = kaiserord(f,a,dev)函数来估计滤波器阶数等，得到凯塞窗滤波器：fcuts = 0.3 0.5; %归一化频率omega/pimags = 1 0;devs = 0.05 10(-2.5);n,Wn,beta,ftype = kaiserord(fcuts,mags,devs); %计算出凯塞窗N，beta的值hh = fir1(n,Wn,ftype,kaiser(n+1,beta),noscale); freqz(hh); 例2 利用雷米兹交替算法设计等波纹滤波器，设

5、计一个线性相位低通FIR数字滤波器，其指标为：通带边界频率fc=800Hz，阻带边界fr=1000Hz，通带波动 阻带最小衰减At=40dB，采样频率fs=4000Hz。 解 在MATLAB中可以用remezord 和remez两个函数设计，其结果如图2，MATLAB程序如下： fedge=800 1000; mval=1 0; dev=0.0559 0.01; fs=4000; N,fpts,mag,wt=remezord(fedge,mval,dev,fs); b=remez(N,fpts,mag,wt); h,w=freqz(b,1,256); plot(w*2000/pi,20*log

6、10(abs(h); grid; xlabel(频率/Hz) ; ylabel(幅度/dB); 参考close all,clear all,clc;N=(40-7.95)/2.286/(0.15*pi);N=round(N)+1 %N取奇数omega=0.15*pi,0.3*pi,0.5*pi,0.65*pi;wc1=(omega(1)+omega(2)/2; %理想带通的截止频率wc2=(omega(3)+omega(4)/2;a=(N-1)/2; %满足线性相位的条件for n=0:100 %此处我只画出理想带通滤波器101个点的单位样值 if n=a %响应，当然还可以取其它值 hd(n

7、+1)=(wc2-wc1)/pi; else hd(n+1)=(sin(n-a)*wc2)-sin(n-a)*wc1)/pi/(n-a); %由频率特性求单位样值相应的公式 endendn=0:100;stem(n,hd);ns=40;b=0.5842*(ns-21)0.4+0.07886*(ns-21)for m=0:N-1 b1=b*sqrt(1-(1-2*m/(N-1)2); wk(m+1)=besl(b1)/besl(b); %调用besl函数，求得kaiser窗函数endm=0:N-1;figure;stem(m,wk); %画出窗函数for n=0:N-1 h(n+1)=hd(n+1)*wk(n+1); %加窗后的带通滤波器的单位样值响应endfigure;stem(m,h);k=2*pi/(pi/40)+1;H=zeros(1,k); %对数字带通滤波器频率特性初始化q=0;for w=-pi:pi/40:pi q=q+1; for n2=0:N-1 H(q)=H(q)+h(n2+1)*exp(-j*w*n2);%计算数字滤波器的幅频特性 endendfigure,hold on,grid on;H=20*log10(abs(H); %转换为对数幅频特性w=-pi:pi/