p2.m

已知：信号中心波长为2,天线阵元的间距为1米,快拍数为2000

%姓名: 吴茜   学号:xxxxxxxxxxxxx    08.6.26
%(1)在不加入噪声的情况下,观察并计算协方差矩阵特征值,并对它的特点加以说明;
%(2)分别采用MUSIC算法,CAPON算法, ESPRIT算法在下面两种情况下，
%   对上述四个信号源的波达方向进行估计,并画出它们的空间谱图；
%     ①在不加入噪声的情况下(注释掉Line35和Line36即可)
%     ②在加入高斯白噪声的情况下,假设信噪比为10dB

clear all; 
clc;
close all;
format short

snap=2000;                  %快拍数
n_x=0:snap-1;
M=7;                       %阵元数     
N=4;                       %信源数
K=0:M-1;
fs=100;                    %采样率
c=300000000;
lamda=2;
d=lamda/2;  
f0=c/lamda ;               %目标信号的中心频率
t=[0:snap-1]/fs;
angle=[45,60,30,-45];        %信号方位角
fre=[0.015,0.05,0.02,0.035];  %信号归一化频率

%============产生阵元输出数据矩阵==============
for ii=1:N
  S(ii,:)=exp(j*2*pi*(f0*fre(ii)*t+0.5*2^(ii-1)*t.^2));%确保不相干
 %S(ii,:)=exp(j*2*pi*f0*fre(ii)*t);
 A(:,ii)=exp(j*pi*K'*sin(angle(ii)/180*pi));
end
  x=A*S;
  
% snr=10;           %定义信噪比
% x=awgn(x,snr);    %加高斯白噪声

%=========求阵列接受数据的协方差矩阵R,查看其特征值D===========
R=(x*x')/snap;
[V,D] = eig(R);
                                
%==========music算法============
for i=1:M
    DD(i)=D(i,i);
end
[h,e]=sort(DD);
for i=1:M-N
Vn(:,i)=V(:,e(i));
end
theta=-90:0.1:90;
for ii=1:length(theta)
   a_theta=exp(j*2*pi*d/lamda*K'*sin(theta(ii)/180*pi));
   Pmusic(ii)=10*log10(1/abs(a_theta'*Vn*Vn'*a_theta));
end
figure(1)
 plot(theta,Pmusic);
 grid on;
 xlabel('角度/degree');
 ylabel('功率谱');
 title('music算法');
%=========capon算法===============
L=length(theta);
Pcapon=capon_sp(x,L,d/lamda);
figure(2)
 plot(theta,Pcapon);
 grid on;
 xlabel('角度/degree');
 ylabel('功率谱');
 title('capon算法');
%=========esprit方法=============
est_theta=esprit_doa(x,N,d/lamda);