%% 多点目标 单站SAR
% 2014/03/17
clc;
clear all;
close all;
%% parameters
c = 3e8;
j = sqrt(-1);
pi = 3.1416;
fc = 1e9;
lamda = c/fc;
D = 4;
Va = 150;
Kr = 20e12;
Tr = 2.5e-6;
sq_ang = 0/180*pi; %斜视角
Br = Kr*Tr;
Frfactor = 1.2;
Fr = Br*Frfactor; %距离向采样率
Ba = 0.886*2*Va*cos(sq_ang)/D;
Fafactor = 1.2;
Fa = Ba*Fafactor; %方位向采样率
R_near = 2e4;
R_far = R_near + 1000;
La_near = 0.886*R_near*lamda/cos(sq_ang)^2/D; % 近距合成孔径长度
La_far = 0.886*R_far *lamda/cos(sq_ang)^2/D;
Tc_near = -R_near*tan(sq_ang)/Va; % 波束中心时间
Tc_far = -R_far *tan(sq_ang)/Va;
fdc = 2*Va*sin(sq_ang)/lamda; % 多普勒中心频率
Y_min = Va*Tc_far;
Y_max = Y_min+100;
Rmin = sqrt(R_near^2+(Tc_near*Va+La_near/2)^2);
Rmax = sqrt(R_far^2 +(Tc_far*Va -La_far/2)^2);
disp('parameters:');
disp('minimal slant range');disp(Rmin);
disp('maximal slant range');disp(Rmax);
disp('range resolution:');disp(0.886*(c/2/Br));
disp('azimuth resolution:');disp(0.886*(Va/Ba));
disp('doppler centroid frequency');disp(fdc);
%% echo model
Nr = (2*Rmax/c - 2*Rmin/c + Tr)*Fr;
Nr = 2^nextpow2(Nr);
tr = linspace(-Tr/2 + 2*Rmin/c,Tr/2 + 2*Rmax/c,Nr); % 快时间t采样点，最早收到回波~最晚收到回波
Fr = (Nr - 1)/(Tr/2 + 2*Rmax/c - (-Tr/2 + 2*Rmin/c));
Na = ((Tc_near + La_near/2/Va)-(Tc_far - La_far/2/Va))*Fa;
Na = 2^nextpow2(Na);
ta = linspace(Tc_far - La_far/2/Va,Tc_near + La_near/2/Va,Na);% 慢时间采样点
Fa = (Na -1)/(Tc_near + La_near/2/Va -(Tc_far - La_far/2/Va));
Rpt = [R_near R_near+500 R_near+1000]; % 点目标距离坐标
Ypt = [0 0 0]; % 点目标方位坐标
La = 0.886*Rpt*lamda/(cos(sq_ang)^2)/D; % 点目标合成孔径长度
Tc = -Rpt*tan(sq_ang)/Va;
Npt = length(Rpt);
Y_high = max(Ypt) + 50;
Y_low = min(Ypt) - 50;
R_left = R_near - 50;
R_right = R_far + 50;
disp('number of point targets:');disp(Npt);
disp('range sample number:');disp(Nr);
disp('azimuth sample number:');disp(Na);
disp('range sample rate:');disp(Fr);
disp('azimuth sample rate:');disp(Fa);
sig = zeros(Na,Nr);
for k = 1:Npt
delay = 2/c*sqrt(Rpt(k)^2+(Ypt(k)-ta*Va).^2);
Dr = ones(Na,1)*tr - delay'*ones(1,Nr); % t-tao
sig = sig + exp(j*pi*Kr*Dr.^2 - j*2*pi*fc*delay'*ones(1,Nr))...
.*(abs((ta - Ypt(k)/Va - Tc(k))'*ones(1,Nr)) <= La(k)/2/Va)...
.*(abs(Dr) <= Tr/2);
end
figure('Name','回波信号幅度');imagesc(abs(sig));
title('回波信号幅度');xlabel('距离向');ylabel('方位向');
figure('Name','回波相位');imagesc(abs(angle(sig)));
title('回波信号相位');xlabel('距离向');ylabel('方位向');
%% BP成像
% 步骤一：距离向匹配滤波
sig_rd = fft(sig,[],2);
fr = -1/2:1/Nr:(1/2-1/Nr);
fr = fftshift(fr*Fr);
frrr=linspace(-Fr/2,Fr/2,Nr);
filter_r = ones(Na,1)*exp(j*pi*fr.^2/Kr);
sig_rd = sig_rd.*filter_r;
nup = 2;
% % nup = 1; %%不进行升采样也是可以的
Nr_up = Nr*nup;
nz = Nr_up - Nr;
dtr = 1/nup/Fr;
sig_rd_up = [sig_rd(:,1:Nr/2),zeros(Na,nz),sig_rd(:,(Nr/2+1):Nr)];
sig_rdt = ifft(sig_rd_up,[],2);
% % sig_rdt = ifft(sig_rd,[],2);
figure('Name','距离向匹配滤波之后');imagesc(abs(sig_rdt));
title('距离向匹配滤波之后');xlabel('距离向');ylabel('方位向');
%步骤二：对成像区域进行网格化
R = zeros(1,Nr);
for ii = 1:Nr
R(1,ii) = R_left + (R_right - R_left)/(Nr-1)*(ii-1);
end
R = ones(Na,1)*R;
Y = zeros(1,Na);
for ii = 1:Na
Y(1,ii) = Y_low + (Y_high - Y_low)/(Na-1)*(ii-1);
end
Y = Y'*ones(1,Nr);
%步骤三：根据时延在回波域寻找相应位置并进行成像
f_back = zeros(Na,Nr);
for ii = 1:Na
R_ij = sqrt(R.^2 + (Y-Va*ta(ii)).^2);
t_ij = 2*R_ij/c;
t_ij = round((t_ij - (2*Rmin/c - Tr/2))/dtr); %t_ij对应中心处，距离压缩后幅度最大的位置
it_ij = (t_ij>0 & t_ij<=Nr_up);
t_ij = t_ij.*it_ij + Nr_up*(1-it_ij);
sig_rdta = sig_rdt(ii,:);
sig_rdta(Nr_up) = 0;
f_back = f_back + sig_rdta(t_ij).*exp(1i*4*pi*R_ij/lamda);
end
figure('Name','BP算法处理后');imagesc(abs(f_back));
title('BP算法处理后');xlabel('距离向');ylabel('方位向');
