p3.m

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

%姓名：吴茜  学号：xxxxxxxxx  08.6.26
%(3)在信噪比变化范围为-20dB－+20dB时，分别采用MUSIC算法,CAPON算法, ESPRIT算法，
%   对上述四个信号源的波达方向进行估计,并画出波达方向估计误差,随着信噪比变化的曲线图

clear all; 
clc;
close all;
format short

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

%=========求阵列接受数据的协方差矩阵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
estmusic_theta(jj,:)=peak_seek(Pmusic,theta,N);
theta_error_music(jj)=sum(abs(estmusic_theta(jj,:)-sort(angle)))/N;   %估计偏差

%=========capon算法===============
L=length(theta);
Pcapon=capon_sp(x,L,d/lamda);
estcapon_theta(jj,:)=peak_seek(Pcapon,theta,N);
theta_error_capon(jj)=sum(abs(estcapon_theta(jj,:)-sort(angle)))/N;   %估计偏差

%=========esprit方法=============
est_esprit_theta=transpose(esprit_doa(x,N,d/lamda));
theta_error_esprit(jj)=sum(abs(est_esprit_theta-[-45 60 45 30]))/N;      %估计偏差

end

figure(1)
plot(snr,theta_error_music,'r-*');grid on;hold on;
plot(snr,theta_error_capon,'g-o');grid on;hold on;
plot(snr,theta_error_esprit,'-s');grid on;hold on;
title('波达方向估计误差,随着信噪比变化的曲线图');
xlabel('信噪比');
ylabel('估计误差');
legend('music','capon','esprit');