MUSIC算法介绍
MUSIC算法介绍
MUSIC算法，全称为Multiple Signal Classification（多重信号分类），是用于波达方向估计的一种高分辨率谱估计算法。它最初由Schmidt在1979年提出，广泛应用于雷达、声纳、无线通信等领域中，以确定多个信号源的方向。该算法利用了 信号子空间和噪声子空间正交性 的特点，构造噪声空间然后通过谱峰搜索来检测信号的波达方向。
1. MUSIC 算法步骤
输入：多快拍接收数据 Y
∈
C
M
×
L
Y\in\mathbb{C}^{M\times L}
1）计算采样协方差矩阵 R
R
；
2）对 R
R
进行特征分解，求出对应特征值及特征矢量；
3）按照特征值大小排序，根据特征值判断信号数 K
K
并确定信号子空间 U
S
U_S
与噪声子空间 U
N
U_N
；
4）通过 P
(
θ
)
=
1
a
(
θ
)
H
U
N
U
N
H
a
(
θ
)
P(\theta)=\frac{1}{a(\theta)^HU_NU_N^Ha(\theta)}
得到空间谱估计；
5）对空间谱进行谱峰搜索，前 K
K
个极大值点对应的角度就是目标入射方向。
输出：目标入射角度
2. MUSIC 算法原理证明
首先，将观测数据向量的协方差矩阵特征分解
R
x
x
=
A
R
s
s
A
H
+
σ
2
I
=
U
Σ
U
H
=
[
U
S
U
N
]
[
Σ
S
O
O
Σ
N
]
[
U
S
H
U
N
H
]
\begin{aligned} R_{xx}=AR_{ss}A^H+\sigma^2I =U\Sigma U^H =\begin{bmatrix}U_S&U_N\end{bmatrix}\begin{bmatrix}\Sigma_S & O \cr O & \Sigma_N \end{bmatrix}\begin{bmatrix}U_S^H\cr U_N^H\end{bmatrix} \end{aligned}
右乘 U
N
U_N
有
R
x
x
U
N
=
[
U
S
,
U
N
]
[
Σ
S
O
O
Σ
N
]
[
U
S
H
U
N
H
]
U
N
=
[
U
S
,
U
N
]
[
Σ
S
O
O
Σ
N
]
[
O
I
]
=
σ
2
U
N
\begin{aligned} R_{xx}U_N =\begin{bmatrix}U_S,U_N\end{bmatrix}\begin{bmatrix}\Sigma_S & O \cr O & \Sigma_N \end{bmatrix}\begin{bmatrix}U_S^H\cr U_N^H\end{bmatrix}U_N =\begin{bmatrix}U_S,U_N\end{bmatrix}\begin{bmatrix}\Sigma_S & O \cr O & \Sigma_N \end{bmatrix}\begin{bmatrix}O\cr I\end{bmatrix} =\sigma^2U_N \end{aligned}
以及
R
x
x
U
N
=
A
R
s
s
A
H
U
N
+
σ
2
U
N
\begin{aligned} R_{xx}U_N=AR_{ss}A^HU_N+\sigma^2U_N \end{aligned}
联立两式得到
A
R
s
s
A
H
U
N
=
O
AR_{ss}A^HU_N=O
左乘 U
N
H
U_N^H
有
U
N
H
A
R
s
s
A
H
U
N
=
O
U_N^HAR_{ss}A^HU_N=O
若 Q 非奇异时 t
H
Q
t
=
0
t^HQt=0
，则 t
=
0
t=0
。故当
R
S
S
R_{SS}
非奇异，也就是入射信号不相干时，有
A
H
U
N
=
O
A^HU_N=O
也就是
[
a
(
ω
1
)
,
a
(
ω
2
)
,
…
,
a
(
ω
N
)
]
H
U
N
=
O
\begin{bmatrix}a(\omega_1),a(\omega_2),\dots,a(\omega_N)\end{bmatrix}^HU_N=O
即
a
(
ω
)
H
U
N
=
0
T
,
ω
=
ω
1
,
ω
2
,
…
,
ω
N
a
(
ω
)
H
U
N
≠
0
T
,
ω
≠
ω
1
,
ω
2
,
…
,
ω
N
\begin{aligned} a(\omega)^HU_N=0^T,\quad\omega=\omega_1,\omega_2,\dots,\omega_N\cr a(\omega)^HU_N\ne0^T,\quad\omega\ne\omega_1,\omega_2,\dots,\omega_N \end{aligned}
写成标量形式，定义一种类似功率谱的函数
P
(
ω
)
=
1
a
(
ω
)
H
U
N
U
N
H
a
(
ω
)
P(\omega)=\frac{1}{a(\omega)^HU_NU_N^Ha(\omega)}
上式的N个峰值对应N个信号的来波方向
3. MATLAB仿真
clear all
close all
theta = [10 30 60]; % 信源入射角度
source_number = length(theta); % 信元数
sensor_number = 8; % 阵元数
w = pi/4; % 信号频率
lambda = (2*pi*3e8)/w; % 信号波长
d = 0.5*lambda; % 阵元间距
snr = 10; % 信噪比为10
n = 500; % 快拍数500
A = exp(-1j*2*pi*(0:d:(sensor_number-1)*d).'*sind(theta)/lambda); % 阵列流形
S = randn(source_number,n); % 随机数矩阵
X = A*S; % 信号矩阵
X1 = awgn(X,snr,'measured'); % 添加噪声
Rxx = X1*X1'/n; % 协方差矩阵
InvS = inv(Rxx);
[EV,D] = eig(Rxx); % 协方差矩阵对角化 EV特征向量矩阵 D对角阵
EVA = diag(D)'; % 提取对角元素
[EVA,I] = sort(EVA); % EVA对角元素从小到大排序 I位置索引
EVA = fliplr(EVA); % 左右翻转
EV = fliplr(EV(:,I)); % 使用索引I对EV列重新排列 左右翻转
for k = 1:361
angle(k)=(k-181)/2;
a = exp(-1j*2*pi*(0:d:(sensor_number-1)*d)*sind((k-181)/2)/lambda).'; % 方向向量
En=EV(:,source_number+1:sensor_number); % 提取source_number+1列到sensor_number列
SP(k)=(a'*a)/(a'*En*En'*a); % 存储计算值
end
SP=abs(SP); % SP的绝对值
SPmax=max(SP);
SP=10*log10(SP/SPmax); % 以dB为单位
h=plot(angle,SP); % 绘制angle与SP的关系图
set(h,'Linewidth',2) % 绘制线条宽度为 2
xlabel('angle(degree)') % 为 x 轴和 y 轴添加标签
ylabel('magnituded(dB)')
title('Music Algorithm For Doa', 'fontsize', 16);
set(gca,'XTick',[-90:30:90]) % x 轴刻度
grid on
运行结果：
这里可以看到，music算法能够识别多个方向，而且具有很高的分辨率。需要注意的是，当入射信号相干时，music算法会失效，在推导过程中也是以入射信号不相干为前提的。对于相干信号，可采用平滑music算法进行处理，这种方法在之后的文章中进行介绍。
4. 不同信噪比下的RMSE
在信源数为1，阵元数为8，间距 λ
2
\frac{\lambda}{2}
，蒙托卡罗次数为1000 时的 RMSE
clear all
close all
source_number = 1; % 信元数
sensor_number = 8; % 阵元数
w = pi/4; % 信号频率
lambda = (2*pi*3e8)/w; % 信号波长
d = 0.5*lambda; % 阵元间距
kps = 100; % 快拍数
theta = 30; % 信源入射角度
mengtimes = 1000; % 蒙特卡洛次数
A = exp(-1j*2*pi*(0:d:(sensor_number-1)*d).'*sind(theta)/lambda); % 阵列流形
S = randn(source_number,kps); % 随机数矩阵
X = A*S; % 信号矩阵
snr = 0:1:25; % 信噪比区间
eror = zeros(length(snr),1); % 估计误差
for n = 1:size(snr,2)
for m = 1:mengtimes
X1 = awgn(X,snr(n),'measured'); % 添加噪声
Rxx = X1*X1'/kps; % 协方差矩阵
InvS = inv(Rxx);
[EV,D] = eig(Rxx); % 协方差矩阵对角化 EV特征向量矩阵 D对角阵
EVA = diag(D)'; % 提取对角元素
[EVA,I] = sort(EVA); % EVA对角元素从小到大排序 I位置索引
EVA = fliplr(EVA); % 左右翻转
EV = fliplr(EV(:,I)); % 使用索引I对EV列重新排列 左右翻转
for k = 1:3601
a = exp(-1j*2*pi*(0:d:(sensor_number-1)*d)*sind((k-1801)/20)/lambda).'; % 方向向量
En = EV(:,source_number+1:sensor_number); % 提取source_number+1列到sensor_number列
SP(k) = (a'*a)/(a'*En*En'*a); % 存储计算值
end
SP = abs(SP); % SP的绝对值
[maxsp,index] = max(SP); % 搜索谱峰最大值坐标
estmt(m) = (index-1801)/20;
eror(n) = (estmt(m)-theta)^2 + eror(n); % 误差平方和
end
rmse(n) = sqrt(eror(n)/mengtimes);
end
figure
plot(snr,rmse,'Linewidth',2)
grid on
xlabel('SNR(dB)','Fontname','Times New Roman','FontSize',17)
ylabel('\fontname{Times New Roman}\fontsize{17}RMSE°')
legend('\fontname{Times New Roman}\fontsize{17}MUSIC','Location','Best')
运行结果：
参考：
Music算法原理部分参考于：《矩阵分析与应用》张贤达