%%========================================================
clear;clc;close all;
%%========================================================
%%Parameter--constant
C=3e8; %propagation speed
%%Parameter--radar characteristics
Fc=1e9; %carrier frequency 1GHz载波频率
lambda=C/Fc; %wavelength
%%Parameter--target area
Xmin=0; %target area in azimuth is within[Xmin,Xmax]
Xmax=50;
Yc=10000; %center of imaged area 成像区域中线
Y0=500; %target area in range is within[Yc-Y0,Yc+Y0]
%imaged width 2*Y0
%%Parameter--orbital information
V=100; %SAR velosity 100 m/s
H=5000; %height 5000 m
R0=sqrt(Yc^2+H^2);
%%Parameter--antenna
D=4; %antenna length in azimuth direction方位向天线长度
Lsar=lambda*R0/D; %SAR integration lengthSAR合成孔径长度
Tsar=Lsar/V; %SAR integration time合成孔径时间 SAR照射时间
%%Parameter--slow-time domain方位向
Ka=-2*V^2/lambda/R0; %doppler frequency modulation rate多普勒调制频率
Ba=abs(Ka*Tsar); %doppler frequency modulation bandwidth多普勒调制带宽
PRF=Ba; %pulse repitition frequency脉冲重复频率 PRF其实为多普勒频率的采样率，又为复频率，所以等于Ba
PRT=1/PRF; %pulse repitition time脉冲重复时间
ds=PRT; %sample spacing in slow-time domain 慢时域的时间步长，所谓步长就是你将需要观测的数值均匀分成若干个区间，每个区间的长度就叫步长。时间步长法应该是以时间轴为主自变量，确定步长后定点测值，简单来说就是每隔一段时间取一个值记录下来。
Nslow=ceil((Xmax-Xmin+Lsar)/V/ds); %sample number in slow-time domain
Nslow=2^nextpow2(Nslow); %for fft
sn=linspace((Xmin-Lsar/2)/V,(Xmax+Lsar/2)/V,Nslow);%discrete time array in slow-time domain 慢时间域的时间矩阵%（linspace用法：linspace(x1,x2,N) 用于产生x1，x2之间的N点行矢量，相邻数据跨度相同。其中x1、x2、N分别为起始值、终止值、元素个数。若缺省N，默认点数为100。）
PRT=(Xmax-Xmin+Lsar)/V/Nslow; %refresh 由于Nslow改变了，所以相应的一些参数也需要更新，周期减小了
PRF=1/PRT;
ds=PRT;
%%Parameter--fast-time domain距离向
Tr=5e-6; %pulse duration 10us
Br=30e6; %chirp frequency modulation bandwidth 30MHz
Kr=Br/Tr; %chirp slope
Fsr=3*Br; %sampling frequency in fast-time domain
dt=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);
Nfast=ceil(2*(Rmax-Rmin)/C/dt+Tr/dt);%sample number in fast-time domain
Nfast=2^nextpow2(Nfast); %for fft
tm=linspace(2*Rmin/C,2*Rmax/C+Tr,Nfast); %discrete time array in fast-time domain 快时域的离散时间矩阵
dt=(2*Rmax/C+Tr-2*Rmin/C)/Nfast; %refresh
Fsr=1/dt;
%%Parameter—resolution分辨率
DY=C/2/Br; %range resolution距离向分辨率
DX=D/2; %cross-range resolution方位向分辨率
%%Parameter--point targets点目标
Ntarget=3; %number of targets
%format [x, y, reflectivity] %[x坐标，y坐标，分辨率] 点目标格式[x,y,反射系数sigma]
Ptarget=[Xmin,Yc,1 %position of targets
Xmin,Yc+10*DY,1
Xmin+20*DX,Yc+50*DY,1];
disp('Parameters:')
disp('Sampling Rate in fast-time domain');disp(Fsr/Br)
disp('Sampling Number in fast-time domain');disp(Nfast)
disp('Sampling Rate in slow-time domain');disp(PRF/Ba)
disp('Sampling Number in slow-time domain');disp(Nslow)
disp('Range Resolution');disp(DY)
disp('Cross-range Resolution');disp(DX)
disp('SAR integration length');disp(Lsar)
disp('Position of targets');disp(Ptarget)
%%========================================================
%%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距离向采样数?
T=Ptarget; %position of targets目标位置
Srnm=zeros(N,M); %生成零矩阵存储回波信号
for k=1:1:K %总共K个目标
sigma=T(k,3); %得到目标散射系数
Dslow=sn*V-T(k,1); %方位向距离，投影到方位向的距离
R=sqrt(Dslow.^2+T(k,2)^2+H^2); %实际距离矩阵
tau=2*R/C; %回波相对于发射波的延时
Dfast=ones(N,1)*tm-tau'*ones(1,M); %%(t-tau)，其实就是时间矩阵，ones(N,1)和ones(1,M)都是为了将其扩展为矩阵
phase=pi*Kr*Dfast.^2-(4*pi/lambda)*(R'*ones(1,M)); %相位
Srnm=Srnm+sigma*exp(j*phase).*(0<Dfast&Dfast<Tr).*((abs(Dslow)<Lsar/2)'*ones(1,M)); %由于是多个目标反射的回波，所以此处进行叠加
end
%%========================================================
%%Range compression距离向压缩 矩阵行做FFT变换
tr=tm-2*Rmin/C;
Refr=exp(j*pi*Kr*tr.^2).*(0<tr&tr<Tr);
Sr=ifty(fty(Srnm).*(ones(N,1)*conj(fty(Refr))));
Gr=abs(Sr);
%%Azimuth compression方位向压缩
ta=sn-Xmin/V;
Refa=exp(j*pi*Ka*ta.^2).*(abs(ta)<Tsar/2);
Sa=iftx(ftx(Sr).*(conj(ftx(Refa)).'*ones(1,M)));
Ga=abs(Sa);
%%========================================================
%%graw the intensity image of signal
colormap(gray); %colormap设置图像的颜色；
figure(1)
subplot(211);
row=tm*C/2-2008;col=sn*V-26; %row距离向 col方位向
imagesc(row,col,255-Gr); %intensity image of Sr imagesc(A) 将矩阵A中的元素数值按大小转化为不同颜色，并在坐标轴对应位置处以这种颜色染色imagesc(x,y,A) x,y决定坐标范围，x,y应是两个二维向量，即x=[x1 x2],y=[y1 y2],matlab会在[x1,x2]*[y1,,y2]的范围内染色。
axis([Yc-Y0,Yc+Y0,Xmin-Lsar/2,Xmax+Lsar/2]);
xlabel('\rightarrow\itRange in meters'),ylabel('\itAzimuth in meters\leftarrow'),
title('Stripmap SAR after range compression'),
subplot(212);
imagesc(row,col,255-Ga); %intensity image of Sa
axis([Yc-Y0,Yc+Y0,Xmin-Lsar/2,Xmax+Lsar/2]);
xlabel('\rightarrow\itRange in meters'),ylabel('\itAzimuth in meters\leftarrow'),
title('Stripmap SAR after range and azimuth compression'),
%%========================================================
%%draw 3D picture 两点目标的回波仿真3D图
figure(2)
waterfall(real(Srnm((200:205),:)));axis tight
xlabel('Range'),ylabel('Azimuth'),
title('Real part of the raw signal'),
figure(3)
waterfall(Gr((200:205),(600:1000)));axis tight
xlabel('Range'),ylabel('Azimuth'),
title('Stripmap SAR after range compression'),
figure(4)
mesh(Ga((200:300),(750:860)));axis tight %mesh —— 三维网线绘图函数 mesh(z) —— z为n×m的矩阵，x与y坐标为元素的下标
xlabel('Range'),ylabel('Azimuth'),
title('Stripmap SAR after range and azimuth compression'),
%%========================================================
%%draw -3dB contour 两点目标压缩后的3dB等高线图
figure(5)
a=max(max(Ga));
contour(row,col,Ga,[0.707*a,a],'r');grid on
axis([9995,10050,-20,20]),
xlabel('\rightarrow\itRange in meters'),ylabel('\itAzimuth in meters\leftarrow'),
title('Resolution Demo: -3dB contour'); %分辨率演示
%%========================================================
