function [cleanvar,Qa]=fix_vbrange(x, settings)
%fix_vbrange attempts to clean up adv and pcadp stats data with spikes
% a test program to evaluate a method of cleaning spiky data-
% this is a two step process, where step 1 is a simple min-max
% excluder, and step 2 uses std to find spikes.  It replaces only those
% points with data from medfilt over 15 points in first part, over 7 points
% in the second.
%
% usage: [cl_dat,Qa]=fix_vbrange(x,settings);
%   inputs are the data and the settings
%   settings - the fields below are needed:
%     settings.min=550;settings.max=700;
%     settings.npts=50, settings.std_threshold=20; settings.nstds=2;
%  outputs are the clean data, and the points where something was replaced

% check for arguments
if nargin < 2; help mfilename; return; end

% set up defaults for settings, if not supplied
if exist('settings', 'var'),
    if ~isfield(settings, 'min'), settings.min=400; end
    if ~isfield(settings, 'max'), settings.max=500; end
    if ~isfield(settings,'nstds'), settings.nsds=2.0; end
    if ~isfield(settings, 'npts'), settings.npts=20; end
    if ~isfield(settings, 'std_threshold'), settings.std_threshold=10; end
end
    Qa.nbad = 0;
    Qa.nfixed_mm = 0;
    Qa.nfixed_filt = 0;
    neg_flg=0;      %flag indicating the presence of negative data


% find spiky data, using min and max values to defing the spikes
loc=find(x< settings.min | x > settings.max);

% now do a median filter over a longish number of points
od_mf=medfilt(x,15);

%finally replace only the spiky locations with the filtered points
clvar=x;
clvar(loc)=od_mf(loc);
cleanvar=clvar;
Qa.nfixed_mm=length(loc);

% the above is all that's needed if the vrange or brange was relatively flat
% but if the  data slopes or has biggish jumps, it won't do it all

% use std as a threshold for doing more checking
if gstd(cleanvar) > settings.std_threshold
    % subset the data and find spikes by chunk
    npcs=length(cleanvar)/settings.npts;
    indx=1; spikes=[];
    for ik=1:npcs
        dat=cleanvar(indx:(indx+settings.npts)-1);
        chunkstd=gstd(dat);
        chunkmed=gmedian(dat);
        if chunkmed < 0
            nbad=find(dat <0)
        else
           nbad=find(abs(dat-chunkmed) > settings.nstds*chunkstd);
        end
        spikes=[spikes; nbad+(indx-1)];
        indx=indx+settings.npts;
        clear nbad chunkstd chunkmed dat
    end
    dat=cleanvar(indx:end);
    chunkstd=gstd(dat);
    chunkmed=gmedian(dat);
    if chunkmed < 10
            nbad=find(dat <=0);
            neg_flg=1;
    else
       nbad=find(abs(dat-chunkmed) > settings.nstds*chunkstd);
    end
    spikes=[spikes; nbad+(indx-1)];
    clear nbad chunkstd chunkmed dat

    % then do a median filter of the whole thing with 7 points
    if neg_flg==1       %negative data is present
        cleanvar(spikes)=ones(size(spikes)).*1e35;
    else
        vr_mf=medfilt(cleanvar,7);
        fvr=cleanvar;
        % substitute the median filtered data in place of the spikes
        fvr(spikes)=vr_mf(spikes);
        cleanvar=fvr;
    end
    
    % do one last check for negative data
      negv=find(cleanvar <=0);
      if ~isempty(negv)
          % when autonNaN is ON during processing, putting Nan in results
          % in 1e35 in data for real (putting in 1e35 is OK too)
       cleanvar(negv)=ones(size(negv)).*1e35;
       spikes=sort([spikes; negv]);
      end
    Qa.nfixed_filt = length(spikes);
    loc=[loc; spikes];
end
Qa.nbad=length(loc);