%% CSA算法 正侧视
%% =================================
clear;clc;close all;
%%=================================
% 参数设置
C=3e8; %光速
Fc=1e9; %载频
lambda=C/Fc; %波长
Xmin=-50; %方位向最小距离
Xmax=50; %方位向最大距离
Yc=10000; %观测带中心距离
Y0=500; %观测带地距宽度
V=100; %方位向速度
H=5000; %飞机高度
R0=sqrt(Yc^2+H^2); %观测带中心到飞机的位置
% RNear=sqrt((Yc-Y0)^2+H^2);
% RFar=sqrt((Yc+Y0)^2+H^2);
% 天线参数
D=4; %实际孔径大小
Lsar=lambda*R0/D; %单孔径方位分辨率
Tsar=Lsar/V; % 一个合成孔径时间
%%=================================
% 方位维参数
Ka=-2*V^2/lambda/R0; %方位调频率
Ba=abs(Ka*Tsar); %方位向带宽
PRF=Ba; %理想方位向采样率
PRT=1/PRF; %理想方位向采样间隔
display(PRF);
ds=PRT;
Nslow=ceil((Xmax-Xmin+Lsar)/V/ds);
Nslow=2^nextpow2(Nslow);%为了方便FFT
sn=linspace((Xmin-Lsar/2)/V,(Xmax+Lsar/2)/V,Nslow);
PRT=(Xmax-Xmin+Lsar)/V/Nslow;%实际采样间隔
PRF=1/PRT;%实际采样率
display(PRF);
ds=PRT;%实际采样间隔
%%=================================
% 距离维参数
Tr=5e-6; %调频周期
Br=30e6; %调频带宽
Kr=Br/Tr; %调频率
Fsr=3*Br; %采样率
dt=1/Fsr; %采样时间
Rmin=sqrt((Yc-Y0)^2+H^2);
Rmax=sqrt((Yc+Y0)^2+H^2+(Lsar/2)^2);
Nfast=ceil(2*(Rmax-Rmin)/C/dt+Tr/dt);
Nfast=2^nextpow2(Nfast);
tm=linspace(2*Rmin/C-Tr/2,2*Rmax/C+Tr/2,Nfast);
dt=(2*Rmax/C+Tr-2*Rmin/C)/Nfast; %实际采样时间
Fsr=1/dt;%实际采样频率
RNear=Rmin-Tr/2*C/2;%近距
RFar=Rmax+Tr/2*C/2;%远距
%%=================================
% 分辨率
DY=C/2/Br;
DX=D/2;
display(DY);
display(DX);
%%=================================
% 目标参数
Ntarget=3;
%format [x, y, reflectivity]
Ptarget=[(Xmin+Xmax)/2,Yc,1 %方位向 距离向 反射系数
(Xmin+Xmax)/2,Yc-Y0/2,1
(Xmin+Xmax)/2,Yc+Y0/2,1];
%%=================================
% 产生回波数据
T=Ptarget;
Srnm=zeros(Nslow,Nfast);
for k=1:1:Ntarget
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(Nslow,1)*tm-tau'*ones(1,Nfast);%tau-2*R/c
phase=pi*Kr*Dfast.^2-(4*pi/lambda)*(R'*ones(1,Nfast));%相位
Srnm=Srnm+sigma*exp(1i*phase).*(-Tr/2<Dfast&Dfast<Tr/2).*((abs(Dslow)<Lsar/2)'*ones(1,Nfast));
end
figure;
mesh(abs(Srnm));
view(2);
xlabel('距离向');ylabel('方位向');title('时域图像')
axis tight
%%
%%=================================
% 傅立叶变换，转换到距离多普勒域(FFT思路值得参考)
SrnmTKx=zeros(Nslow,Nfast);
for k = 1:Nfast
SrnmTKx(:,k) = fftshift(fft(Srnm(:,k)));
end
figure;
mesh(abs(SrnmTKx));
view(2); %%从上面往下看
xlabel('距离向');ylabel('方位向');title('方位向FFT图像')
axis tight
%%
DR = linspace(RNear,RFar,Nfast);
dx = ds*V;
Kx = linspace(-1/2,1/2,Nslow)*PRF;%fn
W = linspace(-1/2,1/2,Nfast)*Fsr;%ftau
D_Kx_V = sqrt(1-(C*Kx).^2/(4*V^2*Fc^2));%徙动因子D
%Km = Kr./(1-(Kr*C*R0*Kx.^2)./(2*V^2*Fc^3*D_Kx_V.^3)); % 新的调频率
Km = Kr./(1-(Kr*C*(ones(Nslow,1)*DR).*((Kx.^2).'*ones(1,Nfast)))./(2*V^2*Fc^3*(D_Kx_V.^3).'*ones(1,Nfast))); % 新的调频率
tmN = ones(Nslow,1)*tm - (2*R0./(C*D_Kx_V)).'*ones(1,Nfast);%参考距离是R0
Kxref = 0;%由于参考距离是R0 因此徙动因子是1 可知方位频率Kxref=0
D_Kxref_V = sqrt(1-(C*Kxref)^2/(4*V^2*Fc^2));
% 与变标方程进行相乘
%SrnmTKx = SrnmTKx.*exp(1i*pi*Km.'*ones(1,Nfast).*(D_Kxref_V./(D_Kx_V.'*ones(1,Nfast))-1).*tmN.^2);
SrnmTKx = SrnmTKx.*exp(1i*pi*Km.*(D_Kxref_V./(D_Kx_V.'*ones(1,Nfast))-1).*tmN.^2);
Hf=exp(-j*pi*Kr*Dfast.^2).*(-Tr/2<Dfast&Dfast<Tr/2).*((abs(Dslow)<Lsar/2)'*ones(1,Nfast));
% sr_cmps=ifty(fty(Hf).*fty(SrnmTKx1));
% sr_cmp=ifty(fty(Hf).*fty(SrnmTKx));
% figure
% mesh(abs(sr_cmps)),view(2);axis tight
% figure
% mesh(abs(sr_cmp)),view(2);axis tight
% mesh(abs(SrnmTKx1)),view(2);axis tight
%%
% 距离维匹配滤波
for n = 1:Nslow
SrnmTKx(n,:) = fftshift(fft(SrnmTKx(n,:)));
end
Refr = exp(1i*pi*(ones(Nslow,1)*W).^2.*(D_Kx_V.'*ones(1,Nfast)./D_Kxref_V./Km ) );
SrnmTKx = SrnmTKx.*Refr;
figure;
mesh(abs(SrnmTKx));
view(2);
xlabel('距离向');ylabel('方位向');title('二维频域图像')
axis tight
%%
% 一致RCM
RCMbulk = exp( 1i*4*pi/C*R0*ones(Nslow,1)*W.*( (1./D_Kx_V-1./D_Kxref_V).'*ones(1,Nfast)) );
%RCMbulk = exp( 1i*4*pi/C*(ones(Nslow,1)*(W.*DR)).*( (1./D_Kx_V-1./D_Kxref_V).'*ones(1,Nfast)) );
SrnmTKx = SrnmTKx.*RCMbulk;
figure;
mesh(abs(SrnmTKx));
view(2);
xlabel('距离向');ylabel('方位向');title('一致RCM图像')
axis tight
%%
for n = 1:Nslow
SrnmTKx(n,:) = ifft(fftshift(SrnmTKx(n,:)));
end
% 方位维相位补偿
Dphase1 = 4*pi*Fc*(ones(Nslow,1)*DR).*(D_Kx_V.'*ones(1,Nfast))./C;%距离相关的方位调制
%Dphase2 = -4*pi/C^2*Km.'*ones(1,Nfast).*(1-(D_Kx_V./D_Kxref_V).'*ones(1,Nfast)).*( (ones(Nslow,1)*DR)./(D_Kx_V.'*ones(1,Nfast)) - R0./(D_Kx_V.'*ones(1,Nfast))).^2;
Dphase2 = -4*pi/C^2*Km.*(1-(D_Kx_V./D_Kxref_V).'*ones(1,Nfast)).*( (ones(Nslow,1)*DR)./(D_Kx_V.'*ones(1,Nfast)) - R0./(D_Kx_V.'*ones(1,Nfast))).^2;
%附加相位补偿
SrnmTKx = SrnmTKx.*exp(1i*Dphase1).*exp(1i*Dphase2);
figure;
mesh(abs(SrnmTKx));
view(2);
xlabel('距离向');ylabel('方位向');title('方位FFT图像')
axis tight
%%
for n = 1:Nfast
SrnmTKx(:,n) = ifft(SrnmTKx(:,n));
end
figure;
mesh(abs(SrnmTKx));
view(2);
view(2);
xlabel('距离向');ylabel('方位向');title('最后图像')
axis tight
