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