Solving the amplitude of interferometer observations

Solve the amplitude

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|G|And baseline-based visibility amplitude|V|The relationship with is as follows.

|\hat{V}|_k = |G|_j |G|_j |V|_k\tag{5.1}

On the right side|V|_kIs the true visibility amplitude, on the left side|\hat{V}|_kIs the observed visibility amplitude. The true visibility amplitude is not known until it is calibrated and determined.|V|_kDivided by|\hat{V}|_k / |\hat{V}|_kTreating as the left side does not lose generality. in this case,|G|Is an appropriately scaled value and has only a meaning as a relative amplitude. In this case the observation equation(5.1)Is the following(5.2)It looks like.

|\hat{V}|_k = |G|_j |G|_j \tag{5.2}

Simple solution by logarithm

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 middlePToki, on the left side\log |V|Vector\boldsymbol{V},On the right side\log |G|Vector\boldsymbol{G}If you say\boldsymbol{V} = P \boldsymbol{G}So\boldsymbol{G}To solve\boldsymbol{G} = (P^T P)^{-1} P^T \boldsymbol{V}is.

Labor saving in calculation

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])
#

Non-linear solution

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.

Labor saving in calculation

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.

Summary code

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

Solving the amplitude of interferometer observations
Solving the phase of interferometer observations
Solving the complex gain of interferometer observations
Solve the delay of interferometer observation
Solving the equation of motion in Python (odeint)
The beginning of cif2cell
The meaning of self
the zen of Python
The story of sys.path.append ()
python chrome driver ver. Solving the problem of difference
AttributeError: The story of solving module'tensorflow' has no attribute'log'.
Revenge of the Types: Revenge of types
Align the version of chromedriver_binary
Scraping the result of "Schedule-kun"
10. Counting the number of lines
The story of building Zabbix 4.4
Towards the retirement of Python2
[Apache] The story of prefork
Compare the fonts of jupyter-themes
About the ease of Python
Get the number of digits
Explain the code of Tensorflow_in_ROS
Reuse the results of clustering
GoPiGo3 of the old man
Change the theme of Jupyter
The popularity of programming languages
Change the style of matplotlib
Visualize the orbit of Hayabusa2
About the components of Luigi
Connected components of the graph
Filter the output of tracemalloc
About the features of Python
Simulation of the contents of the wallet
The Power of Pandas: Python
Solving the Maze with Python-Supplement to Chapter 6 of the Algorithm Quick Reference-
Solving Knapsack Problems with Google's OR-Tools-Practicing the Basics of Combinatorial Optimization Problems