find_rfi

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

  function full_rfi = find_rfi(d, parm, skymodel)

    Seriusly, what do you think this does?


    sjcm


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

0001 function [full_rfi rfiflags sim] = find_rfi(d, parm,thedate_day_portion,skymodel)
0002 
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 %  function full_rfi = find_rfi(d, parm, skymodel)
0006 %
0007 %    Seriusly, what do you think this does?
0008 %
0009 %
0010 %    sjcm
0011 %
0012 %
0013 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0014 
0015 MAXSCAN = 4.35;  %Hz, deg/s
0016 DEFSCALEP= 6.0;   % default scaling to make Polsky match pol data.
0017 DEFSCALEI= 1.05;   % default scaling to make Isky match I data.
0018 DEFOPACITY= 0.008; % default opacity.
0019 
0020 if(nargin>1)
0021   FASTSIG = parm.rfi.fast;
0022   POLSIG  = parm.rfi.polarization;
0023   INTSIG  = parm.rfi.intensity;
0024 else
0025   FASTSIG = 4.0;
0026   POLSIG  = 3.0;
0027   INTSIG  = 5.0;
0028 end
0029 
0030 
0031 [home,installeddir]=where_am_i();
0032 
0033 
0034 scaleTerms_filename = ['scaleTerms_',thedate_day_portion,'.txt'];
0035 
0036 
0037 
0038 % Stephen''s RFI flagging scheme:
0039 % I:  fast RFI
0040 %   This step looks for RFI events which are happening at faster than our
0041 %   scan rate.  It basically:
0042 %    Ia.  applies a filter
0043 %    Ib.  searches for peaks in the hi-pass data which are N-sigma greater
0044 %    than the rms (calculated from the hi-pass data itself)
0045 %    Ic.  checks whether these events are Noise Diode events
0046 %    Id.  checks the lengths of the events, if less than 2, it is a "glitch"
0047 % II:  slow RFI
0048 %   This step then takes the lo-passed data, and compares that to a sky
0049 %   model.  In essence it:
0050 %    IIa.  Given a sky model, calculates what we expect to see in I/Q/U
0051 %    IIb.  Looks at the difference between our low-passed data and this
0052 %    model
0053 %    IIc.  Calculates a metric based on the discrepancy and the noise levels
0054 %    in both the data and the model.
0055 %    IId.  Is it Sun/Moon
0056 % III.  Metric for RFI
0057 %   This is where we use all the information available to us.  There are
0058 %   many schemes we can try where we can weight polarization offsets more,
0059 %   etc etc.  Offhand, I would say that in step I, anything with discrepancy
0060 %   greater than snr of 3 which lasts longer than 3 samples is probably an RFI.
0061 
0062 % setup variables
0063 glitches    = logical(zeros(size(d.antenna0.receiver.data,1),1));
0064 fast_rfi    = logical(zeros(size(d.antenna0.receiver.data,1),1));
0065 full_rfi    = logical(zeros(size(d.antenna0.receiver.data,1),1));
0066 pol_rfi     = logical(zeros(size(d.antenna0.receiver.data,1),1));
0067 int_rfi     = logical(zeros(size(d.antenna0.receiver.data,1),1));
0068 wings_rfi   = logical(zeros(size(d.antenna0.receiver.data,1),1));
0069 noise_trans = logical(zeros(size(d.antenna0.receiver.data,1),1));
0070 
0071 % I: FAST RFI
0072 % Ia. Apply filter:
0073 d = both_filter(d, MAXSCAN, MAXSCAN+0.2, 'cos', 0);
0074 
0075 % Ib. calculate rms of hi
0076 hi = d.antenna0.receiver.hi;
0077 hirms = rms(hi);
0078 hisig = hi./(repmat(hirms, [size(hi,1) 1]));
0079 hihi  = hisig > FASTSIG;
0080 hihi = sum(hihi,2) > 0;
0081 [si ei] = findStartStop(hihi, 10);
0082 
0083 % Ic. Is it noise diode on/off
0084 [on off] = findStartStop(d.index.noise.fast, 5);
0085 noiseTransitions = [on off];
0086 for m=1:length(si)
0087   aa = abs(noiseTransitions - si(m));
0088   %f  = find( aa == min(aa));
0089   %if (length(f) >1)
0090   % keyboard;
0091   %end
0092   %offset_from_ND(m) = aa(f);
0093 
0094   offset_from_ND(m) = min(aa); 
0095 end
0096 ND_events = offset_from_ND < 20;
0097 si(ND_events>0) = [];
0098 ei(ND_events>0) = [];
0099 % actually flag the transitions as well for future reference
0100 for m=1:length(noiseTransitions)
0101   ss = noiseTransitions(m)-4;
0102   ee = noiseTransitions(m)+4;
0103   if(ss<1)
0104     ss = 1;
0105   end
0106   if(ee > length(noise_trans))
0107     ee = length(noise_trans);
0108   end
0109   noise_trans(ss:ee) = 1;
0110 end
0111 
0112 %Id.  check on lengths:
0113 li = ei - si;
0114 fg = find(li < 4);
0115 ffr= find(li >=4);
0116 for m=1:length(fg)
0117   glitches(si(fg(m)):ei(fg(m))) = 1;
0118 end
0119 for m=1:length(ffr)
0120   fast_rfi(si(ffr(m)):ei(ffr(m))) = 1;
0121 end
0122 full_rfi = glitches | fast_rfi;
0123 
0124 % SLOW RFI
0125 % first we get the values form the sky.
0126 [Ivals, Pvals, Ivar, Pvar] = decomposeSky(d, 1024, 2);
0127 
0128 sim.I.val = Ivals;
0129 sim.I.var = Ivar;
0130 sim.P.val = Pvals;
0131 sim.P.var = Pvar;
0132 
0133 if(isempty(Pvar))
0134   Pvar = zeros(size(Pvals));
0135 end
0136 if(isempty(Ivar))
0137   Ivar = zeros(size(Ivals));
0138 end
0139 
0140 % next we need to remove a baseline, and rescale by noise diode temp to get
0141 % the units to match
0142 try
0143   [nV nI] =  getNoiseDiodeTemps(d);
0144   meanNoise = mean(nI);
0145 catch exception
0146   display('find_rfi:: WARNING:  NO NOISE DIODE TEMPS PRESENT');
0147   display('find_rfi:: USING VALUES OF 3K');
0148   meanNoise = [3 3];
0149 end
0150 
0151 % next we rescale the map values
0152 Idata = d.antenna0.receiver.data(:,1)*meanNoise(1) + d.antenna0.receiver.data(:,8)*meanNoise(2);
0153 Pdata = sqrt(prod(meanNoise))*(pyth(d.antenna0.receiver.data(:,6), d.antenna0.receiver.data(:,7)));
0154 
0155 % before we remove the background, let's take out the effect of the changing
0156 % elevation
0157 Iel = (mean(d.array.weather.airTemperature) + 273.15 + ...
0158     2.73)*(DEFOPACITY./sind(d.antenna0.servo.el));
0159 Iel = Iel - min(Iel);
0160 Idata = Idata - Iel;
0161 
0162 % now we remove a baseline
0163 t = (d.antenna0.receiver.utc - d.antenna0.receiver.utc(1))*24*60;
0164 offI = remove_background(t, Idata, 'StepSize', 0.2, 'WindowSize', 2);
0165 offP = remove_background(t, Pdata, 'StepSize', 0.2, 'WindowSize', 2);
0166 
0167 
0168 % next, remove the effect of the noise diode.
0169 % in I, that's pretty straightforward, we just remove the mean
0170 for m=1:size(d.correction.alpha.indices,1)
0171   ssi = d.correction.alpha.indices(m,4)+1;
0172   eei = d.correction.alpha.indices(m,2);
0173   offI(ssi:eei) = offI(ssi:eei) - mean(meanNoise);
0174   rmsI(m) = sqrt(var(offI(ssi:eei)));
0175   
0176   offs1i = d.correction.alpha.indices2(m,1);
0177   offe1i = d.correction.alpha.indices2(m,5);
0178   offs2i = d.correction.alpha.indices2(m,6);
0179   offe2i = d.correction.alpha.indices2(m,3);
0180   
0181   mPoff = mean(offP([offs1i:offe1i offs2i:offe2i]));
0182   rmsP(m) = sqrt(var(offP([offs1i:offe1i offs2i:offe2i])));
0183   mPon  = mean(offP(ssi:eei));
0184   offP(ssi-1:eei) = offP(ssi-1:eei) - (mPon - mPoff);
0185 end
0186 
0187 offI = smooth(offI);
0188 offP = smooth(offP);
0189 
0190 rmsI = median(rmsI);
0191 rmsP = median(rmsP);
0192 
0193 % next we need to correct for a gain factor in the polarization -- but this
0194 % is pretty hard when there's no signal, so we should pick it when it's on a
0195 % source.
0196 unsource = unique(d.antenna0.tracker.source);
0197 hasTauA  = strmatch('TauA', unsource);
0198 if(~isempty(hasTauA))
0199   ftau = strmatch('TauA', d.antenna0.tracker.source);
0200   indTau = zeros(size(d.antenna0.tracker.source));
0201   indTau(ftau) = 1;
0202   
0203   % we want on-source data which does not have the noise diode on
0204   indOff = d.antenna0.tracker.horiz_off(:,1)==0;
0205   indScan= d.antenna0.tracker.scan_off(:,1)==0 & ...
0206       d.antenna0.tracker.scan_off(:,2) == 0;
0207   indSrc = d.antenna0.tracker.offSource == 0;
0208   indNoise = ~d.index.noise.slow;
0209   fullInd = indTau & indOff & indScan & indSrc & indNoise;
0210   [dc, indReg, indServ, indRx] = framecut(d, logical(fullInd));
0211 
0212   gainFac = offP(indRx)./Pvals(indRx);
0213   scaleFacP = median(gainFac);
0214 
0215   gainFac = offI(indRx)./Ivals(indRx);
0216   scaleFacI = median(gainFac);
0217 
0218   % we need to save this for posterity
0219   Pvals = Pvals*scaleFacP;
0220   Ivals = Ivals*scaleFacI;
0221   Pvar  = Pvar*scaleFacP^2;
0222   Ivar  = Ivar*scaleFacI^2;
0223   iselephant = 1;
0224   if(iselephant)
0225     scaleTerms = [d.array.frame.utc(1) scaleFacI scaleFacP];
0226     %allScale = load([home,'/',installeddir,'/constants/scaleTerms.txt']);
0227     %allScale = [allScale; scaleTerms];
0228     save([home,'/',installeddir,'/constants/scaleTerms.txt'], 'scaleTerms', '-ascii', '-double','-append');
0229     save([home,'/',installeddir,'/constants/',scaleTerms_filename], 'scaleTerms', '-ascii', '-double','-append');
0230 
0231 %    system('cvs commit -m "update" constants/scaleTerms.txt');
0232   end
0233   
0234 else
0235   Pvals = Pvals*DEFSCALEP;
0236   Ivals = Ivals*DEFSCALEI;
0237   Pvar  = Pvar*DEFSCALEP^2;
0238   Ivar  = Ivar*DEFSCALEI^2;
0239 end
0240 
0241 % and now, we calculate the significance of each point.
0242 rmsOffP = sqrt(Pvar + rmsP.^2);
0243 rmsOffI = sqrt(Ivar + rmsI.^2);
0244 
0245 Psig  = abs(offP - Pvals)./rmsOffP;
0246 Isig  = abs(offI - Ivals)./rmsOffI;
0247 
0248 % to account for skydips (the fact that the telescope doesn't behave as we
0249 % expect), we should change the significance.  I know, it's horrible.
0250 inddip  = d.index.skydip.fast;
0251 Isig(inddip) = Isig(inddip)/5;
0252 Psig(inddip) = Psig(inddip)/2;
0253 
0254 % lastly, we identify the RFI:
0255 [si ei] = findStartStop(Psig>POLSIG & ~full_rfi, 10);
0256 if (numel(si)>0)
0257     for m=1:length(si)
0258         pol_rfi(si(m):ei(m)) = 1;
0259     end
0260 end
0261 full_rfi = full_rfi | pol_rfi;
0262 
0263 [si ei] = findStartStop(Isig>INTSIG & ~full_rfi, 10);
0264 if (numel(si)>0)
0265     for m=1:length(si)
0266         int_rfi(si(m):ei(m)) = 1;
0267     end
0268 end
0269 full_rfi = full_rfi | int_rfi;
0270 
0271 % lastly, we find the RFI that are longer than a particular time, and then
0272 % add a few extra samples before and after, just for good measure.
0273 [si ei] = findStartStop(full_rfi);
0274 f = find(ei - si > 10);
0275 si(si<51) = 51;
0276 ei(ei>(length(full_rfi)-50)) = length(full_rfi)-50;
0277 for m=1:length(f)
0278   wings_rfi(si(f(m))-50:ei(f(m))+50) = 1;
0279 end
0280 full_rfi = full_rfi | wings_rfi;
0281 
0282 rfiflags.data.fast            = ones(size(full_rfi));
0283 
0284 [gs gm gf] = determineInd(wings_rfi);
0285 rfiflags.wings.slow        = logical(gs);
0286 rfiflags.wings.medium      = logical(gm);
0287 rfiflags.wings.fast        = logical(gf);
0288 
0289 [gs gm gf] = determineInd(noise_trans);
0290 rfiflags.noiseTransition.slow        = logical(gs);
0291 rfiflags.noiseTransition.medium      = logical(gm);
0292 rfiflags.noiseTransition.fast        = logical(gf);
0293 
0294 [gs gm gf] = determineInd(pol_rfi);
0295 rfiflags.polarization.slow        = logical(gs);
0296 rfiflags.polarization.medium      = logical(gm);
0297 rfiflags.polarization.fast        = logical(gf);
0298 
0299 [gs gm gf] = determineInd(int_rfi);
0300 rfiflags.intensity.slow        = logical(gs);
0301 rfiflags.intensity.medium      = logical(gm);
0302 rfiflags.intensity.fast        = logical(gf);
0303 
0304 [gs gm gf] = determineInd(fast_rfi);
0305 rfiflags.fast.slow        = logical(gs);
0306 rfiflags.fast.medium      = logical(gm);
0307 rfiflags.fast.fast        = logical(gf);
0308 
0309 
0310 [gs gm gf] = determineInd(glitches);
0311 rfiflags.glitches.slow        = logical(gs);
0312 rfiflags.glitches.medium      = logical(gm);
0313 rfiflags.glitches.fast        = logical(gf);
0314 
0315 %
0316 %keyboard;
0317 %
0318 %
0319 %
0320 %Idata = d.antenna0.receiver.data(:,1)*meanNoise(1) + d.antenna0.receiver.data(:,8)*meanNoise(2);
0321 %figure(1);
0322 %setwinsize(gcf, 1200, 800);
0323 %clf
0324 %subplot(2,1,1);
0325 %plot(t, Idata);
0326 %hold on
0327 %plot(t(full_rfi), Idata(full_rfi), 'y.');
0328 %plot(t(glitches), Idata(glitches), 'r.');
0329 %plot(t(fast_rfi), Idata(fast_rfi), 'g.');
0330 %plot(t(int_rfi),  Idata(int_rfi), 'c.');
0331 %plot(t(pol_rfi),  Idata(pol_rfi), 'k.');
0332 %legend('Data', 'wings', 'glitches', 'Fast', 'Intens', 'Pol');
0333 %xlim([0 30]);
0334 %title('Total Intensity Data');
0335 %subplot(2,1,2);
0336 %plot(t, Isig);
0337 %hold on
0338 %plot(t(full_rfi), Isig(full_rfi), 'y.');
0339 %plot(t(glitches), Isig(glitches), 'r.');
0340 %plot(t(fast_rfi), Isig(fast_rfi), 'g.');
0341 %plot(t(int_rfi),  Isig(int_rfi), 'c.');
0342 %plot(t(pol_rfi),  Isig(pol_rfi), 'k.');
0343 %legend('Data', 'wings', 'glitches', 'Fast', 'Intens', 'Pol');
0344 %xlim([0 30]);
0345 %title('Total Intensity Metric');
0346 %
0347 %
0348 %figure(2);
0349 %setwinsize(gcf, 1200, 800);
0350 %clf
0351 %subplot(2,1,1);
0352 %plot(t, Pdata);
0353 %hold on
0354 %plot(t(full_rfi), Pdata(full_rfi), 'y.');
0355 %plot(t(glitches), Pdata(glitches), 'r.');
0356 %plot(t(fast_rfi), Pdata(fast_rfi), 'g.');
0357 %plot(t(int_rfi),  Pdata(int_rfi), 'c.');
0358 %plot(t(pol_rfi),  Pdata(pol_rfi), 'k.');
0359 %legend('Data', 'wings', 'glitches', 'Fast', 'Intens', 'Pol');
0360 %xlim([0 30]);
0361 %title('Polarization data');
0362 %subplot(2,1,2);
0363 %plot(t, Psig);
0364 %hold on
0365 %plot(t(full_rfi), Psig(full_rfi), 'y.');
0366 %plot(t(glitches), Psig(glitches), 'r.');
0367 %plot(t(fast_rfi), Psig(fast_rfi), 'g.');
0368 %plot(t(int_rfi),  Psig(int_rfi), 'c.');
0369 %plot(t(pol_rfi),  Psig(pol_rfi), 'k.');
0370 %legend('Data', 'wings', 'glitches', 'Fast', 'Intens', 'Pol');
0371 %xlim([0 30]);
0372 %title('Polarization metric');
0373 %
0374 %figure(3);
0375 %clf
0376 %fData = Pdata./Idata*100;
0377 %plot(t, fData);
0378 %hold on
0379 %plot(t(full_rfi), fData(full_rfi), 'y.');
0380 %plot(t(glitches), fData(glitches), 'r.');
0381 %plot(t(fast_rfi), fData(fast_rfi), 'g.');
0382 %plot(t(int_rfi),  fData(int_rfi), 'c.');
0383 %plot(t(pol_rfi),  fData(pol_rfi), 'k.');
0384 %legend('Data', 'wings', 'glitches', 'Fast', 'Intens', 'Pol');
0385 %xlim([0 30]);
0386 
0387 
0388 return;
0389 
0390 
0391 
0392 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0393 function [indReg, indServ, indRx] = determineInd(ind)
0394 
0395 indRx   = ind;
0396                                 
0397 indReg = reshape(indRx, [100 length(indRx)/100]);
0398 indReg  = mean(indReg)~=0;
0399 indReg  = indReg';
0400 
0401 indServ = reshape(indRx, [20 length(indRx)/20]);
0402 indServ = mean(indServ)~=0;
0403 indServ = indServ';
0404 
0405 return;
0406         
0407 
0408 
0409   
0410 
0411 
0412