kurtclip

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

  function lowRFI = kurtclip(data,binSize, sig)

   Returns an array of same length as data.  Elements of 
   returned array are 0 when RFI is present, 1 otherwise.
   Kurtosis is calculated in bins of length binSize elements
   (100 works well).  RFI is detected if the kurtosis is 
   greater than sig sigma above 3.  If sig is not provided, 
   a default of 5 is used.

   11-May-2010 (MAS)

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

0001 function lowRFI = kurtclip(data,binSize, sig)
0002 
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 %  function lowRFI = kurtclip(data,binSize, sig)
0006 %
0007 %   Returns an array of same length as data.  Elements of
0008 %   returned array are 0 when RFI is present, 1 otherwise.
0009 %   Kurtosis is calculated in bins of length binSize elements
0010 %   (100 works well).  RFI is detected if the kurtosis is
0011 %   greater than sig sigma above 3.  If sig is not provided,
0012 %   a default of 5 is used.
0013 %
0014 %   11-May-2010 (MAS)
0015 %
0016 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0017 
0018 if(nargin<3)
0019      sig = 5;
0020 end
0021 
0022 
0023 % Calculate the actual kurtosis threshold.
0024 thresh = sig * sqrt(24 / binSize) + 3;
0025 
0026 
0027 
0028 % Find the size of the data array:
0029 dataSize = size(data,1);
0030 
0031 % How many time bins does this give us?  Round down.
0032 timeSize = floor(dataSize / binSize);
0033 
0034 
0035 % Clip the data so that we can make equal-sized bins.
0036 clData = data(1:binSize*timeSize);
0037 
0038 % Reshape the array to make the bins.
0039 reData = reshape(clData,binSize,timeSize);
0040 
0041 % Calculate the kurtosis of the bins.
0042 kurtData = kurtosis(reData);
0043 
0044 
0045 % We now want to interpolate the kurtosis, so we need two x-arrays.  One
0046 % for the caluclated kurtosis, and one for the interpolated array.
0047 x = binSize*(1:timeSize) - binSize/2;
0048 xi = 1:dataSize;
0049 
0050 % Do the actual interpolation.  Any elements that were clipped are set to
0051 % 1.
0052 lowRFI = logical(interp1(x,kurtData < thresh,xi,'nearest',1) > 0);
0053 
0054 % And that's it!
0055 return;
0056