transient_duration_plot_disabled.m

function [TD, maximum, transient_norm_square] =....
    transient_duration_plot_disabled(y,Fs,ED,Td,SW_Plot)


%--------------------------------------------------------------------------
%                   Benchmark on Adaptive Regulation:
%
%   Rejection of unknown/time-varying multiple narrow band disturbances
%
%--------------------------------------------------------------------------
%
%   This function computes the transient duration for narrow band disturbance
%   rejection (Application to a Benchmark problem)
%
%   The transient is the time necessary to the output to reach the value
%   +/-2.17sigma , where sigma is the standard deviation for the residual output
%   (once the algo has converged) computed during the 3 last seconds of the experiment.
%
%   sigma (x) = sqrt ( var(x))
%   where "x" is the vector of data corresponding to the 3 last seconds of
%   the experiment.
%
%   In addition the output after has reached for
%   the first time 2.17 sigma should remain within +/- 4 sigma.
%   The test has been developped assuming that the residual output is a gaussian
%   process.(ie 97% of values are between +/-2.17 sigma, and 99.99% of values are
%   between +/- 4 sigma).
%
%   However it is tolerated that 0.1% of the measurements be outside +/- 4
%   sigma. The percentage is denoted Dp. (to take in account that it is not
%   a pure gaussian variable.
%
%
%   Written by M. Alma and I.D. Landau (GIPSA-LAB)
%   Version 1, September 6, 2010
% 
%   Updated by Xu Chen, University of Washington
%   Version 2, 2010-09-26
%              Added residual norm calculation
%              Added a switch to disable the plots
% 
%  TD = transient_duration(y,Fs,ED,Td)
%  TD : transient duration (in seconds)
%
%  y : is vector of data
%  Fs : is the sampling frequency ( Fs = 800 Hz by default)
%  Ed : is the exepriment duration ( ED = 15 sec by default)
%       The experiment duration should be at least equal to 10 sec.
%  Td : Start disturbance time (Td = 5 sec by default)
%  threshold 1 = 2.17 * standard deviation (after convergence)
%  threshold 2 = 4 * standard deviation (after convergence)

if nargin<2, Fs = 800; ED = 15; Td = 5; SW_Plot = 'PlotOff'; end
if nargin<3, ED = 15; Td = 5; SW_Plot = 'PlotOff'; end
if nargin<4, Td = 5; SW_Plot = 'PlotOff'; end
if nargin<5, SW_Plot = 'PlotOff'; end

if strcmp(SW_Plot,'PlotOn');
    SW_PLOT = 1;
else
    SW_PLOT = 0;
end

if Td>=ED, error('Disturbance must be injected before the end of the experiment'),end

y=y-mean(y);

yres = y((ED-3)*Fs+1:ED*Fs);
perc = 99.9;

threshold1 = 2.17*sqrt(var(yres));
threshold2 = 4*sqrt(var(yres));
l=0;
for n = 1:length(y)
    g=0;
    if abs(y(n)) <= threshold1
        for k = n:length(y)
            if abs(y(k)) > threshold2
                g= g+1;
            end
        end
        if g<(length(y)-n+1)*((100-perc)/100)
            l=n;
            break
        end
    end
end
% g;

Dp = g/(length(y)-n+1) *100;
if SW_PLOT
    fprintf ('Percentage of points above +/- 4 sigma is :')
    Dp
end

if l==0, error('Algorithm does not converge'),end

maximum = max(abs(y));
TD=(l+1)/Fs - Td;
transient_norm_square = sum(y( Td*Fs+1 : (Td+TD)*Fs ).^2);

if SW_PLOT
    fprintf ('The transient maximum value is:')
    maximum
    
    fprintf ('Transient duration is :')
    TD
    
    fprintf ('Transient norm square is:')
    transient_norm_square 
end

t= 0:1/800:ED-1/800;

threshold_vec1=threshold1*ones(1,length(y));

threshold_vec2=threshold2*ones(1,length(y));

vec_ver = -max(abs(y)):0.001:max(abs(y));
TD_ver =  (l+1)/Fs*ones(1,length(vec_ver));

if SW_PLOT  % added by Xu Chen, xuchen@cal.berkeley.edu
    figure
    plot(t,y);
    hold on
    plot(t,threshold_vec1,'r','LineWidth',2)
    hold on
    plot(t,-threshold_vec1,'r','LineWidth',2)
    hold on
    plot(t,threshold_vec2,'k','LineWidth',2)
    hold on
    plot(t,-threshold_vec2,'k','LineWidth',2)
    grid on
    xlabel('Time [sec]')
    ylabel('Residual force [V]')
    title('Adaptive disturbance rejection')
    
    
    figure
    %hold on
    plot(t,y);
    hold on
    plot(t,threshold_vec1,'r','LineWidth',2)
    hold on
    plot(t,-threshold_vec1,'r','LineWidth',2)
    hold on
    plot(t,threshold_vec2,'k','LineWidth',2)
    hold on
    plot(t,-threshold_vec2,'k','LineWidth',2)
    hold on
    plot(TD_ver,vec_ver,'-.k','LineWidth',2)
    hold on
    text(Td+TD,max(abs(y))*3/4,...
        ['Transient duration = ',num2str(TD)],...
        'HorizontalAlignment','center',...
        'BackgroundColor',[.7 .9 .7],...
        'FontSize',16);
    hold on
    text(Td+TD,-max(abs(y))*3/4,...
        ['Maximum transient value = ',num2str(maximum)],...
        'HorizontalAlignment','center',...
        'BackgroundColor',[.7 .9 .7],...
        'FontSize',16);
    grid on
    xlabel('Time [sec]')
    ylabel('Residual force [V]')
    title('Adaptive disturbance rejection')
    xlim([Td-0.1 Td+2*TD]);
    ylim([-max(abs(y)) max(abs(y))]);
end