%程序清理段
close all;
clc;
%%
%参数定义段
f0 = 10e9; %载频(本振)
tau = 10e-6; %脉冲宽度
B = 10e6; %带宽
C = 3e8; %光速
fs = 100e6; %采样频率
Rd = 300; %两个雷达之间的距离
R1 = 1000; %节点一雷达到目标的距离
R2 = 1500; %节点二雷达到目标的距离
Gate_echo = 3500; %距离波门
det_f1 = randn; %节点一本振频率误差
det_f2 = randn; %节点二本振频率误差
det_fai1 = randn; %节点一本振相位误差
det_fai2 = randn; %节点二本振相位误差
fai1 = randn; %节点一中频相位误差
fai2 = randn; %节点二中频相位误差
det_t = 1e-7*randn; %节点一时钟和节点二时钟系统误差
K = B/tau; %线性调频系数
Td = Rd/C ; %直达波延时
T1 = R1/C; %节点一到达目标的延时
T2 = R2/C; %节点二到达目标的延时
%%
%信号构建
Gate_signal = C*tau/2; %信号脉宽对应距离
Gate = Gate_echo+Gate_signal; %距离波门加脉宽
Gate_Number = round(Gate/C*fs); %波门内采样点个数
t = 0:1/fs:tau*2; %信号长度
t_1 = t; %节点二的时间
echo_t = linspace(0,Gate/C,Gate_Number); %波门长度
trasmit_signal1 = rectpuls((t_1-tau/2)/(tau)).*exp(1i*2*pi*K*(t_1-tau/2).^2).*exp(1i*2*pi*(det_f1)*t+1i*fai1); %节点一发射信号
echo_direct1_2 = rectpuls((t_1-det_t-Td-tau/2)/(tau)).*exp(1i*2*pi*K*(t_1-det_t-Td-tau/2).^2).*exp(1i*2*pi*(f0+det_f1).*(t_1-det_t-Td)...
+1i*fai1-1i*2*pi*(f0+det_f2).*t_1-1i*det_fai2);
echo_direct1_2 =echo_direct1_2 + sqrt(1e-6)*(randn(size(echo_direct1_2))+1i*randn);%echo_direct1_2; %节点二接收直达波信号下变频
echo_1_2 = rectpuls((t_1-det_t-T1-T2-tau/2)/(tau)).*exp(1i*2*pi*K*(t_1-det_t-T1-T2-tau/2).^2).*exp(1i*2*pi*(f0+det_f1).*(t_1-det_t-T1-T2)...
+1i*fai1-1i*2*pi*(f0+det_f2).*t_1-1i*det_fai2);%+sqrt(1e-4)*(randn(size(echo_1_2))+1i*randn); %节点二接收目标信号下变频
echo_1_1 =rectpuls((t_1-2*T1-tau/2)/(tau)).*exp(1i*2*pi*K*(t_1-2*T1-tau/2).^2).*exp(1i*2*pi*(f0+det_f1).*(t_1-2*T1)...
+1i*fai1-1i*2*pi*(f0+det_f2).*t_1-1i*det_fai2);%+sqrt(1e-4)*(randn(size(echo_direct1_2))+1i*randn); %节点一接收信号目标下变频
trasmit_signal2 = rectpuls((t-tau/2)/(tau)).*exp(1i*2*pi*K*(t-tau/2).^2).*exp(1i*2*pi*(det_f2)*t+1i*fai2); %节点二发射信号
echo_direct2_1 = rectpuls((t+det_t-Td-tau/2)/(tau)).*exp(1i*2*pi*K*(t+det_t-Td-tau/2).^2).*exp(1i*2*pi*(f0+det_f2).*(t+det_t-Td)+1i*fai2...
-1i*2*pi*(f0+det_f1).*t-1i*det_fai1);
echo_direct2_1 = echo_direct2_1+sqrt(1e-6)*(randn(size(echo_direct2_1))+1i*randn); %节点一接收直达波信号下变频
echo_2_1 = rectpuls((t+det_t-T1-T2-tau/2)/(tau)).*exp(1i*2*pi*K*(t+det_t-T1-T2-tau/2).^2).*exp(1i*2*pi*(f0+det_f2).*(t+det_t-T1-T2)+1i*fai2...
-1i*2*pi*(f0+det_f1).*t-1i*det_fai1); %节点一接收目标信号下变频
echo_2_2 = rectpuls((t-2*T2-tau/2)/(tau)).*exp(1i*2*pi*K*(t-2*T2-tau/2).^2).*exp(1i*2*pi*(f0+det_f2).*(t-2*T2)+1i*fai2...
-1i*2*pi*(f0+det_f1).*t-1i*det_fai1); %节点二接收目标信号下变频
%%
%脉压
%直达波脉压
FFT_Number = 2^nextpow2(length(t)+Gate_Number-1); %FFT点数
Srw1 = fft(trasmit_signal1,FFT_Number);
Sw1 = fft(echo_direct1_2,FFT_Number);
Sot1_1 = ifft(Sw1.*conj(Srw1));
Z1 = abs(Sot1_1(1:Gate_Number)); %这里只保留了接受窗内的数
Z = Z1/max(Z1);
Z = 20*log10(Z); %db形式
[maxvalue1,maxposition1] = max(abs(Sot1_1)); %寻找最大值点
%对直达波的信号进行处理
t1=t_1(maxposition1);
times=16;
t_1_up = (0:512*times-1)/fs/times;
up_ind_1 = mod((1:513),length(Sot1_1));
Sot1 = interpft(Sot1_1(up_ind_1),512*times);
[maxvalue_1,maxposition_1] = max(abs(Sot1));
det1=t_1_up(maxposition_1) - Td;
N_1=round((Td+det_t)*fs*16);
max_ideal_1=abs(Sot1(N_1));
% Sot1 = interpft(Sot1_1(1:512),512*16);
% [maxvalue_1,maxposition_1] = max(Sot1);
% r_n_1 = (1:1/16:513)/fs*C;
% r_n_1 = r_n_1(1:end-1);
% rmaxp_1 = r_n_1(maxposition_1)
% t_maxp1 = maxposition_1/C;
% max_fai1 =2*pi*(f0+det_f1).*(t_maxp1-det_t-Td)+fai1-2*pi*(f0+det_f2).*t_maxp1-det_fai2;
Srw2 = fft(trasmit_signal2,FFT_Number);
Sw2 = fft(echo_direct2_1,FFT_Number);
Sot2_1 = ifft(Sw2.*conj(Srw2));
Y1 = abs(Sot2_1(1:Gate_Number)); %这里只保留了接受窗内的数
Y = Y1/max(Y1);
Y = 20*log10(Y); %db形式
[maxvalue2,maxposition2] = max(abs(Sot2_1)); %寻找最大值点
%对直达波的信号进行处理
t2=t(maxposition2);
times=16;
t_2_up = (0:512*times-1)/fs/times;
up_ind_2 = mod((1:513),length(Sot2_1));
Sot2 = interpft(Sot2_1(up_ind_2),512*times);
[maxvalue_2,maxposition_2] = max(abs(Sot2));
det2=t_2_up(maxposition_2) - Td;
N_2=round((Td-det_t)*fs*16);
max_ideal_2=abs(Sot2(N_2));
% Sot2 = interpft(Sot2_1(1:512),512*16);
% [maxvalue_2,maxposition_2] = max(Sot2);
% r_n_2 = (1:1/16:513)/fs*C;
% r_n_2 = r_n_2(1:end-1);
% rmaxp_2 = r_n_2(maxposition_2)
% t_maxp2 = maxposition_2/C;
% max_fai2 =2*pi*(f0+det_f2).*(t_maxp2+det_t-Td)+fai2-2*pi*(f0+det_f1).*t_maxp2-det_fai1;
% det_max_fai1=max_fai2-max_fai1
%计算出带误差的节点一与节点二的系统延时
det_t_j=(1/2)*(N_1-N_2)/16/fs;
%目标信号脉压
Srw3 = fft(trasmit_signal1,FFT_Number);
Sw3 = fft(echo_1_2,FFT_Number);
Sot3_1 = ifft(Sw3.*conj(Srw3));
W1 = abs(Sot3_1(1:Gate_Number)); %这里只保留了接受窗内的数
W = W1/max(W1);
W = 20*log10(W); %db形式
[maxvalue3,maxposition3] = max(abs(Sot3_1)); %寻找最大值点
%对目标的信号进行处理
t3=t_1(maxposition3);
times=16;
t_3_up = t_1(maxposition3)-255/fs+(0:512*times-1)/fs/times;
up_ind_3 = mod(maxposition3 + (-255:256),length(Sot3_1));
Sot3 = interpft(Sot3_1(up_ind_3),512*times);
[maxvalue_3,maxposition_3] = max(abs(Sot3));
det3=t_3_up(maxposition_3) - T1 -T2;
N_3=round(1+(T1+T2+det_t_j-t_1(maxposition3)+255/fs)*fs*16);
max_ideal_3=abs(Sot3(N_3));
% Sot3 = interpft(Sot3_1((maxposition3-255):(maxposition3+256)),512*16);
% [maxvalue_3,maxposition_3] = max(Sot3);
% r_n_3 = (maxposition_3-255:1/16:maxposition_3+256)/fs*C;
% r_n_3 = r_n_3(1:end-1);
% rmaxp_3 = r_n_3(maxposition_3)
% t_maxp3 = maxposition_3/C;
% max_fai3 = 2*pi*(f0+det_f1).*(t_maxp3-det_t-T1-T2)+fai1-2*pi*(f0+det_f2).*t_maxp3-det_fai2;
Srw4 = fft(trasmit_signal2,FFT_Number);
Sw4 = fft(echo_2_1,FFT_Number);
Sot4_1 = ifft(Sw4.*conj(Srw4));
Q2 = abs(Sot4_1(1:Gate_Number)); %这里只保留了接受窗内的数
Q = Q2/max(Q2);
Q = 20*log10(Q); %db形式
[maxvalue4,maxposition4]