function [P,T,th,Qd] = newPT(Qraw, factor, onsetTime, endTime);
% have the respiration signal in a column vector called Q
% call by >> [P,T,th,Qd] = newPT(Q);
% code adapted from Todd & Andrews. The Identification of Peaks in Physiological
% Signals. Computers and Biomedical Research 32, 322-335 (1999).

Qd = decimate(Qraw, 5);  % downsample the signal

d = 1;      % variable to show if signal unknown (1), trough-to-peak (2),
            % or peak-to-trough (3)
a = 1;      % index of maximal element since last trough
b = 1;      % index of minimal element since last peak
S = [];     % indices of maximal elements since last trough if signal rising (d = 2)
            % or minimal elements since last peak if signal falling (d = 3)
P = [];     % peak elements
T = [];     % trough elements
% factor is the number to use when making the threshold; default is 0.1

th = abs(prctile(Qd,75) - prctile(Qd,25)) * factor;
% th = abs(max(Qd) - min(Qd)) * .5;  % threshold

[n] = max(size(Qd));  % get size of the array so can do the for loop
for i = 1:n
    if d == 1
        if Qd(a) >= Qd(i) + th
            d = 3;
        elseif Qd(i) >= Qd(b) + th
            d = 2;
        end;
        if Qd(a) < Qd(i)
            a = i;
        elseif Qd(i) < Qd(b)
            b = i;
        end;
        S = i;
    elseif d == 2  % signal rising, trough-to-peak
        if Qd(a) < Qd(i) % still rising
            S = i; 
            a = i;
        elseif Qd(a) == Qd(i)
            S = [S,i];            
        elseif Qd(a) >= Qd(i) + th
            P = [P,S]; S = i;
            b = i; d = 3;
        end;
    elseif d == 3  % signal falling, peak-to-trough
        if Qd(i) <= Qd(b)
            S = i; 
            b = i;
        elseif Qd(b) == Qd(i)
           S = [S,i];
        elseif Qd(i) >= Qd(b) + th
            T = [T,S]; S = i;
            a = i; d = 2;
        end;
    end;
end;
        
% go through P and T, removing marks for peaks/troughs before the start point
% and after the end point
goodP = [];  goodT = [];
onsetTime = onsetTime/5; endTime = endTime/5;  % to match decimation
n = length(P);
for i = 1:n
     if P(i) > onsetTime
        if P(i) < endTime
           goodP = [goodP, P(i)];   % only add in peaks bigger than onsetTime
	end;
     end;
end;
n = length(T);
for i = 1:n
     if T(i) > onsetTime
        if T(i) < endTime
          goodT = [goodT, T(i)];   % only add in troughs bigger than onsetTime
	end;
     end;
end;

P = []; T = []; P = goodP; T = goodT;  % reset arrays for return values


plot(Qd, 'm'); whitebg([.9 .9 .9]);   % set background color to gray
hold on;  % plot the signal and peaks
plot(P, Qd(P), 'ob', 'MarkerSize',10);
plot(T, Qd(T), 'ob', 'MarkerSize',10);
hold off;