
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                %
%   Phased Array Imaging of Stationary Targets   %
%              Using Beam Steering               %
%                                                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

colormap(gray(256))

cj=sqrt(-1); pi2=2*pi;

c=155000;                 % Propagation speed, cm/sec
fc=3e6;                   % Carrier frequency, Hz
f0=1e6;                   % [fc-f0,fc+f0] is the bandwidth
kc=(pi2*fc)/c;            % Wavenumber at the carrier frequency

x0=4.8;                   % Length of the target area in range is 2 x0
dkx=pi/x0;                % Sample spacing in the target k_x domain
dk=dkx/2;                 % Sample spacing in the wavenumber domain
df=(c*dk)/pi2;            % Sample spacing in the frequency domain
n=ceil((2*f0)/df);        % Number of frequencies used
k=kc+(-n/2:n/2-1)*dk;     % Wavenumber array
lambda=pi2/k(n);          % Wavelength at the highest frequency

dx=pi2/(dkx*n);           % Sample spacing in the target x domain
x=dx*(-n/2:n/2-1);        % Spatial domain range (x) array

r0=1.4*x0;                % Distance of the center of the array from the
                          % center of the target region; in the book "Fourier
                          % Array Imaging," r0 is something else (the depth of
                          % reconstruction)

y=x;                      % Spatial domain cross-range (y) array
y0=-y(1);                 % Length of the target area in cross-range is 2 y0
dky=pi/y0;                % Sample spacing in the target k_y domain
dy=pi2/(dky*n);           % Sample spacing in the target y domain

L=.75*(r0*lambda)/(4*dy); % Length of aperture is 2L. The factor .75
                          % can be removed and/or adjusted to yield the
                          % desired aperture length. However, this factor
                          % should always be less than "one" to avoid aliasing

du=(r0*lambda)/(4*(L+y0));  % Sample spacing in the aperture domain
m=ceil((2*L)/du);           % Number of elements on phased array; it can be
                            % readjusted to have an "odd" number of elements
                            % such that the array is symmeteric; this is done
                            % in the following (NOT IMPORTANT):
if floor(m/2)*2 == m, m=m+1; end; m2=floor(m/2); du=(2*L)/(m-1);

u=du*(-m2:m2);             % Aperture array

dte=lambda/(4*L);          % sample spacing in the beam steering domain
te0=asin(y0/r0);           % Beam steering is done in the interval [-te0,te0]
M=ceil((2*te0)/dte);       % Number of steer angles
                           % Readjusting it to become "odd" (NOT IMPORTANT)
if floor(M/2)*2 == M, M=M+1; end; M2=floor(M/2); dte=(2*te0)/(M-1);

te=dte*(-M2:M2);           % Beam steer array

%%%%%   Simulated targets' parameters

io=4;
xx=zeros(1,io); yy=xx;
xx(1)=68*dx; yy(1)=88*dy;
xx(2)=60*dx; yy(2)=80*dy;
xx(3)=68*dx; yy(3)=80*dy;
xx(4)=60*dx; yy(4)=88*dy;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                            %
%        SIMULATION          %
%                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

s=zeros(n,M);
for ii=1:io;
 ii
 d0=sqrt((r0+xx(ii))^2+(yy(ii)+u(:)*ones(1,M)).^2)+u(:)*sin(te);
      % d0 contains two components; the first one, i.e., "sqrt..."
      % is the target's distance from an element of the phased
      % array; the second component, i.e., "u*sin(te)", is the
      % delay for beam steering at angle "te" for the element at "u"
%
 for i=1:n;
  s(i,:)=s(i,:)+sum(exp(cj*k(i)*d0)).^2;
 end;
end;




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                            %
%      RECONSTRUCTION        %
%                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     Interpolation form "te" domain into "k_u=k_y=2k sin(te)" domain   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dKU=pi/L; KU0=2*k(n)*((L+y0)/r0);
MU=ceil((2*KU0)/dKU);
  if floor(MU/2)*2 ~= MU, MU=MU+1; end;
 KU=dKU*(-MU/2:MU/2-1);

ffs=polte(-KU0,dKU,MU,te,2*k,s);

ny=(2*KU0)/dky; ny=2^(ceil(log(ny)/log(2)));
dky=(2*KU0)/ny;
ky=dky*(-ny/2:ny/2-1);
dY=pi2/(ny*dky); Y=dY*(-ny/2:ny/2-1);

if ny > MU,
 nz=(ny-MU)/2;
 ff=ifty([zeros(n,nz),fty(ffs),zeros(n,nz)]);
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    Addition of the phase function;   see page 230,232   %
%    of the book "Fourier Array Imaging"                  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
KX=4*(k(:).^2)*ones(1,ny)-ones(n,1)*(ky.^2);
KX=sqrt(KX.*(KX > 0)); 
F=ff.*exp(-cj*KX*r0);     % r0 is the center of the target
                          % region (not depth as in the book)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     Interpolation form "k" domain into "k_x" domain   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ff=polk(2*k(1),dkx,n,2*k,ky,F);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Display reconstruction in the spatial frequency domain

kx=2*k(1)+dkx*(0:n-1);
 G=abs(ff)';img;
 image(kx,-ky,256-cg*(G-ng))
 axis('equal')
 title('Spatial Frequency Domain Reconstruction')
 ylabel('Cross-range Spatial Frequency Ky, radians/cm')
 xlabel('Range Spatial Frequency Kx, radians/cm')

%  pause


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    Inverse Fourier transform into spatial domain   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f=ftx(ifty(ff));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%  DISPLAY   %%%%%%%%%%%%%%%%%%

 G=abs(f)';img;
 image(r0+x,-Y,256-cg*(G-ng))
 axis('equal')
 title('Spatial Domain Reconstruction')
 ylabel('Cross-range, cm')                                   
 xlabel('Range, cm')                                         

