predendpoints.m

一种新的时频分析方法的matlab源程序。

function [max_x, max_y, min_x, min_y] = predendpoints(data_x, data_y, ex_max_x, ex_max_y, ex_min_x, ex_min_y)
 
% The function PREDENDPOINTS appends predicted endpoints to data.
%	
% Non MATLAB library routines used are:
%	PREDICT, LOCALMAX, LOCALMIN.
%
% Calling sequence-
% [max_x, max_y, min_x, min_y] = predendpoints(data_x, data_y, 
%                       ex_max_x, ex_max_y, ex_min_x, ex_min_y)
%
% Input-
%	data_x		- vector representing the data x coordinates
%	data_y		- vector representing the data values
%	ex_max_x	- vector representing the max extrema x coordinates
%	ex_max_y	- vector representing the max extrema values
%	ex_min_x	- vector representing the min extrema x coordinates
%	ex_min_y	- vector representing the min extrema values
% Output-
%	max_x		- vector representing the max extrema x coordinates
%				with predicted endpoints appended 
%	max_y		- vector representing the max extrema values
%				with predicted endpoints appended
%	min_x		- vector representing the min extrema x coordinates
%				with predicted endpoints appended
%	min_y		- vector representing the min extrema values
%				with predicted endpoints appended
 
% Jelena Marshak (NASA GSFC)	May 3, 2004 Modified
%	(replaced the names of main and called non MATLAB functions:
% 	'pred_endpoints()' to 'predendpoints()',
% 	'local_min()' to 'localmin',
% 	'local_max()' to 'localmax'.

%primitive endpoint selection
%%max_x = [data_x(1), ex_max_x, data_x(end)];
%%max_y = [data_y(1), ex_max_y, data_y(end)];
%%
%%min_x = [data_x(1), ex_min_x, data_x(end)];
%%min_y = [data_y(1), ex_min_y, data_y(end)];
%%return;

%Dr. Huang's suggestion: reduce NPOLES% NPOLES = 10;
NPOLES = 4;

%set size of predicted data length to 2 times greatest extrema period
if(size(ex_max_x, 2) > size(ex_min_x, 2))
    NFUT = 2 * (ex_max_x(end) - ex_max_x(end-1));
%    NFUTb = 2 * (ex_max_x(2) - ex_max_x(1));
else
    NFUT = 2 * (ex_min_x(end) - ex_min_x(end-1));
%    NFUTb = 2 * (ex_min_x(2) - ex_min_x(1));
end

NFUTbak = NFUT;

%if # extrema < 3, use complete data set for prediction.
%else, use time from smallest x extrama value to data end
if((size(ex_max_x, 2) <= 3) | (size(ex_min_x, 2) <= 3))
    disp('less than three')
    
    %%experiment
   Y = data_y;
   TY = data_x;
    
%    Yb = fliplr(data_y);
    Yb = data_y;
    TYb = data_x;
    
%    Y = data_y(ex_max_x(1):ex_max_x(end));
%    TY = data_x(ex_max_x(1):ex_max_x(end));
%     
%    Yb = fliplr(Y);
%    TYb = TY;

else
    disp('more than three');
    %selection for end extrema
    
    %select last 3 end extrema
    if(ex_max_x(end-3) < ex_min_x(end-3))
        c = ex_max_x(end-3);
    else
        c = ex_min_x(end-3);
    end
    
 
    %for k=c: size(data_x, 2)
    %    Y = data_y(k);
    %    TY = data_x(k);
    %end
    
    Y = data_y(c:end);
    TY = data_x(c:end);
    
   
    %selection for beginning extrema
    
    if(ex_max_x(3) > ex_min_x(3))
        c = ex_max_x(3);
    else
        c = ex_min_x(3);
    end
    
    %for k=1: c
    %    Yb = data_y(k);
    %    TYb = data_x(k);
    %end
    
    Yb = data_y(1:c);
    TYb = data_x(1:c);
    
%    Yb = fliplr(Yb);
    
end

%actual start of emulation of C++ GetEndExtrema code

cdt = 1;
count = 0;

%fix for prediction bug - size Y must be >= NPOLES + 1
if(length(Y) < NPOLES+1)
    NPOLES = length(Y)-1;
end

%generate coefficients
a = lpc(Y, NPOLES);

%r = xcorr(Y);
%r(1:length(Y)-1) = [];       % Remove corr. at negative lags
%a = levinson(r,NPOLES);

b = polystab(a);

a = real(a);
b = real(b);


while ((cdt == 1) & (count < 3))
  
    %predict future values
%    ext_y = filter([0 -a(2:end)], 1, Y);

%    ext_y = real(ext_y);

    future = predict(Y, b, NFUT);    
 
    end_y = [data_y(end-1), data_y(end)];
    end_x = [data_x(end-1), data_x(end)];
    
        
    for k=1:NFUT
        end_y = [end_y, future(k)];
        end_x = [end_x, abs(end_x(1)-end_x(2))+end_x(end)];
    end
    
    %end_y
    %end_x
    
    [max_x, max_y] = localmax(end_x, end_y);
    [min_x, min_y] = localmin(end_x, end_y);

% %  figure(21);plot(end_x, end_y, max_x, max_y, data_x(end-100), data_y(end-100),TY, Y);
   
%     if((size(max_x,2) > 0) & (size(min_x,2) > 0))
    if((size(max_x,2) > 1) & (size(min_x,2) > 1))
        if ~((size(max_x, 2) == 1 & max_y(1)==0) | (size(min_x, 2) == 1 & min_y(1)==0))
            cdt = 0;
%%%test 6/21/02
% % %             maxEY = max_y(1);
% % %             maxET = max_x(1);
% % %             minEY = min_y(1);
% % %             minET = min_x(1);    
            if(length(max_y) >= 3) & (length(min_y) >=3)
                 maxEY = max_y(2:3);
                 maxET = max_x(2:3);
                 minEY = min_y(2:3);
                 minET = min_x(2:3);
             else
                 maxEY = [max_y(2), max_y(2)];
                 maxET = [max_x(2), max_x(2)+1];
                 minEY = [min_y(2), min_y(2)];
                 minET = [min_x(2), min_x(2)+1]; 
             end
%              maxEY = max_y(3:4);
%              maxET = max_x(3:4);
%              minEY = min_y(3:4);
%              minET = min_x(3:4);
        else    
            count = count + 1;
            NFUT = NFUT * 2;
        end
    else
        count = count + 1;
        NFUT = NFUT * 2;
    end
end

if(count == 3)
    disp('Could not find good end point, making a cheesy guess')
%    ext_y = filter([0 -a(2:end)], 1, Yb);

%    ext_y = real(ext_y);
    
%    minET = TY(1);
%    minEY = ext_y(TY(1));
%    maxET = minET;
%    maxEY = maxEY;

 %   [ext_x, ext_y] = localmax(TY, Y);
    
    maxEY = [mean(ex_max_y), mean(ex_max_y)];
    maxET = [ex_max_x(end)+2, ex_max_x(end)+4];
    
 %   [ext_x, ext_y] = localmin(TY, Y)
    
    minEY = [mean(ex_min_y), mean(ex_min_y)];
    minET = [ex_min_x(end)+2, ex_min_x(end)+4];
end

%end of endextrema C++ code

%start of begextrema C++ code emulation
 data = fliplr(Yb);

if(length(Y) < NPOLES+1)
    NPOLES = length(data)-1;
end 

 a = lpc(data, NPOLES);

 %r = xcorr(data);
 %r(1:length(data)-1) = [];       % Remove corr. at negative lags
%a = levinson(r,NPOLES);
 
 %a
 
 b = polystab(a);
  
 a = real(a);
 b = real(b);
 
 count = 0;
 cdt = 1;
 NFUT = NFUTbak;
 %NFUT = NFUTb;
 
 %y_temp = fliplr(data_y);
 while((cdt == 1) & (count < 3))
 
    future2 = predict(data, b, NFUT); 
        
    beg_y = [Yb(2), Yb(1)];
    beg_x = [TYb(2), TYb(1)];
    
    for k=1:NFUT
        beg_y = [beg_y, future2(k)];
 %       beg_x = [beg_x, beg_x(1) - abs(beg_x(2)-beg_x(1))];
        beg_x = [beg_x, -k];
    end
    
     
    
    %%beg_y = fliplr(beg_y);
    %%beg_x = fliplr(beg_x);
    
    %figure(2);plot(beg_x, beg_y, TYb, Yb);
    %figure(5); plot(end_x, end_y, TY, Y);
    %beg_y
    %beg_x
    
%    figure(4); plot(TYb, Yb);
  %  figure(3);plot(beg_x, beg_y , end_x, end_y, data_x, data_y);
    
    
    [max_x, max_y] = localmax(beg_x, beg_y);
    [min_x, min_y] = localmin(beg_x, beg_y);
    
 %   min_y
  %  max_y
    
    %figure(2);plot(max_x, max_y, min_x, min_y);
    
%     if((size(max_x,2) > 0) & (size(min_x, 2) > 0))
    if((size(max_x,2) > 1) & (size(min_x,2) > 1))
       if ~((size(max_x, 2) == 1 & max_y(1)==0) | (size(min_x, 2) == 1 & min_y(1)==0))
            cdt = 0;
%%%test change for 2 endpoints 6/21/02                   
% % %             maxBY = max_y(1);
% % %             maxBT = max_x(1);
% % %             minBY = min_y(1);
% % %             minBT = min_x(1);

    if(length(max_x) > 3) & (length(min_x) > 3)
        maxBY = max_y(2:3);
        maxBT = max_x(2:3);
        minBY = min_y(2:3);
        minBT = min_x(2:3);
    else
        maxBY = [max_y(2), max_y(2)];
        maxBT = [max_x(2), max_x(2)+1];
        minBY = [min_y(2), min_y(2)];
        minBT = [min_x(2), min_x(2)+1];
    end
        else
            disp('No extrema found 1');
            count = count + 1;
            NFUT = NFUT * 2;
        end
    else
        disp('No extrema found 2');
        count = count + 1;
        NFUT = NFUT * 2;
    end
end

if(count == 3)
    disp('Could not find good begin point, making a cheesy guess')
    %ext_y = filter([0 -a(2:end)], 1, Yb);

    %ext_y = real(ext_y);
    
    %minBT = TYb(end);
    %minBY = ext_y(end);
    %maxBT = minBT;
    %maxBY = maxBY;
    
 %   [ext_x, ext_y] = localmax(TYb, Yb);
    
    maxBY =[mean(ex_max_y), mean(ex_max_y)];
    maxBT = [ex_max_x(1) - 2, ex_max_x(1) - 4];
    
   % [ext_x, ext_y] = localmin(TYb, Yb)
    
    minBY = [mean(ex_min_y), mean(ex_min_y)];
    minBT = [ex_min_x(1)-2, ex_min_x(1)-4];
%     minBT = ext_x - 2;
          
end
    

%end of getBegExtrema

%maxBT
%maxBY
%maxET
%maxEY
%minBT
%minBY
%minET
%maxEY

max_x = [maxBT, ex_max_x, maxET];
max_y = [maxBY, ex_max_y, maxEY];
min_x = [minBT, ex_min_x, minET];
min_y = [minBY, ex_min_y, minEY];

% max_x = [maxBT, ex_max_x, maxET];
% max_y = [-ex_max_y(1), ex_max_y, maxEY];
% min_x = [minBT, ex_min_x, minET];
% min_y = [-ex_min_y(1), ex_min_y, minEY]

%max_x = real(max_x);
%max_y = real(max_y);
%min_x = real(min_x);
%min_y = real(min_y);

return;