% 收发多径都考虑的情况
clear all
close all
clc
c=3e8;
f0 = 150e6; % 载波频率
fd = 100; % 多普雷频率
lemda = c/f0;
d =lemda/1; % 阵元间距
snap=10; % 快拍数
P=2; % 信源个数（直射信号与镜像反射信号）
Ne = 11; % 发射阵元个数
Nr = 11; % 接收阵元个数
Te = 100e-6; % 脉宽
Tr = 1000e-6; % 1e-3*3e8 = 3e5 最大测距范围 1000us
T0 = Tr; % 100ms
center_f = (Ne-1)/2;
delta_f = 1/Te;
B = Ne*delta_f;
Fs=2*B; % 采样速率
Ts=1/Fs;
R = 100e3; % 阵列目标的距离
ht=4000; % 目标高度
ha=100; % 天线架高
Rd=sqrt(R^2+(ht-ha)^2); % 直达波程
Ri=sqrt(R^2+(ht+ha)^2); % 反射波程
range = [Rd Ri];
theta1 = atan((ht-ha)/R);
theta2 = -atan((ht+ha)/R);
theta = [theta1;theta2];
%theta = [-1;8]/180*pi;
% theta = [theta1;theta2]/pi*180; % DOA
amplify = [1 .9];
%amplify = [1 1];
% amplify = [1 0.9*exp(j*pi)];
% t=0:Ts:T0-Ts; %
t=0:Ts:T0-Ts; % 时间坐标
tao = 0+Te/2:Tr:T0; % Tr
rect_me1 = pulstran(t,tao,'rectpuls',Te);
% s_num = double(int32(Tr/Ts));
S = zeros(Nr,length(t),Ne);
X = zeros(Nr,length(t),Ne);
%%%%%%%%%%%% 发射、接收直达波和反射波都有,相当于四条路径,产生信号
for loop1 = 0:Nr-1 % 接收
for loop2 = 0:Ne-1 % 发射
for loop = 1:2 % 发射直达和发射反射
fj = f0+(loop2-center_f)*delta_f;
at = amplify(loop)*exp(j*2*pi*loop2*d*sin(theta(loop))./lemda);
tao = 2*range(loop)/c;
% tao = 0;
temp = at*rect_me1.*exp(j*2*pi*fj*(t-tao)); % 暂存发射信号,直接用fj做ML不对，因为是异频
% temp = at*exp(j*2*pi*f0*(t-tao)); % 暂存发射信号，如果用f0可以，
for loop3 = 1:2 % 接收直达和反射
% ar = exp(j*2*pi*loop1*d*sin(theta(loop))./lemda);
ar = amplify(loop3)*exp(j*2*pi*(loop1)*d*sin(theta(loop3))./lemda); % 接收直达和反射
S(loop1+1,:,loop2+1) = S(loop1+1,:,loop2+1) + ar*temp; % 乘上接收导向矢量
end
end
end
end
%%%%%%%%%%%% processing
p=1;
% for SNR = -15:5;
SNR = 20;
snr = 10^(SNR/20); %%信噪比
mentcronumber = 10;
for mentcro = 1:mentcronumber
noise = normrnd(0,1,size(S));
X = snr*S + noise;
% X = snr*S;
% Ne /---------------/|
% / / |
% Nr /---------------/ /
% | | /
% 1 |_______________|/
temp_X = zeros(1,length(t));
temp_XXX = zeros(Nr,Ne);
% temp_ss = zeros(11,1);
temp0 = zeros(Nr,length(t));
Rx_ss = zeros(Nr,Nr);
k1=1;
for loop = 1:Nr
for loop3 = 0:Ne-1 % 发射阵元
tao = 2*range(1)/c;
fj = f0+(loop3-center_f)*delta_f;
complement = exp(j*2*pi*fj*(t-tao));%
temp_X = X(loop,:,loop3+1);
temp_XXX((loop3)*Ne + k1) = temp_X*complement';
% temp_XXX(loop3+1,k1) = temp_X*complement';
end
k1=k1+1;
% %%%
% for loop0 = 1:Nr
% temp0(loop0,:) = X(loop,:,loop0);
% end
% temp_ss = temp0;
% % temp_ss = X(:,:,loop);
% Rx_ss = Rx_ss + temp_ss*temp_ss';
end
%%%%%%%%%%%%%%%%%%%% SS
for loop = 1:Nr
temp_ss(:,1) = temp_XXX((loop-1)*Nr+1:loop*Nr);
Rx_ss = Rx_ss + temp_ss*temp_ss';
end
[V,D]=eigs(Rx_ss/Ne,Ne); % 特征分解
G=V(:,P+1:Ne)*V(:,P+1:Ne)';
hang = [0:Nr-1];
k=1;
for loop = -10:0.01:10
a = exp(j*2*pi*d*hang.'*sin(loop/180*pi)./lemda);
power(k) = 1/abs(a'*G*a);
k=k+1;
end
figure(1)
plot(-10:0.01:10,20*log10(abs(power)),'b-')
grid on
hold on
end
xlabel('DOA(度)')
ylabel('MUSIC谱峰')
title('目标高度4000m')