[PYTHON] Résolution du gain complexe des observations d'interféromètre

Faites correspondre l'amplitude et la phase et résolvez le gain d'antenne sous la forme d'un nombre complexe. 2 points:

1. Nombres complexes, c'est-à-dire "réels" au lieu de "amplitude" et "phase" Résolvez comme "partie et partie imaginaire". Puisque les parties réelle et imaginaire de la visibilité sont censées être normalement distribuées, cela empêche la méthode des moindres carrés d'être biaisée (Amplitude est une distribution de Rice. Par conséquent, la méthode des moindres carrés est biaisée).
2. La phase de refant est mise à 0 et ne perd pas sa généralité. En d'autres termes, la partie imaginaire du gain de référence $ \ boldsymbol {G} _0 $ est fixée à 0.

Comme dans le cas de "Résolution de l'amplitude", l'amplitude utilise une mise à l'échelle appropriée.

\hat{V}_k = \boldsymbol{G}_j \boldsymbol{G}^*_i \tag{6.1}

Résolvons pour $ \ boldsymbol {G} $. Le vecteur inconnu est $ 2N_a, tel que $ \ boldsymbol {G} = (reG_0, reG_1, reG_2, \ dots, reG_ {Na-1}, imG_1, imG_2, \ dots, imG_ {Na-1}) $. - Représenté par un nombre réel de 1 $. Écrivons respectivement les parties réelle et imaginaire de $ comme $ R et I $, comme $ \ boldsymbol {G} = R + iI. Puisque l'équation observée est une équation non linéaire contenant des inconnues dans la matrice $ P $, elle est résolue par la méthode itérative. Si la valeur initiale est $ \ boldsymbol {G ^ 0} $, le vecteur résiduel est $ r_k = \ hat {V} _k --G ^ 0_j G_i ^ {* 0} $, et le vecteur de montant de correction $ \ boldsymbol {c } $ Satisfait $ \ boldsymbol {r} = P \ boldsymbol {c} $ et est obtenu par $ \ boldsymbol {c} = (P ^ TP) ^ {-1} P ^ T \ boldsymbol {r} $ .. Lorsque vous écrivez les composants de la matrice $ 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}

Et est représenté par quatre sous-matrices $ 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) \\

Économie de travail dans le calcul matriciel

Il est exprimé comme une sous-matrice comme suit.

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}

Le calcul de $ A ^ TA $ se fait par la formule 5.7 de "Résoudre l'amplitude", alors utilisez-le tel quel. $ C ^ TC $ est très similaire à $ 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}

Ce sera. Le calcul de $ B ^ TA $ et $ D ^ TC $ est le suivant.

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}

Notez que $ I_0 = 0 $. Il est également calculé par $ A ^ TB = (B ^ TA) ^ T $, $ C ^ TD = (D ^ TC) ^ T $. Pour résumer ce qui précède, si vous écrivez une fonction pour obtenir $ P ^ TP $ en Python, ce sera comme suit.

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
#

L'argument CompSol est le vecteur du gain d'antenne initial (nombre complexe). La matrice de valeur de retour PM a des composants réels et une taille de $ (2N_a --1, 2N_a --1) $.

Puis calculez $ P ^ T r $. La fonction Python suivante PTdotR () prend le vecteur de valeur initiale du gain d'antenne (nombre complexe) CompSol et le vecteur résiduel (nombre complexe) Cresid comme arguments, et calcule et renvoie $ 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)]
#

Code récapitulatif

La fonction gainComplex () qui résout le gain du complexe d'antenne par la méthode itérative utilisant PMatrix () et PTdotR () ci-dessus est illustrée ci-dessous. L'argument bl_vis est une visibilité basée sur la ligne de base et est supposé suivre l'ordre canonique. Le niter est le nombre d'itérations, mais environ 2 fois suffiront pour converger. Puisque $ P ^ {T} P $ est une matrice symétrique réelle, nous résoudrons l'équation en obtenant la matrice triangulaire inférieure L par décomposition de Cholesky.

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
#

c'est tout.

Retour à Mokuji