line source
+ − % 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
+ − %
+ − % <a href="matlab:utils.helper.displayMethodInfo('ao', 'spikecleaning')">Parameters Description</a>
+ − %
+ − % VERSION: $Id: spikecleaning.m,v 1.17 2011/04/08 08:56:16 hewitson Exp $
+ − %
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ −
+ − function varargout=spikecleaning(varargin)
+ −
+ − %%% Check if this is a call for parameters
+ − if utils.helper.isinfocall(varargin{:})
+ − varargout{1} = getInfo(varargin{3});
+ − return
+ − end
+ −
+ − if nargout == 0
+ − error('### cat cannot be used as a modifier. Please give an output variable.');
+ − end
+ −
+ − % Collect input variable names
+ − in_names = cell(size(varargin));
+ − for ii = 1:nargin,in_names{ii} = inputname(ii);end
+ −
+ − % Collect all AOs
+ − [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
+ − pli = utils.helper.collect_objects(varargin(:), 'plist', in_names);
+ −
+ − pls = parse(pli, getDefaultPlist());
+ −
+ − % initialise output array
+ − bo = [];
+ −
+ − % go through each input AO
+ − for i=1:numel(as)
+ − a = as(i);
+ − d = a.data;
+ −
+ − % check this is a time-series object
+ − if ~isa(d, 'tsdata')
+ − error(' ### temperature spike detection requires tsdata (time-series) inputs.')
+ − end
+ −
+ − %--- check input parameters
+ − kspike = find(pls, 'kspike'); % kspike*sigma definition
+ − method = find(pls, 'method'); % method of spike-values substitution
+ − pls.pset('gain', 1); % gain of the filter
+ − pls.pset('type', 'highpass'); % type of the filter
+ − fs = plist();
+ − fs.append('fs', d.fs);
+ − pls.combine(fs); % determination of the sampling frequency of the input AO
+ −
+ − % high-pass filtering data
+ − xfiltered = filtfilt(a, miir(pls));
+ −
+ − % standard deviation of the filtered data is calculated
+ − nxfiltered = find(abs(xfiltered) < kspike*std(xfiltered));
+ −
+ − xfiltered_2 = xfiltered.data.y(nxfiltered);
+ −
+ − std_xfiltered_2 = std(xfiltered_2);
+ −
+ − % spikes vector position is determined
+ − nspike = find(abs(xfiltered) > kspike*std_xfiltered_2);
+ −
+ − % substitution of spike values starts here
+ − xcleaned = a.data.y;
+ − for j=1:length(nspike)
+ − if nspike(j) <=2 % just in case a spike is detected in the 1st or 2nd sample
+ − xcleaned(nspike(j)) = mean(xcleaned(1:50));
+ − else
+ − if strcmp(method, 'random') % spike is substituted by a random value: N(0,std_xfiltered)
+ − xcleaned(nspike(j)) = xcleaned(nspike(j)-1) + randn(1)*std_xfiltered_2;
+ − elseif strcmp(method, 'mean') % spike is substituted by the mean if the two previous values
+ − xcleaned(nspike(j)) = (xcleaned(nspike(j)-1) + xcleaned(nspike(j)-2))/2;
+ − elseif strcmp(method, 'previous') % spike is substituted by the pervious value
+ − xcleaned(nspike(j)) = xcleaned(nspike(j)-1);
+ − end
+ − end
+ − end
+ −
+ − % create new output tsdata
+ − ts = tsdata(xcleaned, d.fs);
+ − ts.setYunits(d.yunits);
+ − ts.setXunits(d.xunits);
+ −
+ − % % create new output history
+ − % h = history(ALGONAME, VERSION, pls, a.hist);
+ − % h = set(h, 'invars', invars);
+ −
+ − % make output analysis object
+ − b = ao(ts);
+ − b.name = sprintf('spikecleaning(%s)', ao_invars{i});
+ − b.addHistory(getInfo('None'), pls, ao_invars(i), as(i).hist);
+ −
+ − % add to output array
+ − bo = [bo b];
+ −
+ − end
+ −
+ − % Set output
+ − if nargout == numel(bo)
+ − % List of outputs
+ − for ii = 1:numel(bo)
+ − varargout{ii} = bo(ii);
+ − end
+ − else
+ − % Single output
+ − varargout{1} = bo;
+ − end
+ −
+ − end
+ −
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ − % Local Functions %
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ −
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ − %
+ − % FUNCTION: getInfo
+ − %
+ − % DESCRIPTION: Get Info Object
+ − %
+ − % HISTORY: 11-07-07 M Hewitson
+ − % Creation.
+ − %
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ −
+ − function ii = getInfo(varargin)
+ − if nargin == 1 && strcmpi(varargin{1}, 'None')
+ − sets = {};
+ − pl = [];
+ − else
+ − sets = {'Default'};
+ − pl = getDefaultPlist();
+ − end
+ − % Build info object
+ − ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: spikecleaning.m,v 1.17 2011/04/08 08:56:16 hewitson Exp $', sets, pl);
+ − ii.setModifier(false);
+ − end
+ −
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ − %
+ − % FUNCTION: getDefaultPlist
+ − %
+ − % DESCRIPTION: Get Default Plist
+ − %
+ − % HISTORY: 11-07-07 M Hewitson
+ − % Creation.
+ − %
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ −
+ − function plout = getDefaultPlist()
+ − persistent pl;
+ − if exist('pl', 'var')==0 || isempty(pl)
+ − pl = buildplist();
+ − end
+ − plout = pl;
+ − end
+ −
+ − function pl = buildplist()
+ −
+ − pl = plist();
+ −
+ − % kspike
+ − p = param({'kspike', 'High values imply no correction of relative low amplitude spikes.'}, paramValue.DOUBLE_VALUE(3.3));
+ − pl.append(p);
+ −
+ − % fc
+ − p = param({'fc', 'Frequency cut-off of the IIR filter.'}, paramValue.DOUBLE_VALUE(0.025));
+ − pl.append(p);
+ −
+ − % Order
+ − p = param({'order', 'The order of the IIR filter.'}, paramValue.DOUBLE_VALUE(2));
+ − pl.append(p);
+ −
+ − % Ripple
+ − p = param({'ripple', 'Specify the pass/stop-band ripple for bandpass/bandreject filters'}, ...
+ − paramValue.DOUBLE_VALUE(0.5));
+ − pl.append(p);
+ −
+ − % Method
+ − p = param({'method', 'The method used to replace the spike value.'}, {1, {'random', 'mean'}, paramValue.SINGLE});
+ − pl.append(p);
+ −
+ − end
+ −
+ − % 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>>
+ − %