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)
Spectrum structure constructor | |
normalizes the spreading function | |
Plot a spectral density | |
' Complex conjugate transpose. | |
Wait for user response. | |
Pseudoinverse. | |
Remove singleton dimensions. |
Estimates the directional wave spectrum from timeseries | |
Extended Maximum Entropy Method | |
Iterated maximum likelihood method for estimating the directional distribution | |
maximum entropy method for estimating the directional distribution |
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
Comments or corrections to the WAFO group