unitary MUSIC 算法
论文 A Unitary Transformation Method for Angle-of-Arrival Estimation 中提出了 unitary MUSIC 的算法，直译就是酉 MUSIC 算法，即酉变换 MUSIC 算法。该算法的目的是简化计算复杂度，将传统 MUSIC 算法中的复数 SVD 和复数网格搜索计算转化为实数计算。在学习 unitary MUSIC 之前需要理解 Hermitian 矩阵及 Persymmetric 矩阵的概念及性质：
- Hermitian 矩阵指的是满足 A H = A \mathbf{A}^H = \mathbf{A} AH=A 的矩阵 A \mathbf{A} A；
- Persymmetric 矩阵指的是满足 A J = J A T \mathbf{A}\mathbf{J} = \mathbf{J}\mathbf{A}^T AJ=JAT 的矩阵 A \mathbf{A} A，其中 J \mathbf{J} J 其对角线从左下至右上，很多地方又称 J \mathbf{J} J 为选择矩阵。
- 假若矩阵 A \mathbf{A} A 既为 Hermitian 矩阵又为 Persymmetric 矩阵，则满足：
J A ∗ J = A \mathbf{J}\mathbf{A}^*\mathbf{J}=\mathbf{A} JA∗J=A
其中 A ∗ \mathbf{A}^* A∗ 为 A \mathbf{A} A 的共轭。
在接下来的讨论中， I \mathbf{I} I 和 J \mathbf{J} J 分别用作表示单位矩阵和选择矩阵，下文中将会出现这两种矩阵的运算，例如 A I \mathbf{A}\mathbf{I} AI 或 J B \mathbf{J}\mathbf{B} JB，设 A \mathbf{A} A 和 B \mathbf{B} B 均为方阵，如果没有特别强调，则说明 I \mathbf{I} I 和 J \mathbf{J} J 分别和 A \mathbf{A} A 和 B \mathbf{B} B 同维度。
算法原理
前面讨论的子空间算法中，复协方差矩阵的特征值分解是至关重要的一步，然而该步的计算量很高。为了降低计算量，unitary MUSIC 算法考虑利用一个酉矩阵将原先的复协方差矩阵 R \mathbf{R} R 转换成实协方差矩阵，同时传统算法中的复空间搜索向量 a ( θ ) \mathbf{a}(\theta) a(θ) 也用实向量来代替。
unitary MUSIC 算法的提出基于一个性质，即若不相关的窄带远场信号源射入均匀线阵中，其协方差矩阵不仅是 Hermitian，且 Persymmetric。通常估计的协方差矩阵 R ≜ R ^ \mathbf{R}\triangleq \hat{\mathbf{R}} R≜R^ 仅仅只是 Hermitian 矩阵但不满足 Persymmetric 性质，需要先获得一个满足 Persymmetric 性质的估计协方差矩阵：
R ≜ 1 2 ( R ^ + J R ^ ∗ J ) \begin{equation*} \mathbf{R}\triangleq \frac{1}{2}(\hat{\mathbf{R}}+ \mathbf{J}\hat{\mathbf{R}}^*\mathbf{J}) \end{equation*} R≜21(R^+JR^∗J)
假设阵元数 M M M 为偶数，unitary MUSIC 算法引入了一个酉矩阵 Q ∈ C M × M \mathbf{Q}\in\mathbb{C}^{M\times M} Q∈CM×M：
Q = 1 2 [ I J j J − j I ] \mathbf{Q} = \frac{1}{\sqrt{2}} \begin{bmatrix} \mathbf{I} & \mathbf{J} \\ \mathrm{j}\mathbf{J} & -\mathrm{j}\mathbf{I} \end{bmatrix} Q=21[IjJJ−jI]
其中 I \mathbf{I} I 和 J \mathbf{J} J 分别为单位矩阵和选择矩阵，且该两个矩阵维度均为 M 2 × M 2 \frac{M}{2}\times \frac{M}{2} 2M×2M。易得 Q \mathbf{Q} Q 为酉矩阵，即 Q − 1 = Q H \mathbf{Q}^{-1} = \mathbf{Q}^H Q−1=QH，同时满足：
Q ∗ J = Q \mathbf{Q}^*\mathbf{J} = \mathbf{Q} Q∗J=Q
至此到了本算法的关键，它在于证明由于 R \mathbf{R} R 是 Hermitian 且 Persymmetric 矩阵， Q R Q H \mathbf{Q}\mathbf{R}\mathbf{Q}^H QRQH 是实对称矩阵：
因为 R \mathbf{R} R 为 Hermitian，易得 Q R Q H \mathbf{Q}\mathbf{R}\mathbf{Q}^H QRQH 为 Hermitian；因此只需要证明 Q R Q H \mathbf{Q}\mathbf{R}\mathbf{Q}^H QRQH 是实矩阵，即证明 ( Q R Q H ) ∗ = Q R Q H (\mathbf{Q}\mathbf{R}\mathbf{Q}^H)^* = \mathbf{Q}\mathbf{R}\mathbf{Q}^H (QRQH)∗=QRQH：
( Q R Q H ) ∗ = Q ∗ R ∗ Q T = ( Q ∗ J ) ( J R ∗ J ) ( J Q T ) = Q R Q H \begin{equation*} \begin{aligned} &(\mathbf{Q}\mathbf{R}\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{Q}\mathbf{R}\mathbf{Q}^H \end{aligned} \end{equation*} ===(QRQH)∗Q∗R∗QT(Q∗J)(JR∗J)(JQT)QRQH
由此得证。
综上所述，unitary MUSIC 算法引入酉矩阵 Q \mathbf{Q} Q 并令 R ≜ Q R Q H \mathbf{R}\triangleq\mathbf{Q}\mathbf{R}\mathbf{Q}^H R≜QRQH，使得酉变换后的协方差矩阵变为实对称矩阵，接着对其特征值分解即可进行后续的搜索步骤。而 ULA 的搜索方向矢量为：
a ( θ ) = [ 1 , e − j 2 π d sin θ / λ , ⋯ , e − j ( M − 1 ) 2 π d sin θ / λ ] T \mathbf{a}(\theta) = \left[1, e^{-\mathrm{j}2\pi d\sin\theta/\lambda},\cdots,e^{-\mathrm{j}(M-1)2\pi d\sin\theta/\lambda}\right]^T a(θ)=[1,e−j2πdsinθ/λ,⋯,e−j(M−1)2πdsinθ/λ]T
则 unitary MUSIC 算法的搜索方向矢量为 a ( θ ) ≜ Q a ( θ ) \mathbf{a}(\theta)\triangleq\mathbf{Q}\mathbf{a}(\theta) a(θ)≜Qa(θ)。
为了进一步降低算法计算复杂度，Unitary MUSIC 算法考虑将搜索方向矢量也用实变量代替，做法如下：
a ( θ ) ≜ e j M − 1 2 2 π d sin θ / λ Q a ( θ ) = Q [ e j M − 1 2 2 π d sin θ / λ , ⋯ , e j 1 2 2 π d sin θ / λ , e − j 1 2 2 π d sin θ / λ , ⋯ , e − j M − 1 2 2 π d sin θ / λ ] T \begin{aligned} \mathbf{a}(\theta) &\triangleq e^{j\frac{M-1}{2}2\pi d\sin\theta/\lambda}\mathbf{Q}\mathbf{a}(\theta) \\ &= \mathbf{Q}\left[e^{\mathrm{j}\frac{M-1}{2}2\pi d\sin\theta/\lambda},\cdots, e^{\mathrm{j}\frac{1}{2}2\pi d\sin\theta/\lambda},e^{-\mathrm{j}\frac{1}{2}2\pi d\sin\theta/\lambda},\cdots,e^{-\mathrm{j}\frac{M-1}{2}2\pi d\sin\theta/\lambda}\right]^T \end{aligned} a(θ)≜ej2M−12πdsinθ/λQa(θ)=Q[ej2M−12πdsinθ/λ,⋯,ej212πdsinθ/λ,e−j212πdsinθ/λ,⋯,e−j2M−12πdsinθ/λ]T
不难看出原方向矢量对应的阵列索引位置为 { 0 , 1 , ⋯ , M − 1 } \{0,1,\cdots,M-1\} {0,1,⋯,M−1}，而更新后方向矢量对应的阵列索引位置为 { − 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}。此时 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−12πdsinθ/λQa(θ)=2 cos(2M−12πdsinθ/λ)⋮cos(212πdsinθ/λ)sin(−212πdsinθ/λ)⋮sin(−2M−12πdsinθ/λ)
当 M M M 为奇数时，酉矩阵 Q \mathbf{Q} 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=21 IOTjJO2OJOT−jI
其中 O \mathbf{O} O 为零矩阵。
进一步理解
在上一小节中，我整理了论文对于 unitary MUSIC 算法的解释及证明，其核心思想在于酉矩阵 Q \mathbf{Q} Q 的提出。在本小节中，我将进一步对 Q \mathbf{Q} Q 的作用谈谈自己的理解。
在论文中，unitary 算法一直强调 R \mathbf{R} R 的 Hermitian 及 Persymmetric 性质，因为假若 R \mathbf{R} R 不满足这两个性质， Q R Q H \mathbf{Q}\mathbf{R}\mathbf{Q}^H QRQH 将不是实矩阵，但是利用 R e ( Q R Q H ) \mathrm{Re}(\mathbf{Q}\mathbf{R}\mathbf{Q}^H) Re(QRQH) 来进行后续的估计其实是可以估计角度的。假设 M = 4 M=4 M=4 和 φ = 2 π d sin θ / λ \varphi = 2\pi d \sin\theta/\lambda φ=2πdsinθ/λ，展开 Q a ( θ ) \mathbf{Q}\mathbf{a}(\theta) Qa(θ)：
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 100j01j001−j0100−j 1e−jφe−j2φe−j3φ =2 1+e−j3φe−jφ+e−j2φj(e−jφ−e−j2φ)j(1−e−j3φ) =2e−j23φ ej23φ+e−j23φej21φ+e−j21φj(ej21φ−e−j21φ)j(ej23φ−e−j23φ) =2e−j23φ cos(23φ)cos(21φ)sin(−21φ)sin(−23φ) =e−j23φa(θ)
至此可以得到 Q a ( θ ) = e − j M − 1 2 φ a ‾ ( θ ) \mathbf{Q}\mathbf{a}(\theta)=e^{-\mathrm{j}\frac{M-1}{2}\varphi}\overline{\mathbf{a}}(\theta) Qa(θ)=e−j2M−1φa(θ)，进一步我们可以得到：
Q A = e − j M − 1 2 φ A ‾ \mathbf{Q}\mathbf{A}=e^{-\mathrm{j}\frac{M-1}{2}\varphi}\overline{\mathbf{A}} QA=e−j2M−1φA
其中 A ‾ ∈ R M × K \overline{\mathbf{A}}\in\mathbb{R}^{M\times K} A∈RM×K 是由 K K K 个形如 a ‾ ( θ ) \overline{\mathbf{a}}(\theta) a(θ) 的实向量组成的矩阵。最后一步，我们可以得到：
Q R Q H = ( e − j M − 1 2 φ A ‾ ) R s ( e j M − 1 2 φ A ‾ T ) = A ‾ R s A ‾ T \begin{aligned} \mathbf{Q}\mathbf{R}\mathbf{Q}^H&=\left(e^{-\mathrm{j}\frac{M-1}{2}\varphi}\overline{\mathbf{A}}\right)\mathbf{R}_s\left(e^{\mathrm{j}\frac{M-1}{2}\varphi}\overline{\mathbf{A}}^T\right)\\ &=\overline{\mathbf{A}}\mathbf{R}_s\overline{\mathbf{A}}^T \end{aligned} QRQH=(e−j2M−1φA)Rs(ej2M−1φAT)=ARsAT
因此有 R e ( Q R Q H ) = A ‾ R e ( R s ) A ‾ T \mathrm{Re}(\mathbf{Q}\mathbf{R}\mathbf{Q}^H) = \overline{\mathbf{A}}\mathrm{Re}(\mathbf{R}_s)\overline{\mathbf{A}}^T Re(QRQH)=ARe(Rs)AT，不难看出即使 R \mathbf{R} R 不满足 Hermitian 及 Persymmetric 性质，仍然不会破坏正交性。
总的来说， Q \mathbf{Q} Q 的作用就是使得方向矢量转为实向量，如此便可以利用协方差矩阵的实部进行后续的计算。
算法步骤
unitary MUSIC 算法步骤如下（输入为阵列接收矩阵 X \mathbf{X} X）：
- 计算协方差矩阵 R = 1 T X X H \mathbf{R} = \frac{1}{T} \mathbf{X}\mathbf{X}^H R=T1XXH 和酉矩阵 Q \mathbf{Q} Q，接着以 R ≜ Q R Q H \mathbf{R}\triangleq\mathbf{Q}\mathbf{R}\mathbf{Q}^H R≜QRQH 更新协方差矩阵。
- 对 R \mathbf{R} R 进行特征值分解，并对特征值进行排序，然后取得 M − K M-K M−K 个较小特征值对应的特征向量来组成噪声子空间 U n \mathbf{U}_n Un。
- 以下式遍历 θ ∈ [ − 9 0 ∘ , 9 0 ∘ ] \theta \in [-90^{\circ}, 90^{\circ}] θ∈[−90∘,90∘]：
P ( θ ) = 1 a ‾ H ( θ ) U n U n T a ‾ ( θ ) \begin{equation*} P(\theta) = \frac{1}{\overline{\mathbf{a}}^H(\theta)\mathbf{U}_n\mathbf{U}_n^T\overline{\mathbf{a}}(\theta)} \end{equation*} P(θ)=aH(θ)UnUnTa(θ)1
此时得到一组 P ( θ ) P(\theta) P(θ)， K K K 个最大值对应的 θ \theta θ 就是需要返回的结果。
代码实现
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 matrixA = exp(-1j*(0:M-1)'*2*pi*d/lambda*sind(theta)); % steering vector matrixN = noise_flag.*sqrt(10.^(-SNR/10))*(1/sqrt(2))*(randn(M,T)+1j*randn(M,T)); % noise matrixX = A*S+N; % received matrixR = X*X'/T; % covariance matrix
endfunction [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) == 0Q = [I,J;1j*J,-1j*I]./sqrt(2);elseQ = [I,O,J;O',sqrt(2),O';1j*J,O,-1j*I]./sqrt(2);endR = 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 searchP = 10*log10(P./max(P));[~, idx] = findpeaks(P,'NPeaks',K,'SortStr','descend'); % find K peakstheta = sort((idx-1)*grid-90);
end
参考内容
- 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.
- 【wikipedia】Persymmetric matrix
- 【wikipedia】Hermitian matrix
-----CityscapesScripts生成自己的标签)
)
)
