delay.m

很多matlab的源代码

function [num,den] = delay(ty,n,r)
% DELAY Rational function expression for analog filter group delay.
%
%       [NUM,DEN] = DELAY(TY,N,R) returns the rational function expression 
%       for the group delay of an LPP filter of order N in NUM and DEN.
%
%	TY can be 'bw' for Butterworth,'c1' for Chebyshev and 'be' for Bessel.
%       For TY = 'bw' or 'c1', R = passband ripple in dB [Default: 3.01 dB]
%
%       RESULT is plotted if NO OUTPUT arguments are specified.
%
%       DELAY (with no input arguments) invokes the following example:
%
%       % Plot the delay of a 1-dB Butterworth LPF with N = 3
%       >>delay('bw',3,1);


% 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 delay,disp('Strike a key to see results of the example')
pause,delay('bw',3,1);return,end

ty=ty(1:2);
if ty=='bw'
if nargin<3,r=10*log10(2);end,esq=(10)^(.1*r)-1;
e=sqrt(esq);nu3=(1/e)^(1/n);
m=n-1:-1:0;num=(1)./(sin((2*m+1)*pi/2/n).*(nu3.^(2*m)));
num=[0*m;num];num=num(:)';
num=[0 num]/esq;
den=[1 zeros(1,2*n-1) 1/esq];
elseif ty=='c1'
if nargin<3,r=10*log10(2);end,esq=(10)^(.1*r)-1;
es=sqrt(esq);nu=asinh(1/es)/n;l=2*n+1;num=zeros(1,l);
for m=0:n-1
add=chcoeff(2*m,2)*sinh((2*n-2*m-1)*nu)/sin((2*m+1)*pi/2/n)/es;
num=num+[zeros(1,l-2*m-1) add];
end
den=0.5*esq*chcoeff(2*n);den(2*n+1)=den(2*n+1)+1+esq/2;
num=num/den(1);den=den/den(1);
elseif ty=='be'
b=ones(1,n+1);k=3:2:2*n-1;b(1)=b(1)*prod(k);ii=1;for i=1:n-1,k=n-i+1:2*n-i;
b(i+1)=b(i+1)*prod(k);ii=i*ii;b(i+1)=b(i+1)/(2^(n-i))/ii;end
b=b(n+1:-1:1);e=b(1:2:n+1);o=b(2:2:n+1);le=length(e);lo=length(o);
e(2:2:le)=-e(2:2:le);o(2:2:lo)=-o(2:2:lo);e2=conv(e,e);o2=conv(o,o);
if rem(n,2)==1,e2=[e2 0];o2=[0 o2];else,o2=[0 o2 0];end
den=[zeros(1,length(e2));e2+o2];den=den(:)';den(1)=[];
num=den;num(1)=0;
else,error('Unrecognized filter type'),return,end
nn=length(num);dd=length(den);
if nn<dd,num=[zeros(1,dd-nn) num];elseif nn>dd,den=[zeros(1,nn-dd) den];
else,end,tf=[num;den];
if nargout==0,if ty=='be',w=0:n/100:n;else,w=0:.01:1;end
dl=polyval(num,w)./polyval(den,w);plot(w,dl),grid,
if ty=='be',t='Bessel';elseif ty=='bw',t='Butterworth';else,t='Chebyshev';end
if ty=='c1'|ty=='bw',rr=[' and  ripple = ' num2str(r) ' dB'];else,rr='';end
title(['Delay of ' t ' LPP of order n = ' int2str(n) rr]),
xlabel('Frequency   rad/s'),
if exist('version')~=5,pause,end
end