Home > classes > @ao > spikecleaning.m

spikecleaning

PURPOSE ^

spikecleaning detects and corrects possible spikes in analysis objects

SYNOPSIS ^

function varargout=spikecleaning(varargin)

DESCRIPTION ^

 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

 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: spikecleaning.m,v 1.6 2008/09/05 11:05:29 ingo Exp $

 HISTORY:     24-12-2007 J Sanjuan
                 Creation

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 % spikecleaning detects and corrects possible spikes in analysis objects
0002 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0003 %
0004 % DESCRIPTION: SPIKECLEANING detects spikes in the temperature data and
0005 %               replaces them by artificial values depending on the method chosen ('random',
0006 %               'mean', 'previous').
0007 %               Spikes are defined as singular samples with an (absolute) value
0008 %               higher than kspike times the standard deviation of the high-pass
0009 %               filtered (IIR filter) input AO.
0010 %
0011 % CALL:        b = spikecleaning(a1, a2, ..., an, pl)
0012 %
0013 % INPUTS:    aN - a list of analysis objects
0014 %               pl - parameter list
0015 %
0016 % OUTPUTS:     b - a list of analysis objects with "spike values" removed
0017 %              and corrected
0018 %
0019 % PARAMETERES: 'kspike' - set the kspike value. High values imply
0020 %                         not correction of relative low amplitude spike
0021 %                         [default: 3.3]
0022 %               'method' - method used to replace the spike value: 'random,
0023 %                         'mean', 'previous' [default:random]
0024 %               'fc' - frequency cut-off of the IIR filter [default: 0.025]
0025 %               'order' - order of the IIR filter [default: 2]
0026 %               'ripple' - specify pass/stop-band ripple for bandpass
0027 %                         and bandreject filters
0028 %                         <<default: 0.5>>
0029 %
0030 % VERSION:     $Id: spikecleaning.m,v 1.6 2008/09/05 11:05:29 ingo Exp $
0031 %
0032 % HISTORY:     24-12-2007 J Sanjuan
0033 %                 Creation
0034 %
0035 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0036 
0037 function varargout=spikecleaning(varargin)
0038 
0039   %%% Check if this is a call for parameters
0040   if utils.helper.isinfocall(varargin{:})
0041     varargout{1} = getInfo(varargin{3});
0042     return
0043   end
0044 
0045   if nargout == 0
0046     error('### cat cannot be used as a modifier. Please give an output variable.');
0047   end
0048 
0049   % Collect input variable names
0050   in_names = cell(size(varargin));
0051   for ii = 1:nargin,in_names{ii} = inputname(ii);end
0052 
0053   % Collect all AOs
0054   [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
0055   pli             = utils.helper.collect_objects(varargin(:), 'plist', in_names);
0056 
0057   pls = combine(pli, getDefaultPlist());
0058 
0059   % initialise output array
0060   bo = [];
0061 
0062   % go through each input AO
0063   for i=1:numel(as)
0064     a = as(i);
0065     d = a.data;
0066 
0067     % check this is a time-series object
0068     if ~isa(d, 'tsdata')
0069       error(' ### temperature spike detection requires tsdata (time-series) inputs.')
0070     end
0071 
0072     %--- check input parameters
0073     kspike = find(pls, 'kspike'); % kspike*sigma definition
0074     method = find(pls, 'method'); % method of spike-values substitution
0075     pls.pset('gain', 1);          % gain of the filter
0076     pls.pset('type', 'highpass'); % type of the filter
0077     fs = plist();
0078     fs.append('fs', d.fs);
0079     pls.combine(fs);              % determination of the sampling frequency of the input AO
0080 
0081     % high-pass filtering data
0082     xfiltered = filtfilt(a, miir(pls));
0083 
0084     % standard deviation of the filtered data is calculated
0085     nxfiltered = find(abs(xfiltered) < kspike*std(xfiltered));
0086 
0087     xfiltered_2 = xfiltered.data.y(nxfiltered);
0088 
0089     std_xfiltered_2 = std(xfiltered_2);
0090 
0091     % spikes vector position is determined
0092     nspike = find(abs(xfiltered) > kspike*std_xfiltered_2);
0093 
0094     % substitution of spike values starts here
0095     xcleaned = a.data.y;
0096     for j=1:length(nspike)
0097       if nspike(j) <=2 % just in case a spike is detected in the 1st or 2nd sample
0098         xcleaned(nspike(j)) = mean(xcleaned(1:50));
0099       else
0100         if strcmp(method, 'random') % spike is substituted by a random value: N(0,std_xfiltered)
0101           xcleaned(nspike(j)) = xcleaned(nspike(j)-1) + randn(1)*std_xfiltered_2;
0102         elseif strcmp(method, 'mean') % spike is substituted by the mean if the two previous values
0103           xcleaned(nspike(j)) = (xcleaned(nspike(j)-1) + xcleaned(nspike(j)-2))/2;
0104         elseif strcmp(method, 'previous') % spike is substituted by the pervious value
0105           xcleaned(nspike(j)) = xcleaned(nspike(j)-1);
0106         end
0107       end
0108     end
0109 
0110     % create new output tsdata
0111     ts = tsdata(xcleaned, d.fs);
0112     ts.setYunits(d.yunits);
0113     ts.setXunits(d.xunits);
0114 
0115     % % create new output history
0116     % h = history(ALGONAME, VERSION, pls, a.hist);
0117     % h = set(h, 'invars', invars);
0118 
0119     % make output analysis object
0120     b = ao(ts);
0121     b.setName(sprintf('spikecleaning(%s)', ao_invars{i}), 'internal');
0122     %warning('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%');
0123     %warning('%%%%%%%%%%%%  Uncomment the following line (add Hsitory)  %%%%%');
0124     %warning('%%%%%%%%%%%%  Uncomment the following line (add Hsitory)  %%%%%');
0125     %warning('%%%%%%%%%%%%  Uncomment the following line (add Hsitory)  %%%%%');
0126     b.addHistory(getInfo(), pls, ao_invars(i), as(i).hist);
0127 
0128     % add to output array
0129     bo = [bo b];
0130 
0131     varargout{1} = bo;
0132   end
0133 end
0134 
0135 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0136 %                               Local Functions                               %
0137 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0138 
0139 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0140 %
0141 % FUNCTION:    getInfo
0142 %
0143 % DESCRIPTION: Get Info Object
0144 %
0145 % HISTORY:     11-07-07 M Hewitson
0146 %                Creation.
0147 %
0148 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0149 
0150 function ii = getInfo(varargin)
0151   if nargin == 1 && strcmpi(varargin{1}, 'None')
0152     sets = {};
0153     pl   = [];
0154   else
0155     sets = {'Default'};
0156     pl   = getDefaultPlist();
0157   end
0158   % Build info object
0159   ii = minfo(mfilename, 'ao', '', utils.const.categories.sigproc, '$Id: spikecleaning.m,v 1.6 2008/09/05 11:05:29 ingo Exp $', sets, pl);
0160 end
0161 
0162 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0163 %
0164 % FUNCTION:    getDefaultPlist
0165 %
0166 % DESCRIPTION: Get Default Plist
0167 %
0168 % HISTORY:     11-07-07 M Hewitson
0169 %                Creation.
0170 %
0171 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0172 
0173 function plo = getDefaultPlist()
0174 
0175   plo = plist(       ...
0176     'kspike', 3.3,   ...
0177     'fc',     0.025, ...
0178     'order',  2,     ...
0179     'ripple', 0.5,   ...
0180     'method', 'random');
0181 end
0182

Generated on Mon 08-Sep-2008 13:18:47 by m2html © 2003