lpp.m

很多matlab的源代码

function [nn,dd] = lpp(ty,n,a)
% LPP Normalized Lowpass Prototype Analog Filter.
%
%	[NN,DD] = LPP(TY,N,A) NORMALIZED ANALOG Lowpass Prototype Filter 
%	TY = 'bw', 'c1', 'c2', or 'el' 
%	N = filter order (Must be <=10 for elliptic)
%	A = Ap = passband attenuation for 'bw' and 'c1'
%	A = As = stopband attenuation for 'c2'
%	A = [Ap As] for 'el'
%	[NN,DD] = the numerator and denominator of the LPP which is: 
%	Normalized with respect to Unit Passband for 'bw', 'c1', and 'el'.
%       NOTE: For 'c2', the LPP is normalized wrto Unit Stopband (if A=As) 
%	or w. r to the  unit passband if A is specified as [Ap As].
%
%     LPP (with no input arguments) invokes the following example:
%
%     % Find the LPP for a 4th order BW filter with passband ripple = 2dB 
%     >>[N,D] = lpp('bw',4,2)


% 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 lpp,disp('Strike a key to see results of the example')
pause,[N,D]=lpp('bw',4,2),return,end


mv=3;if exist('version')==5,mv=4;end
%DESIGN LPP
if n==1,if ty=='el',ty='c2';end,end  %Special case
if ty=='el',if n>10, error('Elliptic requires order <= 10');return,end,end
ty1=0;sc=1;ap=a(1);if length(a)>1,ty1=1;as=a(2);end
esq=(10)^(.1*ap)-1;if ty1==1,esq2=1/((10)^(.1*as)-1);end

if ty=='bw',sc=(1/esq)^(.5/n);end
if ty=='c2',if ty1==1,
sc=cosh(acosh(sqrt(1/esq/esq2))/n);
else,esq2=1/((10)^(.1*ap)-1);
end,end


if ty=='el'
if ty1==0,error('Both Ap and As must be specified for elliptic'),return,end
delsq=((10)^(.1*as)-1)/esq;mm=[1-1/delsq 1/delsq];
if mv==4,m=eval('ellipke(mm)');else,m=eval('ellipk(1,mm)');end
n1=n-0.9;kr=m(1)/m(2);
ax=1.5;wsn=2.2*as/ap; %Increase if necessary
if n<8 & n~=5,wsn=wsn/2;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);
%wsn,pause
end %
end
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;p=cplxpair(p);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
%i=[1:m/2; m:-1:m/2+1];z=cplxpair(z(i(:)));
%c1r=sinh(mu1)*br;c1i=cosh(mu1)*bi;

cr=sinh(mu)*br;ci=cosh(mu)*bi;pc=cr+j*ci;p=(1)./pc;
k=real(prod(-p)/prod(-z));

else   %if ty=='el'     %ellip  %MODIFIED TO C2 IF 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=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));p=cplxpair(p);
k=real(prod(-p)/prod(-z));if r==0,k=k/sqrt(1+esq);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];