Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/gapfilling.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 % GAPFILLING fills possible gaps in data.
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: GAPFILLING interpolated data between two data
5 % segments. This function might be useful for possible
6 % gaps or corrupted data. Two different types of
7 % interpolating are available: linear and spline, the latter
8 % results in a smoother curve between the two data segments.
9 %
10 % CALL: b = gapfilling(a1, a2, pl)
11 %
12 % INPUTS: a1 - data segment previous to the gap
13 % a2 - data segment posterior to the gap
14 % pl - parameter list
15 %
16 % OUTPUTS: b - data segment containing a1, a2 and the filled data
17 % segment, i.e., b=[a1 datare_filled a2].
18 %
19 % <a href="matlab:utils.helper.displayMethodInfo('ao', 'gapfilling')">Parameters Description</a>
20 %
21 % VERSION: $Id: gapfilling.m,v 1.21 2011/11/07 15:34:34 miquel Exp $
22 %
23 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
24
25 function varargout = gapfilling(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 if length(as)~=2
48 error('only two analysis objects are needed!')
49 end
50
51 % go through each input AO
52 for i=1:numel(as)
53 % check this is a time-series object
54 if ~isa(as(i).data, 'tsdata')
55 error(' ### Gap filling requires tsdata (time-series) inputs.')
56 end
57 end
58
59 %--- check input parameters
60 method = find(pls, 'method'); % method definition: linear or 'spline'
61 addnoise = find(pls, 'addnoise'); % decide whether add noise or not in the filled data
62
63 %---------------------------
64 % rename AOs
65 a1 = as(1);
66 a2 = as(2);
67
68 if a1.fs ~= a2.fs
69 warning('Sampling frequencies of the two AOs are different. The sampling frequency of the first AO will be used to reconstruct the gap.')
70 end
71
72 % Different units error
73 if a1.xunits ~= a2.xunits
74 error('!!! Data has different X units !!!');
75 end
76
77 if a1.yunits ~= a2.yunits
78 error('!!! Data has different Y units !!!');
79 end
80
81 a1_length = len(a1);
82 a2_length = len(a2);
83 start_a2 = a2.t0.utc_epoch_milli/1000 + a2.x(1);
84 end_a1 = (a1.t0.utc_epoch_milli/1000 + a1.x(1) + a1.nsecs - 1);
85 gaptime = (start_a2 - end_a1);
86 gapn = gaptime*a1.fs-1;
87 t = (0:1:gapn-1)'/a1.fs;
88
89 %--- gapfilling process itself
90 if strcmp(method,'linear')
91 % linear interpolation method ---xfilled=(deltay/deltax)*t+y1(length(y1))---
92 if len(a1)>10 && len(a2)>10
93 dy = mean(a2.y(1:10))-mean(a1.y(a1_length-10:a1_length));
94
95 filling_data = (dy/gaptime)*t + mean(a1.y(a1_length-10:a1_length));
96 filling_time = (1:1:gapn)'/a1.fs + a1.x(a1_length);
97 else
98 error('!!! Not enough data in the data segments (min=11 for each one for the linear method) !!!');
99 end
100
101 elseif strcmp(method,'spline') % spline method xfilled = a*T^3 + b*T^2 + c*T +d
102
103 if len(a1)>1000 && len(a2)>1000
104
105 % derivatives of the input data are calculated
106 da1 = diff(a1.y(1:100:a1_length))*(a1.fs/100);
107 da1 = tsdata(da1, a1.fs/100);
108 da1 = ao(da1);
109
110 da2 = diff(a2.y(1:100:a2_length))*(a2.fs/100);
111 da2 = tsdata(da2, a2.fs/100);
112 da2 = ao(da2);
113
114 % This filters the previous derivatives
115 % filters parameters are obtained
116 plfa1 = getlpFilter(a1.fs/100);
117 plfa2 = getlpFilter(a2.fs/100);
118
119 lpfa1 = miir(plfa1);
120 lpfpla1 = plist(param('filter', lpfa1));
121
122 lpfa2 = miir(plfa2);
123 lpfpla2 = plist(param('filter', lpfa2));
124
125 % derivatives are low-pass filtered
126 da1filtered = filtfilt(da1, lpfpla1);
127 da2filtered = filtfilt(da2, lpfpla2);
128
129 % coefficients are calculated
130 c = mean(da1filtered.y(len(da1filtered)...
131 -10:len(da1filtered)));
132 d = mean(a1.y(len(a1)-10:len(a1)));
133
134 a=(2*d+(c+mean(da2filtered.y(1:10)))...
135 *gaptime-2*mean(a2.y(1:10)))/(gaptime.^3);
136
137 b=-(3*d+2*c*gaptime+mean(da2filtered.y(1:10))...
138 *gaptime-3*mean(a2.y(1:10)))/(gaptime^2);
139
140 % filling data is calculated with the coefficients a, b, c and d
141 filling_data = a*t.^3+b*t.^2+c*t+d;
142 filling_time = (1:1:gapn)'/a1.fs + a1.x(a1_length);
143 else
144 error('!!! Not enough data in data segments (min=1001 in spline method)');
145 end
146
147 end
148
149 % this add noise (if desired) to the filled gap
150 if strcmp(addnoise,'yes');
151 % calculation of the standard deviation after eliminating the low-frequency component
152 phpf = gethpFilter(a1.fs);
153 ax = tsdata(a1.y, a1.fs);
154 ax = ao(ax);
155 hpf = miir(phpf);
156 hpfpl = plist(param('filter', hpf));
157 xhpf = filter(ax, hpfpl);
158 hfnoise = std(xhpf);
159
160 % noise is added to the filling data
161 filling_data = filling_data + randn(length(filling_data),1)*hfnoise.data.getY;
162 end
163
164 % join data
165 filling_data = [a1.y; filling_data; a2.y];
166 filling_time = [a1.x; filling_time; a2.x];
167
168 % create new output tsdata
169 ts = tsdata(filling_time, filling_data);
170 ts.setYunits(a1.yunits);
171 ts.setXunits(a1.xunits);
172
173 % make output analysis object
174 b = ao(ts);
175 b.setT0(a1.t0)
176 b.name = sprintf('gapfilling(%s,%s)', ao_invars{1}, ao_invars{2});
177 b.addHistory(getInfo('None'), pls, [ao_invars(1) ao_invars(2)], [as(1).hist as(2).hist]);
178
179 % Set output
180 if nargout == numel(b)
181 % List of outputs
182 for ii = 1:numel(b)
183 varargout{ii} = b(ii);
184 end
185 else
186 % Single output
187 varargout{1} = b;
188 end
189
190 end
191
192 function plf = getlpFilter(x)
193
194 plf = plist();
195 plf = append(plf, param('gain', 1));
196 plf = append(plf, param('ripple', 0.5));
197 plf = append(plf, param('type', 'lowpass'));
198 plf = append(plf, param('order', 2));
199 plf = append(plf, param('fs', x));
200 plf = append(plf, param('fc', 0.1/100));
201
202 end
203
204 function phf = gethpFilter(x)
205
206 phf = plist();
207 phf = append(phf, param('gain', 1));
208 phf = append(phf, param('ripple', 0.5));
209 phf = append(phf, param('type', 'highpass'));
210 phf = append(phf, param('order', 2));
211 phf = append(phf, param('fs', x));
212 phf = append(phf, param('fc', 0.1/100));
213
214 end
215
216 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
217 % Local Functions %
218 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
219
220 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
221 %
222 % FUNCTION: getInfo
223 %
224 % DESCRIPTION: Get Info Object
225 %
226 % HISTORY: 11-07-07 M Hewitson
227 % Creation.
228 %
229 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
230
231 function ii = getInfo(varargin)
232
233 if nargin == 1 && strcmpi(varargin{1}, 'None')
234 sets = {};
235 pl = [];
236 else
237 sets = {'Default'};
238 pl = getDefaultPlist();
239 end
240 % Build info object
241 ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: gapfilling.m,v 1.21 2011/11/07 15:34:34 miquel Exp $', sets, pl);
242 ii.setModifier(false);
243 end
244
245 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
246 %
247 % FUNCTION: getDefaultPlist
248 %
249 % DESCRIPTION: Get Default Plist
250 %
251 % HISTORY: 11-07-07 M Hewitson
252 % Creation.
253 %
254 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
255
256 function plout = getDefaultPlist()
257 persistent pl;
258 if exist('pl', 'var')==0 || isempty(pl)
259 pl = buildplist();
260 end
261 plout = pl;
262 end
263
264 function pl = buildplist()
265
266 pl = plist();
267
268 % Method
269 p = param({'method', 'The method used to interpolate data.'}, {1, {'linear', 'spline'}, paramValue.SINGLE});
270 pl.append(p);
271
272 % Add noise
273 p = param({'addnoise', ...
274 ['Noise can be added to the interpolated data.<br>'...
275 'This noise is defined as random variable with<br>'...
276 'zero mean and variance equal to the high-frequency<br>'...
277 'noise of the first input.']}, paramValue.YES_NO);
278 p.val.setValIndex(2);
279 pl.append(p);
280
281 end
282
283 % PARAMETERES: 'method' - method used to interpolate data between a1 and a2.
284 % Two options can be used: 'linear' and 'spline'.
285 % Default values is 'linear'.
286 % 'addnoise' - noise can be added to the interpolated data.
287 % This noise is defined as random variable with
288 % zero mean and variance equal to the high-frequency
289 % noise if a1. 'Yes' adds noise. Default value
290 % is 'no'.