%%產生Stripmap SAR的回波
clear all

thetaT=0;%T平臺波束斜視角
thetaT=thetaT*pi/180;%rad
thetaR=0;%R平臺波束斜視角
thetaR=thetaR*pi/180;
c=3e8;%光速
fc=1.5e9;%載頻
lambda=c/fc;%波長

%%測繪帶區域
X0=200;%方位向[-X0,X0]
Rtc=3000; 
Rrc=3000;
Rc=(Rtc+Rrc)/2;
R0=150;%距離向[Rc-R0,Rc+R0]

%%距離向(Range),r/t domain
Tr=1.5e-6;%LFM信號脈寬 1.5us (200m)
Br=150e6; %LFM信號帶寬 150MHz
Kr=Br/Tr; %調頻斜率
Nr=512;
r=Rc+linspace(-R0,R0,Nr);
t=2*r/c;%t域序列
dt=R0*4/c/Nr;%采樣周期
f=linspace(-1/2/dt,1/2/dt,Nr);%f域序列

%%方位向(Azimuth,Cross-Range),x/u domain
v=100;%SAR 平臺速度
Lsar=300;%合成孔徑長度
Na=1024;
x=linspace(-X0,X0,Na);%u域序列
u=x/v;
du=2*X0/v/Na;
fu=linspace(-1/2/du,1/2/du,Na);%fu域序列
ftdc=v*sin(thetaT);
ftdr=-(v*cos(thetaT))^2/lambda/Rtc;
frdc=v*sin(thetaR);
frdr=-(v*cos(thetaR))^2/lambda/Rrc;
fdc=ftdc+frdc;%Doppler調頻中心頻率
fdr=ftdr+frdr;%Doppler調頻斜率
%％目標位置
Ntar=3;%目標個數
Ptar=[Rrc,0,1 %距離向坐標,方位向坐標,sigma             
          Rrc+50,-50,1
         Rrc+50,50,1];


%%產生回波
s_ut=zeros(Nr,Na);
U=ones(Nr,1)*u;%擴充為矩陣
T=t'*ones(1,Na);
for i=1:1:Ntar
    rn=Ptar(i,1);xn=Ptar(i,2);sigma=Ptar(i,3);
    rtn=rn+Rtc-Rrc;
    RT=sqrt(rtn^2+(rtn*tan(thetaT)+xn-v*U).^2);
    RR=sqrt(rn^2+(rn*tan(thetaT)+xn-v*U).^2);
    R=RT+RR;
    DT=T-R/c;
    phase=pi*Kr*DT.^2-2*pi/lambda*R;
    s_ut=s_ut+sigma*exp(j*phase).*(abs(DT)<Tr/2).*(abs(v*U-xn)<Lsar/2);
end;

%%距離壓縮
p0_t=exp(j*pi*Kr*(t-2*Rc/c).^2).*(abs(t-2*Rc/c)<Tr/2);%距離向LFM信號
p0_f=fftshift(fft(fftshift(p0_t)));
s_uf=fftshift(fft(fftshift(s_ut)));%距離向FFT
src_uf=s_uf.*(conj(p0_f).'*ones(1,Na));%距離壓縮
src_ut=fftshift(ifft(fftshift(src_uf)));%距離壓縮后的信號
src_fut=fftshift(fft(fftshift(src_ut).')).';%距離多普勒域

%%二次距離壓縮,距離遷移校正原理仿真
src_fuf=fftshift(fft(fftshift(src_uf).')).';%距離壓縮后的二維頻譜
F=f'*ones(1,Na);%擴充為矩陣
FU=ones(Nr,1)*fu;
p0_2f=exp(j*pi/fc^2/fdr*(FU.*F).^2+j*pi*fdc^2/fc/fdr*F-j*pi/fc/fdr*FU.^2.*F);
s2rc_fuf=src_fuf.*p0_2f;
s2rc_fut=fftshift(ifft(fftshift(s2rc_fuf)));%距離多普勒域

%%方位壓縮
p0_2fu=exp(j*pi/fdr*(FU-fdc).^2);%方位向壓縮因子
s2rcac_fut=s2rc_fut.*p0_2fu;%方位壓縮
s2rcac_fuf=fftshift(fft(fftshift(s2rcac_fut)));%距離方位壓縮后的二維頻譜
s2rcac_ut=fftshift(ifft(fftshift(s2rcac_fut).')).';%方位向IFFT


subplot(221)
G=20*log10(abs(s_ut)+1e-6);
gm=max(max(G));
gn=gm-40;%顯示動態范圍40dB
G=255/(gm-gn)*(G-gn).*(G>gn);
imagesc(x,r-Rc,-G),colormap(gray)
grid on,axis tight,
xlabel('Azimuth')
ylabel('Range')
title('(a)原始信號')

subplot(222)
G=20*log10(abs(src_fut)+1e-6);
gm=max(max(G));
gn=gm-40;%顯示動態范圍40dB
G=255/(gm-gn)*(G-gn).*(G>gn);
imagesc(fu,r-Rc,-G),colormap(gray)
grid on,axis tight,
xlabel('Azimuth')
ylabel('Range')
title('(b)距離多普勒域頻譜')

subplot(223)
G=20*log10(abs(s2rc_fut)+1e-6);
gm=max(max(G));
gn=gm-40;%顯示動態范圍40dB
G=255/(gm-gn)*(G-gn).*(G>gn);
imagesc(fu,r-Rc,-G),colormap(gray)
grid on,axis tight,
xlabel('Azimuth')
ylabel('Range')
title('(c)RMC后的RD域頻譜')


subplot(224)
G=20*log10(abs(s2rcac_ut)+1e-6);
gm=max(max(G));
gn=gm-60;%顯示動態范圍40dB
G=255/(gm-gn)*(G-gn).*(G>gn);
imagesc(x,r-Rc,G),colormap(gray)
grid on,axis tight,
xlabel('Azimuth')
ylabel('Range')
title('(d)目標圖象')