关于MUSIC算法的知识点和MATLAB程序详细注释—

[该文章是自己没事学习的一些笔记，其内容可能会包含参考书和一些论文上的话语，这只是一个学习过程的笔记。文末会备注参考文献]如有侵权，可让我删除。

MUSIC算法是经典的波达方向（DOA）估计算法

1.MUSIC是啥？中文：多重信号分类（Multiple Signal Classification）

2.基本思想：特征分解任意阵列输出数据的协方差矩阵，从而得到与信号分类相对应的信号子空间和与信号分量相正交的噪声子空间，接着构造两个子空间的正交空间谱函数，从而通过搜索谱峰检测出信号的波达方向。

3.基本原理

注明！！！：参考文献

4.MUSIC计算量分析

5.天线阵元数，阵元间距，快拍数和信噪比都对MUSIC算法的性能影响很大。

6.MATLAB代码部分：

//几乎每句话都给了注释，助于理解啦

运行结果：

[该文章是自己没事学习的一些笔记，其内容可能会包含参考书和一些论文上的话语，这只是一个学习过程的笔记。文末会备注参考文献]如有侵权，可让我删除。

MUSIC算法是经典的波达方向（DOA）估计算法

1.MUSIC是啥？中文：多重信号分类（Multiple Signal Classification）

2.基本思想：特征分解任意阵列输出数据的协方差矩阵，从而得到与信号分类相对应的信号子空间和与信号分量相正交的噪声子空间，接着构造两个子空间的正交空间谱函数，从而通过搜索谱峰检测出信号的波达方向。

3.基本原理

常用典型一维DOA估计数据模型

注明！！！：参考文献

闫锋刚,沈毅,刘帅,等. 高效超分辨波达方向估计算法综述[J]. 系统工程与电子技术,2015(7):1465-1475. DOI:10.3969/j.issn.1001-506X.2015.07.01.

4.MUSIC计算量分析

主要包含三个计算单元：1.利用N快拍接收数据估计协方差矩阵R；2.对R进行特征值分解（EVD）或奇异值分解（SVD）估计信号子空间矩阵和噪声子空间矩阵；3.在空间谱范围内遍历搜索谱峰极值，进而得到信号DOA。

5.天线阵元数，阵元间距，快拍数和信噪比都对MUSIC算法的性能影响很大。

6.MATLAB代码部分：

//几乎每句话都给了注释，助于理解啦

clear all %清除工作空间的所有变量，函数，和MEX文件
close all %关闭所有的Figure窗口

%%%%输入信号模拟%%%%
derad = pi/180; %角度转弧度
radeg = 180/pi; %弧度转角度
twpi = 2*pi;

%%%%输入阵列配置参数：阵元数量、阵元坐标（间距）
kelm = 8; % 阵列数量（阵元个数）
dd = 0.5; % 阵元间距d = lambda/2 （阵元间距与波长的比值）
d=0:dd:(kelm-1)*dd; % 阵元序列
iwave = 3; % 信号源数目

%构建阵元接收信号
theta = [10 30 60]; % 入射信号角度
snr = 10; % 信噪比
n = 500; % 快拍个数（采样点数）
A=exp(-j*twpi*d.'*sin(theta*derad));%构建信号导向矢量矩阵
S=randn(iwave,n); %（randn生成标准正态分布随机数 3*n的） (空间信号源矩阵) %信源信号，射入信号
X=A*S; %构造接收信号
X1=awgn(X,snr,'measured'); %将高斯白噪声加入到信号中
Rxx=X1*X1'/n; %计算协方差矩阵【用的最大似然估计公式】

InvS=inv(Rxx);
[EV,D]=eig(Rxx);%特征值分解计算Rxx的特征值对应的对角阵D和特征向量构成的矩阵EV
EVA=diag(D)';%将特征值矩阵对角线提取并转为一行 diag抽取矩阵对角线元素
[EVA,I]=sort(EVA);%将特征值排序从小到大
EVA=fliplr(EVA);%特征值左右翻转，从小到大排序
EV=fliplr(EV(:,I));%对特征向量排序

% 构造MUSIC谱函数
for iang = 1:361 %遍历每个角度，计算空间谱
angle(iang)=(iang-181)/2; %【范围变成-90-90度】这是线阵的测量范围
phim=derad*angle(iang); %角度变为弧度
a=exp(-j*twpi*d*sin(phim)).';
L=iwave;
En=EV(:,L+1:kelm); % 取矩阵的第L+1到N列组成的噪声子空间
SP(iang)=(a'*a)/(a'*En*En'*a);
end

%作图
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

运行结果：