Contents

%
%   Computational 2D SIMO imaging
%   Thomas Fromenteze
%

clear all
clc

Scene definition

% Parameters
nbf = 1001;                  % number of frequency samples
f = linspace(15e9,20e9,nbf); % frequency vector
c = 3e8;                     % speed of light
k = 2*pi*f/c;                % wavenumber
df = f(2)-f(1);              % frequency resolution
t = linspace(0,1/df,nbf);    % time vector

% Transmit antenna
xT = 0;
yT = -0.05;

% Receive array
dxR = c/f(end)/2;            % receive antenna array sampling
nbR = 51;                   % number of antennas
xR = (1:nbR)*dxR;            % receive array coordinates
xR = xR - mean(xR);          % centered on x = 0;
yR = xR .*0;                 % receive array on the plane y = 0;

% Parameter of the frequency-diverse 1D aperture
tau = 0.25e-7;         % decay time of the aperture
Q = 2*pi*mean(f)*tau;  % composite quality factor
hR = randn(nbR,nbf) .* repmat(exp(-t/(2*tau)),[nbR,1]);
HR = fft(hR,[],2);

% Target
xC = 0;
yC = 0.2;
SigC = 1; % Reflectivity

Propagation

% for loops are used to facilitate understanding. It's a very bad idea if you try to optimize computation times

S = zeros(nbR,numel(f));
for mC = 1:numel(xC)
    rT = sqrt((xT-xC(mC)).^2 + (yT-yC(mC)).^2); % distance from transmit antennas to target(s)
    rR = sqrt((xR-xC(mC)).^2 + (yR-yC(mC)).^2); % distance from targets to the receive antenna

    for mf = 1:numel(f)
        % interaction of the transmit waves with the target(s) and measurement
        S(:,mf) = S(:,mf) + (exp(-1j*k(mf)*rT)./sqrt(rT) .* SigC(mC) .* exp(-1j*k(mf)*rR)./sqrt(rR)).';
    end
end

Rho = sum(HR.*S,1); % Frequency signal on the output of the receive aperture
g = ifft(Rho);     % Time domain counterpart

% Reconstruction grid (should be defined according to the resolution)
x = linspace(-0.4,0.4,41);
y = linspace(0.1,1.1,41);

% Computation of the transfer (Green) matrix
[XR,F,X,Y] = ndgrid(single(xR),single(f),single(x),single(y));
G = exp(-1j*(2*pi*F/c).* ((sqrt((xT-X).^2 + (yT-Y).^2)+(sqrt((XR-X).^2 + (yR(1)-Y).^2)))));
% G = reshape(G,[numel(xR)*numel(f),numel(x)*numel(y)]);

HR_repmat = repmat(single(HR),[1,1,numel(x),numel(y)]);
H = squeeze(sum(HR_repmat.*G,1));
clearvars G HR_repmat XR F X Y

H = reshape(H,[nbf,numel(x)*numel(y)]);

Im = H'*Rho.';
Im = reshape(Im,[numel(x),numel(y)]);

figure(1), clf()
subplot(2,4,[1 2 5 6])
hold on
plot(xT,yT,'kv')
plot(xR,yR,'r.')
plot(xC,yC,'bo')
hold off
legend('Transmit antennas','Receive aperture','Target(s)')
grid on
xlabel('x (m)')
ylabel('y (m)')
daspect([1 1 1])
ylim([-0.4 0.4])
ylim([-0.1 1])

subplot(2,4,[3 4])
plot(t*1e9,real(g)/250)
xlabel('t (ns)')
title('g(t)')
ylim([-1 1])
set(gca,'YTickLabel',[]);

subplot(2,4,7)
pcolor(x,y,abs(Im).')
shading flat
xlabel('x (m)')
ylabel('y (m)')
daspect([1 1 1])

title('Reconstructed image')