


spikecleaning detects and corrects possible spikes in analysis objects
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
DESCRIPTION: SPIKECLEANING detects spikes in the temperature data and
replaces them by artificial values depending on the method chosen ('random',
'mean', 'previous').
Spikes are defined as singular samples with an (absolute) value
higher than kspike times the standard deviation of the high-pass
filtered (IIR filter) input AO.
CALL: b = spikecleaning(a1, a2, ..., an, pl)
INPUTS: aN - a list of analysis objects
pl - parameter list
OUTPUTS: b - a list of analysis objects with "spike values" removed
and corrected
PARAMETERES: 'kspike' - set the kspike value. High values imply
not correction of relative low amplitude spike
[default: 3.3]
'method' - method used to replace the spike value: 'random,
'mean', 'previous' [default:random]
'fc' - frequency cut-off of the IIR filter [default: 0.025]
'order' - order of the IIR filter [default: 2]
'ripple' - specify pass/stop-band ripple for bandpass
and bandreject filters
<<default: 0.5>>
VERSION: $Id: spikecleaning.m,v 1.3 2008/08/08 15:31:38 josep Exp $
HISTORY: 24-12-2007 J Sanjuan
Creation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


0001 function varargout=spikecleaning(varargin) 0002 % spikecleaning detects and corrects possible spikes in analysis objects 0003 % 0004 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0005 % 0006 % DESCRIPTION: SPIKECLEANING detects spikes in the temperature data and 0007 % replaces them by artificial values depending on the method chosen ('random', 0008 % 'mean', 'previous'). 0009 % Spikes are defined as singular samples with an (absolute) value 0010 % higher than kspike times the standard deviation of the high-pass 0011 % filtered (IIR filter) input AO. 0012 % 0013 % CALL: b = spikecleaning(a1, a2, ..., an, pl) 0014 % 0015 % INPUTS: aN - a list of analysis objects 0016 % pl - parameter list 0017 % 0018 % OUTPUTS: b - a list of analysis objects with "spike values" removed 0019 % and corrected 0020 % 0021 % PARAMETERES: 'kspike' - set the kspike value. High values imply 0022 % not correction of relative low amplitude spike 0023 % [default: 3.3] 0024 % 'method' - method used to replace the spike value: 'random, 0025 % 'mean', 'previous' [default:random] 0026 % 'fc' - frequency cut-off of the IIR filter [default: 0.025] 0027 % 'order' - order of the IIR filter [default: 2] 0028 % 'ripple' - specify pass/stop-band ripple for bandpass 0029 % and bandreject filters 0030 % <<default: 0.5>> 0031 % 0032 % VERSION: $Id: spikecleaning.m,v 1.3 2008/08/08 15:31:38 josep Exp $ 0033 % 0034 % HISTORY: 24-12-2007 J Sanjuan 0035 % Creation 0036 % 0037 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0038 0039 %%% Check if this is a call for parameters 0040 if utils.helper.isinfocall(varargin{:}) 0041 varargout{1} = getInfo(varargin{3}); 0042 return 0043 end 0044 0045 if nargout == 0 0046 error('### cat cannot be used as a modifier. Please give an output variable.'); 0047 end 0048 0049 % Collect input variable names 0050 in_names = cell(size(varargin)); 0051 for ii = 1:nargin,in_names{ii} = inputname(ii);end 0052 0053 % Collect all AOs 0054 [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names); 0055 pli = utils.helper.collect_objects(varargin(:), 'plist', in_names); 0056 0057 pls = combine(pli, getDefaultPlist()); 0058 0059 % initialise output array 0060 bo = []; 0061 0062 % go through each input AO 0063 for i=1:numel(as) 0064 a = as(i); 0065 d = a.data; 0066 0067 % check this is a time-series object 0068 if ~isa(d, 'tsdata') 0069 error(' ### temperature spike detection requires tsdata (time-series) inputs.') 0070 end 0071 0072 %--- check input parameters 0073 kspike = find(pls, 'kspike'); % kspike*sigma definition 0074 method = find(pls, 'method'); % method of spike-values substitution 0075 pls.pset('gain', 1); % gain of the filter 0076 pls.pset('type', 'highpass'); % type of the filter 0077 fs = plist(); 0078 fs.append('fs', d.fs); 0079 pls.combine(fs); % determination of the sampling frequency of the input AO 0080 0081 % high-pass filtering data 0082 xfiltered = filtfilt(a, miir(pls)); 0083 0084 % standard deviation of the filtered data is calculated 0085 nxfiltered = find(abs(xfiltered) < kspike*std(xfiltered)); 0086 0087 xfiltered_2 = xfiltered.data.y(nxfiltered); 0088 0089 std_xfiltered_2 = std(xfiltered_2); 0090 0091 % spikes vector position is determined 0092 nspike = find(abs(xfiltered) > kspike*std_xfiltered_2); 0093 0094 % substitution of spike values starts here 0095 xcleaned = a.data.y; 0096 for j=1:length(nspike) 0097 if nspike(j) <=2 % just in case a spike is detected in the 1st or 2nd sample 0098 xcleaned(nspike(j)) = mean(xcleaned(1:50)); 0099 else 0100 if strcmp(method, 'random') % spike is substituted by a random value: N(0,std_xfiltered) 0101 xcleaned(nspike(j)) = xcleaned(nspike(j)-1) + randn(1)*std_xfiltered_2; 0102 elseif strcmp(method, 'mean') % spike is substituted by the mean if the two previous values 0103 xcleaned(nspike(j)) = (xcleaned(nspike(j)-1) + xcleaned(nspike(j)-2))/2; 0104 elseif strcmp(method, 'previous') % spike is substituted by the pervious value 0105 xcleaned(nspike(j)) = xcleaned(nspike(j)-1); 0106 end 0107 end 0108 end 0109 0110 % create new output tsdata 0111 ts = tsdata(xcleaned, d.fs); 0112 ts.setYunits(d.yunits); 0113 ts.setXunits(d.xunits); 0114 0115 % % create new output history 0116 % h = history(ALGONAME, VERSION, pls, a.hist); 0117 % h = set(h, 'invars', invars); 0118 0119 % make output analysis object 0120 b = ao(ts); 0121 b.setName(sprintf('spikecleaning(%s)', ao_invars{i}), 'internal'); 0122 %warning('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%'); 0123 %warning('%%%%%%%%%%%% Uncomment the following line (add Hsitory) %%%%%'); 0124 %warning('%%%%%%%%%%%% Uncomment the following line (add Hsitory) %%%%%'); 0125 %warning('%%%%%%%%%%%% Uncomment the following line (add Hsitory) %%%%%'); 0126 b.addHistory(getInfo(), pli, ao_invars(i), as(i).hist); 0127 0128 % add to output array 0129 bo = [bo b]; 0130 0131 varargout{1} = bo; 0132 end 0133 end 0134 0135 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0136 % Local Functions % 0137 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0138 0139 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0140 % 0141 % FUNCTION: getInfo 0142 % 0143 % DESCRIPTION: Get Info Object 0144 % 0145 % HISTORY: 11-07-07 M Hewitson 0146 % Creation. 0147 % 0148 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0149 0150 function ii = getInfo(varargin) 0151 if nargin == 1 && strcmpi(varargin{1}, 'None') 0152 sets = {}; 0153 pl = []; 0154 else 0155 sets = {'Default'}; 0156 pl = getDefaultPlist(); 0157 end 0158 % Build info object 0159 ii = minfo(mfilename, 'ao', '', 'Signal Processing', '$Id: spikecleaning.m,v 1.3 2008/08/08 15:31:38 josep Exp $', sets, pl); 0160 end 0161 0162 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0163 % 0164 % FUNCTION: getDefaultPlist 0165 % 0166 % DESCRIPTION: Get Default Plist 0167 % 0168 % HISTORY: 11-07-07 M Hewitson 0169 % Creation. 0170 % 0171 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0172 0173 function plo = getDefaultPlist() 0174 0175 plo = plist( ... 0176 'kspike', 3.3, ... 0177 'fc', 0.025, ... 0178 'order', 2, ... 0179 'ripple', 0.5, ... 0180 'method', 'random'); 0181 end 0182