function pv=sigratio(x,win,outwin,outalpha,plfg,plstr)
%
%		function pv=sigratio(x,win,outwin,outalpha,plfg,plstr)
%
%  x = data vector
%  win = half of sigma comparison window
%  outwin = half width of outlier rejection window
%  outalpha = outlier false alarm prob
%  plfg is plotflag, >0 to plot
%  plstr = data title string
%  pv is length(x) x 3  pv(:,1) is for sigsq2 > sigsq1, pv(:,2) is sigsq2 < sigsq1
%  and pv(:,3) is for sigsq2 ~= sigsq1, two sided


out=tout(x,outwin,outalpha);
outind=find(out>0);	% indexes where out = 1
nout=length(outind);
nx=length(x);
pv=ones(nx,3);
if 2*win+1 > nx
   error('not enough data in vector x');
end
nk=nx-2*win+1; % nk=1 if nx=2*win;
for k=1:nk
   k1=k+win; 
   start1=k1-win;stop1=start1+win-1;	% 
   start2=k1;stop2=start2+win-1;			% second window starts at k1
   index1=setdiff([start1:stop1],outind); % excludes ones where out=1
   n1=length(index1);
   index2=setdiff([start2:stop2],outind);
   n2=length(index2);
   sigsq1=cov(x(index1));
   sigsq2=cov(x(index2));
   fval=sigsq2/sigsq1;
   pv(k1,1)=1-fcdf(fval,n2,n1); % note sigsq2 on top and testing for sigsq2 > sigsq1
   pv(k1,2)=1-fcdf(1/fval,n1,n2); % test sigsq1 > sigsq2
   pv(k1,3)=min(min(pv(k1,1),pv(k1,2))*2,1); % pv for not equal test
end
if plfg
   figure;
   subplot(211);
   plot(x);
   title(['original data : ',plstr]);
   subplot(212);
   plot(x.*(1-out));
   title([int2str(nout),' outliers nulled, alpha=',num2str(outalpha),' outwin=',int2str(outwin)]);
   
   figure;
   subplot(211);
   plot(x.*(1-out));
   title([plstr,'  ',int2str(nout),' outliers nulled, alpha=',num2str(outalpha),' outwin=',int2str(outwin)]);
   subplot(212);
   semilogy(pv(:,1));
   title(['P-value for \sigma_2 > \sigma_1 test, win:',int2str(win)]);
   
   figure;
   subplot(211);
   plot(x.*(1-out));
   title([plstr,'  ',int2str(nout),' outliers nulled, alpha=',num2str(outalpha),' outwin=',int2str(outwin)]);
   subplot(212);
   semilogy(pv(:,2));
   title(['P-value for \sigma_2 < \sigma_1 test, win:',int2str(win)]);
   
   figure;
   subplot(211);
   plot(x.*(1-out));
   title([plstr,'  ',int2str(nout),' outliers nulled, alpha=',num2str(outalpha),' outwin=',int2str(outwin)]);

     subplot(212);
   semilogy(pv(:,3));
   title(['P-value for \sigma_2 \neq \sigma_1 test, win:',int2str(win)]);
end