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

这篇具有很好参考价值的文章主要介绍了【学习笔记】【DOA子空间算法】3 root-MUSIC 算法。希望对大家有所帮助。如果存在错误或未考虑完全的地方，请大家不吝赐教，您也可以点击"举报违法"按钮提交疑问。

3 root-MUSIC 算法

3.1 算法原理

root-MUSIC 算法是 MUSIC 算法的一种多项式求根形式，其基本思想是 Pisarenko 分解。相比于 MUSIC 算法，root-MUSIC 算法无须谱峰搜索，降低了复杂度。
由前面可知，ULA 的方向矢量 $\mathbf{a}(\theta)$ 如下：
$\begin{equation*} \mathbf{a}(\mathbf{\theta}) = \begin{bmatrix} \exp(-j2\pi d_1\sin\theta/\lambda) \\ \vdots \\ \exp(-j2\pi d_m\sin\theta/\lambda) \\ \vdots \\ \exp(-j2\pi d_M\sin\theta/\lambda) \end{bmatrix} \end{equation*}$
根据 $d_m = (m-1)d$ ，此时令 $\omega = -2\pi d \sin\theta/\lambda$ ，则有：
$\begin{equation*} \mathbf{a}(\mathbf{\omega}) = \begin{bmatrix} 1 \\ \exp(j\omega) \\ \vdots \\ \exp\left[j(M-1)\omega\right] \end{bmatrix} \end{equation*}$
令 $e^{j\omega}$ ，则有：
$\begin{equation*} \mathbf{p}(z) = [1, z, \cdots, z^{M-1}]^T \end{equation*}$
而根据 $\mathbf{a}(\theta)^H\mathbf{u}_i = 0, i = K+1, \cdots, M$ ，可以得到：
$\begin{equation*} p_i(z) = \mathbf{u}_i^H\mathbf{p}(z) = 0, i = K+1, \cdots, M \end{equation*}$
综合以上，可得到：
$\begin{equation*} f(z) = \mathbf{p}^H(z)\mathbf{U}_N\mathbf{U}_N^H \mathbf{p}(z) \end{equation*}$
求得 $f (z)$ 的根即可得到估计角的信息。然而，上式并不是 $z$ 的多项式，因为它还包含了 $z^*$ 的幂次项。由于我们只对单位圆上的 $z$ 值感兴趣，所以可以用 $\mathbf{p}^{T}(z^{-1})$ 代替 $\mathbf{p}^{H}(z)$ ，这就给出了求根 MUSIC 的多项式：
$\begin{equation*} f(z) = z^{M-1}\mathbf{p}^{T}(z^{-1})\mathbf{U}_N\mathbf{U}_N^H \mathbf{p}(z) \end{equation*}$
令 $\mathbf{G}_N = \mathbf{U}_N\mathbf{U}_N^H$ ，则有：
$\begin{equation*} z^{M-1}\mathbf{p}^{T}(z^{-1})\mathbf{G}_N = \begin{bmatrix} g_{[1,1]}z^{M-1} + \cdots + g_{[M-1,1]} z + g_{[M,1]} \\ \vdots \\ g_{[1,M]}z^{M-1} + \cdots + g_{[M-1,M]} z + g_{[M,M]} \end{bmatrix}^T \end{equation*}$
其中 $g_{[i,j]}$ 表示矩阵 $\mathbf{G}_N$ 的第 $(i, j)$ 元素。此时易得：
$\begin{equation*} \begin{aligned} 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 \\ &\qquad \vdots \\ &+ 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{aligned} \end{equation*}$
由此可以看出 $f (z)$ 中的阶数为 $2 M - 2$ ，因此存在 $M - 1$ 对根，且每对根是相互共轭的关系。在理想情况下有 $K$ 个根 $z_1, \cdots, z_K$ 分布在单位圆上，因此在现实情况中只需要求得上式的 $K$ 个接近于单位圆上的根即可。对于 ULA，有：
$\begin{equation*} \theta_i = \arcsin\left(-\frac{\lambda}{2\pi d} \arg \{ z_i \} \right), i = 1, \cdots, K \end{equation*}$
在 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}$ ）：文章来源地址https://www.toymoban.com/news/detail-516792.html

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

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);