Mercurial > hg > ltpda
view m-toolbox/classes/@ao/spikecleaning.m @ 44:409a22968d5e default
Add unit tests
author | Daniele Nicolodi <nicolodi@science.unitn.it> |
---|---|
date | Tue, 06 Dec 2011 18:42:11 +0100 |
parents | f0afece42f48 |
children |
line wrap: on
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>> %