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
generates conditionally simulated values | |
Simulates a Gaussian process and its derivative | |
Estimate auto covariance function from data. | |
Transforms x using the transformation g. | |
Estimate transformation, g, from data. | |
Removes a trend from data using a moving average | |
Transforms xx using the inverse of transformation g. | |
returns the N-point Parzen window in a column vector. | |
Makes a transformation that is suitable for efficient transforms. | |
Transforms process X and up to four derivatives | |
Create or alter TRANSFORM OPTIONS structure. | |
Plots the surface elevation of timeseries. | |
Clear current figure. | |
Difference and approximate derivative. | |
Display message and abort function. | |
Get structure field contents. | |
Convert integer to string (Fast version). | |
1-D interpolation (table lookup) | |
True if arrays are numerically equal. | |
Display legend. | |
Average or mean value. | |
Not-a-Number. | |
Convert number to string. (Fast version) | |
Wait for user response. | |
Create axes in tiled positions. | |
Start a stopwatch timer. | |
Read the stopwatch timer. |
% CHAPTER2 Modelling random loads and stochastic waves | |
% CHAPTER3 Demonstrates distributions of wave characteristics | |
setup all global variables of the RECDEMO |
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
Comments or corrections to the WAFO group