Some notes on the averaging of bearings.

A bearing is derived after each Fourier transform, for each
frequency bin.   At 192k samples/sec and a transform frame of
8192 samples, there are about 23 bearing estimates generated per
second, and it is necessary to average these over a period of
time (say, 120 seconds) to produce some kind of mean bearing
and also ideally a measure of the spread, such as the standard
deviation.

This is a non trivial matter, since the bearings wrap from 359
degrees around to zero.  Therefore, a naive mean will sometimes
produce incorrect answers.  For example, the average of two
raw readings (350, 10) if calculated as sum/count gives 180,
whereas we would require zero degrees (North) in this case.

Although this is a fairly well known problem, the many ad-hoc
programming solutions found by an internet search nearly all
turn out to be incorrect in some way.

There are two viable one-pass solutions.  In pseudo code,
let the set of N raw bearings be b[1..N].

The Mitsuta method
------------------
   sum = D = b[1];

   for( i from 2 to N)
   {
      delta = b[i] - D;

      if( delta < -180) D = D + delta + 360;
      else
      if( delta < 180) D = D + delta;
      else
         D = D + delta - 360;

      sum = sum + D;
   }
  
   mean = sum/N;

The direction cosine method
---------------------------
   sum_sin = 0;
   sum_cos = 0;

   for( i from 1 to N)
   {
      sum_sin = sum_sin + sin( b[i]);
      sum_cos = sum_cos + cos( b[i]);
   } 

   mean = atan2( sum_sin, sum_cos);

In the case of a 2-axis receiver, bearings are obtained modulo
180 rather than 360.  In this case, change the constants
(180,360) in Mitsuta to (90,180), and in the direction cosine
method, multiply the raw angles by 2 in the arguments to sin
and cos, and divide the arctangent result by 2.

For the remainder of these notes, assume that we are working
mod 180.

Both methods give reasonable results but have some caveats.

The Mitsuta method breaks down if consecutive raw readings
differ by more than 90 degrees.  This can easily occur in
the VLF rx, eg when a raw sample is upset by a strong sferic,
or if the signal is just plain weak.  In this situation, the
algorithm doesn't degrade gracefully at all, and gives wild
averages as soon as the range of raw bearings exceeds +/-45
deg of the true bearing.  The direction cosine method is much
better in this respect, and continues to give a useful average,
even when the range of raw bearings is as high as +/- 75 degrees
of the actual bearing, but then rapidly deteriorates above that.

The Mitsuta method does have two advantages.

a) It can be calculated without trig functions, which may be
important in some situations.

b) The Mitsuta mean is mathematically correct, in that the
total deviation from the mean of all the samples is zero.
This is not so of the direction cosine method.

Standard deviation
------------------
It is useful to produce a measure of by how much the raw
readings are spread around the mean.   This is easy to add to
the Mitsuta method:-

   sum = D = b[1];
   sumsq = D * D;

   for( i from 2 to N)
   {
      ...

      sum = sum + D;
      sumsq = sumsq + D * D;
   }
  
   mean = sum/N;
   std_dev = sqrt( sumsq/N - mean * mean);

The direction cosine method does not offer a proper standard
deviation.   But there is a quite accurate approximation for
the standard deviation which can be used here - the Yamartino
method:

For modulo 360 angles:
   sa = sum_sin/N;
   ca = sum_cos/N;
   e = sqrt( 1 - sa * sa - ca * ca);
   std_dev = asin( e) * (1 + 0.1547 * e * e * e);

For modulo 180 angles, the return value from the arcsine must
be divided by 2.

Here's a comparison of Mitsuta and direction cosine on a VLF
spectrum,

 http://abelian.org/vlf/100107a.mitsuta.png
 http://abelian.org/vlf/100107a.dircos.png
 
Blue crosses are the bearing averages for each bin.  For some
regions of the spectrum containing only noise, the direction
cosine method still manages to extract a fairly decent average
bearing, whereas the Mitsuta mean degenerates into random
bearings.  Direction cosines succeeds on some of the weaker
signals (example near 77.5kHz) where Mitsuta fails.

Here is a comparison table for various widths of raw bearing
distributions.   Bearings are averaged over 2812 raw samples
to simulate the behaviour of the VLF rx, and the true bearing
in this simulation is 90 deg.  Standard deviation in brackets.

Real deviation    Mitsuta          Direction Cosines
----------------------------------------------------
    5            90.061 (5.016)       90.061 (5.015)
   10            90.123 (10.032)      90.122 (10.054)
   20            90.245 (20.064)      90.237 (20.521)
   25            90.306 (25.080)      90.288 (26.044)
   30           134.920 (131.514)     90.327 (31.670)
   40            93.371 (136.339)     90.284 (42.414)
   50           106.231 (126.675)     86.753 (50.987)
   55            77.432 (125.161)    176.870 (50.221)

Mitsuta stops working when the standard deviation of the raw
bearings is above about 25, the direction cosines return valid
results up to a deviation of 50.

A problem with the Yamartino estimate of standard deviation is
that when the actual standard deviation exceeds about 55, the
Yamartino estimate begins to return lower numbers.  Thus, for
example if the receiver outputs these averages,

   50           106.231 (126.675)     86.753 (50.987)
   55            77.432 (125.161)    176.870 (50.221)
   60           101.354 (127.370)      1.489 (47.799)

it is not possible to tell that the first line is a good
estimate of the source and the second two are completely wrong.

The only safe thing to do is to discard averages where the
Yamartino deviation is greater than about 40.  Unfortunately
this means discarding some valid results too.