Home > wafo > multidim > private > imlm.m

# imlm

## PURPOSE

Iterated maximum likelihood method for estimating the directional distribution

## SYNOPSIS

DS = imlm(Sxy,Gwt,thetai,fi,k,opt)

## DESCRIPTION

``` 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)```

## CROSS-REFERENCE INFORMATION

This function calls:
 createspec Spectrum structure constructor mlm maximum likelihood method for estimating the directional distribution normspfn normalizes the spreading function simpson Numerical integration with the Simpson method wspecplot Plot a spectral density close Close figure. num2str Convert number to string. (Fast version) pause Wait for user response. squeeze Remove singleton dimensions. waitbar Display wait bar.
This function is called by:
 dat2dspec Estimates the directional wave spectrum from timeseries

## SOURCE CODE

```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```

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