%%多点运动目标的RD成像试验算法
%% with Range-Doppler Algorithm（采用距离像-多普勒算法）
%%=======================================================================
clear;clc;close all;
%% Parameter--Constant
C=3e8; %propagation speed
%% Parameter--Radar characteristics 参数——雷达特征参数
Fc=1e9; %carrier frequency 1GHz 载频为1GHz
lambda=C/Fc; %wavelength 波长
%% Parameter--Target area 参数——目标域
Xmin=0; %target area in azimuth is within[Xmin,Xmax] 目标域方位向为[Xmin,Xmax]
Xmax=50;
Yc=10000; %center of imaged area 成像域的中心
Y0=500; %target area in range is within[Yc-Y0,Yc+Y0] 目标域距离向为[Yc-Y0,Yc+Y0]
%imaged width 2*Y0 成像宽度 2*Y0
%% Parameter--Orbital information 轨道信息
V=100; %SAR velosity 100 m/s SAR平台速度（水平速度）
H=5000; %height 5000 m 高度 5000m
R0=sqrt(Yc^2+H^2); %目标离SAR较远时，近似后的目标与SAR的垂直斜距（SAR平台与测绘带的垂直斜距）
%% Parameter--Antenna 参数——天线
D=4; %antenna length in azimuth direction 天线等效孔径
Lsar=lambda*R0/D; %SAR integration length 合成孔径长度
Tsar=Lsar/V; %SAR integration time 合成孔径时间
%% Parameter--Slow-time domain 参数——慢时域（即方位域）
Ka=-2*V^2/lambda/R0; %doppler frequency modulation rate 多普勒调频斜率
fd=abs(Ka*Tsar); %doppler frequency modulation bandwidth 多普勒调频带宽
PRF=fd; %pulse repitition frequency PRF相当于方位向上的采样率,
% 而根据复采样定律，采样率不能小于带宽；
PRT=1/PRF;
dt_s=PRT; %pulse repitition time 脉冲重复周期，也即为方位向上的采样间隔, 为方便运算令PRT等于dt_s；
Nslow=ceil((Xmax-Xmin+Lsar)/V/dt_s); %sample number in slow-time domain 慢时域中的采样数
%本步骤算出的Nslow仅有445
Nslow=2^nextpow2(Nslow); %for fft 为了便于做fft 将所得Nslow进行处理 取离该值最近的2的幂
sn=linspace((Xmin-Lsar/2)/V,(Xmax+Lsar/2)/V,Nslow);%discrete time array in slow-time domain
%慢时域中的离散时间序列
%% Parameter--Fast-time domain 参数——快时域（即距离域）
Tr=1e-5; %pulse duration 10us 脉冲持续时间为10us
Br=30e6; %chirp frequency modulation bandwidth 30MHz 线性调频脉冲信号（即chirp信号）的调频带宽为30MHz
Kr=Br/Tr; %chirp slope chirp信号的调频斜率
Fsr=2*Br; %sampling frequency in fast-time domain 快时域中的采样频率
dt_f=1/Fsr; %sample spacing in fast-time domain 快时域中的采样间隔
Rmin=sqrt((Yc-Y0)^2+H^2); %近距点距离
Rmax=sqrt((Yc+Y0)^2+H^2+(Lsar/2)^2); %远距点距离(之所以要考虑进Lsar是因为四面锥顶点到雷达距离才是Rmax)
Nfast=ceil(2*(Rmax-Rmin)/C/dt_f+Tr/dt_f); %sample number in fast-time domain 快时域中的采样数
Nfast=2^nextpow2(Nfast); %for fft 为了便于做fft 将所得Nslow进行处理 取离该值最近的2的幂
tm=linspace(2*Rmin/C,2*Rmax/C+Tr,Nfast); %discrete time array in fast-time domain 快时域中的离散时间序列
%（由于整个波束照射的目标距离不同，即从近距点回波（上升沿）到达时刻开始计算，直到远距点回波（下降沿）到达结束）
%% Parameter--Resolution 参数——分辨率
DY=C/2/Br; %range resolution 距离向分辨率
DX=D/2; %cross-range resolution 方位向分辨率
%由此可以看出X轴向即为航迹向（正测视时为方位向），Y轴为距离向；
%% Parameter--point targets 参数——点目标(目标形状暂时由点轮廓构成)
%format [x, y, reflectivity] 格式 [x坐标,y坐标,目标散射系数]
Vtar_X=5; %目标X与Y方向运动参数设置；
Vtar_Y=0;
delta_X=Vtar_X*PRT; %每个雷达脉冲间隔，运动目标所发生的距离变化；
delta_Y=Vtar_Y*PRT;
x=0; %设初始状态目标静止；
Xtar=[Xmin-50*DX+x*delta_X,Xmin+x*delta_X,Xmin+50*DX+x*delta_X,Xmin-50*DX+x*delta_X,Xmin+x*delta_X,Xmin+50*DX+x*delta_X];
%目标沿X轴方向运动；
Ytar=[Yc-x*delta_Y,Yc-x*delta_Y,Yc-x*delta_Y,Yc+50*DY-x*delta_Y,Yc+50*DY-x*delta_Y,Yc+50*DY-x*delta_Y]; %Y轴方向上点坐标不变；
rcs=ones(1,length(Ytar));
Ntarget=length(Ytar); %number of targets 目标数；
Pos=[Xtar',Ytar',rcs']; %position of targets 目标位置；
%% Generate the raw signal data 生成原始信号数据
K=Ntarget; %number of targets 目标数
N=Nslow; %number of vector in slow-time domain 慢时域中的列向量数
M=Nfast; %number of vector in fast-time domain 快时域中的列向量数
Srnm=zeros(N,M); %生成一个N*M的零矩阵（方位向发射N个脉冲，距离向采样得到M个样值点，则SAR回波为一N*M矩阵）
for k=1:1:K %k从1以1为步长变化到K
for x=0:1:N
sigma=Pos(k,3); %sigma为点目标的雷达散射截面,取矩阵T=Ptarget的第k行第3列元素，即目标位置中的目标散射系数
Dslow=sn*V-Pos(k,1); %一系列方位向上坐标，即X轴上对应各离散时间点的距离
R=sqrt(Dslow.^2+Pos(k,2)^2+H^2); %R也为一个行向量，表示各对应时刻的目标与雷达的斜距
tau=2*R/C; %tau为电磁波在雷达与目标之间传播的双程时间，也为一个行向量，各个不同位置不同时刻所对应的时间
Dfast=ones(N,1)*tm-tau'*ones(1,M); %一个关于时间差的N×M的矩阵，反应在Word文档中公式（3.4）中为[t（m）-tau]
%具体的物理含义即为时间差，当前采样时刻减去来回传播延迟；
phase=pi*Kr*Dfast.^2-(4*pi/lambda)*(R'*ones(1,M)); %表示整个单点回波信号的相位信息
Srnm=Srnm+sigma*exp(1i*phase).*(0<Dfast&Dfast<Tr).*((abs(Dslow)<Lsar/2)'*ones(1,M)); %包含约束条件的K个理想点目标的回波经采样后的表达式
Xtar=[Xmin-50*DX+x*delta_X,Xmin+x*delta_X,Xmin+50*DX+x*delta_X,Xmin-50*DX+x*delta_X,Xmin+x*delta_X,Xmin+50*DX+x*delta_X];
Ytar=[Yc-x*delta_Y,Yc-x*delta_Y,Yc-x*delta_Y,Yc+50*DY-x*delta_Y,Yc+50*DY-x*delta_Y,Yc+50*DY-x*delta_Y];
Pos=[Xtar',Ytar',rcs'];
end
end
%Srnm=Srnm+sigma*exp(1i*phase); 也可得到同样的成像结果；
%整个“.*(0<Dfast&Dfast<Tr).*((abs(Dslow)<Lsar/2)'*ones(1,M))”均为约束条件，不属于算法的必要部分；
%% Range compression 距离压缩
tr=tm-2*Rmin/C; %从0开始计算的相对离散时间阵列
Refr=exp(1i*pi*Kr*tr.^2); %表示参考信号
Sr=ifty(fty(Srnm).*(ones(N,1)*conj(fty(Refr)))); %运用FFT来进行实质为匹配滤波的距离压缩处理
Gr=abs(Sr); %对处理结果取绝对值
%% Azimuth compression 方位压缩
ta=sn-Xmin/V; %相对离散时间阵列，但由于Xmin=0，故ta=sn
Refa=exp(1i*pi*Ka*ta.^2); %(abs(ta)<Tsar/2); %参考信号（仿真中，调频斜率 已知，因此不需要进行Doppler参数估计）
Sa=iftx(ftx(Sr).*(conj(ftx(Refa)).'*ones(1,M))); %运用FFT来进行实质为匹配滤波的方位压缩处理
Ga=abs(Sa); %对处理结果取绝对值
figure(1)
imagesc(Gr); %作距离压缩后的Sr灰度值图像
xlabel('Range'),ylabel('Azimuth'),
figure(2)
imagesc(Ga); %作方位压缩后的Ga灰度值图像
xlabel('Range'),ylabel('Azimuth'),
figure(3)
mesh(Ga); %作Ga的三维网格图
