function [U]=wvfrms2bases(X,varargin)
% Computes an orthogonal basis from a given criteria.  Options are identical
% to those of sortwvfrms: basis are computed by signal width, localization, and
% noise.  The vectors are sorted by the given criteria, projected onto a
% subspace orthogonal to the previously sorted vectors, and sorted in that
% subspace, then projected, etc.  For example, in the case that 'width' is 
% used as a sorting option, this function first finds the most 'compact' 
% waveform among the entire set: that is, the waveform with with the smallest 
% width.  It then projects the matrix of data X onto a subspace orthogonal to 
% this waveform.  The function then finds the waveform with the smallest width 
% in this subspace, and projects the data onto a subspace orthogonal to it, and 
% so on.  The end result is a basis for X formed of vectors with progressively 
% broader signal widths. The pdf option finds the 'least noisy' data, as another
% example.
% 
% [U]=wvfrms2bases(X,varargin)
% 
% INPUT
%      X:  a matrix of whose columns form a set of linearly independent
%          signals/vectors.
% sortopt:  (optional) a string.  Can be either 'width', 'erg' or 'pdf'.
%          Default is 'width'.  This option is what sorts the data.  See
%          above for an explanation or website below.
% plotopt:  (optional) a string. 'plot' produces a plot.  Anything else
%           doesn't.
% 
% OUTPUT
% U: an orthonormal matrix.  The columns are ordered by the user specifed
%    criteria.  In the case that sortopt = 'width', the first column of U,
%    given by U(:,1), has the smallest signal width.  The next column has
%    the smallest width among all other vectors in the subspace of
%    functions orthogonal to U(:,1).
% 
% Examples:
% [U]=wvfrms2bases(W);
% [U]=wvfrms2bases(W,'plot');
% [U]=wvfrms2bases(W,'width','plot');
% [U]=wvfrms2bases(W,'erg','plot');
% [U]=wvfrms2bases(W,'pdf','plot');
% 
% See also http://www.ess.washington.edu/~joshuadc for downloading this.

[M,N]=size(X);
Ndisp=min(25,N);
plotopt='none';

if(nargin==1)
    sortopt='width';
end;

if(nargin==2)
    if(~strcmp(varargin{1},'plot'))
        sortopt=varargin{1};
    end;
    if(strcmp(varargin{1},'plot'))
        plotopt=varargin{1};
        sortopt='width';
    end;
end;

if(nargin==3)
    if(~strcmp(varargin{1},'plot'))
        sortopt=varargin{1};
    end;
    if(~strcmp(varargin{2},'plot'))
        sortopt=varargin{2};
    end;
    if(strcmp(varargin{1},'plot'))
        plotopt=varargin{1};
    end;
    if(strcmp(varargin{2},'plot'))
        plotopt=varargin{2};
    end;
end;

%for the first vector
U=zeros(M,N);
W=X;
[Is]=sortwvfrms(W,sortopt);
Ws=W(:,Is);
uhat=Ws(:,1)/norm(Ws(:,1));
U(:,1)=uhat;
P=eye(M,M)-uhat*uhat';
W=W(:,Is(2:end));
Wproj=P*W;

for k=2:N %not keeping track of Is correctly; 
    [Is]=sortwvfrms(Wproj,sortopt);%%?
    Ws=Wproj(:,Is);
    uhat=Ws(:,1)/norm(Ws(:,1));
    U(:,k)=uhat;
    P=(eye(M,M)-uhat*uhat')*P;
    Wproj=Wproj(:,Is(2:end)); 
    Wproj=P*Wproj;
end

if(strcmp(plotopt,'plot'))
    plotXmatrix(U(:,1:Ndisp));
    title(sprintf('%i Basis Vectors using %s Criteria',Ndisp,char(sortopt)));
    axis tight;
end;