Match the amplitude and phase and solve the antenna gain as a complex number. 2 points:
As in the case of "Solving the amplitude", the amplitude uses appropriate scaling.
\hat{V}_k = \boldsymbol{G}_j \boldsymbol{G}^*_i \tag{6.1}
Let's solve for $ \ boldsymbol {G} $. The unknown vector is $ 2N_a, such as $ \ boldsymbol {G} = (reG_0, reG_1, reG_2, \ dots, reG_ {Na-1}, imG_1, imG_2, \ dots, imG_ {Na-1}) $. --Represented by 1 $ real number. Let's write the real and imaginary parts of $ as $ R and I $, respectively, like $ \ boldsymbol {G} = R + iI. Since the observation equation is a nonlinear equation containing unknowns in the matrix $ P $, it is solved by the iterative method. If the initial value is $ \ boldsymbol {G ^ 0} $, the residual vector is $ r_k = \ hat {V} _k --G ^ 0_j G_i ^ {* 0} $, and the correction amount vector for it $ \ boldsymbol {c } $ Satisfies $ \ boldsymbol {r} = P \ boldsymbol {c} $ and is obtained by $ \ boldsymbol {c} = (P ^ TP) ^ {-1} P ^ T \ boldsymbol {r} $ .. When you write out the components of the matrix $ P $,
P = \left( \begin{array}{ccccccccc}
R_1 & R_0 & 0 & 0 & \cdots & I_0 & 0 & 0 & \cdots \\
R_2 & 0 & R_0 & 0 & \cdots & 0 & I_0 & 0 & \cdots \\
0 & R_2 & R_1 & 0 & \cdots & I_2 & I_0 & 0 & \cdots \\
R_3 & 0 & 0 & R_0 & \cdots & 0 & 0 & I_0 & \cdots \\
0 & R_3 & 0 & R_1 & \cdots & I_3 & 0 & I_1 & \cdots \\
0 & 0 & R_3 & R_2 & \cdots & 0 & I_3 & I_2 & \cdots \\
\vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \ddots \\
I_1 & -I_0 & 0 & 0 & \cdots & R_0 & 0 & 0 & \cdots \\
I_2 & 0 & -I_0 & 0 & \cdots & 0 & R_0 & 0 & \cdots \\
0 & I_2 & -I_1 & 0 & \cdots & -R_2 & R_0 & 0 & \cdots \\
I_3 & 0 & 0 & -I_0 & \cdots & 0 & 0 & R_0 & \cdots \\
0 & I_3 & 0 & -I_1 & \cdots & -R_3 & 0 & R_1 & \cdots \\
0 & 0 & I_3 & -I_2 & \cdots & 0 & -R_3 & R_2 & \cdots \\
\vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \ddots \\
\end{array} \right) = \left( \begin{array}{cc}
A & B \\
C & D \\
\end{array} \right) \tag{6.2}
And is represented by four submatrix $ A, B, C, D $.
A = \left( \begin{array}{ccccc}
R_1 & R_0 & 0 & 0 & \cdots \\
R_2 & 0 & R_0 & 0 & \cdots \\
0 & R_2 & R_1 & 0 & \cdots \\
R_3 & 0 & 0 & R_0 & \cdots \\
0 & R_3 & 0 & R_1 & \cdots \\
0 & 0 & R_3 & R_2 & \cdots \\
\vdots & \vdots & \vdots & \vdots & \ddots \\
\end{array} \right) \\
B = \left( \begin{array}{cccc}
I_0 & 0 & 0 & \cdots \\
0 & I_0 & 0 & \cdots \\
I_2 & I_0 & 0 & \cdots \\
0 & 0 & I_0 & \cdots \\
I_3 & 0 & I_1 & \cdots \\
0 & I_3 & I_2 & \cdots \\
\vdots & \vdots & \vdots & \ddots \\
\end{array} \right) \\
C = \left( \begin{array}{ccccc}
I_1 & -I_0 & 0 & 0 & \cdots \\
I_2 & 0 & -I_0 & 0 & \cdots \\
0 & I_2 & -I_1 & 0 & \cdots \\
I_3 & 0 & 0 & -I_0 & \cdots \\
0 & I_3 & 0 & -I_1 & \cdots \\
0 & 0 & I_3 & -I_2 & \cdots \\
\vdots & \vdots & \vdots & \vdots & \ddots \\
\end{array} \right) \\
D = \left( \begin{array}{cccc}
R_0 & 0 & 0 & \cdots \\
0 & R_0 & 0 & \cdots \\
-R_2 & R_0 & 0 & \cdots \\
0 & 0 & R_0 & \cdots \\
-R_3 & 0 & R_1 & \cdots \\
0 & -R_3 & R_2 & \cdots \\
\vdots & \vdots & \vdots & \ddots \\
\end{array} \right) \\
It is expressed as a submatrix as follows.
P^T P = \left( \begin{array}{cc}
A^T A + C^TC & A^T B + C^T D \\
B^T A + D^T C & B^TB + D^T D
\end{array} \right) \tag{6.3}
The calculation of $ A ^ TA $ is done by the formula 5.7 of "Solving the amplitude", so use it as it is. $ C ^ TC $ is very similar to $ A ^ TA $
C^TC = \left( \begin{array}{ccccc}
\sum_k I^2_k - I^2_0 & -I_0 I_1 & -I_0 I_2 & -I_0 I_3 & \cdots \\
-I_0 I_1 & \sum_k I^2_k - I^2_1 & -I_1 I_2 & -I_1 I_3 & \cdots \\
-I_0 I_2 & -I_1 I_2 & \sum_k I^2_k - I^2_2 & -I_2 I_3 & \cdots \\
-I_0 I_3 & -I_1 I_3 & -I_2 I_3 & \sum_k I^2_k - I^2_3 & \cdots \\
\vdots & \vdots & \vdots & \vdots &\ddots
\end{array} \right) \tag{6.4}
It will be. The calculations for $ B ^ TA $ and $ D ^ TC $ are as follows.
B^TA = \left( \begin{array}{ccccc}
R_1 I_0 & \boldsymbol{R}\cdot\boldsymbol{I} - R_1 I_1 & R_1 I_2 & R_1 I_3 & \cdots \\
R_2 I_0 & R_2 I_1 & \boldsymbol{R}\cdot\boldsymbol{I} - R_2 I_2 & R_2 I_3 & \cdots \\
R_3 I_0 & R_3 I_1 & R_3 I_2 & \boldsymbol{R}\cdot\boldsymbol{I} - R_3 I_3 & \cdots \\
\vdots & \vdots & \vdots & \vdots &\ddots
\end{array} \right) \tag{6.5}
D^TC = \left( \begin{array}{ccccc}
R_0 I_1 & -\boldsymbol{R}\cdot\boldsymbol{I} + R_1 I_1 & R_2 I_1 & R_3 I_1 & \cdots \\
R_0 I_2 & R_1 I_2 & -\boldsymbol{R}\cdot\boldsymbol{I} + R_2 I_2 & R_3 I_2 & \cdots \\
R_0 I_3 & R_1 I_3 & R_2 I_3 & -\boldsymbol{R}\cdot\boldsymbol{I} + R_3 I_3 & \cdots \\
\vdots & \vdots & \vdots & \vdots &\ddots
\end{array} \right) \tag{6.6}
Note that $ I_0 = 0 $. It is also calculated by $ A ^ TB = (B ^ TA) ^ T $, $ C ^ TD = (D ^ TC) ^ T $. Summarizing the above, if you write a function to get $ P ^ TP $ in Python, it will be as follows.
def ATAmatrix(Gain): # Gain is a vector of antenna-based gain amplitude (real)
antNum = len(Gain); normG = Gain.dot(Gain)
PTP = np.zeros([antNum, antNum]) + Gain
for ant_index in range(antNum):
PTP[ant_index,:] *= Gain[ant_index]
PTP[ant_index, ant_index] = normG - Gain[ant_index]**2
#
return PTP
#
def CTCmatrix(Gain): # Gain is a vector of antenna-based gain amplitude (imag)
antNum = len(Gain); normG = Gain.dot(Gain)
PTP = np.zeros([antNum, antNum]) + Gain
for ant_index in range(antNum):
PTP[ant_index,:] *= (-Gain[ant_index])
PTP[ant_index, ant_index] = normG - Gain[ant_index]**2
#
return PTP
#
def ATBmatrix(Gain): # Gain is a vector of antenna-based complex gain
antNum = len(Gain); normG = Gain.real.dot(Gain.imag)
PTP = np.zeros([antNum, antNum]) + Gain.imag
for ant_index in range(antNum):
PTP[ant_index,:] = Gain.real[ant_index]* Gain.imag + Gain.imag[ant_index]* Gain.real
PTP[ant_index, ant_index] = 0.0
#
return PTP[1:antNum]
#
def PMatrix(CompSol):
antNum = len(CompSol); matSize = 2*antNum-1
PM = np.zeros([matSize, matSize])
PM[0:antNum][:,0:antNum] = ATAmatrix(CompSol.real) + CTCmatrix(CompSol.imag) # ATA + CTC
PM[antNum:matSize][:,antNum:matSize] = ATAmatrix(CompSol.imag)[1:antNum][:,1:antNum] + CTCmatrix(CompSol.real)[1:antNum][:,1:antNum] # BTB + DTD
PM[antNum:matSize][:,0:antNum] = ATBmatrix(CompSol) # ATB + CTD
PM[0:antNum][:,antNum:matSize] = PM[antNum:matSize][:,0:antNum].T # BTA + DTC
return PM
#
The argument CompSol is a vector of the initial antenna gain (complex number). The return matrix PM has real components and a size of $ (2N_a ―― 1, 2N_a ―― 1) $.
Then calculate $ P ^ T r $. The following Python function PTdotR () takes the antenna gain (complex number) initial value vector CompSol and the residual vector (complex number) Cresid as arguments, and calculates and returns $ P ^ T r $.
def PTdotR(CompSol, Cresid):
antNum = len(CompSol)
blNum = antNum* (antNum-1) / 2
ant0, ant1= np.array(ANT0[0:blNum]), np.array(ANT1[0:blNum])
PTR = np.zeros(2*antNum)
for ant_index in range(antNum):
index0 = np.where(ant0 == ant_index)[0].tolist()
index1 = np.where(ant1 == ant_index)[0].tolist()
PTR[range(ant_index)] += (CompSol[ant_index].real* Cresid[index0].real + CompSol[ant_index].imag* Cresid[index0].imag)
PTR[range(ant_index+1,antNum)] += (CompSol[ant_index].real* Cresid[index1].real - CompSol[ant_index].imag* Cresid[index1].imag)
PTR[range(antNum, antNum+ant_index)] += (CompSol[ant_index].imag* Cresid[index0].real - CompSol[ant_index].real* Cresid[index0].imag)
PTR[range(antNum+ant_index+1,2*antNum)] += (CompSol[ant_index].imag* Cresid[index1].real + CompSol[ant_index].real* Cresid[index1].imag)
#
return PTR[range(antNum) + range(antNum+1, 2*antNum)]
#
The function gainComplex () that solves the antenna complex gain by iterative method using PMatrix () and PTdotR () above is shown below. The argument bl_vis is a baseline-based visibility (complex number) and is assumed to follow Canonical ordering. The niter is the number of iterations, but about 2 will be enough to converge. Since $ P ^ {T} P $ is a real symmetric matrix, we will solve the equation by obtaining the lower triangular matrix L by Cholesky decomposition.
def gainComplex( bl_vis, niter=2 ):
blNum = len(bl_vis)
antNum = Bl2Ant(blNum)[0]
ant0, ant1, kernelBL = ANT0[0:blNum], ANT1[0:blNum], KERNEL_BL[range(antNum-1)].tolist()
CompSol = np.zeros(antNum, dtype=complex)
#---- Initial solution
CompSol[0] = sqrt(abs(bl_vis[0])) + 0j
CompSol[1:antNum] = bl_vis[kernelBL] / CompSol[0]
#---- Iteration
for iter_index in range(niter):
PTP = PMatrix(CompSol)
L = np.linalg.cholesky(PTP) # Cholesky decomposition
Cresid = bl_vis - CompSol[ant0]* CompSol[ant1].conjugate()
t = np.linalg.solve(L, PTdotR(CompSol, Cresid))
correction = np.linalg.solve(L.T, t)
CompSol = CompSol + correction[range(antNum)] + 1.0j* np.append(0, correction[range(antNum, 2*antNum-1)])
#
return CompSol
#
that's all.
Recommended Posts