基于MUSIC算法的DOA估计

function index = MUSIC(X,Chirp_Num,Target_Num)
%% MUSIC算法：用于方位角的角度估计
%% BY YUXULIANG,ZJU,20211123

    d2rad = pi / 180; % pi rad = 180 ° 1°= pi/180 rad
    N = 8; % 阵元数目
    M = Target_Num; % 信源数目
    lambda = physconst('lightspeed') / 77e9; % 波长
    D = 0.5 * lambda; % 阵元间距
    DD = 0:D:(N-1)*D; % 信号源到每个阵列的距离，维度为1×8
    % 计算协方差矩阵
    RX = X * X' / Chirp_Num;
    % 特征值分解
    [EV,D] = eig(RX); % 特征值分解
    EVA = diag(D)'; % 特征值对角线提取并转为一行
    [~,I] = sort(EVA); % 特征值排序从小到大
    EV = fliplr(EV(:,I)); % 特征矢量排序
    En = EV(:,M+1:N); % 噪声子空间
    for angle_index = 1:181
        angle(angle_index) = angle_index - 91; % 角度区间为-90到90
        phim = d2rad * angle(angle_index);
        a = exp(-1j * 2 * pi * DD * sin(phim)/lambda) .'; % 维度为8×1
        Pmusic(angle_index) = 1 / (a' * En * En' * a);
    end
    Pmusic = abs(Pmusic);
    index = find(Pmusic == max(Pmusic));
end

function [index1,index2,Pmusic] = MUSIC_2D(X1,CHIRP_NUM,TARGET_NUM)
%% MUSIC算法：用于方位角和俯仰角的角度估计

    d2rad = pi / 180; % pi rad = 180 ° 1°= pi/180 rad
    N = 8; % 阵元数目
    M = TARGET_NUM; % 信源数目
    lambda = physconst('lightspeed') / 77e9; % 波长
    D = 0.5 * lambda; % 阵元间距
    D1 = 0:D:(N-1)*D; % TX0和TX2的线阵间距
    D2 = 2 * D : D : 5 * D; % TX1的线阵间距
    D_AZI = [D1,D2]; % TX0 TX2 和 TX1的水平间距
    D_ELE = [zeros(1,8) D*ones(1,4)]; % TX0/TX2和TX1的纵向间距

    %% 计算协方差矩阵
    RX = X1 * X1' / CHIRP_NUM;
    [EV,D] = eig(RX); % 特征值分解
    EVA = diag(D);% 特征值对角线提取并转为一行
    [~,I] = sort(EVA);% 特征值排序从小到大，eig默认有排序
    EV = fliplr(EV(:,I));% 特征矢量排序
    EN = EV(:,M+1:N); % 噪声子空间

    %% 遍历每个方位角和俯仰角 计算空间谱
    for ele = 1:181 % 俯仰角遍历
        phim(ele) = ele - 91; 
        phim1 = d2rad * phim(ele);
        for azi = 1:181 % 方位角遍历
            theta(azi) = azi - 91;
            theta1 = d2rad * theta(azi);
            a = exp(-1j * 2 * pi * (D_AZI * cos(phim1) * sin(theta1) + D_ELE * sin(phim1)) / lambda).'; % 构建导向矢量
            Pmusic(azi,ele) = 1 / (a' * EN *EN' * a); % 空间谱函数
        end
    end
    Pmusic = abs(Pmusic);
    [index1,index2] = find(Pmusic == max(max(Pmusic))); % 找到空间谱谱峰并返回俯仰角和方位角
end