3 root-MUSIC 算法

3.1 算法原理

??root-MUSIC 算法是 MUSIC 算法的一種多項(xiàng)式求根形式，其基本思想是 Pisarenko 分解。相比于 MUSIC 算法，root-MUSIC 算法無(wú)須譜峰搜索，降低了復(fù)雜度。
??由前面可知，ULA 的方向矢量 $\mathbf{a}(\theta )$ 如下：
a ( θ ) = [ exp ? ( ? j 2 π d 1 sin ? θ / λ ) ? exp ? ( ? j 2 π d m sin ? θ / λ ) ? exp ? ( ? j 2 π d M sin ? θ / λ ) ] \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*} a(θ)= ?exp(?j2πd1?sinθ/λ)?exp(?j2πdm?sinθ/λ)?exp(?j2πdM?sinθ/λ)? ??
根據(jù) $d_m=(m?1)d$，此時(shí)令 $\omega =?2\pi d\sin \theta /\lambda$，則有：
a ( ω ) = [ 1 exp ? ( j ω ) ? exp ? [ j ( M ? 1 ) ω ] ] \begin{equation*} \mathbf{a}(\mathbf{\omega}) = \begin{bmatrix} 1 \\ \exp(j\omega) \\ \vdots \\ \exp\left[j(M-1)\omega\right] \end{bmatrix} \end{equation*} a(ω)= ?1exp(jω)?exp[j(M?1)ω]? ??
令 $z=e^{j\omega }$，則有：
$$\begin{array}{c}\mathbf{p}(z)=[1,z,?,z^{M?1}{]}^T\end{array}$$
??而根據(jù) $\mathbf{a}(\theta {)}^H{\mathbf{u}}_i=0,i=K+1,?,M$，可以得到：
$$\begin{array}{c}{p}_i(z)={\mathbf{u}}_i^H\mathbf{p}(z)=0,i=K+1,?,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)$ 的根即可得到估計(jì)角的信息。然而，上式并不是 $z$ 的多項(xiàng)式，因?yàn)樗€包含了 $z^{?}$ 的冪次項(xiàng)。由于我們只對(duì)單位圓上的 $z$ 值感興趣，所以可以用 ${\mathbf{p}}^T(z^{?1})$ 代替 ${\mathbf{p}}^H(z)$，這就給出了求根 MUSIC 的多項(xiàng)式：
$$\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$，則有：
z M ? 1 p T ( z ? 1 ) G N = [ g [ 1 , 1 ] z M ? 1 + ? + g [ M ? 1 , 1 ] z + g [ M , 1 ] ? g [ 1 , M ] z M ? 1 + ? + g [ M ? 1 , M ] z + g [ M , M ] ] T \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*} zM?1pT(z?1)GN?= ?g[1,1]?zM?1+?+g[M?1,1]?z+g[M,1]??g[1,M]?zM?1+?+g[M?1,M]?z+g[M,M]?? ?T?
其中 $g_{[i,j]}$ 表示矩陣 ${\mathbf{G}}_N$ 的第 $(i,j)$ 元素。此時(shí)易得：
$$\begin{array}{c}\begin{array}{rl}f(z)& =g_{[1,1]}{z}^{M?1}+g_{[2,1]}{z}^{M?2}+?+g_{[M?1,1]}z+g_{[M,1]}\\ & +g_{[1,2]}{z}^M+g_{[2,2]}{z}^{M?1}+?+g_{[M?1,2]}{z}^2+g_{[M,2]}z\\ & ?\\ & +g_{[1,M]}{z}^{2M?2}+g_{[2,M]}{z}^{2M?3}?+g_{[M?1,M]}{z}^M+g_{[M,M]}{z}^{M?1}\end{array}\end{array}$$
由此可以看出 $f(z)$ 中的階數(shù)為 $2M?2$，因此存在 $M?1$ 對(duì)根，且每對(duì)根是相互共軛的關(guān)系。在理想情況下有 $K$ 個(gè)根 $z_1,?,z_K$ 分布在單位圓上，因此在現(xiàn)實(shí)情況中只需要求得上式的 $K$ 個(gè)接近于單位圓上的根即可。對(duì)于 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 代碼實(shí)現(xiàn)中，需要利用前面步驟求解得到噪聲子空間 ${\mathbf{U}}_N$ 來(lái)構(gòu)造 ${\mathbf{G}}_N$，然后通過(guò) ${\mathbf{G}}_N$ 中的元素來(lái)得到多項(xiàng)式 $f(x)$ 的系數(shù)，最后利用 roots 函數(shù)對(duì)多項(xiàng)式求根，部分代碼實(shí)現(xiàn)如下：

Gn = Un*Un';
% 提取多項(xiàng)式系數(shù)并對(duì)多項(xiàng)式求根
coe = zeros(1, 2*M-1);
for i = -(M-1):(M-1)
    coe(-i+M) = sum(diag(Gn,i));
end
r = roots(coe);                                       % 利用roots函數(shù)求多項(xiàng)式的根

??補(bǔ)充說(shuō)明：root-MUSIC 算法與譜搜索方式的 MUSIC 算法原理是一樣的，只不過(guò)是用關(guān)于 $z$ 的矢量來(lái)代替導(dǎo)向矢量，從而用求根過(guò)程代替搜索過(guò)程。

3.2 算法步驟

??root-MUSIC 算法步驟如下（輸入為陣列接收矩陣 $\mathbf{X}$）：文章來(lái)源地址http://www.zghlxwxcb.cn/news/detail-516792.html

1. 計(jì)算協(xié)方差矩陣 $\mathbf{R}=\frac{1}{T}\mathbf{X}{\mathbf{X}}^H$。
2. 對(duì) $\mathbf{R}$ 進(jìn)行特征值分解，并對(duì)特征值進(jìn)行排序，然后取得 $M?K$ 個(gè)較小特征值對(duì)應(yīng)的特征向量來(lái)組成噪聲子空間 ${\mathbf{U}}_N$。
3. 根據(jù)下式定義多項(xiàng)式 $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}$$
   并且求出多項(xiàng)式 $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 代碼實(shí)現(xiàn)

% root_music.m
clear;
clc;
close all;

%% 參數(shù)設(shè)定
c = 3e8;                                              % 光速
fc = 500e6;                                           % 載波頻率
lambda = c/fc;                                        % 波長(zhǎng)
d = lambda/2;                                         % 陣元間距，可設(shè) 2*d = lambda
twpi = 2.0*pi;                                        % 2pi
derad = pi/180;                                       % 角度轉(zhuǎn)弧度
theta = [-20, 30]*derad;                              % 待估計(jì)角度
idx = 0:1:7; idx = idx';                              % 陣元位置索引

M = length(idx);                                      % 陣元數(shù)
K = length(theta);                                    % 信源數(shù)
T = 512;                                              % 快拍數(shù)
SNR = 0;                                              % 信噪比

%% 信號(hào)模型建立
S = randn(K,T) + 1j*randn(K,T);                       % 復(fù)信號(hào)矩陣S，維度為K*T
% A = exp(-1j*twpi*d*idx*sin(theta)/lambda);          % 方向矢量矩陣A，維度為M*K
A = exp(-1j*pi*idx*sin(theta));                       % 2d = lambda，直接忽略不寫(xiě)
X = A*S;                                              % 接收矩陣X，維度為M*T
X = awgn(X,SNR,'measured');                           % 添加噪聲

%% root-MUSIC 算法
% 計(jì)算協(xié)方差矩陣
R = X*X'/T;
% 特征值分解并取得噪聲子空間
[U,D] = eig(R);                                       % 特征值分解
[D,I] = sort(diag(D));                                % 將特征值排序從小到大
U = fliplr(U(:, I));                                  % 對(duì)應(yīng)特征矢量排序，fliplr 之后，較大特征值對(duì)應(yīng)的特征矢量在前面
Un = U(:, K+1:M);                                     % 噪聲子空間
Gn = Un*Un';
% 提取多項(xiàng)式系數(shù)并對(duì)多項(xiàng)式求根
coe = zeros(1, 2*M-1);
for i = -(M-1):(M-1)
    coe(-i+M) = sum(diag(Gn,i));
end
r = roots(coe);                                       % 利用roots函數(shù)求多項(xiàng)式的根
r = r(abs(r)<1);                                      % 找出在單位圓里的根
[lamda,I] = sort(abs(abs(r)-1));                      % 挑選出最接近單位圓的K個(gè)根
Theta = r(I(1:K));
Theta = asin(-angle(Theta)/pi)/derad;                 % 計(jì)算信號(hào)到達(dá)方向角
Theta = sort(Theta).';
disp('估計(jì)結(jié)果：');
disp(Theta);

3.4 參考內(nèi)容

1. 張小飛,陳華偉,仇小鋒.陣列信號(hào)處理及MATLAB實(shí)現(xiàn)[M].電子工業(yè)出版社,2015.
2. 王永良.空間譜估計(jì)理論與算法[M].清華大學(xué)出版社,2004.
3. 【CSDN】root-MUSIC算法
4. 【CSDN】DOA算法1：MUSIC算法（二）

到了這里，關(guān)于【學(xué)習(xí)筆記】【DOA子空間算法】3 root-MUSIC 算法的文章就介紹完了。如果您還想了解更多內(nèi)容，請(qǐng)?jiān)谟疑辖撬阉鱐OY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！