%%========================================================
clear;clc;close all;format long;
%%========================================================
%%Parameter--constant
C=3e8; %propagation speed
%%Parameter--radar characteristics
Fc=10e9; %carrier frequency 10GHz
lambda=C/Fc; %wavelength
%%Parameter--target area
Xmin=-100; %target area in azimuth is within[Xmin,Xmax]
Xmax=100;
Yc=10000; %center of imaged area
Y0=1000; %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
%D=0.5;
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
Ba=abs(Ka*Tsar); %doppler frequency modulation bandwidth
PRF=Ba; %pulse repitition frequency
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
PRT=(Xmax-Xmin+Lsar)/V/Nslow; %refresh
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=-Tr/2+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=9; %number of targets
%Ptarget=[0,Yc,255];
%format [x, y, reflectivity]
Ptarget=[Xmin+20,Yc-200,255 %position of targets
Xmin+20,Yc,255
Xmin+20,Yc+400,255
0,Yc-200,255
0,Yc,255
0,Yc+400,255
40,Yc-200,255
40,Yc,255
40,Yc+400,255];
%%========================================================
%%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
temp_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);
phase=pi*Kr*Dfast.^2-(4*pi/lambda)*(R'*ones(1,M));
Srnm=Srnm+temp_sigma*exp(j*phase).*(-Tr/2<=Dfast&Dfast<=Tr/2).*((abs(Dslow)<=Lsar/2)'*ones(1,M));
end
clear Dfast;
clear Dslow;
clear phase;
%%========================================================
%%Range compression
tr=tm-2*Rmin/C;
Refr=exp(j*pi*Kr*tr.^2).*(-Tr/2<=tr&tr<=Tr/2);
temp1=fty(Srnm);
temp2=conj(ones(Nslow,1)*fty(Refr));
Sr=zp_ifty(temp1.*temp2); %Target positioning in range direction is determined by time correlation operation via FFT, IFFT result 0 in range direction versus Rmin
clear temp1;
clear temp2;
Gr=abs(Sr);
%%Azimuth compression
ta=sn-Xmin/V;
Refa=exp(j*pi*Ka*ta.^2).*(abs(ta)<=Tsar/2);
Sa=zp_iftx(ftx(Sr).*(conj(ftx(Refa)).'*ones(1,M))); %Target positioning in azimuth direction is determined by time correlation operation via FFT, IFFT result 0 in azimuth direction versus Xmin
Ga=abs(Sa);
clear Sr;
clear Sa;
%%========================================================
%%graw the intensity image of signal
row1=sqrt((tm*C/2).^2-H^2);
col1=sn*V;
row2=sqrt(((tm+Tr/2)*C/2).^2-H^2);
col2=(sn+Lsar/(2*V))*V;
scene=zeros(Nslow,Nfast);
for k=1:1:K
temp_x=T(k,1);
temp_y=T(k,2);
temp_sigma=T(k,3);
temp_col=abs(col2-temp_x);
temp_index_x=find(temp_col==min(temp_col));
temp_row=abs(row2-temp_y);
temp_index_y=find(temp_row==min(temp_row));
scene(temp_index_x,temp_index_y)=temp_sigma;
end
figure();
colormap(gray(256));
imagesc(row2,col2,scene);
clear scene;
axis([Yc-Y0,Yc+Y0,Xmin,Xmax]);
xlabel('\rightarrow\itRange in meters');
ylabel('\itAzimuth in meters\leftarrow');
title('Target distribution in scene');
figure();
colormap(gray(256));
imagesc(row1,col1,abs(Srnm));
clear Srnm;
axis([Yc-Y0,Yc+Y0,Xmin-Lsar/2,Xmax+Lsar/2]);
xlabel('\rightarrow\itRange in meters');
ylabel('\itAzimuth in meters\leftarrow');
title('target echo');
figure();
colormap(gray(256));
imagesc(row2,col1,Gr); %intensity image of Sr
clear Gr;
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');
figure();
colormap(gray(256));
imagesc(row2,col2,Ga); %intensity image of Sa
clear Ga;
axis([Yc-Y0,Yc+Y0,Xmin,Xmax]);
xlabel('\rightarrow\itRange in meters');
ylabel('\itAzimuth in meters\leftarrow');
title('Stripmap SAR after range and azimuth compression');
%%draw -3dB contour
% figure();
% a=max(max(Ga));
% contour(row2,col2,Ga,[0.5*a,a],'b');grid on
% axis([Yc-Y0,Yc+Y0,Xmin-Lsar/2,Xmax+Lsar/2]);
% xlabel('\rightarrow\itRange in meters');
% ylabel('\itAzimuth in meters\leftarrow');
% title('Resolution Demo: -3dB contour');
%%=======================================================
