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.3 2008/08/08 15:31:38 josep Exp $

 HISTORY:     24-12-2007 J Sanjuan
                 Creation

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function varargout=spikecleaning(varargin)
0002   % spikecleaning detects and corrects possible spikes in analysis objects
0003   %
0004   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0005   %
0006   % DESCRIPTION: 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 = 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: spikecleaning.m,v 1.3 2008/08/08 15:31:38 josep Exp $
0033   %
0034   % HISTORY:     24-12-2007 J Sanjuan
0035   %                 Creation
0036   %
0037   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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(), pli, 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', '', 'Signal Processing', '$Id: spikecleaning.m,v 1.3 2008/08/08 15:31:38 josep 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 25-Aug-2008 22:39:29 by m2html © 2003