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.m,v 1.2 2008/01/29 10:50:26 josep 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.m,v 1.2 2008/01/29 10:50:26 josep Exp $ 0033 % 0034 % HISTORY: 24-12-2007 J Sanjuan 0035 % Creation 0036 % 0037 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0038 0039 ALGONAME = mfilename; 0040 VERSION = '$Id: ltpda_spikecleaning.m,v 1.2 2008/01/29 10:50:26 josep Exp $'; 0041 CATEGORY = 'SIGPROC'; 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 % Capture input variables names 0059 invars = {}; 0060 as = []; 0061 ps = []; 0062 for j=1:nargin 0063 invars = [invars cellstr(inputname(j))]; 0064 if isa(varargin{j}, 'ao') 0065 as = [as varargin{j}]; 0066 end 0067 if isa(varargin{j}, 'plist') 0068 ps = [ps varargin{j}]; 0069 end 0070 end 0071 0072 % check plist 0073 if isempty(ps) 0074 pl = getDefaultPL(); 0075 else 0076 pl = combine(ps, getDefaultPL); 0077 end 0078 0079 % initialise output array 0080 bo = []; 0081 0082 % go through each input AO 0083 na = length(as); 0084 for i=1:na 0085 a = as(i); 0086 d = a.data; 0087 0088 % check this is a time-series object 0089 if ~isa(d, 'tsdata') 0090 error(' ### temperature spike detection requires tsdata (time-series) inputs.') 0091 end 0092 0093 %--- check input parameters 0094 kspike = find(pl, 'kspike'); % kspike*sigma definition 0095 method = find(pl, 'method'); % method of spike-values substitution 0096 pl = pset(pl, param('gain', 1)); % gain of the filter 0097 pl = pset(pl, param('type', 'highpass')); % type of the filter 0098 fs = plist(); 0099 fs = append(fs, param('fs', d.fs)); 0100 pl = combine(pl, fs); % determination of the sampling frequency of the input AO 0101 0102 % High-pass filtering data 0103 xfiltered = filtfilt(a, miir(pl)); 0104 0105 % standard deviation of the filtered data is calculated 0106 nxfiltered = find(abs(xfiltered) < kspike*std(xfiltered)); 0107 0108 xfiltered_2 = xfiltered.data.y(nxfiltered); 0109 0110 std_xfiltered_2 = std(xfiltered_2); 0111 0112 % spikes vector position is determined 0113 nspike = find(abs(xfiltered) > kspike*std_xfiltered_2); 0114 0115 % substitution of spike values starts here 0116 xcleaned = a.data.y; 0117 for j=1:length(nspike) 0118 if nspike(j) <=2 % just in case a spike is detected in the 1st or 2nd sample 0119 xcleaned(nspike(j)) = mean(xcleaned(1:50)); 0120 else 0121 if strcmp(method, 'random') % spike is substituted by a random value: N(0,std_xfiltered) 0122 xcleaned(nspike(j)) = xcleaned(nspike(j)-1) + randn(1)*std_xfiltered_2; 0123 elseif strcmp(method, 'mean') % spike is substituted by the mean if the two previous values 0124 xcleaned(nspike(j)) = (xcleaned(nspike(j)-1) + xcleaned(nspike(j)-2))/2; 0125 elseif strcmp(method, 'previous') % spike is substituted by the pervious value 0126 xcleaned(nspike(j)) = xcleaned(nspike(j)-1); 0127 end 0128 end 0129 end 0130 0131 % create new output tsdata 0132 ts = tsdata(xcleaned, d.fs); 0133 ts = set(ts, 'name', sprintf('cleaned data of %s', d.name)); 0134 ts = set(ts, 'yunits', d.yunits); 0135 ts = set(ts, 'xunits', d.xunits); 0136 0137 % create new output history 0138 h = history(ALGONAME, VERSION, pl, a.hist); 0139 h = set(h, 'invars', invars); 0140 0141 % make output analysis object 0142 b = ao(ts, h); 0143 0144 % set name 0145 if isempty(invars{i}) 0146 n = a.name; 0147 else 0148 n = invars{i}; 0149 end 0150 0151 b = setnh(b, 'name', sprintf('spikecleaning(%s)', n)); 0152 0153 % add to output array 0154 bo = [bo b]; 0155 0156 varargout{1} = bo; 0157 end 0158 0159 %------------------------------- 0160 % Get default parameters 0161 function plo = getDefaultPL(); 0162 0163 disp('* creating default plist...'); 0164 plo = plist(); 0165 plo = append(plo, param('kspike', 3.3)); 0166 plo = append(plo, param('fc', 0.025)); 0167 plo = append(plo, param('order', 2)); 0168 plo = append(plo, param('ripple', 0.5)); 0169 plo = append(plo, param('method', 'random')); 0170 disp('* done.');