[算法]music测向算法matlab仿真

clear all
close all

%%%%输入信号模拟%%%%  
derad = pi/180;         %角度转弧度
theta = [-30 0 30 60];  % 入射信号角度
snr = 10;               % 信噪比
n = 512;                % 快拍个数（采样点数）
iwave = 4;              % 信号源数目

%%%%输入阵列配置参数：阵元数量、阵元坐标（间距）
kelm = 8;               % 阵列数量（阵元个数）
dd = 0.5;               % 阵元间距d = lambda/2  （阵元间距与波长的比值）
d=0:dd:(kelm-1)*dd;     % 阵元序列 % 0 0.5000 1.0000 1.5000 2.0000 2.5000 3.0000 3.5000

%构建阵元接收信号
A=exp(-j*2*pi*d.'*sin(theta*derad));%构建信号导向矢量矩阵 
S=randn(iwave,n);       %（randn生成标准正态分布随机数 3*n的）  (空间信号源矩阵)   %信源信号，射入信号
X=A*S;                  %构造接收信号
X1=awgn(X,snr,'measured');%将高斯白噪声加入到信号中

Rxx=X1*X1'/n;           %计算协方差矩阵  【用的最大似然估计公式】
[EV,D]=eig(Rxx);        %特征值分解 计算Rxx的特征值对应的对角阵D和特征向量构成的矩阵EV
EVA=diag(D)';           %将特征值矩阵对角线提取并转为一行   diag抽取矩阵对角线元素
[EVA,I]=sort(EVA);      %将特征值排序 从小到大
EV=fliplr(EV(:,I));     %对特征向量排序

% MUSIC
for iang = 1:361                     %遍历每个角度，计算空间谱
        angle(iang)=(iang-181)/2;    %【范围变成-90-90度】这是线阵的测量范围 %-90至90分成360等分,每0.5一份
        phim=derad*angle(iang);      %角度变为弧度
        a=exp(-j*2*pi*d*sin(phim)).';%方向矢量
        L=iwave;    
        En=EV(:,L+1:kelm);           % 取矩阵的第L+1到N列组成的噪声子空间
        pmuaic(iang)=1/(a'*En*En'*a);%MUSIC算法步骤归纳，计算谱函数
end
   
%作图
pmuaic=abs(pmuaic);                         %绝对值处理函数   即数值的绝对值和复数的幅值
pmax=max(pmuaic);
pmuaic=10*log10(pmuaic/pmax);              %归一化处理
h=plot(angle,pmuaic);
set(h,'Linewidth',3)
xlabel('入射角 (degree)')
ylabel('空间谱 (dB)')
axis([-90 90 -50 50])
set(gca, 'XTick',[-90:30:90])
grid on