function rsf = trbw(n,d,pc)
% TRBW	Risetime, settling time and bandwidth of lowpass filters 
%
%	RSF = TRBW(N,D,P): Computes risetime, settling time and bandwidth
%	N, D are the num and den of H(s) in descending powers of s
%	P is for P% settling time (DEFAULT: P = 5)
%	NOTE: The risetime and settling time are computed for the STEP RESPONSE 
%	      The bandwidth is the half-power (3-dB) bandwidth of H(f).
%	RSF is a 3 element vector given by
%	RSF =[10-90% risetime; P% settling time; bandwidth (HERTZ)] 
%
%	Results are plotted if NO OUTPUT ARGUMENT is given
%
%	TRBW (with no input arguments) invokes the following example
%
%	% Plot tr, ts(2%) and bw for a 3rd order Butterworth filter whose 
%	% transfer function is H(s) = 1/(s^3 +2s^2 +2s +1)
%	  >>trbw(1,[1 2 2 1],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: akambard@mtu.edu
% Copyright (c) 1998


if nargin==0,help trbw,disp('Strike a key to see results of the example')
pause,trbw(1,[1 2 2 1],2),return,end


if nargin<3,pc=5;end
while n(1)==0,n(1)=[];end
while d(1)==0,d(1)=[];end
ln=length(n);ld=length(d);
if ln>=ld | n(ln)==0 | any(d)==0,error('Cannot find risetime'),return,end
gain=n(ln)/d(ld);
y=sysresp2('s',n,d,1,[1 0]);
rt=roots(d);r=abs(real(rt));
i=find(r==0);if ~isempty(i),r(i)=[];end,rm=min(r);
i=find(abs(r-rm)<=0.1);if rm>=0.99,rm=2*rm/length(i);else,rm=rm*length(i);end

%Find settling time
ys=gain*pc/100;
i=0;while i(1)<2,t=10/rm:-1/rm:0;yt=eval(y);
yt=yt.*(abs(yt)>0.0001);
i=find(abs(gain-yt)>ys);end
t1=t(i(1));dt=6/rm/600;t=t1:dt:t1+6/rm;yt=eval(y);
yt=yt.*(abs(yt)>0.0001);
i=find(abs(gain-yt)>ys);ls=length(i);iset=i(ls);ts=t(iset+1);ys=yt(iset);

%dt=ts/rm/600;
t=0:dt:ts+2/rm;yt=eval(y);yt=yt.*(abs(yt)>0.0001);
t10=find(yt>=0.1*gain);t10=t10(1);
t90=find(yt>=0.9*gain);t90=t90(1);
tr=dt*(t90-t10);

ri=0;i=find(imag(rt)~=0);if ~isempty(i),ri=rt(i);rt(i)=[];end
rf=max([min(abs(rt));min(sqrt(abs(ri)))]);
[m,p,f]=tfplot('s',n,d,[0 rf]);
mx=max(m)/sqrt(2);i=find(m>mx);
bw=[];if ~isempty(i),ll=length(i);fi=i(1);ff=i(ll);bw=f(ff+1)-f(fi);end
rsf=[tr;ts;bw];

if nargout==0
 vx=matverch;
 if vx < 4, eval('clg');else,eval('clf');end
subplot(211),plot(t,yt,'r'),grid
hold on,plot(t([t10 t90]),yt([t10 t90]),'o'),plot(ts,ys,'*'),hold off,
pcn=num2str(pc);trn=num2str(tr);tsn=num2str(ts);
title(['Step response:  tr = ' trn ' s     ts (' pcn '%) = ' tsn ' s']),
xlabel('Time [sec]')
subplot(212),plot(f,m),grid
if f(fi)==0,fp=ff+1;else,fp=[fi-1 ff+1];end
hold on,plot(f(fp),m(fp),'o'),hold off
bwn=num2str(bw);title(['Frequency response:  Bandwidth = ' bwn ' Hz']),
xlabel('frequency  [Hz]'),subplot
if exist('version') ~= 5,pause,end
end