Home > wafo > multidim > private > emem.m

emem

PURPOSE ^

Extended Maximum Entropy Method

SYNOPSIS ^

DS = emem(Sxyn,Gwt,theta,fi,k,opt)

DESCRIPTION ^

 EMEM  Extended Maximum Entropy Method  
  
   CALL: DS = emem(Sxyn,Gwt,thetai,fi,k)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function DS = emem(Sxyn,Gwt,theta,fi,k,opt) 
002 %EMEM  Extended Maximum Entropy Method  
003 % 
004 %  CALL: DS = emem(Sxyn,Gwt,thetai,fi,k) 
005 % 
006  
007 %Tested on matlab 6.0   
008 %History :  
009 %By pab Oct-Nov 2002 
010    
011 %Reference: 
012 % Hashimoto,N. (1997)  
013 % "Analysis of the directional wave spectra from field data.", Advances in  
014 % Coastal and Ocean Engineering, Vol.3., pp.103-143 
015  
016 [m, nt, nf] = size(Gwt); 
017  
018  
019 H   = zeros(m*m,nt,nf); 
020 phi = zeros(m*m,nf); 
021  
022 % Eliminate meaningless equations such as those determined  
023 % from the zero co-spectrum and zero quadrature-spectrum. 
024 M = 0;  % M is the number of independent equations 
025  
026 dtheta = theta(2)-theta(1); 
027 tol    = sqrt(eps); % treshold defining zero for transfer functions 
028 for ix=1:m 
029   for iy=ix:m 
030     Htemp  = Gwt(ix,:,:).*conj(Gwt(iy,:,:));                         
031     if (any(any(abs(diff(real(Htemp),1))>tol*dtheta))) 
032       M        = M+1; 
033       phi(M,:) = real(Sxyn(ix,iy,:)); 
034       H(M,:,:) = real(Htemp);  
035     end 
036     if any(any(abs(diff(imag(Htemp),1))>tol*dtheta) ) 
037       M            = M+1; 
038       phi(M,:) = imag(Sxyn(ix,iy,:)); 
039       H(M,:,:) = imag(Htemp);           
040     end 
041   end 
042 end 
043  
044 % Note: H and Phi here are normalized in another way than described in Hashimoto,N. (1997). 
045 H   = H(1:M,:,:); 
046 phi = phi(1:M,:); 
047  
048 warningState = warning; 
049 warning off; 
050  
051 M2 = floor(M/2+1); 
052  
053 % Constants controlling the calculation 
054 coefAbsTol    = 0.01; 
055 coefAbsTol2   = 500; 
056 errorTol      = 0.0005;%sqrt(eps); 
057 maxIter       = 25; 
058 minModelOrder = 1; 
059 maxModelOrder = M2; 
060 Li            = 1 ; %relaxation parameter 
061 AICTol        = 0.1; 
062 maxCoef       = 10000; 
063  
064 display =0; 
065 if nargin>5 
066   if ~isempty(opt.errortol), errorTol = opt.errortol; end 
067   if ~isempty(opt.maxiter),  maxIter  = opt.maxiter; end 
068   if ~isempty(opt.relax),    Li       = opt.relax; end 
069   if ~isempty(opt.maxcoef),  maxCoef  = opt.maxcoef; end 
070   if ~isempty(opt.message),  display  = opt.message; end 
071   if ~isempty(opt.coefabstol),       coefAbsTol = opt.coefabstol; end 
072   if ~isempty(opt.minmodelorder), minModelOrder = opt.minmodelorder; end 
073   if ~isempty(opt.maxmodelorder), maxModelOrder = opt.maxmodelorder; end 
074 end 
075  
076  
077  
078 cosn  = cos([1:maxModelOrder]'*theta'); 
079 sinn  = sin([1:maxModelOrder]'*theta'); 
080 cosnt = permute(repmat(cosn,[1 1 M]),[2 3 1]); 
081 sinnt = permute(repmat(sinn,[1 1 M]),[2 3 1]); 
082  
083 AIC       = zeros(1,maxModelOrder); 
084 XY        = zeros(2*maxModelOrder,M); 
085 coef      = zeros(1,2*maxModelOrder); % [a,b] 
086 deltaCoef = zeros(1,2*maxModelOrder); % [deltaA,deltaB] 
087  
088 Nbest = 0; %Best model order 
089  
090 DS = repmat(1/(2*pi),nt,nf); % initialize DS 
091       
092 h = waitbar(0,'Please wait...EMEM calculation'); 
093  
094 % compute a fast estimate in order to find a good 
095 % starting guess for the Lagrange multipliers  
096 DS0   = log(mlm(Sxyn,Gwt,theta,fi,k)); 
097 Oneth = ones(nt,1); 
098  
099 for ff=k, % loop over frequencies where S(f)>0 
100   waitbar(ff/k(end),h) 
101      
102   Hj     = H(:,:,ff).';  %H(1:M,1:nt)  
103   Phij = repmat(phi(1:M,ff).',nt,1);  
104      
105   stop = 0; 
106    
107   % Assume model order change smoothly over the frequencies 
108   localMinModelOrder = max(minModelOrder,ceil(Nbest/2)); 
109    
110   N    = localMinModelOrder-1; 
111       
112   while(~stop), % loop until the right model order is found 
113     N         = N+1; %increment modelOrder 
114      
115     % coef0=starting guess for coefs multipliers  
116     coef0 = zeros(1,2*N+1); 
117     %coef0 = [Oneth, cosn(1:N,:).', sinn(1:N,:).'] \ DS0(:,ff);  
118     coef(1:2*N) = coef0(2:2*N+1); 
119     %coef(1:2*N)      = zeros(1,2*N); % [a,b] 
120     deltaCoef(1:2*N) = zeros(1,2*N); % [deltaA,deltaB] 
121      
122     %coef(2*N-1:2*N)      = zeros(1,2); % [a,b] 
123     %deltaCoef(2*N-1:2*N) = zeros(1,2); % [deltaA,deltaB] 
124      
125     count      = 0; 
126     lambda     = Li; %relaxation parameter      
127     modelFound = 0; 
128      
129     exponent  = (coef(1:N)*cosn(1:N,:)+coef(N+1:2*N)*sinn(1:N,:))'; 
130     a0     = -max(exponent); % trick in order to avoid infinities 
131     Fn     = exp(a0+exponent); 
132     Dn     = repmat(Fn/(trapz(Fn)*dtheta),1,M); %Fn(theta|f)/norm(Fn)=Dn(theta|f) 
133     PhiHD  = (Phij-Hj).*Dn; 
134     Z      = trapz(PhiHD,1)*dtheta; % 1xM 
135     error1 = max(abs(Z)); 
136      
137     while ( ~stop & ~modelFound),%Use Newton-Raphson iteration to find model. 
138        
139       count = count+1; 
140        
141       %XY is the negative jacobian in the Newton iteration. 
142       for ix=1:N 
143     %This better than eq. 97 and 98 in Hashimoto (1997). 
144     %X(ix,1:M)= 
145     XY(ix,:)   = (Z.*trapz(Dn.*cosnt(:,:,ix),1) - ... 
146               trapz(PhiHD.*cosnt(:,:,ix),1))*dtheta; 
147     %Y(ix,1:M)= 
148     XY(N+ix,:) = (Z.*trapz(Dn.*sinnt(:,:,ix),1) - ... 
149               trapz(PhiHD.*sinnt(:,:,ix),1))*dtheta ; 
150       end 
151        
152       deltaCoefOld     = deltaCoef(1:2*N); 
153       deltaCoef(1:2*N) = Z/XY(1:2*N,:); % solve eq. 95 by least square 
154                      
155       coef(1:2*N) = coef(1:2*N) + lambda*deltaCoef(1:2*N); 
156        
157       if (maxCoef<inf) %This option is not described in Hashimoto,N. (1997). 
158     %Make sure the coefficients do not diverge to infinity 
159     k0 = find(abs(coef(1:2*N))>maxCoef); 
160     if any(k0) 
161       deltaCoef(k0)=(sign(deltaCoef(k0)).*maxCoef-... 
162              (coef(k0)-lambda*deltaCoef(k0)))/lambda; 
163       coef(k0)     = sign(coef(k0)).*maxCoef; 
164     end 
165       end 
166        
167       exponent  = (coef(1:N)*cosn(1:N,:)+coef(N+1:2*N)*sinn(1:N,:))'; 
168       a0     = -max(exponent); % trick in order to avoid infinities 
169       Fn     = exp(a0+exponent); 
170       Dn     = repmat(Fn/(trapz(Fn)*dtheta),1,M); %Fn(theta|f)/norm(Fn)=Dn(theta|f) 
171       PhiHD  = (Phij-Hj).*Dn; 
172       Z      = trapz(PhiHD,1)*dtheta; % 1xM 
173        
174       error2 = error1; 
175       error1 = max(abs(Z)); 
176        
177       if (~any(abs(deltaCoef(1:2*N))>coefAbsTol)  | (error1 <= errorTol)); 
178      modelFound = 1; 
179       elseif (count>maxIter | ((error2<error1) & ... 
180     (any(abs(deltaCoef(1:2*N)-deltaCoefOld)>coefAbsTol2) ))), 
181      
182     % Coefficients are diverging and error is increasing 
183     % solution: under relax the computation or quit 
184     if(lambda>Li*2^-4),  
185       lambda      = lambda*0.5; 
186       count       = 0; 
187       deltaCoef(1:2*N) = 0; 
188      % coef(1:2*N)      = 0; 
189       coef(1:2*N) = coef0(2:2*N+1); 
190       Dn     = repmat(1/(2*pi),nt,M); %Dn(theta|f) 
191       PhiHD  = (Phij-Hj).*Dn; 
192       Z      = trapz(PhiHD,1)*dtheta; % 1xM 
193       error2 = inf; 
194       error1 = max(abs(Z)); 
195     else 
196       stop = 1; 
197     end 
198       end 
199        
200       if 0 %modelFound 
201       Np = 4; 
202       subplot(Np,1,1), semilogy(abs(Z)+eps,'*'),hline(errorTol), title('error') 
203       subplot(Np,1,2),plot(deltaCoef(1:2*N),'g*'), title('deltaCoef') 
204       subplot(Np,1,3), plot(coef(1:2*N),'r*'), title('coef') 
205       subplot(Np,1,4), plot(theta,Dn) 
206       drawnow,% 
207       disp('Hit any key'),pause 
208       end 
209     end %while Newton-Raphson 
210              
211     AIC(N) = M*(log(2*pi*var(Z))+1)+4*N+2; %Aikakes Information Criterion 
212        
213     if(N>localMinModelOrder), 
214       stop = stop | (AIC(N)+AICTol>AIC(N-1)) ; 
215     end 
216     if isnan(AIC(N)), stop = 1;  end %this should never happen 
217      
218     if (stop), 
219       if (N>localMinModelOrder) 
220     Nbest = N-1; 
221     coef(1:2*Nbest) = coefOld; 
222       else 
223     Nbest = 1; 
224     coef = zeros(1,2*Nbest); 
225       end 
226     else 
227       Nbest   = N; 
228       stop    = (N>=maxModelOrder); 
229       coefOld = coef(1:2*N); 
230     end 
231        
232   end % while 
233    
234   if display>0 
235     disp(sprintf('f = %g \t \t Model order = %d \t \t error = %g',fi(ff),Nbest,max(abs(Z)))) 
236   end 
237   %plot(deltaCoef,'r*'),hold on,plot(coef(1:2*Nbest),'g*'),hold off,pause 
238    
239   exponent  = (coef(1:Nbest)*cosn(1:Nbest,:)+coef(Nbest+1:2*Nbest)*sinn(1:Nbest,:)).'; 
240   a0     = -max(exponent); % trick in order to avoid infinities 
241   DS(:,ff) = exp(a0+exponent); 
242          
243 end % for ff= 
244 close(h) 
245 DS = normspfn(DS,theta); 
246  
247 warning(warningState); 
248  
249 if (all(abs(DS(:)-DS(1))<sqrt(eps))) 
250   warning('No main direction found. Check the estimated spectrum!') 
251 end 
252  
253 return; % EMEM

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