unitary MUSIC 算法

??論文 A Unitary Transformation Method for Angle-of-Arrival Estimation 中提出了 unitary MUSIC 的算法，直譯就是酉 MUSIC 算法，即酉變換 MUSIC 算法。該算法的目的是簡(jiǎn)化計(jì)算復(fù)雜度，將傳統(tǒng) MUSIC 算法中的復(fù)數(shù) SVD 和復(fù)數(shù)網(wǎng)格搜索計(jì)算轉(zhuǎn)化為實(shí)數(shù)計(jì)算。在學(xué)習(xí) unitary MUSIC 之前需要理解 Hermitian 矩陣及 Persymmetric 矩陣的概念及性質(zhì)：

1. Hermitian 矩陣指的是滿足 ${\mathbf{A}}^H=\mathbf{A}$ 的矩陣 $\mathbf{A}$；
2. Persymmetric 矩陣指的是滿足 $\mathbf{AJ}=\mathbf{J}{\mathbf{A}}^T$ 的矩陣 $\mathbf{A}$，其中 $\mathbf{J}$ 其對(duì)角線從左下至右上，很多地方又稱 $\mathbf{J}$ 為選擇矩陣。
3. 假若矩陣 $\mathbf{A}$ 既為 Hermitian 矩陣又為 Persymmetric 矩陣，則滿足：
   $$\mathbf{J}{\mathbf{A}}^{?}\mathbf{J}=\mathbf{A}$$
   其中 ${\mathbf{A}}^{?}$ 為 $\mathbf{A}$ 的共軛。

??在接下來(lái)的討論中，$\mathbf{I}$ 和 $\mathbf{J}$ 分別用作表示單位矩陣和選擇矩陣，下文中將會(huì)出現(xiàn)這兩種矩陣的運(yùn)算，例如 $\mathbf{AI}$ 或 $\mathbf{JB}$，設(shè) $\mathbf{A}$ 和 $\mathbf{B}$ 均為方陣，如果沒(méi)有特別強(qiáng)調(diào)，則說(shuō)明 $\mathbf{I}$ 和 $\mathbf{J}$ 分別和 $\mathbf{A}$ 和 $\mathbf{B}$ 同維度。

算法原理

??前面討論的子空間算法中，復(fù)協(xié)方差矩陣的特征值分解是至關(guān)重要的一步，然而該步的計(jì)算量很高。為了降低計(jì)算量，unitary MUSIC 算法考慮利用一個(gè)酉矩陣將原先的復(fù)協(xié)方差矩陣 $\mathbf{R}$ 轉(zhuǎn)換成實(shí)協(xié)方差矩陣，同時(shí)傳統(tǒng)算法中的復(fù)空間搜索向量 $\mathbf{a}(\theta )$ 也用實(shí)向量來(lái)代替。
??unitary MUSIC 算法的提出基于一個(gè)性質(zhì)，即若不相關(guān)的窄帶遠(yuǎn)場(chǎng)信號(hào)源射入均勻線陣中，其協(xié)方差矩陣不僅是 Hermitian，且 Persymmetric。通常估計(jì)的協(xié)方差矩陣 $\mathbf{R}?\hat{\mathbf{R}}$ 僅僅只是 Hermitian 矩陣但不滿足 Persymmetric 性質(zhì)，需要先獲得一個(gè)滿足 Persymmetric 性質(zhì)的估計(jì)協(xié)方差矩陣：
$$\begin{array}{c}\mathbf{R}?\frac{1}{2}(\hat{\mathbf{R}}+\mathbf{J}{\hat{\mathbf{R}}}^{?}\mathbf{J})\end{array}$$
??假設(shè)陣元數(shù) $M$ 為偶數(shù)，unitary MUSIC 算法引入了一個(gè)酉矩陣 $\mathbf{Q}\in {\mathbb{C}}^{M\times M}$：
$$\mathbf{Q}=\frac{1}{\sqrt{2}}\left[\begin{array}{cc}\mathbf{I}& \mathbf{J}\\ \mathrm{j}\mathbf{J}& ?\mathrm{j}\mathbf{I}\end{array}\right]$$
其中 $\mathbf{I}$ 和 $\mathbf{J}$ 分別為單位矩陣和選擇矩陣，且該兩個(gè)矩陣維度均為 $\frac{M}{2}\times \frac{M}{2}$。易得 $\mathbf{Q}$ 為酉矩陣，即 ${\mathbf{Q}}^{?1}={\mathbf{Q}}^H$，同時(shí)滿足：
$${\mathbf{Q}}^{?}\mathbf{J}=\mathbf{Q}$$
??至此到了本算法的關(guān)鍵，它在于證明由于 $\mathbf{R}$ 是 Hermitian 且 Persymmetric 矩陣，$\mathbf{QR}{\mathbf{Q}}^H$ 是實(shí)對(duì)稱矩陣：

因?yàn)?$\mathbf{R}$ 為 Hermitian，易得 $\mathbf{QR}{\mathbf{Q}}^H$ 為 Hermitian；因此只需要證明 $\mathbf{QR}{\mathbf{Q}}^H$ 是實(shí)矩陣，即證明 $(\mathbf{QR}{\mathbf{Q}}^H{)}^{?}=\mathbf{QR}{\mathbf{Q}}^H$：
$$\begin{array}{c}\begin{array}{rl}& (\mathbf{QR}{\mathbf{Q}}^H{)}^{?}\\ =& {\mathbf{Q}}^{?}{\mathbf{R}}^{?}{\mathbf{Q}}^T\\ =& ({\mathbf{Q}}^{?}\mathbf{J})(\mathbf{J}{\mathbf{R}}^{?}\mathbf{J})(\mathbf{J}{\mathbf{Q}}^T)\\ =& \mathbf{QR}{\mathbf{Q}}^H\end{array}\end{array}$$
由此得證。

??綜上所述，unitary MUSIC 算法引入酉矩陣 $\mathbf{Q}$ 并令 $\mathbf{R}?\mathbf{QR}{\mathbf{Q}}^H$，使得酉變換后的協(xié)方差矩陣變?yōu)閷?shí)對(duì)稱矩陣，接著對(duì)其特征值分解即可進(jìn)行后續(xù)的搜索步驟。而 ULA 的搜索方向矢量為：
$$\mathbf{a}(\theta )={\left[1,e^{?\mathrm{j}2\pi d\sin \theta /\lambda },?,e^{?\mathrm{j}(M?1)2\pi d\sin \theta /\lambda }\right]}^T$$
則 unitary MUSIC 算法的搜索方向矢量為 $\mathbf{a}(\theta )?\mathbf{Qa}(\theta )$。
??為了進(jìn)一步降低算法計(jì)算復(fù)雜度，Unitary MUSIC 算法考慮將搜索方向矢量也用實(shí)變量代替，做法如下：
$$\begin{array}{rl}\mathbf{a}(\theta )& ?e^{j\frac{M?1}{2}2\pi d\sin \theta /\lambda }\mathbf{Qa}(\theta )\\ & =\mathbf{Q}{\left[e^{\mathrm{j}\frac{M?1}{2}2\pi d\sin \theta /\lambda },?,e^{\mathrm{j}\frac{1}{2}2\pi d\sin \theta /\lambda },e^{?\mathrm{j}\frac{1}{2}2\pi d\sin \theta /\lambda },?,e^{?\mathrm{j}\frac{M?1}{2}2\pi d\sin \theta /\lambda }\right]}^T\end{array}$$
不難看出原方向矢量對(duì)應(yīng)的陣列索引位置為 { 0 , 1 , ? ? , M ? 1 } \{0,1,\cdots,M-1\} {0,1,?,M?1}，而更新后方向矢量對(duì)應(yīng)的陣列索引位置為 { ? M ? 1 2 , ? M ? 3 2 , ? ? , ? 1 2 , 1 2 , ? ? , M ? 3 2 , M ? 1 2 } \{-\frac{M-1}{2},-\frac{M-3}{2},\cdots,-\frac{1}{2},\frac{1}{2},\cdots,\frac{M-3}{2},\frac{M-1}{2}\} {?2M?1?,?2M?3?,?,?21?,21?,?,2M?3?,2M?1?}。此時(shí) unitary MUSIC 算法的搜索方向矢量更新為：
a ￣ ( θ ) = e j M ? 1 2 2 π d sin ? θ / λ Q a ( θ ) = 2 [ cos ? ( M ? 1 2 2 π d sin ? θ / λ ) ? cos ? ( 1 2 2 π d sin ? θ / λ ) sin ? ( ? 1 2 2 π d sin ? θ / λ ) ? sin ? ( ? M ? 1 2 2 π d sin ? θ / λ ) ] \begin{aligned} \overline{\mathbf{a}}(\theta)&= e^{\mathrm{j}\frac{M-1}{2}2\pi d\sin\theta/\lambda}\mathbf{Q}\mathbf{a}(\theta) \\ &= \sqrt{2} \begin{bmatrix} \cos\left(\frac{M-1}{2}2\pi d\sin\theta/\lambda\right) \\ \vdots \\ \cos\left(\frac{1}{2}2\pi d\sin\theta/\lambda\right) \\ \sin\left(-\frac{1}{2}2\pi d\sin\theta/\lambda\right) \\ \vdots \\ \sin\left(-\frac{M-1}{2}2\pi d\sin\theta/\lambda\right) \end{bmatrix} \end{aligned} a(θ)?=ej2M?1?2πdsinθ/λQa(θ)=2 ? ?cos(2M?1?2πdsinθ/λ)?cos(21?2πdsinθ/λ)sin(?21?2πdsinθ/λ)?sin(?2M?1?2πdsinθ/λ)? ??
??當(dāng) $M$ 為奇數(shù)時(shí)，酉矩陣 $\mathbf{Q}$ 的形式為：
Q = 1 2 [ I O J O T 2 O T j J O ? j I ] \mathbf{Q} = \frac{1}{\sqrt{2}} \begin{bmatrix} \mathbf{I} & \mathbf{O} & \mathbf{J} \\ \mathbf{O}^T & \sqrt{2} & \mathbf{O}^T \\ \mathrm{j}\mathbf{J} & \mathbf{O} & -\mathrm{j}\mathbf{I} \end{bmatrix} Q=2 ?1? ?IOTjJ?O2 ?O?JOT?jI? ?
其中 $\mathbf{O}$ 為零矩陣。

進(jìn)一步理解

??在上一小節(jié)中，我整理了論文對(duì)于 unitary MUSIC 算法的解釋及證明，其核心思想在于酉矩陣 $\mathbf{Q}$ 的提出。在本小節(jié)中，我將進(jìn)一步對(duì) $\mathbf{Q}$ 的作用談?wù)勛约旱睦斫狻?br> ??在論文中，unitary 算法一直強(qiáng)調(diào) $\mathbf{R}$ 的 Hermitian 及 Persymmetric 性質(zhì)，因?yàn)榧偃?$\mathbf{R}$ 不滿足這兩個(gè)性質(zhì)，$\mathbf{QR}{\mathbf{Q}}^H$ 將不是實(shí)矩陣，但是利用 $\mathrm{Re}(\mathbf{QR}{\mathbf{Q}}^H)$ 來(lái)進(jìn)行后續(xù)的估計(jì)其實(shí)是可以估計(jì)角度的。假設(shè) $M=4$ 和 $\phi =2\pi d\sin \theta /\lambda$，展開(kāi) $\mathbf{Qa}(\theta )$：
Q a ( θ ) = 2 [ 1 0 0 1 0 1 1 0 0 j ? j 0 j 0 0 ? j ] [ 1 e ? j φ e ? j 2 φ e ? j 3 φ ] = 2 [ 1 + e ? j 3 φ e ? j φ + e ? j 2 φ j ( e ? j φ ? e ? j 2 φ ) j ( 1 ? e ? j 3 φ ) ] = 2 e ? j 3 2 φ [ e j 3 2 φ + e ? j 3 2 φ e j 1 2 φ + e ? j 1 2 φ j ( e j 1 2 φ ? e ? j 1 2 φ ) j ( e j 3 2 φ ? e ? j 3 2 φ ) ] = 2 e ? j 3 2 φ [ cos ? ( 3 2 φ ) cos ? ( 1 2 φ ) sin ? ( ? 1 2 φ ) sin ? ( ? 3 2 φ ) ] = e ? j 3 2 φ a ￣ ( θ ) \begin{aligned} \mathbf{Q}\mathbf{a}(\theta) &= \sqrt{2} \begin{bmatrix} 1 & 0 & 0 & 1 \\ 0 & 1 & 1 & 0 \\ 0 & \mathrm{j} & -\mathrm{j} & 0 \\ \mathrm{j} & 0 & 0 & -\mathrm{j} \end{bmatrix} \begin{bmatrix} 1\\ e^{-\mathrm{j}\varphi}\\ e^{-\mathrm{j}2\varphi}\\ e^{-\mathrm{j}3\varphi} \end{bmatrix} \\ &=\sqrt{2}\begin{bmatrix} 1 + e^{-\mathrm{j}3\varphi}\\ e^{-\mathrm{j}\varphi} + e^{-\mathrm{j}2\varphi}\\ \mathrm{j}(e^{-\mathrm{j}\varphi} - e^{-\mathrm{j}2\varphi})\\ \mathrm{j}(1 - e^{-\mathrm{j}3\varphi}) \end{bmatrix}\\ &=\sqrt{2}e^{-\mathrm{j}\frac{3}{2}\varphi}\begin{bmatrix} e^{\mathrm{j}\frac{3}{2}\varphi} + e^{-\mathrm{j}\frac{3}{2}\varphi}\\ e^{\mathrm{j}\frac{1}{2}\varphi} + e^{-\mathrm{j}\frac{1}{2}\varphi}\\ \mathrm{j}(e^{\mathrm{j}\frac{1}{2}\varphi} - e^{-\mathrm{j}\frac{1}{2}\varphi})\\ \mathrm{j}(e^{\mathrm{j}\frac{3}{2}\varphi} - e^{-\mathrm{j}\frac{3}{2}\varphi}) \end{bmatrix}\\ &=\sqrt{2}e^{-\mathrm{j}\frac{3}{2}\varphi}\begin{bmatrix} \cos(\frac{3}{2}\varphi) \\ \cos(\frac{1}{2}\varphi) \\ \sin(-\frac{1}{2}\varphi) \\ \sin(-\frac{3}{2}\varphi) \end{bmatrix}\\ &=e^{-\mathrm{j}\frac{3}{2}\varphi}\overline{\mathbf{a}}(\theta) \end{aligned} Qa(θ)?=2 ? ?100j?01j0?01?j0?100?j? ? ?1e?jφe?j2φe?j3φ? ?=2 ? ?1+e?j3φe?jφ+e?j2φj(e?jφ?e?j2φ)j(1?e?j3φ)? ?=2 ?e?j23?φ ?ej23?φ+e?j23?φej21?φ+e?j21?φj(ej21?φ?e?j21?φ)j(ej23?φ?e?j23?φ)? ?=2 ?e?j23?φ ?cos(23?φ)cos(21?φ)sin(?21?φ)sin(?23?φ)? ?=e?j23?φa(θ)?
??至此可以得到 $\mathbf{Qa}(\theta )=e^{?\mathrm{j}\frac{M?1}{2}\phi }\overline{\mathbf{a}}(\theta )$，進(jìn)一步我們可以得到：
$$\mathbf{QA}=e^{?\mathrm{j}\frac{M?1}{2}\phi }\overline{\mathbf{A}}$$
其中 $\overline{\mathbf{A}}\in {\mathbb{R}}^{M\times K}$ 是由 $K$ 個(gè)形如 $\overline{\mathbf{a}}(\theta )$ 的實(shí)向量組成的矩陣。最后一步，我們可以得到：
$$\begin{array}{rl}\mathbf{QR}{\mathbf{Q}}^H& =\left(e^{?\mathrm{j}\frac{M?1}{2}\phi }\overline{\mathbf{A}}\right){\mathbf{R}}_s\left(e^{\mathrm{j}\frac{M?1}{2}\phi }{\overline{\mathbf{A}}}^T\right)\\ & =\overline{\mathbf{A}}{\mathbf{R}}_s{\overline{\mathbf{A}}}^T\end{array}$$
因此有 $\mathrm{Re}(\mathbf{QR}{\mathbf{Q}}^H)=\overline{\mathbf{A}}\mathrm{Re}({\mathbf{R}}_s){\overline{\mathbf{A}}}^T$，不難看出即使 $\mathbf{R}$ 不滿足 Hermitian 及 Persymmetric 性質(zhì)，仍然不會(huì)破壞正交性。
??總的來(lái)說(shuō)，$\mathbf{Q}$ 的作用就是使得方向矢量轉(zhuǎn)為實(shí)向量，如此便可以利用協(xié)方差矩陣的實(shí)部進(jìn)行后續(xù)的計(jì)算。

算法步驟

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

1. 計(jì)算協(xié)方差矩陣 $\mathbf{R}=\frac{1}{T}\mathbf{X}{\mathbf{X}}^H$ 和酉矩陣 $\mathbf{Q}$，接著以 $\mathbf{R}?\mathbf{QR}{\mathbf{Q}}^H$ 更新協(xié)方差矩陣。
2. 對(duì) $\mathbf{R}$ 進(jìn)行特征值分解，并對(duì)特征值進(jìn)行排序，然后取得 $M?K$ 個(gè)較小特征值對(duì)應(yīng)的特征向量來(lái)組成噪聲子空間 ${\mathbf{U}}_n$。
3. 以下式遍歷 $\theta \in [?90^{°},90^{°}]$：
   $$\begin{array}{c}P(\theta )=\frac{1}{{\overline{\mathbf{a}}}^H(\theta ){\mathbf{U}}_n{\mathbf{U}}_n^T\overline{\mathbf{a}}(\theta )}\end{array}$$
   此時(shí)得到一組 $P(\theta )$，$K$ 個(gè)最大值對(duì)應(yīng)的 $\theta$ 就是需要返回的結(jié)果。

代碼實(shí)現(xiàn)

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
grid       = 0.1;             % search grid

%% Signals
R = generateSignal(M,K,T,theta,lambda,d,noise_flag,SNR);

%% DoA 
% unitary-MUSIC
[theta_unitary_music,P_unitary_music] = unitaryMUSIC(R,M,K,lambda,d,grid);

%% plot
figure;
hold on;
ang_list = -90:grid:90;
plot(ang_list, P_unitary_music);
hold off;

function [R,X,A,S] = generateSignal(M,K,T,theta,lambda,d,noise_flag,SNR)
    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
end

function [theta,P] = unitaryMUSIC(R,M,K,lambda,d,grid)
    M_half = floor(M/2);
    O = zeros(M_half,1);
    I = eye(M_half);
    J = fliplr(I);
    if mod(M,2) == 0
        Q = [I,J;1j*J,-1j*I]./sqrt(2);
    else
        Q = [I,O,J;O',sqrt(2),O';1j*J,O,-1j*I]./sqrt(2);
    end
    R = real(Q*R*Q');
    [U,~] = svd(R);
    Un = U(:, K+1:M);

    a_list = exp(-1j*(0:M-1).'*2*pi*d/lambda*sind(-90:grid:90));
    a_list = a_list.*exp(1j*(M-1)*ones(M,1)*pi*d/lambda*sind(-90:grid:90));
    a_list = real(Q*a_list);
    P = arrayfun(@(i) 1/norm(Un'*a_list(:,i),'fro'),1:size(a_list,2)); % spectral spectrum grid search
    P = 10*log10(P./max(P));

    [~, idx] = findpeaks(P,'NPeaks',K,'SortStr','descend'); % find K peaks
    theta = sort((idx-1)*grid-90);
end

參考內(nèi)容

1. Huarng K C, Yeh C C. A unitary transformation method for angle-of-arrival estimation[J]. IEEE Transactions on Signal Processing, 1991, 39(4): 975-977.
2. 【wikipedia】Persymmetric matrix
3. 【wikipedia】Hermitian matrix

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