1 |
function [S]=IMLM(xps,trm,kx,Ss,pidirs,miter,displ) |
---|
2 |
|
---|
3 |
gamma=0.1; |
---|
4 |
beta=1.0; |
---|
5 |
alpha=0.1; |
---|
6 |
|
---|
7 |
szd=size(xps,1); |
---|
8 |
ffreqs=size(xps,3); |
---|
9 |
ddirs=size(trm,3); |
---|
10 |
|
---|
11 |
ddir=8*atan(1.0)/ddirs; |
---|
12 |
|
---|
13 |
if(displ<2) |
---|
14 |
warning off; |
---|
15 |
end |
---|
16 |
|
---|
17 |
|
---|
18 |
for ff=1:ffreqs |
---|
19 |
if(displ>=1) |
---|
20 |
disp(['calculating for frequency' blanks(1) num2str(ff) ' of' blanks(1) num2str(ffreqs)]); |
---|
21 |
end |
---|
22 |
|
---|
23 |
for m=1:szd |
---|
24 |
for n=1:szd |
---|
25 |
H(1:ddirs,m,n)=trm(n,ff,1:ddirs); |
---|
26 |
Hs(1:ddirs,m,n)=conj(trm(m,ff,1:ddirs)); |
---|
27 |
expx(1:ddirs,m,n)=exp(i*kx(m,n,ff,1:ddirs)); |
---|
28 |
iexpx(1:ddirs,m,n)=exp(-i*kx(m,n,ff,1:ddirs)); |
---|
29 |
Htemp(:,m,n)=H(:,m,n).*Hs(:,m,n).*expx(:,m,n); |
---|
30 |
iHtemp(:,m,n)=H(:,m,n).*Hs(:,m,n).*iexpx(:,m,n); |
---|
31 |
end |
---|
32 |
end |
---|
33 |
|
---|
34 |
|
---|
35 |
invcps=inv(xps(:,:,ff)); |
---|
36 |
Sftmp=zeros(ddirs,1); |
---|
37 |
for m=1:szd |
---|
38 |
for n=1:szd |
---|
39 |
xtemp=invcps(m,n)*Htemp(:,m,n); |
---|
40 |
Sftmp(:)=Sftmp(:)+xtemp; |
---|
41 |
end |
---|
42 |
end |
---|
43 |
Eo=(1./Sftmp(:)); |
---|
44 |
kappa=1./(ddir*sum(Eo)); |
---|
45 |
Eo=kappa*Eo; |
---|
46 |
E=Eo; |
---|
47 |
T=Eo; |
---|
48 |
|
---|
49 |
for it=1:miter |
---|
50 |
for m=1:szd |
---|
51 |
for n=1:szd |
---|
52 |
expG(m,n,:)=iHtemp(:,m,n).*E(:); |
---|
53 |
ixps(m,n)=sum(expG(m,n,:))*ddir; |
---|
54 |
end |
---|
55 |
end |
---|
56 |
invcps=inv(ixps); |
---|
57 |
Sftmp=zeros(ddirs,1); |
---|
58 |
for m=1:szd |
---|
59 |
for n=1:szd |
---|
60 |
xtemp=invcps(m,n)*Htemp(:,m,n); |
---|
61 |
Sftmp(:)=Sftmp(:)+xtemp; |
---|
62 |
end |
---|
63 |
end |
---|
64 |
Told=T; |
---|
65 |
T=(1./Sftmp(:)); |
---|
66 |
kappa=1./(ddir*sum(T)); |
---|
67 |
T=kappa*T; |
---|
68 |
|
---|
69 |
%lambda=ones(size(T))-(T./Eo) |
---|
70 |
%ei=gamma*lambda.*E; |
---|
71 |
ei=gamma*((Eo-T)+alpha*(T-Told)); |
---|
72 |
E=E+ei; |
---|
73 |
kappa=1./(ddir*sum(E)); |
---|
74 |
E=kappa*E; |
---|
75 |
|
---|
76 |
|
---|
77 |
end |
---|
78 |
|
---|
79 |
S(ff,:)=Ss(1,ff)*E'; |
---|
80 |
end |
---|
81 |
|
---|
82 |
warning on; |
---|
83 |
|
---|
84 |
|
---|
85 |
|
---|