%姓名：吳茜  學號：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');