clc;clear;
c=3e8;%光速
fc=1.5e9;%载频
lambda=c/fc;%波长
%%测绘带区域
X0=200;%方位向[-X0,X0]
Rc=3000;
R0=150;%距离向[Rc-R0,Rc+R0]
%%距离向(Range),r/t
Tr=1.5e-6;%脉宽 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),x/u
v=300;%SAR 平台速度
Lsar=300;%合成孔径长度
Na=1024;
x=linspace(-X0,X0,Na);
u=x/v;%u域序列
du=2*X0/v/Na;
fu=linspace(-1/2/du,1/2/du,Na);%fu域序列
fdc=0;%Doppler调频中心
fdr=-2*v^2/lambda/Rc;%Doppler调频斜率
%%目标位置
Ntar=6;%目标个数
Ptar=[Rc , 0 ,1 %距离向坐标,方位向坐标,目标截面积RCS sigma
Rc+50 , -50 ,1
Rc+50 , 50 ,1
Rc+50 , -150,1
Rc+50 , 150,1
Rc+100 , 0 ,1];
%%产生回波
s_ut=zeros(Nr,Na);
U=ones(Nr,1)*u;%扩充为矩阵
T=t'*ones(1,Na);
for i=1:Ntar
rn=Ptar(i,1);xn=Ptar(i,2);sigma=Ptar(i,3);
if i==1
R=sqrt((rn+3*U).^2+(xn+30*U-v*U).^2);
else
R=sqrt(rn^2+(xn-v*U).^2);
end
DT=T-2*R/c;
phase=pi*Kr*DT.^2-2*pi/lambda*R*2;
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);%距离向匹配函数(行向量)
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).')).';%方位向FFT 距离-多普勒域
%距离弯曲校正(二维去耦)
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)));%距离向IFFT 距离-多普勒域
%方位压缩
p0_2fu=exp(j*pi/fdr*FU.^2);%方位向压缩因子
s2rcac_fut=s2rc_fut.*p0_2fu;%方位压缩
s2rcac_fuf=fftshift(fft(fftshift(s2rcac_fut)));%距离方位压缩后的二维频谱
s2rcac_ut=fftshift(ifft(fftshift(s2rcac_fut).')).';%方位向IFFT
figure(1)
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)RD域频谱')
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)目标图象')