%%产生Stripmap SAR的回波
clear all
c=3e8;%光速
%%测绘带区域
X0=204.8;%方位向[-X0,X0],条带长度的“一半”
Rc=10000;%条带中心参考距离
R0=480;%距离向[Rc-R0,Rc+R0],条带宽度的“一半”
lambda=0.02;%波长
fc=c/lambda;%载频
%%距离向(Range),r/t domain
Tr=40e-6;%LFM信号脉宽
Br=10e6; %LFM信号带宽(=Kr*Tr),Hz
Kr=Br/Tr; %调频斜率
Fsr=40e6;%径向快时间采样率,应>Br
%还应尽量使下面的Nr=2若干次幂,以利于FFT。
Nr=ceil(Fsr*4*R0/c);
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
Va=200;%SAR 平台速度
Lsar=150;%合成孔径长度
PRF=1000;%Hz,应>多普勒带宽Ba(k*Tsar),还应尽量使下面的Na=2若干次幂,以利于FFT
Ta=1/PRF;
D=Va*Ta;%两孔径间距
Na=ceil(PRF*2*X0/Va);
x=linspace(-X0,X0,Na);%u域序列
u=x/Va;
du=2*X0/Va/Na;
fu=linspace(-1/2/du,1/2/du,Na);%fu域序列
fdc=0;%Doppler调频中心频率
fdr=-2*(Va)^2/lambda/Rc;%Doppler调频斜率
Ba=abs(fdr*Lsar/Va);%计算多普勒带宽
%目标位置
Ptar=[Rc,0,0,0,1 %距离向坐标,方位向坐标,距离向速度,方位向速度,sigma
Rc,30,10,0,1];
Ntar=size(Ptar,1);%目标个数
%%%%%%%%%%%%%%% 1、产生距压、徙校后的回波 %%%%%%%%%%%
s_ut=zeros(Nr,Na);
s_ut2=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);
vrn=Ptar(i,3);vcn=Ptar(i,4);
sigma=Ptar(i,5);
R=sqrt(((Va-vcn)*(U-xn/Va)).^2+(rn+vrn*(U-xn/Va)).^2);%位于x=0处的动目标波程,
phase=-4*pi*R/lambda;
tmp=zeros(size(T));tmp(128,:)=ones(1,size(T,2));% 距离向上本应是乘Sinc函数,但可用delta函数近似,即tmp,tmp在128处,故目标距离均应设在Rc处
s_ut=s_ut+sigma*exp(j*phase).*(abs(Va*U-xn)<Lsar/2).*tmp; % .*sinc(Br*(T-2*rn/c));
R2=sqrt(((Va-vcn)*(U-xn/Va)+Va*Ta-D-vrn*Ta).^2+(rn+vrn*(U-xn/Va)+vrn*Ta).^2 );
phase2=-4*pi*R2/lambda;
s_ut2=s_ut2+sigma*exp(j*phase2).*(abs(Va*U-xn)<Lsar/2).*tmp;
end;
%%%%%%%%%%%%%%% 第一路成像 %%%%%%%%%%%%%
FU=ones(Nr,1)*fu;
p0_2fu=exp(j*pi/fdr*(FU-fdc).^2);%方位向压缩因子 %%% .*(abs(FU-fdc)<Ba/2)给匹配函数加窗
% p0_2fu=exp(j*2*pi*Rc/Va*(4*Va*Va/lambda/lambda-FU.^2).^0.5);%采用式(5.5.2)
s_fut=fftshift(fft(fftshift(s_ut).')).';%距离多普勒域
scc_fut=s_fut.*p0_2fu;%方位压缩
scc_ut=fftshift(ifft(fftshift(scc_fut).')).';%方位向IFFT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% 第二路成像 %%%%%%%%%%%%%
s_fut2=fftshift(fft(fftshift(s_ut2).')).';%距离多普勒域
scc_fut2=s_fut2.*p0_2fu;%方位压缩
scc_ut2=fftshift(ifft(fftshift(scc_fut2).')).';%方位向IFFT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% DPCA对消 %%%%%%%%%%%%%%%%
s_dpca=abs(scc_ut-scc_ut2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% 画图 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% figure;
% 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)原始信号')
%
% figure;
% G=20*log10(abs(scc_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('(b)第一路成像')
%
% figure;
% G=20*log10(abs(scc_ut2)+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('(c)第二路成像')
figure;
plot(u,real(scc_ut(128,:)),'-',u,real(scc_ut2(128,:)),':',u,real(scc_ut2(128,:)-scc_ut(128,:)),'r');
legend('第一路成像结果','第二路成像结果','DPCA对消后结果');
figure;
subplot(311);
plot(u,real(scc_ut(128,:)));
title('第一路成像结果');
subplot(312);
plot(u,real(scc_ut2(128,:)));
title('第二路成像结果');
subplot(313);
plot(u,real(scc_ut2(128,:)-scc_ut(128,:)));
title('DPCA对消后结果');