function [OS,ts,tr,tp,T,U,Y,index] = StepResponse(t,u,y,tol)
% StepResponse:
%
% This function finds the percent overshoot (OS), settling time (ts), rise
% time (tr) and time to peak (tp) for a step response. This function is
% primarily intended to analyze measured data but can also be used on
% simulation results.
%
% This function is similar to the stepinfo function but does not require
% the Control Systems Toolbox. It also provides index points useful for
% plotting.
%
% useage: [OS,ts,tr,tp,T,U,Y,index] = StepResponse(t,u,y,tol)
%
% arguments (input):
% - t: the time vector
% - u: the step command vector. Must contain only 1 step command that
% starts at 0 and moves to +/-step in 1 time sample. Must be the same
% length as t.
% - y: the response vector. Must be the same length as t
% - tol: tolerance used to compare floating point numbers.
%
% arguments (output):
% OS: percent overshoot of the response.
% ts: settling time of the response. The time when the response settles
% within 2% of the command.
% tr: rise time of the response. The change in time for the response to
% move from 10% to 90% of the command.
% tp: time to peak. The time to reach the the maximum response
% T: truncated time vector, includes only the data after the step command
% is provided
% U: truncated step command vector, includes only the data after the step command
% is provided
% Y: truncated response vector, includes only the data after the step command
% is provided
% index: vector of the T and Y locations were ts, tr (upper), tr (lower)
% and tp values are found
%
% Example: This example uses tf and lsim to generate the data to be
% analyzed and therefor requires the Control Systems Toolbox.
% % classic second order system
% w = 10;
% z = 0.25;
% Hs = tf(w^2,[1,2*z*w,w^2]);
%
% % create input arguments
% t = [0:0.001:3];
% u = ones(size(t)); u(1:200) = zeros(200,1);
% y = lsim(Hs,u,t)+0.03*(rand(size(t'))-0.5);
% tol = 10^-6;
%
% % call StepResponse
% [OS,ts,tr,tp,T,U,Y,index] = StepResponse(t,u,y,tol);
%
% % display results
% disp(strcat('%OS :',num2str(OS)))
% disp(strcat('Ts :',num2str(ts)))
% disp(strcat('Tr :',num2str(tr)))
% disp(strcat('Tp :',num2str(tp)))
% % plot results
% subplot(1,2,1)
% plot(t,u,t,y)
% grid
% xlabel('Time')
% ylabel('Command and Response')
% subplot(1,2,2)
% plot(T,U,T,Y,T(index),Y(index),'*')
% grid
% xlabel('Truncated Time')
% ylabel('Truncated Command and Response')
%
% Acknowledgements:
% John D'Errico - for some many useful comments and examples on the MatLab
% File exchange
%
% Author: Clinton Cathey
% E-mail: clintoncathey@yahoo.com
% Release: 1.0
% Date: 3/17/2009
% error check input
% check number of input arguments
if (nargin < 4) || (nargin > 4)
error('StepResponse must have 4 arguments only')
end
% verify input argments are vectors
if (min(size(t)) > 1) || (min(size(u)) > 1) || (min(size(y)) > 1)
error('StepResponse arguments t, u, and y must be vectors, not matrices')
end
if (max(size(t)) < 2) || (max(size(u)) < 2) || (max(size(y)) < 2)
error('StepResponse arguments t, u, and y must be vectors, not scalars')
end
% verify input argument tol is a scalar
if (length(tol) > 1)
error('StepResponse arguments tol must be a scalars')
end
% verify input arguments are the same length
if (length(t) > length(u)+tol) || (length(t) < length(u)-tol) || (length(t) > length(y)+tol) || (length(t) < length(y)-tol)
error('StepResponse arguments t, u, and y must be vectors of the same length')
end
% find time when step command was given
t0 = 0; % initialize t0, the vector index when the step command is given
for J = 1:length(t),
% search for change in u. the 1st change above the tol is assumed
% to be the step command
if (t0 == 0) && (abs(u(J)) >= tol),
t0 = J;
end
% continue to search for changes in u. if additional changes
% greater than tol are found, generate an error
if (t0 > 0) && ((abs(u(J)) > abs(u(t0)) + tol) || (abs(u(J)) < abs(u(t0)) - tol))
error('StepResponse must have only 1 step command. The step command must start at 0units and move to +/- step size in 1 time sample.')
end
end
% truncate data series to t0 and above. zero the t vector
T = t(t0:length(t))-t(t0);
U = u(t0:length(t));
Y = y(t0:length(t));
% determine step size, ts tolerance, and tr points
StepSize = abs(U(1)); % step command size
Ts_utol = 1.02*StepSize; % upper ts tolerance
Ts_ltol = 0.98*StepSize; % lower ts tolerance
Tr_u = 0.9*StepSize; % upper tr point
Tr_l = 0.1*StepSize; % lower tr point
% initialize index
index = zeros(4,1);
% calculate OS and tp
[Ypeak,Ipeak] = max(abs(Y)); % find max response value and index
tp = T(Ipeak); % tp is time at max response
index(4) = Ipeak; % tp index
if (Ypeak >= StepSize), % if max response is greater than StepSize then calculate OS
OS = 100*(Ypeak/StepSize)-100;
else % if max response is less than StepSize then OS = 0%
OS = 0;
end
% determine ts and tr
ts = 0; % initialize ts
Tr_ucount = 0; % initialize tr counters
Tr_lcount = 0;
% search response vector
for J = length(T):-1:1,
% search for when Y is no longer within 2% of StepSize
if (ts == 0) && ((abs(Y(J)) > Ts_utol) || (abs(Y(J)) < Ts_ltol))
ts = T(J);
index(1) = J;
end
% search for when Y is ~90% of U
if (Tr_ucount == 0) && (abs(Y(J)) <= Tr_u) && (T(J) < T(Ipeak))
Tr_ucount = J;
end
% search for when Y is ~10% of U
if (Tr_lcount == 0) && (abs(Y(J)) <= Tr_l) && (T(J) < T(Ipeak))
Tr_lcount = J;
end
end
% error check ts - ts should have been set during loop and Y should
% settle out before end of data
if (ts == 0) || (ts >= max(T)-tol)
ts = NaN;
index(1) = NaN;
end
% calculate tr and store index values
tr = T(Tr_ucount)-T(Tr_lcount);
index(2) = Tr_ucount;
index(3) = Tr_lcount;
% error check tr - tr should have been set during loop and tr should be
% positive
if tr <= 0,
tr = NaN;
index(2) = NaN;
index(3) = NaN;
end