stripmap.m

SAR雷达基本信号处理算法

G=abs(fs)';
xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
image(k*c/pi2,ku,256-cg*(G-ng));
axis('square');axis('xy')
xlabel('Fast-time Frequency \omega/(2*pi), Hertz')
ylabel('Synthetic Aperture (Slow-time) Frequency k_u, rad/m')
title('Stripmap SAR Signal Spectrum after Digital Spotlighting')
print P6.5.ps
pause(1)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%    SLOW-TIME DOPPLER SUBSAMPLING     %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if Y0 < L,
 ny=2*ceil(1.2*Y0/du);      % Number of samples in y domain
                            % with 20 percent guard band
 ms=floor(m/ny);            % subsampling ratio
 tt=floor(m/(2*ms));
 I=m/2+1-tt*ms:ms:m/2+1+(tt-1)*ms; % subsampled index in ku domain
 [tt,ny]=size(I);           % number of subsamples
 fs=fs(:,I);                % subsampled SAR signal spectrum
 ky=ku(I);                  % subsampled ky array
 dky=dku*ms;                % ky domain sample spacing
else,
 dky=dku;
 ny=m;
 ky=ku;
end;

dy=pi2/(ny*dky);            % y domain sample spacing
y=dy*(-ny/2:ny/2-1);        % cross-range array

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%             RECONSTRUCTION           %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
ky=ones(n,1)*ky;       % ky array
kx=(4*k(:).^2)*ones(1,ny)-ky.^2;
kx=sqrt(kx.*(kx > 0));                  % kx array
%
plot(kx(1:20:n*ny),ky(1:20:n*ny),'.')
xlabel('Spatial Frequency k_x, rad/m')
ylabel('Spatial Frequency k_y, rad/m')
title('Stripmap SAR Spatial Frequency Data Coverage')
axis image; axis xy
print P6.6.ps
pause(1)
%
kxmin=min(min(kx));
kxmax=max(max(kx));
dkx=pi/X0;        % Nyquist sample spacing in kx domain
nx=2*ceil((.5*(kxmax-kxmin))/dkx); % Required number of
                      % samples in kx domain;
                      % This value will be increased slightly
                      % to avoid negative array index
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                                                         %%%
%%%   FIRST TWO OPTIONS FOR RECONSTRUCTION:                 %%%
%%%                                                         %%%
%%%     1. 2D Fourier Matched Filtering and Interpolation   %%%
%%%     2. Range Stacking                                   %%%
%%%                                                         %%%
%%%     Note: For "Range Stacking," make sure that the      %%%
%%%           arrays nx, x, and kx are defined.             %%%
%%%                                                         %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%    2D FOURIER MATCHED FILTERING AND INTERPOLATION    %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Matched Filtering
%
fs0=(kx > 0).*exp(cj*kx*Xc+cj*.25*pi ...
           -cj*2*k(:)*ones(1,ny)*Xc); % reference signal complex conjugate
fsm=fs.*fs0;     % 2D Matched filtering

% Interpolation
%
is=8;       % number of neighbors (side lobes) used for sinc interpolator
I=2*is+1;
kxs=is*dkx; % plus/minus size of interpolation neighborhood in KX domain
%
nx=nx+2*is+4;  % increase number of samples to avoid negative
               %  array index during interpolation in kx domain
KX=kxmin+(-is-2:nx-is-3)*dkx;     % uniformly-spaced kx points where
                                  % interpolation is done
kxc=KX(nx/2+1);                   % carrier frequency in kx domain
KX=KX(:)*ones(1,ny);
%
F=zeros(nx,ny);         % initialize F(kx,ky) array for interpolation

for i=1:n;                       % for each k loop
  i                              % print i to show that it is running
 icKX=round((kx(i,:)-KX(1,1))/dkx)+1; % closest grid point in KX domain
 cKX=KX(1,1)+(icKX-1)*dkx;            % and its KX value
 ikx=ones(I,1)*icKX+[-is:is]'*ones(1,ny);
 ikx=ikx+nx*ones(I,1)*[0:ny-1];
 nKX=KX(ikx);
 SINC=sinc((nKX-ones(I,1)*kx(i,:))/dkx);             % interpolating sinc
 HAM=.54+.46*cos((pi/kxs)*(nKX-ones(I,1)*kx(i,:)));  % Hamming window
         %%%%%   Sinc Convolution (interpolation) follows  %%%%%%%%
 F(ikx)=F(ikx)+(ones(I,1)*fsm(i,:)).*(SINC.*HAM);
end
%
%  DISPLAY interpolated spatial frequency domain image F(kx,ky)

KX=KX(:,1).';
KY=ky(1,:);

G=abs(F)';
xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
image(KX,KY,256-cg*(G-ng));
axis image; axis xy
xlabel('Spatial Frequency k_x, rad/m')
ylabel('Spatial Frequency k_y, rad/m')
title('Wavefront Stripmap SAR Reconstruction Spectrum')
print P6.7.ps
pause(1)

%
f=iftx(ifty(F));     % Inverse 2D FFT for spatial domain image f(x,y)
%
dx=pi2/(nx*dkx);     % range sample spacing in reconstructed image
x=dx*(-nx/2:nx/2-1); % range array
%
% Display SAR reconstructed image

G=abs(f)';
xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
image(Xc+x,y,256-cg*(G-ng));axis([Xc-X0 Xc+X0 -Y0 Y0]);
axis image; axis xy
xlabel('Range X, meters')
ylabel('Cross-range Y, meters')
title('Wavefront Stripmap SAR Reconstruction')
print P6.8.ps
pause(1)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%      RANGE STACK WAVEFRONT RECONSTRUCTION       %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
f_stack=zeros(nx,ny); % Initialize reconstruction array in (x,y) domain
for i=1:nx; i        % Stack's loop for reconstruction at each range
 f_stack(i,:)=ifty(sum(fs.*exp(cj*kx*(Xc+x(i)) ...
  +cj*.25*pi-cj*2*k(:)*ones(1,ny)*Xc)));
end;

% Remove carrier in range domain
f_stack=f_stack.*exp(-cj*x(:)*kxc*ones(1,ny));
%
f_stack=f_stack/nx; % Scale it for comparison with Fourier interpolation
                    % Use "f_stack-f" to display difference of two
                    % reconstructions
           
G=abs(f_stack)';
xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
image(Xc+x,y,256-cg*(G-ng));axis([Xc-X0 Xc+X0 -Y0 Y0]);
axis image; axis xy
xlabel('Range X, meters')
ylabel('Cross-range Y, meters')
title('Range Stack Stripmap SAR Reconstruction')
print P6.9.ps
pause(1)
                
F_stack=ftx(fty(f_stack)); % Reconstruction array in spatial frequency
                           % domain; Use "F_stack-F" to display
                           %  difference of two reconstructions
%
G=abs(F_stack)';
xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
image(KX,KY,256-cg*(G-ng));
axis image; axis xy
xlabel('Spatial Frequency k_x, rad/m')
ylabel('Spatial Frequency k_y, rad/m')
title('Range Stack Stripmap SAR Reconstruction Spectrum')
print P6.10.ps
pause(1)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%     TIME DOMAIN CORRELATION RECONSTRUCTION      %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
f_tdc=zeros(nx,ny); % Initialize reconstruction array in (x,y) domain

for i=1:nx; i
   Bmaxi=(Xc+x(i))*tan(asin(lambda/D));  % Half-beamwidth at range x(i)
   for j=1:ny;
      Iu=(ones(n,1)*abs(y(j)-u) <= Bmaxi(:)*ones(1,m));
      t_ij=(2*sqrt((x(i)+Xc)^2+(y(j)-u).^2))/c;
      f_tdc(i,j)=sum(sum(s_ds.*exp(cj*w(:)*(t_ij-tm(n/2+1))).* ...
         (Iu & ones(n,1)*(t_ij >= Ts & t_ij <= Tf))));
   end;
end;

% Remove carrier in range domain
f_tdc=f_tdc.*exp(-cj*x(:)*kxc*ones(1,ny));
%
           
G=abs(f_tdc)';
xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
image(Xc+x,y,256-cg*(G-ng));axis([Xc-X0 Xc+X0 -Y0 Y0]);
axis image; axis xy
xlabel('Range X, meters')
ylabel('Cross-range Y, meters')
title('TDC Stripmap SAR Reconstruction')
print P6.11.ps
pause(1)
                
F_tdc=ftx(fty(f_tdc));     % Reconstruction array in spatial frequency
%
G=abs(F_tdc)';
xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
image(KX,KY,256-cg*(G-ng));
axis image; axis xy
xlabel('Spatial Frequency k_x, rad/m')
ylabel('Spatial Frequency k_y, rad/m')
title('TDC Stripmap SAR Reconstruction Spectrum')
print P6.12.ps
pause(1)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%           BACKPROJECTION RECONSTRUCTION         %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
f_back=zeros(nx,ny); % Initialize reconstruction array in (x,y) domain
n_ratio=100;              % Upsampling ratio in fast-time domain
nu=n_ratio*n;             % Size of upsampled s(t,u) array in t domain
nz=nu-n;                  % Number of zeros
dtu=(n/nu)*dt;            % Fast-time sample spacing of upsampled array
tu=dtu*(-nu/2:nu/2-1);    % Upsampled reference fast-time array
X=x(:)*ones(1,ny);
Y=ones(nx,1)*y;

for j=1:m; j
   t_ij=(2*sqrt((X+Xc).^2+(Y-u(j)).^2))/c;
   t_ij=round((t_ij-tm(n/2+1))/dtu)+nu/2+1;
   it_ij=(t_ij > 0 & t_ij <= nu);
   t_ij=t_ij.*it_ij+nu*(1-it_ij);
   S=ifty([zeros(1,nz/2),s_ds(:,j).',zeros(1,nz/2)])...
      .*exp(cj*wc*tu);                % Upsampled data
   S(nu)=0;
   Bmax_xy=(Xc+X)*tan(asin(lambda_max/D));   % Maximum half-beamwidth
   Iu=(abs(Y-u(j)) <= Bmax_xy);
   f_back=f_back+S(t_ij).*Iu;
end;

clear X Y

% Remove carrier in range domain
f_back=f_back.*exp(-cj*x(:)*kxc*ones(1,ny));
%
  
G=abs(f_back)';
xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
image(Xc+x,y,256-cg*(G-ng));axis([Xc-X0 Xc+X0 -Y0 Y0]);
axis image; axis xy
xlabel('Range X, meters')
ylabel('Cross-range Y, meters')
title('Backprojection Stripmap SAR Reconstruction')
print P6.13.ps
pause(1)
                
F_back=ftx(fty(f_back));     % Reconstruction array in spatial frequency
%
G=abs(F_back)';
xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
image(KX,KY,256-cg*(G-ng));
axis image; axis xy
xlabel('Spatial Frequency k_x, rad/m')
ylabel('Spatial Frequency k_y, rad/m')
title('Backprojection Stripmap SAR Reconstruction Spectrum')
print P6.14.ps
pause(1)