Home > classes > @ao > gapfilling.m

gapfilling

PURPOSE ^

GAPFILLING fills possible gaps in data.

SYNOPSIS ^

function varargout = gapfilling(varargin)

DESCRIPTION ^

 GAPFILLING fills possible gaps in data.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 DESCRIPTION: GAPFILLING interpolated data between two data
              segments. This function might be useful for possible
              gaps or corrupted data. Two different types of
              interpolating are available: linear and spline, the latter
              results in a smoother curve between the two data segments.

 CALL:        b = gapfilling(a1, a2, pl)

 INPUTS:      a1 - data segment previous to the gap
               a2 - data segment posterior to the gap
               pl - parameter list

 OUTPUTS:     b - data segment containing a1, a2 and the filled data
                  segment, i.e., b=[a1 datare_filled a2].

 PARAMETERES: 'method' - method used to interpolate data between a1 and a2.
                         Two options can be used: 'linear' and 'spline'.
                         Default values is 'linear'.
              'addnoise' - noise can be added to the interpolated data.
                           This noise is defined as random variable with
                           zero mean and variance equal to the high-frequency
                           noise if a1. 'Yes' adds noise. Default value
                           is 'no'.

 VERSION:     $Id: gapfilling.m,v 1.5 2008/09/05 11:05:29 ingo Exp $

 HISTORY:     28-12-2007 J Sanjuan
          Creation

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 % GAPFILLING fills possible gaps in data.
0002 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0003 %
0004 % DESCRIPTION: GAPFILLING interpolated data between two data
0005 %              segments. This function might be useful for possible
0006 %              gaps or corrupted data. Two different types of
0007 %              interpolating are available: linear and spline, the latter
0008 %              results in a smoother curve between the two data segments.
0009 %
0010 % CALL:        b = gapfilling(a1, a2, pl)
0011 %
0012 % INPUTS:      a1 - data segment previous to the gap
0013 %               a2 - data segment posterior to the gap
0014 %               pl - parameter list
0015 %
0016 % OUTPUTS:     b - data segment containing a1, a2 and the filled data
0017 %                  segment, i.e., b=[a1 datare_filled a2].
0018 %
0019 % PARAMETERES: 'method' - method used to interpolate data between a1 and a2.
0020 %                         Two options can be used: 'linear' and 'spline'.
0021 %                         Default values is 'linear'.
0022 %              'addnoise' - noise can be added to the interpolated data.
0023 %                           This noise is defined as random variable with
0024 %                           zero mean and variance equal to the high-frequency
0025 %                           noise if a1. 'Yes' adds noise. Default value
0026 %                           is 'no'.
0027 %
0028 % VERSION:     $Id: gapfilling.m,v 1.5 2008/09/05 11:05:29 ingo Exp $
0029 %
0030 % HISTORY:     28-12-2007 J Sanjuan
0031 %          Creation
0032 %
0033 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0034 
0035 function varargout = gapfilling(varargin)
0036 
0037   %%% Check if this is a call for parameters
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   % Collect input variable names
0048   in_names = cell(size(varargin));
0049   for ii = 1:nargin,in_names{ii} = inputname(ii);end
0050 
0051   % Collect all AOs
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   % go through each input AO
0062   for i=1:numel(as)
0063     a = as(i);
0064     d = a.data;
0065   end
0066 
0067   % check this is a time-series object
0068   if ~isa(d, 'tsdata')
0069     error(' ### temperature spike detection requires tsdata (time-series) inputs.')
0070   end
0071 
0072   %--- check input parameters
0073   method = find(pls, 'method'); % method definition: linear or 'spline'
0074   addnoise = find(pls, 'addnoise'); % decide whether add noise or not in the filled data
0075 
0076   %---------------------------
0077   % Get data type from the first AO
0078   a1    = as(1).data;
0079   dinfo1 = whos('a1');
0080   dtype1 = dinfo1.class;
0081 
0082   % Get data type from the second AO
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   % Different units error
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   %--- gapfilling process itself
0107   if strcmp(method,'linear')
0108     % linear interpolation method ---xfilled=(deltay/deltax)*t+y1(length(y1))---
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') % spline method xfilled = a*T^3 + b*T^2 + c*T +d
0119 
0120     if length(a1.getY)>1000 && length(a2.getY)>1000
0121 
0122       % derivatives of the input data are calculated
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       % This filters the previous derivatives
0132       % filters parameters are obtained
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       % derivatives are low-pass filtered
0143       da1filtered = filtfilt(da1, lpfpla1);
0144       da2filtered = filtfilt(da2, lpfpla2);
0145 
0146       % coefficients are calculated
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       % filling data is calculated with the coefficients a, b, c and d
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   % this add noise (if desired) to the filled gap
0167   if strcmp(addnoise,'yes');
0168     % calculation of the standard deviation after eliminating the low-frequency component
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     % noise is added to the filling data
0178     filling_data = filling_data + randn(length(filling_data),1)*hfnoise.data.getY;
0179   end
0180 
0181   % join data
0182   filling_data = [a1.getY; filling_data; a2.getY];
0183   filling_time = [a1.getX; filling_time; a2.getX];
0184 
0185   % create new output tsdata
0186   ts = tsdata(filling_time, filling_data);
0187   ts.setYunits(a1.yunits);
0188   ts.setXunits(a1.xunits);
0189 
0190   % create new output history
0191   %%histin = [as(1).hist as(2).hist];
0192   %%h = history(ALGONAME, VERSION, pls, histin);
0193   %%h = set(h, 'invars', invars);
0194 
0195   % make output analysis object
0196   b = ao(ts);
0197   b.setName(sprintf('gapfilling(%s&%s)', ao_invars{1}, ao_invars{2}), 'internal');
0198   %warning('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%');
0199   %warning('%%%%%%%%%%%%  Uncomment the following line (add Hsitory)  %%%%%');
0200   %warning('%%%%%%%%%%%%  Uncomment the following line (add Hsitory)  %%%%%');
0201   %warning('%%%%%%%%%%%%  Uncomment the following line (add Hsitory)  %%%%%');
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 %                               Local Functions                               %
0234 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0235 
0236 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0237 %
0238 % FUNCTION:    getInfo
0239 %
0240 % DESCRIPTION: Get Info Object
0241 %
0242 % HISTORY:     11-07-07 M Hewitson
0243 %                Creation.
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   % Build info object
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 % FUNCTION:    getDefaultPlist
0264 %
0265 % DESCRIPTION: Get Default Plist
0266 %
0267 % HISTORY:     11-07-07 M Hewitson
0268 %                Creation.
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

Generated on Mon 08-Sep-2008 13:18:47 by m2html © 2003