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.3 2008/03/05 17:02:25 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.3 2008/03/05 17:02:25 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.3 2008/03/05 17:02:25 josep Exp $';
0041 CATEGORY = 'Signal Processing';
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 % collect input ao's, plist's and ao variable names
0059 in_names = {};
0060 for ii = 1:nargin
0061     in_names{end+1} = inputname(ii);
0062 end
0063 [as, pls, invars] = collect_inputs(varargin, in_names);
0064 na = length(as);
0065 
0066 pls = combine(pls, getDefaultPL);
0067 
0068 % initialise output array
0069 bo = [];
0070 
0071 % go through each input AO
0072 for i=1:na
0073    a = as(i);
0074    d = a.data;
0075 
0076    % check this is a time-series object
0077    if ~isa(d, 'tsdata')
0078      error(' ### temperature spike detection requires tsdata (time-series) inputs.')
0079    end
0080 
0081    %--- check input parameters
0082    kspike = find(pls, 'kspike'); % kspike*sigma definition
0083    method = find(pls, 'method'); % method of spike-values substitution
0084    pls = pset(pls, param('gain', 1)); % gain of the filter
0085    pls = pset(pls, param('type', 'highpass')); % type of the filter
0086    fs = plist();
0087    fs = append(fs, param('fs', d.fs));
0088    pls = combine(pls, fs); % determination of the sampling frequency of the input AO
0089       
0090 % high-pass filtering data
0091 xfiltered = filtfilt(a, miir(pls));
0092 
0093 % standard deviation of the filtered data is calculated
0094 nxfiltered = find(abs(xfiltered) < kspike*std(xfiltered));
0095 
0096 xfiltered_2 = xfiltered.data.y(nxfiltered);
0097 
0098 std_xfiltered_2 = std(xfiltered_2);
0099 
0100 % spikes vector position is determined
0101 nspike = find(abs(xfiltered) > kspike*std_xfiltered_2);
0102 
0103 % substitution of spike values starts here
0104 xcleaned = a.data.y;
0105 for j=1:length(nspike)
0106     if nspike(j) <=2 % just in case a spike is detected in the 1st or 2nd sample
0107         xcleaned(nspike(j)) = mean(xcleaned(1:50)); 
0108     else
0109       if strcmp(method, 'random') % spike is substituted by a random value: N(0,std_xfiltered)
0110         xcleaned(nspike(j)) = xcleaned(nspike(j)-1) + randn(1)*std_xfiltered_2;
0111       elseif strcmp(method, 'mean') % spike is substituted by the mean if the two previous values
0112         xcleaned(nspike(j)) = (xcleaned(nspike(j)-1) + xcleaned(nspike(j)-2))/2;
0113       elseif strcmp(method, 'previous') % spike is substituted by the pervious value
0114         xcleaned(nspike(j)) = xcleaned(nspike(j)-1);       
0115       end
0116     end
0117 end
0118 
0119 % create new output tsdata
0120 ts = tsdata(xcleaned, d.fs);
0121 ts = set(ts, 'name', sprintf('cleaned data of %s', d.name));
0122 ts = set(ts, 'yunits', d.yunits);
0123 ts = set(ts, 'xunits', d.xunits);
0124 
0125 % create new output history
0126 h = history(ALGONAME, VERSION, pls, a.hist);
0127 h = set(h, 'invars', invars);
0128 
0129 % make output analysis object
0130 b = ao(ts, h);
0131 
0132 % set name
0133 if isempty(invars{i})
0134     n = a.name;
0135 else
0136     n = invars{i};
0137 end
0138 
0139 b = setnh(b, 'name', sprintf('spikecleaning(%s)', n));
0140 
0141 % add to output array
0142 bo = [bo b];
0143 
0144 varargout{1} = bo;
0145 end
0146 
0147 %-------------------------------
0148 % Get default parameters
0149 function plo = getDefaultPL();
0150 
0151 disp('* creating default plist...');
0152 plo = plist();
0153 plo = append(plo, param('kspike', 3.3));
0154 plo = append(plo, param('fc', 0.025));
0155 plo = append(plo, param('order', 2));
0156 plo = append(plo, param('ripple', 0.5));
0157 plo = append(plo, param('method', 'random'));
0158 disp('* done.');

Generated on Mon 31-Mar-2008 13:54:54 by m2html © 2003