%回波信号仿真
%parameter
c=3e8; %光速
fc=5.3e9;
lambda=c/fc;
%SAR基本参数
%SAR平台
v=100; %雷达平台水平速度
H=5000;%高度5000m
%测绘带数据
xmin=-100;
xmax=100;
yc=1500;
y0=200;
thita=0;
R0=sqrt(yc^2+H^2);%目标与SAR垂直斜距
D=3; %天线等效孔径
La=3;
Lsar=lambda*R0/D/cos(thita)^2;%SAR长度
Tsar=Lsar/v;
%方位向参数
Ka=2*v^2*cos(thita)^3/lambda/R0;%调频率
BW=Ka*Tsar;%带宽
PRF=2.0*BW;%脉冲频率
PRT=1/PRF;%脉冲时间
ds=PRT;
Nslow=ceil((xmax-xmin+Lsar)/v/ds);
Nslow=2^nextpow2(Nslow);
sn=linspace((xmin-Lsar/2)/v,(xmax+Lsar/2)/v,Nslow);
PRT=(xmax-xmin+Lsar)/v/Nslow;
PRF=1/PRT;
ds=PRT;
%距离向参数
Tr=2e-6;
Br=100e6; %脉冲信号带宽
Kr=Br/Tr;
Fsr=1.2*Br; %距离向采样频率
dt=1/Fsr;
Rmin=sqrt((yc-y0)^2+H^2);
Rmax=sqrt((yc+y0)^2+H^2);
Nfast=ceil(2*(Rmax-Rmin)/c/dt+Tr/dt);
tm=linspace(2*Rmin/c-Tr/2,2*Rmax/c+Tr/2,Nfast);
dt=(2*Rmax/c+Tr-2*Rmin/c)/Nfast;
f0=0;
Fsr=1/dt;
%分辨率
DY=c/2/Br;%距离向
DX=D/2;%方位向
%目标点
Ntarget=1;
Ptarget=[(xmin+xmax)/2, yc, 1 %目标位置
xmin, yc-y0, 1
xmin yc+y0 1
xmax yc-y0 1
xmax yc+y0 1
1/3*xmin yc-y0 1
1/3*xmax yc-y0 1
2/3*xmin yc-y0 1
2/3*xmax yc+y0 1
1/3*xmin yc-1/3*y0 1
1/3*xmin yc-2/3*y0 1
1/3*xmin yc+1/3*y0 1
1/3*xmin yc+2/3*y0 1
1/3*xmax yc-1/3*y0 1
1/3*xmax yc-2/3*y0 1
1/3*xmax yc+1/3*y0 1
1/3*xmax yc+2/3*y0 1
2/3*xmin yc-1/3*y0 1
2/3*xmin yc-2/3*y0 1
2/3*xmin yc+1/3*y0 1
2/3*xmin yc+2/3*y0 1
2/3*xmax yc-1/3*y0 1
2/3*xmax yc-2/3*y0 1
2/3*xmax yc+1/3*y0 1
2/3*xmax yc+2/3*y0 1];
M=Nfast;
N=Nslow;
K=Ntarget;
T=Ptarget;
Srnm=zeros(N,M);
for k=1:1: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);
phase=pi*Br/Tr*Dfast.^2-(4*pi/lambda)*(R'*ones(1,M));%相位
ssssssss=-Tr/2<Dfast&Dfast<Tr/2;
Srnm=Srnm+sigma*exp(1i*phase).*ssssssss.*((abs(Dslow)<Lsar/2)'*ones(1,M));
end
tr=tm-2*Rmin/c;
Refr=exp(j*pi*Kr*tr.^2).*(-Tr/2<tr&tr<Tr/2);
Sr=ifft((fft(Srnm')'.*(ones(N,1)*conj(fft(Refr')')))');
Gr=abs(Sr);
ta=sn-xmin/v;
Refa=exp(j*pi*Ka*ta.^2).*(abs(ta)<Tsar/2);
Sa=ifft(fft(Sr).*((conj(fft(Refa)).'*ones(1,M)))');
Ga=abs(Sa);
%w_fast=zeros(1,512);
%for ii=1:512;
% ss=tm(1,ii);
%w_fast(1,ii)=rectpuls(ss/Tr);
%end
%W_fast=fft(w_fast);
%pa_thita=sinc(La*thita/lambda);
%yitac=0;
%w_slow=(sinc(atan(v*(sn-yitac)/R0))).^2;
%W_slow=fft(w_slow);
uuu=(atan(v*sn/R0)*La/lambda);%时域方位向天线方向图中间变量
pa=sinc(uuu).^2;
Srnm=Srnm.*(pa'*ones(1,Nfast));
figure;
title('点目标仿真回波信号');
imagesc(abs(Srnm));
colormap(gray);
Srnm1=fft2(Srnm);
imagesc(abs(Srnm1));
colormap(gray);
Srnm1=fftshift(Srnm1);
Ka=2*v^2*(cos(thita))^3/lambda/R0;
fyita=linspace(-PRF/2,PRF/2,Nslow); %获得慢时间频率变量;在(-Fslow/2,Fslow/2)范围;
ftao= linspace(-Fsr/2,Fsr/2,Nfast); %获得快时间频率变量;在(-Ffast/2,Ffast/2)范围;
Theta_range_comp = ones(Nslow,1)*(pi/Kr*ftao.^2); %用于距离压缩的相位
Theta_sqrt_1 = ones(Nslow,1) * ((fc+ftao).^2); %根号内第一项
Theta_sqrt_2 = (c^2*fyita.^2/4/v^2).'*ones(1,Nfast); %根号内第二项
Theta_sqrt = 4*pi*R0/c*sqrt(Theta_sqrt_1-Theta_sqrt_2); %根号项
Theta_2d_ref = Theta_range_comp + Theta_sqrt; %获得总相位
Ref_2d= exp(1i*Theta_2d_ref); %获得参考函数
Srnm2= Ref_2d.*Srnm1; %参考函数相乘
imagesc(abs(Srnm2));
Delta_ffast = sqrt(ones(Nslow,1)*(fc + ftao).^2 + (c^2*fyita.^2/4/v^2).'*ones(1,Nfast)) - ones(Nslow,1) * (ftao+fc);
Delta_ffast = Delta_ffast * Nfast/Fsr;%归一化
Nd = max(max(fix(Delta_ffast)));
temp = zeros(Nslow,Nfast); %临时数据矩阵temp;
Ni = 16/2; %sinc插值核长度Ni=8;
for i = 1 : 1 : Nslow
data = Srnm2(i,:);
for m = Ni+Nd : 1 : Nfast-Nd-Ni
Insert_AMT = Delta_ffast(i,m);
AMT = Insert_AMT; %获得插值量(位移量);
AMT_N = fix(AMT); %取插值位移量的整数部分
AMT_F = AMT - AMT_N; %取插值位移量的小数部分
NN = AMT_N + (-Ni+1:Ni); %所取的插值位置
NN1 = AMT_F - (-Ni+1:Ni); %所取的sinc核位置
%% 进行sinc核数据插值计算
sinc_kernal = sinc(NN1); %sinc核函数的值;
temp(i,m) = sinc_kernal * data(1,m+NN).';%重建方程
end
end
Srnm3 = temp;
Wfast=kaiser(Nfast,2.5);
Wfast=ones(Nslow,1)*(Wfast.');
Srnm3=Srnm3.*Wfast;
imagesc(abs(Srnm3));
Srnm4 = ifft2(Srnm3);
colormap(gray);
imagesc(255-abs(Srnm4));
figure;
title('Image after 2-d ifft');
xlabel('距离向');ylabel('方位向');
[a,b]=max(Srnm4);
[x,y]=max(a);
z=b(y);
Nslice=32;
times=16;
slice_2D=zeros(Nslice,Nslice);
slice_2D(1:Nslice,1:Nslice)=Srnm4(z-Nslice/2:1:z+Nslice/2-1,y-Nslice/2:1:y+Nslice/2-1);
slice_2D_fft=fft2(slice_2D);
imagesc(abs(slice_2D_fft));
slice_2D_interp=zeros(Nslice*times,Nslice*times);
slice_2D_interp(1:Nslice,1:Nslice)=slice_2D_fft;
slice_2D_ifft=ifft2(fftshift(slice_2D_interp));
slice=abs(slice_2D_ifft);
figure(3)
[a1,b1]=max(slice);
[x1,y1]=max(a1);
z1=b1(y1);
plot(20*log10(slice(:,y1)/max(slice(:,y1))));
xlabel('方位向');
ylabel('幅度（db）');
title('方位向剖面图');
figure(4)
plot(20*log10(slice(z1,:)/max(slice(z1,:))));
xlabel('距离向');
ylabel('幅值（db）');
title('距离向剖面图');
