sarautofocus.m

本程序是一个基于matlab的雷达SAR成像自聚焦算法

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    　 STRIPMAP SAR　AUTOFOCUS SIMULATION          %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear

%%%%%%%%%%%%%%% constant parameters %%%%%%%%%%%%%%%%%%%

c=3e8;                   % propagation speed,3*10^8m/s

%%%%%%%%%%%%%%% chirp signal parameters %%%%%%%%%%%%%%%
fc=1.5e9;                % carrier frequency,1.5GHz
wc=2*pi*fc;
lambda_c=c/fc;           % Wavelength at carrier frequency
Tr=1.5e-6;               % chirp pusle duration 1.5us
Br=150e6;                % chirp frequebcy bandwidth 150Mhz
Kr=Br/Tr;                % chirp rate 
%%%%%%%%%%%%%% observating strip parameters %%%%%%%%%%%%
Rc=3000;                 % Range distance to center of target area,3000m
R0=150;                  % Target area in range is within [Rc-R0,Rc+R0]
X0=200;                  % Target area in cross-range is within [-X0,X0]

%%%%%%%%%%%% range ,r/t domain parameters %%%%%%%%%%%%%%%
Nr=512;                  % The number of r (fast-time) domain samples
r=Rc+linspace(-R0,R0,Nr);
t=2*r/c;                 % t domain series (fast-time domain) 
dt=4*R0/c/Nr;            % sample spacing in fast-time domain (t domain) 
f=linspace(-1/2/dt,1/2/dt,Nr);% f domain (FT about t)

%%%%%%%%%%%Azimuth, cross range parameters x/u domian %%%%%%
v=100;                   % speed of the SAR aircraft  100m/s
theta=0;                 % squint angle of beam center
Lsar=300;                % synthesize aperture length 300m
Tsar=Lsar/v;             % synthesize aperture time
Na=1024;                 % The number of x (slow-time) domain samples
x=linspace(-X0,X0,Na);   % x domian series (cross-range)
u=x/v;                   % u (slow-time) domain
du=2*X0/v/Na;            % sample spacing in slow-time domain (u domain)
fu=linspace(-1/2/du,1/2/du,Na); % fu series(FT about u domain)
fdc=2*v*sin(theta)/lambda_c;    % the centric Doppler frequency 
fr=-2*v^2/(lambda_c*Rc);        % Doppler chirp rate
Ba=abs(fr)*Tsar;                % Doppler frequency bandwidth Ba

%%%%%%%%%%%%% SAR　Resolution  %%%%%%%%%%%%%%%%%%%%%%%%%
Dr=c/2/Br;                       %range resolution
DX=v/Ba;                         %cross-range resolution

%%%%%%%%%%%%% target parameters ( one point target) %%%%%%%%%%%%%%%%%%%%%%%
% yn: range;            xn= cross-range;   fn: reflectivity
rn=Rc;                  xn=0;              fn=1 ;       

%%%%%%%%%%%%%%  targets echo %%%%%%%%%%%%%%%%%%
s_tu=zeros(Na,Nr);   
U=u'*ones(1,Nr);       % extend to matrix
T=ones(Na,1)*t;        
R=zeros(Na,Nr);
DT=zeros(Na,Nr);             
R=sqrt(rn^2+(xn-v*U).^2);
DT=T-2*R/c;
phase=pi*Kr*DT.^2-4*pi*R/lambda_c;
s_tu=fn*exp(j*phase).*(abs(DT)<Tr/2).*((abs(xn-v*U))<Lsar/2);

%%%%%%%%%%%%%%%  Range Compression %%%%%%%%%%%%%%%
p0_t=exp(j*pi*Kr*(t-2*Rc/c).^2).*(abs(t-2*Rc/c)<Tr/2);   % Range referrence signal
p0_f=fftshift(fft(fftshift(p0_t)));
s_fu=fftshift(fft(fftshift(s_tu.'))).';                  % Range FFT
src_fu=s_fu.*(ones(Na,1)*conj(p0_f));                    % Range Compression
src_tu=fftshift(ifft(fftshift(src_fu.'))).';             % signal after Range Compression
src_tfu=fftshift(fft(fftshift(src_tu)));                 % Azimuth FFT and Range-Doppler domain

%%%%%%%%%%%%%%%  RCM Compensation  %%%%%%%%%%%%%%%
src_ffu=fftshift(fft(fftshift(src_fu)));                 % 2D spectrum after RC
F=ones(Na,1)*f;
FU=fu'*ones(1,Nr);
p0_2f=exp(j*pi/fc^2/fr*(FU.*F).^2+j*pi*fdc^2/fc/fr*F-j*pi/fc/fr*FU.^2.*F); % RCM Compensation
s2rc_ffu=(src_ffu).*p0_2f;
s2rc_tfu=fftshift(ifft(fftshift(s2rc_ffu).')).';        % Range-Doppler domain
s2rc_tu=fftshift(ifft(fftshift(s2rc_tfu)));             % t-u domain

%%%%%%%%%%%%%% Azimuth Compression %%%%%%%%%%%%%%%
%  p0_fu=exp(j*pi/fr*(FU-fdc).^2);                         % Azimuth referrence signal in Doppler domain
%  s2rcac_tfu=s2rc_tfu.*p0_fu;                             % Azimuth Compression
%  s2rcac_tu=fftshift(ifft(fftshift(s2rcac_tfu)));         % Amizuth IFFT
 
% %%%%%%%%%%%%%% AutoFocus by PGA algorithm %%%%%%%%%%%%%%%%%%%%
% threshold=1e-3*pi;
% Nw=127;                                                  % the length of window Nw
% window=zeros(Na,1);                                      % the triangle window
% TriWindow=triang(Nw);                                    
% window((Na/2-(Nw-1)/2):Na/2)=TriWindow(1:(Nw+1)/2);
% window(Na/2+1:Na/2+(Nw+1)/2)=TriWindow((Nw+1)/2:end);
% dfr=pi;                                                      % the variety of Doppler chirp rate
% count=0;
% while abs(dfr)>threshold
%     count=count+1
%     DistinctPoint=max(max(abs(s2rcac_tu)));                  % find the DistinctPoint
%     [index_a index_r]=find(abs(s2rcac_tu)==DistinctPoint);   % the  range unit index_r
%     gn=s2rcac_tu(:,index_r);
%     gn=circshift(gn,Na/2-index_a);                           % shift the DistinctPoint to the center
%     gn=gn.*window;
%     n=[-Na/2:-1 1:Na/2];
%     nk=[Na/2-(Nw-1)/2:Na/2+(Nw-1)/2+1];
%     nkk=-(Nw-1)/2:(Nw-1)/2+1;
%     gdn=-j*n'.*gn;                                            %  the differential  of  signal after RC,RCM
%     Gdk=fftshift(fft(gdn));
%     Gk=fftshift(fft(gn));
%     dphase=imag(Gdk(nk).*conj(Gk(nk)))./(Gk(nk).*conj(Gk(nk))); % phase gradient 
%     [a S]=polyfit(nkk,dphase',1);                              % estimate the slope rate a(1) by  LS
%     dfr=a(1)*du*(fr)^2;                                        % the variance of Doppler chirp rate
%     %%%%%%%%% Azimuth Compression %%%%%%%%%%
%     fr=fr+dfr;
%     p0_fu=exp(j*pi/fr*(FU-fdc).^2);                         % Azimuth referrence signal in Doppler domain
%     s2rcac_tfu=s2rc_tfu.*p0_fu;                             % Azimuth Compression
%     s2rcac_tu=fftshift(ifft(fftshift(s2rcac_tfu)));         % Amizuth IFFT
% end

%%%%%%%%%% Minimun Entropy Algorithm %%%%%%%%%%%%%%
frc=-2*v^2/(lambda_c*Rc);        % intial 
fr0=abs(-2*v^2/(lambda_c*Rc))*0.2;   % search region [frc-fr0 frc+fr0]
Ns=128;                               % number of sampling
frs=linspace(frc-fr0,frc+fr0,Ns);
H=zeros(1,Ns);
for i=1:Ns
    %%%%%%%%%%%%%%% Azimuth Compression %%%%%%%%%%%%%%%
    p0_fu=exp(j*pi/frs(i)*(FU-fdc).^2);                     % Azimuth referrence signal in Doppler domain
    s2rcac_tfu=s2rc_tfu.*p0_fu;                             % Azimuth Compression
    s2rcac_tu=fftshift(ifft(fftshift(s2rcac_tfu)));         % Amizuth IFFT
    s2rcac_X=abs(s2rcac_tu);
    s2rcac_X=s2rcac_X./(ones(Na,1)*sum(s2rcac_X));
    H(i)=-sum(sum(s2rcac_X.*log(s2rcac_X)));
end
 [Hmin index]=min(H);
 p0_fu=exp(j*pi/frs(index)*(FU-fdc).^2);                 % the minumum entropy Doppler chirp rate
 s2rcac_tfu=s2rc_tfu.*p0_fu;                             % Azimuth Compression
 s2rcac_tu=fftshift(ifft(fftshift(s2rcac_tfu))); 
    



%%%%%%%%%%%%%%%%    Image show    %%%%%%%%%%%%%%%%%

% Show Measured Stripmap SAR Signal
figure(1)
G=20*log10(abs(s_tu)+1e-6);
xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
colormap(gray)
image(r-Rc,x,256-cg*(G-ng));
axis([-R0 R0 -X0 X0])
grid on,axis tight
xlabel('Range')
ylabel('Azimuth')
title('Measured Stripmap SAR Signal')

% Image after range compression  
figure(2)
G=20*log10(abs(src_tu)+1e-6);
colormap(gray)
xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
image(r-Rc,x,256-cg*(G-ng));
axis([-R0 R0 -X0 X0])
grid on,axis tight
xlabel('Range')
ylabel('Azimuth')
title('Image after RC ')

% Image after RCM 
figure(3)
G=20*log10(abs(s2rc_tu)+1e-6);
colormap(gray)
xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
image(r-Rc,fu,256-cg*(G-ng));
axis([-R0 R0 -X0 X0])
grid on,axis tight
xlabel('Range')
ylabel('Doppler frequency')
title(' Image after RCM')

% show image after RC,AC 
figure(4)
G=20*log10(abs(s2rcac_tu)+1e-6);
colormap(gray)
xg=max(max(G)); ng=xg-60; cg=255/(xg-ng);
imagesc(r-Rc,x,cg*(G-ng).*(G>ng));
axis([-R0 R0 -X0 X0])
grid on,axis tight
xlabel('Range')
ylabel('Azimuth')
title('Targets Image ')