function [pupper_env, plower_env, pmean_env, supper_env, slower_env, smean_env] = pchiptest(data_y)

% The function PCHIPTEST views envelopes generated using pchip and spline.
% Endpoint prediction is applied to fix the endpoints.
%
% Non MATLAB Library routines used are:
%	LOCALMIN, LOCALMAX, 
%	PREDICTENDPOINTS, PREDENDPOINTS.
%
% Example-
%	[upper_env, lower_env, mean_env] = pchiptest(data_y)
%
% Calling sequence-
% [pupper_env, plower_env, pmean_env, supper_env, slower_env, smean_env] = pchiptest(data_y)
%
% Input-
%	data_y		- the dataset for which the envelope
%			  is being generated
% Output-
%	pupper_env	- the upper envelope via pchip
%	plower_env	- the lower envelope via pchip
%	pmean_env	- the mean envelope via pchip
%	supper_env	- the upper envelope via spline
%	slower_env	- the lower envelope via spline
%	smean_env	- the mean envelope via spline
 
% Karin Blank (NASA GSFC)	    Initial
% Jelena Marshak (NASA GSFC)	May 5 2004 Modified
% (replaced the names of main and called non MATLAB functions as follows:
%	'pchip_test()' to 'pchiptest()'
%	'local_min()' to 'localmin',
%	'local_max()' to 'localmax',
%	'predict_endpoints()' to 'predictendpoints()',
%	'pred_endpoints()' to 'predendpoints()'.	 

data_x = 0:size(data_y, 2)-1;

[local_min_x, local_min_y] = localmin(data_x, data_y);
[local_max_x, local_max_y] = localmax(data_x, data_y);

if(data_y(1)>local_max_y(1))
    local_max_x=[data_x(1) local_max_x];
    local_max_y=[data_y(1) local_max_y];
end
if(data_y(1)<local_min_y(1))
    local_min_x=[data_x(1) local_min_x];
    local_min_y=[data_y(1) local_min_y];
end

if(data_y(end)>local_max_y(end))
    local_max_x=[local_max_x data_x(end)];
    local_max_y=[local_max_y data_y(end)];
end
if(data_y(end)<local_min_y(end))
    local_min_x=[local_min_x data_x(end)];
    local_min_y=[local_min_y data_y(end)];
end


%endpoint prediction
error1 = 0;
error2 = 0; 
[local_max_x, local_max_y, error1] = predictendpoints(local_max_x, local_max_y, length(data_x), data_y(1), data_y(end), 1);
[local_min_x, local_min_y, error2] = predictendpoints(local_min_x, local_min_y, length(data_x), data_y(1), data_y(end), -1);

    
if((error1 > 0) | (error2 > 0))
    disp('Using linear prediction')
    error1
    error2
    [temp_max_x, temp_max_y, temp_min_x, temp_min_y] = predendpoints(data_x, data_y, local_max_x, local_max_y, local_min_x, local_min_y);
end

if(error1 > 0)
    switch error1
    case 1
        local_max_x = [temp_max_x(1:2), local_max_x];
        local_max_y = [temp_max_y(1:2), local_max_y];
    case 2
        local_max_x = [local_max_x, temp_max_x(end-1:end)];
        local_max_y = [local_max_y, temp_max_y(end-1:end)];
    case 3
        local_max_x = [temp_max_x(1:2), local_max_x, temp_max_x(end-1:end)];
        local_max_y = [temp_max_y(1:2), local_max_y, temp_max_y(end-1:end)];
    end
end

if(error2 > 0)
    switch error2
    case 1
        local_min_x = [temp_min_x(1:2), local_min_x];
        local_min_y = [temp_min_y(1:2), local_min_y];
    case 2
        local_min_x = [local_min_x, temp_min_x(end-1:end)];
        local_min_y = [local_min_y, temp_min_y(end-1:end)];
    case 3
        local_min_x = [temp_min_x(1:2), local_min_x, temp_min_x(end-1:end)];
        local_min_y = [temp_min_y(1:2), local_min_y, temp_min_y(end-1:end)];
    end
end


%210 : Construct upper envelope of physical signal
supper_env = spline(local_max_x, local_max_y, data_x);
pupper_env = pchip(local_max_x, local_max_y, data_x);
%230 : Construct lower envelope of physical signal
slower_env = spline(local_min_x, local_min_y, data_x); 
plower_env = pchip(local_min_x, local_min_y, data_x);
%240 : Determine envelope mean of physical signal
temp = [pupper_env; plower_env];
pmean_env = mean(temp);

temp = [supper_env; slower_env];
smean_env = mean(temp);

disp('Figure 1 displays pchip envelopes');
disp('Figure 2 displays spline enveleopes');
figure(1);plot(data_x, data_y, data_x, pupper_env, data_x, plower_env, data_x, pmean_env);
figure(2);plot(data_x, data_y, data_x, supper_env, data_x, slower_env, data_x, smean_env);