stripmap1.m

合成孔径雷达系统MATLALB模拟实现源码。

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %    PULSED STRIPMAP SAR SIMULATION AND RECONSTRUCTION   %
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


colormap(gray(256))
cj=sqrt(-1);
pi2=2*pi;
%
c=3e8;                   % propagation speed
f0=50e6;                 % baseband bandwidth is 2*f0
w0=pi2*f0;
fc=100e6;                % carrier frequency
wc=pi2*fc;
lambda_min=c/(fc+f0);    % Wavelength at highest frequency
lambda_max=c/(fc-f0);    % Wavelength at lowest frequency
kc=(pi2*fc)/c;           % wavenumber at carrier frequency
kmin=(pi2*(fc-f0))/c;    % wavenumber at lowest frequency
kmax=(pi2*(fc+f0))/c;    % wavenumber at highest frequency
%
Xc=500;                  % Range distance to center of target area
X0=200;                  % Target area in range is within [Xc-X0,Xc+X0]
Y0=300;                  % Target area in cross-range is within [-Y0,Y0]
%
D=2*lambda_max;          % Diameter of planar radar
                         % (Beamwidth of a planar radar aperture varies
                         %  with fast-time frequency)
%
Bmin=(Xc-X0)*tan(asin(lambda_min/D));    % Minimum half-beamwidth
Bmax=(Xc+X0)*tan(asin(lambda_max/D));    % Maximum half-beamwidth
L=Bmax+Y0;                          % Synthetic aperture length is 2*L

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%    u domain parameters and arrays for SAR signal     %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
du=D/4;                           % sample spacing in aperture domain
du=du/1.2;                        % 10 percent guard band                    
m=2*ceil(L/du);                   % number of samples on aperture
dku=pi2/(m*du);                   % sample spacing in ku domain
u=du*(-m/2:m/2-1);                % synthetic aperture array
ku=dku*(-m/2:m/2-1);              % ku array

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% u domain parameters and arrays for compressed SAR signal %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
x_t=[Xc-X0:2*X0:Xc+X0];
lambda_t=[lambda_min:lambda_max-lambda_min:lambda_max];
B_t=x_t(:)*tan(asin(lambda_t/D));
L_t=B_t+Y0;
duc=min(min((x_t(:)*lambda_t)./(4*L_t))); % sample spacing in u domain
                                          % for compressed SAR signal
                                          % Less restrictive: use duc=du or
                                          % duc=du/2
clear x_t lambda_t
dkuc=dku;
mc=2*ceil(pi/(duc*dkuc));      % number of samples on aperture
mc=2^ceil(log(mc)/log(2));     % is chosen to be a power of 2
duc=pi2/(mc*dkuc);
uc=duc*(-mc/2:mc/2-1);         % synthetic aperture array
dkuc=pi2/(mc*duc);             % sample spacing in ku domain
kuc=dkuc*(-mc/2:mc/2-1);       % kuc array

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%       Fast-time domain parameters and arrays         %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
Tp=1.5e-7;                     % Chirp pulse duration
alpha=w0/Tp;                   % Chirp rate
wcm=wc-alpha*Tp;               % Modified chirp carrier
%
Rmin=Xc-X0;
Rmax=sqrt((Xc+X0)^2+Bmax^2);
Ts=(2/c)*Rmin;                 % start time of sampling
Tf=(2/c)*Rmax+Tp;              % end time of sampling
T=Tf-Ts;                       % fast-time interval of measurement
Ts=Ts-.1*T;                    % start slightly earlier (10% guard band)
Tf=Tf+.1*T;                    % end slightly later (10% guard band)
T=Tf-Ts;
theta_ax=asin(lambda_min/D);
Tmin=max(T,(4*X0)/ ...
         (c*cos(theta_ax))); % Minimum required fast-time interval
%
dt=1/(4*f0);                 % Time domain sampling (guard band factor 2)
n=2*ceil((.5*Tmin)/dt);      % number of time samples
t=Ts+(0:n-1)*dt;             % time array for data acquisition
dw=pi2/(n*dt);               % Frequency domain sampling
w=wc+dw*(-n/2:n/2-1);        % Frequency array (centered at carrier)
k=w/c;                       % Wavenumber array
Ik=find(k >= kmin & k <= kmax);
[temp,nk]=size(Ik);
lambda=zeros(1,n);
lambda(Ik)=(pi2*ones(1,nk))./k(Ik);      % Wavelength array
phi_d=asin(lambda/D);                    % Divergence angle
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%       Resolution         %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

DX=c/(4*f0);                      % range resolution
DY=D/2;                           % cross-range resolution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%           Parameters of Targets                 %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
ntarget=9;                        % number of targets
% Set ntarget=1 to see "clean" PSF of target at origin
% Try this with other targets

% xn: range;            yn= cross-range;    fn: reflectivity
  xn=zeros(1,ntarget);  yn=xn;              fn=xn;

% Targets within digital spotlight filter
%
  xn(1)=0;              yn(1)=0;            fn(1)=1;
  xn(2)=.7*X0;          yn(2)=-.6*Y0;       fn(2)=1.4;
  xn(3)=0;              yn(3)=-.85*Y0;      fn(3)=.8;
  xn(4)=-.5*X0;         yn(4)=.75*Y0;       fn(4)=1.;
  xn(5)=-.5*X0+DX;      yn(5)=.75*Y0+DY;    fn(5)=1.;

% Targets outside digital spotlight filter  
% (Run the code with and without these targets)
%
  xn(6)=-1.2*X0;        yn(6)=.75*Y0;       fn(6)=1.;
  xn(7)=.5*X0;          yn(7)=1.25*Y0;      fn(7)=1.;
  xn(8)=1.1*X0;         yn(8)=-1.1*Y0;      fn(8)=1.;
  xn(9)=-1.2*X0;        yn(9)=-1.75*Y0;     fn(9)=1.;
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                   SIMULATION                    %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
s=zeros(n,m);        % SAR signal array
%
for i=1:ntarget;
 td=t(:)*ones(1,m)-2*ones(n,1)*sqrt((Xc+xn(i))^2+(yn(i)-u).^2)/c;
 sn=fn(i)*exp(cj*wcm*td+cj*alpha*(td.^2)).*(td >= 0 & td <= Tp & ...
   ones(n,1)*abs(u) <= L & t(:)*ones(1,m) < Tf);

 sn=sn.*exp(-cj*wc*t(:)*ones(1,m));      % Fast-time baseband conversion????

%
% calculate (omega,u)-domain beam pattern
     theta_n=ones(nk,1)*atan((yn(i)-u)/(Xc+xn(i)));
     beam_p=zeros(n,m);
     beam_p(Ik,:)=(.5+.5*cos((pi*theta_n)./(phi_d(Ik).'*ones(1,m)))).* ...
        (abs(theta_n) <= phi_d(Ik).'*ones(1,m));
     
% incorporate beam pattern in (omega,u) domain
%
  s=s+iftx(ftx(sn).*beam_p);
end;
%

G=abs(s)';
xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
image(t,u,256-cg*(G-ng));
axis('square');axis('xy')
xlabel('Fast-time t, sec')
ylabel('Synthetic Aperture (Slow-time) U, meters')
title('Measured Stripmap SAR Signal')
print P6.1.ps
pause(1)
%
td0=t(:)-2*Xc/c;
s0=exp(cj*wcm*td0+cj*alpha*(td0.^2)).*(td0 >= 0 & td0 <= Tp);
s0=s0.*exp(-cj*wc*t(:));            % Baseband reference fast-time signal

s=ftx(s).*(conj(ftx(s0))*ones(1,m));  % Fast-time matched filtering
%
G=abs(iftx(s))';
xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
tm=(2*Xc/c)+dt*(-n/2:n/2-1);   % fast-time array after matched filtering
image(tm,u,256-cg*(G-ng));
axis('square');axis('xy')
xlabel('Fast-time t, sec')
ylabel('Synthetic Aperture (Slow-time) U, meters')
title('Stripmap SAR Signal after Fast-time Matched Filtering')
print P6.2.ps
pause(1)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Digital Spotlighting via Slow-time Compression  %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%% NOTE: This code performs "narrow-bandwidth" digital spotlighting.
%%       To improve results, convert the code to "wide-bandwidth"
%%       digital spotlighting (polar format processing).
%
% Upsample data before compression for digital spotlighting
%
fs=fty(s);
mz=mc-m;        % number of zeros is 2*mz
fs=(mc/m)*[zeros(n,mz/2),fs,zeros(n,mz/2)];
s=ifty(fs);

Ls=100*du;                 % Half-width of a subaperture (200 samples
                           % within each subaperture)
mcs=2*ceil(Ls/duc);        % Number of subaperture slow-time samples
                           % for compressed SAR signal
mcs=2^ceil(log(mcs)/log(2));    % Readjust to power of 2
if mcs > mc, mcs=mc; end;
ns=mc/mcs;                  % Number of subapertures
us=duc*(-mcs/2:mcs/2-1);    % Subaperture u array
dkus=pi2/(mcs*duc);         % Subaperture ku (Doppler) spacing
kus=dkus*(-mcs/2:mcs/2-1);  % Subaperture ku array
 

s_ds=zeros(n,mc);           % Initialize digital spotlight SAR array

PH=asin(kus/(2*kc));        % Angular Doppler domain
R=(c*tm)/2;                 % Transform fast-time to range

for i=1:ns; i                 % loop for each subaperture
    I=(1:mcs)+(i-1)*mcs;
    Yi=uc(I(mcs/2+1));        % Squint cross-range of i-th subaperture
    tei=atan(Yi/Xc);          % Squint angle of i-th subaperture
    Ri=sqrt(Xc^2+Yi^2);       % Squint radial range of i-th subaperture
    tt=(2*sqrt(Xc^2+uc(I).^2)-2*Ri)/c;
    fpi=iftx(fty(s(:,I).*exp(cj*w(:)*tt)));    % slow-time compression and
                                   % Fourier transform to (t,ku) domain
                                   % to obtain polar format reconstruction
                                   % within i-th subaperture
    W_di=(abs(R(:)*cos(PH+tei)-Xc) < X0 & ...
       abs(R(:)*sin(PH+tei)-Yi) < Y0);         % Digital spotlight filter
    G=(abs(fpi)/max(max(abs(fpi)))+.1*W_di)';
    xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
    image(tm,kuc,256-cg*(G-ng));
    xlabel('Fast-time t, sec')
    ylabel('Synthetic Aperture (Slow-time) Doppler Frequency k_u, rad/m')
    title('Polar Format SAR Reconstruction with Digital Spotlight Filter')
    axis square; axis xy;
    print P6.3.ps
    pause(1)
    s_ds(:,I)=ifty(ftx(fpi.*W_di)).*exp(-cj*w(:)*tt); % decompression
                          % and Fourier transforming to (omega,u) domain
end;
 
% Downsample data in ku domain to remove redundancy
% (Remove zeros which were added before slow-time compression)
%
fs=fty(s_ds);              % F.T. of SAR signal w.r.t. u
%
mz=mc-m;        % number is redundant samples is mz
fs=(m/mc)*fs(:,mz/2+1:mz/2+m);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%    2D BEAM PATTERN MATCHED FILTERING     %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                              CAUTION                            %
% Applying 2D beam pattern matched filter worsens the PSF. For    %
% digital reconstruction via spatial frequency interpolation or   %
% range stacking, one may remove this beam pattern matched filter %
% to improve PSF. However, for TDC and backprojection, removal of %
% this matched filter degrades the reconstruction in certain      %
% frequency bands.                                                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
phi=zeros(n,m);
Iku=(ones(nk,1)*abs(ku) < 2*k(Ik).'*ones(1,m));
phi(Ik,:)=asin((Iku.*(ones(nk,1)*ku))./(2*k(Ik).'*ones(1,m)));
beam_p=zeros(n,m);
beam_p(Ik,:)=(.5+.5*cos((pi*phi(Ik,:))./(phi_d(Ik).'*ones(1,m)))).* ...
        Iku.*(abs(phi(Ik,:)) <= phi_d(Ik).'*ones(1,m));
     
fs=fs.*beam_p;       % Beam pattern 2D matched filtering. To remove
                     % this filter, use: beam_p=ones(n,m)     

s_ds=ifty(fs);       % Save digitally-spotlighted (omega,u) domain
                     % data for TDC and projection

G=abs(iftx(s_ds))';
xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
image(tm,u,256-cg*(G-ng));
axis('square');axis('xy')
xlabel('Fast-time t, sec')
ylabel('Synthetic Aperture (Slow-time) U, meters')
title('Stripmap SAR Signal after Digital Spotlighting')
print P6.4.ps
pause(1)