基本信息
源码名称：用于SAR成像的RD算法
源码大小：8.15KB
文件格式：.zip
开发语言：MATLAB
更新时间：2020-10-04
友情提示：（无需注册或充值，赞助后即可获取资源下载链接）
嘿，亲！知识可是无价之宝呢，但咱这精心整理的资料也耗费了不少心血呀。小小地破费一下，绝对物超所值哦！如有下载和支付问题，请联系我们QQ(微信同号)：813200300
本次赞助数额为： 2 元×
微信扫码支付：2 元
×
请留下您的邮箱，我们将在2小时内将文件发到您的邮箱
源码介绍
用于SAR成像的RD算法
clear all; close all;
C = 2.9979e8;
Fc = 9.6e9;
lambda = C/Fc;
Rfactor = 2; % 距离维采样倍数因子 1.5
Afactor = 2; % 方位向采样倍数因子 1.2
V = 150;
H = 6000;
fuyang = 45/180*pi;
Yc = H * tan(fuyang);
R0 = sqrt(H^2 Yc^2); %R0为景中心距离飞机飞行轨迹的最短距离（即为垂线长度）
Yw = 100;
Theta_rc = 45/180*pi; %斜视角
disp(['斜视角 = ',num2str(Theta_rc/pi*180),'度']);
Tr = 0.0000025;
Br = 100e6;
Kr = Br/Tr;
disp(['信号带宽 = ',num2str(Br/1e6),'MHz']);
juli_fenbian = 0.886*C/2/Br;
disp(['距离向分辨率 = ',num2str(juli_fenbian)]);
La=juli_fenbian*2;
%La = 10;
beita = 0.886*lambda/La;
fangwei_fenbian = La/2;
disp(['方位向分辨率 = ',num2str(fangwei_fenbian)]);
Fdop = 0.886*2*V*cos(Theta_rc)/La; %多普勒带宽
disp(['多普勒带宽 = ',num2str(Fdop),'Hz']);
B_truth = Fdop 2*V*Br/C*sin(Theta_rc)*cos(beita/2);
PRF = Afactor*B_truth;
disp(['PRF = ',num2str(PRF)]);
PRT = 1/PRF;
ds = PRT;
Xmin = -50;
Xmax = 50;
Lsar = R0/cos(Theta_rc)*beita/cos(Theta_rc);
disp(['合成孔径长度 = ',num2str(Lsar),'m']);
Tsar = Lsar/V;
disp(['目标照射时间 = ',num2str(Tsar),'s']);
alpha = beita/2;
Ymin = Yc-Yw;
Ymax = Yc Yw;
Rmin = sqrt(H^2 Ymin^2);
Rmax = sqrt(H^2 Ymax^2);
Rw = Yw *sin(fuyang);
Rw2 = Tr*C/2;
R_min = Rmin/cos(Theta_rc-alpha) - 12*Rw2; %1.25Ghz 12 9.6Ghz 2
R_max = Rmax/cos(Theta_rc alpha) 17*Rw2; %1.25Ghz 17 9.6Ghz 3
Fsr = Rfactor * Kr * Tr;
dt = 1/Fsr;
tm = 2*R_min/C:dt:(2*R_max/C-dt);
M = length(tm);
if mod(M,2)==1 %判断M是否是偶数
disp('M是奇数，M=M 1变为偶数')
M=M 1;
tm = 2*R_min/C:dt:2*R_max/C;
end
S_qian = 1*(Rmax*tan(alpha Theta_rc) - Rmax*tan(Theta_rc)); %实际上是一个近似
S_hou = 1*(Rmin*tan(alpha-Theta_rc) Rmin*tan(Theta_rc)); %实际上是一个近似
% S_qian = 1*(R0*tan(alpha Theta_rc) - R0*tan(Theta_rc)); %实际上是一个近似
% S_hou = 1*(R0*tan(alpha-Theta_rc) R0*tan(Theta_rc)); %实际上是一个近似
% S_qian = S_hou;
disp(['S_qian = ',num2str(S_qian)]);
disp(['S_hou = ',num2str(S_hou)]);
eta_c = -R0*tan(Theta_rc)/V;
sn = (Xmin-S_qian-Rmax*tan(Theta_rc))/V:PRT:((Xmax S_hou-Rmin*tan(Theta_rc))/V-PRT);
%sn = ((Xmin-S_qian-R0*tan(Theta_rc))/V:PRT:((Xmax S_hou-R0*tan(Theta_rc))/V-PRT)) (Xmin-S_qian)/V;
N = length(sn);
if mod(N,2)==1 %判断N是否是偶数
disp('N是奇数，N=N 1变为偶数')
N=N 1;
sn = ((Xmin-S_qian-Rmax*tan(Theta_rc))/V:PRT:(Xmax S_hou-Rmin*tan(Theta_rc))/V);
%sn = ((Xmin-S_qian-R0*tan(Theta_rc))/V:PRT:(Xmax S_hou-R0*tan(Theta_rc))/V) (Xmin-S_qian)/V;
end
Ntarget = 1;
Ptarget = [(Xmin Xmax)/2,Yc, 2
Xmin,Ymin,1
(Xmin Xmax)/2,Ymin,1
Xmax,Ymin,1
Xmin,Ymax,1
(Xmin Xmax)/2,Ymax,1
Xmax,Ymax,1
Xmin,Yc,1
Xmax,Yc,1 ];
%% 每个格点间的距离
juli_gedian = dt*C/2;
disp(['距离向单元格点间距离为 ',num2str(juli_gedian),'m']);
fangwei_gedian = PRT*V;
disp(['方位向单元格点间距离为 ',num2str(fangwei_gedian),'m']);
%% 产生需要处理的数据矩阵
K = Ntarget;
T = Ptarget;
Srnm = zeros(N,M);
h = waitbar(0,'SAR回波生成');
for k=1:1:K
T(k,2) = sqrt(H^2 T(k,2)^2);
S_qian = T(k,2)*tan(alpha Theta_rc) - T(k,2)*tan(Theta_rc);
S_hou = T(k,2)*tan(alpha-Theta_rc) T(k,2)*tan(Theta_rc);
% S_qian = S_hou;
sigma =T(k,3);
Dslow = sn*V T(k,2)*tan(Theta_rc)-T(k,1); %波束中心距离目标的距离
% faiweijuli = T(k,2)*tan(Theta_rc) - Dslow; %sn*V-T(k,1) %雷达距离目标的水平距离
% Ra = sqrt(faiweijuli.^2 T(k,2)^2);
Ra = sqrt((sn*V-T(k,1)).^2 T(k,2)^2); %雷达距离目标的距离
tau = 2*Ra/C;
Dfast = ones(N,1)*tm - tau'*ones(1,M);
clear tau
% myphase = pi*Kr*(Dfast-Tr/2).^2- (4*pi/lambda)*(Ra'*ones(1,M));
Srnm = Srnm sigma * exp(1i*(pi*Kr*(Dfast-Tr/2).^2- (4*pi/lambda)*(Ra'*ones(1,M))))...
.*(0<Dfast&Dfast<Tr).*(((Dslow>=(-S_qian))&(Dslow<=S_hou)).'*ones(1,M))...
;%.*(sinc(0.886/beita*(atan(tan(Theta_rc)-Dslow/T(k,2))-Theta_rc)).'*ones(1,M)).^2; %最后一项为模拟sinc^2函数
clear myphase Dfast Dslow Ra
waitbar(k/K);
end
close(h)
clear tm
%% 距离向脉冲压缩(方式1)
tr1 = -Tr/2:dt:Tr/2-dt;
Refr = zeros(1,M);
Refr(1:length(tr1)) = exp(1i*pi*Kr*tr1.^2);
%%% Refr(1:length(tr1)) = kaiser(length(tr1),4).'.*exp(1i*pi*Kr*tr1.^2); %% 加窗
Refr =ones(N,1)* Refr;
Srnm = fftshift(ifty(fty(Srnm).*(conj(fty(Refr)))),2);
f_tao = [-M/2:1:M/2-1]/M*Fsr;
% Refr = exp(1j*pi/Kr*ones(N,1)*(f_tao.^2));
%%% Refr = (ones(N,1)*kaiser(M,5).').*exp(1j*pi/Kr*ones(N,1)*(f_tao.^2)); %加窗
%%% Refr = (ones(N,1)*hamming(M).').*exp(1j*pi/Kr*ones(N,1)*(f_tao.^2)); %加窗
% Srnm = ifty(fty(Srnm).*Refr);
% R_min = R_min - Tr/2*C/2;
clear Refr tr1
% figure
% imagesc(abs(Sr)); axis xy;
%% 距离多普勒域
Sr_RD = ftx(Srnm);
% figure
% imagesc(abs(Sr_RD));axis xy;
% title('距离徙动校正前，无f_c');
%% 多普勒频率估计
f_c = 2*V*sin(Theta_rc)/lambda;
D = sqrt(1 - lambda^2/(4*V^2)*(f_c (-PRF/2:PRF/N:PRF/2-PRF/N)').^2);
Sr_RD = ftx(Srnm.*exp(-1j*2*pi*f_c*((sn).')*ones(1,M)));
mohushu = fix(f_c/PRF);
disp(['模糊数 = ',num2str(mohushu)]);
% tiaozheng = - round(mod(f_c,PRF)/PRF*size(Sr_RD,1)); %反向调整的数
% Sr_RD = circshift(Sr_RD,tiaozheng,1);
% figure
% imagesc(abs(Sr_RD));axis xy;
% title('距离徙动校正前，有f_c');
clear Srnm sn
%clear sn
Sr_f2 = fty(Sr_RD);
H_src2 = exp(1j*4*pi*R0*Fc/C*1/(2*1)*(((D.^2-1)./D.^3)*f_tao.^2)/Fc^2);
Sr_f2 = Sr_f2.*H_src2;
clear H_src2
%% 三次相位补偿
H_xiangwei3 =exp(1j*4*pi*R0*Fc/C*(-1)*(1./(2*(D.^3)) - 1./(2*(D.^5)))*f_tao.^3/Fc^3);
Sr_f2 = Sr_f2.*H_xiangwei3;
clear H_xiangwei3
%% 四次相位补偿
H_xiangwei4 = exp(1j*4*pi*R0*Fc/C*...
(-1)*(1./(8*(D.^3)) - 3./(4*(D.^5)) 5./(8*(D.^7)))*f_tao.^4/Fc^4);
Sr_f2 = Sr_f2.*H_xiangwei4;
clear H_xiangwei4
%% 五次相位补偿
H_xiangwei5 = exp(1j*4*pi*R0*Fc/C*...
(3./(8*(D.^5)) - 5./(4*(D.^7)) 7./(8*(D.^9)))*f_tao.^5/Fc^5);
Sr_f2 = Sr_f2.*H_xiangwei5;
clear H_xiangwei5
%% 六次相位补偿
H_xiangwei6 = exp(1j*4*pi*R0*Fc/C*...
(1./(16*(D.^5)) - 15./(16*(D.^7)) 35./(16*(D.^9)) - 21./(16*(D.^11)))...
*f_tao.^6/Fc^6);
Sr_f2 = Sr_f2.*H_xiangwei6;
clear H_xiangwei6
%% 七次相位补偿
H_xiangwei7 = exp(1j*4*pi*R0*Fc/C*...
(-1)*(5./(16*(D.^7)) - 35./(16*(D.^9)) 63./(16*(D.^11)) - 33./(16*(D.^13)))...
*f_tao.^7/Fc^7);
Sr_f2 = Sr_f2.*H_xiangwei7;
clear H_xiangwei7
%% 八次相位补偿
H_xiangwei8 = exp(1j*4*pi*R0*Fc/C*...
(-1)*(5./(128*(D.^7)) - 35./(32*(D.^9)) 315./(64*(D.^11)) - 231./(32*(D.^13)) 429./(128*(D.^15)))...
*f_tao.^8/Fc^8);
Sr_f2 = Sr_f2.*H_xiangwei8;
clear H_xiangwei8
%% 九次相位补偿
H_xiangwei9 = exp(1j*4*pi*R0*Fc/C*...
(35./(128*(D.^9)) - 105./(32*(D.^11)) 693./(64*(D.^13)) - 429./(32*(D.^15)) 715./(128*(D.^17)))...
*f_tao.^9/Fc^9);
Sr_f2 = Sr_f2.*H_xiangwei9;
clear H_xiangwei9
%% 十次相位补偿
H_xiangwei10 = exp(1j*4*pi*R0*Fc/C*...
(7./(256*(D.^9)) - 315./(256*(D.^11)) 1155./(128*(D.^13)) - 3003./(128*(D.^15)) 6435./(256*(D.^17)) - 2431./(256*(D.^19)))...
*f_tao.^10/Fc^10);
Sr_f2 = Sr_f2.*H_xiangwei10;
clear H_xiangwei10
%% 十一次相位补偿
H_xiangwei11 = exp(1j*4*pi*R0*Fc/C*...
(-1)*(63./(256*(D.^11)) - 1155./(256*(D.^13)) 3003./(128*(D.^15)) - 6435./(128*(D.^17)) 12155./(256*(D.^19)) - 4199./(256*(D.^21)))...
*f_tao.^11/Fc^11);
Sr_f2 = Sr_f2.*H_xiangwei11;
clear H_xiangwei11
%% 十二次相位补偿
H_xiangwei12 = exp(1j*4*pi*R0*Fc/C*...
(-1)*(21./(1024*(D.^11)) - 693./(512*(D.^13)) 15015./(1024*(D.^15)) - 15015./(256*(D.^17)) 109395./(1024*(D.^19)) - 46189./(512*(D.^21)) 29393./(1024*(D.^23)))...
*f_tao.^12/Fc^12);
Sr_f2 = Sr_f2.*H_xiangwei12;
clear H_xiangwei12
%% 十三次相位补偿
H_xiangwei13 = exp(1j*4*pi*R0*Fc/C*...
(231./(1024*(D.^13)) - 3003./(512*(D.^15)) 45045./(1024*(D.^17)) - 36465./(256*(D.^19)) 230945./(1024*(D.^21)) - 88179./(512*(D.^23)) 52003./(1024*(D.^25)))...
*f_tao.^13/Fc^13);
Sr_f2 = Sr_f2.*H_xiangwei13;
clear H_xiangwei13
%% 十四次相位补偿
H_xiangwei14 = exp(1j*4*pi*R0*Fc/C*...
(33./(2048*(D.^13)) - 3003./(2048*(D.^15)) 45045./(2048*(D.^17)) - 255255./(2048*(D.^19)) 692835./(2048*(D.^21)) - 969969./(2048*(D.^23)) 676039./(2048*(D.^25)) - 185725./(2048*(D.^27)))...
*f_tao.^14/Fc^14);
Sr_f2 = Sr_f2.*H_xiangwei14;
clear H_xiangwei14
%Sr_f2_he =Sr_f2_he Sr_f2;
%end
%Sr_RD = ifty(Sr_f2_he);
Sr_RD = ifty(Sr_f2);
clear Sr_f2 Sr_f2_he;
% figure
% imagesc(abs(Sr_RD));axis xy;
% title('SRC后，有f_c');
%
% %%
tic
%Rmin = R_min 2*Rw;
RMC = 2*Fsr/C*((1./D)*(Rmin-Rw2 (0:1:M-1)*C/2/Fsr)-R_min);
%% 单独的两个循环比一个循环的时间还要快，为什么？
for n=1:N
%RMCmaxtix(n,:) = interp1(Sr_RD(n,:),1 RMC(n,:),'spline');
Sr_RD(n,:) = interp1(Sr_RD(n,:),1 RMC(n,:),'spline');%加1的原因是因为矩阵序列是从1开始的
end
toc
tic
%% 一种优化方式
for n=1:N
for m=1:M
if RMC(n,m)<0 %%数据溢出时的处理方式
Sr_RD(n,m) = 0;
elseif (1 RMC(n,m))>M
Sr_RD(n,m:M) = zeros(1,M-m 1);
break
end
end
end
% for n=1:N
% for m=1:M
% if RMC(n,m)<0 || (1 RMC(n,m))>M %%数据溢出时的处理方式
% RMCmaxtix(n,m) = 0;
% end
% end
% end
clear RMC
toc
%
% tic
% RMCmaxtix = zeros(N,M);
% h = waitbar(0,'Sinc插值');
% P = 8; %8点sinc插值.
% P_start = P/2;
% P_chazhi = -P/2 1:1:P/2;
% for n=1:N
% for m=P_start:M
% %delta_R =(Rmin-2*RwC (m-1)*/2/Fsr)*(1/D(n)-1);
% %delta_R =(Rmin-1*Rw (m-1)*C/2/Fsr)*(1/D(n)-1); %这一步必须慎重，R(d)=R(0)/D,注意一一对应.
% %zhengduan_num = R_min - Rmin 2*Rw - (m-1)*C/2/Fsr;
% %delta_R_new = (Rmin-2*Rw (m-1)*C/2/Fsr)*(1/D(n)-1) -(R_min - Rmin 2*Rw - (m-1)*C/2/Fsr);
% %RMC = 2*Fsr/C*((Rmin-2*Rw (m-1)*C/2/Fsr)*(1/D(n)-1) -(R_min - Rmin 2*Rw - (m-1)*C/2/Fsr));
% RMC = 2*Fsr/C*((Rmin-2*Rw (m-1)*C/2/Fsr)*(1/D(n)) - R_min);
% if RMC>=4 && (floor(RMC) P/2)<=M
% %delta_R =(1/8)*(lambda/V)^2*(R_min (m)*C/2/Fsr)*((n-N/2)*PRF/N f_c)^2;
% %delta_R =(1/8)*(lambda/V)^2*(Rmin (m)*C/2/Fsr)*((n-N/2)*PRF/N)^2;
% delta_RMC = RMC - floor(RMC);
% %biaoji(n,m) = delta_RMC;
% for i=-P/2 1:1:P/2
% RMCmaxtix(n,m) = RMCmaxtix(n,m) Sr_RD(n,floor(RMC) i)*sinc(-i delta_RMC);
% end
% % RMCmaxtix(n,m) = RMCmaxtix(n,m) sum(Sr_RD(n,floor(RMC) P_chazhi).*sinc(-P_chazhi delta_RMC));
% end
% end
% waitbar(n/N);
% end
% close(h)
% clear RMC Sr_RD
% toc
% figure
% imagesc(abs(RMCmaxtix));axis xy;
% title('距离徙动校正后');
% figure
% imagesc(abs(fty(RMCmaxtix))); axis xy;
% Sr_RD = RMCmaxtix;
%% 方位向压缩
Fa=PRF;
fn = f_c (-Fa/2:Fa/N:Fa/2-Fa/N)';
%M=floor(M/6)
tau = (0:M-1)*dt 2*(Rmin - Rw2)/C;
Refa = exp(1j*4*pi/lambda*(D)*tau*C/2); %注意D减1的效果.
myRD = Sr_RD(:,1:M).*Refa;
%myRD = Sr_RD.*Refa;
clear Sr_RD tau Ka Refa D
R0 = sqrt(H^2 Yc^2); %R0为景中心距离飞机飞行轨迹的最短距离（即为垂线长度）
lingpin = exp(1j*2*pi*R0*tan(Theta_rc)/V*(fn*ones(1,M)));
Sa = iftx(myRD.*lingpin); %压到零频率位置，由于混叠，我在频域乘以一个线性相位，在时间域还原到零频位置.
clear lingpin myRD f_tao fn
[a,b]=find(abs(Sa)==max(max(abs(Sa))))
figure
imagesc(abs(Sa)); axis xy;
%% 后续处理
n_1=100;
n_2=100;
zonghe = Sa(a-n_1:1:a n_1-1,b-n_2:1:b n_2-1);
% figure
% imagesc(abs(zonghe));axis xy;
zonghe_fft2 = fftshift(fft2(zonghe));
% figure
% imagesc(abs(zonghe_fft2));axis xy;
zonghe_fft2_yidong = findpu_min(zonghe_fft2);
% zonghe_fft2_yidong = zonghe_fft2;
% figure
% imagesc(abs(zonghe_fft2_yidong));axis xy;
%% 16倍数升采样
Da = zeros(2*n_1*16,2*n_2*16);
Da(2*n_1*8-n_1:1:2*n_1*8 n_1-1,2*n_2*8-n_2:1:2*n_2*8 n_2-1)=zonghe_fft2_yidong;
figure
imagesc(abs(Da))
zuihou_result =ifft2(ifftshift(Da));axis xy;
zuihou =abs(zuihou_result)/max(max(abs(zuihou_result)));
figure
imagesc(abs(zuihou));axis xy;
figure
contour(abs(zuihou),35);
zuihou2 = 20*log10(zuihou);
figure
contour(zuihou2,linspace(-35,0,8));
fanwei = 3;
I = zuihou2(2*n_1*8-fanwei*n_1:1:2*n_1*8 fanwei*n_1-1,2*n_2*8-fanwei*n_2:1:2*n_2*8 fanwei*n_2-1); %% -3*n_1
I2 = zuihou_result(2*n_1*8-fanwei*n_1:1:2*n_1*8 fanwei*n_1-1,2*n_2*8-fanwei*n_2:1:2*n_2*8 fanwei*n_2-1);
I2 = I2/max(max(abs(I2)));
figure,set(gcf,'Color','w')
contour(I,linspace(-35,0,8));
axis square
%% 后续处理
%% 将两个方向的格点拉均匀
%%
[m,n] = size(I)
%% 变换后，两个方向上的格点距离均为fangwei_gedian
%n = round(n*juli_gedian/fangwei_gedian)
%gedian_juli = fangwei_gedian;
%% 变换后，两个方向上的格点距离均为juli_gedian
m = round(m*fangwei_gedian/juli_gedian)
gedian = juli_gedian;
I_new = imresize(I2,[m n],'bicubic');
I_new_log = 20*log10(abs(I_new)/max(max(abs(I_new))));
figure
contour(I_new_log,linspace(-35,0,8));
if m<n
if mod(m,2)==1
m=m-1;
end
I_new2 = I_new(:,round(n/2)-m/2 1:1:round(n/2) m/2);
else
if mod(n,2)==1
n=n-1;
end
I_new2 = I_new(round(m/2)-n/2 1:1:round(m/2) n/2,:);
end
I_new2_log = 20*log10(abs(I_new2)/max(max(abs(I_new2))));
figure
contour(I_new2_log,linspace(-35,0,8));axis square
% figure
% imagesc(abs(fftshift(fft2(I_new2))));
I_new3_real = imrotate(real(I_new2),Theta_rc/pi*180,'bicubic','crop');%,'crop'
I_new3_imag = imrotate(imag(I_new2),Theta_rc/pi*180,'bicubic','crop');%,'crop'
I_new3 = I_new3_real 1j*I_new3_imag;
I_new3_log = 20*log10(abs(I_new3)/max(max(abs(I_new3))));
figure
contour(I_new3_log,linspace(-35,0,8));axis square;
[m,n] = find(I_new3_log ==max(max(I_new3_log)))
% figure
% plot(I_new3_log(m,:)); title('距离向剖面图');
% axis([-inf inf -40 0]);grid on;
% figure
% plot(I_new3_log(:,n)); title('方位向剖面图');
% axis([-inf inf -40 0]);grid on;
%
a_pow_new = I_new3(:,n);
r_pow_new = I_new3(m,:);
% figure
% plot(angle(r_pow_new)*180/pi); title('距离向剖面图');
% figure
% plot(angle(a_pow_new)*180/pi); title('方位向剖面图');
[a_PSIR a_ISIR a_IRW]=PSISIRIRW(a_pow_new,gedian/16,fangwei_fenbian) %除以16是因为进行了16倍升采样
[r_PSIR r_ISIR r_IRW]=PSISIRIRW(r_pow_new,gedian/16,juli_fenbian) %除以16是因为进行了16倍升采样
用于SAR成像的RD算法
clear all; close all;
C = 2.9979e8;
Fc = 9.6e9;
lambda = C/Fc;
Rfactor = 2; % 距离维采样倍数因子 1.5
Afactor = 2; % 方位向采样倍数因子 1.2
V = 150;
H = 6000;
fuyang = 45/180*pi;
Yc = H * tan(fuyang);
R0 = sqrt(H^2 Yc^2); %R0为景中心距离飞机飞行轨迹的最短距离（即为垂线长度）
Yw = 100;
Theta_rc = 45/180*pi; %斜视角
disp(['斜视角 = ',num2str(Theta_rc/pi*180),'度']);
Tr = 0.0000025;
Br = 100e6;
Kr = Br/Tr;
disp(['信号带宽 = ',num2str(Br/1e6),'MHz']);
juli_fenbian = 0.886*C/2/Br;
disp(['距离向分辨率 = ',num2str(juli_fenbian)]);
La=juli_fenbian*2;
%La = 10;
beita = 0.886*lambda/La;
fangwei_fenbian = La/2;
disp(['方位向分辨率 = ',num2str(fangwei_fenbian)]);
Fdop = 0.886*2*V*cos(Theta_rc)/La; %多普勒带宽
disp(['多普勒带宽 = ',num2str(Fdop),'Hz']);
B_truth = Fdop 2*V*Br/C*sin(Theta_rc)*cos(beita/2);
PRF = Afactor*B_truth;
disp(['PRF = ',num2str(PRF)]);
PRT = 1/PRF;
ds = PRT;
Xmin = -50;
Xmax = 50;
Lsar = R0/cos(Theta_rc)*beita/cos(Theta_rc);
disp(['合成孔径长度 = ',num2str(Lsar),'m']);
Tsar = Lsar/V;
disp(['目标照射时间 = ',num2str(Tsar),'s']);
alpha = beita/2;
Ymin = Yc-Yw;
Ymax = Yc Yw;
Rmin = sqrt(H^2 Ymin^2);
Rmax = sqrt(H^2 Ymax^2);
Rw = Yw *sin(fuyang);
Rw2 = Tr*C/2;
R_min = Rmin/cos(Theta_rc-alpha) - 12*Rw2; %1.25Ghz 12 9.6Ghz 2
R_max = Rmax/cos(Theta_rc alpha) 17*Rw2; %1.25Ghz 17 9.6Ghz 3
Fsr = Rfactor * Kr * Tr;
dt = 1/Fsr;
tm = 2*R_min/C:dt:(2*R_max/C-dt);
M = length(tm);
if mod(M,2)==1 %判断M是否是偶数
disp('M是奇数，M=M 1变为偶数')
M=M 1;
tm = 2*R_min/C:dt:2*R_max/C;
end
S_qian = 1*(Rmax*tan(alpha Theta_rc) - Rmax*tan(Theta_rc)); %实际上是一个近似
S_hou = 1*(Rmin*tan(alpha-Theta_rc) Rmin*tan(Theta_rc)); %实际上是一个近似
% S_qian = 1*(R0*tan(alpha Theta_rc) - R0*tan(Theta_rc)); %实际上是一个近似
% S_hou = 1*(R0*tan(alpha-Theta_rc) R0*tan(Theta_rc)); %实际上是一个近似
% S_qian = S_hou;
disp(['S_qian = ',num2str(S_qian)]);
disp(['S_hou = ',num2str(S_hou)]);
eta_c = -R0*tan(Theta_rc)/V;
sn = (Xmin-S_qian-Rmax*tan(Theta_rc))/V:PRT:((Xmax S_hou-Rmin*tan(Theta_rc))/V-PRT);
%sn = ((Xmin-S_qian-R0*tan(Theta_rc))/V:PRT:((Xmax S_hou-R0*tan(Theta_rc))/V-PRT)) (Xmin-S_qian)/V;
N = length(sn);
if mod(N,2)==1 %判断N是否是偶数
disp('N是奇数，N=N 1变为偶数')
N=N 1;
sn = ((Xmin-S_qian-Rmax*tan(Theta_rc))/V:PRT:(Xmax S_hou-Rmin*tan(Theta_rc))/V);
%sn = ((Xmin-S_qian-R0*tan(Theta_rc))/V:PRT:(Xmax S_hou-R0*tan(Theta_rc))/V) (Xmin-S_qian)/V;
end
Ntarget = 1;
Ptarget = [(Xmin Xmax)/2,Yc, 2
Xmin,Ymin,1
(Xmin Xmax)/2,Ymin,1
Xmax,Ymin,1
Xmin,Ymax,1
(Xmin Xmax)/2,Ymax,1
Xmax,Ymax,1
Xmin,Yc,1
Xmax,Yc,1 ];
%% 每个格点间的距离
juli_gedian = dt*C/2;
disp(['距离向单元格点间距离为 ',num2str(juli_gedian),'m']);
fangwei_gedian = PRT*V;
disp(['方位向单元格点间距离为 ',num2str(fangwei_gedian),'m']);
%% 产生需要处理的数据矩阵
K = Ntarget;
T = Ptarget;
Srnm = zeros(N,M);
h = waitbar(0,'SAR回波生成');
for k=1:1:K
T(k,2) = sqrt(H^2 T(k,2)^2);
S_qian = T(k,2)*tan(alpha Theta_rc) - T(k,2)*tan(Theta_rc);
S_hou = T(k,2)*tan(alpha-Theta_rc) T(k,2)*tan(Theta_rc);
% S_qian = S_hou;
sigma =T(k,3);
Dslow = sn*V T(k,2)*tan(Theta_rc)-T(k,1); %波束中心距离目标的距离
% faiweijuli = T(k,2)*tan(Theta_rc) - Dslow; %sn*V-T(k,1) %雷达距离目标的水平距离
% Ra = sqrt(faiweijuli.^2 T(k,2)^2);
Ra = sqrt((sn*V-T(k,1)).^2 T(k,2)^2); %雷达距离目标的距离
tau = 2*Ra/C;
Dfast = ones(N,1)*tm - tau'*ones(1,M);
clear tau
% myphase = pi*Kr*(Dfast-Tr/2).^2- (4*pi/lambda)*(Ra'*ones(1,M));
Srnm = Srnm sigma * exp(1i*(pi*Kr*(Dfast-Tr/2).^2- (4*pi/lambda)*(Ra'*ones(1,M))))...
.*(0<Dfast&Dfast<Tr).*(((Dslow>=(-S_qian))&(Dslow<=S_hou)).'*ones(1,M))...
;%.*(sinc(0.886/beita*(atan(tan(Theta_rc)-Dslow/T(k,2))-Theta_rc)).'*ones(1,M)).^2; %最后一项为模拟sinc^2函数
clear myphase Dfast Dslow Ra
waitbar(k/K);
end
close(h)
clear tm
%% 距离向脉冲压缩(方式1)
tr1 = -Tr/2:dt:Tr/2-dt;
Refr = zeros(1,M);
Refr(1:length(tr1)) = exp(1i*pi*Kr*tr1.^2);
%%% Refr(1:length(tr1)) = kaiser(length(tr1),4).'.*exp(1i*pi*Kr*tr1.^2); %% 加窗
Refr =ones(N,1)* Refr;
Srnm = fftshift(ifty(fty(Srnm).*(conj(fty(Refr)))),2);
f_tao = [-M/2:1:M/2-1]/M*Fsr;
% Refr = exp(1j*pi/Kr*ones(N,1)*(f_tao.^2));
%%% Refr = (ones(N,1)*kaiser(M,5).').*exp(1j*pi/Kr*ones(N,1)*(f_tao.^2)); %加窗
%%% Refr = (ones(N,1)*hamming(M).').*exp(1j*pi/Kr*ones(N,1)*(f_tao.^2)); %加窗
% Srnm = ifty(fty(Srnm).*Refr);
% R_min = R_min - Tr/2*C/2;
clear Refr tr1
% figure
% imagesc(abs(Sr)); axis xy;
%% 距离多普勒域
Sr_RD = ftx(Srnm);
% figure
% imagesc(abs(Sr_RD));axis xy;
% title('距离徙动校正前，无f_c');
%% 多普勒频率估计
f_c = 2*V*sin(Theta_rc)/lambda;
D = sqrt(1 - lambda^2/(4*V^2)*(f_c (-PRF/2:PRF/N:PRF/2-PRF/N)').^2);
Sr_RD = ftx(Srnm.*exp(-1j*2*pi*f_c*((sn).')*ones(1,M)));
mohushu = fix(f_c/PRF);
disp(['模糊数 = ',num2str(mohushu)]);
% tiaozheng = - round(mod(f_c,PRF)/PRF*size(Sr_RD,1)); %反向调整的数
% Sr_RD = circshift(Sr_RD,tiaozheng,1);
% figure
% imagesc(abs(Sr_RD));axis xy;
% title('距离徙动校正前，有f_c');
clear Srnm sn
%clear sn
Sr_f2 = fty(Sr_RD);
H_src2 = exp(1j*4*pi*R0*Fc/C*1/(2*1)*(((D.^2-1)./D.^3)*f_tao.^2)/Fc^2);
Sr_f2 = Sr_f2.*H_src2;
clear H_src2
%% 三次相位补偿
H_xiangwei3 =exp(1j*4*pi*R0*Fc/C*(-1)*(1./(2*(D.^3)) - 1./(2*(D.^5)))*f_tao.^3/Fc^3);
Sr_f2 = Sr_f2.*H_xiangwei3;
clear H_xiangwei3
%% 四次相位补偿
H_xiangwei4 = exp(1j*4*pi*R0*Fc/C*...
(-1)*(1./(8*(D.^3)) - 3./(4*(D.^5)) 5./(8*(D.^7)))*f_tao.^4/Fc^4);
Sr_f2 = Sr_f2.*H_xiangwei4;
clear H_xiangwei4
%% 五次相位补偿
H_xiangwei5 = exp(1j*4*pi*R0*Fc/C*...
(3./(8*(D.^5)) - 5./(4*(D.^7)) 7./(8*(D.^9)))*f_tao.^5/Fc^5);
Sr_f2 = Sr_f2.*H_xiangwei5;
clear H_xiangwei5
%% 六次相位补偿
H_xiangwei6 = exp(1j*4*pi*R0*Fc/C*...
(1./(16*(D.^5)) - 15./(16*(D.^7)) 35./(16*(D.^9)) - 21./(16*(D.^11)))...
*f_tao.^6/Fc^6);
Sr_f2 = Sr_f2.*H_xiangwei6;
clear H_xiangwei6
%% 七次相位补偿
H_xiangwei7 = exp(1j*4*pi*R0*Fc/C*...
(-1)*(5./(16*(D.^7)) - 35./(16*(D.^9)) 63./(16*(D.^11)) - 33./(16*(D.^13)))...
*f_tao.^7/Fc^7);
Sr_f2 = Sr_f2.*H_xiangwei7;
clear H_xiangwei7
%% 八次相位补偿
H_xiangwei8 = exp(1j*4*pi*R0*Fc/C*...
(-1)*(5./(128*(D.^7)) - 35./(32*(D.^9)) 315./(64*(D.^11)) - 231./(32*(D.^13)) 429./(128*(D.^15)))...
*f_tao.^8/Fc^8);
Sr_f2 = Sr_f2.*H_xiangwei8;
clear H_xiangwei8
%% 九次相位补偿
H_xiangwei9 = exp(1j*4*pi*R0*Fc/C*...
(35./(128*(D.^9)) - 105./(32*(D.^11)) 693./(64*(D.^13)) - 429./(32*(D.^15)) 715./(128*(D.^17)))...
*f_tao.^9/Fc^9);
Sr_f2 = Sr_f2.*H_xiangwei9;
clear H_xiangwei9
%% 十次相位补偿
H_xiangwei10 = exp(1j*4*pi*R0*Fc/C*...
(7./(256*(D.^9)) - 315./(256*(D.^11)) 1155./(128*(D.^13)) - 3003./(128*(D.^15)) 6435./(256*(D.^17)) - 2431./(256*(D.^19)))...
*f_tao.^10/Fc^10);
Sr_f2 = Sr_f2.*H_xiangwei10;
clear H_xiangwei10
%% 十一次相位补偿
H_xiangwei11 = exp(1j*4*pi*R0*Fc/C*...
(-1)*(63./(256*(D.^11)) - 1155./(256*(D.^13)) 3003./(128*(D.^15)) - 6435./(128*(D.^17)) 12155./(256*(D.^19)) - 4199./(256*(D.^21)))...
*f_tao.^11/Fc^11);
Sr_f2 = Sr_f2.*H_xiangwei11;
clear H_xiangwei11
%% 十二次相位补偿
H_xiangwei12 = exp(1j*4*pi*R0*Fc/C*...
(-1)*(21./(1024*(D.^11)) - 693./(512*(D.^13)) 15015./(1024*(D.^15)) - 15015./(256*(D.^17)) 109395./(1024*(D.^19)) - 46189./(512*(D.^21)) 29393./(1024*(D.^23)))...
*f_tao.^12/Fc^12);
Sr_f2 = Sr_f2.*H_xiangwei12;
clear H_xiangwei12
%% 十三次相位补偿
H_xiangwei13 = exp(1j*4*pi*R0*Fc/C*...
(231./(1024*(D.^13)) - 3003./(512*(D.^15)) 45045./(1024*(D.^17)) - 36465./(256*(D.^19)) 230945./(1024*(D.^21)) - 88179./(512*(D.^23)) 52003./(1024*(D.^25)))...
*f_tao.^13/Fc^13);
Sr_f2 = Sr_f2.*H_xiangwei13;
clear H_xiangwei13
%% 十四次相位补偿
H_xiangwei14 = exp(1j*4*pi*R0*Fc/C*...
(33./(2048*(D.^13)) - 3003./(2048*(D.^15)) 45045./(2048*(D.^17)) - 255255./(2048*(D.^19)) 692835./(2048*(D.^21)) - 969969./(2048*(D.^23)) 676039./(2048*(D.^25)) - 185725./(2048*(D.^27)))...
*f_tao.^14/Fc^14);
Sr_f2 = Sr_f2.*H_xiangwei14;
clear H_xiangwei14
%Sr_f2_he =Sr_f2_he Sr_f2;
%end
%Sr_RD = ifty(Sr_f2_he);
Sr_RD = ifty(Sr_f2);
clear Sr_f2 Sr_f2_he;
% figure
% imagesc(abs(Sr_RD));axis xy;
% title('SRC后，有f_c');
%
% %%
tic
%Rmin = R_min 2*Rw;
RMC = 2*Fsr/C*((1./D)*(Rmin-Rw2 (0:1:M-1)*C/2/Fsr)-R_min);
%% 单独的两个循环比一个循环的时间还要快，为什么？
for n=1:N
%RMCmaxtix(n,:) = interp1(Sr_RD(n,:),1 RMC(n,:),'spline');
Sr_RD(n,:) = interp1(Sr_RD(n,:),1 RMC(n,:),'spline');%加1的原因是因为矩阵序列是从1开始的
end
toc
tic
%% 一种优化方式
for n=1:N
for m=1:M
if RMC(n,m)<0 %%数据溢出时的处理方式
Sr_RD(n,m) = 0;
elseif (1 RMC(n,m))>M
Sr_RD(n,m:M) = zeros(1,M-m 1);
break
end
end
end
% for n=1:N
% for m=1:M
% if RMC(n,m)<0 || (1 RMC(n,m))>M %%数据溢出时的处理方式
% RMCmaxtix(n,m) = 0;
% end
% end
% end
clear RMC
toc
%
% tic
% RMCmaxtix = zeros(N,M);
% h = waitbar(0,'Sinc插值');
% P = 8; %8点sinc插值.
% P_start = P/2;
% P_chazhi = -P/2 1:1:P/2;
% for n=1:N
% for m=P_start:M
% %delta_R =(Rmin-2*RwC (m-1)*/2/Fsr)*(1/D(n)-1);
% %delta_R =(Rmin-1*Rw (m-1)*C/2/Fsr)*(1/D(n)-1); %这一步必须慎重，R(d)=R(0)/D,注意一一对应.
% %zhengduan_num = R_min - Rmin 2*Rw - (m-1)*C/2/Fsr;
% %delta_R_new = (Rmin-2*Rw (m-1)*C/2/Fsr)*(1/D(n)-1) -(R_min - Rmin 2*Rw - (m-1)*C/2/Fsr);
% %RMC = 2*Fsr/C*((Rmin-2*Rw (m-1)*C/2/Fsr)*(1/D(n)-1) -(R_min - Rmin 2*Rw - (m-1)*C/2/Fsr));
% RMC = 2*Fsr/C*((Rmin-2*Rw (m-1)*C/2/Fsr)*(1/D(n)) - R_min);
% if RMC>=4 && (floor(RMC) P/2)<=M
% %delta_R =(1/8)*(lambda/V)^2*(R_min (m)*C/2/Fsr)*((n-N/2)*PRF/N f_c)^2;
% %delta_R =(1/8)*(lambda/V)^2*(Rmin (m)*C/2/Fsr)*((n-N/2)*PRF/N)^2;
% delta_RMC = RMC - floor(RMC);
% %biaoji(n,m) = delta_RMC;
% for i=-P/2 1:1:P/2
% RMCmaxtix(n,m) = RMCmaxtix(n,m) Sr_RD(n,floor(RMC) i)*sinc(-i delta_RMC);
% end
% % RMCmaxtix(n,m) = RMCmaxtix(n,m) sum(Sr_RD(n,floor(RMC) P_chazhi).*sinc(-P_chazhi delta_RMC));
% end
% end
% waitbar(n/N);
% end
% close(h)
% clear RMC Sr_RD
% toc
% figure
% imagesc(abs(RMCmaxtix));axis xy;
% title('距离徙动校正后');
% figure
% imagesc(abs(fty(RMCmaxtix))); axis xy;
% Sr_RD = RMCmaxtix;
%% 方位向压缩
Fa=PRF;
fn = f_c (-Fa/2:Fa/N:Fa/2-Fa/N)';
%M=floor(M/6)
tau = (0:M-1)*dt 2*(Rmin - Rw2)/C;
Refa = exp(1j*4*pi/lambda*(D)*tau*C/2); %注意D减1的效果.
myRD = Sr_RD(:,1:M).*Refa;
%myRD = Sr_RD.*Refa;
clear Sr_RD tau Ka Refa D
R0 = sqrt(H^2 Yc^2); %R0为景中心距离飞机飞行轨迹的最短距离（即为垂线长度）
lingpin = exp(1j*2*pi*R0*tan(Theta_rc)/V*(fn*ones(1,M)));
Sa = iftx(myRD.*lingpin); %压到零频率位置，由于混叠，我在频域乘以一个线性相位，在时间域还原到零频位置.
clear lingpin myRD f_tao fn
[a,b]=find(abs(Sa)==max(max(abs(Sa))))
figure
imagesc(abs(Sa)); axis xy;
%% 后续处理
n_1=100;
n_2=100;
zonghe = Sa(a-n_1:1:a n_1-1,b-n_2:1:b n_2-1);
% figure
% imagesc(abs(zonghe));axis xy;
zonghe_fft2 = fftshift(fft2(zonghe));
% figure
% imagesc(abs(zonghe_fft2));axis xy;
zonghe_fft2_yidong = findpu_min(zonghe_fft2);
% zonghe_fft2_yidong = zonghe_fft2;
% figure
% imagesc(abs(zonghe_fft2_yidong));axis xy;
%% 16倍数升采样
Da = zeros(2*n_1*16,2*n_2*16);
Da(2*n_1*8-n_1:1:2*n_1*8 n_1-1,2*n_2*8-n_2:1:2*n_2*8 n_2-1)=zonghe_fft2_yidong;
figure
imagesc(abs(Da))
zuihou_result =ifft2(ifftshift(Da));axis xy;
zuihou =abs(zuihou_result)/max(max(abs(zuihou_result)));
figure
imagesc(abs(zuihou));axis xy;
figure
contour(abs(zuihou),35);
zuihou2 = 20*log10(zuihou);
figure
contour(zuihou2,linspace(-35,0,8));
fanwei = 3;
I = zuihou2(2*n_1*8-fanwei*n_1:1:2*n_1*8 fanwei*n_1-1,2*n_2*8-fanwei*n_2:1:2*n_2*8 fanwei*n_2-1); %% -3*n_1
I2 = zuihou_result(2*n_1*8-fanwei*n_1:1:2*n_1*8 fanwei*n_1-1,2*n_2*8-fanwei*n_2:1:2*n_2*8 fanwei*n_2-1);
I2 = I2/max(max(abs(I2)));
figure,set(gcf,'Color','w')
contour(I,linspace(-35,0,8));
axis square
%% 后续处理
%% 将两个方向的格点拉均匀
%%
[m,n] = size(I)
%% 变换后，两个方向上的格点距离均为fangwei_gedian
%n = round(n*juli_gedian/fangwei_gedian)
%gedian_juli = fangwei_gedian;
%% 变换后，两个方向上的格点距离均为juli_gedian
m = round(m*fangwei_gedian/juli_gedian)
gedian = juli_gedian;
I_new = imresize(I2,[m n],'bicubic');
I_new_log = 20*log10(abs(I_new)/max(max(abs(I_new))));
figure
contour(I_new_log,linspace(-35,0,8));
if m<n
if mod(m,2)==1
m=m-1;
end
I_new2 = I_new(:,round(n/2)-m/2 1:1:round(n/2) m/2);
else
if mod(n,2)==1
n=n-1;
end
I_new2 = I_new(round(m/2)-n/2 1:1:round(m/2) n/2,:);
end
I_new2_log = 20*log10(abs(I_new2)/max(max(abs(I_new2))));
figure
contour(I_new2_log,linspace(-35,0,8));axis square
% figure
% imagesc(abs(fftshift(fft2(I_new2))));
I_new3_real = imrotate(real(I_new2),Theta_rc/pi*180,'bicubic','crop');%,'crop'
I_new3_imag = imrotate(imag(I_new2),Theta_rc/pi*180,'bicubic','crop');%,'crop'
I_new3 = I_new3_real 1j*I_new3_imag;
I_new3_log = 20*log10(abs(I_new3)/max(max(abs(I_new3))));
figure
contour(I_new3_log,linspace(-35,0,8));axis square;
[m,n] = find(I_new3_log ==max(max(I_new3_log)))
% figure
% plot(I_new3_log(m,:)); title('距离向剖面图');
% axis([-inf inf -40 0]);grid on;
% figure
% plot(I_new3_log(:,n)); title('方位向剖面图');
% axis([-inf inf -40 0]);grid on;
%
a_pow_new = I_new3(:,n);
r_pow_new = I_new3(m,:);
% figure
% plot(angle(r_pow_new)*180/pi); title('距离向剖面图');
% figure
% plot(angle(a_pow_new)*180/pi); title('方位向剖面图');
[a_PSIR a_ISIR a_IRW]=PSISIRIRW(a_pow_new,gedian/16,fangwei_fenbian) %除以16是因为进行了16倍升采样
[r_PSIR r_ISIR r_IRW]=PSISIRIRW(r_pow_new,gedian/16,juli_fenbian) %除以16是因为进行了16倍升采样