[rt, pz_ix, pl_ix, ph_ix]=rise_time(y, Fs, percent_low, percent_high, peak_index, make_plot)

function [rt, pz_ix, pl_ix, ph_ix]=rise_time(y, Fs, percent_low, percent_high, peak_index, make_plot)
% % rise_time: 
% % 
% % Syntax:
% % 
% % [rt]=rise_time(y, Fs, percent_low, percent_high, peak_index, make_plot);
% % 
% % *********************************************************************
% % 
% % Description
% % 
% % Calculates rise time according to Federal Telecommunications Glossary.
% % Default is the time to rise from 10 percent to 90 percent of the 
% % signal peak.  There is an option for the program to output a plot of 
% % the peak 10 percent and 90 percent points and the previous zero.    
% % 
% % 
% % 
% % *********************************************************************
% % 
% % Input Variables
% % 
% % y                   % multichannel input time record in (Pa).  
% %                     % Processsing assumes that y has more channels 
% %                     % than time record samples.
% %                     % default is from analytic impoulse program.
% % 
% % Fs=50000;           % (Hz) sampling rate in Hz.  
% %                     % default is Fs=50000; (Hz).
% % 
% % percent_low=10;     % percentage of the range from 0 to the peak level
% %                     % 10 percent is the typical value from the Federal
% %                     % Telecommunications Glossary on rise time.  
% %                     % default is percent_low=10;
% % 
% % percent_high=90;    % percentage of the range from 0 to the peak level
% %                     % 90 percent is the typical value from the Federal
% %                     % Telecommunications Glossary on rise time.  
% %                     % default is percent_high=90;
% % 
% % peak_index=[];      % peak index for a one peak for each channel. 
% %                     % default is [maxy, peak_index]=max(y, [], 1);
% % 
% % make_plot=1;        % 1 Plots the processing data of the peaks.  
% %                     % Otherwise no plots are made. 
% %                     % 
% %                     % default is make_plot=0;  
% %
% %
% % *********************************************************************
% % 
% % Output Variables
% % 
% % rt is the rise time array (seconds).  Each channel has a single rise
% %         time.  Array size is [num_channels,1];
% % 
% % pz_ix is the array of the indices of the previous zero before the peak
% %         index for each channel.  Array size is [num_channels,1];
% % 
% % pl_ix is the array of the indices of the percent low crossing for each
% %         channel. Array size is [num_channels,1];
% % 
% % ph_ix is the array of the indices of the percent high crossing for each
% %         channel.  Array size is [num_channels,1];
% % 
% % 
% % *********************************************************************
% 
% Example='1';
% 
% % Each of the inputs should be a constant 
% 
% tau=0.001;        % seconds decay constant for modelling the ringing
% td=1;             % seconds time duration of the test signal
% A2=20;            % Pa initial ampitude of the of the sin wave
% A1=1;             % Pa amplitude of the background gaussian noise
% Fs=100000;        % Hz sampling rate of the impulsive noise
% fc=1000;          % Hz ringing frequency of sin wave
% delay=0.1;        % seconds time delay between impulsive noise 
%                   % and test signal recording
% rt=0.0005;        % seconds time to reach peak value from back ground
%                   % noise level
% 
% [y, t]=analytic_impulse(Fs, fc, td, tau, delay, A1, A2, rt);
% 
% percent_low=10;
% percent_high=90; 
% peak_index=[];
% make_plot=1;
% 
% [rt]=rise_time(y, Fs, percent_low, percent_high, peak_index, make_plot);
% 
% 
% 
% Example='2';
% 
% tau=0.001;            % seconds decay constant for modelling the ringing
% td=0.2;               % seconds time duration of the test signal
% A2=20;                % Pa initial ampitude of the of the sin wave
% A1=1;                 % Pa amplitude of the background gaussian noise
% Fs=800000;            % Hz sampling rate of the impulsive noise
% fc=40000;             % Hz ringing frequency of sin wave
% delay=0.05;           % seconds time delay between impulsive noise 
%                       % and test signal recording
% rt=0.00005;           % seconds time to reach peak value from back ground
%                       % noise level (variation of rise time)
% 
% % [y, t]=analytic_impulse(Fs, fc, td, tau, delay, A1, A2, rt);
% 
% 
% percent_low=10;
% percent_high=90; 
% peak_index=[];
% make_plot=1;
% 
% [rt]=rise_time(y, Fs, percent_low, percent_high, peak_index, make_plot);
% 
% % *********************************************************************
% % 
% % References
% % 
% %     From the electronic Federal Telecommunications Glossary.
% %     rise time: In the approximation of a step function, the time 
% %     required for a signal to change from a specified low value to 
% %     a specified high value. Typically, these values are 10% and 90% 
% %     of the step height. (188) 
% % 
% %     http://www.its.bldrdoc.gov/fs-1037/fs-1037c.htm
% % 
% % 
% % *********************************************************************
% % 
% % 
% % Subprograms
% % 
% % 
% % 
% % List of Dependent Subprograms for 
% % rise_time
% % 
% % FEX ID# is the File ID on the Matlab Central File Exchange
% % 
% % 
% % Program Name   Author   FEX ID#
% % 1) analytic_impulse		Edward L. Zechmann			
% % 2) convert_double		Edward L. Zechmann			
% % 3) find_previous_crossing       Edward L. Zechmann	
% % 4) geospace		Edward L. Zechmann			
% % 
% % 
% % 
% % *********************************************************************
% %
% % rise_time is written by Edward L. Zechmann 
% %  
% %     date  4 August      2010    Created program
% %  
% % modified  9 August      2010    Updated Comments
% % 
% % modified 10 August      2010    Updated Comments
% % 
% % modified 30 September   2010    Fixed a bug in plotting the peak 
% %                                 indices.  Updated Comments
% % 
% %
% % 
% % 
% %  
% % *********************************************************************
% % 
% % Please feel free to modify this code.
% % 
% % See also: A_duration, B_mil_1474D_duration, B_Duration, C_Duration,
% %           D_Duration
% % 


if nargin < 1 || isempty(y) || ~isnumeric(y)
    tau=0.001;        % seconds decay constant for modelling the ringing
    td=1;             % seconds time duration of the test signal
    A2=20;            % Pa initial ampitude of the of the sin wave
    A1=1;             % Pa amplitude of the background gaussian noise
    Fs=100000;        % Hz sampling rate of the impulsive noise
    fc=1000;          % Hz ringing frequency of sin wave
    delay=0.1;        % seconds time delay between impulsive noise 
    %                 % and test signal recording
    rt=0.0005;        % seconds time to reach peak value from back ground
    %                 % noise level
    % 
    [y]=analytic_impulse(Fs, fc, td, tau, delay, A1, A2, rt);

end

% Make the data have the correct data type and size
[y]=convert_double(y);

% Make sure the matrix y is oriented correctly.  
% Transpose if necessary.  
[num_samples, num_channels]=size(y);

if num_samples < num_channels
    y=y';
    [num_samples, num_channels]=size(y);
end

if nargin < 2 || isempty(Fs) || ~isnumeric(Fs)
    Fs=50000;
end

if nargin < 3 || isempty(percent_low) || ~isnumeric(percent_low)
    percent_low=10;
end

% Set the settling time of the filters default is 0.1 seconds.
if (nargin < 4 || isempty(percent_high)) || ~isnumeric(percent_high)
    percent_high=90;
end

if (nargin < 5 || isempty(peak_index)) || any(~isnumeric(peak_index))
    dim_size=size(y);
    [buf max_dim]=max(dim_size);
    [maxy, peak_index]=max(abs(y), [], max_dim);
end

if (nargin < 6 || isempty(make_plot)) || ~isnumeric(make_plot)
    make_plot=0;
end


rt=zeros(num_channels, 1);
pz_ix=zeros(num_channels, 1);
pl_ix=zeros(num_channels, 1);
ph_ix=zeros(num_channels, 1);
p_ixa=zeros(num_channels, 1);

for e1=1:num_channels;
    
    if length(peak_index) >= e1
        p_ix=peak_index(e1);
    else
        p_ix=peak_index(end);
    end

    p_ixa(e1)=p_ix;

    if ~isempty(p_ix)
        
        ypeak=y(p_ix, e1);
        

        [last_zero]=find_previous_crossing(y(1:p_ix, e1), p_ix, 0);   
        
        if isempty(last_zero)
            last_zero=1;
        end     

        
        % Percent low index
        [percent_low_ix, at1]=find_previous_crossing(y(1:p_ix, e1), p_ix, percent_low/100*ypeak);  
        
        if isempty(last_zero)
            percent_low_ix=1;
        end     

        % Percent high index
        [percent_high_ix, at2]=find_previous_crossing(y(1:p_ix, e1), p_ix, percent_high/100*ypeak); 
        
        if isempty(last_zero)
            percent_high_ix=1;
        end     
        
        pz_ix(e1)=last_zero;
        pl_ix(e1)=percent_low_ix;
        ph_ix(e1)=percent_high_ix;
        
        % rise time array seconds
        rt(e1)=(at2-at1)/Fs;
        
    else
        rt(e1)=NaN;
    end
    
end


if isequal(make_plot, 1)
    
    Lt=2;
    
    close('all');
    figure(1);
    
    for e1=1:num_channels;
        
        
        subplot(num_channels, 1, e1);
        plot( 1/Fs*(-1+(1:num_samples)), y(:, e1));
        hold on;
        plot( 1/Fs*(-1+pz_ix(e1)), y(pz_ix(e1), e1), 'ok', 'linestyle', 'none', 'markersize', 7, 'LineWidth', Lt);
        plot( 1/Fs*(-1+pl_ix(e1)), y(pl_ix(e1), e1), 'vm', 'linestyle', 'none', 'markersize', 7, 'LineWidth', Lt);
        plot( 1/Fs*(-1+ph_ix(e1)), y(ph_ix(e1), e1), 'sg', 'linestyle', 'none', 'markersize', 7, 'LineWidth', Lt);
        plot( 1/Fs*(-1+p_ixa(e1)), y(p_ixa(e1), e1), 'or', 'linestyle', 'none', 'markersize', 7, 'LineWidth', Lt);
    
    end
    
    hold off;
    
end