1 |
function [S]=BDM(xps,trm,kx,Ss,pidirs,miter,displ) |
---|
2 |
|
---|
3 |
nmod=6; |
---|
4 |
|
---|
5 |
szd=size(xps,1); |
---|
6 |
freqs=size(xps,3); |
---|
7 |
ddirs=size(trm,3); |
---|
8 |
|
---|
9 |
ddir=abs(pidirs(2)-pidirs(1)); |
---|
10 |
pi=4.0*atan(1.0); |
---|
11 |
|
---|
12 |
if(sum(kx)==0) |
---|
13 |
warning('BDM method may not work with three quantity measurements'); |
---|
14 |
disp(' '); |
---|
15 |
end |
---|
16 |
|
---|
17 |
if(displ<2) |
---|
18 |
warning off; |
---|
19 |
end |
---|
20 |
|
---|
21 |
Co=real(xps); |
---|
22 |
Quad=-imag(xps); |
---|
23 |
|
---|
24 |
|
---|
25 |
for ff=1:freqs |
---|
26 |
xpsx(:,:,ff)=diag(xps(:,:,ff))*(diag(xps(:,:,ff))'); |
---|
27 |
sigCo(:,:,ff)=sqrt(0.5*(xpsx(:,:,ff)+Co(:,:,ff).^2-Quad(:,:,ff).^2)); |
---|
28 |
sigQuad(:,:,ff)=sqrt(0.5*(xpsx(:,:,ff)-Co(:,:,ff).^2+Quad(:,:,ff).^2)); |
---|
29 |
end |
---|
30 |
|
---|
31 |
for ff=1:freqs |
---|
32 |
|
---|
33 |
index=0; |
---|
34 |
for m=1:szd |
---|
35 |
for n=m:szd |
---|
36 |
expx(1:ddirs)=exp(-i*kx(m,n,ff,1:ddirs)); |
---|
37 |
Hh(1:ddirs)=trm(m,ff,1:ddirs); |
---|
38 |
Hhs(1:ddirs)=conj(trm(n,ff,1:ddirs)); |
---|
39 |
Htemp=(Hh.*Hhs.*expx); |
---|
40 |
|
---|
41 |
if(Htemp(1)~=Htemp(2)) |
---|
42 |
index=index+1; |
---|
43 |
phi(index,ff)=real(xps(m,n,ff))./(sigCo(m,n,ff)*Ss(1,ff)); |
---|
44 |
H(1:ddirs,index,ff)=real(Htemp)./sigCo(m,n,ff); |
---|
45 |
if(kx(m,n,1,1)+kx(m,n,1,2)~=0) |
---|
46 |
index=index+1; |
---|
47 |
phi(index,ff)=imag(xps(m,n,ff))./(sigQuad(m,n,ff)*Ss(1,ff)); |
---|
48 |
H(1:ddirs,index,ff)=imag(Htemp)./sigQuad(m,n,ff); |
---|
49 |
end |
---|
50 |
end |
---|
51 |
end |
---|
52 |
end |
---|
53 |
end |
---|
54 |
M=index; |
---|
55 |
|
---|
56 |
k=ddirs; |
---|
57 |
|
---|
58 |
dd=diag(ones(k,1))+diag(-2*ones(k-1,1),-1)+diag(ones(k-2,1),-2); |
---|
59 |
dd(1,k-1:k)=[1 -2]; |
---|
60 |
dd(2,k)=1; |
---|
61 |
|
---|
62 |
|
---|
63 |
for ff=1:freqs |
---|
64 |
if(displ>0) |
---|
65 |
disp(['calculating for frequency' blanks(1) num2str(ff) ' of' blanks(1) num2str(freqs)]); |
---|
66 |
end |
---|
67 |
a(1:k,1:M)=H(1:k,1:M,ff)*ddir; |
---|
68 |
A=a'; |
---|
69 |
B=phi(:,ff); |
---|
70 |
n=0; |
---|
71 |
keepgoing=1; |
---|
72 |
|
---|
73 |
while(keepgoing) |
---|
74 |
n=n+1; |
---|
75 |
if(displ>0) |
---|
76 |
disp(strcat('model :',num2str(n))); |
---|
77 |
end |
---|
78 |
if(n<=nmod) |
---|
79 |
u=0.5^n; |
---|
80 |
x(1:k,1)=log(1/(2*pi))*ones(k,1); |
---|
81 |
stddiff=1; |
---|
82 |
rlx=1.0; |
---|
83 |
count=0; |
---|
84 |
while(stddiff>0.001); |
---|
85 |
count=count+1; |
---|
86 |
F=exp(x); |
---|
87 |
E=diag(F); |
---|
88 |
A2=A*E; |
---|
89 |
B2=B-A*F+A*E*x; |
---|
90 |
Z(1:M,1:k)=A2; |
---|
91 |
Z(M+1:M+k,1:k)=u*dd; |
---|
92 |
Z(1:M,k+1)=B2; |
---|
93 |
Z(M+1:M+k,k+1)=zeros(k,1); |
---|
94 |
|
---|
95 |
[Q,U]=qr(Z); |
---|
96 |
|
---|
97 |
UZ=U; |
---|
98 |
TA=UZ(1:k,1:k); |
---|
99 |
Tb=UZ(1:k,k+1); |
---|
100 |
|
---|
101 |
x1=TA\Tb; |
---|
102 |
stddiff=std(x-x1); |
---|
103 |
x=(1-rlx)*x+rlx*x1; |
---|
104 |
if(count>miter|sum(isfinite(x))~=k) |
---|
105 |
if(rlx>0.0625) |
---|
106 |
rlx=rlx*0.5; |
---|
107 |
if(displ==2) |
---|
108 |
disp(['relaxing computation...factor:' num2str(rlx,4)]); |
---|
109 |
end |
---|
110 |
if(sum(isfinite(x))~=k) |
---|
111 |
x(1:k,1)=log(1/(2*pi))*ones(k,1); |
---|
112 |
end |
---|
113 |
count=0; |
---|
114 |
else |
---|
115 |
if(displ==2) |
---|
116 |
warning('computation fully relaxed..bailing out'); |
---|
117 |
end |
---|
118 |
if(n>1) |
---|
119 |
keepgoing=0; |
---|
120 |
end |
---|
121 |
break; |
---|
122 |
end |
---|
123 |
end |
---|
124 |
end |
---|
125 |
sig2=((norm(A2*x-B2)).^2+(u*norm(dd*x)).^2)/M; |
---|
126 |
ABIC(n)=M*(log(2*pi*sig2)+1)-k*log(u*u)+sum(log(diag(TA).^2)); |
---|
127 |
|
---|
128 |
if(n>1) |
---|
129 |
if(ABIC(n)>ABIC(n-1)) |
---|
130 |
keepgoing=0; |
---|
131 |
n=n-1; |
---|
132 |
end |
---|
133 |
end |
---|
134 |
|
---|
135 |
if(keepgoing) |
---|
136 |
xold=x; |
---|
137 |
else |
---|
138 |
x=xold; |
---|
139 |
end |
---|
140 |
|
---|
141 |
else |
---|
142 |
keepgoing=0; |
---|
143 |
end |
---|
144 |
|
---|
145 |
end |
---|
146 |
if(displ==2) |
---|
147 |
disp(['best: ' num2str(n)]); |
---|
148 |
end |
---|
149 |
G=exp(x); |
---|
150 |
S(ff,1:k)=Ss(1,ff)*G'; |
---|
151 |
end |
---|
152 |
|
---|
153 |
warning on; |
---|
154 |
|
---|
155 |
|
---|
156 |
|
---|