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
Estimate frequency spectrum for the surface elevation from the bfs timeseries | |
Spectrum structure constructor | |
Removes a trend from data using a moving average | |
Extended Maximum Entropy Method | |
Returns Fourier coefficients. | |
Compute the cross spectra by integration | |
Estimate absolute value of transfer function H(w) from sensor spectra | |
Iterated maximum likelihood method for estimating the directional distribution | |
maximum entropy method for estimating the directional distribution | |
maximum likelihood method for estimating the directional distribution | |
Create or alter SPECTRUM OPTIONS structure. | |
Computes transfer functions based on linear wave theory | |
Toggle Transform between angular frequency and frequency spectrum | |
Plot a spectral density | |
Cross Spectral Density estimate. | |
Remove a linear trend from a vector, usually for FFT processing. | |
Display message and abort function. | |
Inverse discrete Fourier transform. | |
Input argument name. | |
True if arrays are numerically equal. | |
Linearly spaced vector. | |
Convert string to lowercase. | |
Average or mean value. | |
Wait for user response. | |
Power Spectral Density estimate. | |
Semi-log scale plot. | |
Set structure field contents. | |
Remove singleton dimensions. | |
Standard deviation. | |
Compare strings. | |
Create axes in tiled positions. | |
Set unique. | |
Display warning message; disable or enable warning messages. |
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
Comments or corrections to the WAFO group