sarimgdata.asv

条带式合成孔径雷达的cs。RDA。RMA算法 作者何志华

%%%%%%%%%%%%%%%%%%%%%%%%
%SAR点目标Stripmap 模式下回波生成和成像仿真
%NUDT,Hezhihua
%E-mail:[e-mail]
%last update:8/11/2005

clear all
%%range:x domain
%x=c*t/2,kx=2*f/c
Tx=200;%时宽200m(1.33us)
Bx=1;%带宽1(1/m)(150MHz)
a=Bx/Tx;%调频斜率
kxc=10;%载频10(1/m)(1.5GHz)
Nx=512;
Xc=3000;X0=150;
x=Xc+linspace(-X0,X0,Nx);%x域序列:Xc-X0~Xc+X0
dx=2*X0/Nx;
kx=linspace(-1/dx/2,1/dx/2,Nx);%kx域序列

%%cross-range:y domain
Ty=300;%时宽300m,合成孔径长度
By=1;%带宽1(1/m)
b=kxc/Xc;%调频斜率 b=By/Ty=kxc/Xc
Ny=1024;
Y0=200;
y=linspace(-Y0,Y0,Ny);%y域序列:-Y0~Y0
dy=2*Y0/Ny;
ky=linspace(-1/dy/2,1/dy/2,Ny);%ky域序列

%%target geometry
Ntar=3;%目标个数
Ptar=[Xc,0,1+0j  %x坐标,y坐标,复后向散射系数             
          Xc+50,-50,1+0j
          Xc+50,50,1+0j];
      
%%SAR raw data
s_xy=zeros(Nx,Ny);
for i=1:1:Ntar
    xn=Ptar(i,1);yn=Ptar(i,2);sigma=Ptar(i,3);
    X=x'*ones(1,Ny);%扩充为矩阵
    Y=ones(Nx,1)*y;
    DX=X-sqrt(xn^2+(Y-yn).^2);
    phase=pi*a*DX.^2-2*pi*kxc*sqrt(xn^2+(Y-yn).^2);
    s_xy=s_xy+sigma*exp(j*phase).*(abs(DX)<Tx/2).*(abs(Y-yn)<Ty/2);%SAR回波信号
end
disp('SAR回波信号完毕...')

%%产生相关的数据
p0_x=exp(j*pi*a*(x-Xc).^2).*(abs(x-Xc)<Tx/2);%距离向LFM信号
p0_kx=fftshift(fft(fftshift(p0_x)));
p0_y=exp(-j*pi*b*y.^2).*(abs(y)<Ty/2);%方位向LFM信号
p0_ky=fftshift(fft(fftshift(p0_y)));
 
s_kxy=fftshift(fft(fftshift(s_xy)));%对raw data距离向FFT
sxc_kxy=s_kxy.*(conj(p0_kx).'*ones(1,Ny));
sxc_kxky=fftshift(fft(fftshift(sxc_kxy).')).';%距离压缩后的2D频域信号
sxc_xy=fftshift(ifft(fftshift(sxc_kxy)));%距离压缩后的信号

sxc_xky=fftshift(fft(fftshift(sxc_xy).')).';%距离压缩后，距离-多普勒域
%sxcyc_xky=sxc_xky.*(ones(Nx,1)*conj(p0_ky));
%sxcyc_xy=fftshift(ifft(fftshift(sxcyc_xky).')).';%距离，方位压缩后的信号

%load sarecho.mat x y kx ky Xc dx;
%load sarecho.mat kxc Nx Ny ;
%load sarecho.mat sxc_xky p0_ky;
%%距离迁移校正
RM_xky=(x'*ones(1,Ny)).*(kxc./sqrt(kxc^2-(ones(Nx,1)*ky).^2)-1);%距离迁移量

method=1;%%选取不同的插值方法

if method==1
%%最邻域近似插值
RMC=zeros(Nx,Ny);
RM=round(RM_xky/dx);
for i=1:Nx-max(max(RM))
    for j=1:Ny
        RMC(i,j)=sxc_xky((i+RM(i,j)),j);
    end
end

elseif method==2
%%截断sinc插值
P=8/2;%8点sinc插值
RMC=zeros(Nx,Ny);
RM=floor(RM_xky/dx);
for i=P:Nx-max(max(RM))-P
     for j=1:Ny
       T=RM(i,j)+(-P+1:P);
       f=sxc_xky(i+T,j).';
       RMC(i,j)=sum(f.*sinc(RM_xky(i,j)/dx-T));
     end
     i
end

else
%%线性插值
RMC=zeros(Nx,Ny);
RM=floor(RM_xky/dx);
for i=2:Nx-max(max(RM))-1
    for j=1:Ny
        f1=sxc_xky(i+RM(i,j),j);
        f2=sxc_xky(i+RM(i,j)+1,j);
        RMC(i,j)=f1+(f2-f1)*(RM_xky(i,j)-RM(i,j)*dx)/dx;
    end
end

end

sxcrmc_xky=RMC;
%%方位压缩
sxcrmcyc_xky=sxcrmc_xky.*(ones(Nx,1)*conj(p0_ky));
sxcrmcyc_xy=fftshift(ifft(fftshift(sxcrmcyc_xky).')).';

%save rda_img.mat sxcrmcyc_xy x y Xc Nx Ny%存储变量
%load rda_img.mat sxcrmcyc_xy x y Xc Nx Ny;
%%是否显示数据
display=1;
if display==1
colormap(gray)
G=10*log10(abs(sxcrmcyc_xy));
subplot(221)
imagesc(y,x-Xc,G)
xlabel('y domain(cross-range),m')
ylabel('x domain(range),m')
title('|s_x_c_r_m_c_y_c(x,k_y)|')
axis tight

subplot(222)
imagesc(y,x-Xc,G)
xlabel('y domain(cross-range),m')
ylabel('x domain(range),m')
title('|s_x_c_r_m_c_y_c(x,k_y)|')
axis([-40,40,-20,20])

Range=sxcrmcyc_xy(:,512).';
Azimuth=sxcrmcyc_xy(256,:).';
Range=20*log10(abs(Range)+1e-6);
Azimuth=20*log10(abs(Azimuth)+1e-6);
subplot(223)
plot(x-Xc,Range-max(Range))
xlabel('x domain,m')
ylabel('Magnitude,dB')
title('Range')
axis([-5,5,-30,0])
grid on

subplot(224)
plot(y,Azimuth-max(Azimuth))
xlabel('y domain,m')
ylabel('Magnitude,dB')
title('Azimuth')
axis([-5,5,-30,0])
grid on
end