plot_rfi

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  function  [...] = plot_rfi(xxx)

   Jan 26, 2012 (MAS) - Created.

   Plot the individual RFI candidates; user is responsible for determining
   which are real and which are not.

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

0001 function r = plot_rfi(startMjd,endMjd,basedir)    
0002 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0003 %  function  [...] = plot_rfi(xxx)
0004 %
0005 %   Jan 26, 2012 (MAS) - Created.
0006 %
0007 %   Plot the individual RFI candidates; user is responsible for determining
0008 %   which are real and which are not.
0009 %
0010 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0011 
0012 
0013 
0014 % % Read in the three files.
0015 % statFile = open([basedir 'stat_flag_results.mat']);
0016 %
0017 % statUtc = statFile.statUtcVec;
0018 % statFlag = statFile.statFlagArray;
0019 %
0020 %
0021 %
0022 % foldFile = open([basedir 'fold_flag_results.mat']);
0023 %
0024 % foldFlagInd = round((foldFile.dataUtc - statUtc(1)) * 24 * 3600 * 100 + 1);
0025 % foldFlag = zeros(length(statUtc),1);
0026 % foldFlag(foldFlagInd) = 1;
0027 
0028 
0029 % Read in the candidate file:
0030 candidateFile = open([basedir 'rfi_candidates.mat']);
0031 %disp(candidateFile);
0032 
0033 candidateUtc = candidateFile.candidateUtc;
0034 candidateLen = candidateFile.candidateLen;
0035 candidateSet_f = candidateFile.candidateSet_f;
0036 candidateSetA_f = candidateFile.candidateSetA_f;
0037 candidateSetZ_f = candidateFile.candidateSetZ_f;
0038 candidateSet_p = candidateFile.candidateSet_p;
0039 candidateSetA_p = candidateFile.candidateSetA_p;
0040 candidateSetZ_p = candidateFile.candidateSetZ_p;
0041 candidateFold_f = candidateFile.candidateFold_f;
0042 candidateFold_p = candidateFile.candidateFold_p;
0043 %candidateFold = candidateFile.candidateFold;
0044 candidateFeat = candidateFile.candidateFeat;
0045 
0046 nCandidates = length(candidateUtc);
0047 
0048 candidateVer = zeros(nCandidates,1);
0049 
0050 
0051 
0052 % How many time sections will we use?
0053 nTimes = round(50 * (endMjd - startMjd));
0054 
0055 
0056 % Initialize the inner iterator:
0057 l = 1;
0058 
0059 
0060 % Loop over the time segments.
0061 disp('Looping over the time segments...');
0062 for k=1:nTimes
0063 
0064 
0065     startSec = mjd2string(startMjd + 0.02 * (k-1) - 5/24/60);
0066 %    endSec = mjd2string(tstr2mjd(startSec)+1/6/24);
0067     endSec = mjd2string(startMjd + 0.02 * k + 5/24/60);
0068     
0069     disp(['Data from ' startSec ' to ' endSec]);
0070 
0071     
0072     % Make sure we have got something here, before reading in...
0073     if candidateUtc(l) >= tstr2mjd(endSec)
0074         disp('  No RFI candidates in this range...');
0075         continue;
0076     end
0077     
0078     
0079     % Read in the data, but just the interesting parts...
0080     [d, pointFlag, galacticFlag, sunFlag, moonFlag] = ...
0081         read_posflag(startSec,endSec);
0082 %        read_posflag(startSec,endSec);
0083     
0084     
0085     % Apply the alpha and 1.2hz corrections:
0086     disp('Applying the alpha step.');
0087     d = assembleAlphaStreams(d,'FILTERED');
0088 %     d = applyAlpha(d,'FILTERED');
0089 %     disp('Applying the 1.2Hz correction step.');
0090 %     d = load_template_corr(d,0,0);
0091 
0092 %     % Remove the noise diode events:
0093 %     disp('Removing noise diode events.');
0094 %     d = noiseRemove(d);
0095     
0096     % The load and noise diode removals have not been applied.
0097     loadApply = 0;
0098     noiseApply = 0;
0099 
0100     
0101     utcVec = d.antenna0.receiver.utc;
0102     I1Vec = -d.antenna0.receiver.data(:,1);
0103     QVec = d.antenna0.receiver.data(:,6);
0104     UVec = d.antenna0.receiver.data(:,7);
0105     I2Vec = -d.antenna0.receiver.data(:,8);
0106     
0107     dataLen = length(utcVec);    
0108 
0109 
0110     % Grab the predicted sky model.
0111     disp('Get the predicted sky model.');
0112     sky_predict = read_skymodel(d.antenna0.servo.equa);
0113     
0114 
0115     printhelp;
0116     
0117     % Loop over candidates.
0118     while (l <= nCandidates && candidateUtc(l) < tstr2mjd(endSec))
0119 
0120         disp(['Candidate ' int2str(l) ' out of ' int2str(nCandidates)]);
0121         
0122         % Candidate Starting Index:
0123         canStart = max(1, ...
0124             round((candidateUtc(l) - utcVec(1)) * 24 * 3600 * 100) - ...
0125             5*candidateLen(l));
0126         % Candidate Ending Index:
0127         canEnd = min(dataLen, ...
0128             round((candidateUtc(l) - utcVec(1)) * 24 * 3600 * 100) + ...
0129             5*candidateLen(l));
0130         % Candidate Centre Index:
0131         canCen = round((candidateUtc(l) - utcVec(1)) * 24 * 3600 * 100) + 1;
0132         
0133 
0134         I1Med = median(I1Vec(canStart:canEnd));
0135         I1Range = range(I1Vec(canStart:canEnd));
0136         QMed = median(QVec(canStart:canEnd));
0137         QRange = range(QVec(canStart:canEnd));
0138         UMed = median(UVec(canStart:canEnd));
0139         URange = range(UVec(canStart:canEnd));
0140         I2Med = median(I2Vec(canStart:canEnd));
0141         I2Range = range(I2Vec(canStart:canEnd));
0142         
0143         modelMed = median(sky_predict(canStart:canEnd));
0144         modelRange = range(sky_predict(canStart:canEnd));
0145         if modelRange == 0
0146             modelRange = 1;
0147         end
0148         disp(modelRange);
0149         
0150         maxRange = max([I1Range QRange URange I2Range]);
0151         disp(maxRange);
0152         
0153         
0154         % Initialize this guy.
0155         canY = 1;
0156         
0157         % Create the figure window (if necessary).
0158         figure(gcf);
0159 
0160 
0161         p = 1;
0162         while (p > 0)
0163 
0164             % Clear the figure window.
0165             clf;
0166 
0167 
0168             % The Time vector.
0169             tVec = (utcVec(canStart:canEnd) - candidateUtc(l)) * 24 * 3600;
0170 
0171             % Generate the axes...
0172             axis([min(tVec), max(tVec), -1.5*maxRange*canY, 1.5*maxRange*canY]);
0173             hold on;
0174 
0175             
0176 
0177             % Plot the position flags.
0178             flagplot(tVec,sunFlag(canStart:canEnd),[0.8 0.7 0.7],maxRange*canY);
0179             flagplot(tVec,moonFlag(canStart:canEnd),[0.7 0.7 0.8],maxRange*canY);
0180             flagplot(tVec,galacticFlag(canStart:canEnd),[0.7 0.8 0.7],maxRange*canY);
0181             flagplot(tVec,pointFlag(canStart:canEnd),[0.7 0.7 0.7],maxRange*canY);
0182 
0183             
0184             % Plot the predicted data
0185             % (from de Oliveira-Costa's all-sky model.)
0186             plot(tVec,(sky_predict(canStart:canEnd) - modelMed)/modelRange*maxRange,'color',[0.5, 0.5, 0.5],'linewidth',5);
0187             
0188             
0189             % Plot the candidate.
0190             plot(tVec,I1Vec(canStart:canEnd) - I1Med + maxRange*canY,'b');
0191             plot(tVec,QVec(canStart:canEnd) - QMed + maxRange/3*canY,'r');
0192             plot(tVec,UVec(canStart:canEnd) - UMed - maxRange/3*canY,'g');
0193             plot(tVec,I2Vec(canStart:canEnd) - I2Med - maxRange*canY,'m');
0194 
0195             
0196             % Show the RFI region.
0197             plot([-candidateLen(l)/200, -candidateLen(l)/200, ...
0198                 candidateLen(l)/200, candidateLen(l)/200], ...
0199                 [-1.5*maxRange*canY, 1.5*maxRange*canY, ...
0200                 1.5*maxRange*canY, -1.5*maxRange*canY], 'k');
0201             
0202             
0203             legend('Sun', 'Moon', 'Galaxy', 'Point Sources', ...
0204                 'Sky Model', ...
0205                 'I1', 'Q', 'U', 'I2', 'RFI?');
0206             xlabel('Time (s)');
0207             ylabel('Power (s)');
0208             title(mjd2string(candidateUtc(l)));
0209             
0210             hold off;
0211 
0212             % Now wait for the user to give a command...
0213             w = 0;
0214             while (~w)
0215                 w = waitforbuttonpress;
0216             end
0217 
0218             [p canStart canEnd canY procCommand verdict] = evalkey(canStart,canEnd,canCen,canY,dataLen);
0219 
0220             % If there was a DSP command:
0221             if procCommand ~= 0
0222                 if procCommand == 1 && noiseApply == 0
0223                     disp('Removing noise diode events...');
0224                     d = noiseRemove(d);
0225                     disp('   ...done');
0226 
0227                     noiseApply = 1;
0228                 elseif procCommand == 1
0229                     disp('  Noise diode removed already.');
0230 
0231                 elseif procCommand == 2 && loadApply == 0
0232                     disp('Removing 1.2Hz...');
0233                     d = load_template_corr(d,0,0);
0234                     disp('   ...done');
0235 
0236                     loadApply = 1;
0237                 elseif procCommand == 2
0238                     disp('  1.2Hz removed already.');
0239 
0240                 end
0241                 
0242 
0243                 I1Vec = -d.antenna0.receiver.data(:,1);
0244                 QVec = d.antenna0.receiver.data(:,6);
0245                 UVec = d.antenna0.receiver.data(:,7);
0246                 I2Vec = -d.antenna0.receiver.data(:,8);
0247 
0248                 I1Med = median(I1Vec(canStart:canEnd));
0249                 I1Range = range(I1Vec(canStart:canEnd));
0250                 QMed = median(QVec(canStart:canEnd));
0251                 QRange = range(QVec(canStart:canEnd));
0252                 UMed = median(UVec(canStart:canEnd));
0253                 URange = range(UVec(canStart:canEnd));
0254                 I2Med = median(I2Vec(canStart:canEnd));
0255                 I2Range = range(I2Vec(canStart:canEnd));
0256             
0257                 maxRange = max([I1Range QRange URange I2Range]);
0258 
0259             end            
0260         end
0261 
0262         % Record the verdict.
0263         candidateVer(l) = verdict;
0264 
0265         if p == -1
0266             l = max([1 (l-1)]);
0267         else
0268             l = l + 1;
0269         end
0270 
0271     end
0272     
0273 end
0274 
0275 
0276 % Write the output.
0277 save([basedir 'rfi_candidates.mat'], 'candidateUtc', 'candidateLen', ...
0278     'candidateSet_f', 'candidateSetA_f', 'candidateSetZ_f', ...
0279     'candidateSet_p', 'candidateSetA_p', 'candidateSetZ_p', ...
0280     'candidateFold_f', 'candidateFold_p', 'candidateFeat', 'candidateVer');
0281 %    'candidateFold', 'candidateFeat', 'candidateVer');
0282 
0283 
0284 r = 1;
0285 
0286 end
0287 
0288 
0289 
0290 %%%%%%%%%%%%%%%%%%%%%%%%%
0291 function flagplot(tVec,flagVec,flagColour,maxRange)
0292 
0293 
0294 % Do we even have any such flags...?
0295 if max(flagVec) > 0
0296 
0297 
0298     % Identify the flagged regions...
0299     dFlag = diff(flagVec);
0300 
0301     flagStart = find(dFlag > 0);
0302     flagEnd = find(dFlag < 0);
0303 
0304     if isempty(flagStart)
0305         flagStart = 1;
0306     end
0307     if isempty(flagEnd)
0308         flagEnd = length(flagVec);
0309     end
0310     if flagEnd(1) < flagStart(1)
0311         flagStart = cat(1,1,flagStart);
0312     end
0313     if flagEnd(end) < flagStart(end)
0314         flagEnd = cat(1,flagEnd,length(flagVec));
0315     end
0316 
0317 
0318     nFlag = length(flagStart);
0319 
0320     % Prepare the vertices for the filled region...
0321     xVec = zeros(4*nFlag,1);
0322     yVec = zeros(4*nFlag,1);
0323     for k=1:nFlag
0324         xVec(4*k-3) = tVec(flagStart(k));
0325         xVec(4*k-2) = tVec(flagStart(k));
0326         xVec(4*k-1) = tVec(flagEnd(k));
0327         xVec(4*k) = tVec(flagEnd(k));
0328 
0329         yVec(4*k-3) = -1.5*maxRange;
0330         yVec(4*k-2) = 1.5*maxRange;
0331         yVec(4*k-1) = 1.5*maxRange;
0332         yVec(4*k) = -1.5*maxRange;
0333     end
0334 
0335 else
0336 
0337 
0338     xVec = [0 0];
0339     yVec = [-2*maxRange -2*maxRange];
0340 
0341 end
0342 
0343 
0344 fill(xVec,yVec,flagColour,'EdgeColor',flagColour);
0345 
0346 end
0347 
0348 
0349 
0350 %%%%%%%%%%%%%%%%%%%%%%%%%
0351 function [p canStart canEnd canY procCommand verdict] = evalkey(canStart,canEnd,canCen,canY,dataLen)
0352 
0353 verdict = 0;
0354 procCommand = 0;
0355 
0356 ch=get(gcf,'CurrentCharacter');
0357 switch ch
0358         % Zoom in (X).
0359     case '.'
0360         
0361         canLen = max([(canEnd-canCen) (canCen-canStart)]);
0362         
0363         canStart = max(1,round(canCen-0.5*canLen));
0364         canEnd = min(dataLen,round(canCen+0.5*canLen));
0365         
0366         p = 1;
0367         
0368         % Zoom out (X).
0369     case ','
0370         
0371         canLen = max([(canEnd-canCen) (canCen-canStart)]);
0372         
0373         canStart = max(1,round(canCen-2*canLen));
0374         canEnd = min(dataLen,round(canCen+2*canLen));
0375         
0376         p = 1;
0377         
0378         % Zoom in (Y).
0379     case ']'
0380         
0381         canY = canY / 2;
0382 
0383         p = 1;
0384         
0385         % Zoom out (Y).
0386     case '['
0387         
0388         canY = canY * 2;
0389         
0390         p = 1;
0391         
0392         % Classify as RFI
0393     case 'r'
0394         verdict = 1;
0395         p = 0;
0396         
0397         % Classify as false positive
0398     case 'f'
0399         verdict = 0;
0400         p = 0;
0401         
0402         % Decline to classify
0403     case 'd'
0404         verdict = -1;
0405         p = 0;
0406         
0407         % Previous plot
0408     case 'p'
0409         verdict = -1;
0410         p = -1;
0411         
0412         % Apply the noise diode removal..
0413     case 'n'
0414         procCommand = 1;
0415         p = 1;
0416         
0417         % Apply the 1.2hz correction..
0418     case 'l'
0419         procCommand = 2;
0420         p = 1;
0421         
0422         % Display the help.
0423     case 'h'
0424         printhelp;
0425         p = 1;
0426         
0427         % Other commands not recognized.
0428     otherwise
0429         printcmr(ch);
0430         p = 1;
0431 end
0432 
0433 
0434 end
0435 
0436 
0437 %%%%%%%%%%%%%%%%%%%%%%%%%
0438 function printhelp()
0439 
0440 disp(' ')
0441 disp('******************************')
0442 disp('        RFI Plot Help         ')
0443 disp('******************************')
0444 disp(' ');
0445 disp('The user is expected to judiciously discriminate');
0446 disp('between probable RFI events and probable false ');
0447 disp('positives.  The user may choose to reserve judgement');
0448 disp('for some events, if the nature is ambiguous.');
0449 disp(' ')
0450 disp('h - this help menu')
0451 disp('r - classify as RFI')
0452 disp('f - classify as a false positive')
0453 disp('d - decline to classify')
0454 disp('n - remove the noise diode events from these data')
0455 disp('l - apply the load (1.2Hz) correction to these data')
0456 disp('. - zoom in on current plot  (X)')
0457 disp(', - zoom out on current plot (X)')
0458 disp('] - zoom in on current plot  (Y)')
0459 disp('[ - zoom out on current plot (Y)')
0460 disp(' ')
0461 disp('*******************WARNING*********************')
0462 disp('Do not exit plot by manually killing the window.')
0463 disp('This will cause the program to crash!')
0464 disp('***********************************************')
0465 disp(' ')
0466 
0467 
0468 end
0469 
0470 %%%%%%%%%%%%%%%%%%%%%%%%%
0471 function printcmr(ch)
0472 
0473 disp(' ')
0474 disp(['Oops!  Command not recognized: ' ch]);
0475 disp('Press h for help.');
0476 disp(' ')
0477 
0478 
0479 end
0480