DOA算法3：Matrix Pencil

注：本博文为本人阅读论文、文章后的原创笔记，未经授权不允许任何转载或商用行为，否则一经发现本人保留追责权利。有问题可留言联系，欢迎指摘批评，共同进步！！！

为什么需要Matrix Pencil算法？

注意到，之前所述的DOA经典算法如MUSIC、ESPRIT算法都需要首先估计出相关矩阵$\mathbf{R}$。而要估计出这个矩阵是非常消耗计算量的：对于数据向量$\mathbf{x}$，需要至少$K$次的快照，并且要求$K>2N$。

算法介绍

Matrix Pencil最初是用于确定系统极点的。
在第$n$个时隙所接收到的信号可以建模为：
$$x_n=\sum _{m=1}^M{\mathbf{A}}_m{z}_m^n+n_n$$
其中$z_m=e^{jkdcos{\varphi }_m△t}$就是系统的极点（不用记），$n_n$代表加性高斯白噪声。目标就是给定接收信号$x_n$来估计极点$z_m$。

数理基础——广义特征值

存在$\lambda$，使得式$\mathbf{Ax}=\lambda \mathbf{Bx}$存在非0解向量，那么$\lambda$就称作矩阵$\mathbf{A}$对$\mathbf{B}$的广义特征值；相应的$\mathbf{x}$是广义特征向量。

• 显然，当$\mathbf{B}=\mathbf{E}$时，就是普通的求特征值问题；
• 当$\mathbf{B}\ne \mathbf{E}$时，就是广义特征值问题。

同理，等式可以化为如下形式：
$$(\mathbf{A}-\lambda \mathbf{B})\mathbf{x}=0$$
这是个齐次方程组

齐次方程组要有非零解，那么系数矩阵所形成的的行列式必定为0

根据要求，那么使得$det(\mathbf{A}-\lambda \mathbf{B})=0$即可解出广义特征值。

算法原理

与ESPRIT算法相似，Matrix Pencil也需要分解矩阵，它分解的就是接收信号矩阵$\mathbf{X}$. 首先我们定义两个$(N-L)\times L$大小的矩阵${\mathbf{X}}_0$和${\mathbf{X}}_1$：
$${\mathbf{X}}_0=\left[\begin{array}{cccc}{x}_0& x_1& \cdots & x_{L-1}\\ x_1& x_2& \cdots & x_L\\ ⋮& ⋮& \ddots & ⋮\\ x_{N-L-1}& x_{N-L}& \cdots & x_{N-2}\end{array}\right],{\mathbf{X}}_1=\left[\begin{array}{cccc}{x}_1& x_2& \cdots & x_L\\ x_2& x_3& \cdots & x_{L+1}\\ ⋮& ⋮& \ddots & ⋮\\ x_{N-L}& x_{N-L+1}& \cdots & x_{N-1}\end{array}\right]$$
其中，$\mathbf{L}$必须满足的条件是：
$$\{\begin{array}{l}M\le L\le N-LN是偶数\\ M\le L\le N-L+1N是奇数\end{array}$$
在Matrix Pencil矩阵中，我们可以得到如下矩阵等式：
$$\begin{array}{rl}{\mathbf{X}}_0& ={\mathbf{Z}}_1\mathbf{A}{\mathbf{Z}}_2,\\ {\mathbf{X}}_1& ={\mathbf{Z}}_1\mathbf{A\Phi }{\mathbf{Z}}_2,\end{array}$$
其中矩阵$\mathbf{\Phi }$就是与ESPRIT算法中相同的对角矩阵$\mathbf{\Phi }$.上面等式中的四个矩阵分别给定为：
$$\begin{array}{rl}{\mathbf{Z}}_1& ={\left[\begin{array}{cccc}1& 1& \cdots & 1\\ z_1& z_2& \cdots & z_M\\ ⋮& ⋮& \ddots & ⋮\\ z_1^{(N-L-1)}& z_2^{(N-L-1)}& \cdots & z_M^{(N-L-1)}\end{array}\right]}_{(N-L)\times M}\\ {\mathbf{Z}}_2& ={\left[\begin{array}{cccc}1& z_1& \cdots & z_1^{L-1}\\ 1& z_2& \cdots & z_2^{L-1}\\ ⋮& ⋮& \ddots & ⋮\\ 1& z_M& \cdots & z_M^{L-1}\end{array}\right]}_{M\times L}\\ \mathbf{\Phi }& ={\left[\begin{array}{cccc}{z}_1& 0& \cdots & 0\\ 0& z_2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & z_M\end{array}\right]}_{M\times M}\\ \mathbf{A}& ={\left[\begin{array}{cccc}{\alpha }_1& 0& \cdots & 0\\ 0& {\alpha }_2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {\alpha }_M\end{array}\right]}_{M\times M}\end{array}$$
在无噪声的情况下，$L$也满足约束条件，我们可以看到：矩阵${\mathbf{Z}}_1$和${\mathbf{Z}}_2$都是满秩矩阵，因此矩阵${\mathbf{X}}_0$和${\mathbf{X}}_1$的秩与中间的对角阵$\mathbf{\Phi }$和$\mathbf{A\Phi }$相同，都是$M$.
考虑矩阵${\mathbf{X}}_1-\lambda {\mathbf{X}}_0={\mathbf{Z}}_1\mathbf{A}[\mathbf{\Phi }-\lambda \mathbf{I}]{\mathbf{Z}}_2$。注意到中间的$\mathbf{\Phi }-\lambda \mathbf{I}$部分：

• 当$\lambda \ne z_m(m\in [1,M])$时，$\mathbf{\Phi }-\lambda \mathbf{I}$的秩还是$M$，因此${\mathbf{X}}_1-\lambda {\mathbf{X}}_0$的秩还是$M$；
• 当$\lambda =z_m(m\in [1,M])$时，$\mathbf{\Phi }-\lambda \mathbf{I}$的秩就减小为$M-1$。此时矩阵不满秩，那么行列式$det({\mathbf{X}}_1-\lambda {\mathbf{X}}_0)=0$恒成立。

回顾到广义特征值的性质，我们可以看到，$\lambda =z_m(m\in [1,M])$时，刚好可以作为矩阵${\mathbf{X}}_1$对${\mathbf{X}}_0$的广义特征值.

矩阵${\mathbf{X}}_0$对${\mathbf{X}}_1$的广义特征值就是求 使得式${\mathbf{X}}_0\mathbf{x}=\lambda {\mathbf{X}}_1\mathbf{x}$存在非零解的特征值$\lambda$。由推导可以知道，就是求 使得行列式$det({\mathbf{X}}_1-\lambda {\mathbf{X}}_0)=0$恒成立的$\lambda$。根据上面所述，当$\lambda =z_m(m\in [1,M])$时，可满足条件。因此，我们可以通过找出矩阵${\mathbf{X}}_1$对${\mathbf{X}}_0$的$M$个广义特征值来作为对$z_m(m\in [1,M])$的估计。

算法步骤

1. 给定$N$和$M$，并且选择合适的$L$;
2. 创建矩阵${\mathbf{X}}_0$和${\mathbf{X}}_1$；
3. 求矩阵${\mathbf{X}}_1$对${\mathbf{X}}_0$的$M$个广义特征值${\lambda }_m(m\in [1,M])$，这些特征值就是对$z_m(m\in [1,M])$的估计；
4. 用下式计算DOA:
   

参考文献（仅写文章的标题，以做记录）

[1]. Direction of Arrival Estimation