【MUSIC、最大似然与克拉美-罗下界】MUSIC与ESPRIT 算法来估计到达角(AoA)，并尝试推导克拉美-罗下界(CRLB)以分析其性能研究（Matlab代码实现）

  💥💥💞💞欢迎来到本博客❤️❤️💥💥

🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。

⛳️座右铭：行百里者，半于九十。

📋📋📋本文内容如下：🎁🎁🎁

⛳️赠与读者

👨‍💻做科研，涉及到一个深在的思想系统，需要科研者逻辑缜密，踏实认真，但是不能只是努力，很多时候借力比努力更重要，然后还要有仰望星空的创新点和启发点。建议读者按目录次序逐一浏览，免得骤然跌入幽暗的迷宫找不到来时的路，它不足为你揭示全部问题的答案，但若能解答你胸中升起的一朵朵疑云，也未尝不会酿成晚霞斑斓的别一番景致，万一它给你带来了一场精神世界的苦雨，那就借机洗刷一下原来存放在那儿的“躺平”上的尘埃吧。

    或许，雨过云收，神驰的天地更清朗.......🔎🔎🔎

💥1 概述

MUSIC与ESPRIT算法估计到达角(AoA)及克拉美-罗下界(CRLB)性能分析

一、MUSIC与ESPRIT算法原理及AoA估计实现

1. MUSIC算法

MUSIC（Multiple Signal Classification）算法通过分解接收信号的协方差矩阵，构建信号子空间与噪声子空间，利用两者的正交性构造空间谱函数，通过谱峰搜索确定信号到达角（AoA）。其核心步骤如下：

• 协方差矩阵构建：计算接收信号的协方差矩阵 R=E[X(t)XH(t)]，其中 X(t) 为阵列接收信号向量。
• 特征值分解：对 R 进行特征分解，得到大特征值对应的信号子空间 Us 和小特征值对应的噪声子空间 Un。
• 谱函数构造：定义空间谱函数 PMUSIC(θ)=aH(θ)UnUnHa(θ)1，其中 a(θ) 为阵列流型向量。
• 谱峰搜索：通过搜索谱函数的峰值位置，确定信号源的AoA。

优势：高分辨率，适用于多信号源场景；劣势：计算复杂度高，需全空间搜索，对阵列校准敏感。

2. ESPRIT算法

ESPRIT（Estimation of Signal Parameters via Rotational Invariance Techniques）算法利用阵列流型的旋转不变性，通过子空间分解直接估计信号参数，避免谱峰搜索。其核心步骤如下：

• 数据矩阵构建：将阵列接收数据划分为两个子阵列 X1 和 X2。
• 子空间分解：对子阵列数据进行奇异值分解（SVD），得到信号子空间 Us。
• 旋转矩阵估计：利用子阵列间的旋转不变性关系 X2=ΦX1，估计旋转矩阵 Φ。
• 参数估计：通过 Φ 的特征值分解，确定信号源的AoA。

优势：计算复杂度低，无需全空间搜索；劣势：对阵列几何结构有一定要求，需子阵列间存在旋转不变性。

二、克拉美-罗下界（CRLB）推导与性能分析

1. CRLB定义

CRLB为无偏估计量的方差下限，表示为 Var(θ^)≥CRLB(θ)，其中 θ 为待估参数（如AoA）。CRLB由Fisher信息矩阵 I(θ) 的逆决定：

编辑

2. AoA估计的CRLB推导

假设阵列接收信号为 X(t)=A(θ)S(t)+N(t)，其中 A(θ) 为阵列流型矩阵，S(t) 为信号源信号，N(t) 为高斯白噪声。

• Fisher信息矩阵计算：
  在已知信号源数量 D 和信噪比（SNR）条件下，Fisher信息矩阵可表示为：
  

编辑

• CRLB表达式：
  对于单信号源场景，CRLB可简化为：
  

编辑

3. 性能分析

• MUSIC算法：在小采样快拍数条件下，MUSIC算法的空间谱估计方差近似CRLB；在多信号源场景下，需满足高SNR和大采样快拍数条件，方差方可接近CRLB。
• ESPRIT算法：由于避免了谱峰搜索，ESPRIT算法的计算复杂度低于MUSIC，但在低SNR或小采样快拍数条件下，其估计方差可能高于CRLB。
• 比较：MUSIC算法在理想条件下（高SNR、大采样快拍数）性能更优，但ESPRIT算法在计算效率和鲁棒性方面表现更佳。

三、算法选择与优化策略

1. 算法选择建议

• 高分辨率需求：选择MUSIC算法，适用于雷达、声纳等对精度要求极高的场景。
• 计算效率需求：选择ESPRIT算法，适用于无线通信、实时信号处理等对实时性要求较高的场景。
• 相干信号处理：结合空间平滑技术（如SS、MSS）改进MUSIC算法，或选择TLS-ESPRIT等变体算法处理相干信号。

2. 优化策略

• MUSIC算法优化：

• 采用波束空间MUSIC（B-MUSIC）算法，通过波束形成矩阵降低信噪比门限，提高分辨能力。
• 结合稀疏优化技术，提升小采样快拍数条件下的估计性能。

• ESPRIT算法优化：

• 采用总体最小二乘（TLS-ESPRIT）算法，提高对噪声和模型误差的鲁棒性。
• 推广至二维或三维DOA估计，适应复杂信号环境。

四、未来发展趋势

• MIMO系统集成：随着MIMO技术在无线通信中的广泛应用，MUSIC与ESPRIT算法将进一步优化，以适应大规模阵列和动态场景需求。
• 深度学习融合：结合深度学习技术，提升算法在非理想条件下的性能，如低SNR、多径效应等。
• 实时处理能力提升：通过硬件加速（如FPGA、GPU）和算法简化，实现算法的实时处理，满足工业、医疗等领域的需求。

📚2 运行结果

编辑

编辑

编辑

部分代码：

%% 参数设置
N               =           64;                % 基站处天线
fc              =           100e9;              %100GHz
lambda_c        =           3e8/fc;
d               =           lambda_c/2;
T               =           10;
Nit             =           1e2;
SNRdBs          =           -10:5:20;
%% 
theta   = 0.5;%pi/6
H   = genSteerVector(theta, N, d, lambda_c);
Y   = zeros(N, T);
A   = eye(N);
a1  = genPartialSteerVector(theta, N, d, lambda_c, 1);
a2  = genPartialSteerVector(theta, N, d, lambda_c, 2);
MseMUSIC    = zeros(length(SNRdBs), 1);
MseESPRIT   = zeros(length(SNRdBs), 1);
CRLB        = zeros(length(SNRdBs), 1);
D           = a1;
Cst         = D'*(eye(N)-H*inv(H'*H)*H')*D;
for i_SNR = 1:length(SNRdBs)
    fprintf('==============================================\n');
    fprintf('SNR        = %ddB\n', SNRdBs(i_SNR));
    for i_it = 1:Nit
        
        if mod(i_it,10) == 0
            fprintf('i_it/Nit   = %3d/%3d\n', i_it, Nit);
        end
        X = (sqrt(2)/2*randn(1,T)+1j*sqrt(2)/2*randn(1,T));
        for i_T = 1:T
            Y(:, i_T) = awgn(H*X(i_T), SNRdBs(i_SNR), 'measured');
        end
        
        theta_MUISC     = MUSIC(Y);
        psi = TLS_ESPRIT_Algorithm(Y, 1);
        theta__ESPRIT = log(psi)/(1j*pi); 
        MseMUSIC(i_SNR)     = MseMUSIC(i_SNR) + abs(asin(theta_MUISC) - asin(theta))^2;
        MseESPRIT(i_SNR)    = MseESPRIT(i_SNR) + abs(asin(theta__ESPRIT) - asin(theta))^2;
        
        sigma2  = 10^(-SNRdBs(i_SNR)/10);
% 这个sigma2和上面的值是渐进一致的
%         sigma2  = (norm(Y, 'fro')^2-norm(H*X, 'fro')^2)/...
%                 (size(Y, 1)*size(Y, 2));
%         X_bar   = kron(X.', eye(N));
%         y       = reshape(Y,[],1);
%         CRLB(i_SNR)     = CRLB(i_SNR) + sigma2/2./...
%                             real(-y'*X_bar*a2 + a2'*(X_bar'*X_bar)*H + a1'*(X_bar'*X_bar)*a1);
        CRLB(i_SNR)     = CRLB(i_SNR) + sigma2/2/real((Cst*(X*X')));
    end
end
MseMUSIC     = MseMUSIC/Nit;
MseESPRIT    = MseESPRIT/Nit;
CRLB         = CRLB/Nit;
% fprintf("计算得到的SNRdB = %f\n", 10*log10((norm(H*X, 'fro')^2)/(norm(Y, 'fro')^2-norm(H*X, 'fro')^2)))
%% 绘图
set(0,'defaultfigurecolor','w') 
figure; hold on; grid on; box on;
xlabel('SNR/dB');
ylabel('MSE(rad^2)');
set(gca, 'YScale', 'log');
plot(SNRdBs, MseESPRIT,'r:s', 'LineWidth', 2);
plot(SNRdBs, MseMUSIC,'g:o', 'LineWidth', 2);
plot(SNRdBs, CRLB,'m:<', 'LineWidth', 2);
l = legend({  
        'ESPRIT',...    
        'MUSIC',...
        'CRLB'},...    
        'Interpreter','latex', 'Box','off'); 
l.FontSize = 12;

🎉3 参考文献

文章中一些内容引自网络，会注明出处或引用为参考文献，难免有未尽之处，如有不妥，请随时联系删除。(文章