root MUSIC 算法补充说明

  这篇笔记是上一篇关于 root MUSIC 笔记的补充。

多项式求根

  要理解 root MUSIC 算法，需要理解多项式求根的相关知识。给定多项式 $P(x)$：
$$P(x)=a_0+a_{1}x+\cdots +a_n{x}^n$$
容易看出 $P(x)$ 中只有一个未知数 $x$，且未知数的最高次数为 $n$，因此称 $P(x)$ 为一元 $n$ 次多项式，同时系数 { a i ∈ C : i = 0 , ⋯   , n } \{a_i\in\mathbb{C}:i = 0,\cdots, n\} {ai​∈C:i=0,⋯,n}。而多项式求根就是求得一元 $n$ 次方程式 $P(x)=0$ 的解，这个解被称作根或者零点。
  在进行后续的讨论前，还需要清楚，根据代数基本定理，$n$ 次复系数多项式方程在复数域内有且只有 $n$ 个根（这里的重根按重数计算）。

root MUSIC 算法原理

  root MUSIC 算法是 MUSIC 算法的一种多项式求根形式。回忆一下，传统 MUSIC 算法利用了噪声子空间矩阵 ${\mathbf{U}}_n$ 和搜索方向矢量 $\mathbf{a}(\theta )$ 来构造空间谱：
$$\begin{gathered} P(\theta )=\frac{1}{{\mathbf{a}}^H(\theta ){\mathbf{U}}_n{\mathbf{U}}_n^H\mathbf{a}(\theta )} \\ \mathbf{a}(\theta )={\left[1,e^{-\mathrm{j}2\pi d\sin \theta /\lambda },\cdots ,e^{-\mathrm{j}2\pi (M-1)d\sin \theta /\lambda }\right]}^T \end{gathered}$$
在 { θ = θ k : k = 1 , ⋯   , K } \{\theta = \theta_k:k = 1,\cdots,K\} {θ=θk​:k=1,⋯,K} 时 $P(\theta )$ 将产生峰值，换句话说此时 $P^{-1}(\theta )=0$。
  在接下来的讨论中，我们令 $P^{-1}(\theta )={\mathbf{a}}^H(\theta )\mathbf{Ga}(\theta )$，此时我们可以知道，MUSIC 算法满足 $\mathbf{G}\triangleq {\mathbf{U}}_n{\mathbf{U}}_n^H$，而 Capon 算法满足 $\mathbf{G}\triangleq {\mathbf{R}}^{-1}$。需要注意的是无论是 MUSIC 算法还是 Capon 算法，$\mathbf{G}$ 均是 Hermitian 矩阵。
  令 $\omega =-2\pi d\sin \theta /\lambda$ 以及 $z=e^{\mathrm{j}\omega }$，我们将会得到：
$$\begin{array}{rl}\mathbf{a}(z)& =[1,z,z^2,\cdots ,z^{M-1}{]}^T=\mathbf{a}(\theta )\\ P^{-1}(z)& ={\mathbf{a}}^H(z)\mathbf{Ga}(z)=P^{-1}(\theta )\end{array}$$
接下来我们展开 $P^{-1}(z)$：
$$\begin{array}{rl}{P}^{-1}(z)& ={\mathbf{a}}^H(z)\mathbf{Ga}(z)\\ & =[1,z^{\ast },(z^{\ast }{)}^2,\cdots ,(z^{\ast }{)}^{M-1}]\mathbf{G}[1,z,z^2,\cdots ,z^{M-1}{]}^T\\ & =[1,z^{-1},z^{-2},\cdots ,z^{-M+1}]\mathbf{G}[1,z,z^2,\cdots ,z^{M-1}{]}^T\\ & =\sum _{m=0}^{M-1}\sum _{n=0}^{M-1}{z}^{-m}{\mathbf{G}}_{[m,n]}{z}^n\\ & =\sum _{m=0}^{M-1}\sum _{n=0}^{M-1}{z}^{n-m}{\mathbf{G}}_{[m,n]}\\ & =\sum _{p=-M+1}^{M-1}{a}_p{z}^{-p}\end{array}$$
其中 ${\mathbf{G}}_{[m,n]}$ 表示矩阵 $\mathbf{G}$ 的第 $m$ 行第 $n$ 列元素，$a_p$ 表示矩阵 $\mathbf{G}$ 的第 $p$ 条对角线的求和：
$$a_p\triangleq \sum _{m-n=p}{\mathbf{G}}_{[m,n]}$$
  到这里我们已经可以看出，传统 MUSIC 算法对 $P(\theta )$ 求峰值，其实等价于对 $P^{-1}(z)$ 求根，为了方便大家的理解，我们令 $M=3$，此时会得到一条简单的式子：
$$P^{-1}(z)=a_2{z}^{-2}+a_1{z}^{-1}+a_0{z}^0+a_{-1}{z}^1+a_{-2}{z}^2$$
可以看出，其实 $P^{-1}(\theta )$ 是一个 $2M-1=5$ 项的多项式，但还存在一个问题，$P^{-1}(\theta )$ 中存在负整数次数，我们令 $P^{-1}(z)z^{M-1}$ 将负整数次数消除即可，操作前后，求根的结果是一样的，因此我们可以说 $P^{-1}(z)z^{M-1}$ 是一个一元的 $2M-1$ 项的 $2M-2$ 次的多项式。更进一步地，我们可以说求解 $P^{-1}(z)z^{M-1}=0$ 将会得到 $2M-2$ 个根，从已知条件我们知道，其中 $K$ 个根必定是 $e^{\mathrm{j}{\omega }_k}$（$e^{\mathrm{j}{\omega }_k}$ 的幅值是 $1$，因此该 $K$ 点在单位圆上），在这里 ${\omega }_k=-2\pi d\sin {\theta }_k/\lambda$。
  总结一下，MUSIC 算法的谱峰搜索操作等价于对方程式 $P^{-1}(z)z^{M-1}=0$ 求根，root MUSIC 算法所做的，就是利用 $\mathbf{G}$ 的多条对角线求和得到对应的多项式系数，从而求解得 $2M-2$ 个根，接着筛选得到合适的 $K$ 个根 $z_k$，再通过 $z_k$ 推导得到原先的 ${\theta }_k$。

如何从 $2M-2$ 个根中确定 $K$ 个根

  那么如何从 $2M-2$ 个根中确定 $K$ 个根？这个问题大部分的论文和博客都一笔带过了。从前面的讨论可知，多项式系数是由 $\mathbf{G}$ 的多条对角线求和得到，同时 $\mathbf{G}$ 是 Hermitian 矩阵，因此以下式子可以得到：
$$a_p=a_{-p}^{\ast }$$
这个等式意味着在 $2M-1$ 个系数 { a p : p = − M + 1 , ⋯   , M − 1 } \{a_p:p=-M+1,\cdots,M-1\} {ap​:p=−M+1,⋯,M−1} 中，前 $M-1$ 个和后 $M-1$ 个系数是前后共轭对称，同时正中间的系数是实数。
  我们继续假设 $M=3$，$P^{-1}(z)=0$ 可以进一步表示如下：
$$P^{-1}(z)=a_2{z}^{-2}+a_1{z}^{-1}+a_0+a_1^{\ast }{z}^1+a_2^{\ast }{z}^2=0$$
  如此我们分析 $P^{-1}(1/z^{\ast })$，可得：
$$\begin{array}{rl}{P}^{-1}(1/z^{\ast })& =a_2(z^{\ast }{)}^2+a_1(z^{\ast }{)}^1+a_0+a_1^{\ast }(z^{\ast }{)}^{-1}+a_2^{\ast }(z^{\ast }{)}^{-2}\\ & =[P^{-1}(z){]}^{\ast }=P^{-1}(z)=0\end{array}$$
这意味着假若 $z_1=\rho e^{\mathrm{j}\phi }$ 是 $P^{-1}(z)=0$ 的根，那么 $z_2=1/z_1^{\ast }=1/\rho e^{\mathrm{j}\phi }$ 同样是 $P^{-1}(z)=0$ 的根。观察 $z_1$ 和 $z_2$ 在复平面的位置，将会观察得到 $z_1$ 和 $z_2$ 是关于单位圆有一个类似对称的关系；简单来说，这个现象是因为 $z_1$ 和 $z_2$ 是幅值互为倒数而相位相等的关系，因此它们就像是关于 $e^{\mathrm{j}\phi }$ 对称一样（$e^{\mathrm{j}\phi }$ 的幅值是 $1$，因此该点在单位圆上）。
  综上所述，通过 $P^{-1}(z)$ 得到 $2M-2$ 个根，它们是关于单位圆对称的 $M-1$ 对根，因此一定有 $K$ 对根在单位圆附近，所以我们只需要从 $2M-2$ 个根中找 $M-1$ 个处于单位圆内的根（找 $M-1$ 个处于单位圆外的根同样是可以的，因为角度信息其实只存在于 $z_k$ 的相位中，与幅值无关），最后确定最接近单位圆的 $K$ 个根就可以确定 $z_k$。

从复数域上观察 $2M-2$ 个根的分布

  我们从实验中进一步观察 $2M-2$ 个根的分布，matlab 代码实现如下：

clear; close all; clc;

%% Parameters
lambda     = 3e8/1e9;         % wavelength, c/f
d          = lambda/4;        % distance between sensors
theta      = [10,20];         % true DoAs, 1 times K vector
theta      = sort(theta);
M          = 16;              % # of sensors
T          = 500;             % # of snapshots
K          = length(theta);   % # of signals
noise_flag = 1;
SNR        = 0;               % signal-to-noise ratio

%% Signals
S = exp(1j*2*pi*randn(K,T)); % signal matrix
A = exp(-1j*(0:M-1)'*2*pi*d/lambda*sind(theta)); % steering vector matrix
N = noise_flag.*sqrt(10.^(-SNR/10))*(1/sqrt(2))*(randn(M,T)+1j*randn(M,T)); % noise matrix
X = A*S+N; % received matrix
R = X*X'/T; % covariance matrix

%% DoA:root-MUSIC
[U,~] = svd(R); % SVD
Un = U(:, K+1:end); % noise subspace matrix
Gn = Un*Un';
coe = arrayfun(@(i) sum(diag(Gn, M-i)),(1:2*M-1));
r = roots(coe); % 2M-2 roots

%% plot
dis = sort(abs(r)-1);
disp(dis);
cnt = sum(dis<0);
disp(cnt); % 记录单位圆内的根个数

% 提取实部和虚部
realPart = real(r);
imaginaryPart = imag(r);

% 绘制复平面
figure;
scatter(realPart, imaginaryPart, 'filled');
hold on;

% 绘制单位圆
theta = linspace(0, 2*pi, 100);
unitCircleReal = cos(theta);
unitCircleImag = sin(theta);
plot(unitCircleReal, unitCircleImag, 'r--', 'LineWidth', 1);

xlabel('实部');
ylabel('虚部');
title('复平面上的复数点和单位圆');
grid on;
box on;

%% find K roots
r = r(abs(r)<1);
[~, idx] = sort(abs(abs(r)-1));
z = angle(r(idx));
theta = sort(asin(-z(1:K)/2/pi/d*lambda)/pi*180).';

  设 $M=16$ 和 $K=2$，根的分布如下图所示，可以看到 $2M-2=30$ 个根，其中 $2K=4$ 个接近单位圆的根对应着估计角度：