%************************MUSIC算法**********************************************
%******************************************************************************
clc
clear all
close all
%function MUSI();
%===============================参数设定========================================
clc
clear all
%if nargin==0;
d=1.5;
M=8;//阵元数
lamda=3;
sam=2000;//采样点数
theta1=20;//入射波角度1
theta2=10;//入射波角度2
degrad=pi/180;
theta=[theta1,theta2];
m=[0:M-1]';
A1=exp(-2*pi*j*(d/lamda)*m*sin(theta1*degrad));
A2=exp(-2*pi*j*(d/lamda)*m*sin(theta2*degrad));
A=[A1,A2];%方向矩阵；
n=1:sam;
v1=1e+12;
v2=3e+12;
SNR=30;//信噪比
Am=10^(SNR/20);
signal=Am*[sin(v1*n);sin(v2*n)];%接收信号；
noise=randn(M,sam)+j*randn(M,sam);
A1=A*signal;//模拟天线接收到的信号
A1(3,:)=A1(3,:)/10;
X=A1+noise;//加噪声
%end
mm=size(signal);
nn=size(X);
Rx=X*X';%求协方差矩阵；
Rx1=Rx/sam;
Rxr=real(Rx1);
Rxi=imag(Rx1);
save Rxi(8).txt Rxi -ascii;
save Rxr(8).txt Rxr -ascii;
%dlmwrite('Rxi(8).txt',Rxi,'\t',0,0)
%dlmwrite('Rxr(8).txt',Rxr,'\t',0,0);
[V,D]=eig(Rx);
D=diag(D);
[D,pin]=sort(D,'descend');
Vn= V(:,pin([mm(1)+1:nn(1)],1));%噪声矩阵；
theta0=[0:80]';
tmp=[0:nn(1)-1]';
a2=-j*tmp*pi*sin(theta0'*degrad);
A2=exp(a2);
Ena=Vn'*A2;
den=diag(Ena'*Ena);
p=1./den;
%画图；
figure(1)
semilogy(theta0,abs(p));
xlabel('角度(度)');
ylabel('P(dB)');
title('MUSIC算法空间谱图');
% axis([10 30 0.1 1e6]);
grid on;
hold on;
%end