空间谱估计MUSIC算法(含matlab仿真)

最新推荐文章于 2023-09-19 14:02:23 发布

学习小混子

最新推荐文章于 2023-09-19 14:02:23 发布

阅读量2.7k

点赞数 8

分类专栏：空间谱估计理论与算法读书笔记文章标签：算法 matlab 信号处理

本文链接：https://blog.csdn.net/weixin_48627103/article/details/130782840

版权

空间谱估计理论与算法读书笔记专栏收录该内容

1 篇文章 0 订阅

订阅专栏

经典MUSIC算法，在窄带远场假设的情况下信号的DOA数学模型为

$\boldsymbol{X}(\boldsymbol{t})=\boldsymbol{A}(\theta)\boldsymbol{s}(\boldsymbol{t})+\boldsymbol{N}(\boldsymbol{t})$

则阵列数据的协方差矩阵为

$\begin{aligned}\mathbf{R}&=E\big[\mathbf{XX}^{\mathbb{H}}\big]\\ &=\mathbf{AE}\big[\mathbf{SS}^{\mathbb{H}}\big]\mathbf{A}^{\mathbb{H}}+\boldsymbol{R}_{\mathrm{N}}\\ &=\mathbf{AR}_{\mathbb{S}}\mathbf{A}^{\mathbb{H}}+\sigma^2\mathbf{I}\end{aligned}$

式中 $\boldsymbol{R}_{\mathrm{s}}$ ， $\boldsymbol{R}_{\mathrm{N}}$ 分别为信号协方差矩阵和噪声协方差矩阵，对于理想空间的白噪声且功率为 $\sigma^2$ 。其中 $\mathbf{A}=[\boldsymbol{a}(\theta_1),\boldsymbol{a}(\theta_2),\cdots,\boldsymbol{a}(\theta_K)]\in\mathbb{C}_{M\times K}$ ，K为信号源个数，M为阵列个数。 $\mathbf{A}$ 为列满秩（秩为K）， $\boldsymbol{R}_{\mathrm{s}}$ 为共轭对称矩阵。故 $\mathbf{AR}_{\mathbb{S}}\mathbf{A}^{\mathbb{H}}$ 为M*M维的非奇异Hermitian阵，秩为K。将上述 $\mathbf{R}$ 进行特征分解得

$\boldsymbol{R}=\boldsymbol{U}\boldsymbol{\Sigma}U^{\mathrm{H}}$

式中，U为特征矢量矩阵，其中由特征值组成的对角阵Σ如下：

$\boldsymbol\Sigma=\begin{bmatrix}\lambda_1\\ &\lambda_z\\ &&\ddots\\ &&&\lambda_M\end{bmatrix}$

上式中的特征值满足如下关系： $\lambda_1\geqslant\lambda_2\geqslant\cdots\geqslant\lambda_N\geqslant\lambda_{N+1}=\cdots=\lambda_M=\sigma^2$

进一步 $\mathbf{R}$ 可以特征分解为

$\begin{gathered} \boldsymbol{R}=\sum_{i=1}^{N}\lambda_{i}\boldsymbol{e}_{i}\boldsymbol{e}_{i}^{\mathrm{H}}+\sum_{j=N+1}^{M}\lambda_{j}\boldsymbol{e}_{j}\boldsymbol{e}_{j}\boldsymbol{e}_{j}^{\mathrm{H}} \\ =\begin{bmatrix}\boldsymbol{U}_{\mathrm{S}}&\boldsymbol{U}_{\mathrm{N}}\end{bmatrix}\boldsymbol{\Sigma}\begin{bmatrix}\boldsymbol{U}_{\mathrm{S}}&\boldsymbol{U}_{\mathrm{N}}\end{bmatrix}^{\mathrm{H}} \\ =\boldsymbol{U}_{\mathrm{S}}\boldsymbol{\Sigma}_{\mathrm{S}}\boldsymbol{U}_{\mathrm{S}}^{\mathrm{H}}+\boldsymbol{U}_{\mathrm{N}}\boldsymbol{\Sigma}_{\mathrm{N}}\boldsymbol{U}_{\mathrm{N}}^{\mathrm{H}} \end{gathered}$

式中， $\boldsymbol{\Sigma}_{\mathrm{s}}$ 为大特征值组成的对角阵， $\boldsymbol{\Sigma}_{\mathrm{N}}$ 为小特征值组成的对角阵。 $\boldsymbol{U}_{\boldsymbol{S}}=\begin{bmatrix}\boldsymbol{e}_{1}&{\boldsymbol{e}_{2}}&{\cdots}&{\boldsymbol{e}_{N}}\end{bmatrix}$ 是大特征值对应的信号子空间； $\boldsymbol{U}_{N}=\begin{bmatrix}\boldsymbol{e}_{N+1}&\boldsymbol{e}_{N+2}&\cdots&\boldsymbol{e}_{M}\end{bmatrix}$ 是小特征值对应的信号子空间。 $\boldsymbol{U}_{\mathrm{s}}$ ， $\boldsymbol{U}_N$ 是非奇异Hermitian阵不同特征值对应的特征矩阵，故两个子空间是正交的。

因此信号子空间中的导向矢量也是与噪声子空间正交。

$\pmb{a}^{\mathrm{H}}(\theta)\pmb{U}_{\mathrm{N}}=0$

经典MUSIC算法正是基于上述这个这个性质提出来的，但是考虑到实际接收数据矩阵式有限长的，即数据协方差矩阵的最大似然估计为

$\hat{\boldsymbol{R}}=\frac{1}{L}\sum_{i=1}^{L}\boldsymbol{X}\boldsymbol{X}^{\mathrm{H}}$

对 $\hat{\boldsymbol{R}}$ 进行特征分解可以计算得到噪声子空间特征矢量矩阵 $\hat{\boldsymbol{U}}_{\mathrm{N}}$ 。由于噪声的存在， $\pmb{a}^{\mathrm{H}}(\theta)$ 与 $\hat{\pmb{U}}_{\mathrm{N}}$ 并不能完全正交，式 $\pmb{a}^{\mathrm{H}}(\theta)\pmb{U}_{\mathrm{N}}=0$ 并不成立。因此实际上求DOA是以最小优化搜索实现的，即

$\theta_{\mathrm{MUSIC}}=\operatorname{arg}\min\boldsymbol{a}^{\mathrm{H}}(\theta)\hat{\boldsymbol{U}}_{\mathrm{N}}\hat{\boldsymbol{U}}_{\mathrm{N}}^{\mathrm{H}}\boldsymbol{a}(\theta)$

所以，MUSIC算法的谱估计公式为

$\begin{gathered} P_{\mathrm{MUSIC}} =\frac{1}{\boldsymbol{a}^{\mathrm{H}}(\theta)\hat{\boldsymbol{U}}_{\mathrm{N}}\hat{\boldsymbol{U}}_{\mathrm{N}}^{\mathrm{H}}\boldsymbol{a}(\theta)} \end{gathered}$

下面利用matlab对MUSIC算法进行仿真验证，针对均匀线阵。

参数设置：阵元数8个，信号源3个，入射角度分别为-20，0，30，快拍数为500，讨论不同信噪比、阵元数以及快拍数变化下的MUSIC算法仿真结果。

SNR表示信噪比，snap表示快拍数，M表示阵元数。可以观察到在-20°，0°，30°出现峰值。

仿真代码：

%% 仿真MUSIC算法测向
% 参数：
% 阵元数：8
% 入射信号角度：[-20 0 30]
% 信噪比：10
% 信源数：3
% 快拍数：500
clc
clear
close all
M = [8,32,128];
c = 3e8;
lamad = 0.1;
d = 0.5*lamad;
fc = c/lamad;
fs = 500;
snap = [500,250,50];
dt = 1/fs;
% t = (1:snap(1)).*dt;
SNR = [1/db2mag(10),1,1/db2mag(-10)];
for i = 1:3
    %% 产生三个信号源
    
    t = (1:snap(1)).*dt;
    s1 = exp(1j*2*pi*(fc.*t+0.5*30*t.^2));
    s2 =exp(1j*2*pi*(fc.*t+0.5*20*t.^2));
    s3 = exp(1j*2*pi*(fc.*t+0.5*10*t.^2));
    St = [s1;s2;s3];
    %% 对应的导向矢量
    n = 0:M(i)-1;
    theta1 = -20;
    theta2 = -0;
    theta3 = 30;
    a1 = exp(-1j*pi.*n*sin(theta1/180*pi))';
    a2 = exp(-1j*pi.*n*sin(theta2/180*pi))';
    a3 = exp(-1j*pi.*n*sin(theta3/180*pi))';
    A = [a1,a2,a3];
    %% DOA估计
    Nt = zeros(M(i),length(t));
    for k=1:M(i)
        noise = SNR(1).*randn(1,length(t));
        Nt(k,:) =noise;
    end
    Xt = A*St+Nt;
    R = Xt*Xt'/length(t);
    [U,Sigema] = eig(R);
    ac = 0;
    for theta = -60:1:60
        ac = ac+1;
        a_theta = exp(-1j*pi.*n.*sin(theta./180*pi))';
        Pmusic = 1./(a_theta'*(U(:,1:5)*U(:,1:5)')*a_theta);
        amplitude(ac) = abs(Pmusic);
    end
    f1 = figure(1);
    figure(f1)
    plot(-60:1:60,pow2db(amplitude/max(amplitude)))
    hold on
end
figure(f1)
xlabel('角度(°)');
ylabel('归一化幅度(dB)');
legend('M=8','M=32','M = 128')