PUMA：DOA估计模式的改进实现（Matlab代码实现）
👨🎓个人主页
💥💥💞💞欢迎来到本博客❤️❤️💥💥
🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。
⛳️座右铭：行百里者，半于九十。
📋📋📋本文目录如下：🎁🎁🎁
目录
💥1 概述
📚2 运行结果
🎉3 参考文献
🌈4 Matlab代码、文章下载
💥1 概述
方向估计方法（MODE）具有诸多吸引人的优势，例如渐近效率高、计算复杂度低，以及在处理相干信号时表现出色，这些优势是传统基于子空间的方法所不具备的。然而，MODE在根多项式系数的对称性上引入了额外的假设和约束，这可能会在低信噪比/小样本量的情况下导致性能严重下降，因为任何估计误差都会因对称性而被放大两倍。此外，MODE的标准实现没有用于更新其估计值的闭式解。
在本文中，我们证明了MODE的优化问题等价于主特征向量用于模态分析（PUMA）算法的优化问题。我们指出，PUMA具有闭式解，且不依赖于系数的任何额外假设和约束，因此是比MODE更好的替代方案，用于最小化相同的代价函数。我们通过大量的仿真实验结果来支持我们的观点。
到达方向（DOA）估计是一个基本的阵列处理问题，具有广泛的应用，例如雷达[1]、声呐[2]和无线通信[3]、[4]。在过去四十多年中，DOA估计已经得到了深入研究，产生了许多高效且准确的算法[5]–[22]。其中，最大似然（ML）方法[5]能够提供高效的DOA估计，但代价是巨大的计算复杂度，这主要是由于多维搜索造成的。尽管可以应用诸如牛顿法等优化方法来避免搜索步骤，但由于ML目标函数是非凸的，因此无法保证全局收敛。基于子空间的DOA估计算法能够在分辨率能力和计算复杂度之间提供良好的折衷，并且得到了广泛研究。特别是，作为子空间方法的代表，多信号分类（MUSIC）[6]和通过旋转不变技术估计信号参数（ESPRIT）[7]由于其简单性和高精度而受到广泛关注，随后也开发了大量它们的变体[8]–[18]。例如，在[14]–[18]中设计了计算效率更高的MUSIC算法的改进版本。然而，当出现相干/高相关信号时，这些算法的性能会下降，因此需要去相干技术，例如空间平滑（SS）[19]或前向-后向空间平滑（FBSS）[20]，但代价是损失阵列孔径。特别是，提出了一种实值版本的根MUSIC，称为酉根MUSIC[18]。该算法利用均匀线性阵列（ULA）的中心厄米特性质，将复值数据转换到实空间。这一过程相当于对样本协方差矩阵进行一步前向-后向平滑，因此它最多可以处理两个相干信号。然而，当存在超过两个相干信号时，酉根MUSIC算法无法正常工作。
与MUSIC和ESPRIT相比，方向估计方法（MODE）[21]–[22]引起了广泛关注，因为它在某些方面优于MUSIC和ESPRIT。这些优势包括以下几点：
它的表现类似于ML，但不需要计算密集型的搜索过程，并且在大样本情况下被认为是近似高效的。
它能够处理相干信号。
它本质上不是一种迭代方法，因此不存在收敛问题。
此外，在ULA设置中，MODE可以通过多项式根求解高效地实现。这些特性使得MODE成为DOA估计的一个有力候选方法。然而，为了高效地最小化MODE代价函数，Stoica和Sharman[21]假设多项式系数是对称的，这一假设对于求解根多项式是必要但不充分的，然后他们引入了额外的约束，即固定第一个多项式系数的实部或虚部为1，以避免系数的非唯一性。需要注意的是，这种假设并不完全正确，可能会导致性能损失。在实际应用中，由于噪声干扰，特别是在低信噪比（SNR）或小样本量的情况下，MODE的性能对对称性假设非常敏感，因为系数主要依赖于其一半的元素。当估计值不准确时，在使用共轭对称性补偿另一半后，误差会被放大两倍，从而导致性能不佳。根据我们的经验，对于奇数个源的情况，性能下降非常频繁，甚至在单源情况下也是如此。此外，MODE没有闭式解。
为了弥补MODE在低信噪比区域的性能损失，Gershman和Stoica提出了MODEX（带额外根的MODE）算法。MODEX分别假设源的数量为K和P（P>K），运行两次MODE，生成（P+K）个DOA候选值，然后使用确定性/随机ML代价函数帮助从这些候选值中选择K个作为最终的DOA估计值。如果K是奇数，我们适当选择P=K+1（P为偶数），MODEX将显著改善MODE的阈值性能，其中性能提升主要依赖于P为偶数时的MODE。最近，我们还采用了类似的想法，提出了一种增强型主特征向量用于模态分析（EPUMA）技术来执行DOA估计。由于我们对多项式系数向量没有任何假设，EPUMA在K为奇数的情况下不会出现性能损失。与MODEX不同，EPUMA无论K为何值都能可靠地运行。更详细的比较将在后续的第III节和第IV节中提供。
在本文中，我们专注于介绍MODE的高效实现，即PUMA，并在不引入任何额外假设或约束多项式系数的情况下，建立MODE的方差表达式。我们证明了基于不同理论推导的PUMA和MODE具有等价的代价函数。我们希望对这两种算法及其各自的变体（即EPUMA和MODEX）进行深入比较，以展示PUMA优于MODE，并且应该成为最小化相同代价函数的更好替代方法。仿真结果表明，MODE的性能对源的数量非常敏感。具体来说，MODE在源的数量为偶数时表现良好，但在奇数情况下则不然，甚至在只有一个源时也无法正常工作。然而，PUMA并不存在这样的问题，它不仅可以提供可靠的DOA估计，在奇数源场景中远远优于MODE，而且比MODE更容易实现，计算上也更具吸引力。
V. 结论
在本文中，我们证明了MODE代价函数与最近提出的PUMA方法是等价的。与原始的MODE方法不同，PUMA不需要对多项式系数进行任何假设和约束，并且它具有多项式系数更新的闭式表达式。因此，PUMA比原始方法更容易实现。数值结果展示了PUMA及其变体的优势。
从仿真中可以观察到一些有趣的结果：当源的数量为偶数时，MODE的表现略优于PUMA。然而，当源的数量为奇数时，MODE的性能会严重下降。根据我们的推测，MODE的不稳定性能不仅出现在我们的参数设置中，还可能出现在其他一些设置中——源的数量为偶数（或奇数）有助于（或破坏）MODE的性能。然而，PUMA并不存在这样的问题。由于PUMA和MODE最小化的是同一个无约束优化问题，我们建议使用PUMA而不是MODE来实现稳定性。此外，MODE在某些情况下仍然优于PUMA，我们将致力于挖掘它们的优势并将它们结合起来，以进一步提高在低信噪比和小样本场景下的精度。
📚2 运行结果
部分代码：
M = 10;
N = 50;
DOA = [-5, 2, 12];
K = length(DOA);
SNR = linspace(-10,6,11);
nT = 100;
for iS = 1:length(SNR)
snr = SNR(iS);
for iT = 1:nT
if rem(iT,nT/2) == 0
fprintf( 'n = %d, Trials = %d, total = %d\n',...
iS, iT, (iS-1)*nT+iT );
end
x = StatSigGenerate(M, N, DOA, snr*ones(1,K));
doa1(:,iT) = EPUMA(x, K, K, 3);
doa2(:,iT) = EPUMA(x, K, K+1, 3);
doa3(:,iT) = rMUSIC(x, K, 'FBSS', 2);
[doa4(:,iT),doa5(:,iT)] = MODEX(x, K);
end
RMSE1(iS) = rmse(doa1, DOA);
RMSE2(iS) = rmse(doa2, DOA);
RMSE3(iS) = rmse(doa3, DOA);
RMSE4(iS) = rmse(doa4, DOA);
RMSE5(iS) = rmse(doa5, DOA);
[x, A, R_idl, Rs] = StatSigGenerate(M, N, DOA, snr*ones(1,K));
CRB(iS) = crbdet_w(A,R_idl,Rs,DOA,N,1)*(180/pi)^2;
end
mz = 8;
lw = 2;
figure
semilogy(SNR, RMSE1.^0.5, '-p', 'markersize', mz, 'linewidth', 2); hold on;
semilogy(SNR, RMSE2.^0.5, '-o', 'markersize', mz, 'linewidth', 2);
semilogy(SNR, RMSE3.^0.5, '->', 'markersize', mz, 'linewidth', 2);
semilogy(SNR, RMSE4.^0.5, '-*', 'markersize', mz, 'linewidth', 2)
semilogy(SNR, RMSE5.^0.5, '-*', 'markersize', mz, 'linewidth', 2)
semilogy(SNR, CRB.^0.5, 'k', 'linewidth', 2)
xlabel('SNR (dB)'); ylabel('RMSE (degree)');
legend('PUMA', 'EPUMA', 'root-MUSIC', 'MODEX', 'MODE', 'CRB');
M = 10;
N = 50;
DOA = [-5, 2, 12];
K = length(DOA);
SNR = linspace(-10,6,11);
nT = 100;
for iS = 1:length(SNR)
snr = SNR(iS);
for iT = 1:nT
if rem(iT,nT/2) == 0
fprintf( 'n = %d, Trials = %d, total = %d\n',...
iS, iT, (iS-1)*nT+iT );
end
x = StatSigGenerate(M, N, DOA, snr*ones(1,K));
doa1(:,iT) = EPUMA(x, K, K, 3);
doa2(:,iT) = EPUMA(x, K, K+1, 3);
doa3(:,iT) = rMUSIC(x, K, 'FBSS', 2);
[doa4(:,iT),doa5(:,iT)] = MODEX(x, K);
end
RMSE1(iS) = rmse(doa1, DOA);
RMSE2(iS) = rmse(doa2, DOA);
RMSE3(iS) = rmse(doa3, DOA);
RMSE4(iS) = rmse(doa4, DOA);
RMSE5(iS) = rmse(doa5, DOA);
[x, A, R_idl, Rs] = StatSigGenerate(M, N, DOA, snr*ones(1,K));
CRB(iS) = crbdet_w(A,R_idl,Rs,DOA,N,1)*(180/pi)^2;
end
mz = 8;
lw = 2;
figure
semilogy(SNR, RMSE1.^0.5, '-p', 'markersize', mz, 'linewidth', 2); hold on;
semilogy(SNR, RMSE2.^0.5, '-o', 'markersize', mz, 'linewidth', 2);
semilogy(SNR, RMSE3.^0.5, '->', 'markersize', mz, 'linewidth', 2);
semilogy(SNR, RMSE4.^0.5, '-*', 'markersize', mz, 'linewidth', 2)
semilogy(SNR, RMSE5.^0.5, '-*', 'markersize', mz, 'linewidth', 2)
semilogy(SNR, CRB.^0.5, 'k', 'linewidth', 2)
xlabel('SNR (dB)'); ylabel('RMSE (degree)');
legend('PUMA', 'EPUMA', 'root-MUSIC', 'MODEX', 'MODE', 'CRB');
🎉3参考文献
部分理论来源于网络，如有侵权请联系删除。
[1]C. Qian, L. Huang, M. Cao, H. C. So and J. Xie, "PUMA: An improved realization of MODE for DOA estimation," IEEE Transactions on Aerospace and Electronic Systems, vol. 53, no. 5, pp. 2128-2139, 2017.