clc;
clear all;
close all;
%(1)定义参数
c=3e8;
j=sqrt(-1);
pi=3.14159;
fc=5.3e6;
lamda=c/fc;%波长
D=4; %方位向的天线尺寸
Va=150;%卫星平台速度
Kr=20e12;%距离向脉冲的调频率
Tr=2.5e-6;%脉冲持续时间
sq_ang=3.5/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;%方位复采样率(PRF)
R_near=2e4;%先理解为距目标视野最近距离
R_far=R_near+1000;%先理解为距目标视野最远距离
%R=(R_near+R_far)/2; %assumed a fixed R for implement
%Y_min=-50; %both R_near and R_far are closest approach range
%Y_max=Y_min+100;
La_near=0.886*R_near*lamda/cos(sq_ang)^2/D;%near range的合成孔径长度
La_far=0.886*R_far*lamda/cos(sq_ang)^2/D;%far range的合成孔径长度
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);%shortest range between radar and scene
Rmax=sqrt(R_far^2+(Tc_far*Va-La_far/2)^2);%longest range between radar and scene
%参数定义完了 然后开始disp:用来展示变量的内容,可以是字符串,元胞,矩阵,结构体
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);
%(2)回波模型(怎么用feko来解决导入呢 这里全是点目标的样子。。)
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);
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];%Rpt is there closest approach range
La=0.886*Rpt*lamda/(cos(sq_ang)^2)/D;%每个目标的合成孔径长度
Tc=-Rpt*tan(sq_ang)/Va;%beam center cross time of each target
Npt=length(Rpt);
Y_high=max(Ypt)+50;%%确定成像网格的范围,以使其包括目标点
Y_low=min(Ypt)-50;%%成像区域是(R_right-R_left)*(Y_high-Y_low)的一片区域
R_left=R_near-50;
R_right=R_far+50;
disp('number of point target:');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);
sig=sig+exp(j*pi*Kr*Dr.^2-j*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('回波数据幅度');xlable('距离向(采样点)');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);
filter_r=ones(Na,1)*exp(j*pi*fr.^2/Kr);
sig_rd=sig_rd.*filter_r;
nup=100;%距离向升采样系数为100
Nr_up=Nr*nup;
nz=Nr_up-Nr;
dtr=1/nup/Fr;
sig_rd_up=zeros(Na,Nr_up);
sig_rd_up=[sig_rd(:,1:Nr/2),zeros(Na,nz),sig_rd(:,(Nr/2+1):Nr)];
sig_rdt=ifft(sig_rd_up,[],2);
%figure('Name','距离向匹配滤波后‘);imagesc(abs(sig_rdt));
%title('距离向匹配滤波后’);xlabel('距离向(采样点)');ylabel('方位向(采样点)'0;
%%分别压至(2*sqrt(R_near^2+(Va*Tc_near)^2)/c-(2*Rmin/c-Tr/2))*Fr;
%%----(2*sqrt((R_near+500)^2+(Va*Tc(2))^2)/c-(2*Rmin/c-Tr/2)*Fr;
%%----(2*sqrt((R_near+1000)^2+(Va*Tc(3))^2)/c-(2*Rmin/c-Tr/2))*Fr;
%步骤二 对成像区域进行网格化(Na*Nr)
R=zero(1,Nr);
for ii=1:Nr
%R(1,ii)=R_near+(R_far-R_near)/(Nr-1)*(ii-1);
R(1:ii)=R_left+(R_right-R_left)/(Nr-1)*(ii-1);
end
Y=zeros(1,Na);
for ii=1:Na
Y(1:ii)=Y_low+(Y_high-Y_low)/(Na-1)*(ii-1);
end
R=ones(Na,1)*R;
Y=Y'*ones(1,Nr);
%%步骤三:根据时延在回波域寻找相应的位置并进行成像
f_back=zeros(Na,Nr);
for ii=1:Na
% ii=100;
R_ij=sqrt(R.^2+(Y-Va*ta(ii).^2));
t_ij=2*R_ij/c;
% t_ij=round((t_ij-tr(Nr/2+1))/dtr)+Nr_up/2+1;%%怎么将时域转化为点数
t_ij=round((t_ij-(2*Rmin/c-Tr/2))/dtr);
%t_ij=round(t_ij/dtr);
it_ij=(t_ij>0&t_ij<=Nr_up);
t_ij=t_ij.*it_ij+Nr_up*(1-it_ij);%%将矩阵中为零的元素用Nr_up代替
sig_rdta=sig_rdt(ii,:);
sig_rdta(Nr_up)=0;
f_back=f_back+sig_rdta(t_ij).*exp(li*4*pi*R_ij/lamda);
end
figure('Name','BP算法处理后');imagesc(abs(f_back));
title('BP算法处理后');xlabel('距离向');ylabel('方向位');