用MATLAB设计滤波器

用MATLAB设计滤波器1IIR滤波器的设计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)=A(z)B(z)=n-1--n-1l)za(n....a(2)zll)zb(n....b(2)zb(l)的复频响应H(jω)。[h，w]=freqz(b，a，n)可得到数字滤波器的n点的幅频响应，这n个点均匀地分布在上半单位圆(即0~π)，并将这n点频率记录在w中，相应的频率响应记录在h中。至于n值的选择没有太多的限制，只要n0的整数，但最好能选取2的幂次方，这样就可采用FFT算法进行快速计算。如果缺省，则n=512。[h，f]二freqz(b，a，n，Fs)允许指定采样终止频率Fs(以Hz为单位)，也即在0~Fs/2频率范围内选取n个频率点(记录在f中)，并计算相应的频率响应h。[h，w]=freqz(b，a，n，'whole')表示在0~2π之间均匀选取n个点计算频率响应。[h，f]=freqz(b，a，n，'whole'，Fs)则在O~Fs之间均匀选取n个点计算频率响应。h=freqz(b，a，w)计算在矢量w中指定的频率处的频率响应，但必须注意，指定的频率必须介于0和2π之间。h=freqz(b，a，f，Fs)计算在矢量f中指定的频率处的频率响应，但指定频率必须介于0和Fs之间。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(…)说明:butter函数可设计低通、带通、高通和带阻的数字和模拟滤波器，其特性为使通带内的幅度响应最大限度地平坦，这会损失截止频率处的下降斜度。在期望通带平滑的情况下，可使用butter函数，但在期望下降斜度大的场合，应使用椭圆和Chebyshev(切比雪夫)滤波器。butter函数可设计出数字域和模拟域的Butterworth滤波器。(1)数字域[b，a]=butter(n，Wn)可设计出截止频率为Wn的n阶低通Butterworth滤波器，其滤波器为H(z)=A(z)B(z)=n-1--n-1l)za(n....a(2)zll)zb(n....b(2)zb(l)截止频率是滤波器幅度响应下降至1/2处的频率，Wn∈[0，1],其中1相应于0.5fs(取样频率，即Nyquist频率)。当Wn=[WlW2](WlW2)时，butter函数产生一2n阶的数字带通滤波器，其通带为WlωW2。[b，a]=butter(n，Wn，'ftype')可设计出高通或带阻滤波器:·当ftype=high时，可设计出截止频率为Wn的高通滤波器;·当ftype=stop时，可设计出带阻滤波器，这时Wn=[WlW2]，且阻带为WlωW2利用输出变量个数的不同，可得到滤波器的另外两种表示:零极点增益和状态方程。·[z，p，k]=butter(n，Wn)或[z，p，k]=butter(n，Wn，'ftype')可得到滤波器的零极点增益表示;·[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)=A(s)B(s)=l)a(n....a(2)ssl)b(n....b(2)sb(l)s1-nn1-nn其中截止频率Wn0。模拟域的butter函数说明与数字域的完全相同，它也有六种形式。例1：设数据采样频率为900Hz，现要设计一9阶的高通Butterworth滤波器，截止频率为300Hz。%ex104.m[b,a]=butter(9,300/500,'high');freqz(b,a,128,1000)例2：设计一10阶的带通Butterworth滤波器，带通为100--200Hz，并画出滤波器的冲击响应。%ex105.mn=5;Wn=[100200]/500;[b,a]=butter(n,Wn);figure(1);freqz(b,a,128,1000)figure(2);impz(b,a,101)例3利用双线性变换设计数字Butterworth滤波器2.0pdBRP13.0SdBAS15源程序: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(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);%调用1[b,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('MagnitudeResponse');xlabel('frequencyinpiunits');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('MagnitudeindB');xlabel('frequencyinpiunits');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','YTickLabels',['50';'15';'1';'0']);subplot(2,2,2);plot(w/pi,pha/pi);title('PhaseResponse')xlabel('frequencyinpiunits');ylabel('piunits');axis([0,1,-1,1]);set(gca,'XTickMode','manual','XTick',[0,0.2,0.3,1]);set(gca,'YTickMode','manual','YTick',[-1,0,1]);gridsubplot(2,2,4);plot(w/pi,grd);title('GroupDelay');xlabel('frequencyinpiunits');ylabel('samples');axis([0,1,0,10]);set(gca,'XTickMode','manual','XTick',[0,0.2,0.3,1]);set(gca,'YTickMode','manual','YTick',[0:2:10]);grid子程序1：freqz的修正%freqz_m.mfunction[db,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,A]=sdir2cas(b,a)%DIRECT-formtoCASCADE-formconversinins-plane%-----------------------------------------%[C,B,A]=sdir2cas(b,a)%C=gaincoefficient%B=Kby3matrixofrealcoefficientscontainingbk's%A=Kby3matrixofrealcoefficientscontainingak's%b=numberatorpolynomialcoefficientsofDIRECTform%a=denominatorpolynomialcoefficientsofDIRECTform%Na=length(a)-1;Nb=length(b)-1;%computegaincoefficientCb0=b(1);b=b/b0;a0=a(1);a=a/a0;C=b0/a0;%%Debnominatorsecond-ordersections:p=cplxpair(roots(a));K=floor(Na/2);ifK*2==Na%ComputationwhenNaisevenA=zeros(K,3);forn=1:2:NaArow=p(n:1:n+1,:);Arow=poly(Arow);A(fix((n+1)/2),:)=real(Arow);endelseifNa==1%ComputationwhenNa=1A=[0real(poly(p))];else%ComputationwhenNaisoddand1A=zeros(K+1,3);forn=1:2:2*KArow=p(n:1:n+1,:);Arow=poly(Arow);A(fix(n+1)/2,:)=real(Arow);endA(K+1,:)=[0real(poly(p(Na)))];end%numeratorsecond-ordersections:z=cplxpair(roots(b));K=floor(Nb/2);ifNb==0%computationwhenNb=0B=[00poly(z)];elseifK*2==Nb%computationwhenNbisevenB=zeros(K,3);forn=1:2:NbBrow=z(n:1:n+1,:)Brow=poly(Brow);B(fix((n+1)/2),:)=real(Brow);endelseifNb==1%computationwhenNb=1B=[0real(poly(z))];else%ComputationwhenNbisoddand1B=zeros(K+1,3);forn=1:2:2*KBrow=z(n:1:n+1,:)Brow=poly(Brow);B(fix((n+1)/2),:)=real(Brow);endB(K+1,:)=[0real(poly(z(Nb)))];endchebyl功能：Chebyshev(切比雪夫)I型滤波器设计(通带等波纹)。格式:[b，b]=chebyl(n，Rp，Wn)[b，a]=chebyl(n，Rp，Wn，'ft