Description of flagRFI

0001 function ...
0002     [dataUtc, shortI_A, shortI_B, shortL_A, shortL_B, ...
0003     longI_A, longI_B, longI_B_min, longI_B_max, longL_A, longL_B] = ...
0004     flagRFI_eval(d,preFlag)    
0005 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0006 %  function  [flagShort flagLong] = flagRFI_eval(d,preFlag,params)
0007 %
0008 %   Nov 4, 2011 (MAS) - Branched from April 14, 2011 version of
0009 %   flagRFI_std().
0010 %
0011 %
0012 %   params(1): minimum shortflag 1
0013 %   params(2): minimum shortflag 2
0014 %   params(3): shortflag slope
0015 %   params(4): minimum shortflag (alt)
0016 %   params(5): minimum longflag 1 A
0017 %   params(6): minimum longflag 1 B
0018 %   params(7): minimum longflag 2
0019 %   params(8): minimum I longflag
0020 %   params(9): minimum I difference.
0021 %   params(10): final smoothing (seconds)
0022 %   params(11): How much to smooth out the I long.
0023 %
0024 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0025 
0026 
0027 %%%%%%%%%%%%%%
0028 % Parameters:
0029 %%%%%%%%%%%%%%
0030 
0031 % Number of iterations in 0.5s flagging.
0032 nClip = 5;
0033 
0034 % Two lengths of sub-vectors in first flagging.  50 is 0.5s.
0035 subLength_A = 25;
0036 subLength_B = 2 * subLength_A;
0037 % Smoothing length.  Units of subLength.  ie, 500 corresponds to 250s.
0038 subStdLength = 250; % 100
0039 
0040 % Two lengths of sub-vectors in second flagging.  1000 is 10s.
0041 subLength2_A = 250;
0042 subLength2_B = 4 * subLength2_A;
0043 
0044 
0045 
0046 
0047 %%%%%%%%%%%%%%%%%%%%%%%%%%
0048 % Begin the calculations!
0049 %%%%%%%%%%%%%%%%%%%%%%%%%%
0050 
0051 dataUtc = d.antenna0.receiver.utc;
0052 
0053 % Make the noise diode events look like non-noise diode events.
0054 d = noiseRemove(d);
0055 
0056 % Get its length.
0057 dLength = length(d.antenna0.receiver.data);
0058 
0059 
0060 % The length of the first StdDev vector.
0061 nSub_B = dLength / subLength_B;
0062 nSubR_B = floor(nSub_B);
0063 nSubR_A = nSubR_B * subLength_B / subLength_A;
0064 
0065 % The length of the 10s StdDev vector.
0066 nSub2_B = dLength / subLength2_B;
0067 nSubR2_B = floor(nSub2_B);
0068 nSubR2_A = nSubR2_B * subLength2_B / subLength2_A;
0069 
0070 
0071 
0072 % Now we begin the first flagging.
0073 %%%%%%%%%%%%%%%%%%%%%
0074 
0075 
0076 % There's going to be some smoothing done on the standard deviation vector.
0077 % We should prepare by calculating some length and padding parameters.
0078 sLength_A = 2 ^ ceil(log2(nSubR_A + 2 * 10 * subStdLength));
0079 sBuffStart_A = floor((sLength_A - nSubR_A) / 2);
0080 sBuffEnd_A = ceil((sLength_A - nSubR_A) / 2);
0081 
0082 sLength_B = 2 ^ ceil(log2(nSubR_B + 2 * 10 * subStdLength));
0083 sBuffStart_B = floor((sLength_B - nSubR_B) / 2);
0084 sBuffEnd_B = ceil((sLength_B - nSubR_B) / 2);
0085 
0086 
0087 % Also, let's make the kernel for the convolution.
0088 xKernel_A = (1:sLength_A)' - 1;
0089 gaussKernelStd_A = ...
0090     (exp(-0.5*(xKernel_A / subStdLength).^2) + ...
0091     exp(-0.5*((sLength_A - xKernel_A) / subStdLength).^2)) / ...
0092     (subStdLength*sqrt(2*pi));
0093 
0094 xKernel_B = (1:sLength_B)' - 1;
0095 gaussKernelStd_B = ...
0096     (exp(-0.5*(xKernel_B / subStdLength).^2) + ...
0097     exp(-0.5*((sLength_B - xKernel_B) / subStdLength).^2)) / ...
0098     (subStdLength*sqrt(2*pi));
0099 
0100 % Also, an X array which will later be used for interpolation purposes.
0101 xStd_A = (1:nSubR_A);
0102 xStd_B = (1:nSubR_B);
0103 
0104 
0105 
0106 % Prep the resulting arrays from this endeavor.
0107 dataStdNorm_A = zeros(nSubR_A,6);
0108 dataStdNorm_B = zeros(nSubR_B,6);
0109 
0110 dataStdVec2_A = zeros(nSubR2_A,6);
0111 dataStdVec2_B = zeros(nSubR2_B,6);
0112 
0113 
0114 % This will only proceed with a subvector of the data, determined by
0115 % fitting in an integral number of subLength_B.
0116 dataRStart = floor((dLength - nSubR_B * subLength_B) / 2) + 1;
0117 dataREnd = dataRStart + nSubR_B * subLength_B - 1;
0118 
0119 dataR2Start = floor((dLength - nSubR2_B * subLength2_B) / 2) + 1;
0120 dataR2End = dataR2Start + nSubR2_B * subLength2_B - 1;
0121 
0122 
0123 % Prep the positional flags.
0124 preFlagR = max(reshape(preFlag(dataRStart:dataREnd),subLength_B,nSubR_B));
0125 
0126 
0127 % OK.  Loop over all six polarizations, flagging on the linear
0128 % polarization, and also getting the smoothed overall RMS trends (which
0129 % will be used in the next step).
0130 for k=1:6
0131     
0132     % Calculate the standard deviations!
0133     dataStdVec_A = std(reshape(d.antenna0.receiver.data(dataRStart:dataREnd,k), subLength_A, nSubR_A));
0134     dataStdVec_B = std(reshape(d.antenna0.receiver.data(dataRStart:dataREnd,k), subLength_B, nSubR_B));
0135 
0136     % In case nSub is not an integer...
0137     if nSub_B > nSubR_B
0138         dataStdVec_A(1) = ...
0139             std(d.antenna0.receiver.data(1:dataRStart+subLength_A-1,k));
0140         dataStdVec_A(end) = ...
0141             std(d.antenna0.receiver.data(dataREnd-subLength_A+1:end,k));    
0142 
0143         dataStdVec_B(1) = ...
0144             std(d.antenna0.receiver.data(1:dataRStart+subLength_B-1,k));
0145         dataStdVec_B(end) = ...
0146             std(d.antenna0.receiver.data(dataREnd-subLength_B+1:end,k));    
0147     end    
0148     
0149     
0150     %-----
0151     % So here's the actual short flagging part.  Done in easy steps.
0152     %-----
0153     
0154     % Step 1:  Clip immediate outliers using both arrays.
0155     %-----
0156     flagVec_A = stdClip(dataStdVec_A,3,nClip,1);
0157     flagVec_B = stdClip(dataStdVec_B,3,nClip,1);
0158     
0159             
0160     % Step 2:  Ignoring outliers, calculate running median.
0161     %-----
0162     
0163     % Fold the positional flagging into the flag vector for this part:
0164     flagVec_A_t = repmat(preFlagR == 1,floor(subLength_B/subLength_A),1);
0165     flagVec_A = smooth((flagVec_A == 1) | flagVec_A_t(:)',40)' > 0;
0166     flagVec_B = smooth((flagVec_B == 1) | (preFlagR == 1),20)' > 0;
0167 
0168     
0169     
0170     runningMed_A = ...
0171         padSmooth(xStd_A,dataStdVec_A,gaussKernelStd_A,flagVec_A, ...
0172         [sLength_A sBuffStart_A sBuffEnd_A]);
0173     runningMed_B = ...
0174         padSmooth(xStd_B,dataStdVec_B,gaussKernelStd_B,flagVec_B, ...
0175         [sLength_B sBuffStart_B sBuffEnd_B]);
0176        
0177     
0178         
0179     % Step 3:  Subtract running median.
0180     %-----
0181     dataStdVec_A = dataStdVec_A - runningMed_A';
0182     dataStdVec_B = dataStdVec_B - runningMed_B';
0183     
0184     % Step 4:  Clip current outliers.
0185     %-----
0186     flagVec_A = stdClip(dataStdVec_A,5,nClip,0);
0187     flagVec_B = stdClip(dataStdVec_B,5,nClip,0);
0188 
0189     
0190     % Step 5:  Normalize the standard deviation vectors.
0191     %-----
0192     dataStdNorm_A(:,k) = dataStdVec_A / std(dataStdVec_A(flagVec_A == 0));
0193     dataStdNorm_B(:,k) = dataStdVec_B / std(dataStdVec_B(flagVec_B == 0));
0194 
0195     
0196     % While we're here, let's go ahead and calculate another couple of
0197     % standard deviation vectors.
0198     % These ones are for use with the longer flagging.
0199     % Note that, at this point, we go ahead and divide by the running
0200     % median and subtract by the overall remaining mean (~1).
0201     
0202     % Interpolate the running median a bit more.
0203     runningMedLong = ...
0204         interp1(((1:length(runningMed_A))' - 0.5) * subLength_A, ...
0205         runningMed_A,(1:dLength),'pchip');
0206 
0207     
0208 
0209     dataStdVec2_A(:,k) = ...
0210         std(detrend(reshape(d.antenna0.receiver.data(dataR2Start:dataR2End,k) ./ ...
0211         runningMedLong(dataR2Start:dataR2End)', ...
0212         subLength2_A, nSubR2_A)));
0213     
0214     dataStdVec2_B(:,k) = ...
0215         std(detrend(reshape(d.antenna0.receiver.data(dataR2Start:dataR2End,k) ./ ...
0216         runningMedLong(dataR2Start:dataR2End)', ...
0217         subLength2_B, nSubR2_B)));
0218 
0219     
0220     
0221     flagStdVec2_A = stdClip(dataStdVec2_A(:,k),1,nClip,1);
0222     flagStdVec2_B = stdClip(dataStdVec2_B(:,k),1,nClip,1);
0223     
0224     dataStdVec2_A(:,k) = dataStdVec2_A(:,k) - ...
0225         mean(dataStdVec2_A(flagStdVec2_A == 0,k));
0226     dataStdVec2_B(:,k) = dataStdVec2_B(:,k) - ...
0227         mean(dataStdVec2_B(flagStdVec2_B == 0,k));
0228     
0229     
0230     % In case nSub2_B is not an integer...
0231     if nSub2_B > nSubR2_B
0232 
0233         dataStdVec2_A(1) = ...
0234             std(detrend(d.antenna0.receiver.data(1:dataR2Start+subLength2_A-1,k) ./ ...
0235             runningMedLong(1:dataR2Start+subLength2_A-1)')) - 1;
0236         dataStdVec2_A(end) = ...
0237             std(detrend(d.antenna0.receiver.data(dataR2End-subLength2_A+1:end,k) ./ ...
0238             runningMedLong(dataR2End-subLength2_A+1:end)')) - 1;    
0239 
0240         dataStdVec2_B(1) = ...
0241             std(detrend(d.antenna0.receiver.data(1:dataR2Start+subLength2_B-1,k) ./ ...
0242             runningMedLong(1:dataR2Start+subLength2_B-1)')) - 1;
0243         dataStdVec2_B(end) = ...
0244             std(detrend(d.antenna0.receiver.data(dataR2End-subLength2_B+1:end,k) ./ ...
0245             runningMedLong(dataR2End-subLength2_B+1:end)')) - 1;    
0246     end    
0247 end
0248 
0249 
0250 % Finish off the short flagging.
0251 %%%%%%%%%%%%%%%%
0252 
0253 % First off, we need to logically combine the polarization vectors.
0254 % Further, we'll do a quick interpolation of the B vector to make it the
0255 % same length as A.
0256 dataStdLogic_A = max(min(dataStdNorm_A(:,2:2:4),[],2), ...
0257     min(dataStdNorm_A(:,3:2:5),[],2));
0258 
0259 dataStdI_A = mean(dataStdNorm_A(:,1:5:6),2);
0260 
0261 
0262 dataStdLogic_B = repmat(max(min(dataStdNorm_B(:,2:2:4),[],2), ...
0263     min(dataStdNorm_B(:,3:2:5),[],2))',subLength_B/subLength_A,1);
0264 dataStdLogic_B = dataStdLogic_B(:);
0265 
0266 dataStdI_B = repmat(mean(dataStdNorm_B(:,1:5:6),2)', ...
0267     subLength_B/subLength_A,1);
0268 dataStdI_B = dataStdI_B(:);
0269 
0270 
0271 
0272 % flagLogic = ...
0273 %     ((dataStdLogic_A > params(1)) & ...
0274 %     (dataStdLogic_A > params(2) * dataStdI_A) & ...
0275 %     (dataStdLogic_B > params(2) * dataStdI_B) & ...
0276 %     (dataStdLogic_B < params(3) * dataStdLogic_A)) | ...
0277 %     dataStdLogic_A > params(4);
0278     
0279 
0280 % % Stretch it out.
0281 % flagRep = zeros(dLength,1);
0282 % flagRep(dataRStart:dataREnd) = ...
0283 %     reshape(repmat(flagLogic',subLength_A,1),subLength_A*nSubR_A,1);
0284 %
0285 % if dataRStart > 1
0286 %     flagRep(1:dataRStart-1) = flagLogic(1);
0287 %     flagRep(dataREnd+1:end) = flagLogic(end);
0288 % end
0289 
0290 
0291 
0292 
0293 % Now we do the longer flagging.
0294 %%%%%%%%%%%%%%%%%
0295 
0296 
0297 % First off, we need to smooth the vectors.  At the same time, we'll add
0298 % the matching polarization channels together.
0299 dataStdI2_A = mean(dataStdVec2_A(:,1:5:6),2);
0300 dataStdLogic2_A = max(sum(dataStdVec2_A(:,2:2:4),2) / 2, ...
0301     sum(dataStdVec2_A(:,3:2:5),2) / 2);
0302 
0303 dataStdI2_B = mean(dataStdVec2_B(:,1:5:6),2);
0304 dataStdI2_B_max = max(dataStdVec2_B(:,1:5:6),[],2);
0305 dataStdI2_B_min = min(dataStdVec2_B(:,1:5:6),[],2);
0306 
0307 dataStdLogic2_B = max(sum(dataStdVec2_B(:,2:2:4),2) / 2, ...
0308     sum(dataStdVec2_B(:,3:2:5),2) / 2);
0309 
0310 dataStdLogic2_B = repmat(dataStdLogic2_B',subLength2_B/subLength2_A,1);
0311 dataStdLogic2_B = dataStdLogic2_B(:);
0312 
0313 dataStdI2_B = repmat(dataStdI2_B',subLength2_B/subLength2_A,1);
0314 dataStdI2_B = dataStdI2_B(:);
0315 
0316 dataStdI2_B_max = repmat(dataStdI2_B_max',subLength2_B/subLength2_A,1);
0317 dataStdI2_B_max = dataStdI2_B_max(:);
0318 dataStdI2_B_min = repmat(dataStdI2_B_min',subLength2_B/subLength2_A,1);
0319 dataStdI2_B_min = dataStdI2_B_min(:);
0320 
0321 
0322 
0323 
0324 %if params(11) > 0
0325 %    for k=1:params(11)
0326 %        dataStdI2_B = cat(2, dataStdI2_B, ...
0327 %            cat(1,dataStdI2_B(1+k:end,1),zeros(k,1)), ...
0328 %            cat(1,zeros(k,1),dataStdI2_B(1:end-k,1)));
0329 %    end
0330     
0331 dataStdI2_B = max(dataStdI2_B,[],2);
0332 %end
0333 
0334 
0335 
0336 % % OK, now we flag on these fellows.
0337 % dataStdFlag = ((dataStdLogic2_A > params(5)) & ...
0338 %     (dataStdLogic2_B > params(6)) & ...
0339 %     (dataStdLogic2_B > params(7) * dataStdI2_B)) | ...
0340 %     ((dataStdI2_B_max > params(8)) & ...
0341 %     (dataStdI2_B_max > params(9) * dataStdI2_B_min));
0342 
0343 
0344 % % Stretch these vectors:
0345 % dataStdI_A
0346 % dataStdI_B
0347 % dataStdLogic_A
0348 % dataStdLogic_B
0349 %
0350 % dataStdI2_A
0351 % dataStdI2_B
0352 % dataStdLogic2_A
0353 % dataStdLogic2_B
0354 % dataStdI2_B_max
0355 % dataStdI2_B_min
0356 
0357 
0358 
0359 % Stretch these out:
0360 %%%%%%%%%%%%%%%%%%%%%
0361 shortI_A = zeros(dLength,1);
0362 shortI_B = zeros(dLength,1);
0363 shortL_A = zeros(dLength,1);
0364 shortL_B = zeros(dLength,1);
0365 
0366 longI_A = zeros(dLength,1);
0367 longI_B = zeros(dLength,1);
0368 longI_B_min = zeros(dLength,1);
0369 longI_B_max = zeros(dLength,1);
0370 longL_A = zeros(dLength,1);
0371 longL_B = zeros(dLength,1);
0372 
0373 
0374 
0375 
0376 
0377 % Do it for the short ones:
0378 shortI_A(dataRStart:dataREnd) = ...
0379     reshape(repmat(dataStdI_A',subLength_A,1),subLength_A*nSubR_A,1);
0380 shortI_B(dataRStart:dataREnd) = ...
0381     reshape(repmat(dataStdI_B',subLength_A,1),subLength_A*nSubR_A,1);
0382 shortL_A(dataRStart:dataREnd) = ...
0383     reshape(repmat(dataStdLogic_A',subLength_A,1),subLength_A*nSubR_A,1);
0384 shortL_B(dataRStart:dataREnd) = ...
0385     reshape(repmat(dataStdLogic_B',subLength_A,1),subLength_A*nSubR_A,1);
0386 
0387 if dataRStart > 1
0388     shortI_A(1:dataRStart-1) = shortI_A(1);
0389     shortI_A(dataREnd+1:end) = shortI_A(end);
0390 
0391     shortI_B(1:dataRStart-1) = shortI_B(1);
0392     shortI_B(dataREnd+1:end) = shortI_B(end);
0393 
0394     shortL_A(1:dataRStart-1) = shortL_A(1);
0395     shortL_A(dataREnd+1:end) = shortL_A(end);
0396 
0397     shortL_B(1:dataRStart-1) = shortL_B(1);
0398     shortL_B(dataREnd+1:end) = shortL_B(end);
0399 end    
0400 
0401 
0402 % Now for the long ones:
0403 longI_A(dataR2Start:dataR2End) = ...
0404     reshape(repmat(dataStdI2_A',subLength2_A,1),subLength2_A*nSubR2_A,1);
0405 longI_B(dataR2Start:dataR2End) = ...
0406     reshape(repmat(dataStdI2_B',subLength2_A,1),subLength2_A*nSubR2_A,1);
0407 longI_B_min(dataR2Start:dataR2End) = ...
0408     reshape(repmat(dataStdI2_B_min',subLength2_A,1),subLength2_A*nSubR2_A,1);
0409 longI_B_max(dataR2Start:dataR2End) = ...
0410     reshape(repmat(dataStdI2_B_max',subLength2_A,1),subLength2_A*nSubR2_A,1);
0411 longL_A(dataR2Start:dataR2End) = ...
0412     reshape(repmat(dataStdLogic2_A',subLength2_A,1),subLength2_A*nSubR2_A,1);
0413 longL_B(dataR2Start:dataR2End) = ...
0414     reshape(repmat(dataStdLogic2_B',subLength2_A,1),subLength2_A*nSubR2_A,1);
0415 
0416 if dataR2Start > 1
0417     longI_A(1:dataR2Start-1) = longI_A(1);
0418     longI_A(dataR2End+1:end) = longI_A(end);
0419 
0420     longI_B(1:dataR2Start-1) = longI_B(1);
0421     longI_B(dataR2End+1:end) = longI_B(end);
0422 
0423     longI_B_min(1:dataR2Start-1) = longI_B_min(1);
0424     longI_B_min(dataR2End+1:end) = longI_B_min(end);
0425 
0426     longI_B_max(1:dataR2Start-1) = longI_B_max(1);
0427     longI_B_max(dataR2End+1:end) = longI_B_max(end);
0428 
0429     longL_A(1:dataR2Start-1) = longL_A(1);
0430     longL_A(dataR2End+1:end) = longL_A(end);
0431 
0432     longL_B(1:dataR2Start-1) = longL_B(1);
0433     longL_B(dataR2End+1:end) = longL_B(end);
0434 end    
0435 
0436 
0437 % % Stretch it out.
0438 % flagStdRep = zeros(dLength,1);
0439 % flagStdRep(dataR2Start:dataR2End) = ...
0440 %     reshape(repmat(dataStdFlag',subLength2_A,1),subLength2_A*nSubR2_A,1);
0441 %
0442 % if dataR2Start > 1
0443 %     flagStdRep(1:dataR2Start-1) = dataStdFlag(1);
0444 %     flagStdRep(dataR2End+1:end) = dataStdFlag(end);
0445 % end
0446 
0447 
0448 
0449 % % Alright, let's put these two flags together.
0450 % % First we'll make sure they haven't already been flagged by position and
0451 % % make sure they're longer than 25 samples.
0452 % % Then flag their surroundings.
0453 % flagCombine = (flagRep | flagStdRep) & ~preFlag;
0454 %
0455 %
0456 % flagCombine = flagCombine & ...
0457 %     cat(1,ones(12,1),flagCombine(1:end-12)) & ...
0458 %     cat(1,flagCombine(13:end),ones(12,1));
0459 %
0460 %
0461 % if params(10) > 0
0462 %     flag = (convn(double(flagCombine),ones(100*params(10),1),'same') > 0);
0463 % else
0464 %     flag = (flagCombine > 0);
0465 % end
0466 
0467 
0468 end
0469 
0470 
0471 % The sigma clippling subroutine.
0472 function clipVec = stdClip(dataVec,thresh,nClip,medSubtract)
0473 
0474 % Either we subtract the median or we don't.
0475 if medSubtract ~= 0
0476     medSubtract = 1;
0477 end
0478 
0479 % Initialize the clip vector.
0480 clipVec = zeros(length(dataVec),1);
0481 
0482 % Do the clipping...
0483 for l=1:nClip
0484     
0485     dataStdStd = std(dataVec(clipVec == 0));
0486     dataStdMed = median(dataVec(clipVec == 0));
0487     
0488     clipVec = abs(dataVec - medSubtract * dataStdMed) > ...
0489         thresh * dataStdStd;
0490     
0491 end
0492 
0493 end
0494 
0495 
0496 function smoothVec = padSmooth(xVec,dataVec,kernelVec,flagVec,padParams)
0497 
0498 
0499 % Interpolate over the flagged regions.
0500 dataVecT = interp1(xVec(flagVec == 0), ...
0501     smooth(dataVec(flagVec == 0),ceil(length(dataVec) / 100)), ...
0502     (1:length(dataVec)), 'pchip');
0503 
0504 dataVecF = dataVec .* (1 - flagVec) + flagVec .* dataVecT;
0505 
0506 if flagVec(1) == 1
0507     noFlagPos = find(flagVec == 0,1,'first');
0508     
0509     dataVecF(1:noFlagPos-1) = mean(dataVecF(noFlagPos:noFlagPos+99));
0510 end    
0511 if flagVec(end) == 1
0512     noFlagPos = find(flagVec == 0,1,'last');
0513     
0514     dataVecF(noFlagPos+1:end) = mean(dataVecF(noFlagPos-99:noFlagPos));
0515 end    
0516     
0517 
0518 
0519 % Pad the data up to a power of two.
0520 dataVecPad = zeros(padParams(1),1);
0521 dataVecPad(1:padParams(2)) = ...
0522     repmat(median(dataVecF(1:10)),padParams(2),1);
0523 dataVecPad(padParams(2)+1:padParams(2)+length(dataVec)) = dataVecF;
0524 dataVecPad(end-padParams(3)+1:end) = ...
0525     repmat(median(dataVecF(end-9:end)),padParams(3),1);
0526 
0527 % Perform a Fourier filtering.
0528 dataVecFilt = real(ifft(fft(dataVecPad).*fft(kernelVec)));
0529     
0530 % And then clip off the padding.
0531 smoothVec = dataVecFilt(padParams(2)+1:end - padParams(3));
0532 end
0533
flagRFI_eval

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
