IMLM Iterated maximum likelihood method for estimating the directional distribution CALL DS = imlm(Sxy,Gwt,thetai,fi,k); DS = Directional distribution (spreading function) size nt x nf Sxy = matrix of cross spectral densities size m x m x nf Gwt = matrix of transfer function (abs(Gwt)==1) size m x nt x nf thetai = angle vector length nt fi = frequency vector length nf k = index vector to frequencies where Sf>0 length <= nf (m = number of measurement devices) nf = number frequencies (f or w) nt = number of angles (theta)
Spectrum structure constructor | |
maximum likelihood method for estimating the directional distribution | |
normalizes the spreading function | |
Numerical integration with the Simpson method | |
Plot a spectral density | |
Close figure. | |
Convert number to string. (Fast version) | |
Wait for user response. | |
Remove singleton dimensions. | |
Display wait bar. |
Estimates the directional wave spectrum from timeseries |
001 function DS = imlm(Sxy,Gwt,thetai,fi,k,opt) 002 %IMLM Iterated maximum likelihood method for estimating the directional distribution 003 % 004 % CALL DS = imlm(Sxy,Gwt,thetai,fi,k); 005 % 006 % DS = Directional distribution (spreading function) size nt x nf 007 % Sxy = matrix of cross spectral densities size m x m x nf 008 % Gwt = matrix of transfer function (abs(Gwt)==1) size m x nt x nf 009 % thetai = angle vector length nt 010 % fi = frequency vector length nf 011 % k = index vector to frequencies where Sf>0 length <= nf 012 % 013 % (m = number of measurement devices) 014 % nf = number frequencies (f or w) 015 % nt = number of angles (theta) 016 % 017 018 019 % revised pab jan2005 020 % Added opt to input 021 [m,nt,nf] = size(Gwt); 022 023 DS = mlm(Sxy,Gwt,thetai,fi,k); 024 025 DS0 = DS; 026 DSold = DS; 027 % Parameters controlling the the convergence 028 % Li = relaxation parameter (0< Li<=1.5) 1..1.2 is proposed by Krogstad 029 % Bi = exponent ( 0< Bi) Bi = 1..3 seems appropriate 030 Li = 1.4; 031 Bi = 2; 032 errorTol = 1e-3; % maximum tolerance 033 maxIter = 30; % maximum number of iterations 034 display =0; 035 if nargin>5 036 if ~isempty(opt.errortol), errorTol = opt.errortol; end 037 if ~isempty(opt.maxiter), maxIter = opt.maxiter; end 038 if ~isempty(opt.relax), Li = opt.relax; end 039 if ~isempty(opt.message), display = opt.message; end 040 end 041 tolold = inf; 042 043 h = waitbar(0,'Please wait...IMLM calculation'); 044 for iz = 1:maxIter 045 waitbar(iz/maxIter,h) 046 % Calculation of cross spectra based on DS 047 for ix=1:m 048 Sxy(ix,ix,:) = simpson(thetai,squeeze(Gwt(ix,:,:).*conj(Gwt(ix,:,:))).*DS); 049 for iy=(ix+1):m, 050 Sxy(ix,iy,:) = simpson(thetai,squeeze(Gwt(ix,:,:).*conj(Gwt(iy,:,:))).*DS); 051 Sxy(iy,ix,:) = conj(Sxy(ix,iy,:)); 052 end 053 end 054 tmp = (DS0-mlm(Sxy,Gwt,thetai,fi,k)); 055 DS = DS+Li*sign(tmp).*abs(tmp.^Bi); 056 %Normalize so that int D(theta,f) dtheta = 1 for each f 057 %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 058 DS = normspfn(DS,thetai); 059 tol = max(abs(DS(:)-DSold(:))); 060 disp(['Iteration nr ' num2str(iz),' of ' num2str(maxIter),' Error = ' num2str(tol) ]) 061 if iz>5, 062 if (min(tol,abs(tol-tolold)*3) < errorTol) , 063 disp('Close enough to convergence'),break, 064 end 065 end 066 tolold = tol; 067 DSold = DS; 068 069 if 0, % used for debugging 070 Di = createspec('dir','f'); 071 Di.f = fi; 072 Di.S = DS; 073 Di.theta = thetai; 074 wspecplot(Di,2) 075 pause 076 end 077 end 078 close(h) 079 return; % imlm 080
Comments or corrections to the WAFO group