[X,Y]=maxall(varargin)

function [X,Y]=maxall(varargin)
% MAXimum points in a curve, ALL(default) or manually selected
% with (default) or without parabolic interpolation at extremes
% needs CLINE
%Call: 
%       [X,Y]=maxall(xy  [,INTERSELECT]  )
%       [X,Y]=maxall([INTERPSELECT,]   xy)
%       
%Input:
%   xy = data (coordinates [x, y] or [x;y]) or x+1i*y with distinct x         x
%           or vector of y - values, then x=1:length(y)
%           at least 3 points are necessary
%           i.e.  minimal acceptable length(xy) is 3
%           xy-points can be entered in any x-sequence, NaN ignored
%             
%           INTERPSELECT: switch for interpolation/selection/visualization 
%
%           INTERPSELECT=INTERP +1i*SELECT
%           INTERP: switch for interpolation,
%           if INTERP: parabolic interpolation at internal maxima (default)
%
%           SELECT: switch for manual selection/visualization of maxima
%           if SELECT:  manual positive (SELECT>0) or negative(=deletion) (SELECT<0) selection
%           Visualization only, all maxima: cancel selection by clicking
%           left or right of the plot area (possible at any selection step)
%           
%Output:
%	       X,Y = position(s) and value(s) of maximum point(s), 
%           allowed maxima at the beginning and end. 
%          
%Subfunctions:
%       bip = "0 of sound bit" (pib = bip(1) = "1 of sound bit")
%       o_o = sound "word", expressing pity ("bemitleiden") if error occurs
%       [x,y]=absord(xy) = makes two columns of x,y coordinates
%       nei=neighbor(x,X) = indexes to X elements  closest to x elements   
%       [X,Y]=parextrem(x,y,IND) finds extreme positions and  values of
%       parabolas drawn through every 3 points 
%               [x(IND-1:IND+1),y(IND-1:IND+1)] 
%        IND indicates extremal points: 1 < IND < length(x)
%       [X,Y]=ginselect(X,Y,SELECT) positive or negative selection via
%       clicking nearest abscissas
%
%DEMO 1: smooth data
%x=0:.1:100;
%y=sin(x).*exp(-.01*x);
%xy =[x;y];
%[X,Y]=maxall(y);
%[XX,YY]=maxall(1,x+1i*y);
%norm(Y-YY)
%
%%DEMO 2: noisy data
% clf
% r=savfilt(randn(1,100),10); %c.f. SAVFILT in FEX
% [X,Y]=maxall(r,1i);
% [X,Y]=maxall(1-1i,r);

%Vassili Pastushenko 8th June 2004
%first revision 23-rd  Feb 2006
% second revision 3-rd apr. 2007, back to "num2str" from earlier
% more functional "tvar", thanks to email of A.D. Nair
%=========================================================
ARGS=varargin;
    %Define arguments
    LARG=length(ARGS);
    switch LARG
        case 0
            o_o
            pause(.2)
            error('Missing data');        
        case 1
            xy=ARGS{1};
            if length(xy)<3
              checkpoint
              error('at least 3 points expected')
            end
            
            INTERPSELECT=1;
            
        case 2
            ind=1:2;
            for i=1:LARG
                if length(ARGS{i})>2
                    xy=ARGS{i};
                    break
                end
            end
            
            if exist('xy')
            ind(i)=[];
                INTERPSELECT=ARGS{ind};
            else
                checkpoint
                error('At least 3 points expected')
            end
        otherwise
            o_o
            pause(.3)
            error('Check arguments: Data and INTERPSELECT')
    end
            
            
    %Prepare data: absord = abs_cissa and ord_inate
    [x,y]=absord(xy);
LENDAT=length(x);
%Padding boundaries for unification
x=[2*x(1)-x(2);x;2*x(LENDAT)-x(LENDAT-1)];
y=[y(2);y;y(LENDAT-1)];
LENDAT=LENDAT+2;
%Detecting all maxima
D=diff(y);
D=[0; D];
IND=1:LENDAT;
IND=IND(D(1:LENDAT-1)>0&D(2:LENDAT)<0);

%All Max 
X=x(IND);
Y=y(IND);

%Selection / Interpolation
INTERP=real(INTERPSELECT);
SELECT=imag(INTERPSELECT);

if SELECT
    Index=2:LENDAT-1;
    plot(x(Index),y(Index),'.:')
    [X,Y,IND]=ginselect(X,Y,SELECT,IND);    
    
end

if INTERP
[X,Y]=parextrem(x,y,IND);    
end

if SELECT&INTERP
    line(X,Y,'marker','o','markersize',7,'linestyle','none','color','r','linewidth',3)
end
   
if SELECT
   if INTERP
       RES=' interpolated';
   else
       RES=' rigid';
   end
   LENX=length(X);
   many=LENX>1;
   ad=char(('um'*~many)+'a '*many);
    title({['Detected ',num2str(LENX),RES,' maxim',ad];'  '})
    legend('Data','All Max','rigid Max','interp Max','location','best')
    figure(gcf)
end
   

function  o_o(pib)
%   
%			
%	Vassili Pastushenko	 2005
%==============================
if nargin
         S=bip(1);
else
    S=bip;
end
         S=[S;zeros(1200,1);flipud(S)];
         sound(S);
         
         function  S=bip(pib)
%Sound warning
%Call:
%           bip
%           bip(1)
%Vassili Pastushenko	Jul	2004
%==============================
t=(1:.1:100)';
ONETWO=linspace(1,3,length(t))';
S=cos(t.*ONETWO)./(2+ONETWO);
if nargin>0
    S=flipud(S);
end
if nargout<1
sound(S);
clear S;
end
function checkpoint
    o_o;
    display('Check points');
      
    function nei=neighbor(x,X)
    %x and X are real or complex vectors
    %  nei = index to X, length(nei)=length(x)
    %nei(i): X(nei(i)) is closest of all X to x(i):
    %[dum, nei(i)] = min(abs(x(i)-X));
    % therefore if x = "grooms", then X(nei)= found "neighbors"
    XX=repmat(X(:),1,length(x));
    xx=repmat(x(:)',length(X),1);
    DIST=abs(XX-xx); % distances to points i
    [dum,nei]=min(DIST);
 
  function [x,y]=absord(xy)
 %ABScissa's (x=indep. var.) and ORDinates (y=dep. var.)of a curve
      if isempty(xy)
         checkpoint
         error('Data missing')
     end
      S=size(xy);
      if length(S)>2|min(S)>2
          checkpoint
         error('2D vector expected')
      end
         
      MS=min(S);
    switch MS
        case 1
            if isreal(xy)% e.g. missing x:  xy = y 
                xy=[(1:length(xy))',xy(:)];% x as index to y
            else % complex coordinates entered
                xy=xy(:);
                xy=[real(xy),imag(xy)]; %matrix form
            end
        case 2
            if isreal(xy)
                if S(1)==2
                    xy=xy';
                end
            else  
                checkpoint
                error('Only complex vectors expected')
            end
        otherwise
            checkpoint;
            error('Check data')
    end

xy(any(isnan(xy')),:)=[];
xy=sortrows(xy,1);
x=xy(:,1);
if any(diff(x)==0)
disp('Distinct points expected !')
checkpoint
error('!')
end
y=xy(:,2);
%============================================
function [X,Y,IND]=ginselect(X,Y,SELECT,IND)
    NAR=nargin;
    if NAR<3
        SELECT=1;
    end
    if NAR<4
        IND=1:length(X);
    end
    
    axis tight
        
    h=cline(X+1i*Y,'go3');
    set(h,'linewidth',4);
    set(gca,'Fontsize',15)       
    XLI=get(gca,'XLim');
    SAV=SELECT>0;
    TEX=char(round('retained'*SAV+'deleted '*(1-SAV)));
    title({['Click abscissas of points to be ',TEX];'select all: "ENTER" or click left or right of the plot area'})
    xg=ginput;
    if ~isempty(xg)
    xg=xg(:,1);
    if any(xg<XLI(1)) | any(xg>XLI(2))
        xg=[];
    end
    end
 if ~isempty(xg)
    INDEX=neighbor(xg,X);
    if SELECT>0
        X=X(INDEX);
        Y=Y(INDEX);
        IND=IND(INDEX);
    else
        X(INDEX)=[];
        Y(INDEX)=[];
        IND(INDEX)=[];
    end
 end
  h=cline(X+1i*Y,'9b^');
    set(h,'linewidth',2)  
    %=================================================
    function [X,Y]=parextrem(x,y,IND)
%assumptions:
%x=x(:);y=y(:);isreal(x),length(x)=length(y),
%distinct x, i.e. all(abs(diff(x))>0) is true

        if nargin<3
            IND=2:length(x)-1;
        end
    DELS=x(IND+1)-x(IND-1);
    X=x(IND);
    Y=y(IND);
    xx=(X-x(IND-1))./DELS;
    YY=[y(IND-1)./xx, Y./(xx.*(xx-1)), y(IND+1)./(1-xx)];
    t=(((xx+1).*YY(:,1)+ YY(:,2)+ xx.*YY(:,3))./(sum(YY')'*2));
    Y=(t-xx).*(t-1).*YY(:,1)+ t.*(t-1).*YY(:,2)+ t.*(t-xx).*YY(:,3);
    X=x(IND-1)+t.*DELS;