MUSIC算法估计DOA代码

%%MUSIC算法
clc
clear 
close all

%% 基本参数设置
M = 10;
N = 100;
theta = [10 20 30];
K = length(theta);
d = 0.5;
dd = 0:M-1;
SNR = 20;

%% 构造信号接收模型
s = randn(K,N);
a = exp(1i*2*pi*d*dd'*sind(theta));
x = a*s;
x = awgn(x,SNR,'measured');%  measured  表示使用实际测量的信噪比

%% 计算噪声子空间
R = x*x'/N;
[V,VA] = eig(R);%得到特征向量以及特征值对角矩阵
VA = diag(VA);
[~,I] = sort(VA);
V = fliplr(V(:,I));%按照特征值的降序排列对应特征向量
N_sub = V(:,K+1:M);


%% 空间谱搜索
for i = 1:181
    angle(i) = i-91;
    A = exp(1i*pi*dd'*sind(angle(i)));
    P(i) = 1/(A'*N_sub*N_sub'*A);%(A'* A)结果为M，归一化后对整个估计结果不造成影响
end
P = abs(P);
P_max = max(P);
P = 10*log10(P/P_max); %归一化操作并转化为dB
%% 绘图
figure
h_1 = plot(angle,P,'b','LineWidth',2);
axis([-90 90 -60 0])
set(gca, 'XTick',-90:30:90)
grid on
xlabel('angle(degree)')
ylabel('magnitude spectr(dB)')
legend('MUSIC')