本文首次在公眾號【零妖閣】上發(fā)表，為了方便閱讀和分享，我們將在其他平臺進(jìn)行自動同步。由于不同平臺的排版格式可能存在差異，為了避免影響閱讀體驗(yàn)，建議如有排版問題，可前往公眾號查看原文。感謝您的閱讀和支持！

DoA 估計是指根據(jù)天線陣列的接收信號估計出單個或多個信號源的方位信息。由于激勵信號和方向圖之間存在傅里葉關(guān)系，DoA 估計也可以等效為譜估計問題。

多重信號分類（Mutiple Signal Classification）算法，簡稱 MUSIC 算法，是一種常用的 DoA 估計方法。它的基本思想是將任意陣列輸出數(shù)據(jù)的協(xié)方差矩陣進(jìn)行特征分解，從而得到與信號分量相對應(yīng)的信號子空間和與信號分量相正交的噪聲子空間。

信號模型

假設(shè)空間中存在 $M$ 個不同方向的信號，入射到由 $N$ 個天線單元構(gòu)成的均勻直線陣上。第 $i$ 個信號源的方向?yàn)?${?}_i$（$i=1,\dots ,M$），第 $i$ 個信號源的信號為 ${\alpha }_i(t)$。假設(shè) $M<N$ 。

令第 $n$ 個天線單元的噪聲為 $n_n(t)$。在窄帶遠(yuǎn)場條件下，第 $n$ 個天線單元的輸出信號 $x_n(t)$ 可表示為
$$\begin{array}{rl}{x}_n(t)& =\sum _{i=1}^M{\alpha }_i(t)e^{jk{z}_n\sin {?}_i}+n_n(t)\\ & =\sum _{i=1}^M{\alpha }_i(t)s_n({?}_i)+n_n(t)\end{array}$$

將 $N$ 個天線的輸出信號表示成向量形式 $\mathbf{x}(\mathbf{t})$，則上式可歸納為
$$\mathbf{x}(t)=\mathbf{S}\alpha (t)+\mathbf{n}(t)$$

其中，$\mathbf{S}$ 為陣列的流型矩陣，矩陣規(guī)模為 $N\times M$，具體可表示為 $M$
個不同方向?qū)?yīng)的陣列導(dǎo)向矢量：
$$\mathbf{S}=[\mathbf{s}({?}_1),\mathbf{s}({?}_2),\dots ,\mathbf{s}({?}_M)]$$

由于 $M<N$，流型矩陣 $\mathbf{S}$ 為列滿秩矩陣，$\mathrm{Rank}(\mathbf{S})=M$。

MUSIC 算法思想

假設(shè)不同信號源的信號之間是相互正交的，噪聲與信號之間是正交的，則陣列輸出信號 $\mathbf{x}(t)$ 的協(xié)方差矩陣為
$$\begin{array}{rl}\mathbf{R}& =\mathrm{E}[\mathbf{x}(t)\mathbf{x}(t{)}^H]\\ & =\mathrm{E}[\mathbf{S}\alpha (t)\alpha (t{)}^H{\mathbf{S}}^H+\mathbf{n}(t)\mathbf{n}(t{)}^H]\\ & =\mathbf{SA}{\mathbf{S}}^H+{\sigma }^2\mathbf{I}\\ & ={\mathbf{R}}_S+{\sigma }^2\mathbf{I}\end{array}$$

其中，$\mathbf{A}$ 為不同信號源之間的協(xié)方差矩陣，由于不同信號源之間是相互正交的，$\mathbf{A}$ 為正定對角矩陣：
A = [ E [ ∣ α 1 ( t ) ∣ 2 ] … … … … E [ ∣ α 2 ( t ) ∣ 2 ] … … … … … … … … … E [ ∣ α M ( t ) ∣ 2 ] ] \mathbf A = \left[\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]………?…E[∣α2?(t)∣2]……?…………?………E[∣αM?(t)∣2]? ?

由于信號協(xié)方差矩陣 ${\mathbf{R}}_S$ 的規(guī)模為 $N\times N$，秩為 $M$， ${\mathbf{R}}_S$ 存在 $N?M$ 個特征值為 0 的特征向量，令這種特征向量為 ${\mathbf{q}}_m$，則
$${\mathbf{R}}_S{\mathbf{q}}_m=0$$
$$?\mathbf{SA}{\mathbf{S}}^H{\mathbf{q}}_m=0$$
$$?{\mathbf{q}}_m^H\mathbf{SA}{\mathbf{S}}^H{\mathbf{q}}_m=0$$
$$?{\mathbf{S}}^H{\mathbf{q}}_m=0$$

上述推論說明，${\mathbf{R}}_S$ 的特征值為 0 時對應(yīng)的特征向量 ${\mathbf{q}}_m$ 與信號源對應(yīng)的 $M$ 個導(dǎo)向矢量均正交。

令 ${\mathbf{R}}_S$ 的 $N?M$ 個特征值為 0 時對應(yīng)的特征向量構(gòu)成矩陣 $\mathbf Q_n $，其規(guī)模為 $N\times (N?M)$，則
$${\mathbf{S}}^H{\mathbf{Q}}_n=0$$

則 MUSIC 算法的譜估計公式為
$$P_{\mathrm{MUSIC}}(?)=\frac{1}{||{\mathbf{Q}}_n^H\mathbf{s}(?)|{|}^2}$$

當(dāng)上式中的 $?$ 與信號源方向相同時，分母為零，此時 MUSIC 譜估計為無窮大。因此，MUSIC 譜估計的尖峰數(shù)目與信源數(shù)目相同，尖峰對應(yīng)的方向即為信號源的方向。

如何根據(jù)陣列輸出信號 x \mathbf x x 計算 Q n \mathbf Q_n Qn? ？

通過記錄多組陣列輸出信號快拍，可以計算出輸出信號協(xié)方差矩陣的近似值
$$\mathbf{R}=\frac{1}{K}\sum _{k=1}^K{\mathbf{x}}_k{\mathbf{x}}_k^H$$

那么，如何根據(jù)輸出信號的協(xié)方差矩陣 $\mathbf{R}$ 估計出信號協(xié)方差矩陣 ${\mathbf{R}}_S$ 對應(yīng)的特征值為 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{array}{rl}?\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{array}$$

因此，信號協(xié)方差矩陣 ${\mathbf{R}}_S$ 的特征值 ${\lambda }_m$ 對應(yīng)的特征向量與輸出信號協(xié)方差矩陣 $\mathbf{R}$ 的特征值 ${\lambda }_m+{\sigma }^2$ 對應(yīng)的特征向量相同。

因此，$\mathbf{R}$ 的特征分解可表示為
R = Q ( Λ + σ 2 I ) Q H = Q [ λ 1 + σ 2 0 … 0 0 … 0 0 λ 2 + σ 2 … 0 0 … 0 … … … … … … … 0 0 … λ M + σ 2 0 … 0 0 0 … 0 σ 2 … 0 … … … … … … … 0 0 … 0 0 … σ 2 ] Q H \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(Λ+σ2I)QH=Q ??λ1?+σ20…00…0?0λ2?+σ2…00…0?…………………?00…λM?+σ20…0?00…0σ2…0?…………………?00…00…σ2? ?QH?

上式表明，將輸出信號矩陣 $\mathbf{R}$ 進(jìn)行特征分解，得到的 $N?M$ 個較小且相等的特征值對應(yīng)的特征向量即可構(gòu)成 ${\mathbf{Q}}_n$。

MATLAB 仿真

clc; clear; close all;

%% 參數(shù)設(shè)置
%%% 工作頻率
c = 3e8;
freq = 10e9;
lambda = c/freq;    % 波長
k = 2*pi/lambda;    % 波數(shù)
%%% 陣列參數(shù)
N = 10;                 % 陣元數(shù)量
d = 0.5*lambda;         % 陣元間隔 
z = (0:d:(N-1)*d)';     % 陣元坐標(biāo)分布
%%% 信號源參數(shù)
phi = [-10, -30, 60]'*pi/180;   % 來波方向
M = length(phi);                % 信號源數(shù)目
%%% 仿真參數(shù)
SNR = 10;             % 信噪比(dB)
K = 1000;     % 采樣點(diǎn)數(shù)

%% 陣列接收信號仿真模擬
S = exp(1j*k*z*sin(phi'));          % 流型矩陣
Alpha = randn(M, K);         % 輸入信號
X = S*Alpha;                        % 陣列接收信號
X1 = awgn(X, SNR, 'measured');      % 加載高斯白噪聲

%% MUSIC 算法
%%% 陣列接收信號的協(xié)方差矩陣的特征分解
R = X1*X1'/K;    % 陣列接收信號的協(xié)方差矩陣
[EV, D] = eig(R);       % 特征值分解
EVA = diag(D);          % 提取特征值
[EVA, I] = sort(EVA, 'descend');   % 降序排序
Q = EV(:, I);           % 特征向量構(gòu)成的矩陣
Q_n = Q(:, M+1:N);      % 噪聲子空間
%%% 計算MUSIC譜估計函數(shù)
phi_list = linspace(-pi/2, pi/2, 200)';
S1 = exp(1j*k*z*sin(phi_list'));    % 不同方向?qū)?yīng)的流型矢量構(gòu)成矩陣
P_MUSIC = 1./sum(abs(Q_n'*S1).^2);     % MUSIC 譜估計公式
%%% 轉(zhuǎn)換為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('信號源估計方向?yàn)椋?);
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

參考文獻(xiàn)

[1] 王永良. 空間譜估計理論與算法[M]. 清華大學(xué)出版社, 2004.
[2] 張小飛, 陳華偉, 仇小鋒. 陣列信號處理及MATLAB實(shí)現(xiàn)[M]. 電子工業(yè)出版社, 2015.文章來源地址http://www.zghlxwxcb.cn/news/detail-436088.html

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