Home > wafo > multidim > dat2dspec.m

dat2dspec

PURPOSE ^

Estimates the directional wave spectrum from timeseries

SYNOPSIS ^

[Sd,D,Sw,Fcof,Gwt,Sxy,Sxy1] = dat2dspec2(xn,pos,h,nfft,nt,method,varargin)

DESCRIPTION ^

 DAT2DSPEC Estimates the directional wave spectrum from timeseries  
   
   CALL: [S, D, Sw,Fcof] = dat2dspec(W,pos,h,Nfft,Nt,method,options); 
     
          S = a spectral density structure containing: 
              S      = 2D array of the estimated directional spectrum, size Nt x Nf.  
              w      = frequency vector  0..2*pi*Nyquistfrequency of length Nf 
              theta  = angle vector -pi..pi of length Nt  
                       (theta = 0 -> + x-axis, theta = pi/2 -> + y-axis)  
              h      = water depth  
                       ( S.S=D.S.*Sw.S(:,ones(:,Nt))' ) 
         D  = estimate of spreading function as function of theta and w 
         Sw = angular frequency spectrum 
       Fcof = Fourier coefficients of the spreading function D(w,theta) 
              defined as int D(w,theta)*exp(i*n*theta) dtheta. n=0:nharm 
              size nharm+1 x Nf 
  
         W  = [t w1 w2 ... wM]  a M+1 column data matrix with sampled times  
              and values in column 1 and column 2 to M+1, respectively. 
        pos = [x y z def bfs], matrix characterizing instruments and postions, size M x 5 
              x(j), y(j),z(j) = the coordinate positon of the j-th sensor (j=1:M) 
              def(j) = identifying parameter (1,2,... or 18) for the j-th time series   
                       (See tran for options) 
              bfs(j) = 1, if j-th time series is to be used in final estimate of  
                          "Best" Frequency Spectra. 
                       0, if j-th is to be excluded. (At least one must be set to 1) 
         h  = water depth (default is deep water,i.e., h = inf) 
       Nfft = length of FFT  (determines the smoothing and number of 
              frequencies, Nf = Nfft/2+1) (default 256) 
         Nt = number of angles (default 101) 
     method = 'MLM'  Maximum Likelihood Method (default) 
              'IMLM' Iterative  Maximum Likelihood Method  
              'MEM'  Maximum Entropy Method   (slow) 
              'EMEM' Extended Maximum Entropy Method 
    options = optional options, see specoptset for details. 
  
  NB! requires psd and csd from the signal toolbox 
  
  See also  specoptset, psd, csd, tran, w2k, wspecplot

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [Sd,D,Sw,Fcof,Gwt,Sxy,Sxy1] = dat2dspec2(xn,pos,h,nfft,nt,method,varargin) 
0002 %DAT2DSPEC Estimates the directional wave spectrum from timeseries  
0003 %  
0004 %  CALL: [S, D, Sw,Fcof] = dat2dspec(W,pos,h,Nfft,Nt,method,options); 
0005 %    
0006 %         S = a spectral density structure containing: 
0007 %             S      = 2D array of the estimated directional spectrum, size Nt x Nf.  
0008 %             w      = frequency vector  0..2*pi*Nyquistfrequency of length Nf 
0009 %             theta  = angle vector -pi..pi of length Nt  
0010 %                      (theta = 0 -> + x-axis, theta = pi/2 -> + y-axis)  
0011 %             h      = water depth  
0012 %                      ( S.S=D.S.*Sw.S(:,ones(:,Nt))' ) 
0013 %        D  = estimate of spreading function as function of theta and w 
0014 %        Sw = angular frequency spectrum 
0015 %      Fcof = Fourier coefficients of the spreading function D(w,theta) 
0016 %             defined as int D(w,theta)*exp(i*n*theta) dtheta. n=0:nharm 
0017 %             size nharm+1 x Nf 
0018 % 
0019 %        W  = [t w1 w2 ... wM]  a M+1 column data matrix with sampled times  
0020 %             and values in column 1 and column 2 to M+1, respectively. 
0021 %       pos = [x y z def bfs], matrix characterizing instruments and postions, size M x 5 
0022 %             x(j), y(j),z(j) = the coordinate positon of the j-th sensor (j=1:M) 
0023 %             def(j) = identifying parameter (1,2,... or 18) for the j-th time series   
0024 %                      (See tran for options) 
0025 %             bfs(j) = 1, if j-th time series is to be used in final estimate of  
0026 %                         "Best" Frequency Spectra. 
0027 %                      0, if j-th is to be excluded. (At least one must be set to 1) 
0028 %        h  = water depth (default is deep water,i.e., h = inf) 
0029 %      Nfft = length of FFT  (determines the smoothing and number of 
0030 %             frequencies, Nf = Nfft/2+1) (default 256) 
0031 %        Nt = number of angles (default 101) 
0032 %    method = 'MLM'  Maximum Likelihood Method (default) 
0033 %             'IMLM' Iterative  Maximum Likelihood Method  
0034 %             'MEM'  Maximum Entropy Method   (slow) 
0035 %             'EMEM' Extended Maximum Entropy Method 
0036 %   options = optional options, see specoptset for details. 
0037 % 
0038 % NB! requires psd and csd from the signal toolbox 
0039 % 
0040 % See also  specoptset, psd, csd, tran, w2k, wspecplot 
0041  
0042 %             freq = frequency vector 0..Nyquistfrequency 
0043  
0044 % References: 
0045 % Isobe M, Kondo K, Horikawa K. (1984) 
0046 % "Extension of MLM for estimating directional wave spectrum", 
0047 % In Symposium on Description and modelling of Directional Seas, 
0048 % Copenhagen, pp 1-15 
0049 % 
0050 % Young I.R. (1994) 
0051 % "On the measurement of directional spectra", 
0052 % Applied Ocean Research, Vol 16, pp 283-294 
0053 % 
0054 % Hashimoto,N. (1997)  
0055 % "Analysis of the directional wave spectra from field data.", Advances in  
0056 % Coastal and Ocean Engineering, Vol.3., pp.103-143 
0057 % 
0058 % Massel S.R. and Brinkman R. M. (1998) 
0059 % "On the determination od directional wave spectra  
0060 %   for practical applications", 
0061 %  Applied Ocean Research, Vol 20, pp 357-374 
0062 % 
0063 % Michel K. Ochi (1998), 
0064 % "OCEAN WAVES, The stochastic approach", 
0065 %  OCEAN TECHNOLOGY series 6, Cambridge, chapter 7. 
0066  
0067  
0068 % Tested on: Matlab 5.3, 5.2  
0069 % History: 
0070 % revised pab dec2003   
0071 % revised pab 08.11.2002 
0072 % -changed name to dat2dspec2 in order to incorprate the specoptset option. 
0073 % revised pab 24.10.2002 
0074 % -made it more robust by estimating Hw from data if possible. 
0075 % -added method EMEM and revised the MEM code 
0076 % revised pab 09.05.2001  
0077 %   Due to Vengatesan Venugopal plotting is now handled correctly. 
0078 % revised pab 08.05.2001 
0079 % - added optimset to fsolve for Matlab v 5.3 and higher. 
0080 % revised pab 21.06.2000 
0081 %   - added 'MEM' 
0082 %   - added 'IMLM' and checked with testsurf and testbuoy 
0083 %   - replaced inv with pinv 
0084 %   - changed order of input arguments 
0085 %   - added dflag, ftype to input 
0086 % revised pab 06.01.2000  
0087 %   - tested  narrowbanded cases generated with testsurf and testbuoy -> OK 
0088 %   - changed size of S.S from Nf x Nt   to     Nt x Nf 
0089 %   - added the possibility that W contains timeseries of other  
0090 %     quantities than a surface elevation 
0091 %   - added par 
0092 %   - extended pos with def and bfs  
0093 % revised pab 28.08.99 
0094 %  updated spectral structure 
0095 % By Per A. Brodtkorb 27.03.99 
0096  
0097 % TODO % Add option 'fem'. It is not implemented 
0098 % TODO % measurements should be scaled       
0099 % TODO % all options in specoptset should be implemented  
0100 % TODO % Remove dependence on psd and csd in signal toolbox. 
0101    
0102 opt = specoptset('plotflag','off','dflag','mean','ftype','w',.... 
0103          'thtype','r','nharm',2,'delay',0,'gravity',[],'wdensity',[],... 
0104          'bet',[],'igam',[],'x_axisdir',[],'y_axisdir',[],... 
0105          'maxIter',25,'coefAbsTol',0.01,'errorTol',[],... 
0106          'minModelOrder',1,'relax',[]); 
0107           
0108 % If just 'defaults' passed in, return the default options in Sd 
0109 if nargin==1 & nargout <= 1 & isequal(xn,'defaults') 
0110   Sd = opt;  
0111   return 
0112 end 
0113 error(nargchk(2,inf,nargin))  
0114  
0115  
0116 if nargin>6,  opt  = specoptset(opt,varargin{:}); end 
0117  
0118  
0119 xx = xn; 
0120 if nargin<2, 
0121   error('Must have 2 inputs arguments') 
0122 end 
0123 [n m]= size(xx); 
0124 if n<m, b=m;m=n;n=b; xx=xx'; end 
0125  
0126 if n<4,   error('The vector must have more than 4 elements!'),end 
0127  
0128 switch m 
0129   case {1,2,3} , 
0130     error('Wrong dimension of input! Must at least have 4 columns! ')    
0131  otherwise,   % dimension OK!      
0132 end 
0133 if any([m-1 5]~=size(pos)),   
0134   error('Wrong dimension of pos, must be of size MX5'), 
0135 end 
0136  
0137 Fs = 1/(xx(2,1)-xx(1,1));% sampling frequency 
0138  
0139 % Default values 
0140 %~~~~~~~~~~~~~~~ 
0141 if nargin<3|isempty(h),     h      = inf; end                 % deep water is default 
0142 if nargin<4|isempty(nfft),  nfft   = min(256,n-mod(n,2)); end %make sure nfft<=n and even 
0143 if nargin<5|isempty(nt),    nt     = 101;end 
0144 if nargin<6|isempty(method),method = 'mlm';end 
0145  
0146 nfft = min(nfft+mod(nfft,2),n-mod(n,2)); 
0147  
0148 dflag  = 'mean';   %'linear'; %detrending option 'none','mean','ma' 
0149 ftype  = 'w';    % options are f=S(f,phi),w=S(w,phi) 
0150 nharm  = 2;      % number of harmonics 
0151 par = [];g=[];rho=[];bet=[];igam=[];thx=[];thy=[]; % use default values in tran 
0152 if nargout==0, plotflag = 1;  else,    plotflag = 0;  end 
0153 noverlap = floor(nfft/2); 
0154  
0155 switch opt.dflag; 
0156   case {'mean','ma','linear','none'}, dflag =opt.dflag; 
0157 end 
0158 switch opt.ftype; 
0159   case {'w','f'}, ftype = opt.ftype; 
0160 end 
0161 switch opt.thtype; 
0162   case {'r','d'}, thtype = opt.thtype; 
0163 end 
0164 switch opt.plotflag 
0165   case {'none','off'},   plotflag = 0; 
0166   case 'final', plotflag = 1; 
0167   case 'iter',  plotflag = 2; 
0168     otherwise plotflag = opt.plotflag; 
0169 end 
0170 if ~isempty(opt.noverlap),  noverlap = opt.noverlap;end 
0171 if ~isempty(opt.window),    window   = opt.window; else window = nfft; end 
0172 if ~isempty(opt.message),   message  = opt.message; end 
0173 if ~isempty(opt.delay),     delay    = opt.delay; end 
0174 if ~isempty(opt.nharm),     nharm    = opt.nharm; end 
0175 if ~isempty(opt.gravity),   g        = opt.gravity; end 
0176 if ~isempty(opt.wdensity),  rho      = opt.wdensity; end 
0177 if ~isempty(opt.bet),       bet      = opt.bet; end 
0178 if ~isempty(opt.igam),      igam     = opt.igam; end 
0179 if ~isempty(opt.x_axisdir), thx      = opt.x_axisdir;end 
0180 if ~isempty(opt.y_axisdir), thy      = opt.y_axisdir; end 
0181  
0182  
0183 if nt<30,  warning(' # of angles are low =>  spectral density may be inaccurate'), end 
0184 if nharm>=nt,  
0185   nharm = nt-1; 
0186   disp('Nharm must be less than Nt') 
0187 end 
0188  
0189 bfs = find(pos(:,5)>.5); % indicator if it should be used to estimate spectra 
0190 if isempty(bfs),  error('Not all the elements of bfs can be zero'),end 
0191  
0192 %thetai=(linspace(0,2*pi,nt)); 
0193 thetai   = linspace(-pi,pi,nt)'; 
0194 nf       = nfft/2+1; 
0195  
0196  
0197 %-------------------------------------------- 
0198 % Detrend the measurements before estimation 
0199 %-------------------------------------------- 
0200 switch lower(dflag) 
0201 case 'none', % do nothing 
0202 case 'mean',   mn        = mean(xx(:,2:m)); 
0203                xx(:,2:m) = (xx(:,2:m)-mn(ones(n,1),:));% remove mean 
0204 case 'linear', xx(:,2:m) = detrend(xx(:,2:m));         % remove linear trend 
0205 case 'ma',     xx(:,2:m) = detrendma(xx(:,2:m),2*nfft);% remove moving average 
0206                dflag     = 'linear'; % must change to a valid option for psd and csd 
0207 end 
0208  
0209  
0210  
0211 % Scale measurements so that they have equal standard deviation 
0212 % and zero mean. 
0213 % This is done in order to overcome possible calibration errors 
0214 % between the different measuring devices. 
0215  
0216 stdev = std(xx(:,2:m));  
0217 mn    = mean(xx(:,2:m)); 
0218  
0219 normfact = ones(1,m-1); 
0220 def = unique(pos(:,4).'); 
0221 for ix=def 
0222   k = find(pos(:,4)==ix); 
0223   if any(k) 
0224     normfact(k) = stdev(k)/sqrt(mean(stdev(k).^2));  
0225   end 
0226 end 
0227  
0228 xx(:,2:m) = (xx(:,2:m)-mn(ones(n,1),:))./normfact(ones(n,1),:);% remove mean 
0229  
0230  
0231 %---------------------------------------------------------------------- 
0232 % finding the spectra, matrix of cross-spectra and transfer functions 
0233 %---------------------------------------------------------------------- 
0234 Hw  = zeros(m-1,nf);     % a function of frequency only 
0235 Gwt = zeros(m-1,nt,nf);  % a function of frequency and direction combined. 
0236 EE  = Gwt;              
0237 Sf  = zeros(m-1,nf).'; 
0238 Sxy = zeros(m-1,m-1,nf); 
0239  
0240 kw=[]; 
0241 %kw=w2k(2*pi*fi,0,h); 
0242  
0243 for ix=1:m-1, 
0244   [Sf(:,ix), fi] = psd(xx(:,ix+1),nfft,Fs,window,noverlap,dflag); 
0245   [Hw(ix,:),Gwt(ix,:,:),kw]=tran(2*pi*fi,thetai,pos(ix,1:3),pos(ix,4),h,g,rho,bet,igam,thx,thy,kw); 
0246   Sxy(ix,ix,:)  = Sf(:,ix).'; 
0247   for iy=(ix+1):m-1, 
0248     Sxy(ix,iy,:) = csd(xx(:,ix+1),xx(:,iy+1),nfft,Fs,window,noverlap,dflag); 
0249     Sxy(iy,ix,:) = conj(Sxy(ix,iy,:)); 
0250   end 
0251 end 
0252 if nargout>6,  Sxy1 = Sxy;end 
0253  
0254 Sf  = Sf*2/Fs;         % make sure Sf and Sxy have the correct units 
0255 Sxy = Sxy*2/Fs; 
0256  
0257 Sf = Sf.'; 
0258 kw = kw(:)'; 
0259  
0260 %------------------------------------------------------------------------------- 
0261 %Estimate frequency spectrum for the surface elevation from the bfs timeseries 
0262 %-------------------------------------------------------------------------------- 
0263 SfBest = bfsspec(Sf,Hw,pos,bfs); 
0264 %--------------------------------------------------------------------------------- 
0265 %Estimate the absolute value of the transfer function H(w) from the sensor spectra 
0266 %--------------------------------------------------------------------------------- 
0267 Hw = hwestimate(Sf,SfBest,Hw,pos); 
0268  
0269  
0270 %------------------------------------------ 
0271 %  Normalizing the matrix of cross-spectra 
0272 %----------------------------------------- 
0273  
0274 k = find(SfBest); % avoid division by zero 
0275 Sxyn = Sxy; 
0276 for ix=1:m-1 
0277  Sxyn(ix,ix,k) = squeeze(Sxy(ix,ix,k)).'./(SfBest(k).*Hw(ix,k).^2); %1 
0278  %Sxyn(ix,ix,k) = squeeze(Sxy(ix,ix,k)).'./(Hw(ix,k).^2); %2 
0279  %Sxyn(ix,ix,k) = 1; %3 
0280  %Sxx =  squeeze(Sxy(ix,ix,k)).'; %3 
0281   for iy=(ix+1):m-1, 
0282     Sxyn(ix,iy,k) = squeeze(Sxy(ix,iy,k)).'./(SfBest(k).*Hw(ix,k).*Hw(iy,k)); % 1 
0283     %Sxyn(ix,iy,k) = squeeze(Sxy(ix,iy,k)).'./(Hw(ix,k).*Hw(iy,k)); % 2 
0284     %Syy =  squeeze(Sxy(iy,iy,k)).'; %3 
0285     %Sxyn(ix,iy,k) = squeeze(Sxy(ix,iy,k)).'./(sqrt(Sxx.*Syy)); % 3 
0286     Sxyn(iy,ix,:) = conj(Sxyn(ix,iy,:)); 
0287   end   
0288 end 
0289  
0290 switch lower(method) 
0291   case {'mlm'} % maximum likelihood method 
0292    DS = mlm(Sxyn,Gwt,thetai,fi,k,opt); 
0293  case {'imlm'} % iterative maximum likelihood method 
0294    DS = imlm(Sxyn,Gwt,thetai,fi,k,opt); 
0295   case 'mem' , %Maximum entropy method 
0296     DS = mem(Sxyn,Gwt,thetai,fi,k,opt); 
0297   case 'emem' , %Extended Maximum entropy method       
0298     DS = emem(Sxyn,Gwt,thetai,fi,k,opt); 
0299   case 'bdm', %    
0300     DS = bdm(Sxyn,Gwt,thetai,fi,k,opt);   
0301   case 'fem', % Fourier Expansion method 
0302  
0303   error('FEM method not implemented yet') 
0304 otherwise , error('Unknown method') 
0305 end 
0306  
0307 if 0,%used for debugging 
0308   Sxy2 = getcrossspectra(thetai,Gwt,DS); 
0309   for ix = 1:m-1, 
0310     for iy=ix:m-1, 
0311       %subplot(2,1,1), plot(fi(k),squeeze(real(Sxy2(ix,iy,k))),'r',fi(k),squeeze(real(Sxyn(ix,iy,k))),'g') 
0312       %subplot(2,1,2), plot(fi(k),squeeze(real(Sxy2(ix,iy,k))),'r',fi(k),squeeze(real(Sxyn(ix,iy,k))),'g') 
0313       subplot(2,1,1), semilogy(fi(k),eps+abs(squeeze(real(Sxy2(ix,iy,k)))-squeeze(real(Sxyn(ix,iy,k)))),'g') 
0314       subplot(2,1,2), semilogy(fi(k),eps+abs(squeeze(real(Sxy2(ix,iy,k)))-squeeze(real(Sxyn(ix,iy,k)))),'g') 
0315       disp('hit any key'),pause, disp('hit any key') 
0316     end 
0317   end 
0318 end 
0319  
0320 Sd       = createspec('dir',ftype); 
0321 Sd.theta = thetai(:); 
0322 Sd.h     = h; 
0323 Sd.norm  = 0; 
0324 Sd.note  = ['dat2dspec(' inputname(1) ')']; 
0325  
0326 if strcmp(ftype,'w'), 
0327   fi   = fi*2*pi; 
0328   SfBest   = SfBest/2/pi; 
0329   Sd.w = fi(:); 
0330 else 
0331   Sd.f = fi(:); 
0332 end 
0333 Sd.S = (SfBest(ones(nt,1),:).*DS); 
0334 Sd = ttspec(Sd,thtype); 
0335  
0336 if nargout>1 
0337  D      = Sd; 
0338  D.S    = DS; 
0339  D.note = [D.note ', D(theta,w)']; 
0340  
0341 end 
0342  
0343 if nargout>2 
0344  Sw      = createspec('freq',ftype); 
0345  Sw.S    = SfBest(:); 
0346  Sw.note = [D.note ', S(w)']; 
0347  Sw      = setfield(Sw,ftype,fi(:));  
0348 end 
0349  
0350 if nargout>3 
0351   %   int D(w,theta)*exp(i*n*theta) dtheta. 
0352   if 0, 
0353     Fcof1 = 2*pi*ifft(DS,[],1); % integration by FFT 
0354     Pcor = [1; exp(sqrt(-1)*[1:nharm].'*thetai(1))]; % correction term to get 
0355                                                      % the correct integration limits 
0356     Fcof = Fcof1(1:1+nharm,:).*Pcor(:,ones(1,nf)); 
0357   else 
0358     [a,b]=fourier(thetai,DS,2*pi,nharm); 
0359     Fcof = (a+sqrt(-1)*b)*pi; 
0360   end 
0361 end 
0362  
0363 %-----------------------------------  
0364 %  
0365 %     Plotting the Spectral Density  
0366 % 
0367 %----------------------------------- 
0368 if plotflag>0, 
0369  wspecplot(Sd,plotflag)  
0370 end 
0371 return   
0372  
0373  
0374  
0375  
0376  
0377  
0378  
0379  
0380  
0381  
0382  
0383  
0384  
0385  
0386  
0387

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