boxRejection

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  function rfiFlags = boxRejection(d, nSigma, nSize, upperLimit, lowerLimit)
  
 this function does a boxcar rejection -- calculates a running average and
 rejects data if the mean is above a certain level and the rms as well.
 
  IMPORTANT:  make sure the data you pass in does not have noise source
  observations mixed in with astronomical observations.  This will cause
  more things than you want to be flagged.

  typical timescales for RFI are of 15 seconds, so that's what we're using
  as our default time unit

  d        == data structure
  nSigma   == number of Sigma after which you throw out the data
  nTime    == time value to check on in seconds
  upperLimit == backend value above which we don't believe
  lowerLimit == backend value below which we don't believe

  defaults (5, 15, 8000, -8000) == (nSigma, nTime, upperLimit, lowerLimit);

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

0001 function rfiFlags = boxRejection(d,  nSigma, nTime, upperLimit, lowerLimit)
0002 
0003 
0004 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0005 %  function rfiFlags = boxRejection(d, nSigma, nSize, upperLimit, lowerLimit)
0006 %
0007 % this function does a boxcar rejection -- calculates a running average and
0008 % rejects data if the mean is above a certain level and the rms as well.
0009 %
0010 %  IMPORTANT:  make sure the data you pass in does not have noise source
0011 %  observations mixed in with astronomical observations.  This will cause
0012 %  more things than you want to be flagged.
0013 %
0014 %  typical timescales for RFI are of 15 seconds, so that's what we're using
0015 %  as our default time unit
0016 %
0017 %  d        == data structure
0018 %  nSigma   == number of Sigma after which you throw out the data
0019 %  nTime    == time value to check on in seconds
0020 %  upperLimit == backend value above which we don't believe
0021 %  lowerLimit == backend value below which we don't believe
0022 %
0023 %  defaults (5, 15, 8000, -8000) == (nSigma, nTime, upperLimit, lowerLimit);
0024 %
0025 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0026 
0027 switch(nargin)
0028   
0029   case 1
0030     nSigma = 5;
0031     nTime  = 15;
0032     upperLimit = 8000;
0033     lowerLimit = -8000;
0034 
0035   case 2
0036     nTime  = 15;
0037     upperLimit = 8000;
0038     lowerLimit = -8000;    
0039     
0040   case 3
0041     upperLimit = 8000;
0042     lowerLimit = -8000;    
0043     
0044   case 4    
0045     lowerLimit = -8000;    
0046     
0047   case 5
0048     
0049   otherwise
0050     error('Wrong number of inputs');
0051 end
0052 
0053 % we do have to switch for all allowable bits.  especially if we have a
0054 % noise source turned on.
0055 rfiFlags = zeros(size(d.antenna0.receiver.data,1),1);
0056 
0057 
0058 % first we check the limits, which are easiest
0059 upLimFlags = sum(d.antenna0.receiver.data > upperLimit,2)>0;
0060 loLimFlags = sum(d.antenna0.receiver.data < lowerLimit,2)>0;
0061 
0062 % convert the time value into a number of points
0063 nSize = round(nTime./(mean(diff(d.antenna0.receiver.utc))*24*60*60));
0064 % quick flag on rfi
0065 [flags] = boxCarReject(d.antenna0.receiver.data, nSigma, nSize, upperLimit);
0066 rfiFlags = flags>0 | upLimFlags | loLimFlags;
0067 
0068 return;
0069 
0070 
0071 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0072 %  function [thisData flags] = boxCarReject(d.antenna0.receiver.data, nSigma, nSize, upperLimit);
0073 %
0074 %  this function loops through your bin sizes and just sets things up
0075 %  for the following function.
0076 %
0077 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0078 
0079 
0080 function [fullFlagVec] = boxCarReject(data,nSigma, nSize, upperLimit)
0081 
0082 % check the size of nSize
0083 if(nSize > size(data,1))
0084   nSize = size(data,1);
0085 end
0086 numBox = floor(size(data,1)/nSize);
0087 
0088 finalFlags = [];
0089 finalData = [];
0090 soFar = 0;
0091 meanVals = [];
0092 rmsVals  = [];
0093 indexVec = [];
0094 for m=1:numBox
0095   if(m==numBox)
0096     thisBin = data((m-1)*nSize+1:size(data,1),:,:);
0097     indexVec((m-1)*nSize+1:size(data,1)) = m;
0098   else
0099     thisBin = data((m-1)*nSize+1:m*nSize,:,:);
0100     indexVec((m-1)*nSize+1:m*nSize) = m;  
0101   end
0102 
0103   meanVals(m,:) = nanmean(thisBin);
0104   rmsVals(m,:)  = sqrt(nanvar(thisBin));
0105 end  
0106 indexVec = indexVec';
0107 
0108 % the zeroth order thing is to just flag the data based on it being super
0109 % noisy, i.e., the rmsValue being greater than 1.1x the median value.
0110 medVals = nanmedian(rmsVals);
0111 medVals = repmat(medVals, [size(meanVals,1) 1])*1.1;
0112 
0113 flagged = rmsVals>medVals;
0114 
0115 % now we convert this back to a vector that matches the data.
0116 flagVec = flagged(:,1) + flagged(:,6) > 0;
0117 f = find(flagVec);
0118 fullFlagVec = zeros(size(indexVec));
0119 for m=1:length(f)
0120   fullFlagVec(indexVec==f(m)) = 1;
0121 end
0122 
0123 return;
0124 
0125 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0126 %  function [data saveFlags] = outlier(thisBin, nSigma)
0127 %
0128 %  this is the function that goes through and picks your outliers as
0129 %  determined by the nSigma field.
0130 %
0131 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0132 
0133 function [data saveFlags] = outlier(thisBin, nSigma, upperLimit)
0134 
0135 nSize = size(thisBin,1);
0136 % d is an array
0137 f=1;
0138 saveFlags = zeros(size(thisBin));
0139 if(size(thisBin,1) == 1)
0140   flags = thisBin>upperLimit;
0141   data = thisBin;
0142   return;
0143 end
0144 if(any(nanmean(thisBin,1) > upperLimit))
0145     display(['Mean value greater than upper limit']);
0146     display(['Readjusting upper limit            ']);
0147     upperLimit = max(nanmean(thisBin,1)) + 5;
0148 end
0149 while(~isempty(f))
0150   thisMean = repmat(nanmean(thisBin,1), [nSize 1 1]);
0151   thisStd  = repmat(nanstd(thisBin,1), [nSize 1 1]);
0152   flags = (thisBin > thisMean + nSigma*thisStd | thisBin < thisMean - ...
0153       nSigma*thisStd | thisBin > upperLimit);
0154   f = find(flags);
0155   saveFlags = saveFlags + flags;
0156   thisBin = ~flags.*thisBin + thisMean.*flags;
0157   saveFlags = saveFlags>0;
0158 end
0159 
0160 data = thisBin;
0161 
0162 return;