afd.m

很多matlab的源代码

function [nal,dal,nn,dd] = afd(ty,ty1,a,pb,sb,w3d,pl)
%AFD	Classical analog filter design
%
%	[N,D,NP.DP] = afd(ty,ty1,a,fp,fs,f3,pl) ANALOG FILTER DESIGN
%	ty is  'bw','c1','c2',or'el'        ty1 is  'lp','bp','bs','hp',
%	a=[Ap,As] is the array of passband and stopband attenuation in dB.
%	fp,fs,f3= passband, stopband,3dB edge(s) IN HERTZ. f3 is OPTIONAL
%	NOTE:if f3 has only one element, it is taken as center freq for bp & bs.
%	NOTE: f3 is ignored for elliptic filter design
%	pl=STRING argument for plotting 'p' for plot, 'o' for omit (DEF: 'o')
%
%	The design EXACTLY meets specs in the following order of priority:
%	at f0 (if given) or at f3 (if given) or at the passband edge(s).
%	N,D=[num,den] of actual TF descending order, 
%	NP,DP=[num,den] of Lowpass Prototype
%	NOTE: The LPP is normalized with respect to the PASSBAND EDGE
%
%       NOTE: Cannot design ELLIPTIC filters if N>10 (lp,hp) or N>20 (bp,bs) 
%
%	AFD (with no input arguments) invokes the following example
%
%	%Design a Chebyshev BP filter with passband/stopband attn = [2 30]dB
%	%passband [200 300]Hz and stopband [100 400]Hz and plot results
%	>>[nbp,dbp]=afd('c1','bp',[2 30],[200 300],[100 400],'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


if nargin==0,help afd,disp('Strike a key to see results of the example')
pause,[nbp,dbp]=afd('c1','bp',[2 30],[200 300],[100 400],'p'),return,end

%CHECK ARGUMENTS
mv=3;if exist('version')==5,mv=4;end
if nargin==6,if isstr(w3d),pl=w3d;w3d=0;else,pl='o';end,end  %%'o'=OMIT
if nargin<6,w3d=0;pl='o';end
if ty=='el',if length(w3d)>1,w3d=0;end,end  %DO NOT USE 3DB freq FOR ELLIPTIC
%CONVERT TO LPP
pbw=pb*2*pi;sbw=sb*2*pi;w3w=w3d*2*pi;[w,wsc,wbw]=conv2lpp(ty1,pbw,sbw,w3w);

%DESIGN LPP
a=sort(abs(a));ap=a(1);as=a(2);if ap==0,ap=10*log10(2);end,wp=w(1);ws=w(2);
esq=(10)^(.1*ap)-1;esq2=1/((10)^(.1*as)-1);w3=0;if length(w)>2,w3=w(3);end

if ty=='bw'
if w3~=0,wpn=wp/w3;wsn=ws/w3;w3n=1;sc=w3;np=0;
if wpn~=1,np=ceil(log(sqrt(esq))/log(wpn));end
ns=ceil(log(sqrt(1/esq2))/log(wsn));n=max(np,ns);
else
wpn=1;wsn=ws/wp;n=ceil(log(sqrt(1/esq2/esq))/log(wsn));
w3n=(1/esq)^(.5/n);sc=w3n*wp;
end

elseif ty=='c1'
n=ceil(acosh(sqrt(1/esq2/esq))/acosh(ws/wp));
w3n=cosh(acosh(1/sqrt(esq))/n);sc=wp;
if w3~=0,n=n-1;asn=0;while asn<as
n=n+1;w3n=cosh(acosh(1/sqrt(esq))/n);sc=w3/w3n;wsn=ws/sc;
asn=10*log10(1+esq*(chebyfun(n,wsn))^2);
end,end

elseif ty=='c2'
sc=ws;n=ceil(acosh(sqrt(1/esq2/esq))/acosh(ws/wp));
w3n=1/(cosh(acosh(1/sqrt(esq2))/n));
if w3~=0,n=n-1;asn=0;while asn<as
n=n+1;w3n=1/(cosh(acosh(1/sqrt(esq2))/n));sc=w3/w3n;wsn=ws/sc;
asn=10*log10(1+1/(esq2*(chebyfun(n,1/wsn))^2));
end,end
if w3==0,wpa=1/cosh(acosh(sqrt(1/esq/esq2))/n);sc=wp/wpa;end

elseif ty=='el'
sc=wp;wpn=1;wsn=ws/wp;wsq=wsn*wsn;delsq=((10)^(.1*as)-1)/esq;
mm=[1/wsq 1-1/delsq 1-1/wsq 1/delsq];
if mv==4,m=eval('ellipke(mm)');else,m=eval('ellipk(1,mm)');end
n1=m(1)*m(2)/m(3)/m(4);
n=ceil(n1);kr=m(2)/m(4);n1=n-0.9;
ax=1.5;wsn=2.2*as/ap;if n<8 & n~=5,wsn=wsn/2;end

if n>10
mess1='Order exceeds ';mess3='Cannot design elliptic filter';
if ty1=='bp'|ty1=='bs', mess2='20. ';else, mess2='10. ';end
error([mess1,mess2,mess3]),return,
end 

while n-n1>1e-6
w2=wsn;n2=n1;wsn=wsn-wsn*(n-n1)/ax;q=wsn*wsn;
if q<=1,wsn=(1-1e-3)*w2;q=wsn*wsn;end
if mv==4,ma=eval('ellipke([1/q 1-1/q])');
else,ma=eval('ellipk(1,[1/q 1-1/q])');end
n1=kr*ma(1)/ma(2);
%if wsn<=1 | n1>n,wsn=w2;n1=n2;ax=2*ax;end
if wsn<=1 | n1>n, wsn=w2;n1=n2;
%if n~=6,ax=3*ax;else,ax=2*ax;end
ax=3*ax-ax*(n==6);
end
end

else
disp('Unlnown filter type: Use bw, c1, c2 or  el'),return
end

%FILTER ROOTS
n1=fix(n/2);ak=[];wk=[];j=sqrt(-1);
if n1>0,k=(1:n1)';tk=.5*pi*(2*k-1)/n;ak=-sin(tk);wk=cos(tk);
ak=[ak;ak];wk=[wk;-wk];end,if rem(n,2)==1,ak=[ak;-1];wk=[wk;0];end
p=ak+j*wk;
br=real(p);bi=imag(p);

if ty=='bw'
z=[];k=1;

elseif ty=='c1'     %Cheb1
es=sqrt(esq);mu=asinh(1/es)/n;c1r=sinh(mu)*br;c1i=cosh(mu)*bi;p=c1r+j*c1i;
z=[];k=real(prod(-p));if rem(n,2)==0,k=k/sqrt(1+esq);end

elseif ty=='c2'     %Cheb2
mu=asinh(1/sqrt(esq2))/n;mu1=asinh(1/sqrt(esq))/n;
if (rem(n,2)),m=n-1;z=j./cos([1:2:n-2 n+2:2:2*n-1]*pi/(2*n))';
else,m=n;z=j./cos((1:2:2*n-1)*pi/(2*n))';end
cr=sinh(mu)*br;ci=cosh(mu)*bi;pc=cr+j*ci;p=(1)./pc;k=real(prod(-p)/prod(-z));

else     %ellip
if n==1,z=[];k=1;sc=wp*(1/esq)^(.5);else  %Special case, n=1
wsq=wsn*wsn;
if mv==4,m=eval('ellipke(1/wsq)');else,m=eval('ellipk(1,1/wsq)');end
n1=fix(n/2);r=rem(n,2);
if r==1,zm=ellipj(2*(1:n1)*m(1)/n,ones(1,n1)/wsq);
else,zm=ellipj((2*(1:n1)-1)*m(1)/n,ones(1,n1)/wsq);end
pm=wsn./zm;z=[j*pm -j*pm];
z=z(:); %%z=cplxpair(z(:));
cee=prod((1-pm.*pm)./(1-zm.*zm));
np=poly([pm.*pm]);nz=poly([zm.*zm]);npm=conv(np,np);nzm=conv(nz,nz);
if r==1,nzm=[nzm 0];npm=[0 npm];end
den=esq*cee*cee*nzm+npm;p=-sqrt(-roots(den));
k=real(prod(-p)/prod(-z));if r==0,k=k/sqrt(1+esq);end
end,end

%NORMALIZED LPP TRANSFER FUNCTION
dd=real(poly(p));if ty=='bw' | ty=='c1',nn=k;else,nn=k*real(poly(z));end
az=length(dd)-length(nn);if az>0,nn=[zeros(1,az) nn];end
if sc ~= 1,aa=sc.^(0:n);nn=nn.*aa;dd=dd.*aa;end %This is for c2 filters
tfn=[nn;dd];nn1=nn;dd1=dd;

%DENORMALIZATION
if ty1=='lp',aa=wsc.^(0:n);nn1=nn1.*aa;dd1=dd1.*aa;end
if ty1=='hp',aa=wsc.^(0:n);nn1=nn1(n+1:-1:1).*aa;dd1=dd1(n+1:-1:1).*aa;
bb=nn1(1);nn1=nn1/bb;dd1=dd1/bb;end

%START OF NEW
if ty1=='bp',[nn1,dd1]=polymap(nn1,dd1,[1 0 1],[0 wbw/wsc 0]);tol=1e-10;
nn1=round(nn1).*(abs(nn1-round(nn1))<=tol)+nn1.*(abs(nn1-round(nn1))>tol);
dd1=round(dd1).*(abs(dd1-round(dd1))<=tol)+dd1.*(abs(dd1-round(dd1))>tol);
if abs(dd1(1)-round(dd1(1)))<1e-8,nn1=nn1/dd1(1);dd1=dd1/dd1(1);end 
aa=wsc.^(0:2*n);nn1=nn1.*aa;dd1=dd1.*aa;end

if ty1=='bs',[nn1,dd1]=polymap(nn1,dd1,[0 wbw/wsc 0],[1 0 1]);tol=1e-10;
nn1=round(nn1).*(abs(nn1-round(nn1))<=tol)+nn1.*(abs(nn1-round(nn1))>tol);
dd1=round(dd1).*(abs(dd1-round(dd1))<=tol)+dd1.*(abs(dd1-round(dd1))>tol);
if abs(dd1(1)-round(dd1(1)))<1e-8,nn1=nn1/dd1(1);dd1=dd1/dd1(1);end 
aa=wsc.^(0:2*n);nn1=nn1.*aa;dd1=dd1.*aa;end
nal=nn1/dd1(1);dal=dd1/dd1(1);tf=[nal;dal];

if pl == 'o',return,end


%PLOTS: NORMALIZED LPP
 vx=matverch;
 if vx < 4, eval('clg');else,eval('clf');end
subplot(221),s=1.2*max([abs(p);abs(z)]);
if ty=='c1',s=max([s 1.2*cosh(mu) 1.2*sinh(mu)]);end
s=[-s s -s s];
if mv==3,axis('square'),axis(s);end
ca=real(p)';cb=imag(p)';plot([ca;-ca],[cb;cb],'x'),
if mv==4,axis('square'),axis(s);end
hold on
z1=zeros(1,length(ca));plot([z1;ca],[z1;cb],':')
if ty=='bw',ellipse(1,':');end,
if ty=='c1'
ellipse([sinh(mu),cosh(mu)],':');end,
if ty=='c2'|ty=='el',ca=real(z)';cb=imag(z)';plot([ca;-ca],[cb;cb],'o'),end
title([ty ' PZ Plot, n = ' num2str(n)]),grid,hold off
if mv==3,axis('normal');end
l=100;if ty=='el'|ty=='c2',w0=[0:.01:6];else,w0=[0:.01:2];end
ss=sqrt(-1)*w0;h=polyval(nn,ss)./polyval(dd,ss);
hm=abs(h);hp=unwrap(angle(h(1:l)));

hpd=180*hp/pi;

dl=-diff(hp)/.01;dl=[dl(1) dl];
subplot(222),plot(w0,hm/max(hm)),grid,title('Magnitude'),
subplot(223),plot(w0(1:l),hpd),grid,title('Phase Deg')
subplot(224),plot(w0(1:l),dl),grid,title('Delay Sec'),
subplot,if mv==3,xlabel(['RESULTS FOR NORMALIZED ' ty ' LP PROTOTYPE']),end
pause,hold off

%PLOTS: DENORMALIZED FILTER
ii=find(nn1~=0);jj=ii(1);kk=nn1(jj)/dd1(1);zz=roots(nn1);pp=roots(dd1);
p0=max(abs(pp));s=1.2*p0;
if ty=='c1',s=max([s 1.2*wsc*cosh(mu) 1.2*wsc*sinh(mu)]);end
s=[-s s -s s];
if mv==3,axis('square'),axis(s);end
ca=real(pp)';cb=imag(pp)';plot([ca;-ca],[cb;cb],'x')
if mv==4,axis('square'),axis(s);end
hold on,z1=zeros(1,length(ca));
if ty1=='lp'
plot([z1;ca],[z1;cb],':c3')
if ty=='bw',ellipse(p0,':');end
if ty=='c1',%if n==1,ellipse(p0,':');else,
ellipse([wsc*sinh(mu),wsc *cosh(mu)],':');end,%end
end
title(['denormalized ' ty1 ' ' ty ' poles']),grid,hold off,pause
if mv==3,axis('normal'),end

%PLOT ACTUAL FILTER SPECTRUM
w1=min([pb sb])*.3;w2=max([pb sb])*1.5;if ty1=='lp',w1=0;end
dw=(w2-w1)/500;w0=w1:dw:w2;
ss=sqrt(-1)*w0*2*pi;h=polyval(nn1,ss)./polyval(dd1,ss);

hm=abs(h);hm=hm/max(hm);hp=unwrap(angle(h));
hpd=180*hp/pi;

ass=10*ceil(as/10)+10;ss=[0 w2 -ass 0];axis(ss)
plot(w0,20*log10(hm)),axis(ss),grid,hold on,ak=-a;lp=length(pb);ls=length(sb);
l3=0;if w3d>0,l3=length(w3d);end,if (ty1=='bp' | ty1=='bs') & l3==1,l3=0;end
ap=ak(1)*ones(1,lp);as=ak(2)*ones(1,ls);if l3>0,a3=-3.01*ones(1,l3);end
plot(pb,ap,'o',sb,as,'o');if l3>0,hold on,plot(w3d,a3,'o'),hold off,end
title([ty ' ' ty1 ' Magnitude in dB vs Hz']),pause,hold off


dl=-diff(hp)/(dw*2*pi);dl=[dl(1) dl];
i=find(dl<0);if ~isempty(i),dl(i)=.5*[dl(i-1)+dl(i+1)];end
plot(w0,hpd),grid,title('Phase Deg'),pause,plot(w0,dl),grid,title('Delay Sec')
pause,plot(w0,hm),grid,hold on,ak=1 ./((10).^(a/20));
lp=length(pb);ls=length(sb);l3=0;if w3d>0,l3=length(w3d);end,
if (ty1=='bp' | ty1=='bs') & l3==1,l3=0;end
ap=ak(1)*ones(1,lp);as=ak(2)*ones(1,ls);if l3>0,a3=0.7071*ones(1,l3);end
plot(pb,ap,'o',sb,as,'o');if l3>0,hold on,plot(w3d,a3,'o'),end
title([ty ' ' ty1 ' Magnitude vs Hz'])
hold off,if mv==3,pause,end
%return,end