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