【学习笔记】【DOA子空间算法】3 root-MUSIC 算法

3 root-MUSIC 算法

3.1 算法原理

  root-MUSIC 算法是 MUSIC 算法的一种多项式求根形式，其基本思想是 Pisarenko 分解。相比于 MUSIC 算法，root-MUSIC 算法无须谱峰搜索，降低了复杂度。
  由前面可知，ULA 的方向矢量 $\mathbf{a}(\theta )$ 如下：
$$\begin{array}{c}\mathbf{a}(\theta )=\left[\begin{array}{c}\exp (-j2\pi d_1\sin \theta /\lambda )\\ ⋮\\ \exp (-j2\pi d_m\sin \theta /\lambda )\\ ⋮\\ \exp (-j2\pi d_M\sin \theta /\lambda )\end{array}\right]\end{array}$$
根据 $d_m=(m-1)d$，此时令 $\omega =-2\pi d\sin \theta /\lambda$，则有：
$$\begin{array}{c}\mathbf{a}(\omega )=\left[\begin{array}{c}1\\ \exp (j\omega )\\ ⋮\\ \exp \left[j(M-1)\omega \right]\end{array}\right]\end{array}$$
令 $z=e^{j\omega }$，则有：
$$\begin{array}{c}\mathbf{p}(z)=[1,z,\cdots ,z^{M-1}{]}^T\end{array}$$
  而根据 $\mathbf{a}(\theta {)}^H{\mathbf{u}}_i=0,i=K+1,\cdots ,M$，可以得到：
$$\begin{array}{c}{p}_i(z)={\mathbf{u}}_i^H\mathbf{p}(z)=0,i=K+1,\cdots ,M\end{array}$$
  综合以上，可得到：
$$\begin{array}{c}f(z)={\mathbf{p}}^H(z){\mathbf{U}}_N{\mathbf{U}}_N^H\mathbf{p}(z)\end{array}$$
求得 $f(z)$ 的根即可得到估计角的信息。然而，上式并不是 $z$ 的多项式，因为它还包含了 $z^{\ast }$ 的幂次项。由于我们只对单位圆上的 $z$ 值感兴趣，所以可以用 ${\mathbf{p}}^T(z^{-1})$ 代替 ${\mathbf{p}}^H(z)$，这就给出了求根 MUSIC 的多项式：
$$\begin{array}{c}f(z)=z^{M-1}{\mathbf{p}}^T(z^{-1}){\mathbf{U}}_N{\mathbf{U}}_N^H\mathbf{p}(z)\end{array}$$
令 ${\mathbf{G}}_N={\mathbf{U}}_N{\mathbf{U}}_N^H$，则有：
$$\begin{array}{c}{z}^{M-1}{\mathbf{p}}^T(z^{-1}){\mathbf{G}}_N={\left[\begin{array}{c}{g}_{[1,1]}{z}^{M-1}+\cdots +g_{[M-1,1]}z+g_{[M,1]}\\ ⋮\\ g_{[1,M]}{z}^{M-1}+\cdots +g_{[M-1,M]}z+g_{[M,M]}\end{array}\right]}^T\end{array}$$
其中 $g_{[i,j]}$ 表示矩阵 ${\mathbf{G}}_N$ 的第 $(i,j)$ 元素。此时易得：
$$\begin{array}{c}\begin{array}{rl}f(z)& =g_{[1,1]}{z}^{M-1}+g_{[2,1]}{z}^{M-2}+\cdots +g_{[M-1,1]}z+g_{[M,1]}\\ & +g_{[1,2]}{z}^M+g_{[2,2]}{z}^{M-1}+\cdots +g_{[M-1,2]}{z}^2+g_{[M,2]}z\\ & ⋮\\ & +g_{[1,M]}{z}^{2M-2}+g_{[2,M]}{z}^{2M-3}\cdots +g_{[M-1,M]}{z}^M+g_{[M,M]}{z}^{M-1}\end{array}\end{array}$$
由此可以看出 $f(z)$ 中的阶数为 $2M-2$，因此存在 $M-1$ 对根，且每对根是相互共轭的关系。在理想情况下有 $K$ 个根 $z_1,\cdots ,z_K$ 分布在单位圆上，因此在现实情况中只需要求得上式的 $K$ 个接近于单位圆上的根即可。对于 ULA，有：
θ i = arcsin ⁡ ( − λ 2 π d arg ⁡ { z i } ) , i = 1 , ⋯   , K \begin{equation*} \theta_i = \arcsin\left(-\frac{\lambda}{2\pi d} \arg \{ z_i \} \right), i = 1, \cdots, K \end{equation*} θi​=arcsin(−2πdλ​arg{zi​}),i=1,⋯,K​
  在 matlab 代码实现中，需要利用前面步骤求解得到噪声子空间 ${\mathbf{U}}_N$ 来构造 ${\mathbf{G}}_N$，然后通过 ${\mathbf{G}}_N$ 中的元素来得到多项式 $f(x)$ 的系数，最后利用 roots 函数对多项式求根，部分代码实现如下：

Gn = Un*Un';
% 提取多项式系数并对多项式求根
coe = zeros(1, 2*M-1);
for i = -(M-1):(M-1)
    coe(-i+M) = sum(diag(Gn,i));
end
r = roots(coe);                                       % 利用roots函数求多项式的根

  补充说明：root-MUSIC 算法与谱搜索方式的 MUSIC 算法原理是一样的，只不过是用关于 $z$ 的矢量来代替导向矢量，从而用求根过程代替搜索过程。

3.2 算法步骤

  root-MUSIC 算法步骤如下（输入为阵列接收矩阵 $\mathbf{X}$）：

1. 计算协方差矩阵 $\mathbf{R}=\frac{1}{T}\mathbf{X}{\mathbf{X}}^H$。
2. 对 $\mathbf{R}$ 进行特征值分解，并对特征值进行排序，然后取得 $M-K$ 个较小特征值对应的特征向量来组成噪声子空间 ${\mathbf{U}}_N$。
3. 根据下式定义多项式 $f(z)$：
   $$\begin{array}{c}f(z)=z^{M-1}{\mathbf{p}}^T(z^{-1}){\mathbf{U}}_N{\mathbf{U}}_N^H\mathbf{p}(z)\end{array}$$
   并且求出多项式 $f(z)$ 的根 { z i : i = 1 , ⋯   , K } \{z_i:i = 1,\cdots, K\} {zi​:i=1,⋯,K}。
4. 利用下式求得角度 { θ i : i = 1 , ⋯   , K } \{\theta_i: i = 1,\cdots, K\} {θi​:i=1,⋯,K}：
   θ i = arcsin ⁡ ( − λ 2 π d arg ⁡ { z i } ) , i = 1 , ⋯   , K \begin{equation*} \theta_i = \arcsin\left(-\frac{\lambda}{2\pi d} \arg \{ z_i \} \right), i = 1, \cdots, K \end{equation*} θi​=arcsin(−2πdλ​arg{zi​}),i=1,⋯,K​

3.3 代码实现

% root_music.m
clear;
clc;
close all;

%% 参数设定
c = 3e8;                                              % 光速
fc = 500e6;                                           % 载波频率
lambda = c/fc;                                        % 波长
d = lambda/2;                                         % 阵元间距，可设 2*d = lambda
twpi = 2.0*pi;                                        % 2pi
derad = pi/180;                                       % 角度转弧度
theta = [-20, 30]*derad;                              % 待估计角度
idx = 0:1:7; idx = idx';                              % 阵元位置索引

M = length(idx);                                      % 阵元数
K = length(theta);                                    % 信源数
T = 512;                                              % 快拍数
SNR = 0;                                              % 信噪比

%% 信号模型建立
S = randn(K,T) + 1j*randn(K,T);                       % 复信号矩阵S，维度为K*T
% A = exp(-1j*twpi*d*idx*sin(theta)/lambda);          % 方向矢量矩阵A，维度为M*K
A = exp(-1j*pi*idx*sin(theta));                       % 2d = lambda，直接忽略不写
X = A*S;                                              % 接收矩阵X，维度为M*T
X = awgn(X,SNR,'measured');                           % 添加噪声

%% root-MUSIC 算法
% 计算协方差矩阵
R = X*X'/T;
% 特征值分解并取得噪声子空间
[U,D] = eig(R);                                       % 特征值分解
[D,I] = sort(diag(D));                                % 将特征值排序从小到大
U = fliplr(U(:, I));                                  % 对应特征矢量排序，fliplr 之后，较大特征值对应的特征矢量在前面
Un = U(:, K+1:M);                                     % 噪声子空间
Gn = Un*Un';
% 提取多项式系数并对多项式求根
coe = zeros(1, 2*M-1);
for i = -(M-1):(M-1)
    coe(-i+M) = sum(diag(Gn,i));
end
r = roots(coe);                                       % 利用roots函数求多项式的根
r = r(abs(r)<1);                                      % 找出在单位圆里的根
[lamda,I] = sort(abs(abs(r)-1));                      % 挑选出最接近单位圆的K个根
Theta = r(I(1:K));
Theta = asin(-angle(Theta)/pi)/derad;                 % 计算信号到达方向角
Theta = sort(Theta).';
disp('估计结果：');
disp(Theta);

3.4 参考内容

1. 张小飞,陈华伟,仇小锋.阵列信号处理及MATLAB实现[M].电子工业出版社,2015.
2. 王永良.空间谱估计理论与算法[M].清华大学出版社,2004.
3. 【CSDN】root-MUSIC算法
4. 【CSDN】DOA算法1：MUSIC算法（二）