%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                %
%  Ambiguity Function for Phased Array Imaging   %
%  with Randomized/Nonrandomized 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
dtau=(2*(r0+x0))/c;        % Sample spacing in the slow-time domain which
                           % is restricted by the round trip delay from the
                           % largest depth of the target region, i.e., r0+x0
tau=dtau*(-M2:M2);         % Nonrandomized slow-time array
TAU=tau;                   % This arry will contain the randomized version of
                           % the slow-time array if the flag "ir=1" is set
                           % NOTE: For computational purposes, it is more
                           % convenient to randomize the slow-time with
                           % respect to the beam steer angle. This yields the
                           % same mathematical model as randomizing the beam
                           % steer angle with respect to the slow-time.
                           % 
%%%%%%%% Slow-time randomizing

ir=1;                      % This flag is ON for randomized beam steering
if ir == 1,
 Z=rand(1,M);              % Uniform random generator
 [Z,I]=sort(Z);            % Sorted random data
                           % What follows is the generation of random slow
                           % time array
 for i=1:M;
  TAU(i)=tau(I(i));
 end;
end;

Dtheta=atan(L/r0);         % Beam steer angular interval over which a target
                           % at the origin, i.e., (r0,0), is "observable" is
                           % [-Dtheta,Dtheta].

Dtau=dtau*(Dtheta/dte);    % Slow-time interval over which a target at the
                           % origin, i.e., (r0,0), is "observable" is
                           % [-Dtau,Dtau].

i1=ceil(Dtheta/dte);       
II=M2+1-i1:M2+1+i1;        % II represents the indices in the beam steer
                           % array which represent the angular interval of
                           % the target observability, i.e., [-Dtheta,Dtheta].
                           % The ambiuity function integral is obtained only
                           % over these indices. NOTE: The ambiguity function
                           % looks better if II is chosen to include ALL beam
                           % steer angles, i.e., II=1:M. However, this
                           % implies that the target is obsevable at ALL
                           % beam steer angles which is not realistic. 

dv=c/(4*Dtau*fc);          % Speed resolution
nn=64;                     % Size of the ambiguity array is nn X nn.

aa=.2*dv*(0:nn-1); bb=aa;  % (aa,bb) are the velocities for which the ambiguity
                           % function is computed.

h=zeros(nn,nn);            % Initialize the ambiuity array

%%%%%%%%%%%%% Computing the ambiguity function

for i=1:nn;
 for j=1:nn;
  h(i,j)=sum(sum(exp(cj*2*k(:)*(TAU(II).*(aa(i)* ...
   cos(te(II))+bb(j)*sin(te(II)))))));
 end;
end

G=abs(h)';img; image(aa,bb(nn:-1:1),256-cg*(G-ng))
axis('square')             
ylabel('Perpendicular to Line of Sight Speed, cm/sec')
xlabel('Line of Sight Speed, cm/sec')                           
if ir == 1,
 title('Velocity Ambiguity Function: Randomized Beam Steering');
 else,
 title('Velocity Ambiguity Function: Conventional Beam Steering');
end;