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.m,v 1.2 2008/01/29 10:55:27 josep 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.m,v 1.2 2008/01/29 10:55:27 josep Exp $
0031 %
0032 % HISTORY:     28-12-2007 J Sanjuan
0033 %          Creation
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 % 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 % Capture input variables names
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 % check plist
0075 if isempty(ps)
0076   pl = getDefaultPL();
0077 else
0078   pl = combine(ps, getDefaultPL);
0079 end
0080 
0081 % check input parameters
0082 method = find(pl, 'method'); % method definition: linear or 'spline'
0083 addnoise = find(pl, 'addnoise'); % decide whether add noise or not in
0084                                  % the filled data
0085 
0086 %---------------------------
0087 % Get data type from the first AO
0088 a1    = as(1).data;
0089 dinfo1 = whos('a1');
0090 dtype1 = dinfo1.class;
0091 
0092 % Get data type from the second AO
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 % Different units error
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 % gapfilling process itself
0117 if strcmp(method,'linear') 
0118 % linear interpolation method ---xfilled=(deltay/deltax)*t+y1(length(y1))---
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') % spline method
0129 % xfilled = a*T^3 + b*T^2 + c*T +d
0130 
0131    if length(a1.y)>1000 && length(a2.y)>1000
0132 
0133     % derivatives of the input data are calculated
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     % This filters the previous derivatives
0143     % filters parameters are obtained
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     % derivatives are low-pass filtered
0154     da1filtered = filtfilt(da1, lpfpla1);
0155     da2filtered = filtfilt(da2, lpfpla2);
0156 
0157     % coefficients are calculated
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     % filling data is calculated with the coefficients a, b, c and d
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 % this add noise (if desired) to the filled gap
0177 if strcmp(addnoise,'yes');
0178   % calculation of the standard deviation after eliminating
0179   % the low-frequency component
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   % noise is added to the filling data
0189   filling_data = filling_data + randn(length(filling_data),1)'*hfnoise.data.y;
0190 
0191 end
0192 
0193 % join data
0194 filling_data = [a1.y; filling_data'; a2.y];
0195 filling_time = [a1.x; filling_time'; a2.x];
0196 
0197 % create new output tsdata
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 % create new output history
0204 histin = [as(1).hist as(2).hist];
0205 h = history(ALGONAME, VERSION, pl, histin);
0206 h = set(h, 'invars', invars);
0207 
0208 % make output analysis object
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 % Get default params
0217 function plo = getDefaultPL()
0218 
0219 plo = plist();
0220 plo = append(plo, param('method', 'linear'));
0221 plo = append(plo, param('addnoise', 'no'));
0222 
0223 % Get low-pass filter
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 % Get high-pass filter
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));

Generated on Fri 07-Mar-2008 15:46:43 by m2html © 2003