dfdiir.m

很多matlab的源代码

function [dtfn,dtfd,tfnn,tfnd,tfan,tfad] = dfdiir(ty,ty1,ty2,a,sf,pb,sb,w3d,pl)
% DFDIIR IIR digital filter design.
%
%	[N,D,NP,DP,NU,DU] = DFDIIR(TY,TY1,TY2,A,SF,FP,FS,F3,PL) 
%	TY = 'bw', 'c1', 'c2', or 'el'    TY1 = 'lp', 'bp', 'bs', 'hp'
%	TY2 = 'back','forw','cent'  : backward/forward/central diff. algorithm
%       'simp','tick','rect'        : Simpsons/Ticks/Rectangular
%       'trap','bilinear','tustin'  : Trapezoidal (Bilinear or Tustin),
%       'ams2', 'ams3', 'ams4',ams5 : Adams/Moulton/Schnider's order 2,3,4,5
%       'impulse','step','ramp'     : Response Invariance 
%       'mat0','mat1','mat2'        : Matched z-transform or modifications
%	A = [Ap,As] = attenuation in dB    SF = sampling frequency in HERTZ
%	FP,FS,F3 = passband, stopband and (optional) 3dB edge(s) IN HERTZ.
%	PL = 'o' for no plots, any other string to plot. (DEFAULT: PL ='o')
%	N,D    =Num,Den of digital filter H(z) IN DESCENDING ORDER
%	NP, DP =Num,Den of normalized analog LPP Hn(s)
%	NU, DU =Num,Den of UNWARPED ANALOG FILTER Hu(s) whose ORDER MAY DIFFER 
%
%	DFDIIR (with no input arguments) invokes the following example:
%
%	% Design a Chebyshev BS IIR filter if A = [2 30]dB, Fp = [10 40]Hz,
%	% Fs = [20 30]Hz SF = 160Hz using Bilinear transformation.
%	>>[ni,di] = dfdiir('c1','bs','trap',[2 30],160,[10 40],[20 30],'p')


% ADSP Toolbox: Version 2.0 
% For use with "Analog and Digital Signal Processing", 2nd Ed.
% Published by PWS Publishing Co.
%
% Ashok Ambardar, EE Dept. MTU, Houghton, MI 49931, USA
% http://www.ee.mtu/faculty/akambard.html
% e-mail: [e-mail]
% Copyright (c) 1998


%NOTE: 	For bilinear design, we convert the analog LPP using A2D transformation
%NOTE:	For all others, we convert the analog LPP to a digital LPP using the 
%	given s2z mapping, match the DC gain, and convert it using a D2D 
%	transformation. This approach avoids aliasing effects.

if nargin==0,help dfdiir,disp('Strike a key to see results of the example')
pause,[ni,di]=dfdiir('c1','bs','trap',[2 30],160,[10 40],[20 30],'p'),
return,end

%CHECK ARGUMENTS
if nargin==8,if isstr(w3d),pl=w3d;w3d=0;else,pl='o';end,end  %%'o'=OMIT
if nargin<8,w3d=0;pl='o';end
if ty=='el',if length(w3d)>1,w3d=0;end,end  %IGNORE 3bB freq FOR ELLIPTIC
ty2=ty2(1:4);


%FIND DIGITAL FREQUENCIES AND DESIGN LPP
[y1,y2,y3,y4]=conv2lpp(ty1,pb,sb,w3d);
th3=0;thp=pb*2*pi/sf;ths=sb*2*pi/sf;if w3d>0,th3=w3d*2*pi/sf;end
%[tfa,tfn]=afd(ty,ty1,a,thp*.5/pi,ths*.5/pi,th3*.5/pi,'o');
[tfan,tfad,tfnn,tfnd]=afd(ty,ty1,a,thp*.5/pi,ths*.5/pi,th3*.5/pi,'o');

%FOR BILINEAR ONLY
if ty2=='trap'|ty2=='tust'|ty2=='bili'
%PREWARP ANALOG FREQUENCIES
pw3=0;ppb=2*tan(thp/2);psb=2*tan(ths/2);if w3d>0,pw3=2*tan(th3/2);end
%DESIGN PREWARPED NORMALIZED ANALOG LPP
[tffn,tffd,tfnn,tfnd]=afd(ty,ty1,a,ppb,psb,pw3);
%APPLY DIRECT A2D TRANSFORMATION
[dtfn,dtfd]=lp2iir(ty1,'a',tfnn,tfnd,1,y4/sf);

else

%NOT BILINEAR
%APPLY S2Z MAPPING
%t=tfn;
if ty2=='impu'|ty2=='step'|ty2=='ramp',
[dtfn,dtfd]=s2zinvar(tfnn,tfnd,1,ty2,0);
elseif ty2(1:3)=='mat',tym=ty2(4);[dtfn,dtfd]=s2zmatch(tfnn,tfnd,1,tym,0);
else,[dtfn,dtfd]=s2zni(tfnn,tfnd,1,ty2,0);end
%tfplot('z',dtfn,dtfd,[0 0.5],0,1);pause  %%Try this
%USE D2D TRANSFORMATION
[dtfn,dtfd]=lp2iir(ty1,'d',dtfn,dtfd,1,y4/sf,1/2/pi);

end

%tfd=[dtfn;dtfd];
if pl=='o',return,end

%PLOT RESULTS
f=0:sf/800:sf/2;th=0:pi/400:pi;ss=sqrt(-1)*th;zz=exp(ss);
hd=polyval(dtfn,zz)./polyval(dtfd,zz);hd=abs(hd(:));
%ANALOG
ha=polyval(tfan,ss)./polyval(tfad,ss);ha=abs(ha(:));
 vx=matverch;
 if vx < 4, eval('clg');else,eval('clf');end
if ty2 == 'trap'|ty2=='tust'|ty2=='bili'
%UNWARPED
[nu,du]=s2zni(tfan,tfad,1,ty2);
hu=polyval(nu,zz)./polyval(du,zz);hu=abs(hu(:));
plot(f,hd,'-',f,ha,':',f,hu,'--')
title('Analog(:) Prewarped(-) and Unwarped(--) Spectrum vs f in Hz')
else,plot(f,hd,'-',f,ha,':')
title('Digital(-) & Analog(:) Magnitude Response vs f in Hz')
end
grid,if exist('version') ~= 5,pause,end