I was looking at doing some experiments with non-uniform arrays for direction finding or triangulation. I was thinking about some simple ways to compute the direction considering all elements.
I came to two major methods. The first examines the wave as planar and the second method looks at the wave as spherical with some defined point of origin. I then looked at combining the two methods together.
Now, with two elements a wave both planar and spherical arriving has infinite solutions for the possible origin. This is because the solution space, in 3D and not 2D, of possible origins exists between the forms of a line, cone, and plane.
It is a line with the wave arrives along the axis of the two elements, a plane with it arrives perpendicular to the axis of the two elements, and a cone at all other times.
To examine this visually think of a 3D point and the angle a wave coming from it, planar or spherical, would arrive onto the two elements. It would create a phase difference as a function of its wavelength and with this phase difference one can compute the angle of incidence.
Now, rotate this point around the two-element axis and you should note how the angle of incidence stays the same, but the point/source changes location. This is the cone. You can also move the point further away while maintaining the same angle of incidence and this should further illustrate the shape of a cone.
With a non-uniform multi-element array the solution for the source when the wave is spherical is the intersection of the shapes. In some cases, it would be possible to have multiple solutions or multiple points of intersection.
The algorithm I have been working with is meant to compute these intersections in reasonable time and it is meant to do so for both a planar wave which has no source and a spherical wave which does.
I figured it would be nice to have an algorithm that could find both because in some cases the wave is planar simply because of the numerical precision. In other words, the angle of incidence is either too small to be measurable or because of noise it ends up having no definite source point that is computable. Therefore, having an algorithm that could handle this kind of situation was desirable.
I want to discuss the two algorithms in parts. I want to explain a method to find the direction of a planar wave that produces a 3D vector and then look at a method that finds the origin and it too produces a 3D vector.
Planar Wave
The algorithm is very simple really. You start with a 3D discrete uniform grid with its points all relative to the antenna/elements. Then you make each point a vector of magnitude one that points in the direction of the planar wave for its angle of incident. You create one grid per pairs of elements.
At this point you have a rough approximation of the cone, line, or plane shapes except they are for a planar wave which means they just simply rotated about the element pair axis.
Then you compute the squared error for each point across all grids. Finally, you find the point with the least error. There may be many of them with the same error. Finally, sum the vectors for the least error point for the final vector.
Spherical Wave
Instead of vectors you use scalar magnitudes, and you use a magnitude of one along the shape and then you create a gradient that diminishes to zero as the angle moves away from the cone’s surface. If you are dealing with a plane, then it diminishes to zero as you move toward the ends and for a line it reaches zero at 180-degrees.
Once again you can compute the distance for each point across the grids and look for the point with the highest distance. This point represents the most likely intersection.
Results
The result is that the planar wave idea failed to work because there were cases that depending on the configuration of the array and the source location it would produce completely incorrect approximations of the plane wave direction.
The spherical wave idea worked, and it also seems to approximate the planar wave direction. However, the algorithm can be quite slow because it must be computed over the grid per grid point. Therefore, the run-time increases in proportion to the number of grid points. Yet, it can be computed over a low-resolution grid and while that introduces considerable error it does produce reasonable outputs which means depending on the usage it might be a viable method.
At the end is the code to produce the two different grids. I only demonstrated the spherical. I call it intercept in the code under the function nua_compute_intersection_grid_based. The nua_compute_plane_fronts computes the planar vectors but as I said I wasn’t able to find a valid working way to use those vectors to produce something useful in the general case.
You can edit the test_point variable to produce the signal at different places in space. The other parameters are also editable such as the antenna locations. I have the antennas setup as a four-element array in a 3D configuration.
import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt
import pickle
import scipy.fft as fft
class AntSim:
"""Simulates discrete signal propogation.
"""
def __init__(self, sps, samps):
self.wave_velocity = 299792458.0
self.sps = sps
self.samps = samps
self.ant_locs = []
self.ant_sig = []
def add_ant(self, x, y, z):
self.ant_locs.append(np.array([x, y, z]).astype(np.float64))
self.ant_sig.append(np.zeros(self.samps, np.complex128))
def do_signal(self, signal_pos, center_freq, signal):
#signal_pos = np.array([x, y, z]).astype(np.float64)
assert hasattr(signal_pos, 'dtype')
sample_dist = self.wave_velocity / self.sps
for ndx, ant_loc in enumerate(self.ant_locs):
pos_vec = signal_pos - ant_loc
print('ant_loc', ant_loc)
print('signal_pos', signal_pos)
dist = np.sqrt(np.sum(pos_vec ** 2))
sample_delay = dist / sample_dist
sample_whole_delay = int(sample_delay)
sample_frac_delay = sample_delay - sample_whole_delay
# How much room is left in the buffer. The rest of the
# signal just flies on past.
left = len(self.ant_sig[ndx]) - sample_whole_delay
# This is weird looking but I really feel like it
# follows the laws of physics. Each sinsusodial is
# going to arrive at a different phase.
new_signal = phase_shift_signal_by_distance(
self.sps, center_freq, signal, dist, self.wave_velocity
)
# The attenuation if radiated from an isotropic source.
attn = 1 / dist ** 2
if left > len(new_signal):
left = len(new_signal)
self.ant_sig[ndx][sample_whole_delay:sample_whole_delay+left] += new_signal[0:left] * attn
def phase_shift_signal_by_distance(sps, center_freq, signal, dist, wave_velocity):
"""Simulate the phase shifts a complex signal undergoes as it travels over a distance.
"""
fdom = fft.fft(signal)
fdom_freqs = fft.fftfreq(len(signal), 1 / sps)
phase_distance = dist / (wave_velocity / (fdom_freqs + center_freq)) * np.pi * 2
cycles = phase_distance // (np.pi * 2)
fdom *= np.exp(1j * phase_distance)
return fft.ifft(fdom)
def nua_compute_plane_fronts(phase_thetas, ant_locs, gdx, gdy, gdz, gus, grid_center):
"""Non-Uniform Array compute plane fronts.
Returns 3D grid for each pair of elements with possible plane fronts at each grid
point. Each point is a 3D vector representing possible direction of plane wave.
"""
grid_corner = np.array([
grid_center[0] - (gdx * gus) * 0.5,
grid_center[1] - (gdy * gus) * 0.5,
grid_center[2] - (gdz * gus) * 0.5,
])
grid = np.zeros((
len(phase_thetas) * (len(phase_thetas) - 1), gdx, gdy, gdz, 3
), np.float64)
q = -1
for a in range(0, len(phase_thetas)):
for b in range(a + 1, len(phase_thetas)):
q += 1
theta = phase_thetas[a] * np.conjugate(phase_thetas[b])
a_pos = ant_locs[a]
b_pos = ant_locs[b]
vec_ab = b_pos - a_pos
vec_ab_dist = vec_dist(vec_ab)
for x in range(grid.shape[1]):
for y in range(grid.shape[2]):
for z in range(grid.shape[3]):
g_pos = grid_corner + np.array([x, y, z])
hypto = g_pos - a_pos
hypto_dist = vec_dist(hypto)
if hypto_dist == 0:
grid[q, x, y, z, :] = 0
continue
adj = np.dot(hypto / hypto_dist, vec_ab) * hypto_dist
# Find 90 degree corner of right triangle.
corner = vec_ab * adj + a_pos
# Get axis towards A from corner.
vec_ba = a_pos - corner
vec_ba /= vec_dist(vec_ba)
# Get axis from corner toward grid point.
vec_cg = g_pos - corner
vec_cg /= vec_dist(vec_cg)
# Rotate along the plane defined by vec_ba and vec_cg.
if norm_dot(corner - a_pos, vec_ab) < 0.98:
# If we go behind pos_a we need to reverse theta because
# the right triangle gets flipped.
vec_ca = a_pos - corner
vec_ca /= vec_dist(vec_ca)
the_planar_vec = vec_ca * np.sin(-theta) + vec_cg * np.cos(-theta)
else:
the_planar_vec = vec_ba * np.sin(theta) + vec_cg * np.cos(theta)
grid[q, x, y, z, :] = the_planar_vec
return grid
def nua_compute_intersection_grid_based(
phase_thetas,
ant_locs,
gdx, gdy, gdz,
gus,
grid_center,
freq,
wave_velocity,
test_point):
"""Non-Uniform Array computation of intersection using a grid based method.
Returns 3D grid containing error points. Lowest error points represent intersections
or near intersections of the conic shapes representing angle of incidences of planar
of spherical waves onto the array.
"""
grid_corner = np.array([
grid_center[0] - (gdx * gus) * 0.5,
grid_center[1] - (gdy * gus) * 0.5,
grid_center[2] - (gdz * gus) * 0.5,
])
grid = np.zeros((
len(phase_thetas) * (len(phase_thetas) - 1), gdx, gdy, gdz
), np.float64)
q = -1
for a in range(0, len(phase_thetas)):
for b in range(a + 1, len(phase_thetas)):
q += 1
theta = np.angle(phase_thetas[a] * np.conjugate(phase_thetas[b]))
a_pos = ant_locs[a]
b_pos = ant_locs[b]
vec_ab = b_pos - a_pos
vec_ab_dist = vec_dist(vec_ab)
vec_ab /= vec_ab_dist
wl = wave_velocity / freq
pl = vec_ab_dist / wl * np.pi * 2
# WRONG
# theta = theta / pl * np.pi * 0.5
# RIGHT
theta = np.arcsin(theta / pl)
theta = np.pi * 0.5 - theta
array_center = vec_ab_dist * 0.5 * vec_ab + a_pos
dbg_theta = np.arccos(norm_dot(test_point - array_center, b_pos - array_center))
assert abs(theta - dbg_theta) < 0.01
qq = -1
k = np.zeros((
grid.shape[1] * grid.shape[2] * grid.shape[3], 6
), np.float64)
for x in range(grid.shape[1]):
for y in range(grid.shape[2]):
for z in range(grid.shape[3]):
qq += 1
g_pos = grid_corner + np.array([x * gus, y * gus, z * gus])
vec_gac = array_center - g_pos
vec_gac_dist = vec_dist(vec_gac)
# toa cah soh
v_theta = np.arccos(norm_dot(b_pos - array_center, g_pos - array_center))
if v_theta < theta:
n_theta = theta - v_theta
else:
n_theta = v_theta - theta
grid[q, y, x, z] = np.abs(np.sin(n_theta) * vec_gac_dist)
if np.all(g_pos == test_point):
print('g_pos', g_pos)
print('g_xyz', x, y, z)
print('@theta', theta / np.pi * 180)
print('@v_theta', v_theta / np.pi * 180)
print('@n_theta', n_theta / np.pi * 180)
print('@grid', grid[q, y, x, z])
# end of xyz iteration
plt.title('View of Antenna %s And %s (View of XY plane [minimum across Z])' % (a, b))
plt.imshow(np.min(grid[q, :, :, :], axis=2))
plt.show()
return np.sqrt(np.sum(grid ** 2, axis=0))
def norm_dot(a, b):
ad = vec_dist(a)
if ad > 0:
a = a / ad
bd = vec_dist(b)
if bd > 0:
b = b / bd
return np.dot(a, b)
def vec_dist(v):
return np.sqrt(np.sum(v ** 2))
def main2():
sps = 24000
center_freq = 1000e6
dist = 0 #50000000
fbeg = -5000
fend = 5000
wave_velocity = 299792458.0
ant_sim = AntSim(sps, sps * 10)
#thetas = np.linspace(fbeg, fend, 20000) / (sps * 0.5) * np.pi * 2
#thetas = np.add.accumulate(thetas)
#chirp = np.exp(1j * thetas)
f = 1e3
t = 20000 / sps
chirp = np.exp(1j * np.linspace(0, np.pi * 2 * f * t, 20000))
# Four antennas. Three as a somewhat equlateral triangle.
# The forth behind them centered.
ant_sim.add_ant(-.05, 0, 0)
ant_sim.add_ant(0.05, 0, 0)
ant_sim.add_ant(0, 0.1, 0)
ant_sim.add_ant(0, 0.05, -0.1)
test_point = np.array([-2, 0, 0])
ant_sim.do_signal(test_point, center_freq, chirp)
ant_sig = np.array(ant_sim.ant_sig)
gd = 21
thetas = []
for a in range(ant_sig.shape[0]):
thetas.append(ant_sig[a, 0])
out = nua_compute_intersection_grid_based(
thetas, ant_sim.ant_locs, gd, gd, gd, 1.0, np.array([0., 0., 0.]),
center_freq + 1e3, wave_velocity, test_point
)
yi, xi, zi = np.unravel_index(np.nanargmin(out), out.shape)
plt.title('Error Across All Views (XY view)')
plt.imshow(np.nanmin(out, axis=2))
plt.show()
print('@', xi, yi, zi, out[yi, xi, zi])
if __name__ == '__main__':
main2()