Visualize MIMO radar antenna array radiation pattern.

radar_mimo.py

#!/usr/bin/python3
"""
Visualize MIMO radar antenna array radiation pattern.
"""
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(123)
plt.style.use('ggplot')

# Target angle
phasing = [0, 0]
# 2D plot threshold
threshold = 20
# Number of points in the plot
plot_points = 360

# Antennas. Change 'if 0' to 'if 1' to select the antenna.
if 0:
    # 2D split TX
    Nrx = 8
    Ntx = 8

    de = 0.5

    drx = np.array([de, 0])
    dtx = np.array([Nrx * de / 2, 0])

    tx_rx_offset = np.array([Nrx * de / 2, 1])

    p_rx = np.array([drx * i for i in range(Nrx//2)])
    rx2_offset = np.array([de * Ntx * Nrx/2, 0])
    p_rx2 = np.array([drx * i + rx2_offset for i in range(Nrx//2)])
    p_rx = np.append(p_rx, p_rx2, axis=0)
    p_tx = np.array([dtx * i + tx_rx_offset for i in range(Ntx)])
elif 0:
    # 2D twice-split TX
    Nrx = 8
    Ntx = 8

    de = 0.5

    drx = np.array([de, 0])
    dtx = np.array([Nrx * de / 4, 0])

    tx_rx_offset = np.array([Nrx * de / 4, 1])

    p_rx = np.array([drx * i for i in range(Nrx//4)])
    rx2_offset = np.array([de * Ntx * Nrx/8, 0])
    p_rx2 = np.array([drx * i + rx2_offset for i in range(Nrx//4)])
    p_rx = np.append(p_rx, p_rx2, axis=0)
    p_tx = np.array([dtx * i + tx_rx_offset for i in range(Ntx//2)])

    p_rx3 = p_rx + np.array([(Ntx//2) * (Nrx//2+4) * de, 0])
    p_tx3 = p_tx + np.array([(Nrx//2) * (Ntx//2+0) * de, 0])
    p_rx = np.append(p_rx, p_rx3, axis=0)
    p_tx = np.append(p_tx, p_tx3, axis=0)
elif 0:
    # 2D
    Nrx = 8
    Ntx = 8

    de = 0.5

    drx = np.array([de, 0])
    dtx = np.array([Nrx * de, 0])

    tx_rx_offset = np.array([0, 1])

    p_rx = np.array([drx * i for i in range(Nrx)])
    p_tx = np.array([dtx * i + tx_rx_offset for i in range(Ntx)])
elif 0:
    # Dense 3D
    Nrx = 64
    Ntx = 1

    de = 0.5

    assert int(Nrx**0.5)**2 == Nrx
    p_tx = np.array([[1.75, 1.75]])

    p_rx = np.array([[i*de, j*de] for i in range(int(Nrx**0.5)) for j in range(int(Nrx**0.5))])
elif 0:
    # Linear 3D
    Nrx = 8
    Ntx = 8

    de = 0.5

    drx = np.array([de, 0])
    dtx = np.array([0, de])

    tx_rx_offset = np.array([-de/2, de/2])

    p_rx = np.array([drx * i for i in range(Nrx)])
    p_tx = np.array([dtx * i + tx_rx_offset for i in range(Ntx)])
elif 0:
    # Box 3D
    Nrx = 8
    Ntx = 8

    drx = np.array([0.5, 0])
    dtx = np.array([0, 0.5])

    offset = np.array([-0.5/2, 0.5/2])

    p_rx = np.array([drx * i for i in range(Nrx//2)])
    p_tx = np.array([dtx * i + offset for i in range(Ntx//2)])

    p_rx2 = p_rx + np.array([0, Ntx * 0.5/2])
    p_tx2 = p_tx + np.array([Nrx * 0.5/2, 0])

    p_rx = np.append(p_rx, p_rx2, axis=0)
    p_tx = np.append(p_tx, p_tx2, axis=0)
elif 1:
    # Tiled Box 3D
    Nrx = 8
    Ntx = 8

    de = 0.5

    drx = np.array([de, 0])
    dtx = np.array([0, de])

    tx_rx_offset = np.array([-de/2, de/2])

    p_rx = np.array([drx * i for i in range(Nrx)])
    p_tx = np.array([dtx * i + tx_rx_offset for i in range(Ntx)])

    p_rx2 = p_rx + np.array([0, Nrx*de])
    p_tx2 = p_tx + np.array([Ntx*de, 0])

    p_rx = np.append(p_rx, p_rx2, axis=0)
    p_tx = np.append(p_tx, p_tx2, axis=0)

    p_rx3 = p_rx + np.array([0, 2 * Nrx*de])
    p_tx3 = p_tx + np.array([0, 2 * Ntx*de])

    p_rx = np.append(p_rx, p_rx3, axis=0)
    p_tx = np.append(p_tx, p_tx3, axis=0)

    p_rx4 = p_rx + np.array([2 * Nrx*de, 0])
    p_tx4 = p_tx + np.array([2 * Ntx*de, 0])

    p_rx = np.append(p_rx, p_rx4, axis=0)
    p_tx = np.append(p_tx, p_tx4, axis=0)
else:
    Nrx = 8
    Ntx = 8
    derx = 0.25
    detx = 0.25
    p_rx = np.random.uniform(-Nrx * derx, Nrx * derx, [Nrx, 2])
    p_tx = np.random.uniform(-Ntx * detx, Ntx * detx, [Ntx, 2])

# Plot antenna array
if 1:
    plt.figure(figsize=(5, 5))
    plt.title('{}TX-{}RX element positions'.format(Ntx, Nrx))
    plt.xlabel('Distance (Wavelenghts)')
    plt.ylabel('Distance (Wavelenghts)')

    plt.scatter(p_rx.T[0], p_rx.T[1], marker="+", label='RX')
    plt.scatter(p_tx.T[0], p_tx.T[1], marker="x", label='TX')
    p_virtual = np.array([0.5*(r+t) for r in p_rx for t in p_tx])
    v_spacings = np.diff(sorted(p_virtual[:,0]))
    if len(set(v_spacings)) == 1:
        print('Equal virtual element spacing', set(v_spacings))
    else:
        print('Un-equal virtual element spacing', set(v_spacings))
    plt.scatter(p_virtual.T[0], p_virtual.T[1], s=25, marker=".", label='Virtual')
    x0 = np.min([np.min(p_rx[:,0]), np.min(p_tx[:,0])])
    x1 = np.max([np.max(p_rx[:,0]), np.max(p_tx[:,0])])
    y0 = np.min([np.min(p_rx[:,1]), np.min(p_tx[:,1])])
    y1 = np.max([np.max(p_rx[:,1]), np.max(p_tx[:,1])])
    print('X-size {} lambda, Y-size {} lambda'.format(x1 - x0, y1 - y0))

    vx0 = np.min(p_virtual[:,0])
    vx1 = np.max(p_virtual[:,0])
    vy0 = np.min(p_virtual[:,1])
    vy1 = np.max(p_virtual[:,1])
    print('Virtual array X-size {} lambda, Y-size {} lambda'.format(vx1 - vx0, vy1 - vy0))
    plt.legend(loc='best')
    plt.savefig('virtual.png')

# Plot points
d1d = np.linspace(-np.pi/2, np.pi/2, plot_points)
# Cartesian product to make 2D grid
d = np.array([(i, j) for i in d1d for j in d1d])

phasing = np.array(phasing) * np.pi/180

# Coordinate system
if 1:
    # Spherical
    r0 = np.sin(d[:,0]) * np.cos(d[:, 1])
    r1 = np.sin(d[:, 1])
    phasing = np.array([np.sin(phasing[0]) * np.cos(phasing[1]), np.sin(phasing[1])])
elif 0:
    # Whatever this is called. 1D in both directions
    r0 = np.sin(d[:,0])
    r1 = np.sin(d[:, 1])
    phasing = np.array([np.sin(phasing[0]), np.sin(phasing[1])])

if 1:
    # Decompose to TX and RX patterns. Much quicker to calculate
    tx_af = np.sum([np.exp(1j*2*np.pi*(-(r0 * p[0] + r1 * p[1]) + (phasing[0] * p[0] + phasing[1] * p[1]))) for p in p_tx], axis=0).reshape([plot_points, plot_points]).T
    rx_af = np.sum([np.exp(1j*2*np.pi*(-(r0 * p[0] + r1 * p[1]) + (phasing[0] * p[0] + phasing[1] * p[1]))) for p in p_rx], axis=0).reshape([plot_points, plot_points]).T

    txrx_af = np.abs(tx_af) * np.abs(rx_af)
else:
    p_txrx = [ptx - prx for prx in p_rx for ptx in p_tx]
    txrx_af = np.sum([np.exp(1j*2*np.pi*(-(r0 * p[0] + r1 * p[1]) + (phasing[0] * p[0] + phasing[1] * p[1]))) for p in p_txrx], axis=0).reshape([plot_points, plot_points]).T

txrx_af = 20 * np.log10(np.abs(txrx_af))
txrx_af -= np.max(txrx_af)

plt.figure(figsize=(7, 5))
plt.title('Target response')
plt.plot((180/np.pi)*d1d, txrx_af[txrx_af.shape[0]//2, :])
plt.xlabel('Angle (deg)')
plt.ylabel('Magnitude (dB)')
plt.ylim([-40 ,0])
plt.savefig('target_response_1d.png')

plt.figure(figsize=(5, 5))
plt.grid(False)
plt.title('2D target response')
imgplot = plt.imshow(txrx_af, extent=[-90, 90, -90, 90], origin='lower')
cbar = plt.colorbar(imgplot, orientation='vertical')
cbar.ax.set_ylabel('J [dB]')
m = np.max(txrx_af)
plt.xlabel('Angle (deg)')
plt.ylabel('Angle (deg)')
imgplot.set_clim(m - threshold, m)
plt.savefig('target_response_2d.png')
plt.show()