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.2 2008/08/15 11:45:19 ingo Exp $

 HISTORY:     28-12-2007 J Sanjuan
          Creation

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

Generated on Mon 25-Aug-2008 22:39:29 by m2html © 2003