DoA 估计：多重信号分类 MUSIC 算法（附 MATLAB 代码）

2023-05-06

amp 信号 dots

本文首次在公众号【零妖阁】上发表，为了方便阅读和分享，我们将在其他平台进行自动同步。由于不同平台的排版格式可能存在差异，为了避免影响阅读体验，建议如有排版问题，可前往公众号查看原文。感谢您的阅读和支持！DoA估计是指根据天线阵列的接收信号估计出单个或多个信号源的方位信息。由于激励信号和方向图之间存在

本文首次在公众号【零妖阁】上发表，为了方便阅读和分享，我们将在其他平台进行自动同步。由于不同平台的排版格式可能存在差异，为了避免影响阅读体验，建议如有排版问题，可前往公众号查看原文。感谢您的阅读和支持！

DoA 估计是指根据天线阵列的接收信号估计出单个或多个信号源的方位信息。由于激励信号和方向图之间存在傅里叶关系，DoA 估计也可以等效为谱估计问题。

多重信号分类（Mutiple Signal Classification）算法，简称 MUSIC 算法，是一种常用的 DoA 估计方法。它的基本思想是将任意阵列输出数据的协方差矩阵进行特征分解，从而得到与信号分量相对应的信号子空间和与信号分量相正交的噪声子空间。

信号模型

假设空间中存在 $M$ 个不同方向的信号，入射到由 $N$ 个天线单元构成的均匀直线阵上。第 $i$ 个信号源的方向为 $\phi_i$ （ $1,\dots,M$ ），第 $i$ 个信号源的信号为 $\alpha_i(t)$ 。假设 $M < N$ 。

令第 $n$ 个天线单元的噪声为 $n_n(t)$ 。在窄带远场条件下，第 $n$ 个天线单元的输出信号 $x_n(t)$ 可表示为

\begin{aligned} x_{n} (t) & = \sum_{i = 1}^{M} α_{i} (t) e^{j k z_{n} \sin ϕ_{i}} + n_{n} (t) \\ = \sum_{i = 1}^{M} α_{i} (t) s_{n} (ϕ_{i}) + n_{n} (t) \end{aligned}

$\begin{aligned} x_n(t) &= \sum\limits_{i=1}^{M} \alpha_i(t) e^{jkz_n\sin\phi_i} + n_n(t) \\ &= \sum\limits_{i=1}^{M} \alpha_i(t) s_n(\phi_i) + n_n(t) \end{aligned}$

x_{n} (t) = i = 1 \sum M α_{i} (t) e^{jk z_{n} s i n ϕ_{i}} + n_{n} (t) = i = 1 \sum M α_{i} (t) s_{n} (ϕ_{i}) + n_{n} (t)

将 $N$ 个天线的输出信号表示成向量形式 $\bf x (t)$ ，则上式可归纳为
$\mathbf{x}(t) = \mathbf{S} \mathbf{\alpha}(t) + \mathbf{n}(t)$

其中， $\mathbf{S}$ 为阵列的流型矩阵，矩阵规模为 $N\times M$ ，具体可表示为 $M$
个不同方向对应的阵列导向矢量：
$\mathbf{S} = [\mathbf{s}(\phi_1), \mathbf{s}(\phi_2), \dots, \mathbf{s}(\phi_M) ]$

由于 $M < N$ ，流型矩阵 $\mathbf{S}$ 为列满秩矩阵， $\mathrm{Rank}(\mathbf{S}) = M$ 。

MUSIC 算法思想

假设不同信号源的信号之间是相互正交的，噪声与信号之间是正交的，则阵列输出信号 $\mathbf x(t)$ 的协方差矩阵为

\begin{aligned} R & = E [x (t) x (t)^{H}] \\ = E [S α (t) α (t)^{H} S^{H} + n (t) n (t)^{H}] \\ = S A S^{H} + σ^{2} I \\ = R_{S} + σ^{2} I \end{aligned}

$\begin{aligned} \mathbf R &= \mathrm E [\mathbf x(t) \mathbf x(t)^H] \\ &= \mathrm E [\mathbf{S} \mathbf{\alpha}(t) \mathbf{\alpha}(t)^H \mathbf{S}^H + \mathbf{n}(t)\mathbf{n}(t)^H ] \\ &= \mathbf{S} \mathbf A \mathbf{S}^H + \sigma^2 \mathbf I \\ &= \mathbf{R}_S + \sigma^2 \mathbf I \\ \end{aligned}$

R = E [x (t) x (t)^{H}] = E [S α (t) α (t)^{H} S^{H} + n (t) n (t)^{H}] = SA S^{H} + σ^{2} I = R_{S} + σ^{2} I

其中， $\mathbf A$ 为不同信号源之间的协方差矩阵，由于不同信号源之间是相互正交的， $\mathbf A$ 为正定对角矩阵：
$\mathbf A = \left[$

\begin{matrix} E [| α_{1} (t) |^{2}] & \dots & \dots & \dots \\ \dots & E [| α_{2} (t) |^{2}] & \dots & \dots \\ \dots & \dots & \dots & \dots \\ \dots & \dots & \dots & E [| α_{M} (t) |^{2}] \end{matrix}

$\begin{matrix} &\mathrm E [\mathbf |\alpha_1(t)|^2] &\dots &\dots &\dots \\ &\dots &\mathrm E [\mathbf |\alpha_2(t)|^2] &\dots &\dots \\ &\dots &\dots &\dots &\dots \\ &\dots &\dots &\dots &\mathrm E [\mathbf |\alpha_M(t)|^2] \\ \end{matrix}$ \right]

A = E [∣ α_{1} (t) ∣^{2}] \dots \dots \dots \dots E [∣ α_{2} (t) ∣^{2}] \dots \dots \dots \dots \dots \dots \dots \dots \dots E [∣ α_{M} (t) ∣^{2}]

由于信号协方差矩阵 $\mathbf R_S$ 的规模为 $N\times N$ ，秩为 $M$ ， $\mathbf R_S$ 存在 $N - M$ 个特征值为 0 的特征向量，令这种特征向量为 $\mathbf q_m$ ，则
$\mathbf R_S \mathbf q_m = 0$
$\Rightarrow \mathbf{S} \mathbf A \mathbf{S}^H \mathbf q_m = 0$
$\Rightarrow \mathbf q_m^H \mathbf{S} \mathbf A \mathbf{S}^H \mathbf q_m = 0$
$\Rightarrow \mathbf{S}^H \mathbf q_m = 0$

上述推论说明， $\mathbf R_S$ 的特征值为 0 时对应的特征向量 $\mathbf q_m$ 与信号源对应的 $M$ 个导向矢量均正交。

令 $\mathbf R_S$ 的 $N - M$ 个特征值为 0 时对应的特征向量构成矩阵 $\mathbf Q_n $，其规模为 $\times (N-M)$ ，则
$\mathbf{S}^H \mathbf Q_n = 0$

则 MUSIC 算法的谱估计公式为
$P_{\mathrm{MUSIC}}(\phi) = \frac{1}{|| \mathbf Q_n^H \mathbf s(\phi) ||^2}$

当上式中的 $\phi$ 与信号源方向相同时，分母为零，此时 MUSIC 谱估计为无穷大。因此，MUSIC 谱估计的尖峰数目与信源数目相同，尖峰对应的方向即为信号源的方向。

如何根据阵列输出信号 $\mathbf x$ 计算 $\mathbf Q_n$ ？

通过记录多组阵列输出信号快拍，可以计算出输出信号协方差矩阵的近似值
$\mathbf R = \frac{1}{K} \sum\limits_{k=1}^K \mathbf x_k \mathbf x_k^H$

那么，如何根据输出信号的协方差矩阵 $\mathbf R$ 估计出信号协方差矩阵 $\mathbf R_S$ 对应的特征值为 0 的特征向量矩阵 $\mathbf Q_n$ 呢？

对于 $\mathbf R_S$ 的任意特征向量 $\mathbf q_m \in \mathbf Q$ ，有
$\mathbf R_S \mathbf q_m = \lambda_m \mathbf q_m$

\begin{aligned} \Rightarrow R q_{m} & = R_{S} q_{m} + σ^{2} I q_{m} \\ = (λ_{m} + σ^{2}) q_{m} \end{aligned}

$\begin{aligned} \Rightarrow \mathbf R \mathbf q_m &= \mathbf R_S \mathbf q_m + \sigma^2 \mathbf I \mathbf q_m \\ &= (\lambda_m + \sigma^2)\mathbf q_m \end{aligned}$

\Rightarrow R q_{m} = R_{S} q_{m} + σ^{2} I q_{m} = (λ_{m} + σ^{2}) q_{m}

因此，信号协方差矩阵 $\mathbf R_S$ 的特征值 $\lambda_m$ 对应的特征向量与输出信号协方差矩阵 $\mathbf R$ 的特征值 $\lambda_m+\sigma^2$ 对应的特征向量相同。

因此， $\mathbf R$ 的特征分解可表示为

\begin{aligned} R & = Q (Λ + σ^{2} I) Q^{H} \\ = Q [\begin{matrix} λ_{1} + σ^{2} & 0 & \dots & 0 & 0 & \dots & 0 \\ 0 & λ_{2} + σ^{2} & \dots & 0 & 0 & \dots & 0 \\ \dots & \dots & \dots & \dots & \dots & \dots & \dots \\ 0 & 0 & \dots & λ_{M} + σ^{2} & 0 & \dots & 0 \\ 0 & 0 & \dots & 0 & σ^{2} & \dots & 0 \\ \dots & \dots & \dots & \dots & \dots & \dots & \dots \\ 0 & 0 & \dots & 0 & 0 & \dots & σ^{2} \end{matrix}] Q^{H} \end{aligned}

$\begin{aligned} \mathbf R &= \mathbf Q (\mathbf \Lambda + \sigma^2 \mathbf I) \mathbf Q^H \\ &= \mathbf Q \left[\begin{matrix} &\lambda_1 + \sigma^2 &0 &\dots &0 &0 &\dots &0 \\ &0 &\lambda_2 + \sigma^2 &\dots &0 &0 &\dots &0 \\ &\dots &\dots &\dots &\dots &\dots &\dots &\dots \\ &0 &0 &\dots &\lambda_M + \sigma^2 &0 &\dots &0 \\ &0 &0 &\dots &0 &\sigma^2 &\dots &0 \\ &\dots &\dots &\dots &\dots &\dots &\dots &\dots \\ &0 &0 &\dots &0 &0 &\dots &\sigma^2 \\ \end{matrix}\right] \mathbf Q^H \end{aligned}$

R = Q (Λ + σ^{2} I) Q^{H} = Q λ_{1} + σ^{2} 0 \dots 00 \dots 0 0 λ_{2} + σ^{2} \dots 00 \dots 0 \dots \dots \dots \dots \dots \dots \dots 00 \dots λ_{M} + σ^{2} 0 \dots 0 00 \dots 0 σ^{2} \dots 0 \dots \dots \dots \dots \dots \dots \dots 00 \dots 00 \dots σ^{2} Q^{H}

上式表明，将输出信号矩阵 $\mathbf R$ 进行特征分解，得到的 $N - M$ 个较小且相等的特征值对应的特征向量即可构成 $\mathbf Q_n$ 。

MATLAB 仿真

clc; clear; close all;

%% 参数设置
%%% 工作频率
c = 3e8;
freq = 10e9;
lambda = c/freq;    % 波长
k = 2*pi/lambda;    % 波数
%%% 阵列参数
N = 10;                 % 阵元数量
d = 0.5*lambda;         % 阵元间隔 
z = (0:d:(N-1)*d)';     % 阵元坐标分布
%%% 信号源参数
phi = [-10, -30, 60]'*pi/180;   % 来波方向
M = length(phi);                % 信号源数目
%%% 仿真参数
SNR = 10;             % 信噪比(dB)
K = 1000;     % 采样点数

%% 阵列接收信号仿真模拟
S = exp(1j*k*z*sin(phi'));          % 流型矩阵
Alpha = randn(M, K);         % 输入信号
X = S*Alpha;                        % 阵列接收信号
X1 = awgn(X, SNR, 'measured');      % 加载高斯白噪声

%% MUSIC 算法
%%% 阵列接收信号的协方差矩阵的特征分解
R = X1*X1'/K;    % 阵列接收信号的协方差矩阵
[EV, D] = eig(R);       % 特征值分解
EVA = diag(D);          % 提取特征值
[EVA, I] = sort(EVA, 'descend');   % 降序排序
Q = EV(:, I);           % 特征向量构成的矩阵
Q_n = Q(:, M+1:N);      % 噪声子空间
%%% 计算MUSIC谱估计函数
phi_list = linspace(-pi/2, pi/2, 200)';
S1 = exp(1j*k*z*sin(phi_list'));    % 不同方向对应的流型矢量构成矩阵
P_MUSIC = 1./sum(abs(Q_n'*S1).^2);     % MUSIC 谱估计公式
%%% 转换为dB
P_MUSIC = abs(P_MUSIC);
P_MUSIC_max = max(P_MUSIC);
P_MUSIC_dB = 10*log10(P_MUSIC/P_MUSIC_max);
%%% 提取信号源方向
[P_peaks, P_peaks_idx] = findpeaks(P_MUSIC_dB);     % 提取峰值
[P_peaks, I] = sort(P_peaks, 'descend');    % 峰值降序排序
P_peaks_idx = P_peaks_idx(I);
P_peaks = P_peaks(1:M);             % 提取前M个
P_peaks_idx = P_peaks_idx(1:M);
phi_e = phi_list(P_peaks_idx)*180/pi;   % 估计方向
disp('信号源估计方向为：');
disp(phi_e);
%%% 绘图
figure;
plot(phi_list*180/pi, P_MUSIC_dB, 'k', 'Linewidth', 2);
xlabel('\phi (deg)');
ylabel('Pseudo-spectrum (dB)');
axis([-90, 90, -40, 0]);
grid on;
hold on;
plot(phi_e, P_peaks, 'r.', 'MarkerSize', 25);
hold on;
for idx = 1:M
    text(phi_e(idx)+3, P_peaks(idx), sprintf('%0.1f°', phi_e(idx)));
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63