Home > m > sigproc > time_domain > ltpda_gapfilling.m

ltpda_gapfilling

PURPOSE ^

LTPDA_GAPFILLING fills possible gaps in data.

SYNOPSIS ^

function varargout=ltpda_gapfilling(varargin)

DESCRIPTION ^

 LTPDA_GAPFILLING fills possible gaps in data.

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

 DESCRIPTION: LTPDA_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 = ltpda_spikecleaning(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: ltpda_gapfilling.html,v 1.7 2008/03/31 10:27:44 hewitson 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=ltpda_gapfilling(varargin)
0002 % LTPDA_GAPFILLING fills possible gaps in data.
0003 %
0004 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0005 %
0006 % DESCRIPTION: LTPDA_GAPFILLING interpolated data between two data
0007 %              segments. This function might be useful for possible
0008 %              gaps or corrupted data. Two different types of
0009 %              interpolating are available: linear and spline, the latter
0010 %              results in a smoother curve between the two data segments.
0011 %
0012 % CALL:        b = ltpda_spikecleaning(a1, a2, pl)
0013 %
0014 % INPUTS:      a1 - data segment previous to the gap
0015 %               a2 - data segment posterior to the gap
0016 %               pl - parameter list
0017 %
0018 % OUTPUTS:     b - data segment containing a1, a2 and the filled data
0019 %                   segment, i.e., b=[a1 datare_filled a2].
0020 %
0021 % PARAMETERES: 'method' - method used to interpolate data between a1 and a2.
0022 %                         Two options can be used: 'linear' and 'spline'.
0023 %                         Default values is 'linear'.
0024 %              'addnoise' - noise can be added to the interpolated data.
0025 %                           This noise is defined as random variable with
0026 %                           zero mean and variance equal to the high-frequency
0027 %                           noise if a1. 'Yes' adds noise. Default value
0028 %                           is 'no'.
0029 %
0030 % VERSION:     $Id: ltpda_gapfilling.html,v 1.7 2008/03/31 10:27:44 hewitson Exp $
0031 %
0032 % HISTORY:     28-12-2007 J Sanjuan
0033 %          Creation
0034 %
0035 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0036 
0037 ALGONAME = mfilename;
0038 VERSION  = '$Id: ltpda_gapfilling.html,v 1.7 2008/03/31 10:27:44 hewitson Exp $';
0039 CATEGORY = 'Signal Processing';
0040 
0041 % check if this is a call for parameters
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 % collect input ao's, plist's and ao variable names
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 % Capture input variables names
0065 %invars = {};
0066 %as     = [];
0067 %ps     = [];
0068 %for j=1:nargin
0069 %   invars = [invars cellstr(inputname(j))];
0070 %   if isa(varargin{j}, 'ao')
0071 %     as = [as varargin{j}];
0072 %   end
0073 %   if isa(varargin{j}, 'plist')
0074 %     ps = [ps varargin{j}];
0075 %   end
0076 %end
0077 
0078 if length(as)~=2
0079    error('only two analysis objects are needed!')
0080 end
0081 
0082 % check plist
0083 if isempty(pls)
0084   pls = getDefaultPL();
0085 else
0086   pls = combine(pls, getDefaultPL);
0087 end
0088 
0089 % check input parameters
0090 method = find(pls, 'method'); % method definition: linear or 'spline'
0091 addnoise = find(pls, 'addnoise'); % decide whether add noise or not in
0092                                  % the filled data
0093 
0094 %---------------------------
0095 % Get data type from the first AO
0096 a1    = as(1).data;
0097 dinfo1 = whos('a1');
0098 dtype1 = dinfo1.class;
0099 
0100 % Get data type from the second AO
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 % Different units error
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 % gapfilling process itself
0125 if strcmp(method,'linear') 
0126 % linear interpolation method ---xfilled=(deltay/deltax)*t+y1(length(y1))---
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') % spline method
0137 % xfilled = a*T^3 + b*T^2 + c*T +d
0138 
0139    if length(a1.y)>1000 && length(a2.y)>1000
0140 
0141     % derivatives of the input data are calculated
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     % This filters the previous derivatives
0151     % filters parameters are obtained
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     % derivatives are low-pass filtered
0162     da1filtered = filtfilt(da1, lpfpla1);
0163     da2filtered = filtfilt(da2, lpfpla2);
0164 
0165     % coefficients are calculated
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     % filling data is calculated with the coefficients a, b, c and d
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 % this add noise (if desired) to the filled gap
0185 if strcmp(addnoise,'yes');
0186   % calculation of the standard deviation after eliminating
0187   % the low-frequency component
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   % noise is added to the filling data
0197   filling_data = filling_data + randn(length(filling_data),1)'*hfnoise.data.y;
0198 
0199 end
0200 
0201 % join data
0202 filling_data = [a1.y; filling_data'; a2.y];
0203 filling_time = [a1.x; filling_time'; a2.x];
0204 
0205 % create new output tsdata
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 % create new output history
0212 histin = [as(1).hist as(2).hist];
0213 h = history(ALGONAME, VERSION, pls, histin);
0214 h = set(h, 'invars', invars);
0215 
0216 % make output analysis object
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 % Get default params
0225 function plo = getDefaultPL()
0226 
0227 plo = plist();
0228 plo = append(plo, param('method', 'linear'));
0229 plo = append(plo, param('addnoise', 'no'));
0230 
0231 % Get low-pass filter
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 % Get high-pass filter
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));

Generated on Mon 31-Mar-2008 12:20:24 by m2html © 2003