MUSIC做DOA估计问题是从谱估计思路出发,对接收信号做自相关处理得到R,接着对R做奇异值分解根据目标个数得到信号子空间和噪声子空间,接着利用分解的噪声特征向量得到谱函数,对空域搜索得到不同角度值下的谱函数。
MUSIC算法解释网上解释很多,我原来奇怪的一点是为什么自相关矩阵R分解后的大特征值对应的就是信号子空间,小特征值对应的就是噪声子空间。后来在知乎上看到一个有关协方差的解释 https://zhuanlan.zhihu.com/p/338335181 帮助大家理解
最近在看压缩感知和凸优化相关算法,看到有用正则化做DOA的方法,其实就是在模型s=Ag+n中加了一个约束条件
其中L(g)为代价函数,通常代价函数为g的二范数平方值。lamda为正则化系数。不同的场景下是不同的约束条件。
下边贴一下matlab代码和仿真结果
% 作者:gxy
% 功能:分别用MUSIC算法和正则化方法进行角度估计
% 时间:2021年6月8日
clc
clear all
close all
% 回波构造
N = 8; % 阵列数
dd = 0.5; % 阵列间距与波长比值
d=0:dd:(N-1)*dd; % 阵列相对波长位置
theta = [-10 5]; % 目标角度
NumK = length(theta);
snr = 20; % 输入信噪比dB
Sample = 1; % 快拍数
A=exp(-1i*2*pi*d.'*sind(theta)); % 导向矢量
S=randn(NumK,Sample); % 目标回波
X=A*S; % 阵列接收信号
X1=awgn(X,snr,'measured'); % 加噪信号
%% MUSIC一般处理步骤
Rxx=X1*X1'/Sample; % 接收信号自相关
InvS=pinv(Rxx);
[EV,D]=eig(Rxx);
EVA=diag(D)';
[EVA,I]=sort(EVA);
EVA=fliplr(EVA);
EV=fliplr(EV(:,I));
for iang = 1:361
angle(iang)=(iang-181)/2;
a=exp(-1i*2*pi*d*sind(angle(iang))).';
L=NumK;
En=EV(:,L+1:N);
SP(iang)=(a'*a)/(a'*En*En'*a);
end
% MUSIC结果绘图
SP=abs(SP);
SPmax=max(SP);
SP=10*log10(SP/SPmax);
h=plot(angle,SP);
set(h,'Linewidth',2)
% xlabel('angle (degree)')
% ylabel('magnitude (dB)')
axis([-90 90 -60 0])
set(gca, 'XTick',[-90:30:90])
grid on
%% 用约束条件
d_lamda = 0.5;
AA=[]; %构造过完备基
theta=-90:0.5:90; % 与MUSIC一样的搜索间距及步长
for kkk= 1:length(theta)
g=exp(-1i*2*pi*[0:N-1]'*d_lamda*sind(theta(kkk)));
g = repelem(g,Sample);
AA=[AA,g];
end
% 多快拍
y = reshape(X1.',1,[]);
y = y.';
cvx_begin
variable x(length(theta));
minimize(square_pos(norm(y-AA*x,2)) + 2 *(norm(x,1))); % 约束条件
cvx_end
searching_doa = theta;
% 正则化结果绘图
hold on ;
plot(searching_doa, 10*log10(abs(x)./max(abs(x))),'Linewidth',2);
xlabel('DOA/degree');
ylabel('PowerdB');
ylim([-140 0])
legend('MUSIC' , '正则化方法')
仿真在20dB信噪比下,阵元数为8,单快拍,两目标分别为-10°和5°条件下的参数估计结果。可以看出约束正则化针对低快拍,低信噪比时能准确估计到达角。
推荐一本《凸优化》,作者Stephen Boyd,王书宁等翻译,清华大学出版社
【初步认识优化方法,等过段时间深入了解之后再分享】