clc;
clear;
close all;
tic
M=40;%陣元數(shù)目
N=3;%信源數(shù)目
snap=1000;%快拍數(shù)目
C=3e8;
f0=10e6;
lamda=C/f0;
d=0.5*lamda;
% k=d/lamda;
theta0=22;
theta1=20;
theta2=-5;
fs=1000;
ts=1/fs;
t=(0:snap-1)*ts;
a=[0:M-1]';%陣列矢量
u0=5;
u1=10;
u2=20;


s0=exp(j*2*pi*(f0*t+0.5*u0*t.^2));
s1=exp(j*2*pi*(f0*t+0.5*u1*t.^2));
s2=exp(j*2*pi*(f0*t+0.5*u2*t.^2));

%陣列流行矢量

a_theta0=exp(j*2*pi*d/lamda*a*sin(theta0/180*pi));
a_theta1=exp(j*2*pi*d/lamda*a*sin(theta1/180*pi));
a_theta2=exp(j*2*pi*d/lamda*a*sin(theta2/180*pi));

A=[a_theta0 a_theta1 a_theta2];
S=[s0;s1;s2];
X0=A*S;
SNR=10;
randn('state',0);
real_noise=randn(size(X0));
randn('state',3);
imag_noise=randn(size(X0));
noise0=(real_noise+j*imag_noise)/2^0.5;
noise=10^(-SNR/20)*noise0;
X=X0+noise;%%%%%%%%%%------------=========完整的基帶信號(hào)


R_x=X*X'/length(t);
[V,D]=eig(R_x);
Us=V(:,M-N+1:M);
Un=V(:,1:M-N);




theta=-90:0.1:90;
for ii=1:length(theta)
    a_theta=exp(j*2*pi*d/lamda*a*sin(theta(ii)*pi/180));
    P_A=a_theta*inv(a_theta'*a_theta)*a_theta';
    Pssf1(ii)=abs(trace(P_A*Us*Us'));%信號(hào)子空間擬合
    Pssf2(ii)=1/abs(trace(Un'*a_theta*a_theta'*Un));%噪聲子空間擬合
end
figure(1);
plot(theta,Pssf1);title('信號(hào)子空間擬合');
grid on;
figure(2);
plot(theta,Pssf2);title('噪聲子空間擬合');
grid on;
figure(3);
plot(theta,[Pssf1/max(Pssf1);Pssf2/max(Pssf2)]);title('空間擬合技術(shù)');
legend('信號(hào)子空間擬合','噪聲子空間擬合');
grid on;