paper_calculate_period( t , y , whichTime , whichEndTime)

function [ period , deflectionMean , deflectionStd ] = paper_calculate_period( t , y , whichTime  , whichEndTime)
% PAPER_CALCULATE_PERIOD  Calculates period of a signal
%
% [ PERIOD DEFLECTIONMEAN DEFLECTIONSTD ] = PAPER_CALCULATE_PERIOD( t , y , whichTime )
% calculates the PERIOD of the (t,y) signal based on the time at which the
% WHICHTIMEth peak occurs. The mean value and standard deviation of the
% steady-state deflection (DEFLECTIONMEAN and DEFLECTIONSTD, respectively) are
% calculated based on the four cycles that begin at the WHICHTIMEth peak.
%
% [ PERIOD DEFLECTIONMEAN DEFLECTIONSTD ] = PAPER_CALCULATE_PERIOD( ... , whichEndTime )
% calculates the deflection based on the first four cycles of the last
% WHICHTIME peaks right before WHICHENDTIME

% Copyright 2006, The MathWorks, Inc.

    % Set the following flag to true to plot the signal and the limits used for the
    % period and steady-state calculations
    PLOT_ANALYSIS = false; 
    
    if PLOT_ANALYSIS
        figure; plot( t , y , 'b-' );hold on;title('Period Analysis');
        ybot=min(y);ytop=max(y);ydiff=(ytop-ybot)/10;
        ybotpd=ybot-2*ydiff;ytoppd=ytop+2*ydiff;
        ybot=ybot-ydiff;ytop=ytop+ydiff;
    end
    
    if nargin<4
        whichEndTime = [];
    end
    
    dy = diff(y);
    fdy = dy>0;
    dfdy = diff(fdy);
    wdfdy = find( dfdy==-1);
    
    %
    % Extract period
    %
    Twf=t( wdfdy( whichTime ) );
    Twb=t( wdfdy( 1 ) );
    period = (Twf - Twb) / ( whichTime - 1 );
    if PLOT_ANALYSIS
        plot( [Twf Twf] , [ybotpd ytoppd] , 'g-' );
        plot( [Twb Twb] , [ybotpd ytoppd] , 'g-' );
    end
    
    %
    % Extract steady-state deflection
    %    
    if isempty( whichEndTime )
        endIdx = length ( wdfdy );
    else

        % Find the peak that occurs on or before  whichEndTime
        [ dummy , timeIdx ] = min( abs( t - whichEndTime ) );
        whichEndPeaksIdx = find( wdfdy < timeIdx );
        endIdx = whichEndPeaksIdx( end );
        
    end
    
    % Get peaks near the end
    steadyStatePeakIdx = [
        wdfdy( endIdx-whichTime+4 )
        wdfdy( endIdx-whichTime+3 )
        wdfdy( endIdx-whichTime+2 )
        wdfdy( endIdx-whichTime+1 )
        wdfdy( endIdx-whichTime   )
                   ];
    
    % Average from approximately the midpoint of the oscillation
    steadyStateZeroIdx = round( diff( steadyStatePeakIdx )/4 ) + ...
                         steadyStatePeakIdx( 1:end-1 );
    for j=1:length( steadyStateZeroIdx )-1
        deflection( j ) = mean( y( steadyStateZeroIdx( j+1 ) : ...
                                   steadyStateZeroIdx( j )  ) ) ;
        if PLOT_ANALYSIS
            plot( [ t(steadyStateZeroIdx( j ))  t(steadyStateZeroIdx( j )) ] , ...
                  [ybot ytop] , 'r-' );
        end
    end
    
    deflectionMean = mean( deflection );
    deflectionStd  = std(  deflection );
    
    
    if PLOT_ANALYSIS
        plot( [ t(steadyStateZeroIdx( end ))  t(steadyStateZeroIdx( end )) ] , ...
              [ybot ytop] , 'r:' );
        keyboard; 
    end