Home > wafo > multidim > private > mem.m

mem

PURPOSE ^

maximum entropy method for estimating the directional distribution

SYNOPSIS ^

DS = mem(Sxyn,Gwt,thetai,fi,k)

DESCRIPTION ^

 MEM  maximum entropy method for estimating the directional distribution 
  
  CALL:  DS = mem(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 = mem(Sxyn,Gwt,thetai,fi,k) 
002 %MEM  maximum entropy method for estimating the directional distribution 
003 % 
004 % CALL:  DS = mem(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 %Tested on matlab 6.0   
019 %History :  
020 %revised pab Oct-Nov 2002   
021 %By pab 1999? 
022  
023 [m, nt, nf] = size(Gwt); 
024  
025  
026 H   = zeros(m*m,nt,nf); 
027 phi = zeros(m*m,nf); 
028  
029 % Eliminate meaningless equations such as those determined  
030 % from the zero co-spectrum and zero quadrature-spectrum. 
031 M = 0;  % M is the number of independent equations 
032  
033 dtheta=thetai(2)-thetai(1); 
034 tol = sqrt(eps); % treshold defining zero for transfer functions 
035 for ix=1:m 
036   for iy=ix:m 
037     Htemp  = Gwt(ix,:,:).*conj(Gwt(iy,:,:));                         
038     if (any(any(abs(diff(real(Htemp),1))>tol*dtheta))) 
039       M        = M+1; 
040       phi(M,:) = real(Sxyn(ix,iy,:)); 
041       H(M,:,:) = real(Htemp);  
042     end 
043     if any(any(abs(diff(imag(Htemp),1))>tol*dtheta) ) 
044       M            = M+1; 
045       phi(M,:) = imag(Sxyn(ix,iy,:)); 
046       H(M,:,:) = imag(Htemp);           
047     end 
048   end 
049 end 
050 M 
051 H   = H(1:M,:,:); 
052 phi = phi(1:M,:); 
053  
054 warningState = warning; 
055 warning off; 
056  
057 % Constants controlling the calculation 
058 coefAbsTol    = 0.1; 
059 coefAbsTol2   = 100; 
060  
061 errorTol      = 0.1;%sqrt(eps); 
062 maxIter       = 25; 
063 Li            = 1 ; %relaxation parameter 
064 maxCoef       = 5000; 
065  
066 display =1; 
067  
068 La(1:M)  = zeros(1,M); % Lagrange multiplier 
069 dLa(1:M) = zeros(1,M); % [deltaA,deltaB] 
070      
071  
072 DS = repmat(1/(2*pi),nt,nf); % initialize DS 
073       
074 h = waitbar(0,'Please wait...MEM calculation'); 
075  
076 % compute a fast estimate in order to find a good 
077 % starting guess for the Lagrange multipliers  
078 DS0 = log(mlm(Sxyn,Gwt,thetai,fi,k)); 
079 Oneth = ones(nt,1); 
080  
081 for ff=k, % loop over frequencies where S(f)>0 
082   waitbar(ff/k(end),h) 
083      
084   Hj   = H(:,:,ff).';  %H(1:M,1:nt)  
085   Phij = repmat(phi(1:M,ff).',nt,1);  
086    
087   stop = 0; 
088   count    = 0; 
089   lambda   = Li; %relaxation parameter   
090   La0 = -[Oneth Hj] \ DS0(:,ff); % starting guess for the Lagrange multipliers  
091   La(1:M) = La0(2:M+1); 
092   %La(1:M)  = zeros(1,M); % Lagrange multiplier 
093    
094   dLa(1:M) = zeros(1,M); % [deltaA,deltaB] 
095      
096   %coef(2*N-1:2*N)      = zeros(1,2); % [a,b] 
097   %deltaCoef(2*N-1:2*N) = zeros(1,2); % [deltaA,deltaB] 
098    
099   %size(La),size(Hj)   
100   exponent  = -(Hj*(La.')); 
101   L0     = -max(exponent); % trick in order to avoid infinities 
102   Fn     = exp(L0+exponent); 
103   Dn     = repmat(Fn/(trapz(Fn)*dtheta),1,M); %Fn(theta|f)/norm(Fn)=Dn(theta|f) 
104   PhiHD  = (Phij-Hj).*Dn; 
105   B      = trapz(PhiHD,1)*dtheta; % 1xM 
106      
107   error1 = max(abs(B)); 
108      
109   while(~stop), %%Use Newton-Raphson iteration to find model. 
110     count = count+1; 
111             
112     %A is the jacobian in the Newton iteration. 
113     for ix=1:M 
114       A(ix,:)  = trapz(PhiHD.*repmat(Hj(:,ix),1,M),1)*dtheta; 
115     end 
116        
117     %A, disp('Hit a key'),pause 
118     dLaOld(1:M) = La; 
119     %dLa(1:M)    = B/(A.'); % solve eq. 27  
120     dLa(1:M)   = B*pinv(A.'); % solve eq. 27  
121     La(1:M) = La + lambda*dLa; 
122        
123     if (maxCoef<inf)  
124       %This option is not described in Hashimoto,N. (1997). 
125       %Make sure the coefficients do not diverge to infinity 
126       k0 = find(abs(La)>maxCoef); 
127       if any(k0) 
128     dLa(k0)=(sign(dLa(k0)).*maxCoef-... 
129          (La(k0)-lambda*dLa(k0)))/lambda; 
130     La(k0)     = sign(La(k0)).*maxCoef; 
131       end 
132     end 
133     exponent  = -(Hj*(La.')); 
134     L0     = -max(exponent); % trick in order to avoid infinities 
135     Fn     = exp(L0+exponent); 
136     Dn     = repmat(Fn/(trapz(Fn)*dtheta),1,M); %Fn(theta|f)/norm(Fn)=Dn(theta|f) 
137     PhiHD  = (Phij-Hj).*Dn; 
138     B      = trapz(PhiHD,1)*dtheta; % 1xM 
139      
140      
141     error2 = error1; 
142     error1 = max(abs(B)); 
143        
144     if (~any(abs(dLa(1:M))>coefAbsTol)  | (error1 <= errorTol)); 
145       stop = 1; 
146     elseif (count>maxIter | ((error2<error1) & ... 
147                  (any(abs(dLa-dLaOld)>coefAbsTol2) ))), 
148     dLa(1:M) = 0; 
149     La(1:M)  = La0(2:M+1); 
150      
151     % Coefficients are diverging and error is increasing 
152     % solution: under relax the computation or quit 
153     if(lambda>Li*2^-4),  
154       lambda   = lambda*0.5; 
155       count    = 0; 
156        
157        
158       Dn     = repmat(1/(2*pi),nt,M); %Dn(theta|f) 
159       PhiHD  = (Phij-Hj).*Dn; 
160       B      = trapz(PhiHD,1)*dtheta; % 1xM 
161       error2 = inf; 
162       error1 = max(abs(B)); 
163     else 
164       stop = 1; 
165     end 
166       end 
167        
168       if 0 %stop, 
169       N = 4; 
170       subplot(N,1,1), semilogy(abs(B)+eps,'*'),hline(errorTol), title('error') 
171       subplot(N,1,2),plot(dLa,'g*'), hline(coefAbsTol*[-1,1]),title('deltaCoef') 
172       subplot(N,1,3), plot(La,'r*'), title('coef') 
173       subplot(N,1,4), plot(thetai,Dn), title('D(theta|f)') 
174       drawnow,% 
175       disp('Hit any key'),pause 
176       end 
177     end %while Newton-Raphson 
178              
179    
180   if display>0 
181     disp(sprintf('f = %g \t \t Error = %g',fi(ff),error1)) 
182   end 
183    
184   exponent  = -(Hj*(La.')); 
185   L0        = -max(exponent); % trick in order to avoid infinities 
186   DS(:,ff)  = exp(L0+exponent); 
187          
188 end % for ff= 
189 close(h) 
190 DS = normspfn(DS,thetai); 
191  
192 warning(warningState); 
193  
194 if (all(abs(DS(:)-DS(1))<sqrt(eps))) 
195   warning('No main direction found. Check the estimated spectrum!') 
196 end 
197  
198 return; % MEM 
199  
200  
201 return 
202  
203 % old call kept just in case 
204 [m, nt, nf] = size(Gwt); 
205 % This method assumes  a directional spreading of the form 
206 %   D(theta|w) = C*exp( qj*La.')  
207 % where C is a  normalization constant, La is vector of the unknown Lagrange 
208 % multipliers and qj is a cosine/sinus function. 
209  
210 % compute a fast estimate in order to find a good 
211 % starting guess for the Lagrange multipliers  
212 DS0 = log(mlm(Sxy,Gwt,thetai,fi,k)); 
213  
214 DS = ones(nt,nf)/(2*pi); % If S(f)==0 then set D(theta,f)=1/(2*pi); 
215  
216 ii = [1:m]'*ones(1,m); 
217 jj = ones(m,1)*[1:m]; 
218  
219 [indx, indy] = find(ii < jj);   % indices to upper triangular array 
220 ind = sub2ind([m m],indx,indy); % convert to linear index 
221  
222 Sxy = permute(Sxy,[3,1,2]); 
223 Gwt = permute(Gwt,[2 1 3 ]); 
224  
225  
226  
227 % Set up the non linear equation to solve or to minimize:  
228 % int (Pj(ones(nt,1),:)-qj).*repmat(exp(qj*La.'),1,M2) dtheta^2 =0; 
229 % simpson(thetai, (Pj(ones(nt,1),:)-qj).*repmat(exp(qj*La.'),1,M2)).^2=0; 
230 eqstr=inline('((simpson(P1,P2.*repmat(exp(P3*transpose(x)),1,P4))))',4); 
231  
232 M2 = m*(m-1); M1 = M2/2; 
233  
234 qj = zeros(nt,M2); 
235 La = rand(1,M2)*10-5; %zeros(1,M2); % starting guess 
236 Pj = La; 
237 if 0 
238   tmp = thetai(:) * (1:M1); 
239   qj2 = [cos(tmp) sin(tmp)]; 
240   size(qj2) 
241 end 
242  
243  
244 % find matlab version 
245 vstr = version; vind = find(vstr=='.'); 
246 vstr = str2num(vstr(1:vind(2)-1)); 
247  
248 Oneth = ones(nt,1); 
249  
250  
251 for ix = k, % looping over non-zero values of S(f) only 
252   Pj(1:M1)       = real(Sxy(ix,ind)); 
253   Pj(M1+1:M2)    = imag(Sxy(ix,ind)); 
254   tmp            = Gwt(:,indx,ix).*conj(Gwt(:,indy,ix)); 
255   qj(:,1:M1)     = real(tmp); 
256   qj(:,M1+1:M2)  = imag(tmp); 
257   %qj2 = qj; 
258   %plot(qj) 
259   %pause 
260   % Find a good starting guess for La based on log(DS_mlm). 
261   if 1, 
262      La0 = [Oneth qj] \ DS0(:,ix); 
263   else 
264     La0 = [Oneth qj(:,[1:2 M1+1:M1+2] )] \ DS0(:,ix); 
265     La  = [La0(2:3).' zeros(1,M1-2) La0(4:5).' zeros(1,M1-2) ] 
266   end 
267   %mans = memfun(La,thetai,Pj(ones(nt,1),:)-qj,qj,M2); 
268   %if any(isnan(mans)),  La = rand(1,M2)*10-5;  end 
269    
270   % find the Lagrange multipliers La  
271   if vstr>5.2, % matlab version 5.3 or higher 
272     %La = fsolve('memfun',La,optimset('fsolve'),thetai,Pj(ones(nt,1),:)-qj,qj,M2); 
273     La = fminsearch('memfun',La,[],thetai,Pj(ones(nt,1),:)-qj,qj(:,1:M2),M2); 
274     %La = fminunc('memfun',La,[],thetai,Pj(ones(nt,1),:)-qj,qj,M2) 
275   else  % matlab version 5.2 or higher 
276     %La = fsolve('memfun',La,[],[],thetai,Pj(ones(nt,1),:)-qj,qj,M2); 
277     La = fmins('memfun',La,[],[],thetai,Pj(ones(nt,1),:)-qj,qj(:,1:M2),M2); 
278     %La = fminu('memfun',La,[],[],thetai,Pj(ones(nt,1),:)-qj,qj,M2) 
279     %La = fmins(eqstr,La,[],[],thetai,Pj(ones(nt,1),:)-qj,qj,M2) 
280   end 
281    
282   La(isnan(La)) = 0; % make sure that isnans is zero 
283   %mans2 = memfun(La,thetai,Pj(ones(nt,1),:)-qj,qj,M2); 
284    
285   
286   disp([' Finished ' num2str(ix) ' of ' num2str(nf)]) 
287   %disp(num2str([mans;mans2])) 
288   %disp(num2str(La)) 
289   DS(:,ix) = exp(qj*La.'); 
290   %if sum(abs(mans2)) > 1  La = rand(1,M2)*10-5;  end  
291 end 
292 %Normalize so that int D(theta,f) dtheta = 1 for each f  
293 %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
294 DS = normspfn(DS,thetai); 
295  
296 return; 
297

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