用MATLAB设计滤波器

1、用MATLAB设计滤波器1 IIR滤波器的设计l freqz功能：数字滤波器的频率响应。格式： h，w=freqz(b，a，n) h，f=freqz(b，a，n，Fs) h，w=freqz(b，a，n，'whole') h，f=freqz(b，a，n，'whole'，Fs) h=freqz(b，a，w) h=freqz(b，a，f，Fs) freqz(b，a)说明:freqz用于计算由矢量"和b构成的数字滤波器 H(z)= 的复频响应H(j)。 h，w=freqz(b，a，n)可得到数字滤波器的n点的幅频响应，这n个点均匀地分布在上半单位圆(即0)，并

2、将这n点频率记录在w中，相应的频率响应记录在h中。至于n值的选择没有太多的限制，只要n>0的整数，但最好能选取2的幂次方，这样就可采用FFT算法进行快速计算。如果缺省，则n=512。 h，f二freqz(b，a，n，Fs)允许指定采样终止频率Fs(以Hz为单位)，也即在0Fs/2频率范围内选取n个频率点(记录在f中)，并计算相应的频率响应h。 h，w=freqz(b，a，n，'whole')表示在02之间均匀选取n个点计算频率响应。h，f=freqz(b，a，n，'whole'，Fs)则在OFs之间均匀选取n个点计算频率响应。 h=freqz(b，a，w)

3、计算在矢量w中指定的频率处的频率响应，但必须注意，指定的频率必须介于0和2之间。 h=freqz(b，a，f，Fs)计算在矢量f中指定的频率处的频率响应，但指定频率必须介于0和Fs之间。l butter功能：Butterworth(比特沃思)模拟和数字滤波器设计。格式：有六种 b，a=butter(n，Wn) b，a=butter(n，Wn，'ftype') b，a=butter(n，Wn，'s') b，a=butter(n，Wn，'ftyPe'，'s') z，p，k=butter() A，B，C，D=butter() 说明: b

4、utter函数可设计低通、带通、高通和带阻的数字和模拟滤波器，其特性为使通带内的幅度响应最大限度地平坦，这会损失截止频率处的下降斜度。在期望通带平滑的情况下，可使用butter函数，但在期望下降斜度大的场合，应使用椭圆和Chebyshev(切比雪夫)滤波器。 butter函数可设计出数字域和模拟域的Butterworth滤波器。 (1)数字域 b，a=butter(n，Wn)可设计出截止频率为Wn的n阶低通Butterworth滤波器，其滤波器为 H(z)= =截止频率是滤波器幅度响应下降至1/处的频率，Wn0，1,其中1相应于0.5fs(取样频率，即Nyquist频率)。 当Wn=Wl W2

5、(Wl<W2)时，butter函数产生一2n阶的数字带通滤波器，其通带为Wl<<W2。 b，a=butter(n，Wn，'ftype')可设计出高通或带阻滤波器: ·当ftype=high时，可设计出截止频率为Wn的高通滤波器; ·当ftype=stop时，可设计出带阻滤波器，这时Wn=Wl W2，且阻带为Wl<<W2 利用输出变量个数的不同，可得到滤波器的另外两种表示:零极点增益和状态方程。 ·z，p，k=butter(n，Wn)或z，p，k=butter(n，Wn，'ftype')可得到滤波器的零极

6、点增益表示; ·A，B，C，D=butter(n，Wn)或A，B，C，D=butter(n，Wn，'ftype')可得到滤波器的状态空间表示。 (2)模拟域 b，a=butter(n，Wn，'s')可设计截止频率为Wn的n阶低通模拟Butterworth滤波器为 H(s)= =其中截止频率Wn>0。 模拟域的butter函数说明与数字域的完全相同，它也有六种形式。例1：设数据采样频率为900Hz，现要设计一9阶的高通Butterworth滤波器，截止频率为300Hz。%ex104.mb,a=butter(9,300/500,'high&#

7、39;);freqz(b,a,128,1000)例2：设计一10阶的带通Butterworth滤波器，带通为100-200 Hz，并画出滤波器的冲击响应。%ex105.mn=5;Wn=100 200/500;b,a=butter(n,Wn);figure(1);freqz(b,a,128,1000)figure(2);impz(b,a,101)例3利用双线性变换设计数字Butterworth滤波器 源程序:Wp=0.2*pi;Ws=0.3*pi;Rp=1;As=15;T=1;Fs=1/T;OmegaP=(2/T)*tan(Wp/2);OmegaS=(2/T)*tan(Ws/2);ep=sqrt

8、(10(Rp/10)-1);Ripple=sqrt(1/(1+ep*ep);Attn=1/10(As/20);cs,ds=adf_butt(OmegaP,OmegaS,Rp,As);b,a=bilinear(cs,ds,T);C,B,A=sdir2cas(b,a)cs,ds=s521(OmegaP,OmegaS,Rp,As);%调用1b,a=bilinear(cs,ds,T);b0,B,A=s522(b,a);%调用2%plotfigure(1);db,mag,pha,grd,w=freqz_m(b,a);subplot(2,2,1);plot(w/pi,mag);title('Mag

9、nitude Response');xlabel('frequency in pi units');ylabel('|H|');axis(0,1,0,1);set(gca,'XTickMode','manual','XTick',0,0.2,0.3,1);set(gca,'YTickMode','manual','YTick',0,Attn,Ripple,1);grid;subplot(2,2,3);plot(w/pi,db);title('Magni

10、tude in dB');xlabel('frequency in pi units');ylabel('decibels');axis(0,1,-40,5);set(gca,'XTickMode','manual','XTick',0,0.2,0.3,1);set(gca,'YTickMode','manual','YTick',-50,-15,-1,0);grid ;%set(gca,'YTickLabelMode','manual

11、','YTickLabels','50''15''1''0');subplot(2,2,2);plot(w/pi,pha/pi);title('Phase Response')xlabel('frequency in pi units');ylabel('pi units');axis(0,1,-1,1);set(gca,'XTickMode','manual','XTick',0,0.2,0.3,1);set(

12、gca,'YTickMode','manual','YTick',-1,0,1);gridsubplot(2,2,4);plot(w/pi,grd);title('Group Delay');xlabel('frequency in pi units');ylabel('samples');axis(0,1,0,10);set(gca,'XTickMode','manual','XTick',0,0.2,0.3,1);set(gca,'YTic

13、kMode','manual','YTick',0:2:10);grid子程序1： freqz的修正%freqz_m.mfunctiondb,mag,pha,grd,w=freqz_m(b,a);H,w=freqz(b,a,1000,'whole');H=(H(1:501)'w=(w(1:501)'mag=abs(H);db=20*log10(mag+eps)/max(mag);pha=angle(H);grd=grpdelay(b,a,w);子程序2变滤波器直接形式为级联形式%sdir2cas.mfunction C,B

14、,A=sdir2cas(b,a)%DIRECT-form to CASCADE-form conversin in s-plane%-%C,B,A=sdir2cas(b,a)%C=gain coefficient%B=K by 3 matrix of real coefficients containing bk's%A=K by 3 matrix of real coefficients containing ak's%b=numberator polynomial coefficients of DIRECT form%a=denominator polynomial co

15、efficients of DIRECT form%Na=length(a)-1;Nb=length(b)-1;%compute gain coefficient Cb0=b(1);b=b/b0;a0=a(1);a=a/a0;C=b0/a0;%Debnominator second-order sections:p=cplxpair(roots(a);K=floor(Na/2);if K*2=Na %Computation when Na is even A=zeros(K,3); for n=1:2:Na Arow=p(n:1:n+1, : ); Arow=poly(Arow); A(fix

16、(n+1)/2), : )=real(Arow); endelseif Na=1 %Computation when Na =1 A=0 real(poly(p);else %Computation when Na is odd and>1 A=zeros(K+1,3); for n=1:2:2*K Arow=p(n:1:n+1, : ); Arow=poly(Arow); A(fix(n+1)/2, : )=real(Arow); end A(K+1,:)=0 real(poly(p(Na);end%numerator second-order sections:z=cplxpair(

17、roots(b);K=floor(Nb/2);if Nb=0 %computation when Nb=0 B=0 0 poly(z);elseif K*2=Nb %computation when Nb is even B=zeros(K,3); for n=1:2:Nb Brow=z(n:1:n+1, : ) Brow=poly(Brow); B(fix(n+1)/2), : )=real(Brow); endelseif Nb=1 %computation when Nb=1 B=0 real(poly(z);else %Computation when Nb is odd and>

18、;1 B=zeros(K+1,3); for n=1:2:2*K Brow=z(n:1:n+1, : ) Brow=poly(Brow); B(fix(n+1)/2), : )=real(Brow); end B(K+1, : )=0 real(poly(z(Nb);endl chebyl功能：Chebyshev(切比雪夫)I型滤波器设计(通带等波纹)。格式: b，b=chebyl(n，Rp，Wn) b，a=chebyl(n，Rp，Wn，'ftype') b，a=chebyl(n，Rp，Wn，'s') b，a=chebyl(n，Rp，Wn，'ftype&

19、#39;，'s') z，p，k=chebyl() A，B，C，D=chebyl() 说明: chebyl函数可设计低通、带通、高通和带阻的数字和模拟ChebyshevI型滤波器，其通带内为等波纹，阻带内为单调。ChebyshevI型滤波器的下降斜度比型大，但其代价是在通带内波纹较大。 与butter函数类似，chebyI函数可设计出数字域和模拟域的ChebyshevI型滤波器。 (1)数字域 b，a=chebyI (n，Rp，Wn)可设计出n阶低通数字ChebyshevI型滤波器。其截止频率由Wn确定，通带内的波纹由Rp(分贝)确定，滤波器为 H(z)= 截止频率是滤波器幅度下

20、降至-Rp分贝处的频率，对chebyI函数，Wn0，1，Wn=l时相应于0.5fs。通带波纹Rp越小，可得到更宽的变换宽度。 当Wn=Wl W2时，chebyI函数可产生一2n阶的数字带通滤波器，其通带为Wl<<W2 b,a=chebyI (n，Rp，Wn，'ftype')可设计高通或带阻滤波器: ·当ftype=high时，可设计出截止频率为Wn的高通滤波器; ·当ftype=stop时，可设计出带阻滤波器，这时Wn=Wl W2，且阻带为Wl<<W2。 利用输出变量个数的不同，可得到滤波器的另外两种表示：零极点增益和状态方程。 &#

21、183;z，p，k=chebyI (n，Rp，Wn)或z，p，k=chebyI (n，Rp，Wn，'ftype')可得到滤波器的零极点增益表示； ·A，B，C，D=chebyI (n，Rp，Wn)或A，B，C，D=chebyI (n，Rp，Wn，'ftype')可得到滤波器的状态空间表示。 (2)模拟域 b，a=chebyI (n，Rp，Wn，'s')可设计出截止频率为Wn的n阶低通模拟ChebyshevI型滤波器H(s)= =其中截止频率Wn>0。 模拟域的chebyl函数说明与数字域完全相同，它也有六种形式。例4：设数据采样频

22、率为900Hz，现要设计一9阶的低通Butterworth滤波器，截止频率为300Hz。并设定波纹系数为0.5dB。其程序如下%ex106.mb,a=butter(9,300/500); %butterworthfigure(1);freqz(b,a,128,1000)b,a=cheby1(9,0.5,300/500); %chebyshevfigure(2);freqz(b,a,512,1000) 2 椭圆滤波器的设计 若信号由5 Hz、15 Hz、30 Hz三个正弦频率成分构成。设计一个椭圆滤波器，滤除5 Hz和30 Hz频率成分。% hz15.mFs=100; t=(0:99)/Fs;

23、s1=sin(2*pi*t*5);s2=sin(2*pi*t*15);s3=sin(2*pi*t*30); s=s1+s2+s3; b,a=ellip(4,0.1,40,10 20*2/Fs); H,w=freqz(b,a,512); S=fft(s,512); sf=filter(b,a,s); SF=fft(sf,512); subplot(3,2,1);plot(t,s);title('Time domain'); grid ylabel('original signal'); subplot(3,2,5);plot(t,sf); xlabel('

24、Time(s)'); ylabel('filted signal'); w1=(0:255)/256*(Fs/2); subplot(3,2,2);plot(w1,abs(S(1:256);title('Frequency domain'); grid ylabel('original spectrum'); subplot(3,2,4);plot(w*Fs/(2*pi),abs(H); ylabel('filter'); grid; subplot(3,2,6);plot(w1,abs(SF(1:256); xlabel

25、('Frequency (Hz)'); ylabel('filted spectrum'); grid3 利用kaiser窗函数法设计FIR DF% 利用kaiser窗函数,设计长度为45的带阻滤波器,使阻带衰减为60dB.% Hd=1 0<=w<pi/3,2*pi/3<w<=pi% Hd=0 pi/3<=w<=2*pi/3% aa.mAs=60;M=45; %或M=ceil(As-7.95)/(2.286*deltaw)+1n=0:1:M-1;beta=0.1102*(As-8.7) %若验证不满足指标要求，要增大beta或M，修正设计w_kai=(kaiser(M,beta)'wc1=pi/3;wc2=2*pi/3;hd=ideal_lp(wc1,M)+ideal_lp(pi,M)-ideal_lp(wc2,M);h=hd.*w_kai;db,mag,pha,grd,w=afreq(h,1);