Home > m > sigproc > time_domain > ltpda_spikecleaning.m

ltpda_spikecleaning

PURPOSE ^

LTPDA_spikecleaning detects and corrects possible spikes in analysis objects

SYNOPSIS ^

function varargout=ltpda_spikecleaning(varargin)

DESCRIPTION ^

 LTPDA_spikecleaning detects and corrects possible spikes in analysis objects

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

 DESCRIPTION: LTPDA_ 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 = ltpda_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: ltpda_spikecleaning.m,v 1.2 2008/01/29 10:50:26 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=ltpda_spikecleaning(varargin)
0002 % LTPDA_spikecleaning detects and corrects possible spikes in analysis objects
0003 %
0004 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0005 %
0006 % DESCRIPTION: LTPDA_ 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 = ltpda_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: ltpda_spikecleaning.m,v 1.2 2008/01/29 10:50:26 josep Exp $
0033 %
0034 % HISTORY:     24-12-2007 J Sanjuan
0035 %                 Creation
0036 %
0037 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0038 
0039 ALGONAME = mfilename;
0040 VERSION  = '$Id: ltpda_spikecleaning.m,v 1.2 2008/01/29 10:50:26 josep Exp $';
0041 CATEGORY = 'SIGPROC';
0042 
0043 % Check if this is a call for parameters
0044 if nargin == 1
0045   in = char(varargin{1});
0046   if strcmp(in, 'Params')
0047     varargout{1} = getDefaultPL();
0048     return
0049   elseif strcmp(in, 'Version')
0050       varargout{1} = VERSION;
0051       return
0052   elseif strcmp(in, 'Category');
0053       varargout{1} = CATEGORY;
0054       return
0055   end
0056 end
0057 
0058 % Capture input variables names
0059 invars = {};
0060 as     = [];
0061 ps     = [];
0062 for j=1:nargin
0063    invars = [invars cellstr(inputname(j))];
0064    if isa(varargin{j}, 'ao')
0065      as = [as varargin{j}];
0066    end
0067    if isa(varargin{j}, 'plist')
0068      ps = [ps varargin{j}];
0069    end
0070 end
0071 
0072 % check plist
0073 if isempty(ps)
0074  pl = getDefaultPL();
0075 else
0076   pl = combine(ps, getDefaultPL);
0077 end
0078 
0079 % initialise output array
0080 bo = [];
0081 
0082 % go through each input AO
0083 na = length(as);
0084 for i=1:na
0085    a = as(i);
0086    d = a.data;
0087 
0088    % check this is a time-series object
0089    if ~isa(d, 'tsdata')
0090      error(' ### temperature spike detection requires tsdata (time-series) inputs.')
0091    end
0092 
0093    %--- check input parameters
0094    kspike = find(pl, 'kspike'); % kspike*sigma definition
0095    method = find(pl, 'method'); % method of spike-values substitution
0096    pl = pset(pl, param('gain', 1)); % gain of the filter
0097    pl = pset(pl, param('type', 'highpass')); % type of the filter
0098    fs = plist();
0099    fs = append(fs, param('fs', d.fs));
0100    pl = combine(pl, fs); % determination of the sampling frequency of the input AO
0101       
0102 % High-pass filtering data
0103 xfiltered = filtfilt(a, miir(pl));
0104 
0105 % standard deviation of the filtered data is calculated
0106 nxfiltered = find(abs(xfiltered) < kspike*std(xfiltered));
0107 
0108 xfiltered_2 = xfiltered.data.y(nxfiltered);
0109 
0110 std_xfiltered_2 = std(xfiltered_2);
0111 
0112 % spikes vector position is determined
0113 nspike = find(abs(xfiltered) > kspike*std_xfiltered_2);
0114 
0115 % substitution of spike values starts here
0116 xcleaned = a.data.y;
0117 for j=1:length(nspike)
0118     if nspike(j) <=2 % just in case a spike is detected in the 1st or 2nd sample
0119         xcleaned(nspike(j)) = mean(xcleaned(1:50)); 
0120     else
0121       if strcmp(method, 'random') % spike is substituted by a random value: N(0,std_xfiltered)
0122         xcleaned(nspike(j)) = xcleaned(nspike(j)-1) + randn(1)*std_xfiltered_2;
0123       elseif strcmp(method, 'mean') % spike is substituted by the mean if the two previous values
0124         xcleaned(nspike(j)) = (xcleaned(nspike(j)-1) + xcleaned(nspike(j)-2))/2;
0125       elseif strcmp(method, 'previous') % spike is substituted by the pervious value
0126         xcleaned(nspike(j)) = xcleaned(nspike(j)-1);       
0127       end
0128     end
0129 end
0130 
0131 % create new output tsdata
0132 ts = tsdata(xcleaned, d.fs);
0133 ts = set(ts, 'name', sprintf('cleaned data of %s', d.name));
0134 ts = set(ts, 'yunits', d.yunits);
0135 ts = set(ts, 'xunits', d.xunits);
0136 
0137 % create new output history
0138 h = history(ALGONAME, VERSION, pl, a.hist);
0139 h = set(h, 'invars', invars);
0140 
0141 % make output analysis object
0142 b = ao(ts, h);
0143 
0144 % set name
0145 if isempty(invars{i})
0146     n = a.name;
0147 else
0148     n = invars{i};
0149 end
0150 
0151 b = setnh(b, 'name', sprintf('spikecleaning(%s)', n));
0152 
0153 % add to output array
0154 bo = [bo b];
0155 
0156 varargout{1} = bo;
0157 end
0158 
0159 %-------------------------------
0160 % Get default parameters
0161 function plo = getDefaultPL();
0162 
0163 disp('* creating default plist...');
0164 plo = plist();
0165 plo = append(plo, param('kspike', 3.3));
0166 plo = append(plo, param('fc', 0.025));
0167 plo = append(plo, param('order', 2));
0168 plo = append(plo, param('ripple', 0.5));
0169 plo = append(plo, param('method', 'random'));
0170 disp('* done.');

Generated on Fri 07-Mar-2008 15:46:43 by m2html © 2003