Home > wafo > onedim > reconstruct.m

reconstruct

PURPOSE ^

reconstruct the spurious/missing points of timeseries

SYNOPSIS ^

[y,g,g2,test,tobs,mu1o, mu1oStd]=reconstruct(x,inds,Nsim,L,def,varargin)

DESCRIPTION ^

 RECONSTRUCT reconstruct the spurious/missing points of timeseries
 
  CALL: [y,g,g2,test,tobs,mu1o,mu1oStd]=reconstruct(x,inds,Nsim,L,def,options);
    
        y   = reconstructed signal
       g,g2 = smoothed and empirical transformation, respectively
  test,tobs = test observator int(g(u)-u)^2 du and int(g_new(u)-g_old(u))^2 du,
              respectively, where int limits is given by param in lc2tr. 
              Test is a measure of departure from the Gaussian model for 
              the data. Tobs is a measure of the convergence of the 
              estimation of g.
       mu1o = expected surface elevation of the Gaussian model process.
    mu1oStd = standarddeviation of mu1o.
 
        x   = 2 column timeseries 
              first column sampling times [sec]
              second column surface elevation [m]
       inds = indices to spurious points of x
       Nsim = the maximum # of iterations before we stop
 
          L = lag size of the Parzen window function. 
              If no value is given the lag size is set to
              be the lag where the auto correlation is less than 
              2 standard deviations. (maximum 200) 
        def = 'nonlinear' : transform based on smoothed crossing intensity (default)
              'mnonlinear': transform based on smoothed marginal distribution
              'linear'    : identity.
    options = options structure defining how the estimation of g is
              done, see troptset.
    
  In order to reconstruct the data a transformed Gaussian random process is
  used for modelling and simulation of the missing/removed data conditioned
  on the other known observations.
 
  Estimates of standarddeviations of y is obtained by a call to tranproc
         Std = tranproc(mu1o+/-mu1oStd,fliplr(g));
 
  See also  troptset, findoutliers, cov2csdat, dat2cov, dat2tr, detrendma

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [y,g,g2,test,tobs,mu1o, mu1oStd]=reconstruct(x,inds,Nsim,L,def,varargin)
0002 %RECONSTRUCT reconstruct the spurious/missing points of timeseries
0003 %
0004 % CALL: [y,g,g2,test,tobs,mu1o,mu1oStd]=reconstruct(x,inds,Nsim,L,def,options);
0005 %   
0006 %       y   = reconstructed signal
0007 %      g,g2 = smoothed and empirical transformation, respectively
0008 % test,tobs = test observator int(g(u)-u)^2 du and int(g_new(u)-g_old(u))^2 du,
0009 %             respectively, where int limits is given by param in lc2tr. 
0010 %             Test is a measure of departure from the Gaussian model for 
0011 %             the data. Tobs is a measure of the convergence of the 
0012 %             estimation of g.
0013 %      mu1o = expected surface elevation of the Gaussian model process.
0014 %   mu1oStd = standarddeviation of mu1o.
0015 %
0016 %       x   = 2 column timeseries 
0017 %             first column sampling times [sec]
0018 %             second column surface elevation [m]
0019 %      inds = indices to spurious points of x
0020 %      Nsim = the maximum # of iterations before we stop
0021 %
0022 %         L = lag size of the Parzen window function. 
0023 %             If no value is given the lag size is set to
0024 %             be the lag where the auto correlation is less than 
0025 %             2 standard deviations. (maximum 200) 
0026 %       def = 'nonlinear' : transform based on smoothed crossing intensity (default)
0027 %             'mnonlinear': transform based on smoothed marginal distribution
0028 %             'linear'    : identity.
0029 %   options = options structure defining how the estimation of g is
0030 %             done, see troptset.
0031 %   
0032 % In order to reconstruct the data a transformed Gaussian random process is
0033 % used for modelling and simulation of the missing/removed data conditioned
0034 % on the other known observations.
0035 %
0036 % Estimates of standarddeviations of y is obtained by a call to tranproc
0037 %        Std = tranproc(mu1o+/-mu1oStd,fliplr(g));
0038 %
0039 % See also  troptset, findoutliers, cov2csdat, dat2cov, dat2tr, detrendma
0040 
0041 % Reference
0042 % Brodtkorb, P, Myrhaug, D, and Rue, H (2001)
0043 % "Joint distribution of wave height and wave crest velocity from
0044 % reconstructed data with application to ringing"
0045 % Int. Journal of Offshore and Polar Engineering, Vol 11, No. 1, pp 23--32 
0046 %
0047 % Brodtkorb, P, Myrhaug, D, and Rue, H (1999)
0048 % "Joint distribution of wave height and wave crest velocity from
0049 % reconstructed data"
0050 % in Proceedings of 9th ISOPE Conference, Vol III, pp 66-73
0051 
0052 % tested on: Matlab 5.3, 5.1
0053 % History:
0054 % revised pab 18.04.2001 
0055 % - updated help header by adding def to function call
0056 % - fixed a bug: param was missing
0057 % revised pab 29.12.2000
0058 % - replaced csm1,...., param, with a options structure
0059 % - monitor replaced with plotflag==2
0060 % revised pab 09.10.2000
0061 % - updated call to dat2cov
0062 % revised pab 22.05.2000
0063 % - found a bug concerning cvmax and indr
0064 % - updated call to waveplot
0065 % revised pab 27.01.2000
0066 % - added L,csm1,csm2,param, monitor to the input arguments 
0067 % revised pab 17.12.1999
0068 % -added reference
0069 % revised pab 12.10.1999
0070 %   updated arg. list to dat2tr
0071 % last modified by Per A. Brodtkorb 01.10.98 
0072 
0073 opt = troptset('dat2tr');
0074 if nargin==1 & nargout <= 1 & isequal(x,'defaults')
0075   y = opt; 
0076   return
0077 end
0078 error(nargchk(1,inf,nargin))
0079 tic
0080 xn=x;% remember the old file
0081 [n m]= size(xn);
0082 if n<m
0083  b=m;m=n;n=b; 
0084  xn=xn.';
0085 end
0086 
0087 if n<2, 
0088   error('The vector must have more than 2 elements!')
0089 end
0090 
0091 switch m
0092  case 1, xn=[(1:n)' xn(:)];
0093  case 2, % dimension OK.
0094  otherwise, 
0095    error('Wrong dimension of input! dim must be 2xN, 1xN, Nx2 or Nx1 ')        
0096 end
0097 
0098 %%%%%%%%%%%%%%%%%%
0099 %                %
0100 %  initializing  %
0101 %                %
0102 %%%%%%%%%%%%%%%%%%
0103 if nargin<2|isempty(inds),  inds=isnan(xn(:,2));end
0104 if nargin<3|isempty(Nsim),  Nsim=20; end
0105 if nargin<4|isempty(L),  L=[]; end, % lagsize
0106 if nargin<5|isempty(def),  def='nonlinear';end
0107 if nargin>=6,  opt=troptset(opt,varargin{:});end
0108 
0109 param = getfield(opt,'param');
0110 
0111 switch opt.plotflag
0112   case {'none','off'},   plotflag = 0;
0113   case 'final', plotflag = 1;
0114   case 'iter',  plotflag = 2;
0115   otherwise,    plotflag = opt.plotflag;
0116 end
0117 
0118 clf
0119 olddef = def;
0120 method = 'approx';  %'approx';%'dec2'; % 'dec2' 'approx';
0121 ptime  = opt.delay; % pause for ptime sec if plotflag=2
0122 pause on
0123 expect1 = 1;     % first reconstruction by expectation? 1=yes 0=no
0124 expect  = 1;     % reconstruct by expectation? 1=yes 0=no
0125 tol     = 0.001; % absolute tolerance of e(g_new-g_old)
0126 
0127 cmvmax = 100; % if number of consecutive missing values (cmv) are longer they
0128              % are not used in estimation of g, due to the fact that the
0129              % conditional expectation approaches zero as the length to
0130              % the closest known points increases, see below in the for loop 
0131 
0132 dT=xn(2,1)-xn(1,1);
0133 Lm=min([n,200,floor(200/dT)]); % Lagmax 200 seconds
0134 if ~isempty(L), Lm=max([L,Lm]);end
0135 Lma = 1500;                    % size of the moving average window used
0136                                % for detrending the reconstructed signal 
0137 
0138 
0139 
0140 if any(inds>1),
0141   xn(inds,2)=NaN;
0142   inds=isnan(xn(:,2));
0143 elseif sum(inds)==0,
0144   error('No spurious data given')
0145 end
0146 endpos  = diff([ inds ]);
0147 strtpos = find(endpos>0);
0148 endpos  = find(endpos<0);
0149 
0150 indg=find(~inds); % indices to good points
0151 inds=find(inds);  % indices to spurous points
0152 
0153 
0154 indNaN = []; % indices to points omitted in the covariance estimation
0155 indr  = 1:n; % indices to point used in the estimation of g
0156 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0157 %
0158 % Finding more than cmvmax consecutive spurios points. 
0159 % They will not be used in the estimation of g and are thus removed 
0160 % from indr.
0161 %
0162 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0163 if ~isempty(strtpos) & (isempty(endpos)|endpos(end)<strtpos(end)),
0164   if (n-strtpos(end))>cmvmax 
0165     indNaN= indr(strtpos(end)+1:n);
0166     indr(strtpos(end)+1:n)=[];
0167   end
0168   strtpos(end)=[];
0169 end
0170 if ~isempty(endpos) & (isempty(strtpos)|(endpos(1)<strtpos(1))),
0171   if (endpos(1))>cmvmax 
0172     indNaN=[indNaN,indr(1:endpos(1))];
0173     indr(1:endpos(1))=[];
0174   end
0175   strtpos   = strtpos-endpos(1);
0176   endpos    = endpos-endpos(1);
0177   endpos(1) = [];
0178 end
0179 %length(endpos)
0180 %length(strtpos)
0181 for ix=length(strtpos):-1:1
0182   if (endpos(ix)-strtpos(ix)>cmvmax)
0183     indNaN=[indNaN, indr(strtpos(ix)+1:endpos(ix))];
0184     indr(strtpos(ix)+1:endpos(ix))=[]; % remove this when estimating the transform
0185   end
0186 end
0187 if length(indr)<0.1*n,
0188   error('Not possible to reconstruct signal')
0189 end
0190 
0191 if any(indNaN),
0192   indNaN=sort(indNaN);
0193 end
0194 
0195 switch 1, % initial reconstruction attempt
0196 case  0,% spline
0197   xn(:,2) = interp1(xn(indg,1),xn(indg,2),xn(:,1),'*spline'); 
0198   y=xn;
0199   return
0200   
0201 case 1,% 
0202  % xn(indg,2)=detrendma(xn(indg,2),1500);
0203   
0204   [g test cmax irr g2]  = dat2tr(xn(indg,:),def,opt);
0205   xnt=xn;
0206   xnt(indg,:)=dat2gaus(xn(indg,:),g);
0207   xnt(inds,2)=NaN;
0208   rwin=findrwin(xnt,Lm,L);
0209   disp(['First reconstruction attempt,    e(g-u)=', num2str(test)] )
0210   [samp ,mu1o, mu1oStd]  = cov2csdat(xnt(:,2),rwin,1,method,inds); % old simcgauss
0211   if expect1,% reconstruction by expectation
0212     xnt(inds,2) =mu1o;
0213   else
0214     xnt(inds,2) =samp;
0215   end
0216   xn=gaus2dat(xnt,g);
0217   xn(:,2)=detrendma(xn(:,2),Lma); % detrends the signal with a moving
0218                                    % average of size Lma
0219   g_old=g;
0220   
0221 end
0222 
0223 bias = mean(xn(:,2));
0224 xn(:,2)=xn(:,2)-bias; % bias correction
0225 
0226 if plotflag==2
0227   clf
0228   mind=1:min(1500,n);
0229   waveplot(xn(mind,:),x(inds(mind),:), 6,1)
0230   subplot(111)
0231   pause(ptime)
0232 end
0233 
0234 test0=0;
0235 for ix=1:Nsim,
0236   if 0,%ix==2,
0237     rwin=findrwin(xn,Lm,L);
0238     xs=cov2sdat(rwin,[n 100 dT]);
0239     [g0 test0 cmax irr g2]  = dat2tr(xs,def,opt);
0240     [test0 ind0]=sort(test0);
0241   end
0242   
0243    if 1, %test>test0(end-5),
0244      % 95% sure the data comes from a non-Gaussian process
0245      def = olddef; %Non Gaussian process
0246    else
0247      def = 'linear'; % Gaussian process
0248    end
0249    % used for isope article
0250    % indr =[1:27000 30000:39000];
0251    % Too many consecutive missing values will influence the estimation of
0252    % g. By default do not use consecutive missing values if there are more 
0253    % than cmvmax. 
0254    
0255    [g test cmax irr g2]  = dat2tr(xn(indr,:),def,opt);
0256   if plotflag==2,
0257     pause(ptime)
0258   end
0259   
0260   
0261   %tobs=sqrt((param(2)-param(1))/(param(3)-1)*sum((g_old(:,2)-g(:,2)).^2))
0262   % new call
0263   tobs=sqrt((param(2)-param(1))/(param(3)-1)....
0264     *sum((g(:,2)-interp1(g_old(:,1)-bias, g_old(:,2),g(:,1),'spline')).^2));
0265   
0266   if ix>1 
0267     if tol>tobs2 & tol>tobs,    
0268       break, %estimation of g converged break out of for loop    
0269     end
0270   end
0271  
0272   tobs2=tobs;
0273   
0274   xnt=dat2gaus(xn,g);
0275   if ~isempty(indNaN),    xnt(indNaN,2)=NaN;  end
0276   rwin=findrwin(xnt,Lm,L);    
0277   disp(['Simulation nr: ', int2str(ix), ' of ' num2str(Nsim),'   e(g-g_old)=', num2str(tobs), ',  e(g-u)=', num2str(test)])
0278   [samp ,mu1o, mu1oStd]  = cov2csdat(xnt(:,2),rwin,1,method,inds);
0279   
0280   if expect,
0281     xnt(inds,2) =mu1o;
0282   else
0283     xnt(inds,2) =samp;
0284   end
0285   
0286   xn=gaus2dat(xnt,g);
0287   if ix<Nsim
0288     bias=mean(xn(:,2));
0289     xn(:,2) = (xn(:,2)-bias); % bias correction
0290   end
0291   g_old=g;% saving the last transform
0292   if plotflag==2
0293     waveplot(xn(mind,:),x(inds(mind),:),6,1,[])
0294     subplot(111)
0295     pause(ptime)
0296   end
0297 end % for loop
0298 
0299 if 1, %test>test0(end-5) 
0300   xnt=dat2gaus(xn,g);
0301   [samp ,mu1o, mu1oStd]  = cov2csdat(xnt(:,2),rwin,1,method,inds);
0302   xnt(inds,2) =samp;
0303   xn=gaus2dat(xnt,g);
0304   bias=mean(xn(:,2));
0305   xn(:,2) = (xn(:,2)-bias); % bias correction
0306   g(:,1)=g(:,1)-bias;
0307   g2(:,1)=g2(:,1)-bias;
0308   gn=trangood(g);
0309  
0310   %mu1o=mu1o-tranproc(bias,gn);
0311   muUStd=tranproc(mu1o+2*mu1oStd,fliplr(gn));%
0312   muLStd=tranproc(mu1o-2*mu1oStd,fliplr(gn));%
0313 else
0314   muLStd=mu1o-2*mu1oStd;
0315   muUStd=mu1o+2*mu1oStd;
0316 end
0317 
0318 if  plotflag==2 & length(xn)<10000,
0319   waveplot(xn,[xn(inds,1) muLStd ;xn(inds,1) muUStd ], 6,round(n/3000),[])
0320   legend('reconstructed','2 stdev')
0321   %axis([770 850 -1 1])
0322   %axis([1300 1325 -1 1])
0323 end
0324 y=xn;
0325 toc
0326 
0327 return
0328 
0329 function r=findrwin(xnt,Lm,L)
0330   r=dat2cov(xnt,Lm);%computes  ACF
0331   %finding where ACF is less than 2 st. deviations .
0332   % in order to find a better L  value
0333   if nargin<3|isempty(L)
0334     L=find(abs(r.R)>2*r.stdev)+1;
0335     if isempty(L), % pab added this check 09.10.2000
0336       L = Lm;
0337     else
0338       L = min([floor(4/3*L(end)) Lm]);
0339     end
0340   end
0341   win=parzen(2*L-1);
0342   r.R(1:L)=win(L:2*L-1).*r.R(1:L);
0343   r.R(L+1:end)=0;
0344   return
0345

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group


Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003