雷达信号到达角估计算法MUSIC

kissgoodbye2012

已于 2022-07-13 17:59:12 修改

阅读量1.8k

点赞数 2

分类专栏： 005雷达电磁波信号文章标签：信号到达角 MUSIC

于 2022-07-08 17:33:56 首次发布

本文链接：https://blog.csdn.net/kissgoodbye2012/article/details/125682290

版权

005雷达电磁波信号专栏收录该内容

21 篇文章 13 订阅

订阅专栏

MUSIC算法，在1979年由美国人RALPH O. SCHMIDT提出

参考文献：Multiple Emitter Location and Signal Parameter Estimation

圆阵MUSIC算法带Matlab代码

一、均匀线阵的MUSIC算法

1.1 算法步骤

步骤1：获得协方差矩阵Rx

步骤2：进行特征值分解

步骤3：按照特征值排序，获得噪声矩阵En

步骤4：计算谱函数，寻找峰值得波达方向的估计值

1.2 算法仿真（Matlab）

clc;clear; close all;
%%%%%%%% MUSIC for 标准线形阵列天线%%%%%%%%

%% ======================产生模拟信号========================
derad = pi/180;      % 角度->弧度
N = 8;               % 阵元个数        
M = 4;               % 信源数目
theta = [-60 -12 20 50];  % 待估计角度
snr = 10;            % 信噪比
K = 512;             % 快拍数
 
dd = 0.5;            % 阵元间距 

%各阵元相对于第一个阵元的距离
d=0:dd:(N-1)*dd;
A=exp(-1i*2*pi*d.'*sin(theta*derad));  %方向矢量

%%%%构建信号模型%%%%%
S=randn(M,K);             %信源信号，入射信号
X=A*S;                    %构造接收信号
X1=awgn(X,snr,'measured'); %将白色高斯噪声添加到信号中

%% ======================开始MUSIC算法========================
M = 4;               % 估计信源数目

%% 步骤1：获得协方差矩阵
Rxx=X1*X1'/K;

%% 步骤2：进行矩阵特征值分解
[EV,D]=eig(Rxx);                   %特征值分解
EVA=diag(D)';                      %将特征值矩阵对角线提取并转为一行

%% 步骤3：排序特征值，获得噪声矩阵
[EVA,I]=sort(EVA);                 %将特征值排序 从小到大
EV=fliplr(EV(:,I));                % 对应特征矢量排序                 
 
%% 步骤4：遍历每个角度，计算空间谱
AngleDelta=0.1;%角度间隔
PointNum=(180/AngleDelta)+1;

for iang = 1:PointNum
    angle(iang)=180/(PointNum-1)*(iang-1)-90;%角度遍历    
    phim=derad*angle(iang);%角度值转换为弧度单位
    
    a=exp(-1i*2*pi*d*sin(phim)).'; 
    En=EV(:,M+1:N);                   % 取矩阵的第M+1到N列组成噪声子空间
    Pmusic(iang)=1/(a'*En*En'*a);
end

%% 归一化处理
Pmusic=abs(Pmusic);
Pmmax=max(Pmusic)
Pmusic=10*log10(Pmusic/Pmmax); 

%% 绘制空间谱
h=plot(angle,Pmusic);
set(h,'Linewidth',2);
xlabel('入射角/(degree)');
ylabel('空间谱/(dB)');
set(gca, 'XTick',[-90:30:90]);
grid on;

%% 峰值搜索
[mm_pks,mm_locs] = findpeaks(Pmusic);
mm_locs=angle(mm_locs);
hold on
plot(mm_locs,mm_pks,'ro');

%% 差分
dPmusic=diff(Pmusic);
dPmusic=[0 dPmusic];
figure
plot(angle,dPmusic);