afdbess.m

很多matlab的源代码

function [hnum,hden,nlp,dlp,d] = afdbess(d0,fp,a0,pl)
% AFDBESS Design of Bessel filters.
%
%	[N,D,NP,DP,DATA] = AFDBESS(D0,P,Ap,PL) Designs LOWPASS bessel filter
%	D0 =Td is the required delay in SECONDS.
%	D0 = [Td tol] = includes the tolerence [Default: tol =0.01]
%	P = passband edge in HERTZ
%	If P is a 2 element vector, FP(1) is taken as the RISE TIME.
%	Ap = Passband attenuation in dB [Default: Ap = 3dB] 
%       PL = 'o' to omit plots, PL = 'p' to plot [Default: PL = 'p'].
%
%       N,D= [num,den] of transfer function 
%	NP,DP=[num,den] of prototype filter
%	The first row is the denominator of the actual TF HA(s)
%	The second row is the denominator of the normalized TF HN(s)
%	The numerator of HA is HA(N+1) (the last coefficient).
%	The numerator of HN is HN(N+1) (the last coefficient).
%       DATA is a vector containing [actual delay;actual attn] 
%
%        AFDBESS (with no input arguments) invokes the following example:
%
%	% Design a Bessel filter to have a constant delay of 0.01s and a 
%	% maximum attenuation of 1.5dB up to 25Hz 
%        >>[N,D,NP,DP] = afdbess(0.01,25,1.5)  


% 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 afdbess,disp('Strike a key to see results of the example')
pause,[N,D,NP,DP]=afdbess(0.01,25,1.5),return,end

if nargin<4,pl='p';end,
if nargin==3,if isstr(a0),pl=a0;a0=3;else,pl='p';end,end
if nargin<3,a0=3;end,
w=fp;
if length(w)==2,w=2.2*2*pi/w(1);end,
fp=w(1);
if length(d0)>1,tol=d0(2);d0=d0(1);else,tol=.01;end
wn=w(1)*d0(1)*2*pi;dn=1;n=floor(5*wn*wn/a0/log(10)+.5);d=0;a=2*a0;wc=wn;

dis=0;while (abs(dn-d)/dn>=tol) | a-a0>.005
n=n+1;b=ones(1,n+1);h=zeros(1,n+1);
if n>30 & dis==0,disp([' Expected Filter Order is > ' num2str(n)]),
st=input('Enter 0 to stop or 1 to continue : ');
if isempty(st),st=0;end,if st ~= 1,h=[];return,else,dis=1;end
end

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,
h=b/b(1);
if rem(n,2)==0,ne=n/2;no=.5*n-1;else ne=(n-1)/2;no=ne;end
k=0:ne;y1=2*k;x1=sum(((-1).^k).*h(y1+1).*(wn.^y1));
k=0:no;y1=2*k+1;x4=sum(((-1).^k).*h(y1+1).*(wn.^y1));a1=x1*x1+x4*x4;
if a1<=0,error('Impossible design'),disp(['Order is > ' num2str(n)]),end
a=10*log10(a1);d=1-(wn^(2*n))/b(1)/b(1)/a1;      % OR d=(x1*x2-x3*x4)/a1;
end

hnn=b(n+1:-1:1);d=d*d0;h0=d0.^(n:-1:0);
hn=h(n+1:-1:1);h=hn.*h0;
hden=h/h(1);hnum=hden(n+1);
dlp=hnn;nlp=dlp(n+1);
hh=[hden;hnn];
w0=0:wn:20*wn;s=sqrt(-1)*w0;
h1=abs(1 ./polyval(hn,s));d=[d;a];

if pl=='o',return,end
disp('STRIKE A KEY BETWEEN PLOTS')
n0=400;wn=[wn(1);(0:2*wn/n0:2*wn)'];w=wn/d0;
s=sqrt(-1)*wn;h1=1 ./polyval(hn,s);mag=abs(h1);db=20*log10(mag);
m1=mag(1);db1=db(1);w(1)=[];wn(1)=[];mag(1)=[];db(1)=[];h1(1)=[];
ph=unwrap(angle(h1));
ph1=ph*180/pi;dl=-diff(ph)/(w(2)-w(1));  %dl=-n0*diff(ph)/w(n0+1);
dl=[dl(1);dl];
subplot,hold off,plot(w/2/pi,db,fp,db1,'o'),
grid,title(['Order=' num2str(n) '     magnitude db']),pause
plot(w/2/pi,ph1),grid,title('phase in degrees'),pause
plot(w/2/pi,dl),grid,title('delay response'),pause
t=0:5*d0/200:5*d0;yt=sysresp2('s',1,h,1,[1 0]);y=eval(yt);
plot(t,y),grid,title('step response'),pause
subplot(221),plot(w/2/pi,db,fp,db1,'o'),grid,
title(['n = ' num2str(n) ' magnitude in dB'])
subplot(222),plot(w/2/pi,ph1),grid,title('unwrapped phase deg')
subplot(223),plot(t,y),grid,title('step response')
subplot(224),plot(w/2/pi,dl),grid,title('unwrapped delay'),subplot
if exist('version')~=5,pause,end