Description of stat

0001 function r = stat_rfi(startMjd,endMjd,pArray_f,pArray_p,basedir)    
0002 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0003 %  function  [...] = stat_rfi(d)
0004 %
0005 %   Jan 26, 2012 (MAS) - Created.
0006 %
0007 %   Call stat_eval() a number of times, and run the values through the
0008 %   suggested flagging parameter values.  Write out results to a file in
0009 %   rfidir.
0010 %
0011 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0012 
0013 
0014 % How many time sections will we use?
0015 nTimes = round(10 * (endMjd - startMjd));
0016 
0017 
0018 
0019 % How many parameter sets to we have?
0020 nSets = size(pArray_f,2);
0021 
0022 % Calculate the 10% error parameters.
0023 pArrayA_f = 0.9 * pArray_f;
0024 pArrayA_f(3,:) = 1.1 * pArray_f(3,:);
0025 
0026 pArrayZ_f = 1.1 * pArray_f;
0027 pArrayZ_f(3,:) = 0.9 * pArray_f(3,:);
0028 
0029 pArrayA_p = 0.9 * pArray_p;
0030 pArrayA_p(3,:) = 1.1 * pArray_p(3,:);
0031 
0032 pArrayZ_p = 1.1 * pArray_p;
0033 pArrayZ_p(3,:) = 0.9 * pArray_p(3,:);
0034 
0035 
0036 
0037 outFile = [basedir 'stat_flag_results.mat'];
0038 if exist(outFile,'file')
0039     % Grab the existing stat flagging arrays (and helpers) from file:
0040     load(outFile);
0041 else
0042     % Prep the stat flagging arrays (and helpers) from scratch:
0043     statFlagArray_f = boolean([]);
0044     statFlagArrayA_f = boolean([]);
0045     statFlagArrayZ_f = boolean([]);
0046     statFlagArray_p = boolean([]);
0047     statFlagArrayA_p = boolean([]);
0048     statFlagArrayZ_p = boolean([]);
0049     statUtcVec = [];
0050     statFeatVec = int32([]);
0051 end
0052 
0053 % Loop over the time segments.
0054 disp('Looping over the time segments...');
0055 for k=1:nTimes
0056 
0057 
0058     startSec = mjd2string(startMjd + 0.1 * (k-1) - 1/24/60);
0059     endSec = mjd2string(startMjd + 0.1 * k + 1/24/60);
0060     
0061     disp(['Data from ' startSec ' to ' endSec]);
0062 
0063     % Give 5 minute leeway.
0064     if max(statUtcVec) > (tstr2mjd(endSec) - 0.0035)
0065         disp('Time period already ran... on to next one.');
0066         continue;
0067     end
0068         
0069     % Read in the requested data.
0070     %----------------------------------
0071     d = read_arc(startSec,endSec); 
0072 
0073 
0074     % Add to these output arrays:
0075     %----------------------------------
0076     statUtcVec = cat(1,statUtcVec,d.antenna0.receiver.utc);
0077     statFeatVec = cat(1,statFeatVec, ...
0078         interp1(d.array.frame.utc,double(d.array.frame.features),d.antenna0.receiver.utc,'nearest'));
0079     
0080     
0081     
0082     for m=1:2
0083 
0084         % Apply alpha and load corrections:
0085         %----------------------------------
0086         switch m
0087             case 1
0088                 d1 = assembleAlphaStreams(d,'FILTERED'); 
0089                 d1 = applyAlpha(d1,'FILTERED'); 
0090                 d1 = load_template_corr(d1,0,0);
0091             case 2
0092                 d1 = assembleAlphaStreams(d,'POLONLY'); 
0093                 d1 = applyAlpha(d1,'POLONLY'); 
0094                 d1 = load_template_corr(d1,2,0);
0095         end
0096 
0097         % OK.  Now let us compare with RFI flagging stats...
0098         %----------------------------------
0099         disp('Running statistical flagging...');
0100         [dataUtc, shortI_A, shortI_B, shortL_A, shortL_B, ...
0101             ~, longI_B, longI_B_min, longI_B_max, longL_A, longL_B] = ...
0102             stat_eval(d1.antenna0.receiver.utc,d1.antenna0.receiver.data(:,[1 6 7 8]));    
0103 
0104         
0105         % Grab the wanted parameter sets:
0106         %----------------------------------
0107         switch m
0108             case 1
0109                 pArray = pArray_f;
0110                 pArrayA = pArrayA_f;
0111                 pArrayZ = pArrayZ_f;
0112             case 2
0113                 pArray = pArray_p;
0114                 pArrayA = pArrayA_p;
0115                 pArrayZ = pArrayZ_p;
0116         end
0117         
0118         % Add to the output arrays.
0119         statFlagArray = boolean(zeros(length(dataUtc),nSets));
0120         statFlagArrayA = boolean(zeros(length(dataUtc),nSets));
0121         statFlagArrayZ = boolean(zeros(length(dataUtc),nSets));
0122 
0123         % Loop over the parameter sets, flagging on each:
0124         disp('Attempting the parameter cuts...');
0125         for l=1:nSets
0126 
0127             % Calculate the flags...
0128 
0129             % First flagging set...
0130             if min(pArray(1:3,l)) > 0
0131 
0132                 % Main parameters.
0133                 flagSA = (shortL_A > pArray(1,l)) & ...
0134                     (shortL_A > pArray(2,l) * shortI_A) & ...
0135                     (shortL_B > pArray(2,l) * shortI_B) & ...
0136                     (shortL_B < pArray(3,l) * shortL_A);
0137 
0138                 % Less strict parameters.
0139                 flagSA_A = (shortL_A > pArrayA(1,l)) & ...
0140                     (shortL_A > pArrayA(2,l) * shortI_A) & ...
0141                     (shortL_B > pArrayA(2,l) * shortI_B) & ...
0142                     (shortL_B < pArrayA(3,l) * shortL_A);
0143 
0144                 % More strict parameters.
0145                 flagSA_Z = (shortL_A > pArrayZ(1,l)) & ...
0146                     (shortL_A > pArrayZ(2,l) * shortI_A) & ...
0147                     (shortL_B > pArrayZ(2,l) * shortI_B) & ...
0148                     (shortL_B < pArrayZ(3,l) * shortL_A);
0149             else
0150                 flagSA = zeros(length(shortL_A),1);
0151                 flagSA_A = zeros(length(shortL_A),1);
0152                 flagSA_Z = zeros(length(shortL_A),1);
0153             end
0154 
0155 
0156             % Second flagging set...
0157             if pArray(4,l) > 0
0158 
0159                 % Main parameters.
0160                 flagSB = (shortL_A > pArray(4,l));
0161 
0162                 % Less strict parameters.
0163                 flagSB_A = (shortL_A > pArrayA(4,l));
0164 
0165                 % More strict parameters.
0166                 flagSB_Z = (shortL_A > pArrayZ(4,l));
0167             else
0168                 flagSB = zeros(length(shortL_A),1);
0169                 flagSB_A = zeros(length(shortL_A),1);
0170                 flagSB_Z = zeros(length(shortL_A),1);
0171             end
0172 
0173             % Third flagging set...
0174             if min(pArray(5:7,l)) > 0
0175 
0176                 % Main parameters.
0177                 flagLA = (longL_A > pArray(5,l)) & ...
0178                     (longL_B > pArray(6,l)) & ...
0179                     (longL_B > pArray(7,l) * longI_B);
0180 
0181                 % Less strict parameters.
0182                 flagLA_A = (longL_A > pArrayA(5,l)) & ...
0183                     (longL_B > pArrayA(6,l)) & ...
0184                     (longL_B > pArrayA(7,l) * longI_B);
0185 
0186                 % More strict parameters.
0187                 flagLA_Z = (longL_A > pArrayZ(5,l)) & ...
0188                     (longL_B > pArrayZ(6,l)) & ...
0189                     (longL_B > pArrayZ(7,l) * longI_B);        
0190             else
0191                 flagLA = zeros(length(shortL_A),1);
0192                 flagLA_A = zeros(length(shortL_A),1);
0193                 flagLA_Z = zeros(length(shortL_A),1);
0194             end
0195 
0196 
0197             % Fourth flagging set...
0198             if min(pArray(8:9,l)) > 0
0199 
0200                 % Main parameters.
0201                 flagLB = (longI_B_max > pArray(8,l)) & ...
0202                     (longI_B_max > pArray(9,l) * longI_B_min);
0203 
0204                 % Less strict parameters.
0205                 flagLB_A = (longI_B_max > pArrayA(8,l)) & ...
0206                     (longI_B_max > pArrayA(9,l) * longI_B_min);
0207 
0208                 % More strict parameters.
0209                 flagLB_Z = (longI_B_max > pArrayZ(8,l)) & ...
0210                     (longI_B_max > pArrayZ(9,l) * longI_B_min);
0211             else
0212                 flagLB = zeros(length(shortL_A),1);
0213                 flagLB_A = zeros(length(shortL_A),1);
0214                 flagLB_Z = zeros(length(shortL_A),1);
0215             end
0216 
0217 
0218 
0219             % Combine the flags...
0220             flagVec = flagSA | flagSB | flagLA | flagLB;
0221             flagVec_A = flagSA_A | flagSB_A | flagLA_A | flagLB_A;
0222             flagVec_Z = flagSA_Z | flagSB_Z | flagLA_Z | flagLB_Z;
0223 
0224 
0225             % Add the flags to the flag arrays.
0226             statFlagArray(:,l) = flagVec;
0227             statFlagArrayA(:,l) = flagVec_A;
0228             statFlagArrayZ(:,l) = flagVec_Z;
0229 
0230         end
0231     
0232         clear d1;
0233         
0234         % Add to these output arrays:
0235         %----------------------------------
0236         switch m
0237             case 1
0238                 statFlagArray_f = cat(1,statFlagArray_f,statFlagArray);
0239                 statFlagArrayA_f = cat(1,statFlagArrayA_f,statFlagArrayA);
0240                 statFlagArrayZ_f = cat(1,statFlagArrayZ_f,statFlagArrayZ);
0241             case 2
0242                 statFlagArray_p = cat(1,statFlagArray_p,statFlagArray);
0243                 statFlagArrayA_p = cat(1,statFlagArrayA_p,statFlagArrayA);
0244                 statFlagArrayZ_p = cat(1,statFlagArrayZ_p,statFlagArrayZ);
0245         end
0246         
0247     end
0248     
0249     clear d;
0250 
0251     
0252     %  Save intermediate results...
0253     disp('Saving intermediate statistical flagging results...');
0254 
0255     save(outFile, 'statUtcVec', ...
0256         'statFlagArray_f', 'statFlagArrayA_f', 'statFlagArrayZ_f', ...
0257         'statFlagArray_p', 'statFlagArrayA_p', 'statFlagArrayZ_p', ...
0258         'statFeatVec');
0259 
0260 end
0261 
0262 disp('Finished looping.');
0263 
0264 
0265 %  Remove duplicate time samples.
0266 disp('Now removing duplicate samples...');
0267 
0268 [statUtcVec, I, ~] = unique(statUtcVec,'first');
0269 
0270 statFlagArray_f = statFlagArray_f(I,:);
0271 statFlagArrayA_f = statFlagArrayA_f(I,:);
0272 statFlagArrayZ_f = statFlagArrayZ_f(I,:);
0273 statFlagArray_p = statFlagArray_p(I,:);
0274 statFlagArrayA_p = statFlagArrayA_p(I,:);
0275 statFlagArrayZ_p = statFlagArrayZ_p(I,:);
0276 statFeatVec = statFeatVec(I);
0277 
0278 
0279 
0280 % Write the stat output file.
0281 disp('Writing out the statistical flagging results...');
0282 
0283 save(outFile, 'statUtcVec', ...
0284     'statFlagArray_f', 'statFlagArrayA_f', 'statFlagArrayZ_f', ...
0285     'statFlagArray_p', 'statFlagArrayA_p', 'statFlagArrayZ_p', ...
0286     'statFeatVec');
0287 
0288 
0289 
0290 % Now we identify the candidate RFI events in preparation for writing to file:
0291 disp('Identifying RFI candidates...');
0292 
0293 % First, read in the fold candidates:
0294 foldFile_f = open([basedir 'fold_flag_results_filtered.mat']);
0295 %foldFile_p = open([basedir 'fold_flag_results_polonly.mat']);
0296 
0297 % Convert into a flagging vector of full length.  ; Both foldFlags should
0298 % be the saaaaame...
0299 foldFlagInd = round((foldFile_f.dataUtc - statUtcVec(1)) * 24 * 3600 * 100 + 1);
0300 foldFlagVec = zeros(length(statUtcVec),1);
0301 foldFlagVec(foldFlagInd) = 1;
0302 % foldFlagInd_p = round((foldFile_p.dataUtc - statUtcVec(1)) * 24 * 3600 * 100 + 1);
0303 % foldFlagVec_p = zeros(length(statUtcVec),1);
0304 % foldFlagVec_p(foldFlagInd_p) = 1;
0305 
0306 
0307 % Combine the fold and stat flagging to identify the candidates.
0308 totFlagArray = cat(2,foldFlagVec,statFlagArray_f,statFlagArray_p);
0309 totFlagVec = sum(totFlagArray,2) > 0;
0310 
0311 
0312 % Only keep those flag events which are not isolated.
0313 totFlagVec = totFlagVec & (smooth(totFlagVec,51) > 0.1);
0314 totFlagVec = smooth(totFlagVec,25) > 0;
0315 
0316 
0317 
0318 totFlagDiff = cat(1,0,diff(totFlagVec));
0319 
0320 totStart = find(totFlagDiff > 0);
0321 totEnd = find(totFlagDiff < 0);
0322 
0323 if totEnd(1) < totStart(1)
0324     totStart = cat(1,1,totStart);
0325 end
0326 if totEnd(end) < totStart(end)
0327     totEnd = cat(1,totEnd,length(totFlagVec));
0328 end
0329 
0330 nCandidates = length(totStart);
0331 
0332 % Loop over the RFI events, adding them to the output data...
0333 candidateUtc = zeros(nCandidates,1);
0334 candidateLen = zeros(nCandidates,1);
0335 candidateSet_f = zeros(nCandidates,nSets);
0336 candidateSetA_f = zeros(nCandidates,nSets);
0337 candidateSetZ_f = zeros(nCandidates,nSets);
0338 candidateSet_p = zeros(nCandidates,nSets);
0339 candidateSetA_p = zeros(nCandidates,nSets);
0340 candidateSetZ_p = zeros(nCandidates,nSets);
0341 candidateFold = zeros(nCandidates,1);
0342 candidateFeat = zeros(nCandidates,1);
0343 
0344 for k=1:nCandidates
0345 
0346     iStart = totStart(k);
0347     iEnd = totEnd(k);
0348     
0349     candidateUtc(k) = mean(statUtcVec(iStart:iEnd));
0350     candidateLen(k) = iEnd-iStart;
0351 
0352     candidateSet_f(k,:) = max(statFlagArray_f(iStart:iEnd,:));
0353     candidateSetA_f(k,:) = max(statFlagArrayA_f(iStart:iEnd,:));
0354     candidateSetZ_f(k,:) = max(statFlagArrayZ_f(iStart:iEnd,:));
0355     candidateSet_p(k,:) = max(statFlagArray_p(iStart:iEnd,:));
0356     candidateSetA_p(k,:) = max(statFlagArrayA_p(iStart:iEnd,:));
0357     candidateSetZ_p(k,:) = max(statFlagArrayZ_p(iStart:iEnd,:));
0358 
0359     candidateFold(k) = max(foldFlagVec(iStart:iEnd));
0360     
0361     candidateFeat(k) = mode(double(statFeatVec(iStart:iEnd)));
0362     
0363 end
0364 
0365 % Write the output file.
0366 disp('Writing out the RFI candidates...');
0367 
0368 outFile = [basedir 'rfi_candidates.mat'];
0369 %outFile = [basedir output];
0370 save(outFile, 'candidateUtc', 'candidateLen', ...
0371     'candidateSet_f', 'candidateSetA_f', 'candidateSetZ_f', ...
0372     'candidateSet_p', 'candidateSetA_p', 'candidateSetZ_p', ...
0373     'candidateFold', 'candidateFeat');
0374     
0375 
0376 r = 1;
0377 end
stat_rfi

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
