0001 function varargout = gapfilling(varargin)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039 if utils.helper.isinfocall(varargin{:})
0040 varargout{1} = getInfo(varargin{3});
0041 return
0042 end
0043
0044 if nargout == 0
0045 error('### cat cannot be used as a modifier. Please give an output variable.');
0046 end
0047
0048
0049 in_names = cell(size(varargin));
0050 for ii = 1:nargin,in_names{ii} = inputname(ii);end
0051
0052
0053 [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
0054 pli = utils.helper.collect_objects(varargin(:), 'plist', in_names);
0055
0056 pls = combine(pli, getDefaultPlist());
0057
0058 if length(as)~=2
0059 error('only two analysis objects are needed!')
0060 end
0061
0062
0063 for i=1:numel(as)
0064 a = as(i);
0065 d = a.data;
0066 end
0067
0068
0069 if ~isa(d, 'tsdata')
0070 error(' ### temperature spike detection requires tsdata (time-series) inputs.')
0071 end
0072
0073
0074 method = find(pls, 'method');
0075 addnoise = find(pls, 'addnoise');
0076
0077
0078
0079 a1 = as(1).data;
0080 dinfo1 = whos('a1');
0081 dtype1 = dinfo1.class;
0082
0083
0084 a2 = as(2).data;
0085 dinfo2 = whos('a2');
0086 dtype2 = dinfo2.class;
0087
0088 if a1.fs ~= a2.fs
0089 warning('Sampling frequencies of the two AOs are different. The sampling frequency of the first AO will be used to reconstruct the gap.')
0090 end
0091
0092
0093 if a1.xunits ~= a2.xunits
0094 error('!!! Data has different X units !!!');
0095 end
0096
0097 if a1.yunits ~= a2.yunits
0098 error('!!! Data has different Y units !!!');
0099 end
0100
0101 a1_length = length(a1.getX);
0102 a2_length = length(a2.getX);
0103 gaptime = a2.getX(1) - a1.getX(a1_length);
0104 gapn = gaptime*a1.fs-1;
0105 t = (0:1:gapn-1)/a1.fs;
0106
0107
0108 if strcmp(method,'linear')
0109
0110 if length(a1.getY)>10 && length(a2.getY)>10
0111 dy = mean(a2.getY(1:10))-mean(a1.getY(a1_length-10:a1_length));
0112
0113 filling_data = (dy/gaptime)*t + mean(a1.getY(a1_length-10:a1_length));
0114 filling_time = (1:1:gapn)/a1.fs + a1.getX(a1_length);
0115 else
0116 error('!!! Not enough data in the data segments (min=11 for each one for the linear method) !!!');
0117 end
0118
0119 elseif strcmp(method,'spline')
0120
0121 if length(a1.y)>1000 && length(a2.y)>1000
0122
0123
0124 da1 = diff(a1.y(1:100:a1_length))*(a1.fs/100);
0125 da1 = tsdata(da1, a1.fs/100);
0126 da1 = ao(da1);
0127
0128 da2 = diff(a2.y(1:100:a2_length))*(a2.fs/100);
0129 da2 = tsdata(da2, a2.fs/100);
0130 da2 = ao(da2);
0131
0132
0133
0134 plfa1 = getlpFilter(a1.fs/100);
0135 plfa2 = getlpFilter(a2.fs/100);
0136
0137 lpfa1 = miir(plfa1);
0138 lpfpla1 = plist(param('filter', lpfa1));
0139
0140 lpfa2 = miir(plfa2);
0141 lpfpla2 = plist(param('filter', lpfa2));
0142
0143
0144 da1filtered = filtfilt(da1, lpfpla1);
0145 da2filtered = filtfilt(da2, lpfpla2);
0146
0147
0148 c = mean(da1filtered.data.getY(length(da1filtered.data.getX)...
0149 -10:length(da1filtered.data.getX)));
0150 d = mean(a1.getY(length(a1.getY)-10:length(a1.getY)));
0151
0152 a=(2*d+(c+mean(da2filtered.data.getY(1:10)))...
0153 *gaptime-2*mean(a2.getY(1:10)))/(gaptime.^3);
0154
0155 b=-(3*d+2*c*gaptime+mean(da2filtered.data.getY(1:10))...
0156 *gaptime-3*mean(a2.getY(1:10)))/(gaptime^2);
0157
0158
0159 filling_data = a*t.^3+b*t.^2+c*t+d;
0160 filling_time = (1:1:gapn)/a1.fs + a1.getX(a1_length);
0161 else
0162 error('!!! Not enough data in data segments (min=1001 in spline method)');
0163 end
0164
0165 end
0166
0167
0168 if strcmp(addnoise,'yes');
0169
0170 phpf = gethpFilter(a1.fs);
0171 ax = tsdata(a1.y, a1.fs);
0172 ax = ao(ax);
0173 hpf = miir(phpf);
0174 hpfpl = plist(param('filter', hpf));
0175 xhpf = filter(ax, hpfpl);
0176 hfnoise = std(xhpf);
0177
0178
0179 filling_data = filling_data + randn(length(filling_data),1)'*hfnoise.data.y;
0180 end
0181
0182
0183 filling_data = [a1.getY filling_data a2.getY];
0184 filling_time = [a1.getX filling_time a2.getX];
0185
0186
0187 ts = tsdata(filling_time, filling_data);
0188 ts.setYunits(a1.yunits);
0189 ts.setXunits(a1.xunits);
0190
0191
0192
0193
0194
0195
0196
0197 b = ao(ts);
0198 b.setName(sprintf('gapfilling(%s&%s)', ao_invars{1}, ao_invars{2}), 'internal');
0199
0200
0201
0202
0203 b.addHistory(getInfo(), pli, [ao_invars(1) ao_invars(2)], [as(1).hist as(2).hist]);
0204
0205 varargout{1} = b;
0206
0207 end
0208
0209 function plf = getlpFilter(x)
0210
0211 plf = plist();
0212 plf = append(plf, param('gain', 1));
0213 plf = append(plf, param('ripple', 0.5));
0214 plf = append(plf, param('type', 'lowpass'));
0215 plf = append(plf, param('order', 2));
0216 plf = append(plf, param('fs', x));
0217 plf = append(plf, param('fc', 0.1/100));
0218
0219 end
0220
0221 function phf = gethpFilter(x)
0222
0223 phf = plist();
0224 phf = append(phf, param('gain', 1));
0225 phf = append(phf, param('ripple', 0.5));
0226 phf = append(phf, param('type', 'highpass'));
0227 phf = append(phf, param('order', 2));
0228 phf = append(phf, param('fs', x));
0229 phf = append(phf, param('fc', 0.1/100));
0230
0231 end
0232
0233
0234
0235
0236
0237
0238
0239
0240
0241
0242
0243
0244
0245
0246
0247
0248 function ii = getInfo(varargin)
0249
0250 if nargin == 1 && strcmpi(varargin{1}, 'None')
0251 sets = {};
0252 pl = [];
0253 else
0254 sets = {'Default'};
0255 pl = getDefaultPlist();
0256 end
0257
0258 ii = minfo(mfilename, 'ao', '', 'Signal Processing', '$Id: gapfilling.m,v 1.2 2008/08/15 11:45:19 ingo Exp $', sets, pl);
0259
0260 end
0261
0262
0263
0264
0265
0266
0267
0268
0269
0270
0271
0272
0273 function plo = getDefaultPlist()
0274
0275 plo = plist();
0276 plo = append(plo, param('method', 'linear'));
0277 plo = append(plo, param('addnoise', 'no'));
0278
0279 end
0280