first_cut

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

  function d = first_cut(d,sigma)

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

  Inputs:
    sigma - no. of sigma to cut on 
  sjcm & act

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

0001 function d = first_cut(d,sigma)
0002 
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 %  function d = first_cut(d,sigma)
0006 %
0007 %  rough function to find and flag glitches in the data, focusing on
0008 %  spurious points.  Needs a lot of improvement.  Right now
0009 %  might flag too much data.
0010 %
0011 %  Inputs:
0012 %    sigma - no. of sigma to cut on
0013 %  sjcm & act
0014 %
0015 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0016 
0017 % first we only want to flag on data is not the noise source.
0018 dnonoise = cutObs(d, 'noise', 'not');
0019 
0020 data = dnonoise.antenna0.receiver.data;
0021 
0022 % first we check to see if there spurious points.
0023 indSpurious = zeros(size(data,1),1);
0024 stop = 0;
0025 while(stop==0)
0026   for m=1:6
0027     a(:,m) = smooth(data(:,m));
0028   end
0029   diffVals = a-dnonoise.antenna0.receiver.data;
0030   stdVals  = std(diffVals)
0031   diffNorm = abs(diffVals)./repmat(stdVals, [size(diffVals,1) 1]);
0032  % plot(diffNorm)
0033  
0034  
0035   badVals     = diffNorm > sigma;
0036   newIndSpurious = sum(badVals,2)>1;
0037     
0038   indSpurious(newIndSpurious) = 1;
0039   data(newIndSpurious,:) = a(newIndSpurious,:);
0040 
0041   l = length(find(newIndSpurious));
0042   newOnes = length(find(indSpurious)) - l;
0043   
0044   if(newOnes<1)
0045     stop = 1;
0046   end
0047 end
0048 
0049 %indDeglitchFast = indSpurious | indDeriv;
0050 indDeglitchFast = indSpurious ;
0051 
0052 % convert to size of d, (instead of dnonoise)
0053 indFinal = logical(zeros(size(d.flags.fast,1),1));
0054 indFinal(~d.index.noise.fast) = indDeglitchFast;
0055 
0056 % convert to other timescales.
0057 flags = applyFlag(d, [], [], indFinal, 4);
0058 
0059 %Final sorting out of flags - check this is needed
0060 flags.bit.fast = bitor(flags.bit.fast, d.flags.bit.fast);
0061 flags.bit.medium = bitor(flags.bit.medium, d.flags.bit.medium);
0062 flags.bit.slow = bitor(flags.bit.slow, d.flags.bit.slow);
0063 display('Plotting data');
0064 gen = 3; % 1 if you want to plot for the web
0065 d = packd(d, flags, 'none', 'deglitch', gen, 'Deglitching Plots',[]);
0066 d = updateFlags(d);
0067 d = logcal(d, 'deglitch');
0068 return;
0069 
0070 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0071 function flags = applyFlag(d, flagRegTemp, flagServTemp, ...
0072     flagRxTemp, bitNum)
0073 
0074 if(~isempty(flagRegTemp))
0075   flagServTemp = repmat(flagRegTemp, [1 5]);
0076   flagServTemp = flagServTemp';
0077   flagServTemp = flagServTemp(:);
0078 
0079   flagRxTemp   = repmat(flagRegTemp, [1 100]);
0080   flagRxTemp   = flagRxTemp';
0081   flagRxTemp   = flagRxTemp(:);
0082   
0083 elseif(~isempty(flagServTemp))
0084   flagRegTemp = reshape(flagServTemp, [5, length(flagServTemp)/5]);
0085   flagRegTemp = mean(flagRegTemp) == 1;
0086   flagRegTemp = flagRegTemp';
0087 
0088   flagRxTemp   = repmat(flagServTemp, [1 20]);
0089   flagRxTemp   = flagRxTemp';
0090   flagRxTemp   = flagRxTemp(:);
0091 elseif(~isempty(flagRxTemp))
0092   flagRegTemp = reshape(flagRxTemp, [100, length(flagRxTemp)/100]);
0093   flagRegTemp = mean(flagRegTemp) == 1;
0094   flagRegTemp = flagRegTemp';
0095   
0096   flagServTemp = reshape(flagRxTemp, [20, length(flagRxTemp)/20]);
0097   flagServTemp = mean(flagServTemp) == 1;
0098   flagServTemp = flagServTemp';
0099   else
0100   error('No flags passed');
0101 end
0102 
0103 
0104 flags.bit.slow = uint8(2.^(bitNum)*double(flagRegTemp));
0105 
0106 flags.bit.medium = uint8(2.^(bitNum)*double(flagServTemp));
0107 
0108 flags.bit.fast = repmat(uint8(2.^(bitNum)*double(flagRxTemp)), [1 6]);
0109 
0110 display('DEGLITCHING FLAGGING REPORT');
0111 display(sprintf('%d of %d samples were flagged', length(find(flagRxTemp)), ...
0112     size(d.flags.fast,1)));
0113 
0114 
0115 
0116 return;