为什么需要Matrix Pencil算法？

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

算法介绍

Matrix Pencil最初是用于确定系统极点的。
在第nnn个时隙所接收到的信号可以建模为：
xn=∑m=1MAmzmn+nnx_n=\sum^{M}_{m=1}{\mathbf{A}_mz^n_m+n_n} xn=m=1∑MAmzmn+nn
其中zm=ejkdcosϕm△tz_m=e^{jkd\ cos{\phi_m}\bigtriangleup t}zm=ejkd cosϕm△t就是系统的极点（不用记），nnn_nnn代表加性高斯白噪声。目标就是给定接收信号xnx_nxn来估计极点zmz_mzm。

数理基础——广义特征值

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

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

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

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

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

算法原理

与ESPRIT算法相似，Matrix Pencil也需要分解矩阵，它分解的就是接收信号矩阵X\mathbf{X}X. 首先我们定义两个(N−L)×L(N-L) \times L(N−L)×L大小的矩阵X0\mathbf{X_0}X0和X1\mathbf{X_1}X1：
X0=[x0x1⋯xL−1x1x2⋯xL⋮⋮⋱⋮xN−L−1xN−L⋯xN−2],X1=[x1x2⋯xLx2x3⋯xL+1⋮⋮⋱⋮xN−LxN−L+1⋯xN−1]\mathbf{X_0}= \begin{bmatrix} x_0 & x_1 & \cdots & x_{L-1}\\ x_1 & x_2 & \cdots & x_L\\ \vdots & \vdots & \ddots & \vdots\\ x_{N-L-1} & x_{N-L} & \cdots & x_{N-2} \end{bmatrix}, \mathbf{X_1}= \begin{bmatrix} x_1 & x_2 & \cdots & x_L\\ x_2 & x_3 & \cdots & x_{L+1}\\ \vdots & \vdots & \ddots & \vdots\\ x_{N-L} & x_{N-L+1} & \cdots & x_{N-1} \end{bmatrix} X0=⎣⎡x0x1⋮xN−L−1x1x2⋮xN−L⋯⋯⋱⋯xL−1xL⋮xN−2⎦⎤,X1=⎣⎡x1x2⋮xN−Lx2x3⋮xN−L+1⋯⋯⋱⋯xLxL+1⋮xN−1⎦⎤
其中，L\mathbf{L}L必须满足的条件是：
{M≤L≤N−LN是偶数M≤L≤N−L+1N是奇数\begin{cases} M \le L \le N-L \qquad N是偶数\\ M \le L \le N-L+1 \quad N是奇数 \end{cases} {M≤L≤N−LN是偶数M≤L≤N−L+1N是奇数
在Matrix Pencil矩阵中，我们可以得到如下矩阵等式：
X0=Z1AZ2,X1=Z1AΦZ2,\begin{aligned} \mathbf{X_0} &= \mathbf{Z_1AZ_2}, \\ \mathbf{X_1} &= \mathbf{Z_1A\Phi Z_2}, \end{aligned} X0X1=Z1AZ2,=Z1AΦZ2,
其中矩阵Φ\mathbf{\Phi}Φ就是与ESPRIT算法中相同的对角矩阵Φ\mathbf{\Phi}Φ.上面等式中的四个矩阵分别给定为：
Z1=[11⋯1z1z2⋯zM⋮⋮⋱⋮z1(N−L−1)z2(N−L−1)⋯zM(N−L−1)](N−L)×MZ2=[1z1⋯z1L−11z2⋯z2L−1⋮⋮⋱⋮1zM⋯zML−1]M×LΦ=[z10⋯00z2⋯0⋮⋮⋱⋮00⋯zM]M×MA=[α10⋯00α2⋯0⋮⋮⋱⋮00⋯αM]M×M\begin{aligned} \mathbf{Z_1} &= \begin{bmatrix} 1 & 1 & \cdots & 1\\ z_1 & z_2 & \cdots & z_M\\ \vdots & \vdots & \ddots & \vdots\\ z^{(N-L-1)}_1 & z^{(N-L-1)}_2 & \cdots & z^{(N-L-1)}_M \end{bmatrix}_{(N-L) \times M}\\ \mathbf{Z_2} &= \begin{bmatrix} 1 & z_1 & \cdots & z^{L-1}_1\\ 1 & z_2 & \cdots & z^{L-1}_2\\ \vdots & \vdots & \ddots & \vdots\\ 1 & z_M & \cdots & z^{L-1}_M \end{bmatrix}_{M \times L}\\ \mathbf{\Phi} &= \begin{bmatrix} z_1 & 0 & \cdots & 0\\ 0 & z_2 & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & z_M \end{bmatrix}_{M \times M}\\ \mathbf{A} &= \begin{bmatrix} \alpha_1 & 0 & \cdots & 0\\ 0 & \alpha_2 & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & \alpha_M \end{bmatrix}_{M \times M} \end{aligned} Z1Z2ΦA=⎣⎡1z1⋮z1(N−L−1)1z2⋮z2(N−L−1)⋯⋯⋱⋯1zM⋮zM(N−L−1)⎦⎤(N−L)×M=⎣⎡11⋮1z1z2⋮zM⋯⋯⋱⋯z1L−1z2L−1⋮zML−1⎦⎤M×L=⎣⎡z10⋮00z2⋮0⋯⋯⋱⋯00⋮zM⎦⎤M×M=⎣⎡α10⋮00α2⋮0⋯⋯⋱⋯00⋮αM⎦⎤M×M
在无噪声的情况下，LLL也满足约束条件，我们可以看到：矩阵Z1\mathbf{Z_1}Z1和Z2\mathbf{Z_2}Z2都是满秩矩阵，因此矩阵X0\mathbf{X_0}X0和X1\mathbf{X_1}X1的秩与中间的对角阵Φ\mathbf{\Phi}Φ和AΦ\mathbf{A\Phi}AΦ相同，都是MMM.
考虑矩阵X1−λX0=Z1A[Φ−λI]Z2\mathbf{X_1- \lambda X_0} = \mathbf{Z_1A[\Phi-\lambda I]Z_2}X1−λX0=Z1A[Φ−λI]Z2。注意到中间的Φ−λI\mathbf{\Phi-\lambda I}Φ−λI部分：

当λ≠zm(m∈[1,M])\lambda \neq z_m(m \in [1, M])λ=zm(m∈[1,M])时，Φ−λI\mathbf{\Phi-\lambda I}Φ−λI的秩还是MMM，因此X1−λX0\mathbf{X_1-\lambda X_0}X1−λX0的秩还是MMM；
当λ=zm(m∈[1,M])\lambda =z_m(m \in [1, M])λ=zm(m∈[1,M])时，Φ−λI\mathbf{\Phi-\lambda I}Φ−λI的秩就减小为M−1M-1M−1。此时矩阵不满秩，那么行列式det(X1−λX0)=0det(\mathbf{X_1-\lambda X_0})=0det(X1−λX0)=0恒成立。

回顾到广义特征值的性质，我们可以看到，λ=zm(m∈[1,M])\lambda=z_m(m \in [1, M])λ=zm(m∈[1,M])时，刚好可以作为矩阵X1\mathbf{X_1}X1对X0\mathbf{X_0}X0的广义特征值.

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

算法步骤

给定NNN和MMM，并且选择合适的LLL;
创建矩阵X0\mathbf{X_0}X0和X1\mathbf{X_1}X1；
求矩阵X1\mathbf{X_1}X1对X0\mathbf{X_0}X0的MMM个广义特征值λm(m∈[1,M])\lambda_m(m \in [1, M])λm(m∈[1,M])，这些特征值就是对zm(m∈[1,M])z_m(m \in [1, M])zm(m∈[1,M])的估计；
用下式计算DOA:

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

[1]. Direction of Arrival Estimation

DOA算法3：Matrix Pencil相关推荐

【数理知识】矩阵束 Matrix Pencil
介绍若 A0,A1,⋯,AlA_0, A_1, \cdots, A_lA0,A1,⋯,Al 为 lll 个 n×nn \times nn×n 复矩阵,l≥0,Al≠0l \ge 0, A_l ...
推荐系统算法_CUCS | 推荐系统算法简介Matrix Factorization（矩阵分解）
你是否还沉浸于春节收到各路微信红包的喜悦却发现midterm季悄然来临每当你考完这一周的midterm准备松一口气时却发现下一周又有新的midterm在等着你与其复习到头秃,不如打开Amazon或 ...
DOA算法2：ESPRIT算法
目录算法简述算法要得到什么如何得到矩阵 Φ \Phi Φ ESPRIT算法的步骤参考文献(仅写文章的标题,以做记录) 注:本博文为本人阅读论文.文章后的原创笔记,未经授权不允许任何转载或商用行 ...
读书笔记 | 自动驾驶中的雷达信号处理（第6章到达方向(DOA)估计算法）
本文编辑:调皮哥的小助理 6.1介绍 DOA 估计算法在汽车应用雷达处理中非常重要,它构成了雷达数据立方体(距离.速度和角度)的第三个部分--角度.实际上,会有多个未知数量的源信号同时被接收阵列接收, ...
doa的matlab算法,基于MATLAB的DOA估计算法的二维仿真建模
基本描述: DOA算法使用T形或L形天线接收微波信号,计算到达角,然后对二维坐标建模. 要解决的问题: 实际环境中存在多个相干源(例如多径效应,信号反射等),并给出了最佳DOA估计算法. 功能描述: ...
《3D ISP核心技术算法》系列（1）——ToF测量中的多径现象与消除方法
<3D ISP核心技术算法>系列文章是上海数迹智能科技有限公司结合自身多年来在3D ISP技术方向上的积累,推出的系列型科普文章.系列文章将对影响3D成像质量的主要问题逐一进行讲解,深入浅 ...
KerberosSDR代码笔记(5) 信号处理（采样时间延迟计算、相位差计算的2种方法、MUSIC算法）
https://github.com/rtlsdrblog/kerberossdr/blob/master/_signalProcessing/hydra_signal_processor.py 这部 ...
车载毫米波雷达DOA估计综述
谨以此文献给墨の宝祝天天开心! 说明雷达的波达方向估计(direction of arrival, DOA)足够写一本书了,比如 <阵列信号处理及MATLAB实现>.<空 ...
基于黄金分割的修正Powell算法
1.目标函数 function y= TextF(x) %目标函数 %输入参数x[x1,x2]; %输出参数y y = x(1)*x(1)+2*x(2)*x(2)-4*x(1)-2*x(1)*x(2) ...
阵列信号DOA估计系列(二).导向矢量与空间FFT(附代码)
阵列信号DOA估计系列(二).导向矢量在DOA估计里面,经常会看到导向矢量这个名词,也有的地方叫方向矢量,方向矩阵,基本上都是array steering vector 的翻译. 本文首先对均匀线阵 ...

DOA算法3：Matrix Pencil

目录