## 6 Modified Unitary R-D Root MUSIC
Normally, the standard way to multiply two matrices, particularly \(\mathbf{\hat{R}}_{xx}=\mathbf{X}\mathbf{V}\), is evidenced by Equation (37),  
\[\hat{r}_{xx}(i,j)=(\frac{1}{N})\sum_{k=1}^{N}x(i,k)v(k,j),\forall i,j=1,\dots,M. \tag{37}\]  
Since \(\mathbf{V}=\mathbf{X}^{H}\), then \(v(k,j)=\overline{x}(j,k)\); therefore, Equation (37) could be written as indicated by Equation (38),  
\[\hat{r}_{xx}(i,j)=(\frac{1}{N})\sum_{k=1}^{N}x(i,k)\overline{x}(j,k),\forall i,j=1,\dots,M. \tag{38}\]Moreover, since \(\mathbf{\hat{R}}_{xx}\) is Hermitian, which means \(\hat{r}_{xx}(j,i)=\overline{\hat{r}_{xx}}(i,j)\), thus the solution applies matrix multiplication only on its upper triangular part; therefore, Equation (38) is remodeled as the set of relations in (39),  
\[\hat{r}_{xx}(i,j) =(\frac{1}{N})\sum_{k=1}^{N}x(i,k)*\overline{x}(j,k), \tag{39}\] \[\hat{r}_{xx}(j,i) =\overline{\hat{r}_{xx}}(i,j),\] \[\forall i,j =i,\ldots,M.\]  
From Equation (39), the implemented solution estimates the covariance matrix using the same matrix \(\mathbf{X}\) twice by applying an element-wise conjugate transpose operation, thus it does not need to store \(\mathbf{X}^{H}\). Furthermore, it only computes the upper triangular part of \(\mathbf{\hat{R}}_{xx}\). To sum up, that approach saves execution time by half and memory usage in the order of \(MN\). That is an important improvement, since calculating the sample covariance matrix is the second most time-consuming operation, as shown in Section 7. The same optimization is carried out in Equation (20) for the y-axis ULA. Another minor optimization concerns the computation of the RF switch compensation matrix, which requires the inverse matrix calculation of \(\mathbf{O}\) defined in Equation (33). The Gauss-Jordan elimination is a popular algorithm to calculate the inverse of a matrix whose complex is \(\mathcal{O}(n^{3})\)[35]. However, since \(\mathbf{O}\) is a diagonal matrix in which the generic form of its elements is known, we can avoid applying the time-consuming inverse matrix calculation. From complex arithmetic, we know that \((e^{i\theta})e^{-j\theta}=1\), thus, the inverse of \(\mathbf{O}\) is in line with Equation (40),  
\[\mathbf{O}^{-1}=\text{diag}(1,e^{-j2\pi\Delta f_{c}\Delta T_{2}},\cdots,e^{-j 2\pi\Delta f_{c}\Delta T_{2M-1}})\in\mathbb{C}^{(2M-1)\times(2M-1)}. \tag{40}\]  
Since \(e^{-j2\pi\Delta f_{c}\Delta T_{2}}=e^{-j2\pi\Delta f_{c}(k-1)\Delta T_{2}}\), then \(e^{-j2\pi\Delta f_{c}\Delta T_{2}}=(e^{-j2\pi\Delta f_{c}\Delta T_{2}})^{k-1}\) for \(2\leq k\leq 2M-1\). By defining, \(z=e^{-j2\pi\Delta f_{c}\Delta T_{2}}\), the inverse of \(\mathbf{O}\) can be redefined in a more compact form as evidenced by Equation (41),  
\[\mathbf{O}^{-1}=\text{diag}(z^{0},z^{1},\ldots,z^{2M-2})\in\mathbb{C}^{(2M-1) \times(2M-1)}. \tag{41}\]  
From the redefinition in Equation (41), we can derive the Relation (42),  
\[\mathbf{O}^{-1}(k,k)=\mathbf{O}^{-1}(k-1,k-1)z,\quad 2\leq k\leq 2M-1. \tag{42}\]  
Thus, to calculate \(\mathbf{O}^{-1}\), the implemented solution only needs to set \(\mathbf{O}^{-1}(1,1)=1\), compute \(z=e^{-j2\pi\Delta f_{c}\Delta T_{2}}\), and apply Equation (42) successively for \(k=2\ldots 2M-1\) to take advantage of the previous computation. As a result, the implemented solution does not need to compute each element of \(\mathbf{O}^{-1}\) explicitly, that is \(e^{-j2\pi\Delta f_{c}\Delta T_{k}}\), which may require computing the finite Maclaurin series \(2M-2\) times to calculate all of \(2M-2\) elements, except the first, which is \(1\). Notably, the finite Maclaurin series is a well-known method to evaluate complex numbers by computers, as illustrated in Equation (43),  
\[e^{jx}=\underbrace{\left(1-\frac{x^{2}}{2!}+\frac{x^{4}}{4!}-\cdots\frac{(-1 )^{n}}{(2n)!}x^{2n}\right)}_{\cos(x)}+j\underbrace{\left(x-\frac{x^{3}}{3!}+ \frac{x^{5}}{5!}-\cdots\frac{(-1)^{n}}{(2n+1)!}x^{2n+1}\right)}_{\sin(x)}, \tag{43}\]  
where \(n\) is the number of elements.