clc
clear
close all

参数设置
距离向参数设置
R_eta_c = 20e3;             % 景中心斜距
Tr = 2.5e-6;                % 发射脉冲时宽
Kr = 20e12;                 % 距离向调频率
alpha_os_r = 1.2;           % 距离过采样率
Nrg = 320;                  % 距离线采样数
距离向参数计算
Bw = abs(Kr)*Tr;            % 距离信号带宽
Fr = alpha_os_r*Bw;         % 距离向采样率
Nr = round(Fr*Tr);          % 距离采样点数（脉冲序列长度）
方位向参数设置
c = 3e8;                    % 光速
Vr = 150;                   % 等效雷达速度
Vs = Vr;                    % 卫星平台速度
Vg = Vr;                    % 波束扫描速度
f0 = 5.3e9;                 % 雷达工作频率
Delta_f_dop = 80;           % 多普勒带宽
alpha_os_a = 1.25;          % 方位过采样率
Naz = 256;                  % 距离线数
theta_r_c = 1.6;            % 波束斜视角
方位向参数计算
lambda = c/f0;              % 雷达工作波长
eta_c = -R_eta_c*sind(theta_r_c)/Vr;
                            % 波束中心偏移时间
f_eta_c = 2*Vr*sind(theta_r_c)/lambda;
                            % 多普勒中心频率
La = 0.886*2*Vs*cosd(theta_r_c)/Delta_f_dop;
                            % 实际天线长度
Fa = alpha_os_a*Delta_f_dop;% 方位向采样率
Ta = 0.886*lambda*R_eta_c/(La*Vg*cosd(theta_r_c));
                            % 目标照射时间
R0 = R_eta_c*cosd(theta_r_c);
                            % 最短斜距
Ka = 2*Vr^2*cosd(theta_r_c)^3/(lambda*R0);
                            % 方位向调频率
theta_bw = 0.886*lambda/La; % 方位向3dB波束宽度
theta_syn = Vs/Vg*theta_bw; % 合成角
Ls = R_eta_c*theta_syn;     % 合成孔径
其他参数计算
rho_r = c/2/Bw;             % 距离向分辨率 
rho_a = La/2;               % 方位向分辨率
Trg = Nrg/Fr;               % 发射脉冲宽度
Taz = Naz/Fa;               % 目标照射时间
d_t_tau = 1/Fr;             % 距离向采样时间间隔
d_t_eta = 1/Fa;             % 方位向采样时间间隔
d_f_tau = Fr/Nrg;           % 距离向采样频率间隔
d_f_eta = Fa/Naz;           % 方位向采样频率间隔

目标设置
设置目标点距离景中心的距离
A_r =   0;A_a =   0;
B_r = -50;B_a = +50;
C_r = +50;C_a = B_a+(C_r-B_r)*tand(theta_r_c);
D_r = -50,D_a = -50;
坐标
A_x = R0+A_r;A_y = A_a;
B_x = R0+B_r;B_y = B_a;
C_x = R0+C_r;C_y = C_a;
D_x = R0+D_r;D_y = D_a;
N_position = [A_x,A_y;B_x,B_y;C_x,C_y;D_x,D_y];
波束中心穿越时刻
N_target = 4;
Target_eta_c = zeros(1,N_target);
for i = 1:N_target
    Delta_Y = N_position(i,2)-N_position(i,1)*tand(theta_r_c);
    Target_eta_c(i) = Delta_Y/Vs;
end
绝对零多普勒时刻
Target_eta_0 = zeros(1,N_target);
for i = 1:N_target
    Target_eta_0(i) = N_position(i,2)/Vs; 
end

变量设置
时间变量：以景中心绝对零多普勒时刻作为方位向零点
t_tau = (-Trg/2:d_t_tau:Trg/2-d_t_tau)+2*R_eta_c/c;     % 距离时间变量
t_eta = (-Taz/2:d_t_eta:Taz/2-d_t_eta)+eta_c;           % 方位时间变量
r_tau = (t_tau*c/2)*cosd(theta_r_c);                    % 最近距离变量
频率变量
f_tau = fftshift(-Fr/2:d_f_tau:Fr/2-d_f_tau);           % 距离频率变量
f_tau = f_tau-round((f_tau-0)/Fr)*Fr;                   % 将频率折叠入(-Fr/2,Fr/2)，距离可观测频率变量
f_eta = fftshift(-Fa/2:d_f_eta:Fa/2-d_f_eta);           % 方位频率变量
f_eta = f_eta-round((f_eta-f_eta_c)/Fa)*Fa;             % 将频率折叠入f_eta_c附近(-Fa/2,Fa/2)范围，方位可观测频率变量
坐标设置
[t_tauX,t_etaY] = meshgrid(t_tau,t_eta);                % 距离时间X轴，方位时间Y轴
[f_tauX,f_etaY] = meshgrid(f_tau,f_eta);                % 距离频域X轴，方位频域Y轴
[r_tauX,f_eta_Y] = meshgrid(r_tau,f_eta);               % 距离长度X轴，方位频域Y轴

信号设置，原始回波生成
tic                                                     % 计时，与toc搭配使用
wait_title = waitbar(0,'开始生成回波数据 ...'); 
pause(1);
st_tt = zeros(Naz,Nrg);
for i = 1:1
    R_eta = sqrt(N_position(i,1)^2+Vs^2*(t_etaY-Target_eta_0(i)).^2);
                                                        % 瞬时斜距，还有近似公式可以尝试
    A0 = [1,1,1,1]*exp(+1j*0);                          % 后向散射系数
    wr = (abs(t_tauX-2*R_eta/c)<=Tr/2);                 % 距离向包络
    wa = sinc(0.886*atan(Vs*(t_etaY-Target_eta_c(i))/N_position(i,1))/theta_bw).^2;
                                                        % 方位向包络，用波束穿越时刻
%     wa = sinc(0.886*(atan(Vs*(t_etaY-Target_eta_0(i))/N_position(i,1))+theta_r_c)/theta_bw).^2;
    st_tt_target = A0(i)*wr.*wa.*exp(-1j*4*pi*f0*R_eta/c)...
                               .*exp(1j*pi*Kr*(t_tauX-2*R_eta/c).^2);
    st_tt = st_tt+st_tt_target;
    pause(0.001);
    time = toc;
    Display_Data = num2str(roundn(i/N_target*100,-1));
    Display_Str  = ['Computation Progress',Display_Data,'%',' --- ',...
                    'Using Time: ',num2str(time)];
    waitbar(i/N_target,wait_title,Display_Str);         % 三参数：进度，句柄，展示的话
end
pause(1);
close(wait_title);
toc
% figure('Name','原始数据回波'),subplot(221)
% imagesc(real(st_tt))
% xlabel('距离时间(采样点)'),ylabel('方位时间(采样点)'),title('(a)实部')
% subplot(222)
% imagesc(imag(st_tt))
% xlabel('距离时间(采样点)'),ylabel('方位时间(采样点)'),title('(b)虚部')
% subplot(223)
% imagesc(abs(st_tt))
% xlabel('距离时间(采样点)'),ylabel('方位时间(采样点)'),title('(c)幅度')
% subplot(224)
% imagesc(angle(st_tt))
% xlabel('距离时间(采样点)'),ylabel('方位时间(采样点)'),title('(d)相位')

一次距离压缩
方式三：根据脉冲频谱特性直接在频域生成频域匹配滤波器
window = kaiser(Nrg,2.5)';              % 时域窗
Window = fftshift(window);              % 频域窗
% figure,plot(window)
% figure,plot(Window)
Hrf = (abs(f_tauX)<=Bw/2).*Window.*exp(1j*pi*f_tauX.^2/Kr);
Sf_ft = fft(st_tt,Nrg,2);
Srf_tf = Sf_ft.*Hrf;
srt_tt = ifft(Srf_tf,Nrg,2);
% figure('Name','一次距离压缩'),subplot(121)
% imagesc(real(srt_tt))
% xlabel('距离时间(采样点)'),ylabel('方位时间(采样点)'),title('(a)实部')
% subplot(122)
% imagesc(abs(srt_tt))
% xlabel('距离时间(采样点)'),ylabel('方位时间(采样点)'),title('(b)虚部')

方位向FFT
Saf_tf = fft(srt_tt,Naz,1);
% figure('Name','方位FFT'),subplot(121)
% imagesc(real(Saf_tf)),set(gca,'YDir','normal')
% xlabel('距离时间(采样点)'),ylabel('方位时间(采样点)'),title('(a)实部')
% subplot(122)
% imagesc(abs(Saf_tf)),set(gca,'YDir','normal')
% xlabel('距离时间(采样点)'),ylabel('方位时间(采样点)'),title('(b)幅度')

距离徙动校正——8点插值
RCM = lambda^2*r_tauX.*f_etaY.^2/(8*Vr^2);
RCM = R0+RCM-R_eta_c;                       % 将距离徙动量转换到原图坐标系下
offset = RCM/rho_r;                         % 将距离徙动量转换为距离单元偏移量
计算插值表
x_tmp = repmat(-4:3,[16,1]);                % 插值长度
x_tmp = x_tmp+repmat(((1:16)/16).',[1,8]);   % 量化位移
% figure,imagesc(repmat(((1:16)/16)',[1,8])),colorbar
% figure,imagesc(repmat(-4:3,[16,1])),colorbar
% figure,imagesc(repmat(((1:16)/16)',[1,8])+repmat(-4:3,[16,1])),colorbar
hx = sinc(x_tmp);                           % 生成插值核
% % figure,imagesc(hx)
hx = kaiser(8,2.5)'.*hx;
hx = hx./sum(hx,2);                         % 归一化
插值表校正
tic
wait_title = waitbar(0,'开始进行距离徙动校正');
pause(1)
Srcmf_tf_8 = zeros(Naz,Nrg);
for a_tmp = 1:Naz
    for r_tmp = 1:Nrg
        offset_ceil = ceil(offset(a_tmp,r_tmp));
        offset_frac = round((offset_ceil-offset(a_tmp,r_tmp))*16);
        if offset_frac == 0
            Srcmf_tf_8(a_tmp,r_tmp) = Saf_tf(a_tmp,ceil(mod(r_tmp+offset_ceil-0.1,Nrg)));
        else
            Srcmf_tf_8(a_tmp,r_tmp) = Saf_tf(a_tmp,ceil(mod((r_tmp+offset_ceil-4:r_tmp+offset_ceil+3)-0.1,Nrg)))*hx(offset_frac,:).';
        end
    end
    pause(0.001)
    time = toc;
    Display_Data = num2str(roundn(a_tmp/Naz*100,-1));
    Display_Str  = ['Computation Progress ',Display_Data,'%',' --- ',...
                    'Using Time: ',num2str(time)];
    waitbar(a_tmp/Naz,wait_title,Display_Str)
end
pause(1)
close(wait_title)
toc
% figure('Name','8点距离徙动校正'),subplot(121)
% imagesc(real(Srcmf_tf_8)),set(gca,'YDir','normal')
% xlabel('距离时间(采样点)'),ylabel('方位时间(采样点)'),title('(a)实部')
% subplot(122)
% imagesc(abs(Srcmf_tf_8)),set(gca,'YDir','normal')
% xlabel('距离时间(采样点)'),ylabel('方位时间(采样点)'),title('(b)幅度')

距离徙动校正——最近邻域
x_tmp = repmat(-4:3,[16,1]);
x_tmp = x_tmp+repmat(((1:16)/16).',[1,8]);
hx = sinc(x_tmp);
hx = kaiser(8,2.5)'.*hx;
hx = hx./sum(hx,2);
插值表校正
tic
wait_title = waitbar(0,'开始进行距离徙动校正(最近邻域插值)');
pause(1)
Srcmf_tf_0 = zeros(Naz,Nrg);
for a_tmp = 1:Naz
    for r_tmp = 1:Nrg
        offset_ceil = ceil(offset(a_tmp,r_tmp));
        offset_frac = round(offset_ceil-offset(a_tmp,r_tmp));
        if offset_frac == 0
            point = r_tmp+offset_ceil;
        else
            point = r_tmp+offset_ceil-1;
        end
        if point < 1
            point = 1;
        elseif point > Nrg
            point = Nrg;
        end
        Srcmf_tf_0(a_tmp,r_tmp) = Saf_tf(a_tmp,point);
    end
    pause(0.001)
    time = toc;
    Display_Data = num2str(roundn(a_tmp/Naz*100,-1));
    Display_Str  = ['Computation Progress ',Display_Data,'%',' --- ',...
                    'Using Time: ',num2str(time)];
    waitbar(a_tmp/Naz,wait_title,Display_Str)
end
pause(1)
close(wait_title)
toc

距离徙动校正——4点插值
计算插值表
x_tmp = repmat(-2:1,[16,1]);                % 插值长度
x_tmp = x_tmp+repmat(((1:16)/16).',[1,4]);  % 量化位移
hx = sinc(x_tmp);                           % 生成插值核
hx = kaiser(4,2.5)'.*hx;
hx = hx./sum(hx,2);                         % 归一化
插值表校正
tic
wait_title = waitbar(0,'开始进行距离徙动校正(4点)');
pause(1)
Srcmf_tf_4 = zeros(Naz,Nrg);
for a_tmp = 1:Naz
    for r_tmp = 1:Nrg
        offset_ceil = ceil(offset(a_tmp,r_tmp));
        offset_frac = round((offset_ceil-offset(a_tmp,r_tmp))*16);
        if offset_frac == 0
            Srcmf_tf_4(a_tmp,r_tmp) = Saf_tf(a_tmp,ceil(mod(r_tmp+offset_ceil-0.1,Nrg)));
        else
            Srcmf_tf_4(a_tmp,r_tmp) = Saf_tf(a_tmp,ceil(mod((r_tmp+offset_ceil-2:r_tmp+offset_ceil+1)-0.1,Nrg)))*hx(offset_frac,:).';
        end
    end
    pause(0.001)
    time = toc;
    Display_Data = num2str(roundn(a_tmp/Naz*100,-1));
    Display_Str  = ['Computation Progress ',Display_Data,'%',' --- ',...
                    'Using Time: ',num2str(time)];
    waitbar(a_tmp/Naz,wait_title,Display_Str)
end
pause(1)
close(wait_title)
toc

方位压缩
Ka = 2*Vr^2*cosd(theta_r_c)^3./(lambda*r_tauX);
Haf = exp(-1j*pi*f_etaY.^2./Ka);                    % 匹配滤波器
Haf_offset = exp(-1j*2*pi*f_etaY*eta_c);            % 时间补偿项
Soutf_tf_8 = Srcmf_tf_8.*Haf.*Haf_offset;
soutt_tt_8 = ifft(Soutf_tf_8,Naz,1);
Soutf_tf_0 = Srcmf_tf_0.*Haf.*Haf_offset;
soutt_tt_0 = ifft(Soutf_tf_0,Naz,1);
Soutf_tf_4 = Srcmf_tf_4.*Haf.*Haf_offset;
soutt_tt_4 = ifft(Soutf_tf_4,Naz,1);
% figure('Name','方位压缩8'),subplot(121)
% imagesc(real(soutt_tt_8))
% subplot(122)
% imagesc(abs(soutt_tt_8))
% figure('Name','方位压缩0'),subplot(121)
% imagesc(real(soutt_tt_0))
% subplot(122)
% imagesc(abs(soutt_tt_0))
% figure('Name','方位压缩4'),subplot(121)
% imagesc(real(soutt_tt_4))
% subplot(122)
% imagesc(abs(soutt_tt_4))

点目标分析
srcmt_tt_8 = ifft(Srcmf_tf_8,Naz,1);
Arcm_8 = srcmt_tt_8(:,round(Nrg/2+1+2*(N_position(1,1)-R0)/c*Fr)); % Nrg/2+1为距离维度的中心点，剩下的计算目标点A相对于最短斜距R0的偏移量
Arcm_8 = abs(Arcm_8)/max(abs(Arcm_8));

srcmt_tt_0 = ifft(Srcmf_tf_0,Naz,1);
Arcm_0 = srcmt_tt_0(:,round(Nrg/2+1+2*(N_position(1,1)-R0)/c*Fr)); % Nrg/2+1为距离维度的中心点，剩下的计算目标点A相对于最短斜距R0的偏移量
Arcm_0 = abs(Arcm_0)/max(abs(Arcm_0));

srcmt_tt_4 = ifft(Srcmf_tf_4,Naz,1);
Arcm_4 = srcmt_tt_4(:,round(Nrg/2+1+2*(N_position(1,1)-R0)/c*Fr)); % Nrg/2+1为距离维度的中心点，剩下的计算目标点A相对于最短斜距R0的偏移量
Arcm_4 = abs(Arcm_4)/max(abs(Arcm_4));
方位切片
len_az = 16;
cut_az = -len_az/2:len_az/2-1;
out_az_8 = soutt_tt_8(round(Naz/2+1+N_position(1,2)/Vr*Fa)+cut_az, ...
                      round(Nrg/2+1+2*(N_position(1,1)-R0)/c*Fr)).';
spec_az_8 = fft(out_az_8);
spec_up_az_8 = ifft(spec_az_8,16*len_az);                           % 频域补零
% spec_az_8 = fft(out_az_8,16*len_az);                              % 时域补零，注意区别
% spec_up_az_8 = ifft(spec_az_8);
spec_up_az_8 = db(abs(spec_up_az_8)/max(abs(spec_up_az_8)));

out_az_0 = soutt_tt_0(round(Naz/2+1+N_position(1,2)/Vr*Fa)+cut_az, ...
                      round(Nrg/2+1+2*(N_position(1,1)-R0)/c*Fr)).';
spec_az_0 = fft(out_az_0);
spec_up_az_0 = ifft(spec_az_0,16*len_az);
spec_up_az_0 = db(abs(spec_up_az_0)/max(abs(spec_up_az_0)));

out_az_4 = soutt_tt_4(round(Naz/2+1+N_position(1,2)/Vr*Fa)+cut_az, ...
                      round(Nrg/2+1+2*(N_position(1,1)-R0)/c*Fr)).';
spec_az_4 = fft(out_az_4);
spec_up_az_4 = ifft(spec_az_4,16*len_az);
spec_up_az_4 = db(abs(spec_up_az_4)/max(abs(spec_up_az_4)));

绘图
H = figure();
set(H,'position',[100,100,900,600]);
subplot(231),plot(abs(wa(:,1))),axis([0 Naz,0 1])
xlabel('方位向(采样点)'),ylabel('幅度'),title('(a)理想的距离徙动校正(8点插值)')
subplot(232),plot(abs(Arcm_0)),axis([0 Naz,0 1])
xlabel('方位向(采样点)'),ylabel('幅度'),title('(b)最近邻域')
subplot(233),plot(abs(Arcm_4)),axis([0 Naz,0 1])
xlabel('方位向(采样点)'),ylabel('幅度'),title('(c)4点插值')
subplot(234),plot(spec_up_az_8),axis([0 Naz,-30 0])
xlabel('方位向(采样点)'),ylabel('幅度(dB)')
subplot(235),plot(spec_up_az_0),axis([0 Naz,-30 0])
xlabel('方位向(采样点)'),ylabel('幅度(dB)')
subplot(236),plot(spec_up_az_4),axis([0 Naz,-30 0])
xlabel('方位向(采样点)'),ylabel('幅度(dB)')