001 function DS = mem(Sxyn,Gwt,thetai,fi,k)
002
003
004
005
006
007
008
009
010
011
012
013
014
015
016
017
018
019
020
021
022
023 [m, nt, nf] = size(Gwt);
024
025
026 H = zeros(m*m,nt,nf);
027 phi = zeros(m*m,nf);
028
029
030
031 M = 0;
032
033 dtheta=thetai(2)-thetai(1);
034 tol = sqrt(eps);
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
058 coefAbsTol = 0.1;
059 coefAbsTol2 = 100;
060
061 errorTol = 0.1;
062 maxIter = 25;
063 Li = 1 ;
064 maxCoef = 5000;
065
066 display =1;
067
068 La(1:M) = zeros(1,M);
069 dLa(1:M) = zeros(1,M);
070
071
072 DS = repmat(1/(2*pi),nt,nf);
073
074 h = waitbar(0,'Please wait...MEM calculation');
075
076
077
078 DS0 = log(mlm(Sxyn,Gwt,thetai,fi,k));
079 Oneth = ones(nt,1);
080
081 for ff=k,
082 waitbar(ff/k(end),h)
083
084 Hj = H(:,:,ff).';
085 Phij = repmat(phi(1:M,ff).',nt,1);
086
087 stop = 0;
088 count = 0;
089 lambda = Li;
090 La0 = -[Oneth Hj] \ DS0(:,ff);
091 La(1:M) = La0(2:M+1);
092
093
094 dLa(1:M) = zeros(1,M);
095
096
097
098
099
100 exponent = -(Hj*(La.'));
101 L0 = -max(exponent);
102 Fn = exp(L0+exponent);
103 Dn = repmat(Fn/(trapz(Fn)*dtheta),1,M);
104 PhiHD = (Phij-Hj).*Dn;
105 B = trapz(PhiHD,1)*dtheta;
106
107 error1 = max(abs(B));
108
109 while(~stop),
110 count = count+1;
111
112
113 for ix=1:M
114 A(ix,:) = trapz(PhiHD.*repmat(Hj(:,ix),1,M),1)*dtheta;
115 end
116
117
118 dLaOld(1:M) = La;
119
120 dLa(1:M) = B*pinv(A.');
121 La(1:M) = La + lambda*dLa;
122
123 if (maxCoef<inf)
124
125
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);
135 Fn = exp(L0+exponent);
136 Dn = repmat(Fn/(trapz(Fn)*dtheta),1,M);
137 PhiHD = (Phij-Hj).*Dn;
138 B = trapz(PhiHD,1)*dtheta;
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
152
153 if(lambda>Li*2^-4),
154 lambda = lambda*0.5;
155 count = 0;
156
157
158 Dn = repmat(1/(2*pi),nt,M);
159 PhiHD = (Phij-Hj).*Dn;
160 B = trapz(PhiHD,1)*dtheta;
161 error2 = inf;
162 error1 = max(abs(B));
163 else
164 stop = 1;
165 end
166 end
167
168 if 0
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
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);
186 DS(:,ff) = exp(L0+exponent);
187
188 end
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;
199
200
201 return
202
203
204 [m, nt, nf] = size(Gwt);
205
206
207
208
209
210
211
212 DS0 = log(mlm(Sxy,Gwt,thetai,fi,k));
213
214 DS = ones(nt,nf)/(2*pi);
215
216 ii = [1:m]'*ones(1,m);
217 jj = ones(m,1)*[1:m];
218
219 [indx, indy] = find(ii < jj);
220 ind = sub2ind([m m],indx,indy);
221
222 Sxy = permute(Sxy,[3,1,2]);
223 Gwt = permute(Gwt,[2 1 3 ]);
224
225
226
227
228
229
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;
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
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,
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
258
259
260
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
268
269
270
271 if vstr>5.2,
272
273 La = fminsearch('memfun',La,[],thetai,Pj(ones(nt,1),:)-qj,qj(:,1:M2),M2);
274
275 else
276
277 La = fmins('memfun',La,[],[],thetai,Pj(ones(nt,1),:)-qj,qj(:,1:M2),M2);
278
279
280 end
281
282 La(isnan(La)) = 0;
283
284
285
286 disp([' Finished ' num2str(ix) ' of ' num2str(nf)])
287
288
289 DS(:,ix) = exp(qj*La.');
290
291 end
292
293
294 DS = normspfn(DS,thetai);
295
296 return;
297