Home > RFI > RFI_Char.m

RFI_Char

PURPOSE ^

% This code will use two methods to characterise the RFI events in

SYNOPSIS ^

function RFI_Char()

DESCRIPTION ^

% This code will use two methods to characterise the RFI events in
% the cbass data timestream. It will save the processed data and
% plot it up.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 %% This code will use two methods to characterise the RFI events in
0002 %% the cbass data timestream. It will save the processed data and
0003 %% plot it up.
0004 
0005 function RFI_Char()
0006 cbass_startup
0007 
0008 start_time = '13-Jun-2014:09:04:00';
0009 %stop_time = '13-Jun-2014:23:12:25';
0010 stop_time = '13-Jun-2014:09:12:25';
0011 
0012 d = read_arcSouth_nojunk(start_time, stop_time, {'array.frame.utc double', ...
0013     'array.frame.features double', 'array.frame.received', ...
0014     'antenna0.roach1.utc double',  'antenna0.roach2.utc double', ...
0015     'antenna0.roach1.LLfreq double',  'antenna0.roach2.LLfreq double', ...
0016     'antenna0.roach1.RRfreq double',  'antenna0.roach2.RRfreq double', ...
0017     'antenna0.roach1.Qfreq double',  'antenna0.roach2.Qfreq double', ...
0018     'antenna0.roach1.Ufreq double',  'antenna0.roach2.Ufreq double'});
0019 
0020 LL1 = d.antenna0.roach1.LLfreq;
0021 LL2 = d.antenna0.roach2.LLfreq;
0022 RR1 = d.antenna0.roach1.RRfreq;
0023 RR2 = d.antenna0.roach2.RRfreq;
0024 
0025 Q1 = d.antenna0.roach1.Qfreq;
0026 Q2 = d.antenna0.roach2.Qfreq;
0027 U1 = d.antenna0.roach1.Ufreq;
0028 U2 = d.antenna0.roach2.Ufreq;
0029 
0030 
0031 % Determine the intensity:
0032 I1 = (LL1 + RR1)/2; 
0033 I2 = (LL2 + RR2)/2; 
0034 
0035 % Determine the polarisation:
0036 P1 = sqrt(Q1.*Q1 + U1.*U1);
0037 P2 = sqrt(Q2.*Q2 + U2.*U2);
0038 
0039 
0040 %% Method 1: RMS
0041 % Establish a bin size:
0042 
0043 % Note: a smaller binsize results in the top 10 rms hits all being part of
0044 % the same RFI event. Bigger bin size gets more events, but smoothes the
0045 % rms value across a larger sample value ie possibly including more RFI
0046 % events in one sample
0047 
0048 binsize = 100;
0049 nbins = lower(length(I1)/binsize);
0050 
0051 % Determine rms for each bin:
0052 for m = 1:nbins - 1
0053     rmsval_I1(m) = sqrt(var(I1((m-1)*100+1:m*100)));
0054     rmsval_I2(m) = sqrt(var(I2((m-1)*100+1:m*100)));
0055     rmsval_P1(m) = sqrt(var(P1((m-1)*100+1:m*100)));
0056     rmsval_P2(m) = sqrt(var(P2((m-1)*100+1:m*100)));
0057 end % for loop
0058 
0059 % Pick a brightness limit:
0060 brightRmsLimit_I1 = 2*median(rmsval_I1);
0061 brightRmsLimit_I2 = 4*median(rmsval_I2);
0062 brightRmsLimit_P1 = 4*median(rmsval_P1);
0063 brightRmsLimit_P2 = 4*median(rmsval_P2);
0064 
0065 % Pick RFI events above limit:
0066 f_I1 = find(rmsval_I1 > brightRmsLimit_I1)
0067 f_I2 = find(rmsval_I2 > brightRmsLimit_I2);
0068 f_P1 = find(rmsval_P1 > brightRmsLimit_P1);
0069 f_P2 = find(rmsval_P2 > brightRmsLimit_P2);
0070 
0071 for m=1:size(f_I1,2)
0072     if(f_I1(m)==1)
0073         tstart = d.antenna0.roach1.utc(1)-1/24/60/60;
0074     else
0075         tstart = d.antenna0.roach1.utc((f_I1(m)-1)*100)-1/24/60/60;
0076     end % if statement
0077     tstop  = d.antenna0.roach1.utc(f_I1(m)*100)+1/24/60/60;
0078     drfi_I1{m} = read_arcSouth(mjd2string(tstart), mjd2string(tstop));  
0079 end % for loop
0080 
0081 for m=1:size(f_I2,2)
0082     if(f_I2(m)==1)
0083         tstart = d.antenna0.roach1.utc(1)-1/24/60/60;
0084     else
0085         tstart = d.antenna0.roach1.utc((f_I2(m)-1)*100)-1/24/60/60;
0086     end % if statement
0087     tstop  = d.antenna0.roach1.utc(f_I2(m)*100)+1/24/60/60;
0088     drfi_I2{m} = read_arcSouth(mjd2string(tstart), mjd2string(tstop));  
0089 end % for loop
0090 
0091 for m=1:size(f_P1,2)
0092     if(f_P1(m)==1)
0093         tstart = d.antenna0.roach1.utc(1)-1/24/60/60;
0094     else
0095         tstart = d.antenna0.roach1.utc((f_P1(m)-1)*100)-1/24/60/60;
0096     end % if statement
0097     tstop  = d.antenna0.roach1.utc(f_P1(m)*100)+1/24/60/60;
0098     drfi_P1{m} = read_arcSouth(mjd2string(tstart), mjd2string(tstop));  
0099 end % for loop
0100 
0101 for m=1:size(f_P2,2)
0102     if(f_P2(m)==1)
0103         tstart = d.antenna0.roach1.utc(1)-1/24/60/60;
0104     else
0105         tstart = d.antenna0.roach1.utc((f_P2(m)-1)*100)-1/24/60/60;
0106     end % if statement
0107     tstop  = d.antenna0.roach1.utc(f_P2(m)*100)+1/24/60/60;
0108     drfi_P2{m} = read_arcSouth(mjd2string(tstart), mjd2string(tstop));  
0109 end % for loop
0110 
0111 plotfigs_I(drfi_I1)
0112 plotfigs_I(drfi_I2)
0113 %plotfigs_P(drfi_P1)
0114 %plotfigs_P(drfi_P2)
0115 
0116 end % RFI_Char
0117 
0118 
0119 %%%  plotting figures
0120 function plotfigs_I(drfi)
0121 %clc
0122 %close all
0123 %cbass_startup
0124 
0125 % first, let's just find the peak and see what its spectrum looks like
0126 figure()
0127 allspectra = [];
0128 Intensity_freq1 = [];
0129 Intensity1 = [];
0130 Intensity_freq2 = [];
0131 Intensity2 = [];
0132 
0133 for m=1:size(drfi,2)
0134 
0135     LL1freq = drfi{m}.antenna0.roach1.LLfreq;
0136     LL2freq = drfi{m}.antenna0.roach2.LLfreq;
0137     RR1freq = drfi{m}.antenna0.roach1.RRfreq;
0138     RR2freq = drfi{m}.antenna0.roach2.RRfreq;
0139     Intensity_freq1{m} = (LL1freq + RR1freq)/2;
0140     Intensity_freq2{m} = (LL2freq + RR2freq)/2;
0141     
0142     LL1 = drfi{m}.antenna0.roach1.LL;
0143     LL2 = drfi{m}.antenna0.roach2.LL;
0144     RR1 = drfi{m}.antenna0.roach1.RR;
0145     RR2 = drfi{m}.antenna0.roach2.RR;
0146     Intensity1{m} = (LL1 + RR1)/2;
0147     Intensity2{m} = (LL2 + RR2)/2;
0148     
0149     aa1 = find(Intensity_freq1{m} == max(Intensity_freq1{m}));
0150     aa2 = find(Intensity_freq2{m} == max(Intensity_freq2{m}));
0151     
0152     
0153     % to visualize them (First step)
0154     subplot(2,1,1); plot(Intensity1{m}(aa1,:));hold on
0155     title('Peak Intensity Values Across All Channels: Roach1')
0156     xlabel('Channels')
0157     ylabel('Intensity')
0158     subplot(2,1,2); plot(Intensity2{m}(aa2,:)); hold on
0159     title('Peak Intensity Values Across All Channels: Roach2')
0160     xlabel('Channels')
0161     ylabel('Intensity')
0162     allspectra(m,:) = [Intensity1{m}(aa1,:), Intensity2{m}(aa2,:)];
0163     %allspectra(m,:) = [drfi{m}.Intensity(aa,:)];
0164     %hold on
0165 end
0166 hold off
0167 %{
0168 keyboard;
0169 
0170 figure()
0171 for m=1:length(drfi)
0172   disp('New Event');
0173   f = find(drfi{m}.antenna0.roach1.LLfreq == max(drfi{m}.antenna0.roach1.LLfreq));
0174   minInt = f-200;
0175   maxInt = f+200;
0176   if(minInt < 1)
0177     minInt = 1;
0178   end
0179   if(maxInt > length(drfi{m}.antenna0.roach1.LLfreq))
0180     maxInt = length(drfi{m}.antenna0.roach1.LLfreq);
0181   end
0182   
0183   maxval = max(drfi{m}.antenna0.roach1.LL(:));
0184   cc = [drfi{m}.antenna0.roach1.LL(f,:) , drfi{m}.antenna0.roach2.LL(f,:)];
0185   for mm=minInt:maxInt
0186      bb = [drfi{m}.antenna0.roach1.LL(mm,:) , drfi{m}.antenna0.roach2.LL(mm,:)];
0187      plot(bb);
0188      ylim([0,maxval]);
0189      hold on
0190      plot(cc, 'r-');
0191      txt = sprintf('Event # %d', m);
0192      title(txt);
0193      pause(0.025);
0194      hold off     
0195   end
0196   pause(1)
0197   
0198 end
0199 
0200 keyboard;
0201 %}
0202 figure()
0203 plot(allspectra');
0204 title('Each RFI Event Across All Channels')
0205 ylabel('Intensity')
0206 xlabel('Channels')
0207 end
0208 
0209 
0210 
0211 
0212 
0213
Generated on Sun 14-Jun-2015 17:12:45 by m2html © 2005