0001 function varargout=ltpda_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 ALGONAME = mfilename;
0038 VERSION = '$Id: ltpda_gapfilling.m,v 1.3 2008/03/05 17:02:25 josep Exp $';
0039 CATEGORY = 'Signal Processing';
0040
0041
0042 if nargin == 1
0043 in = char(varargin{1});
0044 if strcmp(in, 'Params')
0045 varargout{1} = getDefaultPL();
0046 return
0047 elseif strcmp(in, 'Version')
0048 varargout{1} = VERSION;
0049 return
0050 elseif strcmp(in, 'Category')
0051 varargout{1} = CATEGORY;
0052 return
0053 end
0054 end
0055
0056
0057 in_names = {};
0058 for ii = 1:nargin
0059 in_names{end+1} = inputname(ii);
0060 end
0061 [as, pls, invars] = collect_inputs(varargin, in_names);
0062 na = length(as);
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078 if length(as)~=2
0079 error('only two analysis objects are needed!')
0080 end
0081
0082
0083 if isempty(pls)
0084 pls = getDefaultPL();
0085 else
0086 pls = combine(pls, getDefaultPL);
0087 end
0088
0089
0090 method = find(pls, 'method');
0091 addnoise = find(pls, 'addnoise');
0092
0093
0094
0095
0096 a1 = as(1).data;
0097 dinfo1 = whos('a1');
0098 dtype1 = dinfo1.class;
0099
0100
0101 a2 = as(2).data;
0102 dinfo2 = whos('a2');
0103 dtype2 = dinfo2.class;
0104
0105 if a1.fs ~= a2.fs
0106 warning('Sampling frequencies of the two AOs are different. The sampling frequency of the first AO will be used to reconstruct the gap.')
0107 end
0108
0109
0110 if ~strcmp(a1.xunits,a2.xunits)
0111 error('!!! Data has different X units !!!');
0112 end
0113
0114 if ~strcmp(a1.yunits,a2.yunits)
0115 error('!!! Data has different Y units !!!');
0116 end
0117
0118 a1_length = length(a1.x);
0119 a2_length = length(a2.x);
0120 gaptime = a2.x(1) - a1.x(a1_length);
0121 gapn = gaptime*a1.fs-1;
0122 t = (0:1:gapn-1)/a1.fs;
0123
0124
0125 if strcmp(method,'linear')
0126
0127 if length(a1.y)>10 && length(a2.y)>10
0128 dy = mean(a2.y(1:10))-mean(a1.y(a1_length-10:a1_length));
0129
0130 filling_data = (dy/gaptime)*t + mean(a1.y(a1_length-10:a1_length));
0131 filling_time = (1:1:gapn)/a1.fs + a1.x(a1_length);
0132 else
0133 error('!!! Not enough data in the data segments (min=11 for each one for the linear method) !!!');
0134 end
0135
0136 elseif strcmp(method,'spline')
0137
0138
0139 if length(a1.y)>1000 && length(a2.y)>1000
0140
0141
0142 da1 = diff(a1.y(1:100:a1_length))*(a1.fs/100);
0143 da1 = tsdata(da1, a1.fs/100);
0144 da1 = ao(da1);
0145
0146 da2 = diff(a2.y(1:100:a2_length))*(a2.fs/100);
0147 da2 = tsdata(da2, a2.fs/100);
0148 da2 = ao(da2);
0149
0150
0151
0152 plfa1 = getlpFilter(a1.fs/100);
0153 plfa2 = getlpFilter(a2.fs/100);
0154
0155 lpfa1 = miir(plfa1);
0156 lpfpla1 = plist(param('filter', lpfa1));
0157
0158 lpfa2 = miir(plfa2);
0159 lpfpla2 = plist(param('filter', lpfa2));
0160
0161
0162 da1filtered = filtfilt(da1, lpfpla1);
0163 da2filtered = filtfilt(da2, lpfpla2);
0164
0165
0166 c = mean(da1filtered.data.y(length(da1filtered.data.x)...
0167 -10:length(da1filtered.data.x)));
0168 d = mean(a1.y(length(a1.y)-10:length(a1.y)));
0169
0170 a=(2*d+(c+mean(da2filtered.data.y(1:10)))...
0171 *gaptime-2*mean(a2.y(1:10)))/(gaptime.^3);
0172 b=-(3*d+2*c*gaptime+mean(da2filtered.data.y(1:10))...
0173 *gaptime-3*mean(a2.y(1:10)))/(gaptime^2);
0174
0175
0176 filling_data = a*t.^3+b*t.^2+c*t+d;
0177 filling_time = (1:1:gapn)/a1.fs + a1.x(a1_length);
0178 else
0179 error('!!! Not enough data in data segments (min=1001 in spline method)');
0180 end
0181
0182 end
0183
0184
0185 if strcmp(addnoise,'yes');
0186
0187
0188 phpf = gethpFilter(a1.fs);
0189 ax = tsdata(a1.y, a1.fs);
0190 ax = ao(ax);
0191 hpf = miir(phpf);
0192 hpfpl = plist(param('filter', hpf));
0193 xhpf = filter(ax, hpfpl);
0194 hfnoise = std(xhpf);
0195
0196
0197 filling_data = filling_data + randn(length(filling_data),1)'*hfnoise.data.y;
0198
0199 end
0200
0201
0202 filling_data = [a1.y; filling_data'; a2.y];
0203 filling_time = [a1.x; filling_time'; a2.x];
0204
0205
0206 ts = tsdata(filling_time, filling_data)
0207 ts = set(ts, 'name', sprintf('filled data between %s and %s', as(1).name, as(2).name));
0208 ts = set(ts, 'yunits', a1.yunits);
0209 ts = set(ts, 'xunits', a1.xunits);
0210
0211
0212 histin = [as(1).hist as(2).hist];
0213 h = history(ALGONAME, VERSION, pls, histin);
0214 h = set(h, 'invars', invars);
0215
0216
0217 b = ao(ts, h);
0218
0219 b = setnh(b, 'name', sprintf('gapfilling(%s, %s)', as(1).name, as(2).name));
0220
0221 varargout{1} = b;
0222
0223
0224
0225 function plo = getDefaultPL()
0226
0227 plo = plist();
0228 plo = append(plo, param('method', 'linear'));
0229 plo = append(plo, param('addnoise', 'no'));
0230
0231
0232 function plf = getlpFilter(x)
0233
0234 plf = plist();
0235 plf = append(plf, param('gain', 1));
0236 plf = append(plf, param('ripple', 0.5));
0237 plf = append(plf, param('type', 'lowpass'));
0238 plf = append(plf, param('order', 2));
0239 plf = append(plf, param('fs', x));
0240 plf = append(plf, param('fc', 0.1/100));
0241
0242
0243 function phf = gethpFilter(x)
0244
0245 phf = plist();
0246 phf = append(phf, param('gain', 1));
0247 phf = append(phf, param('ripple', 0.5));
0248 phf = append(phf, param('type', 'highpass'));
0249 phf = append(phf, param('order', 2));
0250 phf = append(phf, param('fs', x));
0251 phf = append(phf, param('fc', 0.1/100));