predictendpoints.m

这是一个关于hht变换很有用的工具箱

function [x, y, error] = predictendpoints(extrema_x, extrema_y, max_length, beg_inter, end_inter, env_type)

% The function PREDICTENDPOINTS predicts 2 endpoints based on extrema.
%
% Note: predicted extrema are not added to the end in question.
%
% Calling sequence-
% [x, y, error] = predictendpoints(extrema_x, extrema_y,
%                 max_length,beg_inter, end_inter, env_type)
%
% Input-
%	extrema_x	- the points to base prediction on
%	extrema_y	- the points to base prediction on	
%	max_length	- the maximum length of the complete dataset
%	beg_inter	- beginning point of the dataset 
%	end_inter	- last point of the dataset
%	env_type	- indicator: -1 if calculating for minimum extrema
%			  and 1 if calculating for maximum extrema
% Output-
%	x		    - the set of extrema with the end and begining
%			    predicted endpoints tacked on.
%	y		    - the set of extrema with the end and begining
%			    predicted endpoints tacked on.
%	error		- an indicator of the algorithm being unable
%			    to calculate an endpoint without adding
%			    an additional extrema to the dataset.
%			    It is 0 if no error, 
%			    1 if error in begin points, 
%			    2 if error in future end points,
%			    3 if in both.
 
% Jelena Marshak (NASA GSFC)	May 3, 2004 Modified
%	(replaced the name of the function as follows:
% 	'predict_endpoints()' to 'predictendpoints()' ).

disp('.Predicting endpoints')

error = 0;

if (size(extrema_x, 2) <= 10)
    if(size(extrema_x, 2) > 4)
        disp('...Abbreviated Prediction');
        end_x = extrema_x;
        end_y = extrema_y;
        
        beg_x = extrema_x;
        beg_y = extrema_y;
    else
        disp('...Ugly guess');
        x = [-5, 0, extrema_x, max_length, max_length+5];
        y = [extrema_y(1), extrema_y(1), extrema_y, extrema_y(end), extrema_y(end)];
   
        return;
    end
else    
    end_x = extrema_x(end-10:end);
    end_y = extrema_y(end-10:end);
    
    beg_x = extrema_x(1:10);
    beg_y = extrema_y(1:10);
end

%start begin point prediction

%calculate slopes and intercepts

%find slope/intercept relationships between all points
for i=1:size(beg_x, 2)
    for j=1:size(beg_x, 2)
        if(~(i == j) & ~(j+1 == i))
            slope(i, j) = (beg_y(j) - beg_y(i))/(beg_x(j) - beg_x(i));
            intercept(i, j) = beg_y(j) - (slope(i, j) * beg_x(j));
         end
    end
end

%calculate the standard deviation in relation to first point
for j=1:size(beg_x, 2)
        std_slope(j) = std([sum(slope(1, :)); sum(slope(j, :))]);
        std_intercept(j) = std([sum(intercept(1, :)); sum(intercept(j, :))]);
end
temp = length(std_intercept);
%score points (add absolute values of stds together)
std_all = abs(std_slope); 

%sort by score
[sorted_std, all_match] = sort(std_all, 2);

%set to 2nd point (first is prediction base)
current_match = 1;
y_test = beg_inter - env_type;
match = 0;

while( ((y_test > beg_inter) & (env_type == -1)) | ((y_test < beg_inter) & (env_type == 1)))
    
    if(((y_test > beg_inter) & (env_type == -1)) | ((y_test < beg_inter) & (env_type == 1)))
        if(std([y_test, beg_inter]) < .05)
            break;
        end
    end
    
    current_match = current_match + 1;
    
    if(current_match > length(all_match))
        match = all_match(2);
    else
        match = all_match(current_match);
    end
        
    %calculate change x
    delta_x = beg_x(match) - beg_x(match-1);
    
    final_slope = (beg_y(match-1)-beg_y(match))/(beg_x(match-1) - beg_x(match));
    
    %calculate intercept
    B_intercept = beg_y(match) - (final_slope*beg_x(match));
    tB = (B_intercept-beg_y(match))/beg_x(match);
    final_intercept = (tB * beg_x(1)) + beg_y(1);%%%%TEST%%%%+ beg_y(match);
    
    %calculate x and y values
    x_val = beg_x(1) - delta_x;
    y_val = x_val * final_slope + final_intercept;
    
    if(current_match > length(all_match))
       break;
    end
    
    
    y_test = 0 * final_slope + final_intercept;
end

beg_x = [x_val, beg_x];
beg_y = [y_val, beg_y];

if(match < 1)
    match = 2;
end

    
%calculate change 2nd x
delta_x = beg_x(match) - beg_x(match-1);

%final_slope = (beg_y(match-1)-beg_y(match))/(beg_x(match-1) - beg_x(match));
final_slope = (beg_y(match)-beg_y(match-1))/(beg_x(match) - beg_x(match-1));

%calculate 2nd intercept
B_intercept = beg_y(match) - (final_slope*beg_x(match));
tB = (B_intercept-beg_y(match))/beg_x(match);
final_intercept = (tB * beg_x(1)) + beg_y(1);%%%%TEST%%%%+ beg_y(match);

%calculate 2nd x and y values
x_val2 = beg_x(1) - delta_x;
y_val2 = x_val2 * final_slope + final_intercept;

if((x_val2 <= 0)) 
    x = [x_val2, x_val, extrema_x];
    y = [y_val2, y_val, extrema_y];
else
    disp('...Begin point prediction error');
    x_val2
    error = 1;
    x = extrema_x;
    y = extrema_y;
end

%end begin point prediction

%start end point prediction
%find slope/intercept relationships between all points
for i=1:size(end_x, 2)
    for j=1:size(end_x, 2)
         if(~(i == j) & ~(j-1 == i))
%        if(~(i == j) & ~(j+1 == i))
            slope(i, j) = (end_y(j) - end_y(i))/(end_x(j) - end_x(i));
            intercept(i, j) = end_y(j) - (slope(i, j) * end_x(j));
        end
    end
end

%calculate the standard deviation in relation to end point
for j=1:size(end_x, 2)
    std_slope(j) = std([sum(slope(end, :)); sum(slope(j, :))]);
    std_intercept(j) = std([sum(intercept(end, :)); sum(intercept(j, :))]);
end

%score points (add absolute values of stds together)
std_all = abs(std_slope);% + abs(std_intercept);

%sort by score
[sorted_std, all_match] = sort(std_all, 2);

%set to 2nd point (first is prediction base)

current_match = 1;
y_test = end_inter - env_type;
match = 0;
%error_temp = 0;

while( ((y_test > end_inter) & (env_type == -1)) | ((y_test < end_inter) & (env_type == 1)))  
    
    %tests that it's not repeating for minor difference
    if(((y_test > end_inter) & (env_type == -1)) | ((y_test < end_inter) & (env_type == 1)))
        if(std([y_test, end_inter]) < .05)
            break;
        end
    end
    
    
    current_match = current_match+1;
     if(current_match+1 > length(all_match))
         match = all_match(2);
     else
        match = all_match(current_match);
    end
 
    if(match + 1 > length(all_match))
        match = all_match(current_match-1);
    end
    
    %calculate change x
    delta_x = end_x(match+1) - end_x(match);
    
    final_slope = (end_y(match+1) - end_y(match))/(end_x(match+1) - end_x(match));
      
    %calculate intercept
    B_intercept = end_y(match) - (final_slope*end_x(match));
    B = atan((B_intercept-end_y(match))/end_x(match));
    final_intercept = (tan(B) * end_x(end)) + end_y(end);%%%%TEST%%%%end_y(match);
    
    %calculate x and y values
    x_val = end_x(end) + delta_x;
    y_val = x_val * final_slope + final_intercept;
    
    if(current_match+1 > length(all_match))
        break;
    end
    
    y_test = max_length * final_slope + final_intercept;
end

end_x = [end_x, x_val];
end_y = [end_y, y_val];

%calculate change x for 2nd point
match = match + 1;

delta_x = end_x(match+1) - end_x(match);

%calculate slope for 2nd point
%final_slope = (end_y(match+1) - end_y(match))/(end_x(match+1) - end_x(match));
final_slope = (end_y(match) - end_y(match+1))/(end_x(match) - end_x(match+1));

%calculate intercept for 2nd point
B_intercept = end_y(match) - (final_slope*end_x(match));
B = atan((B_intercept-end_y(match))/end_x(match));
final_intercept = (tan(B) * end_x(end)) +  end_y(end);%%%%TEST%%%%end_y(match);

x_val2 = x_val + delta_x;
y_val2 = x_val2 * final_slope + final_intercept;

if((x_val2 < max_length)) 
    disp('...End point prediction error');
    error = error + 2;
else
    x = [x, x_val, x_val2];
    y = [y, y_val, y_val2];
end

%end endpoint prediction

return;