Home > wafo > multidim > private > mlm.m

mlm

PURPOSE ^

maximum likelihood method for estimating the directional distribution

SYNOPSIS ^

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

DESCRIPTION ^

  MLM  maximum likelihood method for estimating the directional distribution 
  
  CALL  DS = mlm(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 = mlm(Sxy,Gwt,thetai,fi,k,opt) 
002 % MLM  maximum likelihood method for estimating the directional distribution 
003 % 
004 % CALL  DS = mlm(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 [m nt nf] = size(Gwt); 
019  
020   
021 % size(Sxy),nf,nt 
022 %----------------------------------------------------- 
023 %   inverting matrix of cross-spectra at every frequency 
024 %------------------------------------------------------ 
025  
026 I=eye(m)*sqrt(eps); 
027 if 1, % New call, slightly faster  
028   % initialize DS 
029   DS = repmat(1/(2*pi),nt,nf); % If S(f)==0 then set D(theta,f)=1/(2*pi); 
030   for ix = k(:)', % looping over non-zero values of S(f) only 
031     %Gmat = Gwt(1:m,1:nt,ix).'; % = transpose(Gwt(1:m,1:nt,ix)) 
032     %H    =  real((conj(Gmat)*pinv(Sxy(:,:,ix)+I)).*Gmat); 
033       
034     H   = real((ctranspose(Gwt(1:m,1:nt,ix))*pinv(Sxy(:,:,ix)+I)).*(Gwt(1:m,1:nt,ix)).'); 
035     tmp = sum(H,2); 
036     if any(tmp==0) 
037       if all(tmp==0) 
038     DS(:,ix)= 1/(2*pi*nt); 
039       else 
040     tmp(tmp==0)=eps; 
041     DS(:,ix) = 1./tmp; 
042       end 
043     else 
044       DS(:,ix) = 1./tmp; %  size(H) = nt x m 
045     end 
046   end 
047 else % Old call: 
048      
049   DS  = zeros(nt,nf); %initialize 
050   for ix=k, % 1:nf 
051     Sxy(:,:,ix) = pinv(Sxy(:,:,ix)); %Gmn^-1 
052   end 
053   for ix=1:m, %    m-1, 
054     Sm1 = real(squeeze(Sxy(ix,ix,:))).'; 
055     Gm1 = conj(squeeze(Gwt(ix,:,:))); 
056     DS = DS+Sm1(ones(nt,1),:).*abs(Gm1).^2 ; 
057     for iy = (ix+1):m, %    m-1, 
058       Sm1 = squeeze(Sxy(ix,iy,:)).'; 
059       Gm2 = (squeeze(Gwt(iy,:,:))); 
060       DS  = DS+2*real(Sm1(ones(nt,1),:).*Gm1.*Gm2 ); 
061     end 
062   end 
063   DS  = 1./real(DS); 
064 end 
065  
066  
067 if 0, 
068   Di   = createspec('dir','f'); 
069   Di.f = fi; 
070   Di.S = DS; 
071   Di.theta = thetai; 
072   wspecplot(Di,2) 
073   pause 
074 end 
075  
076 %Normalize so that int D(theta,f) dtheta = 1 for each f  
077 %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
078 DS = normspfn(DS,thetai); 
079  
080 if 0, 
081   Di = createspec('dir','f'); 
082   Di.f = fi; 
083   Di.S = DS; 
084   Di.theta=thetai; 
085   wspecplot(Di,2) 
086   pause 
087 end 
088 return 
089

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