MUSIC算法原理分析与实现

很多人私信博主本人关于阵列信号处理的相关问题，顾在此开一个专栏，写一写个人关于阵列信号处理相关算法的一些文章，因一些内容为本人已发表和待发表论文（注：本文内容已收入知网），为了方便，其中有发杂公式的部分会采用截图的形式说明。

正文

学模型，如图2.1所示。

仿真实现代码：

% DOA estimation by MUSIC 
% Developed by Megustas (DLUT）
% Blog ：http://blog.csdn.net/megustas_jjc
clear all
close all
derad = pi/180;        % deg -> rad
radeg = 180/pi;
twpi = 2*pi;
kelm = 8;               % number
dd = 0.5;               % space 
d=0:dd:(kelm-1)*dd;     % 
iwave = 3;              % number of DOA
theta = [10 30 60];     % angle
snr = 10;               % input SNR (dB)
n = 500;                 % 
A=exp(-j*twpi*d'*sin(theta*derad));%%%% direction matrix
S=randn(iwave,n);
X=A*S;
X1=awgn(X,snr,'measured');
Rxx=X1*X1'/n;
InvS=inv(Rxx); %%%%
[EV,D]=eig(Rxx);%%%% 
EVA=diag(D)';
[EVA,I]=sort(EVA);
EVA=fliplr(EVA);
EV=fliplr(EV(:,I));

% MUSIC
for iang = 1:361
        angle(iang)=(iang-181)/2;
        phim=derad*angle(iang);
        a=exp(-j*twpi*d*sin(phim)).';
        L=iwave;    
        En=EV(:,L+1:kelm);
        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  

关于ESPRIT算法及均匀圆阵和多数据集、宽带信号处理等问题会在后续博客讲到。