noiseRemove

%%%%%%%%%%%%%%%%%%%%%%
 Correct the noise events to match the rest of the times series.
 Optional parameter switchMode determines whether to apply this to
 switchData or data channels.  0 is data channels, any other value is
 switchData.
%%%%%%%%%%%%%%%%%%%%%%

0001 function d = noiseRemove(d,switchMode)
0002 
0003 %%%%%%%%%%%%%%%%%%%%%%%
0004 % Correct the noise events to match the rest of the times series.
0005 % Optional parameter switchMode determines whether to apply this to
0006 % switchData or data channels.  0 is data channels, any other value is
0007 % switchData.
0008 %%%%%%%%%%%%%%%%%%%%%%%
0009 
0010 if nargin < 2
0011     switchMode = 0;
0012 end
0013 
0014 % How many samples on either side of noise event to be used for noise event
0015 % subtraction.
0016 noiseOff = 200;
0017 
0018 
0019 dLength = length(d.antenna0.receiver.utc);
0020 
0021 
0022 % Identify the noise-on times.
0023 if isfield(d,'index')
0024     noiseOn = d.index.noise.fast;
0025 else
0026     noiseOn = bitand(d.antenna0.receiver.flags,4) > 0;
0027 end
0028 
0029 
0030 % Begin by identifying the steps:
0031 dataStateDiff = diff(noiseOn);
0032 dataStateEdges = cat(1,0,find(dataStateDiff ~= 0),dLength);
0033 nStates = length(dataStateEdges) - 1;
0034 
0035 
0036 % Now we loop over the states.
0037 for k=1:nStates
0038     % check if this is a noise_on state.
0039     if noiseOn(dataStateEdges(k)+1) == 0
0040         continue;
0041     end
0042 
0043     % If the noise diode is on for the entire data vector, then leave the
0044     % loop.
0045     if (dataStateEdges(k) == 0) && (dataStateEdges(k+1) == dLength)
0046         continue;
0047     end
0048     
0049     
0050     % Calculate the noise_on median:
0051     if switchMode == 0
0052         noiseOnMedian = median( ...
0053             d.antenna0.receiver.data(dataStateEdges(k)+1:dataStateEdges(k+1),:));
0054     else
0055         noiseOnMedian = median( ...
0056             d.antenna0.receiver.switchData(dataStateEdges(k)+1:dataStateEdges(k+1),:));
0057     end
0058     
0059     % Calculate the noise_off median:
0060     if dataStateEdges(k) < noiseOff
0061         noiseOffA = dataStateEdges(k+1) + 1;
0062         noiseOffB = dataStateEdges(k+1) + noiseOff;        
0063         noiseOffC = dataStateEdges(k+1) + noiseOff + 1;
0064         noiseOffD = dataStateEdges(k+1) + 2 * noiseOff;        
0065     elseif dLength - dataStateEdges(k+1) <= noiseOff
0066         noiseOffA = dataStateEdges(k) - 2* noiseOff + 1;
0067         noiseOffB = dataStateEdges(k) - noiseOff;
0068         noiseOffC = dataStateEdges(k) - noiseOff + 1;
0069         noiseOffD = dataStateEdges(k);
0070     else
0071         noiseOffA = dataStateEdges(k) - noiseOff + 1;
0072         noiseOffB = dataStateEdges(k);
0073         noiseOffC = dataStateEdges(k+1) + 1;
0074         noiseOffD = dataStateEdges(k+1) + noiseOff;
0075     end
0076 
0077 
0078     if switchMode == 0
0079         noiseOffMedian = ...
0080             median(cat(1, ...
0081             d.antenna0.receiver.data(noiseOffA:noiseOffB,:), ...
0082             d.antenna0.receiver.data(noiseOffC:noiseOffD,:)));
0083     else
0084         noiseOffMedian = ...
0085             median(cat(1, ...
0086             d.antenna0.receiver.switchData(noiseOffA:noiseOffB,:), ...
0087             d.antenna0.receiver.switchData(noiseOffC:noiseOffD,:)));
0088     end
0089 
0090     
0091     % Subtract off the difference of the ON and OFF medians.
0092     if switchMode == 0
0093         d.antenna0.receiver.data(dataStateEdges(k)+1:dataStateEdges(k+1),:) = ...
0094             d.antenna0.receiver.data(dataStateEdges(k)+1:dataStateEdges(k+1),:) - ...
0095             repmat(noiseOnMedian - noiseOffMedian, ...
0096             dataStateEdges(k+1) - dataStateEdges(k),1);
0097     else
0098         d.antenna0.receiver.switchData(dataStateEdges(k)+1:dataStateEdges(k+1),:) = ...
0099             d.antenna0.receiver.switchData(dataStateEdges(k)+1:dataStateEdges(k+1),:) - ...
0100             repmat(noiseOnMedian - noiseOffMedian, ...
0101             dataStateEdges(k+1) - dataStateEdges(k),1);
0102     end
0103     
0104     
0105     % Also blank the transition regions to the local median.
0106 %    if dataStateEdges(k) > 0
0107     minInd = max(dataStateEdges(k)-11,1);
0108     maxInd = min(dataStateEdges(k)+14,dLength);
0109     if switchMode == 0
0110         d.antenna0.receiver.data(minInd:maxInd,:) = ...
0111             repmat(noiseOffMedian,maxInd-minInd+1,1);
0112     else
0113         d.antenna0.receiver.switchData(minInd:maxInd,:) = ...
0114             repmat(noiseOffMedian,maxInd-minInd+1,1);
0115     end
0116 %    end
0117 %    if dataStateEdges(k+1) < dLength
0118     minInd = max(dataStateEdges(k+1)-12,1);
0119     maxInd = min(dataStateEdges(k+1)+13,dLength);
0120     if switchMode == 0
0121         d.antenna0.receiver.data(minInd:maxInd,:) = ...
0122             repmat(noiseOffMedian,maxInd-minInd+1,1);
0123     else
0124         d.antenna0.receiver.switchData(minInd:maxInd,:) = ...
0125             repmat(noiseOffMedian,maxInd-minInd+1,1);
0126     end
0127 %        d.antenna0.receiver.data(dataStateEdges(k+1)+1,:) = noiseOffMedian;
0128 %    end
0129     
0130 end
0131 
0132 end