flagRFI_std

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  function  [flagShort flagLong] = flagRFI_std(d,preFlag,params)

   Dec 20, 2010 (MAS) - Time series RFI flagging for C-BASS data.  This is
   a two step process.  First, data are flagged based upon sub-second RMS 
   values from the polarization channels.  Next, flagging is performed on
   10s RMS values from the polarization channels.  These flagging steps
   largely overlap, but do catch some different types.

   March 2, 2011 (MAS) - Converted to eval version;  allows evaluation of
   flag values.

   March 11, 2011 (MAS) - Ignores those regions flagged by position.

   March 23, 2011 (MAS) - Trial final version; using strategy tested on
   late October 2010 data.

   April 14, 2011 (MAS) - Trial final version; comparison between I1 and
   I2 on long scales.

   April 12, 2012 (MAS) - Updated to work with new pipeline order.
   Beginning some simplifying modifications.
   Renamed variables so that their names aren't nonsense.
   Removed the last parameter because it was no good.

   params(1): minimum shortflag 1
   params(2): minimum shortflag 2
   params(3): shortflag slope
   params(4): minimum shortflag (alt)
   params(5): minimum longflag 1 A
   params(6): minimum longflag 1 B
   params(7): minimum longflag 2
   params(8): minimum I longflag
   params(9): minimum I difference.
   params(10): final smoothing (seconds)

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

0001 function flag = flagRFI_std(d,preFlag,params)
0002 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0003 %  function  [flagShort flagLong] = flagRFI_std(d,preFlag,params)
0004 %
0005 %   Dec 20, 2010 (MAS) - Time series RFI flagging for C-BASS data.  This is
0006 %   a two step process.  First, data are flagged based upon sub-second RMS
0007 %   values from the polarization channels.  Next, flagging is performed on
0008 %   10s RMS values from the polarization channels.  These flagging steps
0009 %   largely overlap, but do catch some different types.
0010 %
0011 %   March 2, 2011 (MAS) - Converted to eval version;  allows evaluation of
0012 %   flag values.
0013 %
0014 %   March 11, 2011 (MAS) - Ignores those regions flagged by position.
0015 %
0016 %   March 23, 2011 (MAS) - Trial final version; using strategy tested on
0017 %   late October 2010 data.
0018 %
0019 %   April 14, 2011 (MAS) - Trial final version; comparison between I1 and
0020 %   I2 on long scales.
0021 %
0022 %   April 12, 2012 (MAS) - Updated to work with new pipeline order.
0023 %   Beginning some simplifying modifications.
0024 %   Renamed variables so that their names aren't nonsense.
0025 %   Removed the last parameter because it was no good.
0026 %
0027 %   params(1): minimum shortflag 1
0028 %   params(2): minimum shortflag 2
0029 %   params(3): shortflag slope
0030 %   params(4): minimum shortflag (alt)
0031 %   params(5): minimum longflag 1 A
0032 %   params(6): minimum longflag 1 B
0033 %   params(7): minimum longflag 2
0034 %   params(8): minimum I longflag
0035 %   params(9): minimum I difference.
0036 %   params(10): final smoothing (seconds)
0037 %
0038 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0039 
0040 
0041 %%%%%%%%%%%%%%
0042 % Parameters:
0043 %%%%%%%%%%%%%%
0044 
0045 % Number of iterations when removing outliers from vectors.
0046 nClip = 5;
0047 
0048 % Two lengthscales used for the short flagging stage.  In units of 10ms
0049 % samples.
0050 shortLength_A = 25;
0051 shortLength_B = 2 * shortLength_A;
0052 
0053 % Smoothing length.  Units of shortLength.  ie, 500 corresponds to 250s.
0054 subStdLength = 250; % 100
0055 
0056 % Two lengthscales used for the long flagging stage.  In units of 10ms
0057 % samples.
0058 longLength_A = 250;
0059 longLength_B = 4 * longLength_A;
0060 
0061 
0062 
0063 
0064 %%%%%%%%%%%%%%%%%%%%%%%%%%
0065 % Preliminary bits...
0066 %%%%%%%%%%%%%%%%%%%%%%%%%%
0067 
0068 % How long is our data time series...
0069 dLength = length(d.antenna0.receiver.data);
0070 
0071 
0072 % Calculate the lengths of the short flagging standard deviation vectors.
0073 nShort_B = floor(dLength / shortLength_B);
0074 nShort_A = nShort_B * shortLength_B / shortLength_A;
0075 
0076 % Calculate the lengths of the short flagging standard deviation vectors.
0077 nLong_B = floor(dLength / longLength_B);
0078 nLong_A = nLong_B * longLength_B / longLength_A;
0079 
0080 
0081 % This will only proceed with a subvector of the data, determined by
0082 % fitting in an integral number of shortLength_B.
0083 shortStart = floor((dLength - nShort_B * shortLength_B) / 2) + 1;
0084 shortEnd = shortStart + nShort_B * shortLength_B - 1;
0085 
0086 longStart = floor((dLength - nLong_B * longLength_B) / 2) + 1;
0087 longEnd = longStart + nLong_B * longLength_B - 1;
0088 
0089 
0090 
0091 
0092 %%%%%%%%%%%%%%%%%%%%%%%%%%
0093 % Smoothing parameters...
0094 %%%%%%%%%%%%%%%%%%%%%%%%%%
0095 
0096 % There's going to be some smoothing done on the standard deviation vector.
0097 % We should prepare by calculating some length and padding parameters.
0098 sLength_A = 2 ^ ceil(log2(nShort_A + 2 * 10 * subStdLength));
0099 sBuffStart_A = floor((sLength_A - nShort_A) / 2);
0100 sBuffEnd_A = ceil((sLength_A - nShort_A) / 2);
0101 
0102 sLength_B = 2 ^ ceil(log2(nShort_B + 2 * 10 * subStdLength));
0103 sBuffStart_B = floor((sLength_B - nShort_B) / 2);
0104 sBuffEnd_B = ceil((sLength_B - nShort_B) / 2);
0105 
0106 
0107 % Also, let's make the kernel for the convolution.
0108 xKernel_A = (1:sLength_A)' - 1;
0109 gaussKernelStd_A = ...
0110     (exp(-0.5*(xKernel_A / subStdLength).^2) + ...
0111     exp(-0.5*((sLength_A - xKernel_A) / subStdLength).^2)) / ...
0112     (subStdLength*sqrt(2*pi));
0113 
0114 xKernel_B = (1:sLength_B)' - 1;
0115 gaussKernelStd_B = ...
0116     (exp(-0.5*(xKernel_B / subStdLength).^2) + ...
0117     exp(-0.5*((sLength_B - xKernel_B) / subStdLength).^2)) / ...
0118     (subStdLength*sqrt(2*pi));
0119 
0120 % Also, an X array which will later be used for interpolation purposes.
0121 xStd_A = (1:nShort_A)';
0122 xStd_B = (1:nShort_B)';
0123 
0124 
0125 
0126 %%%%%%%%%%%%%%%%%%%%%%%%%%
0127 % Some more prep before looping...
0128 %%%%%%%%%%%%%%%%%%%%%%%%%%
0129 
0130 % Make the noise diode events look like non-noise diode events.
0131 d = noiseRemove(d);
0132 
0133 
0134 % Prep the positional flags for use with short flagging.
0135 preFlagShort_A = ...
0136     max(reshape(preFlag(shortStart:shortEnd),shortLength_A,nShort_A))';
0137 preFlagShort_B = ...
0138     max(reshape(preFlag(shortStart:shortEnd),shortLength_B,nShort_B))';
0139 
0140 
0141 % Prep the loop's output arrays.
0142 shortRMS_A = zeros(nShort_A,4);
0143 shortRMS_B = zeros(nShort_B,4);
0144 
0145 longRMS_A = zeros(nLong_A,4);
0146 longRMS_B = zeros(nLong_B,4);
0147 
0148 
0149 
0150 
0151 %%%%%%%%%%%%%%%%%%%%%%%%%%
0152 % Loop over the polarizations,
0153 % calculating the normalized
0154 % RMS over the four timescales.
0155 %%%%%%%%%%%%%%%%%%%%%%%%%%
0156 
0157 % Polarization iterator.
0158 polNum = 0;
0159 for k=[1 6 7 8]
0160     
0161     % Iterate!
0162     polNum = polNum + 1;
0163 
0164     
0165 
0166     %-----
0167     % So here's the actual short flagging part.  Done in easy steps.
0168     %-----
0169     
0170     
0171     % Step 0:  Calculate the raw RMS values.
0172     %-----
0173 
0174     % Calculate the standard deviations!
0175     shortRMS_A(:,polNum) = std( ...
0176         reshape(d.antenna0.receiver.data(shortStart:shortEnd,k), ...
0177         shortLength_A, nShort_A));
0178     shortRMS_B(:,polNum) = std( ...
0179         reshape(d.antenna0.receiver.data(shortStart:shortEnd,k), ...
0180         shortLength_B, nShort_B));
0181 
0182     % In case we missed some data points at the beginning and end.
0183     if shortEnd ~= dLength
0184         shortRMS_A(1,polNum) = ...
0185             std(d.antenna0.receiver.data(1:shortStart+shortLength_A-1,k));
0186         shortRMS_A(end,polNum) = ...
0187             std(d.antenna0.receiver.data(shortEnd-shortLength_A+1:end,k));    
0188 
0189         shortRMS_B(1,polNum) = ...
0190             std(d.antenna0.receiver.data(1:shortStart+shortLength_B-1,k));
0191         shortRMS_B(end,polNum) = ...
0192             std(d.antenna0.receiver.data(shortEnd-shortLength_B+1:end,k));    
0193     end    
0194         
0195     
0196     % Step 1:  Clip immediate outliers using both arrays.
0197     %-----
0198     shortClip_A = stdClip(shortRMS_A(:,polNum),3,nClip,1);
0199     shortClip_B = stdClip(shortRMS_B(:,polNum),3,nClip,1);
0200     
0201             
0202     % Step 2:  Ignoring outliers, calculate running median.
0203     %-----
0204     
0205     % Fold the positional flagging into the flag vector for this part:
0206     shortClip_A = smooth((shortClip_A == 1) | (preFlagShort_A == 1),40) > 0;
0207     shortClip_B = smooth((shortClip_B == 1) | (preFlagShort_B == 1),20) > 0;
0208 
0209     % Get the actual running median
0210     runningMed_A = padSmooth(xStd_A, shortRMS_A(:,polNum), ...
0211         gaussKernelStd_A, shortClip_A, ...
0212         [sLength_A sBuffStart_A sBuffEnd_A]);
0213     runningMed_B = padSmooth(xStd_B, shortRMS_B(:,polNum), ...
0214         gaussKernelStd_B, shortClip_B, ...
0215         [sLength_B sBuffStart_B sBuffEnd_B]);
0216        
0217             
0218     % Step 3:  Subtract running median.
0219     %-----
0220     shortRMS_A(:,polNum) = shortRMS_A(:,polNum) - runningMed_A;
0221     shortRMS_B(:,polNum) = shortRMS_B(:,polNum) - runningMed_B;
0222 
0223     
0224     % Step 4:  Clip current outliers.
0225     %-----
0226     shortClip_A = stdClip(shortRMS_A(:,polNum),5,nClip,0);
0227     shortClip_B = stdClip(shortRMS_B(:,polNum),5,nClip,0);
0228 
0229     
0230     % Step 5:  Normalize the standard deviation vectors.
0231     %-----
0232     shortRMS_A(:,polNum) = ...
0233         shortRMS_A(:,polNum) / std(shortRMS_A(shortClip_A == 0,polNum));
0234     shortRMS_B(:,polNum) = ...
0235         shortRMS_B(:,polNum) / std(shortRMS_B(shortClip_B == 0,polNum));
0236 
0237 
0238     
0239     %-----
0240     % Now we calculate vectors for the long flagging part.  Easy steps.
0241     %-----
0242 
0243     % Step 1:  Interpolate the running median some more.
0244     %-----
0245     runningMedLong = ...
0246         interp1(((1:length(runningMed_A))' - 0.5) * shortLength_A, ...
0247         runningMed_A,(1:dLength),'pchip');
0248 
0249     
0250     % Step 2:  Calculate and normalize the standard deviation vectors.
0251     %-----
0252     longRMS_A(:,polNum) = ...
0253         std(detrend(reshape(d.antenna0.receiver.data(longStart:longEnd,k) ./ ...
0254         runningMedLong(longStart:longEnd)', ...
0255         longLength_A, nLong_A)));
0256     
0257     longRMS_B(:,polNum) = ...
0258         std(detrend(reshape(d.antenna0.receiver.data(longStart:longEnd,k) ./ ...
0259         runningMedLong(longStart:longEnd)', ...
0260         longLength_B, nLong_B)));
0261 
0262 
0263     
0264     % Step 3:  Clip current outliers.
0265     %-----
0266     flagStdVec2_A = stdClip(longRMS_A(:,polNum),1,nClip,1);
0267     flagStdVec2_B = stdClip(longRMS_B(:,polNum),1,nClip,1);
0268     
0269 
0270     % Step 4:  Subtract the mean.
0271     %-----
0272     longRMS_A(:,polNum) = longRMS_A(:,polNum) - ...
0273         mean(longRMS_A(flagStdVec2_A == 0,polNum));
0274     longRMS_B(:,polNum) = longRMS_B(:,polNum) - ...
0275         mean(longRMS_B(flagStdVec2_B == 0,polNum));
0276     
0277     
0278     % In case we missed some data points at the beginning and end.
0279     if longEnd ~= dLength
0280         longRMS_A(1,polNum) = ...
0281             std(detrend(d.antenna0.receiver.data(1:longStart+longLength_A-1,k) ./ ...
0282             runningMedLong(1:longStart+longLength_A-1)')) - 1;
0283         longRMS_A(end,polNum) = ...
0284             std(detrend(d.antenna0.receiver.data(longEnd-longLength_A+1:end,k) ./ ...
0285             runningMedLong(longEnd-longLength_A+1:end)')) - 1;    
0286 
0287         longRMS_B(1,polNum) = ...
0288             std(detrend(d.antenna0.receiver.data(1:longStart+longLength_B-1,k) ./ ...
0289             runningMedLong(1:longStart+longLength_B-1)')) - 1;
0290         longRMS_B(end,polNum) = ...
0291             std(detrend(d.antenna0.receiver.data(longEnd-longLength_B+1:end,k) ./ ...
0292             runningMedLong(longEnd-longLength_B+1:end)')) - 1;    
0293     end    
0294 end
0295 
0296 
0297 % Finish off the short flagging.
0298 %%%%%%%%%%%%%%%%
0299 
0300 % Combine the polarizations.
0301 shortI_A = mean(shortRMS_A(:,[1 4]),2);
0302 shortPol_A = max(shortRMS_A(:,2), shortRMS_A(:,3));
0303 
0304 % Do the same for length B.  Need to stretch it, though, to match A.
0305 shortI_B = reshape(repmat(mean(shortRMS_B(:,[1 4]),2)', ...
0306     shortLength_B/shortLength_A,1), nShort_A, 1);
0307 shortPol_B = reshape(repmat(max(shortRMS_B(:,2), ...
0308     shortRMS_B(:,3))',shortLength_B/shortLength_A,1), nShort_A, 1);
0309 
0310 
0311 % Do the actual short-flagging comparison.
0312 shortFlag_1 = ...
0313     (shortPol_A > params(1)) & ...
0314     (shortPol_A > params(2) * shortI_A) & ...
0315     (shortPol_B > params(2) * shortI_B) & ...
0316     (shortPol_B < params(3) * shortPol_A);
0317 
0318 shortFlag_2 = shortPol_A > params(4);
0319     
0320 shortFlag = shortFlag_1 | shortFlag_2;
0321 
0322 
0323 % Stretch it out.
0324 shortFlagRep = zeros(dLength,1);
0325 shortFlagRep(shortStart:shortEnd) = ...
0326     reshape(repmat(shortFlag',shortLength_A,1),shortLength_A*nShort_A,1);
0327 
0328 if shortEnd ~= dLength
0329     shortFlagRep(1:shortStart) = shortFlag(1);
0330     shortFlagRep(shortEnd:end) = shortFlag(end);
0331 end    
0332 
0333 
0334 
0335 % Now we do the longer flagging.
0336 %%%%%%%%%%%%%%%%%
0337 
0338 % Combine the polarizations.
0339 %longI_A = mean(longRMS_A(:,[1 4]),2);
0340 longPol_A = max(longRMS_A(:,2), longRMS_A(:,3));
0341 
0342 % Do the same for length B.  Need to stretch it, though, to match A.
0343 longI_B = reshape(repmat(mean(longRMS_B(:,[1 4]),2)', ...
0344     longLength_B/longLength_A,1), nLong_A, 1);
0345 longI_B_max = reshape(repmat(max(longRMS_B(:,[1 4]),[],2)', ...
0346     longLength_B/longLength_A,1), nLong_A, 1);
0347 longI_B_min = reshape(repmat(min(longRMS_B(:,[1 4]),[],2)', ...
0348     longLength_B/longLength_A,1), nLong_A, 1);
0349 
0350 longPol_B = reshape(repmat(max(longRMS_B(:,2), ...
0351     longRMS_B(:,3))', longLength_B/longLength_A,1), nLong_A, 1);
0352 
0353 
0354 
0355 % Do the actual long-flagging comparison.
0356 longFlag_1 = ...
0357     (longPol_A > params(5)) & ...
0358     (longPol_B > params(6)) & ...
0359     (longPol_B > params(7) * longI_B);
0360 
0361 longFlag_2 = ...
0362     (longI_B_max > params(8)) & ...
0363     (longI_B_max > params(9) * longI_B_min);
0364 
0365 longFlag = longFlag_1 | longFlag_2;
0366 
0367 
0368 %plot(((1:length(longI_B)) - 0.5)*longLength_A,longFlag,'k');
0369 %hold off;
0370 
0371 % Stretch it out.
0372 longFlagRep = zeros(dLength,1);
0373 longFlagRep(longStart:longEnd) = ...
0374     reshape(repmat(longFlag',longLength_A,1),longLength_A*nLong_A,1);
0375 
0376 if longEnd ~= dLength
0377     longFlagRep(1:longStart) = longFlag(1);
0378     longFlagRep(longEnd:end) = longFlag(end);
0379 end    
0380 
0381 
0382 % Now we combine the short and long flags.
0383 % Don't flag pre-flagged regions.
0384 %%%%%%%%%%%%%%%%%
0385 
0386 % Here is the combination.
0387 flagCombine = (shortFlagRep | longFlagRep) & ~preFlag;
0388 
0389 % Make sure that flags aren't isolated.
0390 flagCombine = flagCombine & ...
0391     cat(1,ones(12,1),flagCombine(1:end-12)) & ...
0392     cat(1,flagCombine(13:end),ones(12,1));
0393 
0394 
0395 % Finally, flag a specified number of seconds before and after a given
0396 % flag.
0397 if params(10) > 0
0398     flag = (convn(double(flagCombine),ones(100*params(10),1),'same') > 0);
0399 else
0400     flag = (flagCombine > 0);
0401 end
0402 
0403 
0404 end
0405 
0406 
0407 % The sigma clippling subroutine.
0408 function clipVec = stdClip(dataVec,thresh,nClip,medSubtract)
0409 
0410 % Either we subtract the median or we don't.
0411 if medSubtract ~= 0
0412     medSubtract = 1;
0413 end
0414 
0415 % Initialize the clip vector.
0416 clipVec = zeros(length(dataVec),1);
0417 
0418 % Do the clipping...
0419 for l=1:nClip
0420     
0421     dataStdStd = std(dataVec(clipVec == 0));
0422     dataStdMed = median(dataVec(clipVec == 0));
0423     
0424     clipVec = abs(dataVec - medSubtract * dataStdMed) > ...
0425         thresh * dataStdStd;
0426     
0427 end
0428 
0429 end
0430 
0431 
0432 function smoothVec = padSmooth(xVec,dataVec,kernelVec,flagVec,padParams)
0433 
0434 
0435 % Interpolate over the flagged regions.
0436 dataVecT = interp1(xVec(flagVec == 0), ...
0437     smooth(dataVec(flagVec == 0),ceil(length(dataVec) / 100)), ...
0438     (1:length(dataVec))', 'pchip');
0439 
0440 dataVecF = dataVec .* (1 - flagVec) + flagVec .* dataVecT;
0441 
0442 if flagVec(1) == 1
0443     noFlagPos = find(flagVec == 0,1,'first');
0444     
0445     dataVecF(1:noFlagPos-1) = mean(dataVecF(noFlagPos:noFlagPos+99));
0446 end    
0447 if flagVec(end) == 1
0448     noFlagPos = find(flagVec == 0,1,'last');
0449     
0450     dataVecF(noFlagPos+1:end) = mean(dataVecF(noFlagPos-99:noFlagPos));
0451 end    
0452     
0453 
0454 
0455 % Pad the data up to a power of two.
0456 dataVecPad = zeros(padParams(1),1);
0457 dataVecPad(1:padParams(2)) = ...
0458     repmat(median(dataVecF(1:10)),padParams(2),1);
0459 dataVecPad(padParams(2)+1:padParams(2)+length(dataVec)) = dataVecF;
0460 dataVecPad(end-padParams(3)+1:end) = ...
0461     repmat(median(dataVecF(end-9:end)),padParams(3),1);
0462 
0463 % Perform a Fourier filtering.
0464 dataVecFilt = real(ifft(fft(dataVecPad).*fft(kernelVec)));
0465     
0466 % And then clip off the padding.
0467 smoothVec = dataVecFilt(padParams(2)+1:end - padParams(3));
0468 end
0469