Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/spikecleaning.m @ 0:f0afece42f48

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Import.
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 23 Nov 2011 19:22:13 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f0afece42f48
1 % spikecleaning detects and corrects possible spikes in analysis objects
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: SPIKECLEANING detects spikes in the temperature data and
5 % replaces them by artificial values depending on the method chosen ('random',
6 % 'mean', 'previous').
7 % Spikes are defined as singular samples with an (absolute) value
8 % higher than kspike times the standard deviation of the high-pass
9 % filtered (IIR filter) input AO.
10 %
11 % CALL: b = spikecleaning(a1, a2, ..., an, pl)
12 %
13 % INPUTS: aN - a list of analysis objects
14 % pl - parameter list
15 %
16 % OUTPUTS: b - a list of analysis objects with "spike values" removed
17 % and corrected
18 %
19 % <a href="matlab:utils.helper.displayMethodInfo('ao', 'spikecleaning')">Parameters Description</a>
20 %
21 % VERSION: $Id: spikecleaning.m,v 1.17 2011/04/08 08:56:16 hewitson Exp $
22 %
23 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
24
25 function varargout=spikecleaning(varargin)
26
27 %%% Check if this is a call for parameters
28 if utils.helper.isinfocall(varargin{:})
29 varargout{1} = getInfo(varargin{3});
30 return
31 end
32
33 if nargout == 0
34 error('### cat cannot be used as a modifier. Please give an output variable.');
35 end
36
37 % Collect input variable names
38 in_names = cell(size(varargin));
39 for ii = 1:nargin,in_names{ii} = inputname(ii);end
40
41 % Collect all AOs
42 [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
43 pli = utils.helper.collect_objects(varargin(:), 'plist', in_names);
44
45 pls = parse(pli, getDefaultPlist());
46
47 % initialise output array
48 bo = [];
49
50 % go through each input AO
51 for i=1:numel(as)
52 a = as(i);
53 d = a.data;
54
55 % check this is a time-series object
56 if ~isa(d, 'tsdata')
57 error(' ### temperature spike detection requires tsdata (time-series) inputs.')
58 end
59
60 %--- check input parameters
61 kspike = find(pls, 'kspike'); % kspike*sigma definition
62 method = find(pls, 'method'); % method of spike-values substitution
63 pls.pset('gain', 1); % gain of the filter
64 pls.pset('type', 'highpass'); % type of the filter
65 fs = plist();
66 fs.append('fs', d.fs);
67 pls.combine(fs); % determination of the sampling frequency of the input AO
68
69 % high-pass filtering data
70 xfiltered = filtfilt(a, miir(pls));
71
72 % standard deviation of the filtered data is calculated
73 nxfiltered = find(abs(xfiltered) < kspike*std(xfiltered));
74
75 xfiltered_2 = xfiltered.data.y(nxfiltered);
76
77 std_xfiltered_2 = std(xfiltered_2);
78
79 % spikes vector position is determined
80 nspike = find(abs(xfiltered) > kspike*std_xfiltered_2);
81
82 % substitution of spike values starts here
83 xcleaned = a.data.y;
84 for j=1:length(nspike)
85 if nspike(j) <=2 % just in case a spike is detected in the 1st or 2nd sample
86 xcleaned(nspike(j)) = mean(xcleaned(1:50));
87 else
88 if strcmp(method, 'random') % spike is substituted by a random value: N(0,std_xfiltered)
89 xcleaned(nspike(j)) = xcleaned(nspike(j)-1) + randn(1)*std_xfiltered_2;
90 elseif strcmp(method, 'mean') % spike is substituted by the mean if the two previous values
91 xcleaned(nspike(j)) = (xcleaned(nspike(j)-1) + xcleaned(nspike(j)-2))/2;
92 elseif strcmp(method, 'previous') % spike is substituted by the pervious value
93 xcleaned(nspike(j)) = xcleaned(nspike(j)-1);
94 end
95 end
96 end
97
98 % create new output tsdata
99 ts = tsdata(xcleaned, d.fs);
100 ts.setYunits(d.yunits);
101 ts.setXunits(d.xunits);
102
103 % % create new output history
104 % h = history(ALGONAME, VERSION, pls, a.hist);
105 % h = set(h, 'invars', invars);
106
107 % make output analysis object
108 b = ao(ts);
109 b.name = sprintf('spikecleaning(%s)', ao_invars{i});
110 b.addHistory(getInfo('None'), pls, ao_invars(i), as(i).hist);
111
112 % add to output array
113 bo = [bo b];
114
115 end
116
117 % Set output
118 if nargout == numel(bo)
119 % List of outputs
120 for ii = 1:numel(bo)
121 varargout{ii} = bo(ii);
122 end
123 else
124 % Single output
125 varargout{1} = bo;
126 end
127
128 end
129
130 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
131 % Local Functions %
132 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
133
134 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
135 %
136 % FUNCTION: getInfo
137 %
138 % DESCRIPTION: Get Info Object
139 %
140 % HISTORY: 11-07-07 M Hewitson
141 % Creation.
142 %
143 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
144
145 function ii = getInfo(varargin)
146 if nargin == 1 && strcmpi(varargin{1}, 'None')
147 sets = {};
148 pl = [];
149 else
150 sets = {'Default'};
151 pl = getDefaultPlist();
152 end
153 % Build info object
154 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);
155 ii.setModifier(false);
156 end
157
158 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
159 %
160 % FUNCTION: getDefaultPlist
161 %
162 % DESCRIPTION: Get Default Plist
163 %
164 % HISTORY: 11-07-07 M Hewitson
165 % Creation.
166 %
167 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
168
169 function plout = getDefaultPlist()
170 persistent pl;
171 if exist('pl', 'var')==0 || isempty(pl)
172 pl = buildplist();
173 end
174 plout = pl;
175 end
176
177 function pl = buildplist()
178
179 pl = plist();
180
181 % kspike
182 p = param({'kspike', 'High values imply no correction of relative low amplitude spikes.'}, paramValue.DOUBLE_VALUE(3.3));
183 pl.append(p);
184
185 % fc
186 p = param({'fc', 'Frequency cut-off of the IIR filter.'}, paramValue.DOUBLE_VALUE(0.025));
187 pl.append(p);
188
189 % Order
190 p = param({'order', 'The order of the IIR filter.'}, paramValue.DOUBLE_VALUE(2));
191 pl.append(p);
192
193 % Ripple
194 p = param({'ripple', 'Specify the pass/stop-band ripple for bandpass/bandreject filters'}, ...
195 paramValue.DOUBLE_VALUE(0.5));
196 pl.append(p);
197
198 % Method
199 p = param({'method', 'The method used to replace the spike value.'}, {1, {'random', 'mean'}, paramValue.SINGLE});
200 pl.append(p);
201
202 end
203
204 % PARAMETERES: 'kspike' - set the kspike value. High values imply
205 % not correction of relative low amplitude spike
206 % [default: 3.3]
207 % 'method' - method used to replace the spike value: 'random,
208 % 'mean', 'previous' [default:random]
209 % 'fc' - frequency cut-off of the IIR filter [default: 0.025]
210 % 'order' - order of the IIR filter [default: 2]
211 % 'ripple' - specify pass/stop-band ripple for bandpass
212 % and bandreject filters
213 % <<default: 0.5>>
214 %