characterizeRfi_P

0001 function characterizeRfi_P(starttime, stoptime, savefilename);
0002 cbass_startup
0003 
0004 % This code bins the timestream and determines the RMS of each bin. This in
0005 % turn allows us to look at the top RFI events across the entire
0006 % timestream. If we didn't bin, the top RFI events could all be part of the
0007 % same overall event.
0008 close all
0009 clc
0010 
0011 %starttime = '08-jun-2014:21:27:00';
0012 %stoptime  = '09-jun-2014:10:52:00';
0013 %stoptime  = '08-jun-2014:21:52:00';
0014 
0015 
0016 d = read_arcSouth_nojunk(starttime, stoptime, {'array.frame.utc double', ...
0017     'array.frame.features double', 'array.frame.received', ...
0018     'antenna0.roach1.utc double',  'antenna0.roach2.utc double', ...
0019     'antenna0.roach1.LLfreq double',  'antenna0.roach2.LLfreq double', ...
0020     'antenna0.roach1.RRfreq double',  'antenna0.roach2.RRfreq double', ...
0021     'antenna0.roach1.Qfreq double',  'antenna0.roach2.Qfreq double', ...
0022     'antenna0.roach1.Ufreq double',  'antenna0.roach2.Ufreq double'});
0023 
0024 % try removing the steps
0025 % This didn't work, and caused issues on cyclops, so I removed it.
0026 
0027 % Determine the intensity:
0028 % I = sqrt(((Q1+Q2)/2)^2 + ((U1+U2)/2)^2)
0029 Q1 = d.antenna0.roach1.Qfreq;
0030 Q2 = d.antenna0.roach2.Qfreq;
0031 U1 = d.antenna0.roach1.Ufreq;
0032 U2 = d.antenna0.roach2.Ufreq;
0033 
0034 Q = (Q1 + Q2)./2;
0035 U = (U1 + U2)./2;
0036 
0037 Polarization = sqrt(Q.*Q + U.*U);
0038 
0039 
0040 % Establish bin size:
0041 
0042 % Note: a smaller binsize results in the top 10 rms hits all being part of
0043 % the same RFI event. Bigger bin size gets more events, but smoothes the
0044 % rms value across a larger sample value ie possibly including more RFI
0045 % events in one sample
0046 
0047 binsize = 100;
0048 nbins = lower(length(Polarization)/binsize);
0049 
0050 % Determine rms for each bin:
0051 
0052 for m = 1:nbins - 1
0053     rmsval(m) = sqrt(var(Polarization((m-1)*100+1:m*100)));
0054 end
0055 
0056 % Sort the rms values into highest to lowest:
0057 
0058 [rmssort I] = sort(rmsval,'descend');
0059 
0060 
0061 top_ten = rmssort(1,1:10); % top 10 RFI rms values
0062 top_ten_bin = I(1,1:10); % bin numbers relating to top 10 RFI rms values
0063 
0064 
0065 
0066 %pick the RMS which are 5 sigma above the median rms
0067 brightRmsLimit = 4*median(rmsval);
0068 
0069 f = find(rmsval > brightRmsLimit);
0070 % with those N points, we want to read in the full data, and look at a
0071 % spectrum for each
0072 
0073 for m=1:size(f,2)
0074     if(f(m)==1)
0075         tstart = d.antenna0.roach1.utc(1)-1/24/60/60;
0076     else
0077         tstart = d.antenna0.roach1.utc((f(m)-1)*100)-1/24/60/60;
0078     end   
0079     
0080     tstop  = d.antenna0.roach1.utc(f(m)*100)+1/24/60/60;
0081     drfi{m} = read_arcSouth(mjd2string(tstart), mjd2string(tstop));
0082    
0083 end
0084 
0085 % Save the variables:
0086 drfisave = drfi;
0087 numEvents = length(drfi);
0088 if(numEvents >100)
0089     disp('splitting it up'); % splitting up into sets of 100 events (else too big)
0090     numFiles = ceil(numEvents/100);
0091     for m=1:numFiles
0092         if(m==numFiles)
0093             dd = drfisave((m-1)*100+1:numEvents);
0094         else
0095             dd = drfisave((m-1)*100+1:(m-1)*100);
0096         end
0097         newsavefile = strcat(savefilename, '_', num2str(m));
0098         drfi = dd;
0099         
0100         txt = sprintf('save %s drfi', newsavefile);        
0101         eval(txt);
0102     end
0103 else
0104     txt = sprintf('save %s drfi', savefilename);
0105     eval(txt);
0106 end
0107 
0108 
0109 end % function