function [varargout]=sortwvfrms(W,varargin)
% Sorts input matrix of data so that each column is re-ordered by
% user-specified criteria.  Options are (1) sort data by signal width (2)
% sort data by energy localization (3) sort data by dissimilarity of
% data to a normal probability distribution.  The default option is to 
% return all 3 options, as vectors of indices so that W(:,Iw), for example, 
% gives the data ordered by signal width.
% 
% [varargout]=sortwvfrms(W,varargin)
% 
% INPUT 
%      W: A matrix whose columns are time series to be sorted.
% sortopt: (Optional) A string.  Can be either 'width', 'erg' or 
%         'pdf'.  Default is all three.
% plotopt: (Optional) A string.  'plot' produces a plot.  Anything
%          else doesn't. 
% Inputting 'all' returns a plot of all 3 sorting results.
% 
% OUTPUT
% I: a vector of indices that sorts the time series, so that W(:,I) has it's
%   columns ordered by the specified criteria.
% 
% Examples:
% [Iw]=sortwvfrms(W,'width');
% Wsort=W(:,Iw); %now contains the sorted signals.
% 
% [Iw,Ie,Ip]=sortwvfrms(W,'all');
% %Returns indices that sort the input data by signal width, energy
% localization, and 'least noisy'.
% 
% [Ip]=sortwvfrms(W,'pdf','plot');
% %Now W(:,Ip) is ordered from least noisy to most noisy signal.  A plot is
% produced too.
% 
% NOTE: Code requires plotXmatrix as well.  Download from file exchange, or
% http://www.ess.washington.edu/~joshuadc


[M,N]=size(W);
tau=zeros(N,1);  %initialize signal width vector
t0=floor(M/2); %mean time value
wvth=zeros(t0,1); %waveform width function to be tested using Wbar
Wbar=0.90; %thresh-hold signal width to be set
Ndisp=min(25,N); %Number of waveforms to display on plot
Erg=zeros(N,1); %initialize signal energy record
%time=linspace(0,M*sampInt,M);
plotopt='none';

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

if(nargin==2)
    if(~strcmp(varargin{1},'plot'))
        sortopt=varargin{1};
    end;
    if(strcmp(varargin{1},'plot'))
        plotopt=varargin{1};
        sortopt='all';
    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;

if(strcmp(sortopt,'all')||(strcmp(sortopt,'width')...
        ||strcmp(sortopt,'erg')))

    for k=1:N
        %%%% COMPUTE SIGNAL WIDTH
        V=W(:,k).^2/(norm(W(:,k).^2));%normalize energy
        [maxwav,tbar]=max(V);%max time value is tbar
        Erg(k)=maxwav;%max value of V stored
        V=circshift(V, (t0-tbar) );%wave shifted to middle of plot.

        for n=0:t0-1 %compute wavewidth as function of n
            wvth(n+1)=trapz(V(t0-n:t0+n));
        end;

        %find value of wavewidth that exceeds Wbar, and the value of tau
        [minv,Itmp]=min(abs(wvth - Wbar));
        tau(k)=Itmp; %value of tau giving Wbar% of signal energy
    end;

    %The minimum values of tau give the most 'compact' waveforms; the
    %waveform list can be sorted from most compact waveform to least
    %compact waveform using sort and permute.
    [sortau,Iw]=sort(tau);
    varargout{1}=Iw;
end;

%Plot results if option invoked.
if((strcmp(plotopt,'plot')&&(strcmp(sortopt,'width')))||strcmp(sortopt,'all'))

    Wcmp=W(:,Iw); %Waveforms sorted by waveform width.
    figure;
    plotXmatrix(Wcmp(:,1:Ndisp),'align');
    title(sprintf('%i of %i Waveforms Sorted By Signal Width',Ndisp,N));
    ylabel(sprintf('%i Waveforms with Smallest Signal Width',Ndisp));
    xlabel('Sample Number');
    axis tight;

end;
%%%% END SIGNAL WIDTH COMPUTATION

if(strcmp(sortopt,'all')||strcmp(sortopt,'erg'))
    %%%% SORT SIGNALS BY AMOUNT OF ENERGY AND LOCALIZATION
    Erg=Erg./tau;
    [sortErg,Ierg]=sort(Erg,1);
    Ierg=flipud(Ierg);
    if(strcmp(sortopt,'erg'))
        varargout{1}=Ierg;
    end;
    if(strcmp(sortopt,'all'))
        varargout{2}=Ierg;
    end;
end;

if((strcmp(plotopt,'plot')&&(strcmp(sortopt,'erg')))||strcmp(sortopt,'all'))

    Werg=W(:,Ierg);
    figure;
    plotXmatrix(Werg(:,1:Ndisp),'align');
    title(sprintf('%i of %i Waveforms Sorted By Energy Localization',Ndisp,N));
    ylabel(sprintf('%i Waveforms most Energetically Localized',Ndisp));
    xlabel('Sample Number');
    axis tight;

end;
%%%% END ENERGY LOCALIZATION SORTING
%% COMPUTE WAVEFORM HISTOGRAM FOR EACH WAVEFORM
%Compute the fit of a normal distribution to the data, then perform a
%chi-square test
if(strcmp(sortopt,'all')||strcmp(sortopt,'pdf'))

    numbins=floor(M*0.06);
    mj=zeros(numbins,N);
    Hemp=zeros(size(mj));
    pdfdev=zeros(N,1);
    for k=1:N
        %define bins to compute fit to normal distribution
        [H,X]=hist(W(:,k),numbins);
        Hemp(:,k)=H;
        t=linspace(min(W(:,k)),max(W(:,k)),M);
        sigma=std(W(:,k));
        mu=mean(W(:,k));
        pdf=1/(sigma*sqrt(2*pi))*exp(-0.5*((t-mu)/(sigma)).^2);
        for j=1:(length(X)-1)
            Ij=find((t > X(j))&(t < X(j+1)));
            bj=trapz(t(Ij),pdf(Ij));
            mj(j,k)=M*bj;
        end;
        pdfdev(k)=norm(Hemp(:,k)-mj(:,k));
    end;
    [sortPdf,Ipdf]=sort(pdfdev);
    Ipdf=flipud(Ipdf);
    if(strcmp(sortopt,'pdf'))
        varargout{1}=Ipdf;
    end;
    if(strcmp(sortopt,'all'))
        varargout{3}=Ipdf;
    end;
end;

if((strcmp(plotopt,'plot')&&(strcmp(sortopt,'pdf')))||(strcmp(sortopt,'all')))
    
    Wpdf=W(:,Ipdf);
    figure;
    plotXmatrix(Wpdf(:,1:Ndisp),'align');
    title(sprintf('%i of %i Waveforms Sorted by Randomness',Ndisp,N));
    ylabel(sprintf('%i Waveforms most dissimilar from Gaussian Distribution',Ndisp));
    xlabel('Sample Number');
    axis tight;

end;