%
% PURPOSE: Bessel and Butterworth Filter Analysis
% AUTHOR:  Rob McEwen
% DATE:    98/11/9
%
% NOTE: The corner frequency = bandwidth (-3dB) for Butterworth only.
%
 if(exist('fc')==0) fc = 10; end %Corner Frequency 
 Tf = 2;      %Max time in sec for plot
 Ff = 100;    %Max freq in Hz for plot
 order = 4;

 [bb,aa] = besself(order,fc);
 besselsys = tf(bb,aa);
 [mag,phase,freq]= bode(besselsys,{.1,Ff});
 [zz,pp,kk] = besself(order,fc);
 disp 'BESSEL'
 disp ' f_n , damping ratio, time constant '
 [abs(pp)/2/pi -real(pp)./abs(pp) -1./real(pp)]
 Tf = max(-10./real(pp))
 [mmin,ifbw] = min(abs(mag-1/sqrt(2)));
 disp([' Half Power Point at ' num2str(freq(ifbw)) 'Hz'])
%
 figure(1);
 clf;
 step(besselsys,Tf);
 figure(2);
 clf;
 bode(besselsys,{.1,Ff});
%
% Butterworth Filter Analysis
%
 [bb,aa] = butter(order,fc,'s');
 buttersys = tf(bb,aa);
 figure(3);
 clf;
 step(buttersys,Tf);
 figure(4);
 clf;
 bode(buttersys,{.1,Ff});
 [zz,pp,kk] = butter(order,fc,'s');
 disp 'BUTTERWORTH'
 disp ' f_n , damping ratio, time constant '
 [abs(pp)/2/pi -real(pp)./abs(pp) -1./real(pp)]
%
% You need to run plotData, with steps_28, first, for the desired plot below.
%
 figure(5)
 clf;
  treg = [0:1/30:max(data(1:N,time))]';
  sztreg = size(treg);
  NN = min([N sztreg(1)]);
  [ybes tt] = lsim(besselsys,  -sgauge(1:NN), treg(1:NN));
  [ybutter tt] = lsim(buttersys,  -sgauge(1:NN), treg(1:NN));
  plot(data(1:NN,time),[data(1:NN,thrust) -sgauge(1:NN) ybes ybutter])
  ylabel('Thrust, Newtons')
  xlabel([' Seconds.  10 Sec Noise = ', num2str(Noise), ' Newtons, rms.'])
  grid on
  axis([310 320 0 500])
  title([myDate,' File: ',nameStr,' Bias: ', num2str(bias)])

