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: This function is called by:

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