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:
 createspec Spectrum structure constructor normspfn normalizes the spreading function wspecplot Plot a spectral density ctranspose ' Complex conjugate transpose. pause Wait for user response. pinv Pseudoinverse. squeeze Remove singleton dimensions.
This function is called by:
 dat2dspec Estimates the directional wave spectrum from timeseries emem Extended Maximum Entropy Method imlm Iterated maximum likelihood method for estimating the directional distribution mem maximum entropy method for estimating the directional distribution

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