Skip to content

Instantly share code, notes, and snippets.

@Ttl
Created May 15, 2021 05:26
Show Gist options
  • Save Ttl/3bb3e7c3213374f0dcb056b256795a76 to your computer and use it in GitHub Desktop.
Save Ttl/3bb3e7c3213374f0dcb056b256795a76 to your computer and use it in GitHub Desktop.
Visualize MIMO radar antenna array radiation pattern.
#!/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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment