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

## CROSS-REFERENCE INFORMATION

This function calls:
 bfsspec Estimate frequency spectrum for the surface elevation from the bfs timeseries createspec Spectrum structure constructor detrendma Removes a trend from data using a moving average emem Extended Maximum Entropy Method fourier Returns Fourier coefficients. getcrossspectra Compute the cross spectra by integration hwestimate Estimate absolute value of transfer function H(w) from sensor spectra imlm Iterated maximum likelihood method for estimating the directional distribution mem maximum entropy method for estimating the directional distribution mlm maximum likelihood method for estimating the directional distribution specoptset Create or alter SPECTRUM OPTIONS structure. tran Computes transfer functions based on linear wave theory ttspec Toggle Transform between angular frequency and frequency spectrum wspecplot Plot a spectral density bdm csd Cross Spectral Density estimate. detrend Remove a linear trend from a vector, usually for FFT processing. error Display message and abort function. ifft Inverse discrete Fourier transform. inputname Input argument name. isequal True if arrays are numerically equal. linspace Linearly spaced vector. lower Convert string to lowercase. mean Average or mean value. pause Wait for user response. psd Power Spectral Density estimate. semilogy Semi-log scale plot. setfield Set structure field contents. squeeze Remove singleton dimensions. std Standard deviation. strcmp Compare strings. subplot Create axes in tiled positions. unique Set unique. warning Display warning message; disable or enable warning messages.
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 %
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
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
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