flagAlpha

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

  function flags = flagAlpha(d, flagParams)

  function to automatically flag results in the alpha correction
 
  sjcm

 ogk - Dec 2010: modified to work with updateAlphaDatabase_DD script.
 ogk - 21 Mar 2011: modified to take the alternative structures of
 d.correction.alpha.indices into account.

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

0001 function [flags,slopeLoad,slopeSky,interceptDiffSky] = flagAlpha(d, flagParams)
0002 
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 %  function flags = flagAlpha(d, flagParams)
0006 %
0007 %  function to automatically flag results in the alpha correction
0008 %
0009 %  sjcm
0010 %
0011 % ogk - Dec 2010: modified to work with updateAlphaDatabase_DD script.
0012 % ogk - 21 Mar 2011: modified to take the alternative structures of
0013 % d.correction.alpha.indices into account.
0014 %
0015 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0016 
0017 % what we want to flag are the noise diode events where we don't see any
0018 % improvement.
0019 
0020 %%%%%%%%%%%%%%%%%%%%%%%%%
0021 % OGK modified this, 21 Mar 2011
0022 % d.correction.alpha.indices = [offStPos onEndPos offEndPos onStPos];
0023 if size(d.correction.alpha.indices,2) == 4
0024     % Before 8-Sep
0025     si = d.correction.alpha.indices(:,1); % when the first off noise diode event starts
0026     ei = d.correction.alpha.indices(:,2); % when the on noise diode event ends
0027     esi = d.correction.alpha.indices(:,4); % when the on noise diode event starts
0028 else
0029     % After 8-Sep
0030     si = d.correction.alpha.indices(:,1); % when the first off noise diode event starts
0031     ei = d.correction.alpha.indices(:,3); % when the on noise diode event ends
0032     esi = d.correction.alpha.indices(:,6); % when the on noise diode event starts
0033 end
0034 %%%%%%%%%%%%%%%%%%%%%%%%%
0035 numEvents = size(d.correction.alpha.indices,1);
0036 
0037 maxInterceptDiff = flagParams(1);
0038 minInterceptDiff = flagParams(2);
0039 maxSlope = flagParams(3);
0040 
0041 % split up all the noise diode events
0042 for m=1:numEvents
0043   x{1,m} = d.antenna0.receiver.data(si(m):ei(m),[1,2]);
0044   x{2,m} = d.antenna0.receiver.data(si(m):ei(m),[9,10]);
0045 end
0046 
0047 % index 2 in x{1:2,m} should be flat
0048 % index 1 in x{1:2,m} should have a step in them
0049 
0050 % first we fit a line to the flat ones, and make sure there's no slope
0051 % greater than that specified in flagParams(2);
0052 
0053 for m=1:numEvents
0054   thisX = [1:1:size(x{1,m},1)]';
0055   thisY = x{1,m}(:,2);
0056   thisX = thisX/100;
0057 
0058   % Channel 1 slope and intercept
0059   [slopeVal(m,1), interceptVal(m,1)] = linfit(thisX, thisY);
0060   
0061   % Channel 2 slope and intercept
0062   thisY = x{2,m}(:,2);
0063   [slopeVal(m,2), interceptVal(m,2)] = linfit(thisX, thisY);  
0064 end
0065 % find those slopeVals that are greater than the parameters
0066 slopeVal = abs(slopeVal);
0067 slopeLoad = slopeVal;
0068 badSlope1 = slopeVal>maxSlope;
0069 badSlope1 = sum(badSlope1,2)>0; % flag if either of the channels has a bad slope val
0070 
0071 % next we fit a line to both sides of the step ones, and
0072 % A) make sure there's no slope greater than that specified in flagParams(2)
0073 % B) the y-intercept of the horizontal lines are no more than flagParams(1)
0074 % apart
0075 clear slopeVal;
0076 clear interceptVal;
0077 for m=1:numEvents
0078   for nchan=1:2
0079     numSamples = size(x{nchan,m},1);
0080     midSample  = esi(m)-si(m); % the official start of the on noise diode event
0081     slack = 50; % don't trust it
0082     
0083     thisX1 = [1:(midSample-slack)]';
0084     thisY1 = x{nchan,m}(thisX1,1);
0085     thisX1 = thisX1/100;
0086     
0087     thisX2 = [midSample+slack:numSamples]';
0088     thisY2 = x{nchan,m}(thisX2,1);
0089     thisX2 = thisX2/100;
0090 
0091     [slopeA intA] = linfit(thisX1, thisY1);
0092     [slopeB intB] = linfit(thisX2, thisY2);
0093     
0094     thisSlope = [slopeA slopeB];
0095     thisInt   = [intA intB];
0096     
0097     slopeVal(m,nchan,:) = thisSlope;
0098     interceptVal(m,nchan,:) = thisInt;
0099   end
0100 end
0101 
0102 % first we check for criteria A
0103 slopeVal = abs(slopeVal);
0104 slopeSky = slopeVal;
0105 badSlope2 = slopeVal>maxSlope;
0106 badSlope2 = sum(sum(badSlope2,3),2)>0;
0107 
0108 % next we check for criteria B
0109 interceptDiff = interceptVal(:,:,1) - interceptVal(:,:,2);
0110 interceptDiff = abs(interceptDiff);
0111 interceptDiffSky = interceptDiff;
0112 badIntercept  = sum(sum(interceptDiff>maxInterceptDiff | ...
0113     interceptDiff<minInterceptDiff,3),2)>0;
0114 
0115 % now we combine it all to find out which events should be flagged
0116 badEvents = badSlope1 | badSlope2 | badIntercept;
0117 
0118 flags = badEvents;
0119 
0120 return;
0121 
0122 % % plotting for debug purposes
0123 % for m=1:numEvents
0124 %   subplot(2,1,1)
0125 %   plot(x{1, m});
0126 %   subplot(2,1,2);
0127 %   plot(x{2, m});
0128 %
0129 %   if(flags(m))
0130 %     gtitle('FLAGGED');
0131 %     keyboard;
0132 %   else
0133 %     gtitle('CLEAR');
0134 %   end
0135 %   pause(1);
0136 % end
0137 %
0138 %
0139 %