boxCarRejection

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  function [fullFlagVec] = boxCarRejection(data,nSigma, nSize, upperLimit)

  this function loops through your bin sizes and just sets things up
  for the following function.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0001 function [fullFlagVec] = boxCarRejection(data,nSigma, nSize, upperLimit)
0002 
0003 
0004 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0005 %  function [fullFlagVec] = boxCarRejection(data,nSigma, nSize, upperLimit)
0006 %
0007 %  this function loops through your bin sizes and just sets things up
0008 %  for the following function.
0009 %
0010 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0011 
0012 
0013 % check the size of nSize
0014 if(nSize > size(data,1))
0015   nSize = size(data,1);
0016 end
0017 numBox = floor(size(data,1)/nSize);
0018 
0019 finalFlags = [];
0020 finalData = [];
0021 soFar = 0;
0022 meanVals = [];
0023 rmsVals  = [];
0024 indexVec = [];
0025 for m=1:numBox
0026   if(m==numBox)
0027     thisBin = data((m-1)*nSize+1:size(data,1),:,:);
0028     indexVec((m-1)*nSize+1:size(data,1)) = m;
0029   else
0030     thisBin = data((m-1)*nSize+1:m*nSize,:,:);
0031     indexVec((m-1)*nSize+1:m*nSize) = m;  
0032   end
0033 
0034   meanVals(m,:) = nanmean(thisBin);
0035   rmsVals(m,:)  = sqrt(nanvar(thisBin));
0036 end  
0037 indexVec = indexVec';
0038 
0039 % the zeroth order thing is to just flag the data based on it being super
0040 % noisy, i.e., the rmsValue being greater than 1.1x the median value.
0041 medVals = nanmedian(rmsVals);
0042 medVals = repmat(medVals, [size(meanVals,1) 1])*1.1;
0043 
0044 flagged = rmsVals>medVals;
0045 
0046 % now we convert this back to a vector that matches the data.
0047 flagVec = flagged(:,1);
0048 f = find(flagVec);
0049 fullFlagVec = zeros(size(indexVec));
0050 for m=1:length(f)
0051   fullFlagVec(indexVec==f(m)) = 1;
0052 end
0053 
0054 return;
0055 
0056 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0057 %  function [data saveFlags] = outlier(thisBin, nSigma)
0058 %
0059 %  this is the function that goes through and picks your outliers as
0060 %  determined by the nSigma field.
0061 %
0062 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0063 
0064 function [data saveFlags] = outlier(thisBin, nSigma, upperLimit)
0065 
0066 nSize = size(thisBin,1);
0067 % d is an array
0068 f=1;
0069 saveFlags = zeros(size(thisBin));
0070 if(size(thisBin,1) == 1)
0071   flags = thisBin>upperLimit;
0072   data = thisBin;
0073   return;
0074 end
0075 if(any(nanmean(thisBin,1) > upperLimit))
0076     display(['Mean value greater than upper limit']);
0077     display(['Readjusting upper limit            ']);
0078     upperLimit = max(nanmean(thisBin,1)) + 5;
0079 end
0080 while(~isempty(f))
0081   thisMean = repmat(nanmean(thisBin,1), [nSize 1 1]);
0082   thisStd  = repmat(nanstd(thisBin,1), [nSize 1 1]);
0083   flags = (thisBin > thisMean + nSigma*thisStd | thisBin < thisMean - ...
0084       nSigma*thisStd | thisBin > upperLimit);
0085   f = find(flags);
0086   saveFlags = saveFlags + flags;
0087   thisBin = ~flags.*thisBin + thisMean.*flags;
0088   saveFlags = saveFlags>0;
0089 end
0090 
0091 data = thisBin;
0092 
0093 return;