sidelobes_cw_sm

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
  function sidelobes_cw(d)

    function ot analyze the sidelobes when viewing a cw source.  

  began as script by matthew stevenson, stephen muchovej modified to clean
  it up a bit.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0001 function sidelobes_cw(d)
0002 
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 %  function sidelobes_cw(d)
0006 %
0007 %    function ot analyze the sidelobes when viewing a cw source.
0008 %
0009 %  began as script by matthew stevenson, stephen muchovej modified to clean
0010 %  it up a bit.
0011 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0012 
0013 % first let's identify the data that was not taken with the noise source.
0014 [d_noise, d_off, indNoise] = cutNoise(d);
0015 indSource = bitsearch(d.array.frame.features, 0, 'only');
0016 indSource = interp1(d.array.frame.utc, indSource, d.antenna0.receiver.utc, 'closest');
0017 indSource = indSource & ~indNoise;
0018 indSource(isnan(indSource)) = 0;
0019 
0020 % next we take our data, and we have to figure out the period and offset of
0021 % the source
0022 data = d.antenna0.receiver.data(:,1);
0023 elevation = d.antenna0.servo.el;
0024 azimuth   = d.antenna0.servo.az;
0025 
0026 %d2 = read_arc('11-Mar-2010:22:07:17', '11-Mar-2010:22:57:40');
0027 %d2_length = length(d2.antenna0.receiver.utc);
0028 %data = d2.antenna0.receiver.data(9500:d2_length,1);
0029 %elevation = d2.antenna0.servo.el(9500:d2_length);
0030 %azimuth   = d2.antenna0.servo.az(9500:d2_length);
0031 
0032 % March 11, 2010.  CW without secondary baffle.
0033 %switch_freq = 0.70072;
0034 %switch_phase = 38;
0035 
0036 % March 11, 2010.  CW with secondary baffle.
0037 %switch_freq = 0.70088;
0038 %switch_phase = 61;
0039 
0040 % quick and dirty way of calculating the switch_freq and phase from data.
0041 % need to find when we're pointed closest to source
0042 f = find(d.antenna0.servo.az<150 & d.antenna0.servo.el <18);
0043 a = d.antenna0.receiver.data(f(1):f(1)+10000,1);
0044 
0045 da = diff(a);
0046 fa = find(abs(da)>30);
0047 dfa = diff(fa);
0048 fa(dfa>100 | dfa<30) = [];
0049 dfa(dfa>100) = [];
0050 dfa(dfa<30) = [];
0051 switch1 = mode(dfa);
0052 ff = find(dfa==switch1);
0053 if(mean(da(fa(ff)))>0)
0054   % switch to a positive higher value, given that the values in these
0055   % channels are negative, this means it's going from off to on
0056   off_time = floor(switch1/2)-5;
0057   offPeriod = switch1;
0058 else
0059   on_time  = floor(switch1/2)-5;
0060   onPeriod = switch1;
0061 end
0062 
0063 fa(dfa==mode(dfa)) = [];
0064 dfa(dfa==mode(dfa)) = [];
0065 switch2 = mode(dfa);
0066 ff = find(dfa==switch2);
0067 if(mean(da(fa(ff)))>0)
0068   % switch to a positive higher value, given that the values in these
0069   % channels are negative, this means it's going from off to on
0070   off_time = floor(switch2/2)-5;
0071   offPeriod = switch2;
0072 else
0073   on_time  = floor(switch2/2)-5;
0074   onPeriod = switch2;
0075 end
0076 
0077 switchFreq = 100./(switch1 + switch2+1);
0078 
0079 % now we calculate the offset.
0080 % I'll do this the quick and dirty way:  with a square wave of frequency
0081 % switch_freq of values 1 and 0, i'll just multiply the data and add it up.
0082 % the one that returns the greatest value should be the offset.
0083 
0084 % construct the square wave
0085 switchPeriod = round(100/switchFreq);
0086 singleSquare = ones(switchPeriod,1);
0087 singleSquare(1:offPeriod) = -1;
0088 numDataPoints = size(d.antenna0.receiver.data,1);
0089 totSquares = ceil(numDataPoints/switchPeriod)+1;
0090 
0091 onMidIndex = round(offPeriod + onPeriod/2);
0092 
0093 dataVals = abs(d.antenna0.receiver.data);
0094 
0095 sumVal = zeros(switchPeriod,6);
0096 for phase=0:switchPeriod-1
0097   % construct the square wave
0098   squareWave = repmat(singleSquare, [totSquares,1]);
0099   if(phase==0)
0100     squareWave = squareWave(1:numDataPoints);
0101   else
0102     squareWave(1:phase) = [];
0103     squareWave = squareWave(1:numDataPoints);
0104   end
0105   toSum = repmat(squareWave, [1 6]).*dataVals;
0106   sumVal(phase+1,:) = nansum(toSum);
0107 end
0108 switchPhase(1) = find(sumVal(:,1)==min(sumVal(:,1)));
0109 switchPhase(2) = find(sumVal(:,6)==min(sumVal(:,6)));
0110 
0111 
0112 
0113 %switch_phase = round((mean(switchPhase))*switchFreq);
0114 %switch_freq = switchFreq;
0115   
0116   
0117 %display('Calculated switch frequency and phase');
0118 %display(sprintf('frequency: %3.2f', switch_freq));
0119 %display(sprintf('phase:  %i', switch_phase));
0120 
0121 % I verified that this matches what matthew found by hand.
0122 % now back to matthew's stuff.
0123 fullPeriod = onPeriod + offPeriod+1;
0124 n_meas = floor(length(data) / fullPeriod);
0125 
0126 % sjcm: this doesn't work. don't know where the 25 and 50 are coming from.
0127 % it results in the negative values for the differences
0128 % On is defined as where (x-2)*0.70072/100 is equal to an integer.
0129 %x_on = floor(100*(1:n_meas) / switch_freq + (switch_phase - 25 / switch_freq));
0130 %x_off = x_on + floor(50 / switch_freq);
0131 
0132 % find the midpoints of each on/off period (approximate)
0133 x_on = round(mean(switchPhase)):fullPeriod:length(data);
0134 x_off = round(mean(switchPhase))- fullPeriod/2:fullPeriod:length(data);
0135 x_off(x_off<0) = [];
0136 
0137 
0138 diff_array = zeros(n_meas-2,1);
0139 elevation_array = zeros(n_meas-2,1);
0140 azimuth_array   = zeros(n_meas-2,1);
0141 on_source       = zeros(n_meas-2,1);
0142 
0143 
0144 on_time  = on_time-5;
0145 off_time = off_time-5;
0146 for i=1:n_meas-2
0147      
0148     on_mean = mean(data((x_on(i+1)-on_time):(x_on(i+1)+on_time)));
0149     off_mean_1 = mean(data(x_off(i):x_off(i)+off_time));
0150     off_mean_2 = mean(data(x_off(i+1)-off_time:x_off(i+1)));
0151     
0152 %    off_mean_1 = mean(data((x_off(i)+1):(x_off(i)+16)));
0153 %    off_mean_2 = mean(data((x_off(i+1)-16):(x_off(i+1)-1)));
0154    
0155     
0156     difference = -(on_mean - 0.5*off_mean_1 - 0.5*off_mean_2);
0157 
0158     if(difference < -10)
0159       display('problem');
0160       keyboard;
0161     end
0162     
0163     
0164     diff_array(i) = difference;
0165     
0166     
0167     elevation_array(i) = mean(elevation(x_on(i+1)-on_time:x_on(i+1)+on_time));
0168     azimuth_array(i)   = mean(azimuth(x_on(i+1)-on_time:x_on(i+1)+on_time));
0169     
0170     on_source(i) = mean(indSource(x_on(i+1)-on_time:x_on(i+1)+on_time))>0.8;
0171 end
0172 
0173 keyboard;
0174 on_source = on_source > 0;  % turn it into a logical array
0175 f = find(azimuth_array > 180);
0176 elevation_array(f) = 180 - elevation_array(f);
0177 
0178 plot(elevation_array(on_source), diff_array(on_source));
0179 xlabel('Elevation angle from horizon');
0180 ylabel('Backend Units');
0181 
0182