+
−% 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>>
+
−%