In this chapter, we will describe how to solve only the amplitude of the antenna gain using the amplitude of the visibility obtained by the interferometer observation. However, instead of actually using it, we recommend Solving complex gain by matching amplitude and phase. This is because the amplitude follows a Rice distribution rather than a normal distribution, so when the signal-to-noise ratio is small, the amplitude tends to be biased to a large extent (while complex gains are normally distributed). If $ SNR> 5 $, it will be close to a normal distribution, so you can outline the method in this chapter. The significance of this chapter is to make it easier to understand How to solve antenna gain as a complex number.
Antenna gain amplitude
|\hat{V}|_k = |G|_j |G|_j |V|_k\tag{5.1}
On the right side
|\hat{V}|_k = |G|_j |G|_j \tag{5.2}
Taking the logarithm of equation (5.2)
\log |\hat{V}|_k = \log |G|_j + \log |G|_j \tag{5.3}
And can be linearized. That is,
\left(
\begin{array}{c}
\log |\hat{V}_0| \\
\log |\hat{V}_1| \\
\log |\hat{V}_2| \\
\log |\hat{V}_3| \\
\log |\hat{V}_4| \\
\log |\hat{V}_5| \\
\vdots
\end{array}
\right) = \left( \begin{array}{ccccc}
1 & 1 & 0 & 0 & \cdots \\
1 & 0 & 1 & 0 & \cdots \\
0 & 1 & 1 & 0 & \cdots \\
1 & 0 & 0 & 1 & \cdots \\
0 & 1 & 0 & 1 & \cdots \\
0 & 0 & 1 & 1 & \cdots \\
\vdots & \vdots & \vdots & \ddots
\end{array} \right)
\left(
\begin{array}{c}
\log |{G}_0| \\
\log |{G}_1| \\
\log |{G}_2| \\
\log |{G}_3| \\
\vdots
\end{array}
\right) \tag{5.4}
is. The procession in the middle
Calculate $ P ^ T P $ by hand.
P^TP = \left( \begin{array}{ccccc}
N_a - 1 & 1 & 1 & 1 & \cdots \\
1 & N_a - 1 & 1 & 1 & \cdots \\
1 & 1 & N_a-1 & 1 & \cdots \\
1 & 1 & 1 & N_a-1 & \cdots \\
\vdots & \vdots & \vdots & \vdots &\ddots
\end{array} \right) \tag{5.5}
And its inverse matrix
(P^TP)^{-1} = \frac{1}{2(N_a-1)(N_a - 2)} \left( \begin{array}{ccccc}
2N_a -3 & -1 & -1 & -1 & \cdots \\
-1 & 2N_a -3 & -1 & -1 & \cdots \\
-1 & -1 & 2N_a -3 & -1 & \cdots \\
-1 & -1 & -1 & 2N_a -3 & \cdots \\
\vdots & \vdots & \vdots & \vdots &\ddots
\end{array} \right) \tag{5.6}
It will be. In short, the diagonal component is $ \ frac {2N_a -3} {2 (N_a-1) (N_a --2)} $, and the off-diagonal component is $-\ frac {1} {2 (N_a-1) (N_a-) 2)} $. This matrix can also be coded in one line in Python.
import numpy as np
PTP_inv = ((2.0* antNum - 2.0)* np.diag(np.ones(antNum)) - 1.0) / (2.0* (antNum - 1.0)* (antNum - 2.0))
The calculation of $ P ^ T \ boldsymbol {V} $ can be done in much the same way as the code when solving the delay.
ant0, ant1 = np.array(ANT0[0:blNum]), np.array(ANT1[0:blNum])
PTV = np.zeros(antNum)
for ant_index in range(antNum):
index0, index1 = np.where(ant0 == ant_index)[0].tolist(), np.where(ant1 == ant_index)[0].tolist()
PTV[ant_index] += np.sum(log_bl[index0]) + np.sum(log_bl[index1])
#
When solving equation (5.2) without taking the logarithm, it becomes a non-linear equation. Assuming that the initial value of the antenna gain amplitude is $ G ^ 0 $, the baseline-based residual is $ r_k = \ hat {V} _k --G ^ 0_j G_i ^ 0 $. The correction vector $ \ boldsymbol {c} $ given to $ G ^ 0 $ to minimize the residual norm $ \ displaystyle \ sum_k r_kr ^ * _k $ is
\left(
\begin{array}{c}
r_0 \\
r_1 \\
r_2 \\
r_3 \\
r_4 \\
r_5 \\
\vdots
\end{array}
\right) = \left( \begin{array}{ccccc}
G_1 & G_0 & 0 & 0 & \cdots \\
G_2 & 0 & G_0 & 0 & \cdots \\
0 & G_2 & G_1 & 0 & \cdots \\
G_3 & 0 & 0 & G_0 & \cdots \\
0 & G_3 & 0 & G_1 & \cdots \\
0 & 0 & G_3 & G_2 & \cdots \\
\vdots & \vdots & \vdots & \ddots
\end{array} \right)
\left(
\begin{array}{c}
c_0 \\
c_1 \\
c_2 \\
c_3 \\
\vdots
\end{array}
\right) \tag{5.7}
Satisfies the equation. The matrix in the middle is non-linear because the $ P $ element contains $ G $ to be solved. Solving equation (5.7) for $ c $ yields $ \ boldsymbol {c} = (P ^ T P) ^ {-1} P ^ T \ boldsymbol {r} $. The $ \ boldsymbol {c} $ thus obtained is used to improve $ \ boldsymbol {G} $.
\boldsymbol{G} \leftarrow \boldsymbol{G} + \boldsymbol{c}
Repeat this until $ \ boldsymbol {G} $ converges to get a solution. Convergence can be determined by making the norm of $ \ boldsymbol {c} $ small enough.
As usual, $ P ^ T P $ is calculated manually.
P^TP = \left( \begin{array}{ccccc}
\sum_k G^2_k - G^2_0 & G_0 G_1 & G_0 G_2 & G_0 G_3 & \cdots \\
G_0 G_1 & \sum_k G^2_k - G^2_1 & G_1 G_2 & G_1 G_3 & \cdots \\
G_0 G_2 & G_1 G_2 & \sum_k G^2_k - G^2_2 & G_2 G_3 & \cdots \\
G_0 G_3 & G_1 G_3 & G_2 G_3 & \sum_k G^2_k - G^2_3 & \cdots \\
\vdots & \vdots & \vdots & \vdots &\ddots
\end{array} \right) \tag{5.8}
Implement this in Python.
def PTPmatrix(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
#
The inverse matrix $ (P ^ T P) ^ {-1} $ cannot be written in a simple form. Since it is a real symmetric matrix, we will use the Cholesky decomposition to solve the equation. Calculate $ \ boldsymbol {y} = P ^ {T} \ boldsymbol {r} $ in advance.
y = np.zeros(antNum)
for ant_index in range(antNum):
index0, index1 = np.where(ant0 == ant_index)[0].tolist(), np.where(ant1 == ant_index)[0].tolist()
y[ant_index] += antGain[ant1[index0]].dot(resid[index0])
y[ant_index] += antGain[ant0[index1]].dot(resid[index1])
#
Then use the Cholesky decomposition to solve $ P ^ T P \ boldsymbol {c} = \ boldsymbol {y} $ for $ \ boldsymbol {c} $.
L = np.linalg.cholesky(PTP)
t = np.linalg.solve(L, y)
correction = np.linalg.solve(L.T, t)
Add the $ \ boldsymbol {c} $ (correction) thus obtained to $ \ boldsymbol {G} $ (antGain), improve it, and iterate until it converges.
Summarizing the above, the Python code for obtaining the antenna-based gain amplitude from the baseline-based visibility amplitude is shown below. The function logamp_solve () is a logarithmically linearized solution that takes a vector of visibility amplitude as an argument and returns a vector of antenna gain amplitude. The function clamp_solve () is a non-linear solution and has the same specifications as logamp_solve () for its arguments and return values. Internally, logamp_solve (bl_amp) is called to get the initial value of the antenna gain amplitude. It also calls the function PTPmatrix () to calculate $ (P ^ T P) ^ {-1} $.
def logamp_solve(bl_amp):
blnum = len(bl_amp); log_bl = np.log(bl_amp)
antNum = Bl2Ant(blnum)[0]
ant0, ant1 = np.array(ANT0[0:blNum]), np.array(ANT1[0:blNum])
PTP_inv = ((2.0* antNum - 2.0)* np.diag(np.ones(antNum)) - 1.0) / (2.0* (antNum - 1.0)* (antNum - 2.0))
PTV = np.zeros(antNum)
for ant_index in range(antNum):
index0, index1 = np.where(ant0 == ant_index)[0].tolist(), np.where(ant1 == ant_index)[0].tolist()
PTV[ant_index] += np.sum(log_bl[index0]) + np.sum(log_bl[index1])
#
return np.exp(PTP_inv.dot(PTV))
#
def PTPmatrix(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 clamp_solve(bl_amp, niter=2):
blnum = len(bl_amp)
antGain = logamp_solve(bl_amp); antNum = len(antGain)
ant0, ant1 = np.array(ANT0[0:blNum]), np.array(ANT1[0:blNum])
for iter_index in range(niter):
resid = bl_amp - antGain[ant0]* antGain[ant1]
y = np.zeros(antNum)
for ant_index in range(antNum):
index0, index1 = np.where(ant0 == ant_index)[0].tolist(), np.where(ant1 == ant_index)[0].tolist()
y[ant_index] += antGain[ant1[index0]].dot(resid[index0])
y[ant_index] += antGain[ant0[index1]].dot(resid[index1])
#
L = np.linalg.cholesky(PTPmatrix(antGain))
t = np.linalg.solve(L, y)
correction = np.linalg.solve(L.T, t)
antGain += correction; antGain = abs(antGain)
#
return antGain
#
that's all. Actually, it is better to solve amplitude and phase as complex numbers at the same time rather than solving only amplitude. Return to the lottery
Recommended Posts