LTPDA_spikecleaning detects and corrects possible spikes in analysis objects %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DESCRIPTION: LTPDA_ 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 = ltpda_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: ltpda_spikecleaning.html,v 1.7 2008/03/31 10:27:44 hewitson Exp $ HISTORY: 24-12-2007 J Sanjuan Creation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0001 function varargout=ltpda_spikecleaning(varargin) 0002 % LTPDA_spikecleaning detects and corrects possible spikes in analysis objects 0003 % 0004 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0005 % 0006 % DESCRIPTION: LTPDA_ 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 = ltpda_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: ltpda_spikecleaning.html,v 1.7 2008/03/31 10:27:44 hewitson Exp $ 0033 % 0034 % HISTORY: 24-12-2007 J Sanjuan 0035 % Creation 0036 % 0037 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0038 0039 ALGONAME = mfilename; 0040 VERSION = '$Id: ltpda_spikecleaning.html,v 1.7 2008/03/31 10:27:44 hewitson Exp $'; 0041 CATEGORY = 'Signal Processing'; 0042 0043 % Check if this is a call for parameters 0044 if nargin == 1 0045 in = char(varargin{1}); 0046 if strcmp(in, 'Params') 0047 varargout{1} = getDefaultPL(); 0048 return 0049 elseif strcmp(in, 'Version') 0050 varargout{1} = VERSION; 0051 return 0052 elseif strcmp(in, 'Category'); 0053 varargout{1} = CATEGORY; 0054 return 0055 end 0056 end 0057 0058 % collect input ao's, plist's and ao variable names 0059 in_names = {}; 0060 for ii = 1:nargin 0061 in_names{end+1} = inputname(ii); 0062 end 0063 [as, pls, invars] = collect_inputs(varargin, in_names); 0064 na = length(as); 0065 0066 pls = combine(pls, getDefaultPL); 0067 0068 % initialise output array 0069 bo = []; 0070 0071 % go through each input AO 0072 for i=1:na 0073 a = as(i); 0074 d = a.data; 0075 0076 % check this is a time-series object 0077 if ~isa(d, 'tsdata') 0078 error(' ### temperature spike detection requires tsdata (time-series) inputs.') 0079 end 0080 0081 %--- check input parameters 0082 kspike = find(pls, 'kspike'); % kspike*sigma definition 0083 method = find(pls, 'method'); % method of spike-values substitution 0084 pls = pset(pls, param('gain', 1)); % gain of the filter 0085 pls = pset(pls, param('type', 'highpass')); % type of the filter 0086 fs = plist(); 0087 fs = append(fs, param('fs', d.fs)); 0088 pls = combine(pls, fs); % determination of the sampling frequency of the input AO 0089 0090 % high-pass filtering data 0091 xfiltered = filtfilt(a, miir(pls)); 0092 0093 % standard deviation of the filtered data is calculated 0094 nxfiltered = find(abs(xfiltered) < kspike*std(xfiltered)); 0095 0096 xfiltered_2 = xfiltered.data.y(nxfiltered); 0097 0098 std_xfiltered_2 = std(xfiltered_2); 0099 0100 % spikes vector position is determined 0101 nspike = find(abs(xfiltered) > kspike*std_xfiltered_2); 0102 0103 % substitution of spike values starts here 0104 xcleaned = a.data.y; 0105 for j=1:length(nspike) 0106 if nspike(j) <=2 % just in case a spike is detected in the 1st or 2nd sample 0107 xcleaned(nspike(j)) = mean(xcleaned(1:50)); 0108 else 0109 if strcmp(method, 'random') % spike is substituted by a random value: N(0,std_xfiltered) 0110 xcleaned(nspike(j)) = xcleaned(nspike(j)-1) + randn(1)*std_xfiltered_2; 0111 elseif strcmp(method, 'mean') % spike is substituted by the mean if the two previous values 0112 xcleaned(nspike(j)) = (xcleaned(nspike(j)-1) + xcleaned(nspike(j)-2))/2; 0113 elseif strcmp(method, 'previous') % spike is substituted by the pervious value 0114 xcleaned(nspike(j)) = xcleaned(nspike(j)-1); 0115 end 0116 end 0117 end 0118 0119 % create new output tsdata 0120 ts = tsdata(xcleaned, d.fs); 0121 ts = set(ts, 'name', sprintf('cleaned data of %s', d.name)); 0122 ts = set(ts, 'yunits', d.yunits); 0123 ts = set(ts, 'xunits', d.xunits); 0124 0125 % create new output history 0126 h = history(ALGONAME, VERSION, pls, a.hist); 0127 h = set(h, 'invars', invars); 0128 0129 % make output analysis object 0130 b = ao(ts, h); 0131 0132 % set name 0133 if isempty(invars{i}) 0134 n = a.name; 0135 else 0136 n = invars{i}; 0137 end 0138 0139 b = setnh(b, 'name', sprintf('spikecleaning(%s)', n)); 0140 0141 % add to output array 0142 bo = [bo b]; 0143 0144 varargout{1} = bo; 0145 end 0146 0147 %------------------------------- 0148 % Get default parameters 0149 function plo = getDefaultPL(); 0150 0151 disp('* creating default plist...'); 0152 plo = plist(); 0153 plo = append(plo, param('kspike', 3.3)); 0154 plo = append(plo, param('fc', 0.025)); 0155 plo = append(plo, param('order', 2)); 0156 plo = append(plo, param('ripple', 0.5)); 0157 plo = append(plo, param('method', 'random')); 0158 disp('* done.');