[w]=signalwidth(x,varargin)

function [w]=signalwidth(x,varargin)
% [w]=signalwidth(x,varargin)
% Estimates the signal width of data input, using definition:
%   
%   rho = cumtrapz(t,x.^2)/trapz(t,x.^2)
% 
% The integration limits go from mean(t)-w/2 and mean(t)+w/2; w is defined
% as the signal width, and gives the time at which rho % of the signal is
% captured.  It is an informal measure of the localization of a signal.
% Rho, in most applications, is ~90%.
%--------------------------------------------------------------------------
% INPUT
% x: a vector or matrix of data.  If a matrix, data are taken column wise.
% rho: (optional) gives the % of the signal defining the signal width.
%      Default is 90%.
% t: (optional) a time-vector, for computing the signal width of time
%     series.
% plotopt: (optional) enter string 'plot' to see the waveforms with signal
%          widths plotted on top of them.
%
% OUTPUT
% w: the width of the signal(s) given as mean(t)-w/2 and mean(t)+w/2.
%
%Example: 
% %Find t-values giving 75% of the signal energy (default is 90%)
% rho=0.85
% t=[0.01:0.01:5]';
% g=zeros(length(t),20);
% for k=1:20
%   g(:,k)=gausswin(length(t),k+1.5);
% end;
% w=signalwidth(g,t,rho,'plot');
%--------------------------------------------------------------------------

%Default input variables
plotopt='none';
rho=0.9;
t=[1:length(x)]';

%Loop over variable inputs and assign values
if(nargin>1)
    for k=1:length(varargin);
        if(isnumeric(varargin{k}))
            if(length(varargin{k})==1)
                rho=varargin{k};
            else
                t=varargin{k};
            end;
        else
            plotopt=varargin{k};
        end;
    end;
end;

%Get time-series indices at from first sample to mid-sample and time
%series indices from mid-sample to end sample.
Il=fliplr(1:floor(length(x)/2));
Ir=ceil(length(x)/2):length(x);
Ir=Ir(1:length(Il));
indefx=abs(cumtrapz(t(Il),x(Il,:).^2))+abs(cumtrapz(t(Ir),x(Ir,:).^2));
defx=trapz(t,x.^2);

w=zeros(size(x,2),2);
for k=1:size(x,2)
    I=find(indefx(:,k)./defx(k)>=rho,1);
    w(k,1)=t(Ir(1)-I);
    w(k,2)=t(I+Ir(1));
end;

if(strcmp(plotopt,'plot'));
    plotXmatrix(maxalign(x),t); hold on;
    for k=1:size(x,2)
        plot(w(k,1),k-1,'c+','linewidth',2); 
        plot(w(k,2),k-1,'c+','linewidth',2);
    end;
    s=title(sprintf('%i Input Signals and %i Percent Signal Width Times',size(x,2),floor(100*rho)));
    set(s,'fontsize',14); set(gca,'fontsize',14);
end;