% 收發多徑都考慮的情況
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')
————————————————
聯系：highspeedlogic
QQ ：1224848052
微信：HuangL1121
郵箱：1224848052@qq.com
網站：http://www.mat7lab.com/
網站：http://www.hslogic.com/
微信掃一掃：
