dat2dspec

 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

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