stat_eval

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  function  [...] = stat_eval(d)

   May 3, 2012 (MAS) - Branched from April 12, 2012 version of
   flagRFI_std.m


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

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