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.6 2008/09/05 11:05:29 ingo Exp $ HISTORY: 24-12-2007 J Sanjuan Creation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0001 % spikecleaning detects and corrects possible spikes in analysis objects 0002 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0003 % 0004 % DESCRIPTION: SPIKECLEANING detects spikes in the temperature data and 0005 % replaces them by artificial values depending on the method chosen ('random', 0006 % 'mean', 'previous'). 0007 % Spikes are defined as singular samples with an (absolute) value 0008 % higher than kspike times the standard deviation of the high-pass 0009 % filtered (IIR filter) input AO. 0010 % 0011 % CALL: b = spikecleaning(a1, a2, ..., an, pl) 0012 % 0013 % INPUTS: aN - a list of analysis objects 0014 % pl - parameter list 0015 % 0016 % OUTPUTS: b - a list of analysis objects with "spike values" removed 0017 % and corrected 0018 % 0019 % PARAMETERES: 'kspike' - set the kspike value. High values imply 0020 % not correction of relative low amplitude spike 0021 % [default: 3.3] 0022 % 'method' - method used to replace the spike value: 'random, 0023 % 'mean', 'previous' [default:random] 0024 % 'fc' - frequency cut-off of the IIR filter [default: 0.025] 0025 % 'order' - order of the IIR filter [default: 2] 0026 % 'ripple' - specify pass/stop-band ripple for bandpass 0027 % and bandreject filters 0028 % <<default: 0.5>> 0029 % 0030 % VERSION: $Id: spikecleaning.m,v 1.6 2008/09/05 11:05:29 ingo Exp $ 0031 % 0032 % HISTORY: 24-12-2007 J Sanjuan 0033 % Creation 0034 % 0035 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0036 0037 function varargout=spikecleaning(varargin) 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(), pls, 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', '', utils.const.categories.sigproc, '$Id: spikecleaning.m,v 1.6 2008/09/05 11:05:29 ingo 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