CS（Chirp Scaling）算法通过相位相乘的形式避免了插值操作，相比于RD算法极大的提高了SAR的成像效率。
仿真代码如下
close all; clear; clc;
c=3e8;%光速
%% 仿真参数
theta = 45; %下视角
thetac = 0; %斜视角
pd = 3; %地距分辨力为3mxy均为3m
pr0 = 3; %方位向分辨力
pr = pd*cosd(theta); %距离分辨力c/2B=>B>c/2pr
h = 40e3; %飞行高度h=10e3;
gr = h*tand(theta); %地距
Rc = h/cosd(theta)/cosd(thetac); %中心穿越时刻距离（距离场景中心距离，对于波束内部的其他点有所差异）
% widthr = 1000; %测绘带宽为1000m,y方向
% widtha = 1000; %方位向带宽x方向1000m
vs = 7500;
vr = vs/1.06; %卫星比等效速度大6% 参考第90页vr=250;%雷达有效速度
vg = vr/1.06;
vr=vs;
% Rc=30e3;%景中心斜距
Tp=10e-6;%脉冲宽度
B=150e6;%信号带宽
Kr=B/Tp;%距离脉冲调频率
fs=1.2*B;%距离采样率
fc=10e9;%载波频率
lambda=c/fc;%波长
% Ka=2*vr^2/lambda/Rc;%方位向调频率
la=2*pr0;%天线真实孔径
fdwidth=0.886*2*vr/la;%多普勒带宽
prf = 5e3; %脉冲重复频率
im=sqrt(-1);%虚数单位
thetac=0;%斜视角度，正侧视
fd0=2*vr*sin(thetac*pi/180)/lambda;%参考中心频率，取多普勒中心频率f_etac 4.34
ls=0.886*Rc*lambda/la;%合成孔径长度
Ta=ls/vr;%目标照射时间
%% 成像区域[Xc-X0,Xc+X0; Yc-Y0,Yc+Y0]
Xc = sqrt(Rc^2-h^2);
Yc = 0;
Xo = 500;
Yo = 500;
Rmin=sqrt(h^2+(Xc-Xo)^2);%观测场景距飞机的最近距离
Rmax=sqrt(h^2+(Xc+Xo)^2);%观测场景距飞机的最远距离
Ra=ls+2*Yo;%正侧视时雷达在方位向行走距离
centerpos=[gr,0,0];
%% 目标位置
target=[centerpos;gr+50,centerpos(2)+50,0;gr-30,centerpos(2)-20,0;gr+50,centerpos(2),0;gr+150,centerpos(2),0]; %点目标位置(第一个为方位向，第二个为 距离向)
widthr = 1000;widtha=1000;
xylim = [gr-widthr/2,gr+widthr/2;-widtha/2,widtha/2;0,0];%成像区域范围
figure;plot(target(:,1),target(:,2),'r.','Markersize',10);
xlim([xylim(1,1),xylim(1,2)]);ylabel('方位向距离/m');xlabel('距离向距离/m'); ylim([xylim(2,1),xylim(2,2)]);title('目标位置')
%% 生成回波
tm=-Ra/vr/2:1/prf:Ra/vr/2-1/prf;%慢时间轴
tk=2*Rmin/c-Tp/2:1/fs:2*Rmax/c+Tp/2-1/fs;%快时间轴
r=sqrt((tk*c/2).^2-h^2);%距离向横坐标
x=r;
nm=length(tm);%方位向采样点数
nk=length(tk);%距离向采样点数
echoall=zeros(nm,nk);%回波
y=vr*tm;%飞机的位置(注意：慢时间轴不同时表达式不同)
Rtm=zeros(size(target,1),nm);
A0=1;%幅度
bar = waitbar(0,"回波生成...");
radar = [zeros(length(tm),1),y',h*ones(length(tm),1)]; %卫星坐标
% Rtm1 = sqrt(sum((target(1,:)-[0,y(280),h]).^2)); plot(rectpuls(tk-2*Rtm1/c-Tp/2,Tp));
for k=1:size(target,1) %目标数
echo=[];
Rtm = sqrt((target(k,1).*ones(length(tm),1)-radar(:,1)).^2+...
(target(k,2).*ones(length(tm),1)-radar(:,2)).^2+...
(target(k,3).*ones(length(tm),1)-radar(:,3)).^2)*ones(1,length(tk));
echo(:,:,k)=(abs(target(k,2)-y')*ones(1,length(tk))<ls/2).*(abs((ones(length(tm),1)*tk-2.*Rtm/c)/2)<=Tp/2).*...
exp(1j*2*pi*fc.*(-2*Rtm./c)+1j*pi*Kr*(tk-2.*Rtm./c).^2); %回波模型(停走模型)
waitbar(k/size(target,1),bar)
echoall = echo(:,:,k) + echoall;
end
close(bar)
% for i=1:size(target,1)
tk_mat = ones(nm,1)*tk; %距离时间轴矩阵
ta_mtx = tm.'*ones(1,nk); %方位时间轴矩阵
ft_mat = ones(nm,1)*linspace(-fs/2,fs/2,nk); %距离频率轴矩阵
fdc=2*vr*sin(thetac*pi/180)/lambda;%参考中心频率，取多普勒中心频率f_etac 4.34
fd1=fdc+linspace(-prf/2,prf/2,nm).'*ones(1,nk);%方位频率轴矩阵
f=linspace(-fs/2,fs/2,nk);
HFk = rectpuls(f,Kr*Tp).*exp(-1j*pi*(f+0).^2./Kr); %PSP计算已经下变频不需要-fc,若以及移动Tp/2要注意移动+B/2(移动Tp/2)
fd = fdc+linspace(-prf/2,prf/2,nm);
figure,imagesc(tk,tm,abs(echoall)) ;title('原始回波幅度')
figure,imagesc(tk,tm,angle(echoall)) ;title('原始回波相位')
caxis([0 1]);axis xy
%% 方位向傅里叶变换
echoallrd=zeros(nm,nk);
% for i=1:nk
% echoallrd(:,i)= fftshift(fft((echoall(:,i)),nm));%转换到多普勒域
% end
echoallrd= fftshift(fft(echoall,nm,1),1);%转换到多普勒域
figure;imagesc(x,fd,abs(echoallrd));title('原始回波幅度距离多普勒域');xlim([xylim(1,1),xylim(1,2)]);
% figure;
% imagesc(abs(echoallrd));
%% Chirp Scaling 操作
Dfd=sqrt(1-(c*fd1).^2/(4*vr^2*fc^2));%徙动参数7.17
Km=Kr./(1-(Kr*c*Rc*fd1.^2)./(2*vr^2*fc^3*Dfd.^3));%改变后的距离向调频率7.18
Dfdc=sqrt(1-(c*fdc).^2/(4*vr^2*fc^2));%7.18
% tao_pie=tk_mat-2*Rc./(c*Dfd);%新的距离时间7.27
srcm=exp(im*pi*Km.*(Dfdc./Dfd-1).*(tk_mat-2*Rc./(c*Dfd)).^2);%变标方程7.30
echorcm1rd=echoallrd.*srcm;%变标后
%% 距离向傅里叶变换
echorcm1ff=zeros(nm,nk);
echorcm1ff = fftshift(fft(echorcm1rd,[],2),2);%二维频域
%% 距离压缩，SRC，一致RCMC
f_tao=ones(nm,1)*linspace(-fs/2,fs/2,nk);
echorcm2ff=zeros(nm,nk);
echorcm2ff=echorcm1ff.*exp(im*pi.*Dfd./(Dfdc.*Km).*ft_mat.^2).*...
exp(im*4*pi/c*(1./Dfd-1./Dfdc)*Rc.*ft_mat);%7.32
%% 距离向傅里叶逆变换
echorcm2=zeros(nm,nk);
echorcm2=ifft(fftshift(echorcm2ff,2),[],2);%多普勒域
%% 方位向匹配滤波以及附加相位校正
echorcmc1=zeros(nm,nk);
echorcmc1= echorcm2.*exp(im*4*pi*Rc*fc*Dfd/c).*...
exp(im*4*pi*Km/c^2.*(1-Dfd./Dfdc).*(Rc./...
Dfd-Rc./Dfdc).^2);%方位压缩以及附加相位补偿7.32
%% SAR成像结果
res_image=zeros(nm,nk);
% for i = 1:nk
% res_image(:,i) = ifft(echorcmc1(:,i)); %变回SAR图像域
% end
res_image = ifft(echorcmc1,[],1); %变回SAR图像域
figure;[R,Y] = meshgrid(r,y);
mesh(R,Y,abs(( res_image)));view(0,90);xlim([xylim(1,1),xylim(1,2)]);
ylabel('方位向距离/m');xlabel('距离向距离/m'); ylim([xylim(2,1),xylim(2,2)]);title('成像结果')
参考资料《合成孔径雷达成像算法与实现》——lan G. Cumming, et al.