The first problem we would like to attack is the case when we have an
array of antennas randomly spread out in an area of about 1 to 3 acres.
In this case we want to use the array to find the azimuth and elevation
of incoming signals in the HF band.

What we need is an algorithm that can quickly construct a correlation
table for the expected phase difference between the 3 to 10 different
antennas (they are mounted on vehicles, one per vehicle). This table
has to be usable from 1 to 30 MHz and for angles up to 85 degrees
elevation. All the vehicles have GPS so the relative positions are known
to within a couple of feet.  The receive equipment can measure the
phase difference fast between all the different antennas.

With help of the correlation table the system can then look up the best
fit and output the result (which is of course a line of bearing as well
as the elevation of the incoming signal). If any vehicle moves, then a
new table has to be constructed.  So the algorithm must be able to take
in new vehicle locations, and continually calculate the new correlation

 - Lars

Hello Lars,

It never occurred to me that you could determine the direction,
including elevation, of an incoming wave by observing a set of phase
relationships among the received signals at several omnidirectional
antennas.  I've enjoyed thinking about this problem, which is a new one
for me.

Since it is somewhat awkward to express math and figures via email,
I've scanned in a couple sheets of notes I made while thinking about
this problem:
As I understand the problem, the core issue is to compute the expected
phase relationship between two antennas, given that a signal comes from
a given direction.  Then you want to access a table of such data, for several
antennas, to find a direction given the measured phase relationships.

I assume that the signal is sufficiently far away for the incoming
wavefront to be thought of as purely planar, with no significant
curvature.  The motion of the wavefront would then be perpendicular to
that plane.  This in turn implies that all the azimuth/elevation
headings will be the same, at each of our receivers.

For a given pair of antennas A and B with known positions, let d be
the distance between them.  Let us express the heading to the
wavefront using azimuth and elevation from antenna A.  The azimuth
will be measured from this axis.  The elevation angle from the ground
plane.  Call these angles "a" and "e" respectively.

In the notes, I worked out the expression for the quantity "L", which
is the distance the wavefront must travel between when it encounters
one antenna (B say), until it encounters the other antenna (A say).
Using these local coordinates, the expression for L is simply:

   L = d cos(a) cos(e)

And then the phase angle between A and B is just

   (L/lambda) * 360 degrees,

where lambda is the wavelength of the incoming signal.  (Omitting sign

So it would be a simple matter to compute the table that you seek, as I understand it.
You can just take the position data and desired heading, convert to the above-sketched
local coordinates, do the simple calculation, and fill up your table.

Of course I might have muffed the calculation.  If so, tell me!  (I don't mind having
my early mistakes caught so that I don't go too far with bad results.)


I would like to know the desired accuracy of the system.  Then I might be able to work out
how big these tables need to be.  Also, note that these tables I described depend
on the wavelength.  So if you have 50 transmitters on different freqs, you'll need
50 tables.

Note also that the stated gps accuracy of a couple feet might start to become
problematical up around 30 MHz, where the wavelength is only 10 meters, or around
30 feet.

Also, I'd really like to know how accurate the phase measurements are,
and what the resolution is.


Instead of using these tables and then looking for a good fit to solve the problem,
it would be better if we could just take the phase measurements, wavelenth, and
position data, and then just compute the heading to the wavefront.  It might turn
out to be feasible, and I have just begun to think about this.

When you do the problem in that order, I think you have to come to terms
with some issues that you get to ignore with the table approach.  One such
issue is that these randomly placed antennas will often be more than a wavelength
away from each other.  But you know where they are, so that should be okay.

In general, I think what you'll end up with, at the core, is a set of equations
that look something like, for example

   cos(a) cos(b) = 0.345

   cos(a) cos(b) = 0.567

where the two equations are using different local coordinates.  So you need
to unify them, and then determine  an "a" and "b" that make both true.  Unless
you're unlucky, there will be only one solution, and that is your azimuth
and elevation.  That would be your answer.  (In these equations, the phase
phase angle measurment and distance at one pair goes into the 0.345, and
a different pair goes into the 0.567.)

This is just a sketch, and I need to think about it more.  But it
might be quite feasible to do it that way, and if we can, that would
be good, wouldn't it?


So that's where I am.  Any feedback would be appreciated.


Hi Lars,

I worked out how to find an exact unique solution for (az, el) in a large
variety of cases, based on phase measurements between antennas.

Recall from the previous email that a central relationship is

  L = d cos(a) cos(e),

where d is the distance between the antennas, "a" azimuth, and "e" elevation.
L is the distance the wavefront travels between the time it encounters
the two antennas.  The expected phase measurement is then

  th = (L / lambda) * 360 degrees,

where lambda is the wavelength of the signal of interest.

Now, if d is not too large, you can reverse the computation and compute
L based on the phase measurement.  I assume that you can measure phase
accurately, but only between -180 and 180.  You can measure d accurately,
and th accurately.  But you really cannot tell if one wave is, say, 1050 degrees
ahead of another.  The reading would say -30 degrees.

So for now, let us assume that d is not too large, and in particular that
it is not more than half a wavelength of interest.  Then we can easily
determine L from a pure phase measurement.  The quantity L is really
what's of interest.

So given 

  L = d cos(a) cos(e),

we get 

  cos(a) cos(e) = (L/d) =: k.

I think of k as the "shrinkage factor".  It is a number between zero and
one that indicates how much shorter than the "d" distance the wave actually

The above equation expresses the relationship between azimuth and elevation
that must hold in order to be consistent with the phase measurement that
was made.


Here's an example.  The distance between the antennas, d, is 10 meters,
and we measure the measured phase angle between the signals at those antennas
to be 60 degrees, and the frequency is 10 MHz, or 30 meters.

Since d is less than a half wavelength, we easily solve

  60 = (L / 30) * 360,

getting L = 5 meters.  In turn, the ratio L/d is 5/10, or 0.5.  So the set of
feasible azimuth elevation pairs is given by

  cos(a) cos(e) = 1/2.

Here is a plot of the set of feasible pairs:

That general shape is always what such a solution set looks like.  Here is a set of
feasibilities for various shrinkage factors between 0 and 1:

Now suppose we make another phase measurement between a different pair
of antennas, also not too far away.  Really, the only other
requirement is that the new axis not be parallel to the original
antenna axis.  One antenna can even be shared between the two pairs.

We come up with an additional equation for the new antenna, and the azimuth/elevation
pairs are now constrained to be

  cos(a) cos(e) = k2,

where k2 is a potentially different shrinkage factor.  It is important to note,
however, that "a" and "e" are in coordinates local to the new axis.  We want
to express these azimuth/elevation pairs in the same coordinate system
as the first ones, and this is easy to do.  Recall that the actual direction
to the wavefront is the same as measured from any antenna, because the source
is far away.  Everybody reads the actual same elevation; it's only their
azimuths that are different in the local coordinate systems.  So all we have
to do to convert one coordinate system to another is to adjust the angle
of the azimuths by the angle between the antenna pairs.  Or alternately,
we can convert each to global az/el by adjusting each az according to the
axis deviation from north.

So what happens is that we end up with two equations like this:

  cos(a) cos(e) = k1
  cos(a+phi) cos(e) = k2

where phi is a constant we know, and k1 and k2 are each shrinkage
factors.  Each of those is a curve like we saw above (the second one
is shifted on the a-axis), and it fairly clearly intersects at a point
(or all the points).  It cannot intersect at two points, for example.
So that point is the only az/el pair that is consistent with our phase
measurements, and is the answer.


For example, suppose we have

  cos(a) cos(e) = 0.70
  cos(a - 30) cos(e) = 0.53

Here they are plotted together:
We see graphically that they intersect at one point.


It is easy and efficient to actually compute that intersection point.  All
you do is expand cos(a+ph) out according to the well-known relation

  cos(a+phi) = cos(a)cos(phi) - sin(a)sin(phi),

noting that cos(phi) and sin(phi) are constants.  When you simplify
appropriately, you get

  tan(a) = (cos(phi) - (k2/k1)) / sin(phi)

So an appropriate arctan of that is "a".  Given "a", it is easy to get

  cos(e) = k1 / cos(a),

so "e" is an appropriate arccos of that.  (It turns out that the "principle value"
arctan and arccos are the ones you want).


These are really efficient calculations.  Just a few mults and divides, a few
cosines, an arccos, and an arctan.  Even slow hardware should be able to do 
many thousands of these per second.


I've attached a simple
Perl script
that works through some example calculations.  It is also the code
I used to generate the plots.  I used Perl because it was expedient
for investigation.  Of course, I can do the core calculation very
efficiently in C.  It's probably workable even if you don't have
floats (but I assume you do these days).

The code plops down a few antennas and specifies a frequency.  Then
you specify an actual (az, el) to the oncoming wavefront.  The code
computes the expected phase measurements, and then given only this
computes the intersection of the two curves to re-derive the specified
direction to the wavefront.  Of course, that part of the code doesn't
"look at" the answer.

So this gives some degree of confidence.  (Note that it does depend, critically,
on the validity of L = d cos(a) cos(e), however.)

One can use this as the beginnings of a robust investigation of precision
issues.  For example, we could perturb the actual phase measurments, adding
errors that we believe model what we would get in practice, and then observe
what happens with different oncoming wavefronts.  We could vary the antenna
placement, frequency, etc.  This could be a useful thing to do, even if we
also analyze the error in more analytical ways.

Real data would be nice, too.


So consider, the following simplified algorithm could work much of the time.
You have a field full of antennas, and a signal of interest.  Pick a pair of antennas
that are within a half wavelength, and measure the phase.  Then pick another
pair, and measure that phase. (Or just a third antenna not in line and also not
too far away).  Do the simple calculation, and get the answer.

Interestingly, you can then use all your other antennas, even ones that are
"too far" from each other, to reinforce the answer you got.  You just
compute the expected phase between an *arbitrary* pair of antennas, and
verify that's what you measured.

But notice that the core result needs only three, or maybe four antennas.


If you cannot find appropriately placed antennas (for the frequency of interest),
that means that either they are all too far from each other, or they are all lined
up in a row.  If they are all in a line, you lose--move them!

If the problem is, say, that you're trying to measure for 30 MHz or 10 meters,
and you don't have any antennas within 5 meters of each other (or only one
such pair), then the problem is a little more involved.  But would this
happen in practice?

If it does, I have ideas on how to proceed.  It would be pretty much
like the above, but the difference would relate to the th-to-L
calculation.  You don't get a unique value for L, but you get several
possible values, the exact number of which depends on lamdbda and d.

I think it's just a matter of looking at a fairly small number of discrete
cases and figuring out what is going on.  There might also be resolution issues that
are harder in this case.


So given my assumptions and understanding of the problem, I think this line
is potentially winning.  Please let me know what you think.