function [TRresp RVT] = RVTcalc(Resp, tstep,TR)
%Tyler B Jones
%NIH/NIMH/LBC/FIM
%10/3/2006
%Resp = load('resppulse/respinterp.1D');
%tstep = 0.02;
%TR=2;
%Program takes respiration timecourse, calculates upper and lower envelopes
%and calculates the RVT function

%clc
%clear
% Offset = 5000;                              %set a nonzero offset
% Resp = Offset - load('Resp_4.txt');         %load ascii data

% tstep = .025;                               %define sampling period
time = tstep*[0:(length(Resp)-1)];          %define time vector

Max1 = max(Resp(1:(10/tstep)));                    %find first maximum
Max1Ind = min(find(Resp(1:(10/tstep)) == Max1));   %record first max index
Min1 = min(Resp(1:(10/tstep)));                    %find first minimum
Min1Ind = min(find(Resp(1:(10/tstep)) == Min1));   %record first min index

%Define Variables
MaxMin = 0;     %defines if current search is for maximum (0) or minimum (1)
Minima = [];    %defines minima vector as empty
Maxima = [];    %defines maxima vector as empty
CurrentInd = 0; %defines current extremum index
MinimaInd = []; %defines minima index vector
MaximaInd = []; %defines maxima index vector
CurrentMax = 0; %defines current max value
CurrentMin = 0; %defines current min value
TempInd = 0;    %defines temporary search index
SearchLen = 4/tstep;%defines how far to search from previous point (40Hz), Max allowable: 160 for 30bpm
Finished = 0;   %stop variable if indices are out of bounds

%Determine if first extremum is max or min
if Max1Ind > Min1Ind
    Minima = Min1;          %save first min in Minima vector
    CurrentInd = Min1Ind;   %update CurrentInd
    TempInd = CurrentInd;   %update TempInd
    MinimaInd = Min1Ind;    %update MinimaInd
    MaxMin = 1;             %set MaxMin to search for next maxima
else
    Maxima = Max1;          %save first min in Minima vector
    CurrentInd = Max1Ind;   %update CurrentInd
    TempInd = CurrentInd;   %update TempInd
    MaximaInd = Max1Ind;    %update MaximaInd
    MaxMin = 0;             %set MaxMin to search for next minima
end

%Search for Extremum Values and Indices
while CurrentInd < length(Resp)
    if MaxMin == 0                  %search for next minima
        for i =1:SearchLen          %search 3s into future
            if (CurrentInd + i) > length(Resp)
                Finished = 1;                   %break if indices out of bounds
                break
            end
            
            if Resp(CurrentInd + i) < Resp(TempInd)
                CurrentMin = Resp(CurrentInd + i);      %retrieve current min if less than previous
                TempInd = CurrentInd + i;               %update TempInd
            end
        end
        CurrentInd = TempInd;                           %update CurrentInd
        Minima = [Minima CurrentMin];                   %update minima vector
        MinimaInd = [MinimaInd CurrentInd];             %update minima index vector
        MaxMin = 1;                                     %change current search to maxima
    end

    if Finished == 1    %break search if indices out of bounds
        break
    end
    
    if MaxMin == 1
        for i = 1:SearchLen
            if (CurrentInd + i) > length(Resp)
                Finished = 1;
                break
            end
      
            if Resp(CurrentInd + i) > Resp(TempInd)
                CurrentMax = Resp(CurrentInd + i);      %retrieve current max if less than previous
                TempInd = CurrentInd + i;               %update current index
            end
        end
        CurrentInd = TempInd;                           %update CurrentInd
        Maxima = [Maxima CurrentMax];                   %update maxima vector
        MaximaInd = [MaximaInd CurrentInd];             %update maxima index vector
        MaxMin = 0;                                     %change current search to minima
    end

    if Finished == 1    %break search if indices out of bounds
        break
    end
end

%Calculate upper and lower envelope time vectors from index vectors
MaxTime = (MaximaInd - 1) * tstep;
MinTime = (MinimaInd - 1) * tstep;

%Calculate RVT
for i = 1:(length(MaxTime)-1)
    RVT(i) = (Maxima(i) - Minima(i))/abs(MaximaInd(i) - MinimaInd(i));
end

%Plot respiration timecourse with upper and lower envelopes
figure(1)
plot(time,Resp,'k',MaxTime,Maxima,'.r--',MinTime,Minima,'.g--')
xlabel('Time (sec)')
ylabel('Respiration (a.u.)')


RVTtime = MaxTime(1:(length(MaxTime)-1));   %Crop MaxTime to fit RVT length

%Plot RVT
figure(2)
plot(RVTtime,RVT,'r')
xlabel('Time (sec)')
ylabel('RVT (a.u.)')

%Interpolate RVT for a specific TR
TRtime = [0:TR:(length(time)*tstep)];
TRresp = interp1(RVTtime,RVT,TRtime);

%Plot interpolated RVT
figure(3)
plot(TRtime,TRresp,'r')
xlabel('Time (sec)')
ylabel('TR Interpolated RVT (a.u.)')