deglitch

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

  function d = deglitch(d, params)

  rough function to find and flag glitches in the data, focusing on
  spurious points and on RFIs.  Needs a lot of improvement.  Right now
  might flag too much data.

  params = [threshld on channel 1, threshold on channel 2];

  sjcm & act

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

0001 function flags = deglitch(d, params)
0002 
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 %  function d = deglitch(d, params)
0006 %
0007 %  rough function to find and flag glitches in the data, focusing on
0008 %  spurious points and on RFIs.  Needs a lot of improvement.  Right now
0009 %  might flag too much data.
0010 %
0011 %  params = [threshld on channel 1, threshold on channel 2];
0012 %
0013 %  sjcm & act
0014 %
0015 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0016 
0017 % first we only want to flag on data is not the noise source.
0018 % swap the orders so we don't reproduce the data structure that we don't
0019 % need.
0020 data = d.antenna0.receiver.data;
0021 % only want the ones that are not the noise source
0022 [s e] = findStartStop(d.index.noise.fast);
0023 % need a buffer for noise going on and off
0024 s = s-5;
0025 e = e+5;
0026 s(s<1) = 1;
0027 e(e>size(data,1)) = size(data,1);
0028 indNoise = zeros(size(data,1),1);
0029 for m=1:length(s)
0030   indNoise(s(m):e(m)) = 1;
0031 end
0032 indNoNoise = ~indNoise;
0033 
0034 data = data(indNoNoise,:);
0035 
0036 
0037 %dnonoise = cutObs(d, 'noise', 'not');
0038 %data = dnonoise.antenna0.receiver.data;
0039 
0040 % first we check to see if there spurious points.
0041 indSpurious = logical(zeros(size(data,1),1));
0042 stop = 0;
0043 while(stop==0)
0044   for m=1:size(data,2)
0045     a(:,m) = smooth(data(:,m));
0046   end
0047   diffVals = a-data;
0048   stdVals  = std(diffVals);
0049   diffNorm = abs(diffVals)./repmat(stdVals, [size(diffVals,1) 1]);
0050   badVals     = diffNorm>5;
0051   newIndSpurious = sum(badVals,2)>0;
0052 
0053   indSpurious(newIndSpurious) = 1;
0054   data(newIndSpurious,:) = a(newIndSpurious,:);
0055 
0056   l = length(find(newIndSpurious));
0057   newOnes = length(find(indSpurious)) - l;
0058   
0059   if(newOnes<1)
0060     stop = 1;
0061   end
0062 end
0063 clear diffNorm;
0064 clear diffVals;
0065 clear a;
0066 
0067 
0068 
0069 % next we look for values where the data spikes up over a small amount of
0070 % time.
0071 da = deriv(data);
0072 indDeriv = abs(da(:,1))>params(1) | abs(da(:,6)>params(2));
0073 
0074 
0075 % next we check for the boxCarRejection on the "pseudo-good" data
0076 % NOT DOING UNTIL WE FIND THE RIGHT PARAMETERS
0077 %dnonoise.antenna0.receiver.data(indSpurious | indDeriv,:) = nan;
0078 %[rfiFlag] = boxRejection(dnonoise,5,10);
0079 %t = (dnonoise.antenna0.receiver.utc - dnonoise.antenna0.receiver.utc(1))*24*60*60;
0080 %plot(t, log(dnonoise.antenna0.receiver.data(:,1)));
0081 %hold on
0082 %plot(t(rfiFlag), log(dnonoise.antenna0.receiver.data(rfiFlag,1)), '.r');
0083 %hold off
0084 
0085 % combine the ones we flagged on
0086 indDeglitchFast = indSpurious | indDeriv;
0087 
0088 % convert to size of d, (instead of dnonoise)
0089 indFinal = logical(zeros(size(d.flags.fast,1),1));
0090 indFinal(indNoNoise) = indDeglitchFast;
0091 
0092 %% flag the transitional data as well.
0093 %indFinal = indFinal | d.index.transition.fast;
0094 
0095 % convert to other timescales.
0096 [d, flags] = setNewFlag(d, indFinal, 'deglitch');
0097 
0098 clear d;
0099 % now we have to set the proper flags in the flag field
0100 return;
0101 
0102 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0103 function flags = applyFlag(d, flagRegTemp, flagServTemp, ...
0104     flagRxTemp, bitNum)
0105 
0106 if(~isempty(flagRegTemp))
0107   flagServTemp = repmat(flagRegTemp, [1 5]);
0108   flagServTemp = flagServTemp';
0109   flagServTemp = flagServTemp(:);
0110 
0111   flagRxTemp   = repmat(flagRegTemp, [1 100]);
0112   flagRxTemp   = flagRxTemp';
0113   flagRxTemp   = flagRxTemp(:);
0114   
0115 elseif(~isempty(flagServTemp))
0116   flagRegTemp = reshape(flagServTemp, [5, length(flagServTemp)/5]);
0117   flagRegTemp = mean(flagRegTemp) == 1;
0118   flagRegTemp = flagRegTemp';
0119 
0120   flagRxTemp   = repmat(flagServTemp, [1 20]);
0121   flagRxTemp   = flagRxTemp';
0122   flagRxTemp   = flagRxTemp(:);
0123 elseif(~isempty(flagRxTemp))
0124   flagRegTemp = reshape(flagRxTemp, [100, length(flagRxTemp)/100]);
0125   flagRegTemp = mean(flagRegTemp) == 1;
0126   flagRegTemp = flagRegTemp';
0127   
0128   flagServTemp = reshape(flagRxTemp, [20, length(flagRxTemp)/20]);
0129   flagServTemp = mean(flagServTemp) == 1;
0130   flagServTemp = flagServTemp';
0131   else
0132   error('No flags passed');
0133 end
0134 
0135 
0136 flags.bit.slow = uint8(2.^(bitNum)*double(flagRegTemp));
0137 
0138 flags.bit.medium = uint8(2.^(bitNum)*double(flagServTemp));
0139 
0140 flags.bit.fast = repmat(uint8(2.^(bitNum)*double(flagRxTemp)), [1 3]);
0141 
0142 display('DEGLITCHING FLAGGING REPORT');
0143 display(sprintf('%d of %d samples were flagged', length(find(flagRxTemp)), ...
0144     size(d.flags.fast,1)));
0145 
0146 
0147 return;