001 function DS = emem(Sxyn,Gwt,theta,fi,k,opt)
002
003
004
005
006
007
008
009
010
011
012
013
014
015
016 [m, nt, nf] = size(Gwt);
017
018
019 H = zeros(m*m,nt,nf);
020 phi = zeros(m*m,nf);
021
022
023
024 M = 0;
025
026 dtheta = theta(2)-theta(1);
027 tol = sqrt(eps);
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
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
054 coefAbsTol = 0.01;
055 coefAbsTol2 = 500;
056 errorTol = 0.0005;
057 maxIter = 25;
058 minModelOrder = 1;
059 maxModelOrder = M2;
060 Li = 1 ;
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);
086 deltaCoef = zeros(1,2*maxModelOrder);
087
088 Nbest = 0;
089
090 DS = repmat(1/(2*pi),nt,nf);
091
092 h = waitbar(0,'Please wait...EMEM calculation');
093
094
095
096 DS0 = log(mlm(Sxyn,Gwt,theta,fi,k));
097 Oneth = ones(nt,1);
098
099 for ff=k,
100 waitbar(ff/k(end),h)
101
102 Hj = H(:,:,ff).';
103 Phij = repmat(phi(1:M,ff).',nt,1);
104
105 stop = 0;
106
107
108 localMinModelOrder = max(minModelOrder,ceil(Nbest/2));
109
110 N = localMinModelOrder-1;
111
112 while(~stop),
113 N = N+1;
114
115
116 coef0 = zeros(1,2*N+1);
117
118 coef(1:2*N) = coef0(2:2*N+1);
119
120 deltaCoef(1:2*N) = zeros(1,2*N);
121
122
123
124
125 count = 0;
126 lambda = Li;
127 modelFound = 0;
128
129 exponent = (coef(1:N)*cosn(1:N,:)+coef(N+1:2*N)*sinn(1:N,:))';
130 a0 = -max(exponent);
131 Fn = exp(a0+exponent);
132 Dn = repmat(Fn/(trapz(Fn)*dtheta),1,M);
133 PhiHD = (Phij-Hj).*Dn;
134 Z = trapz(PhiHD,1)*dtheta;
135 error1 = max(abs(Z));
136
137 while ( ~stop & ~modelFound),
138
139 count = count+1;
140
141
142 for ix=1:N
143
144
145 XY(ix,:) = (Z.*trapz(Dn.*cosnt(:,:,ix),1) - ...
146 trapz(PhiHD.*cosnt(:,:,ix),1))*dtheta;
147
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,:);
154
155 coef(1:2*N) = coef(1:2*N) + lambda*deltaCoef(1:2*N);
156
157 if (maxCoef<inf)
158
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);
169 Fn = exp(a0+exponent);
170 Dn = repmat(Fn/(trapz(Fn)*dtheta),1,M);
171 PhiHD = (Phij-Hj).*Dn;
172 Z = trapz(PhiHD,1)*dtheta;
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
183
184 if(lambda>Li*2^-4),
185 lambda = lambda*0.5;
186 count = 0;
187 deltaCoef(1:2*N) = 0;
188
189 coef(1:2*N) = coef0(2:2*N+1);
190 Dn = repmat(1/(2*pi),nt,M);
191 PhiHD = (Phij-Hj).*Dn;
192 Z = trapz(PhiHD,1)*dtheta;
193 error2 = inf;
194 error1 = max(abs(Z));
195 else
196 stop = 1;
197 end
198 end
199
200 if 0
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
210
211 AIC(N) = M*(log(2*pi*var(Z))+1)+4*N+2;
212
213 if(N>localMinModelOrder),
214 stop = stop | (AIC(N)+AICTol>AIC(N-1)) ;
215 end
216 if isnan(AIC(N)), stop = 1; end
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
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
238
239 exponent = (coef(1:Nbest)*cosn(1:Nbest,:)+coef(Nbest+1:2*Nbest)*sinn(1:Nbest,:)).';
240 a0 = -max(exponent);
241 DS(:,ff) = exp(a0+exponent);
242
243 end
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;