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 的作用就是使得方向矢量转为实向量,如此便可以利用协方差矩阵的实部进行后续的计算。文章来源:https://www.toymoban.com/news/detail-828436.html
算法步骤
unitary MUSIC 算法步骤如下(输入为阵列接收矩阵 X \mathbf{X} X):文章来源地址https://www.toymoban.com/news/detail-828436.html
- 计算协方差矩阵 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 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
参考内容
- 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
到了这里,关于unitary MUSIC 算法的文章就介绍完了。如果您还想了解更多内容,请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章,希望大家以后多多支持TOY模板网!