1 |
function [S]=EMEP(xps,trm,kx,Ss,pidirs,miter,displ) |
---|
2 |
|
---|
3 |
szd=size(xps,1); |
---|
4 |
freqs=size(xps,3); |
---|
5 |
ddirs=size(trm,3); |
---|
6 |
|
---|
7 |
ddir=abs(pidirs(2)-pidirs(1)); |
---|
8 |
pi=4.0*atan(1.0); |
---|
9 |
|
---|
10 |
if(displ<2) |
---|
11 |
warning off; |
---|
12 |
end |
---|
13 |
|
---|
14 |
Co=real(xps); |
---|
15 |
Quad=-imag(xps); |
---|
16 |
|
---|
17 |
for ff=1:freqs |
---|
18 |
|
---|
19 |
xpsx(:,:,ff)=diag(xps(:,:,ff))*(diag(xps(:,:,ff))'); |
---|
20 |
sigCo(:,:,ff)=sqrt(0.5*(xpsx(:,:,ff)+Co(:,:,ff).^2-Quad(:,:,ff).^2)); |
---|
21 |
sigQuad(:,:,ff)=sqrt(0.5*(xpsx(:,:,ff)-Co(:,:,ff).^2+Quad(:,:,ff).^2)); |
---|
22 |
end |
---|
23 |
|
---|
24 |
|
---|
25 |
for ff=1:freqs |
---|
26 |
|
---|
27 |
index=0; |
---|
28 |
for m=1:szd |
---|
29 |
for n=m:szd |
---|
30 |
expx(1:ddirs)=exp(-i*kx(m,n,ff,1:ddirs)); |
---|
31 |
Hh(1:ddirs)=trm(m,ff,1:ddirs); |
---|
32 |
Hhs(1:ddirs)=conj(trm(n,ff,1:ddirs)); |
---|
33 |
Htemp=(Hh.*Hhs.*expx); |
---|
34 |
|
---|
35 |
|
---|
36 |
if(Htemp(1)~=Htemp(2)) |
---|
37 |
index=index+1; |
---|
38 |
phi(index,ff)=real(xps(m,n,ff))./(sigCo(m,n,ff)*Ss(1,ff)); |
---|
39 |
H(1:ddirs,index,ff)=real(Htemp)./sigCo(m,n,ff); |
---|
40 |
|
---|
41 |
if(kx(m,n,1,1)+kx(m,n,1,2)~=0) |
---|
42 |
index=index+1; |
---|
43 |
phi(index,ff)=imag(xps(m,n,ff))./(sigQuad(m,n,ff)*Ss(1,ff)); |
---|
44 |
H(1:ddirs,index,ff)=imag(Htemp)./sigQuad(m,n,ff); |
---|
45 |
end |
---|
46 |
end |
---|
47 |
end |
---|
48 |
end |
---|
49 |
end |
---|
50 |
M=index; |
---|
51 |
|
---|
52 |
|
---|
53 |
for eni=1:M/2+1 |
---|
54 |
cosnt(1:ddirs,1:M,eni)=cos(eni*pidirs')*ones(1,M); |
---|
55 |
sinnt(1:ddirs,1:M,eni)=sin(eni*pidirs')*ones(1,M); |
---|
56 |
end |
---|
57 |
|
---|
58 |
cosn=cos([1:M/2+1]'*pidirs); |
---|
59 |
sinn=sin([1:M/2+1]'*pidirs); |
---|
60 |
|
---|
61 |
|
---|
62 |
for ff=1:freqs |
---|
63 |
if (displ>=1) |
---|
64 |
disp(['calculating for frequency' blanks(1) num2str(ff) ' of' blanks(1) num2str(freqs)]); |
---|
65 |
end |
---|
66 |
Hi(1:ddirs,1:M)=H(1:ddirs,1:M,ff); |
---|
67 |
Phione=(ones(size(pidirs'))*phi(1:M,ff)'); |
---|
68 |
|
---|
69 |
keepgoing=1; |
---|
70 |
n=0; |
---|
71 |
AIC=[]; |
---|
72 |
|
---|
73 |
|
---|
74 |
while(keepgoing==1) |
---|
75 |
n=n+1; |
---|
76 |
|
---|
77 |
if(n<=M/2+1) |
---|
78 |
if(displ>0) |
---|
79 |
disp(strcat('model :',num2str(n))); |
---|
80 |
end |
---|
81 |
a1(1:n)=0.0; |
---|
82 |
b1(1:n)=0.0; |
---|
83 |
|
---|
84 |
a2(1:n)=100.0; |
---|
85 |
b2(1:n)=100.0; |
---|
86 |
|
---|
87 |
count=0; |
---|
88 |
rlx=1.0; |
---|
89 |
while(max(abs(a2(1:n)))>0.01 | max(abs(b2(1:n)))>0.01) |
---|
90 |
count=count+1; |
---|
91 |
Fn=(a1(1:n)*cosn(1:n,:)+b1(1:n)*sinn(1:n,:))'; |
---|
92 |
|
---|
93 |
Fnexp=exp(Fn)*ones(1,M); |
---|
94 |
PhiHF=(Phione-Hi).*Fnexp; |
---|
95 |
Z(1:M)=sum(PhiHF)./sum(Fnexp); |
---|
96 |
for eni=1:n |
---|
97 |
X(eni,1:M)=Z.*(... |
---|
98 |
( sum(Fnexp.*cosnt(:,:,eni))./sum(Fnexp) ) -... |
---|
99 |
( sum(PhiHF.*cosnt(:,:,eni))./sum(PhiHF) )... |
---|
100 |
); |
---|
101 |
Y(eni,1:M)=Z.*(... |
---|
102 |
( sum(Fnexp.*sinnt(:,:,eni))./sum(Fnexp) )-... |
---|
103 |
( sum(PhiHF.*sinnt(:,:,eni))./sum(PhiHF) )... |
---|
104 |
); |
---|
105 |
end |
---|
106 |
C(:,1:n)=(X(1:n,1:M))'; |
---|
107 |
C(:,n+1:2*n)=(Y(1:n,1:M))'; |
---|
108 |
|
---|
109 |
out=C(:,1:n*2)\Z'; |
---|
110 |
out=out'; |
---|
111 |
|
---|
112 |
a2old=a2(1:n); |
---|
113 |
b2old=b2(1:n); |
---|
114 |
a2=out(1:n); |
---|
115 |
b2=out(n+1:2*n); |
---|
116 |
if sum((abs(a2)-abs(a2old))>100) | sum((abs(b2)-abs(b2old))>100 |count>miter) |
---|
117 |
if(rlx>0.0625) |
---|
118 |
rlx=rlx*0.5; |
---|
119 |
if(displ==2) |
---|
120 |
disp(['relaxing computation...factor:' num2str(rlx,4)]); |
---|
121 |
end |
---|
122 |
|
---|
123 |
count=0; |
---|
124 |
a1(1:n)=0.0; |
---|
125 |
b1(1:n)=0.0; |
---|
126 |
else |
---|
127 |
if(displ==2) |
---|
128 |
warning('computation fully relaxed..bailing out'); |
---|
129 |
end |
---|
130 |
keepgoing=0; |
---|
131 |
break; |
---|
132 |
end |
---|
133 |
else |
---|
134 |
a1=a1(1:n)+rlx*a2; |
---|
135 |
b1=b1(1:n)+rlx*b2; |
---|
136 |
end |
---|
137 |
end |
---|
138 |
|
---|
139 |
error=Z-a2(1:n)*X(1:n,:)-b2(1:n)*Y(1:n,:); |
---|
140 |
|
---|
141 |
AIC(n)=M*(log(2*pi*var(error))+1)+4*n+2; |
---|
142 |
|
---|
143 |
if(n>1) |
---|
144 |
if((AIC(n)>AIC(n-1))| isnan(AIC(n))) |
---|
145 |
keepgoing=0; |
---|
146 |
end |
---|
147 |
end |
---|
148 |
|
---|
149 |
|
---|
150 |
a1held(n,1:n)=a1(1:n); |
---|
151 |
b1held(n,1:n)=b1(1:n); |
---|
152 |
best=n; |
---|
153 |
|
---|
154 |
if~(keepgoing) |
---|
155 |
if(n>1) |
---|
156 |
a1=a1held(n-1,1:n-1); |
---|
157 |
b1=b1held(n-1,1:n-1); |
---|
158 |
best=n-1; |
---|
159 |
else |
---|
160 |
a1=0.0; |
---|
161 |
b1=0.0; |
---|
162 |
end |
---|
163 |
end |
---|
164 |
|
---|
165 |
else |
---|
166 |
keepgoing=0; |
---|
167 |
end |
---|
168 |
|
---|
169 |
end |
---|
170 |
if(displ==2) |
---|
171 |
disp(['best: ' num2str(best)]); |
---|
172 |
end |
---|
173 |
|
---|
174 |
G=exp(a1*cosn(1:best,:)+b1*sinn(1:best,:))'; |
---|
175 |
|
---|
176 |
SG=G/(sum(G)*ddir); |
---|
177 |
|
---|
178 |
S(ff,1:ddirs)=Ss(1,ff)*SG'; |
---|
179 |
|
---|
180 |
end |
---|
181 |
|
---|
182 |
warning on; |
---|