%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% 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')