Click here to Skip to main content
15,867,330 members
Articles / Programming Languages / C++

Circular Values Math and Statistics with C++11

Rate me:
Please Sign up or sign in to vote.
4.90/5 (80 votes)
27 Apr 2013CPOL30 min read 118.2K   1.4K   129   29
A C++11 infrastructure for circular values (angles, time-of-day, etc.) mathematics and statistics

1. Introduction

Many scientific and engineering problems involve circular values. Most of the times these values are orientation (angles) or cyclical timing (time-of-day), but it may be any other circular quantity. Such problems are very common in physics, geodesics and navigation, but also appear in other fields such as psychology, criminology (time-of-day statistics), bird-watching and biology (directional statistics and time-of-year statistics).

In contrast to linear values (values on the real line), circular values have a limited range, and a wrap-around property.

Mathematics and Statistics of circular values is very tricky and error prone, both for simple operations such as addition and subtraction, and for more complex operators such as calculating average and median, estimations and interpolations.

For angles, it is even more error prone, since sometimes the range Circular-Values/1.png is used, sometimes Circular-Values/2.png. Sometimes it is Circular-Values/3.png and sometimes Circular-Values/4.png - even in the same application. For time-of-day computations, the range may be Image 5 or something else. We need a scheme for using several representations in the same application, and to safely perform operations between them.

2. Features

In this article I am going to present:

  • Definition of circular value type and its associated operators
  • A C++11 infrastructure for doing mathematics with circular values
  • C++11 statistical distribution classes for circular values: Wrapped Normal and Wrapped Truncated Normal
  • Average, Weighted average and Median of circular values
  • Circular parameter estimation based on noisy independent measurements
  • Interpolation and average estimation of sampled continuous-time circular signal
This is an original work. To the best of my knowledge, some of the ideas mentioned here were not published before.

3. Requirements

To compile the code, you’ll need Visual C++ 2012.
The given code relies heavily of C++11 features.
The code may be easily converted to any other C++11 compiler.
The code may be converted to plain old C++. However, it will lose some of its charm.

4. Some helper functions

We’ll define the following helper functions:

  • Square function: Sqr(r1). This is simple:
    C++
    // square (x*x)
    template <typename T>
    T Sqr(const T& x)
    {
        return x*x;
    }
  • Floating-point modulo function: reminder = mod(dividend, divisor)

    The dividend, divisor and the resulting reminder are all floating-point values.

    The quotient and the reminder must satisfy:
    - quotient is integer
    - |reminder| < |divisor|
    - dividend = quotientdivisor + reminder

    Still, many definitions are possible:
    A. The reminder has the same sign as the dividend
    - Implies that the quotient rounds towards zero
    B. The reminder has the same sign as the divisor
    - Implies that the quotient rounds toward negative infinity
    C. The reminder has a different sign than the dividend
    - Implies that the quotient rounds away from zero
    D. The reminder has a different sign than the divisor
    - Implies that the quotient rounds toward positive infinity
    E. The sign of the reminder is always non-negative
    F. The sign of the reminder is always non-positive
    G. The reminder is closest to zero
    - Implies that the quotient is the integer nearest to the exact value of division's result
    IEEE 754 names this REM operator, and disambiguate it for a boundary case.
    When the fractional part of the division's result is exactly 0.5, the quotient is even.
    H. The reminder is farthest from zero

    Values for (quotient, reminder) according to the various implementations:

    Image 6

    In addition, we should define the behavior when the divisor is 0:
    X. When the divisor is 0, the reminder is undefined
    Y. When the divisor is 0, the reminder is defined as 0
    Z. When the divisor is 0, the reminder is defined as the dividend

    Most programming languages, numerical computing environment, and spreadsheets, use definitions A, B or G. For example:
    - C (ISO/IEC 9899) and C++ (ISO/IEC 14882) fmod() function implements the AX or AY definition
    (the standard allow either)
    - C (starting from ISO/IEC 9899-1999) and C++ (starting from ISO/IEC 14882-2011,
    but not yet supported by many compiles) defines an additional function - reminder(),
    which implements definition GX or GY (the standard allow either)
    - Microsoft Excel's MOD() function implements the BX definition.
    - MATLAB's mod() and rem() functions implement the BZ and AX definitions respectively.

    We will use the following definition:

    Definition 1

    Image 7

    This definition was given by Knuth in The Art of Computer Programming, Vol.1 (p.39 of the 3rd edition), and is equivalent to definition BZ.

    A straightforward implementation would be:

    C++
    template<typename T>
    T Mod(T x, T y)
    {
        static_assert(!std::numeric_limits<T>::is_exact , "Mod: floating-point type expected");
     
        if (0 == y)
            return x;
     
        return x - y * floor(x/y);
    } 

    However, such implementation is not resilient to boundary cases resulting from the noncontinuity of the modulo function and the floating-point limited accuracy. Here are two examples for double-precision values:
    - Mod(-1e-16, 360.) = 360. (should be 0)
    The value 360,-1e-16 = 359.9999999999999999 cannot be represented by double-precision
    - Mod(106.81415022205296, _TWO_PI) = -1.421e-14 (should be _TWO_PI - 1.01e-14)

    Therefore, we'll use the following implementation:

    C++
    // Floating-point modulo
    // The result (the remainder) has the same sign as the divisor.
    // Similar to matlab's mod(); Not similar to fmod() -   Mod(-3,4)= 1   fmod(-3,4)= -3
    template<typename T>
    T Mod(T x, T y)
    {
        static_assert(!std::numeric_limits<T>::is_exact , "Mod: floating-point type expected");
    
        if (0 == y)
            return x;
    
        double m= x - y * floor(x/y);
    
        // handle boundary cases resulting from floating-point limited accuracy:
    
        if (y > 0)              // modulo range: [0..y)
        {
            if (m>=y)           // Mod(-1e-16             , 360.    ): m= 360.
                return 0;
    
            if (m<0 )
            {
                if (y+m == y)
                    return 0  ; // just in case...
                else
                    return y+m; // Mod(106.81415022205296 , _TWO_PI ): m= -1.421e-14 
            }
        }
        else                    // modulo range: (y..0]
        {
            if (m<=y)           // Mod(1e-16              , -360.   ): m= -360.
                return 0;
    
            if (m>0 )
            {
                if (y+m == y)
                    return 0  ; // just in case...
                else
                    return y+m; // Mod(-106.81415022205296, -_TWO_PI): m= 1.421e-14 
            }
        }
    
        return m;
     }

    Theorem 1

    Circular-Values/16.png

    Proof

    Circular-Values/17.png
  • Almost-equal function: IsAlmostEqual(r1, r2)

    When implementing unit-testing for floating-point values computations, it is a common need to verify that 2 computations give equal results. Due to the finite-accuracy of floating-point types, many times the values won’t be equal, but ‘almost equal’. For example double d= 300.; is only almost-equal to double e= exp(log(300.));

    To test for almost-equality, we need to consider not the absolute magnitude of the error, but its magnitude relative to the result. Moreover, we need to take care of marginal situations such as INF/NAN values.

    Google Test is a framework for writing C++ tests. I’ve extracted the relevant functionality for testing almost equality into a single file - FPCompare.h, which I tweaked a little. Based on it, I implemented the following functions:

C++
// check if two floating-points are almost equal
template<typename T>
static bool IsAlmostEq(T x, T y) 
{ 
    static_assert(!std::numeric_limits<T>::is_exact , "IsAlmostEq: floating-point type expected"); 
    FloatingPoint<T> f(x);
    FloatingPoint<T> g(y);
    return f.AlmostEquals(g);
}
 
// assert that 2 floating-points are almost equal
static void AssertAlmostEq(const double f, const double g)
{
    assert(IsAlmostEq(f, g));
} 

For more information, look at Google test at http://code.google.com/p/googletest/

5. A circular value type

Definition 2

A Circular value type Circular-Values/18.png is defined by three constants:

  • The range Circular-Values/19.png
  • The zero-value Circular-Values/20.png.

Note the use of a right-open interval range.

For example, for angles in the range Circular-Values/21.png, it is natural to define the zero-value in the middle, while for angles in the range Circular-Values/22.png it is natural to define zero-value on the edge. Generally, the zero-value may be any value in the range – not necessarily 0.

The zero-value has some special properties which will be defined below.

Circular-Values/23.png

We will use the following macro to define a circular value type:

C++
// macro for defining a circular-value type
#define CircValTypeDef(_Name, _L, _H, _Z)             \
    struct _Name                                      \
    {                                                 \
        static const double L  ; /* range: [L,H) */   \
        static const double H  ;                      \
        static const double Z  ; /* zero-value   */   \
        static const double R  ; /* range        */   \
        static const double R_2; /* half range   */   \
                                                      \
        static_assert((_H>_L) && (_Z>=_L) && (_Z<_H), \
            #_Name##": Range not valid");             \
    };                                                \
                                                      \
    const double _Name::L  = (_L)          ;          \
    const double _Name::H  = (_H)          ;          \
    const double _Name::Z  = (_Z)          ;          \
    const double _Name::R  = ((_H)-(_L))   ;          \
    const double _Name::R_2= ((_H)-(_L))/2.;

And use it to define the following circular value types:

C++
// basic circular-value types
CircValTypeDef(SignedDegRange  , -180.,   180.,  0. )
CircValTypeDef(UnsignedDegRange,    0.,   360.,  0. )
CircValTypeDef(SignedRadRange  , -M_PI,   M_PI,  0. )
CircValTypeDef(UnsignedRadRange,    0., 2*M_PI,  0. )

And some additional circular value types - for testing:

C++
// some additional circular-value types - for testing
CircValTypeDef(TestRange0      ,    3.,    10.,  5.3)
CircValTypeDef(TestRange1      ,   -3.,    10., -3.0)
CircValTypeDef(TestRange2      ,   -3.,    10.,  9.9)
CircValTypeDef(TestRange3      ,  -13.,    -3., -5.3)

6. Defining a circular value class with a given type

The templated class CircVal is used for storing a single circular value.
The template parameter should be defined using the CircValTypeDef macro.

C++
// circular-value
// Type should be defined using the CircValTypeDef macro
template <typename Type> 
class CircVal 
{
…
}

Here is an example of defining a circular value variable with a given type:

C++
CircVal<UnsignedDegRange> c1;

7. Range checking and wrapping

CircVal::IsInRange static function is used for testing that a given linear value is within the range of a CircVal class.
C++
static bool IsInRange(double r)
{
    return (r>=Type::L && r<Type::H);
}

The wrap(r) function is used to ‘wrap around’ a linear value to the range of a given CircVal class. For example: for the range Circular-Values/4.png, the value 360 would be wrapper to 0, and 370 would be wrapped to 10. -350 would be wrapped to 10 as well.

Definition 3

Image 17

We will use this for faster implementation.

The CircVal::Wrap static function calculates wrap(r):

C++
// 'wraps' circular-value to [L,H)
static double Wrap(double r)
{
    // the next lines are for optimization and improved accuracy only
    if (r>=Type::L)
    {
             if (r< Type::H        ) return r        ;
        else if (r< Type::H+Type::R) return r-Type::R;
    }
    else
             if (r>=Type::L-Type::R) return r+Type::R;
 
    // general case
    return Mod(r - Type::L, Type::R) + Type::L;
}

8. A Walk

Definition 4

For linear values, a walk Image 18 is defined by:

  • A start point Image 19
  • A directed length Circular-Values/32.png

A walk is a movement along the linear axis, from the start point, with the given directed length.

When the length is positive, we'll call it ‘an increasing walk’ (toward positive infinity).
When the length is negative, we'll call it ‘a decreasing walk’ (toward negative infinity).

The end point Circular-Values/33.png is the point reached at the end of the walk.
We’ll call a walk with a start point s and an end point e “A walk from s to e”.

Circular-Values/34.png

There is only one possible walk from s to e: <s, e-s >

For any circular value type, a walk Circular-Values/30.png is defined by:

  • A start point Circular-Values/31.png
  • A directed length Circular-Values/32.png

A walk is a movement along the circular axis, from the start point, with the given directed length.

When the length is positive, we'll call it ‘an increasing walk’ (depicted clockwise in the picture above).
When the length is negative, we'll call it ‘a decreasing walk’ (depicted counterclockwise in the picture above).

The end point Circular-Values/35.png is the point reached at the end of the walk.
We’ll call a walk with a start point s and an end point e “A walk from s to e”.
Since the walk is wrapped around the circle:

Circular-Values/36.png

There are infinitely many walks from Circular-Values/37.png:

Circular-Values/38.png

The shortest walk from s to e is the one that minimizes |d|. The walk may be increasing or decreasing.

9. Distance between two values

Definition 5

The directed distance from linear value Circular-Values/40.png:

Image 31

is the directed length of the walk from Circular-Values/40.png.

For example, Circular-Values/43.png

Similarly, the directed distance from circular value Circular-Values/45.png:

Circular-Values/46.png

is the directed length of the shortest walk from Circular-Values/45.png:

Circular-Values/47.png

It can easily be deduced that:

Circular-Values/48.png

Since

Circular-Values/49.png

And since

Circular-Values/50.png

We can write:

Circular-Values/51.png

which we will use for faster implementation.

The CircVal::Sdist static function calculates sdist(c1,c2):

C++
// the length of the directed walk from c1 to c2, with the lowest absolute-value length
// return value is in [-R/2, R/2)
static double Sdist(const CircVal& c1, const CircVal& c2)
{
    double d= c2.val-c1.val;
    if (d <  -Type::R_2) return d + Type::R;
    if (d >=  Type::R_2) return d - Type::R;
                         return d          ;
}

Note: To define a metric space, we can use |sdist(c1,c2)| as our metric (we won't prove this here).

For circular values, we'll define an additional distance function:

Circular-Values/52.png, The increasing distance from circular value Image 43, is the length of the shortest increasing-walk (depicted clockwise in the picture below) from Image 44.

It is always in the range Image 45.

Image 46

There is no similar distance function for linear values.

Definition 6

Circular-Values/54.png

which may also be expressed as:

Circular-Values/55.png

It can easily be deduced that:

Circular-Values/56.png Hence, it cannot be used as a metric

The CircVal::Pdist static function calculates pdist(c1,c2):

C++
// the length of the increasing walk from c1 to c2 with the lowest length
// return value is in [0, R)
static double Pdist(const CircVal& c1, const CircVal& c2)
{
    return c2.val>=c1.val ? c2.val-c1.val : Type::R-c1.val+c2.val;
}

Theorem 2

Circular-Values/57.png

Proof

Rephrasing the definitions:

Circular-Values/58.png

Theorem 3

Circular-Values/60.png (Very intuitive)

Proof

According to the definition above:

Circular-Values/61.png

And also

Circular-Values/62.png

, which are equal in all cases.

Theorem 4

For any two circular values: Circular-Values/63.png

Proof

Circular-Values/64.png

Theorem 5

For any two circular values: Circular-Values/65.png

Proof

Circular-Values/66.png

Theorem 6

For any two circular values: Circular-Values/67.png

Proof

Trivial, based on theorems 4, 5

10. Conversion between different types of circular values

We want to be able to perform operations between different types of circular values.

For example:

C++
CircVal<UnsignedDegRange> d1=  10.;
CircVal<SignedDegRange  > d2= -10.;
CircVal<SignedRadRange  > d3      ;

d3= d1+d2;

To convert a circular value Circular-Values/70.png from circular type Circular-Values/68.png to circular type Circular-Values/69.png, we’ll retain the arc-length between Circular-Values/70.png and Circular-Values/71.png. For that, we can use either of these formulae:

Circular-Values/72.png

Theorem 7

I and II are equivalent

Proof

Circular-Values/73.png

Therefore, in both cases I =II

The conversion can be accomplished easily be implementing appropriate constructor and assignment operator:

C++
// construction based on a circular value of another type
// sample use: CircVal<SignedRadRange> c= c2;   -or-   CircVal<SignedRadRange> c(c2);
template<typename CircVal2>
CircVal(const CircVal2& c2)
{
    double val2= c2.Pdist(c2.GetZ(), c2);
    val= Wrap(val2*Type::R/c2.GetR() + Type::Z);
}
 
// assignment from another type of circular value
template<typename CircVal2>
CircVal& operator= (const CircVal2& c2)
{
    double val2= c2.Pdist(c2.GetZ(), c2);
    val= Wrap(val2*Type::R/c2.GetR() + Type::Z);
    return *this;
}

11. Constructing circular value / Fetching circular value

Construction a CircVal objects based on linear value-representation is implemented by the appropriate constructor and assignment operator. Construction of a linear value based on a CircVal object is implemented with operator double() const overloading.

C++
// construction based on a floating-point value
// should only be called when the floating-point is a value in the range!
// to translate a floating-point such that 0 is mapped to Type::Z, call ToC()
CircVal(double r)
{
    val= Wrap(r);
}
 
// assignment from a floating-point value
// should only be called when the floating-point is a value in the range!
// to translate a floating-point such that 0 is mapped to Type::Z, call ToC()
CircVal& operator= (double r)
{
    val= Wrap(r);
    return *this;
}

operator double() const
{
    return val;
} 

12. Circular value operators – Intuitive description

First, here is an intuitive description for each operator (formal definition will be given later):
  • Negation operator (-c) is defined for any circular value. The result is a circular value.
    Negation of a circular value is the symmetric value relative to z (the zero-value).

    Circular-Values/74.png

  • Opposite operator (~c) is defined for any circular value. The result is a circular value.
    Opposite of a circular value is the circular value pointing to the opposite direction.

    Circular-Values/75.png

  • Addition and subtraction operators Circular-Values/76.png are defined for any two circular values.
    The result is a circular value - the point reached and the end of the walk:

    Addition: Starting from c1, walking with length pdist(z,c2).
    Subtraction: Starting from c1, walking with length -pdist(z,c2).

    Image 70

  • Multiplication operator (c·r) is defined for any circular value c and any linear value r.
    The result is a circular value - the point reached and the end of the walk:

    Starting from z, walking distance r·pdist(c,z).

    Circular-Values/77.png

  • Division operator (c/r) is defined for any circular value c and any non-zero linear value r.
    The result is a circular value - the point reached and the end of the walk:
    Starting from z, walking Image 72 times distance pdist(c,z).
  • Trigonometric functions sin, cos, tan are defined for any circular value c, such that the result will be a linear value, in accordance with the common functions usage.
  • Inverse trigonometric functions asin, acose, atan, atan2 shall be defined for any linear number (-1 <= r < 1, for asin, acos), such that the result will be a circular value of our class in accordance with the common functions usage.

13. Conversion between circular values and linear values

We want a conversion function ToR(c), such that any circular value c in our range [l,h) can be mapped to a linear value, and a conversion function ToR(r), such that any linear value r can be mapped to a circular value in our in our range [l,h). Note that these functions are different from the circular value construction/fetching described above.

Requirements to these conversion functions:

Circular-Values/78.png

We’ll use these conversion functions to simplify the requirements and the definitions of several circular value operators.

14. Requirements for circular value operators

We want to define the circular value operators such that the following conditions are satisfied:

Circular-Values/79.png
Circular-Values/80.png

The set of real numbers in the range [l,h) with the + binary operator is a bounded group (z is the identity element, and -c is the inverse element of c).

15. Comparison operators

Defining operators = and Image 76 is straightforward. The other 4 comparison operators {<, >, =<, >=} can be defined in several different ways. Requirements to these comparison operators (For any circular values c1, c2, c3) should satisfy:

Circular-Values/83.png

All these requirements, however, tell us nothing about the meaning of the > operator.

For linear values we can define:

Circular-Values/84.png

Based on this, we can define the circular values operator in different ways:

Circular-Values/85.png

<>Hence, for two circular values c1,c2, for each comparison operator Circular-Values/86.png, we can use:

Circular-Values/87.png

Though all definitions satisfy all the requirements, they mean different things. The 1st compares the pdist from l, the 2nd compares the pdisht from z, and the 3rd compares the sdist from z.

In many practical situations z = 0.

For C = < -180,180,0 > it seems more appropriate to use (i) or (iii) which are equivalent, while for C = < 0,360,0 >, (i) and (ii) which are equivalent fits better. This is why our selected definition is (i), but it may be different for specific use cases.

C++
bool operator==(const CircVal& c) const { return val == c.val; }
bool operator!=(const CircVal& c) const { return val != c.val; }
 
// note that two circular values can be compared in several different ways.
// check carefully if this is really what you need!
bool operator> (const CircVal& c) const { return val >  c.val; }
bool operator>=(const CircVal& c) const { return val >= c.val; }
bool operator< (const CircVal& c) const { return val <  c.val; }
bool operator<=(const CircVal& c) const { return val <= c.val; }

16. Circular value operators – Formal definition

Now, here is the definition that satisfies all the requirements

Definition 7

Circular-Values/88.png

Theorem 8

The - unary operator can be equivalently defined by each of these two equivalent formulae:

Circular-Values/89.png

The + and - binary operators can be equivalently defined by each of these two equivalent formulae:

Circular-Values/90.png

Proof

Image 85 Here is the implementation of these operators:

C++
template <typename Type> 
class CircVal 
{
    double val; // actual value [L,H)
 
public:
    // ---------------------------------------------
    // convert circular-value c to real-value [L-Z,H-Z). .Z is converted to 0
    friend double ToR(const CircVal& c) { return c.val - Type::Z; }
 
    // ---------------------------------------------
    const CircVal  operator+ (                ) const { return val;                              }
    const CircVal  operator- (                ) const { return Wrap(Type::Z-Sdist(Type::Z,val)); } // return negative circular value
    const CircVal  operator~ (                ) const { return Wrap(val+Type::R_2             ); } // return opposite circular-value
 
    const CircVal  operator+ (const CircVal& c) const { return Wrap(val+c.val        - Type::Z); }
    const CircVal  operator- (const CircVal& c) const { return Wrap(val-c.val        + Type::Z); }
    const CircVal  operator* (const double&  r) const { return Wrap((val-Type::Z)*r  + Type::Z); }
    const CircVal  operator/ (const double&  r) const { return Wrap((val-Type::Z)/r  + Type::Z); }
 
          CircVal& operator+=(const CircVal& c)       { val= Wrap(val+c.val          - Type::Z); return *this; }
          CircVal& operator-=(const CircVal& c)       { val= Wrap(val-c.val          + Type::Z); return *this; }
          CircVal& operator*=(const double&  r)       { val= Wrap((val-Type::Z)*r    + Type::Z); return *this; }
          CircVal& operator/=(const double&  r)       { val= Wrap((val-Type::Z)/r    + Type::Z); return *this; }
 
          CircVal& operator =(const CircVal& c)       { val= c.val                             ; return *this; }
    …
}
 
// ==========================================================================
template <typename Type> static double        sin  (const CircVal<Type>& c) { return std::sin(ToR(CircVal<SignedRadRange>(c)));  }

template <typename Type> static double        cos  (const CircVal<Type>& c) { return std::cos(ToR(CircVal<SignedRadRange>(c)));  }

template <typename Type> static double        tan  (const CircVal<Type>& c) { return std::tan(ToR(CircVal<SignedRadRange>(c)));  }

template <typename Type> static CircVal<Type> asin (double r              ) { return CircVal<SignedRadRange>(std::asin (r    )); }

template <typename Type> static CircVal<Type> acos (double r              ) { return CircVal<SignedRadRange>(std::acos (r    )); }

template <typename Type> static CircVal<Type> atan (double r              ) { return CircVal<SignedRadRange>(std::atan (r    )); }

template <typename Type> static CircVal<Type> atan2(double r1, double r2  ) { return CircVal<SignedRadRange>(std::atan2(r1,r2)); }

template <typename Type> static CircVal<Type> ToC  (double r              ) { return CircVal<Type>::Wrap(r + Type::Z);           }

Note that for asin, acos, atan and atan2 functions a template parameter should be used. For example:
CircVal<SignedDegRange> d1= asin<SignedDegRange>(0.5);d1= asin(0.5) won’t give us the expected result unless our range is SignedRadRange.

17. Testing for near-equality of circular values

Based on the previously describe test for near-equality for two linear values, we will construct a near-equality test for two circular values:
C++
template <typename Type>
class CircValTester
{
    // check if 2 circular-values are almost equal
    static bool IsCircAlmostEq(const CircVal<Type>& _f, const CircVal<Type>& _g)
    {
        double f= _f;
        double g= _g;
 
        if (::IsAlmostEq(f, g))
            return true;
 
        if (f < g)
            return IsAlmostEq(f, g - Type::R);
        else
            return IsAlmostEq(f, g + Type::R);
    }
 
    // assert that 2 circular-values are almost equal
    static void AssertCircAlmostEq(const CircVal<Type>& f, const CircVal<Type>& g)
    {
        assert(IsCircAlmostEq(f, g));
    }
    …
}

18. Testing correctness of circular value class implementation

Using test-driven design, we’re check that our implementation fulfills our requirements.
The Test() function generates 10,000 random test-cases, and verify that all the requirements holds.
We will use C++11 random generation, which is a very useful addition to the language:

C++
// tester for CircVal class
template <typename Type>
class CircValTester
{
    static void Test()
    {
        CircVal<Type> ZeroVal= Type::Z;
 
        // --------------------------------------------------------
        AssertCircAlmostEq(ZeroVal       , -ZeroVal);
 
        AssertAlmostEq    (sin(ZeroVal)  , 0.      );
        AssertAlmostEq    (cos(ZeroVal)  , 1.      );
        AssertAlmostEq    (tan(ZeroVal)  , 0.      );
 
        AssertCircAlmostEq(asin<Type>(0.), ZeroVal );
        AssertCircAlmostEq(acos<Type>(1.), ZeroVal );
        AssertCircAlmostEq(atan<Type>(0.), ZeroVal );
 
        AssertCircAlmostEq(ToC<Type>(0)  , ZeroVal );
        AssertAlmostEq    (ToR(ZeroVal)  , 0.      );
 
        // --------------------------------------------------------
        std::default_random_engine             rand_engine                 ;
        std::uniform_real_distribution<double> c_uni_dist(Type::L, Type::H);
        std::uniform_real_distribution<double> r_uni_dist(0.     , 1000.  ); // for multiplication,division by real-value
        std::uniform_real_distribution<double> t_uni_dist(-1.    , 1.     ); // for inverse-trigonometric functions
        
        std::random_device rnd_device;
        rand_engine.seed(rnd_device()); // reseed engine
 
        for (unsigned i= 10000; i--;)
        {
            CircVal<Type> c1(c_uni_dist(rand_engine)); // random circular value
            CircVal<Type> c2(c_uni_dist(rand_engine)); // random circular value
            CircVal<Type> c3(c_uni_dist(rand_engine)); // random circular value
            double        r (r_uni_dist(rand_engine)); // random real     value [    0, 1000) - for testing *,/ operators
            double        a1(t_uni_dist(rand_engine)); // random real     value [   -1,    1) - for testing asin,acos
            double        a2(t_uni_dist(rand_engine)); // random real     value [-1000, 1000) - for testing atan
 
            assert            (c1                                 == CircVal<Type>((double)c1)         );
 
            AssertCircAlmostEq(+c1                                  , c1                               ); // +c         = c
            AssertCircAlmostEq(-(-c1)                               , c1                               ); // -(-c)      = c
            AssertCircAlmostEq(c1 + c2                              , c2 + c1                          ); // c1+c2      = c2+c1
            AssertCircAlmostEq(c1 + (c2 +c3)                        , (c1 + c2) + c3                   ); // c1+(c2+c3) = (c1+c2)+c3
            AssertCircAlmostEq(c1 + -c1                             , ZeroVal                          ); // c+(-c)     = z
            AssertCircAlmostEq(c1 + ZeroVal                         , c1                               ); // c+z        = c
 
            AssertCircAlmostEq(c1      -  c1                        , ZeroVal                          ); // c-c        = z
            AssertCircAlmostEq(c1      - ZeroVal                    , c1                               ); // c-z        = c
            AssertCircAlmostEq(ZeroVal - c1                         , -c1                              ); // z-c        = -c
            AssertCircAlmostEq(c1      - c2                         , -(c2 - c1)                       ); // c1-c2      = -(c2-c1)
 
            AssertCircAlmostEq(c1 * 0.                              , ZeroVal                          ); // c*0        = 0
            AssertCircAlmostEq(c1 * 1.                              , c1                               ); // c*1        = c
            AssertCircAlmostEq(c1 / 1.                              , c1                               ); // c/1        = c
 
            AssertCircAlmostEq((c1 * (1./(r+1.))) / (1./(r+1.))     , c1                               ); // (c*r)/r    = c, 0<r<=1
            AssertCircAlmostEq((c1 / (    r+1.) ) * (    r+1. )     , c1                               ); // (c/r)*r    = c,   r>=1
 
            // --------------------------------------------------------
            AssertCircAlmostEq(~(~c1)                               , c1                               ); // opposite(opposite(c) = c
            AssertCircAlmostEq(c1 - (~c1)                           , ToC<Type>(Type::R/2.)            ); // c - ~c               = r/2+z


            // --------------------------------------------------------
            AssertAlmostEq    (sin(ToR(CircVal<SignedRadRange>(c1))),  sin(c1)                         ); // member func sin
            AssertAlmostEq    (cos(ToR(CircVal<SignedRadRange>(c1))),  cos(c1)                         ); // member func cos
            AssertAlmostEq    (tan(ToR(CircVal<SignedRadRange>(c1))),  tan(c1)                         ); // member func tan
 
            AssertAlmostEq    (sin(-c1)                             , -sin(c1)                         ); // sin(-c)    = -sin(c)
            AssertAlmostEq    (cos(-c1)                             ,  cos(c1)                         ); // cos(-c)    =  cos(c)
            AssertAlmostEq    (tan(-c1)                             , -tan(c1)                         ); // tan(-c1)   = -tan(c) error may be large

            AssertAlmostEq    (sin(c1+ToC<Type>(Type::R/4.))        ,  cos(c1)                         ); // sin(c+r/4) =  cos(c)
            AssertAlmostEq    (cos(c1+ToC<Type>(Type::R/4.))        , -sin(c1)                         ); // cos(c+r/4) = -sin(c)
            AssertAlmostEq    (sin(c1+ToC<Type>(Type::R/2.))        , -sin(c1)                         ); // sin(c+r/2) = -sin(c)
            AssertAlmostEq    (cos(c1+ToC<Type>(Type::R/2.))        , -cos(c1)                         ); // cos(c+r/2) = -cos(c)
 
            AssertAlmostEq    (Sqr(sin(c1))+Sqr(cos(c1))            , 1.                               ); // sin(x)^2+cos(x)^2 = 1
 
            AssertAlmostEq    (sin(c1)/cos(c1)                      , tan(c1)                          ); // sin(x)/cos(x) = tan(x)
 
            // --------------------------------------------------------
            AssertCircAlmostEq(asin<Type>(a1)                       , CircVal<SignedRadRange>(asin(a1))); // member func asin
            AssertCircAlmostEq(acos<Type>(a1)                       , CircVal<SignedRadRange>(acos(a1))); // member func acos
            AssertCircAlmostEq(atan<Type>(a2)                       , CircVal<SignedRadRange>(atan(a2))); // member func atan
 
            AssertCircAlmostEq(asin<Type>(a1) + asin<Type>(-a1)     , ZeroVal                          ); // asin(r)+asin(-r) = z
            AssertCircAlmostEq(acos<Type>(a1) + acos<Type>(-a1)     , ToC<Type>(Type::R/2.)            ); // acos(r)+acos(-r) = r/2+z
            AssertCircAlmostEq(asin<Type>(a1) + acos<Type>( a1)     , ToC<Type>(Type::R/4.)            ); // asin(r)+acos( r) = r/4+z
            AssertCircAlmostEq(atan<Type>(a2) + atan<Type>(-a2)     , ZeroVal                          ); // atan(r)+atan(-r) = z
 
            // --------------------------------------------------------
            assert            (c1 >  c2                           ==    (c2 <  c1)                     ); // c1> c2 <==>   c2< c1
            assert            (c1 >= c2                           ==    (c2 <= c1)                     ); // c1>=c2 <==>   c2<=c1
            assert            (c1 >= c2                           ==  ( (c1 >  c2) ||  (c1 == c2))     ); // c1>=c2 <==>  (c1> c2)|| (c1==c2)
            assert            (c1 <= c2                           ==  ( (c1 <  c2) ||  (c1 == c2))     ); // c1<=c2 <==>  (c1< c2)|| (c1==c2)
            assert            (c1 >  c2                           ==  (!(c1 == c2) && !(c1 <  c2))     ); // c1> c2 <==> !(c1==c2)&&!(c1< c2)
            assert            (c1 == c2                           ==  (!(c1 >  c2) && !(c1 <  c2))     ); // c1= c2 <==> !(c1> c2)&&!(c1< c2)
            assert            (c1 <  c2                           ==  (!(c1 == c2) && !(c1 >  c2))     ); // c1< c2 <==> !(c1==c2)&&!(c1> c2)
            assert            (!(c1>c2) || !(c2>c3) || (c1>c3)                                         ); // (c1>c2)&&(c2>c3) ==> c1>c3

            // --------------------------------------------------------
            AssertCircAlmostEq(c1                                   , ToC<Type>(ToR( c1)       )       ); //  c1        = ToC(ToR( c1)
            AssertCircAlmostEq(-c1                                  , ToC<Type>(ToR(-c1)       )       ); // -c1        = ToC(ToR(-c1)
            AssertCircAlmostEq(c1 + c2                              , ToC<Type>(ToR(c1)+ToR(c2))       ); // c1+c2      = ToC(ToR(c1)+ToR(c2))
            AssertCircAlmostEq(c1 - c2                              , ToC<Type>(ToR(c1)-ToR(c2))       ); // c1-c2      = ToC(ToR(c1)-ToR(c2))
            AssertCircAlmostEq(c1 * r                               , ToC<Type>(ToR(c1)*r      )       ); // c1*r       = ToC(ToR(c1)*r      )
            AssertCircAlmostEq(c1 / r                               , ToC<Type>(ToR(c1)/r      )       ); // c1/r       = ToC(ToR(c1)/r      )

            // --------------------------------------------------------
        }
    }

public:
    CircValTester()
    {
        Test();
    }
};

The code below executes these tests for all 8 circular-value defined types:

C++
CircValTester<SignedDegRange  > testA;
CircValTester<UnsignedDegRange> testB;
CircValTester<SignedRadRange  > testC;
CircValTester<UnsignedRadRange> testD;
 
CircValTester<TestRange0      > test0;
CircValTester<TestRange1      > test1;
CircValTester<TestRange2      > test2;
CircValTester<TestRange3      > test3;

19. Using the CircVal class

This sample code below speaks for itself:

C++
CircVal<SignedDegRange  > d1=  10. ;
CircVal<UnsignedRadRange> d2=   0.2;
CircVal<SignedDegRange  > d3= d1+d2;
 
d1+= 355.;
double d= d1;
 
d = sin(d1) / cos(d2) + tan(d3);
d1= asin<SignedDegRange>(0.5);

The code is simple, yet powerful.

20. C++11 Distribution classes for circular values

One of the additions to the C++ standard is the <random> header. It is now possible to generate pseudo-random values with a given distribution. Many distributions classes are provided: Bernoulli, Binomial, Cauchy, Chi-Square, Exponential, Extreme Value, Fisher F, Gamma, Geometric, Log-Normal, Negative-Binomial, Normal, Poisson, Student T, Uniform, Weibull, and some interval-related distributions.

Circular distribution classes, however, are not implemented. Some Circular distributions are von Mises, Wrapped Normal and Wrapped Cauchy.

Another distribution, which is very useful when working with circular values, is the truncated normal distribution.

We have implemented three distribution classes: wrapped normal, truncated normal and wrapped truncated normal.

21. Wrapped Normal Distribution class

A wrapped normal distribution results from the "wrapping" of the normal distribution around the range Circular-Values/91.png. Therefore, the probability density function Circular-Values/92.png of the wrapped normal distribution is:

Circular-Values/93.png

The PDF is described in the following image, for same Circular-Values/94.png, but different values of Circular-Values/95.png:

Circular-Values/96.png

This is the distribution of the modulo of a normally-distributed value.

Example: Assume that we have clocks (date + time-of-day), which have a normal-distributed error. It’s time-of-day [00:00, 24:00) distribution would be wrapped-normal.

Wrapped normal distributed values can be generated as follows:

C++
#include "WrappedNormalDist.h"

std::default_random_engine rand_engine;
std::random_device         rnd_device ;
rand_engine.seed(rnd_device()); // reseed engine
 
double fAvrg =    0.;
double fSigma=   45.;
double fL    = -180.; // wrapping-range lower-bound
double fH    =  180.; // wrapping-range upper-bound
 
wrapped_normal_distribution<double> r_wrp(fAvrg, fSigma, fL, fH);
double r1= r_wrp(rand_engine); // random value

The implementation is based on VC 2012 implementation of the normal distribution, added wrapping.

22. Truncated Normal Distribution class

This class has nothing to do with circular values; however it is a base for the wrapped truncated normal distribution class described in the next section.

A truncated normal distribution is the probability distribution of a normally distributed random variable whose value is bounded to the range Image 92. The PDF is described in the following image, for different values of Image 93:

Circular-Values/97.png

Many phenomena behave, or modeled by a truncated normal distribution. One type of examples is the estimation of quantities by “ordinary” humans, such as today’s temperature: if it is 25° Celsius, no “ordinary” human will estimate it below 10°, or above 40°. Another example is the life span of many species (based on data from Shiro Horiuchi - Interspecies Differences in the Life Span Distribution: Humans versus Invertebrates, The Population Council 2003).

Truncated normal distribution is also useful for generating pseudo-random values such as the described above.

Truncated normal distributed values can be generated as follows:

C++
#include "TruncNormalDist.h"

std::default_random_engine rand_engine;
std::random_device         rnd_device ;
rand_engine.seed(rnd_device()); // reseed engine
 
double fAvrg =    0.;
double fSigma=   45.;
double fA    =  -40.; // truncation-range lower-bound
double fB    =   40.; // truncation-range upper-bound
 
truncated_normal_distribution<double> r_trn(fAvrg, fSigma, fA, fB);
double r2= r_trn(rand_engine); // random value

The implementation is based on VC 2012 implementation of the normal distribution as a skeleton, and on C. H. Jackson's R's implementation of the following paper: “Robert, C. P. - Simulation of truncated normal variables. Statistics and Computing (1995) 5, 121–125”

23. Wrapped Truncated Normal Distribution class

Same as for linear values, circular values wrapped normal distribution may be truncated as well. Note that the normal distribution is first truncated to the range Image 95, and then wrapped around the range Circular-Values/91.png

Many phenomena behave, or modeled by a wrapped truncated normal distribution. One type of examples is the estimation of circular quantities by “ordinary” humans, such as the time of sunrise tomorrow, or the azimuth between two drawn points on a piece of paper: If it is 45°, no “ordinary” human will estimate it below 10° or above 80°. Another example may be the exact time-of-day at which John Doe wakes up on an “ordinary” working day.

Wrapped truncated normal distribution is also useful for generating pseudo-random values such as the described above.

Usually, the truncation range is such that b-a < h-l, however, this is not a necessity. For example, we may model the time-of-day error distribution of a given clock that was set 1 year ago, and its absolute error now is modeled by a truncated normal distribution of [-50,50) hours.

After wrapping, it is possible that wrap(b) < wrap(a) - See picture below, which may depict human estimation of an azimuth between two points, where the real azimuth is 0°.

Circular-Values/98.png

Wrapped truncated normal distributed values can be generated as follows:

C++
#include "WrappedTruncNormalDist.h"
 
std::default_random_engine rand_engine;
std::random_device         rnd_device ;
rand_engine.seed(rnd_device()); // reseed engine
 
// normal distribution is first truncated, and then wrapped
 
double fAvrg =    0.;
double fSigma=  100.;
double fA    = -500.; // truncation-range lower-bound
double fB    =  500.; // truncation-range upper-bound
double fL    =    0.; // wrapping  -range lower-bound
double fH    =  360.; // wrapping  -range upper-bound
 
wrapped_truncated_normal_distribution<double> r_ wrp_trn(fAvrg, fSigma, fA, fB, fL, fH);
double d= r_wrp_trn(rand_engine); // random value

24. 3 types of statistical problems

There are 3 distinct types of problems, which tends to be confused one with the other, especially when it comes to circular values. Although these problems may seem similar, the mathematical methods required to solve each of them should be considered separately.

For each problem type, we’ll give a description, and sample problem for linear values and for circular values.

Problem Type 1

Given [circular] values x1..xN – calculate their mean

-or-

Given a random sample x1..xN from a [circular] random-variable X - calculate the sample mean

Usage: Averaging values

Examples:

  • Linear: Given the number of births occurred in the US for each day in the year 2000 -
    Calculate the mean number of births per day
  • Circular: Given time-of-day [00:00-24:00) for each birth occurred in US in the year 2000 -
    Calculate the mean time-of-day

Problem Type 2

Given a multiset of measurements/observations with a [wrapped-/wrapped-truncated-]normal-distributed error of a [circular] parameter – Estimate the parameter

Usage: Estimation of an unknown constant value, based on a multiset of noisy measurements, where the noise has a [wrapped-/wrapped-truncated-]normal-distribution

Examples:

  • Linear: Given a multiset of distance measurements from a stationary transmitter to a stationary receiver, using a measurement technique with a normal distributed error – Estimate the distance.
  • Circular: Given a multiset of direction measurements from a stationary transmitter to a stationary receiver, using a measurement technique with a wrapped normal distributed error – Estimate the direction.
  • Linear: Given a multiset of distance estimates between two points, made by “ordinary” humans (assuming to subject to truncated normal distributed error) - Estimate the distance.
  • Circular: Given a multiset of azimuth estimates between two points, made by “ordinary” humans (assuming to subject to a wrapped truncated normal distributed error) – Estimate the direction.

Problem Type 3

Given a multiset of samples x1..xN of a continuous-time [circular] signal, sampled at times t1..tN respectively, where the measurement technique/instrument is accurate - Estimate the mean value of the signal at period [t1, tN]

Usage: Estimating the average of a continuous-signal, in a given period, based on a set of time-tagged samples

Examples:

  • Linear: Given airplane velocity measured each 1 second [MPH] over a period of one hour – Estimate the mean velocity.
  • Circular: Given airplane heading measured each 1 second over a period of one hour - Estimate the mean heading.

25. Averaging 2 circular values

To average two linear values, Circular-Values/39.png, we can start at Circular-Values/7.png and take a walk with length Circular-Values/99.png. The end-point of the walk is the average. We can, by symmetry, start at Circular-Values/8.png and take a walk with length Image 102.

Circular-Values/100.png

For circular values, given the range [0,360), the average of 330 and 30 is {0} and not {(330+30)/2} = {180}. The average is in the middle of the walk with the lowest absolute-value length, from 330 to 30 (or from 30 to 330).

The average of circular values is a not a single circular value, but a set of circular values.
average(0,180) = {90, 270} - There is no reason to prefer one of these values over the other.

We will use the following formula:

Circular-Values/101.png

26. Averaging n circular values

Circular values average calculation is really tricky. Observing algorithms described in many textbooks, and code written by many programmers, I've almost never seen a correct method used.

Let us start with linear values. Here is our reference problem:

  • Given the number of births occurred in the US for each day in the year 2000 -
    Calculate the mean number of births per day

The solution here is quite simple: Just sum the number of births, and divide it by the number of days.

Why does it work?

The general definition of the mean (arithmetic-average) of values Image 105 for any mathematical field F, is the value of x∈F that minimizes Image 106. Formally: Image 107.

 

For linear values, we can calculate the average as follow:
Image 108
Let y denote the expression we want to minimize.

Image 109
To find the value of x that minimizes y, we need the derivative Image 110 to be 0:

Circular-Values/108.png

Image 112 

We got the well-known arithmetic mean formula for linear values.

 

Now let’s talk about circular values. Here is our reference problem:

  • Given time-of-day [00:00-24:00) for each birth occurred in US in the year 2000 - Calculate the mean time-of-day

For making the following explanation simple, let us consider the range [0,360).

Note that the true average of 0, 0, and 90 should be 30, and not atan Circular-Values/110.png, as suggested in most references.

Based on theorem 3, the absolute value of the distance between 2 circular values in the range [0,360) is
Image 114

This formula is suitable, since we want to minimize the difference between a value and the average, regardless of which ‘direction’ of the average the value is.

Let y denote the expression we want to minimize.

Circular-Values/112.png

The values of x that minimizes y is our average.

Here are some examples:

Example 1: For the input {0,30,60,90}, we’ll plot Circular-Values/113.png as a function of x:

Circular-Values/114.png

We can see that the average is {45} which matches out intuition.

Example 2: For the input {0,90,180,270}, we’ll plot Circular-Values/115.png as a function of x:

Circular-Values/116.png

In this case the average is {45.135, 225, 315}
(For all these values, Circular-Values/117.png is equal).

Example 3: For the input {30,130,230,330}, we’ll plot Image 121 as a function of x:

Circular-Values/118.png

In this case, the average is {0}. Again, it matches our intuition.

The following code was used to collect the data for the above graph:

C++
vector<CircVal<UnsignedDegRange>> Angles2;
Angles2.push_back(CircVal<UnsignedDegRange>( 30.));
Angles2.push_back(CircVal<UnsignedDegRange>(130.));
Angles2.push_back(CircVal<UnsignedDegRange>(230.));
Angles2.push_back(CircVal<UnsignedDegRange>(330.));
auto y= CircAverage(Angles2); 
 
ofstream f0("log0.txt"); 
 
for (double x= 0.; x<=360.; x+= 0.1)
{
    double fSum= 0; 
    for (const auto& a : Angles2)
        fSum+= Sqr(__min(abs(x-a), 360.-abs(x)));
    f0 << x << "\t" << fSum << endl;
} 

To find the values of that x minimizes y, we need the derivative Circular-Values/119.png to be 0.
However, due to the min and abs operators, the expression Circular-Values/120.png is not derivable.

Let us try to rephrase the expression.

We’ll divide the multiset into 3 distinct multi-subsets:

Circular-Values/121.png

We’ll denote with a, b, c, d – the number of elements in A, B, C, D respectively.
Of course Circular-Values/122.png

Circular-Values/123.png

The minimum of y may be where Circular-Values/124.png, or where y is not derivable.

Since in each sector y(x) is a 2nd order polynomial with a positive coefficient for the x² term, it has a single minima, which is when Circular-Values/125.png. In the sample graphs above, it is easy to see the sectors - each sector is a trimmed 2nd order curve. The points where y(x) is not derivable, is on the maximum of each sector.

So, by dividing the input into 3 multisets, we can find the values of x that minimizes y. This is our average.

However, the division of the input elements into 3 multisets is dependent on the value of x.
Here is a graphic visualization of these multisets, for different values of x:

Image 130

If x = 180, all values are in multiset B.
If x is in the range [0,180), multiset C range is (x + 180,360) and multiset D is empty.
If x is in the range (180,360), multiset D range is [0, x - 180) and multiset C is empty.
Multisets C and D cannot coexist.

Here is the average calculation algorithm:

  • Divide the circle into sectors, such that if is within a sector - the elements of multisets B, C and D does not change. The maximal number of sectors equals to the number of elements in A.
  • For each sector: Assume that Circular-Values/127.png (the value of x that minimizes y) is within this sector, calculate its value, and check that the value is indeed within this sector. If so, calculate y as follows:
    1. Assuming x=180, all values are in multiset B:

      Circular-Values/128.png

    2. Assuming x is in the range [0,180), multiset C range is (x+180,360) and multiset D is empty:

      Circular-Values/129.png

    3. Assuming x is in the range (180,360), multiset D range is [0,x-180) and multiset C is empty:

      Circular-Values/130.png

  • Find the sectors(s) that has the lowest y; return their x.

The Function CircAverage calculates the average set of a given vector of circular values.
The return value is a set of all average values:

C++
// calculate average set of circular-values
// return set of average values
// T is a circular value type defined with the CircValTypeDef macro
template<typename T>
set<CircVal<T>> CircAverage(vector<CircVal<T>> const& A)
{
    set<CircVal<T>> MinAvrgVals      ; // results set
 
    // ----------------------------------------------
    // all vars: UnsignedDegRange [0,360)
    double          fSum         = 0.; // of all elements of A
    double          fSumSqr      = 0.; // of all elements of A
    double          fMinSumSqrDiff   ; // minimal sum of squares of differences
    vector<double>  LowerAngles      ; // ascending   [  0,180)
    vector<double>  UpperAngles      ; // descending  (360,180)
    double          fTestAvrg        ;
 
    // ----------------------------------------------
    // local functions - implemented as lambdas
    // ----------------------------------------------
 
    // calc sum(dist(180, Bi)^2) - all values are in set B
    // dist(180,Bi)= |180-Bi|
    // sum(dist(x, Bi)^2) = sum((180-Bi)^2) = sum(180^2-2*180*Bi + Bi^2) = 180^2*A.size - 360*sum(Ai) + sum(Ai^2)
    auto SumSqr = [&]() -> double
    {
        return 32400.*A.size() - 360.*fSum + fSumSqr;
    };
 
    // calc sum(dist(x, Ai)^2). A=B+C; set D is empty
    // dist(x,Bi)= |x-Bi|
    // dist(x,Ci)= 360-(Ci-x)
    // sum(dist(x, Bi)^2)= sum(     (x-Bi) ^2)= sum(        Bi^2 + x^2                      - 2*Bi*x)
    // sum(dist(x, Ci)^2)= sum((360-(Ci-x))^2)= sum(360^2 + Ci^2 + x^2 - 2*360*Ci + 2*360*x - 2*Ci*x)
    // sum(dist(x, Bi)^2) + sum(dist(x, Ci)^2) = nCountC*360^2 + sum(Ai^2) + nCountA*x^2 - 2*360*sum(Ci) + nCountC*2*360*x - 2*x*sum(Ai)
    auto SumSqrC= [&](double x, size_t nCountC, double fSumC) -> double
    {
        return x*(A.size()*x - 2*fSum) + fSumSqr - 2*360.*fSumC + nCountC*( 2*360.*x + 360.*360.);
    };
 
    // calc sum(dist(x, Ai)^2). A=B+D; set C is empty
    // dist(x,Bi)= |x-Bi|
    // dist(x,Di)= 360-(x-Di)
    // sum(dist(x,Bi)^2)= sum(    (x-Bi)^2)= sum(        Bi^2 + x^2                      - 2*Bi*x)
    // sum(dist(x,Di)^2)= sum(360-(x-Di)^2)= sum(360^2 + Di^2 + x^2 + 2*360*Di - 2*360*x - 2*Di*x)
    // sum(dist(x, Bi)^2) + sum(dist(x, Di)^2) = nCountD*360^2 + sum(Ai^2) + nCountA*x^2 + 2*360*sum(Di) - nCountD*2*360*x - 2*x*sum(Ai)
    auto SumSqrD= [&](double x, size_t nCountD, double fSumD) -> double
    {
        return x*(A.size()*x - 2*fSum) + fSumSqr + 2*360.*fSumD + nCountD*(-2*360.*x + 360.*360.);
    };
 
    // update MinAvrgAngles if lower/equal fMinSumSqrDiff found
    auto TestSum= [&](double fTestAvrg, double fTestSumDiffSqr) -> void
    {
        if (fTestSumDiffSqr < fMinSumSqrDiff)
        {
            MinAvrgVals.clear();
            MinAvrgVals.insert(CircVal<UnsignedDegRange>(fTestAvrg));
            fMinSumSqrDiff= fTestSumDiffSqr;
        }
        else if (fTestSumDiffSqr == fMinSumSqrDiff)
            MinAvrgVals.insert(CircVal<UnsignedDegRange>(fTestAvrg));
    };
 
    // ----------------------------------------------
    for (const auto& a : A)
    {
        double v= CircVal<UnsignedDegRange>(a); // convert to [0.360)
        fSum   +=     v ;
        fSumSqr+= Sqr(v);
             if (v < 180.) LowerAngles.push_back(v);
        else if (v > 180.) UpperAngles.push_back(v);
    }
 
    sort(LowerAngles.begin(), LowerAngles.end()                   ); // ascending   [  0,180)
    sort(UpperAngles.begin(), UpperAngles.end(), greater<double>()); // descending  (360,180)
 
    // ----------------------------------------------
    // start with avrg= 180, sets c,d are empty
    // ----------------------------------------------
    MinAvrgVals.clear();
    MinAvrgVals.insert(CircVal<UnsignedDegRange>(180.));
    fMinSumSqrDiff= SumSqr();
 
    // ----------------------------------------------
    // average in (180,360), set D: values in range [0,avrg-180)
    // ----------------------------------------------
    double fLowerBound= 0.; // of current sector
    double fSumD      = 0.; // of elements of set D
 
    auto iter= LowerAngles.begin();
    for (size_t d= 0; d < LowerAngles.size(); ++d)
    {
        // 1st  iteration : average in (                 180, lowerAngles[0]+180]
        // next iterations: average in (lowerAngles[i-1]+180, lowerAngles[i]+180]
        // set D          : lowerAngles[0..d]
 
        fTestAvrg= (fSum + 360.*d)/A.size(); // average for sector, that minimizes SumDiffSqr
 
        if ((fTestAvrg > fLowerBound+180.) && (fTestAvrg <= *iter+180.))  // if fTestAvrg is within sector
            TestSum(fTestAvrg, SumSqrD(fTestAvrg, d, fSumD));             // if fTestAvrg generates lower SumSqr
 
        fLowerBound= *iter      ;
        fSumD     += fLowerBound;
        ++iter;
    } 
 
    // last sector : average in [lowerAngles[lastIdx]+180, 360)
    fTestAvrg= (fSum + 360.*LowerAngles.size())/A.size(); // average for sector, that minimizes SumDiffSqr
 
    if ((fTestAvrg < 360.) && (fTestAvrg > fLowerBound))                   // if fTestAvrg is within sector
        TestSum(fTestAvrg, SumSqrD(fTestAvrg, LowerAngles.size(), fSumD)); // if fTestAvrg generates lower SumSqr
 
    // ----------------------------------------------
    // average in [0,180); set C: values in range (avrg+180, 360)
    // ----------------------------------------------
    double fUpperBound= 360.; // of current sector
    double fSumC      =   0.; // of elements of set C
 
    iter= UpperAngles.begin();
    for (size_t c= 0; c < UpperAngles.size(); ++c)
    {
        // 1st  iteration : average in [upperAngles[0]-180, 360                 )
        // next iterations: average in [upperAngles[i]-180, upperAngles[i-1]-180)
        // set C          : upperAngles[0..c]  (descendingly sorted)
 
        fTestAvrg= (fSum - 360.*c)/A.size(); // average for sector, that minimizes SumDiffSqr
 
        if ((fTestAvrg >= *iter-180.) && (fTestAvrg < fUpperBound-180.))   // if fTestAvrg is within sector
            TestSum(fTestAvrg, SumSqrC(fTestAvrg, c, fSumC));              // if fTestAvrg generates lower SumSqr
 
        fUpperBound= *iter      ; 
        fSumC     += fUpperBound; 
        ++iter;
    } 
 
    // last sector : average in [0, upperAngles[lastIdx]-180)
    fTestAvrg= (fSum - 360.*UpperAngles.size())/A.size(); // average for sector, that minimizes SumDiffSqr
 
    if ((fTestAvrg >= 0.) && (fTestAvrg < fUpperBound))                    // if fTestAvrg is within sector
        TestSum(fTestAvrg, SumSqrC(fTestAvrg, UpperAngles.size(), fSumC)); // if fTestAvrg generates lower SumSqr
 
    // ----------------------------------------------
    return MinAvrgVals;
}

Note that the complexity of this function is the complexity of the sort operation, hence O(a log a). Given sorted input, the complexity can be O(a).

Here is a sample use:

C++
vector<CircVal<UnsignedDegRange>> angles;

angles.push_back( 90.);
angles.push_back(180.);
angles.push_back(270.);

auto avg= CircAverage(angles);

27. Weighted Average of circular values

The general definition of the weighted mean (arithmetic-average) of values Circular-Values/131.png with weights Circular-Values/132.png respectively, for any mathematical field is the value x that minimizes the expression Circular-Values/133.png.
Formally: Circular-Values/134.png.

wi Is a positive real weight associated with value ai, and dist(v1,v2) is a function that measures the minimal distance between v1 and v2.

For linear values, we can calculate the average as follow:
Image 139
Let y denote the expression we want to minimize:

Circular-Values/136.png

To find the value of x that minimizes y, we need the derivative Circular-Values/137.png to be 0:

Circular-Values/138.png
Circular-Values/139.png

We got the well-known weighted arithmetic mean formula for linear values.

For circular values, we showed that
Image 144

Let y denote the expression we want to minimize:

Circular-Values/141.png

Same as for non-weighted average, we will break the input into 3 distinct multi-subsets: B, C and D.

Circular-Values/142.png

And the algorithm for weighted average is:

  • Divide the circle into sectors, such that if is x within a sector - the elements of multisets B, C and D does not change. The maximal number of sectors equals to the number of elements we want to average.
  • For each sector: Assume that Circular-Values/143.png (the value of x that minimizes y) is within this sector, calculate its value, and check that the value is indeed within this sector. If so, calculate y as follows:
    1. Assuming x=180, all values are in multiset B:

      Circular-Values/144.png

    2. Assuming x is in the range [0,180), multiset C range is (x+180,360) and multiset D is empty:

      Circular-Values/145.png

    3. Assuming x is in the range (180,360), multiset D range is [0,x-180) and multiset C is empty:

      Circular-Values/146.png

  • Find the sectors(s) that has the lowest y; return their x.

The Function WeightedCircAverage calculates the weighted average set of a given vector of circular values.
The return value is a set of all average values:

C++
// calculate weighted-average set of circular-values
// return set of average values
// T is a circular value type defined with the CircValTypeDef macro
template<typename T>
set<CircVal<T>> WeightedCircAverage(vector<pair<CircVal<T>,double>> const& A) // vector <value,weight>
{
    set<CircVal<T>>              MinAvrgVals      ; // results set
 
    // ----------------------------------------------
    // all vars: UnsignedDegRange [0,360)
    double                       fASumW        = 0.; // sum(Wi     ) of all elements of A
    double                       fASumWA       = 0.; // sum(Wi*Ai  ) of all elements of A
    double                       fASumWA2      = 0.; // sum(Wi*Ai^2) of all elements of A
    double                       fMinSumSqrDiff    ; // minimal sum of squares of differences
    vector<pair<double, double>> LowerAngles       ; // ascending   [  0,180)  <angle,weight>
    vector<pair<double, double>> UpperAngles       ; // descending  (360,180)  <angle,weight>
    double                       fTestAvrg         ;
 
    // ----------------------------------------------
    // local functions - implemented as lambdas
    // ----------------------------------------------
 
    // calc sum(Wi*dist(180, Bi)^2) - all values are in set B
    // dist(180,Bi)= |180-Bi|
    // sum(Wi*dist(x, Bi)^2) = sum(Wi*(180-Bi)^2) = sum(Wi*(180^2-2*180*Bi + Bi^2)) = 180^2*fSumW - 360*sum(Wi*Ai) + sum(Wi*Ai^2)
    auto SumSqr = [&]() -> double
    {
        return 32400.*fASumW - 360.*fASumWA + fASumWA2;
    };
 
    // calc sum(Wi*dist(x, Ai)^2). A=B+C; set D is empty
    // dist(x,Bi)= |x-Bi|
    // dist(x,Ci)= 360-(Ci-x)
    // sum(Wi*dist(x,Bi)^2)= sum(Wi*(     (x-Bi) ^2))= sum(Wi*(        Bi^2 + x^2                      - 2*Bi*x)) +
    // sum(Wi*dist(x,Ci)^2)= sum(Wi*((360-(Ci-x))^2))= sum(Wi*(360^2 + Ci^2 + x^2 - 2*360*Ci + 2*360*x - 2*Ci*x))
    //                                                 ==========================================================
    //                                                 sum(Wi*(        Ai^2 + x^2                      - 2*Ai*x))
    auto SumSqrC= [&](double x      , 
                      double fCSumW ,            // sum(Wi   ) of all elements of C
                      double fCSumWC ) -> double // sum(Wi*Ci) of all elements of C
    {
        return fASumWA2 + x*x*fASumW -2*x*fASumWA - 720*fCSumWC + (129600+720*x)*fCSumW;
    };
 
    // calc sum(Wi*dist(x, Ai)^2). A=B+D; set C is empty
    // dist(x,Bi)= |x-Bi|
    // dist(x,Di)= 360-(x-Di)
    // sum(Wi*dist(x,Bi)^2)= sum(Wi*(     (x-Bi) ^2))= sum(Wi*(        Bi^2 + x^2                      - 2*Bi*x))
    // sum(Wi*dist(x,Di)^2)= sum(Wi*((360-(x-Di))^2))= sum(Wi*(360^2 + Di^2 + x^2 + 2*360*Di - 2*360*x - 2*Di*x))
    //                                                 ==========================================================
    //                                                 sum(Wi*(        Ai^2 + x^2                      - 2*Ai*x))
    auto SumSqrD= [&](double x      ,
                      double fDSumW ,            // sum(Wi   ) of all elements of D
                      double fDSumWD ) -> double // sum(Wi*Di) of all elements of D
    {
        return fASumWA2 + x*x*fASumW -2*x*fASumWA + 720*fDSumWD + (129600-720*x)*fDSumW;
    };
 
    // update MinAvrgAngles if lower/equal fMinSumSqrDiff found
    auto TestSum= [&](double fTestAvrg, double fTestSumDiffSqr) -> void
    {
        if (fTestSumDiffSqr < fMinSumSqrDiff)
        {
            MinAvrgVals.clear();
            MinAvrgVals.insert(CircVal<UnsignedDegRange>(fTestAvrg));
            fMinSumSqrDiff= fTestSumDiffSqr;
        }
        else if (fTestSumDiffSqr == fMinSumSqrDiff)
            MinAvrgVals.insert(CircVal<UnsignedDegRange>(fTestAvrg));
    };
 
    // ----------------------------------------------
    for (const auto& a : A)
    {
        double v= CircVal<UnsignedDegRange>(a.first); // convert to [0.360)
        double w= a.second;                           // weight
        fASumW += w    ;
        fASumWA+= w*v  ;
        fASumWA2= w*v*v;
 
             if (v < 180.) LowerAngles.push_back(pair<double,double>(v,w));
        else if (v > 180.) UpperAngles.push_back(pair<double,double>(v,w));
    }
 
    sort(LowerAngles.begin(), LowerAngles.end()                                ); // ascending   [  0,180)
    sort(UpperAngles.begin(), UpperAngles.end(), greater<pair<double,double>>()); // descending  (360,180)
 
    // ----------------------------------------------
    // start with avrg= 180, sets c,d are empty
    // ----------------------------------------------
    MinAvrgVals.clear();
    MinAvrgVals.insert(CircVal<UnsignedDegRange>(180.));
    fMinSumSqrDiff= SumSqr();
 
    // ----------------------------------------------
    // average in (180,360), set D: values in range [0,avrg-180)
    // ----------------------------------------------
    double fLowerBound=   0.; // of current sector
    double fDSumW     =   0.; // sum(Wi   ) of all elements of D
    double fDSumWD    =   0.; // sum(Wi*Di) of all elements of D
 
    auto iter= LowerAngles.begin();
    for (size_t d= 0; d < LowerAngles.size(); ++d)
    {
        // 1st  iteration : average in (                 180, lowerAngles[0]+180]
        // next iterations: average in (lowerAngles[i-1]+180, lowerAngles[i]+180]
        // set D          : lowerAngles[0..d]
 
        fTestAvrg= (fASumWA + 360.*fDSumW)/fASumW; // average for sector, that minimizes SumDiffSqr
 
        if ((fTestAvrg > fLowerBound+180.) && (fTestAvrg <= (*iter).first+180.)) // if fTestAvrg is within sector
            TestSum(fTestAvrg, SumSqrD(fTestAvrg, fDSumW, fDSumWD));             // check if fTestAvrg generates lower SumSqr
 
        fLowerBound= (*iter).first                 ;
        fDSumW    += (*iter).second                ;
        fDSumWD   += (*iter).second * (*iter).first;
        ++iter;
    } 
 
    // last sector : average in [lowerAngles[lastIdx]+180, 360)
    fTestAvrg= (fASumWA + 360.*fDSumW)/fASumW; // average for sector, that minimizes SumDiffSqr
 
    if ((fTestAvrg < 360.) && (fTestAvrg > fLowerBound))                         // if fTestAvrg is within sector
        TestSum(fTestAvrg, SumSqrD(fTestAvrg, fDSumW, fDSumWD));                 // check if fTestAvrg generates lower SumSqr
 
    // ----------------------------------------------
    // average in [0,180); set C: values in range (avrg+180, 360)
    // ----------------------------------------------
    double fUpperBound= 360.; // of current sector
    double fCSumW     =   0.; // sum(Wi   ) of all elements of C
    double fCSumWC    =   0.; // sum(Wi*Ci) of all elements of C
 
    iter= UpperAngles.begin();
    for (size_t c= 0; c < UpperAngles.size(); ++c)
    {
        // 1st  iteration : average in [upperAngles[0]-180, 360                 )
        // next iterations: average in [upperAngles[i]-180, upperAngles[i-1]-180)
        // set C          : upperAngles[0..c]  (descendingly sorted)
 
        fTestAvrg= (fASumWA - 360.*fCSumW)/fASumW; // average for sector, that minimizes SumDiffSqr
 
        if ((fTestAvrg >= (*iter).first-180.) && (fTestAvrg < fUpperBound-180.)) // if fTestAvrg is within sector
            TestSum(fTestAvrg, SumSqrC(fTestAvrg, fCSumW, fCSumWC));             // check if fTestAvrg generates lower SumSqr
 
        fUpperBound= (*iter).first                 ;
        fCSumW    += (*iter).second                ;
        fCSumWC   += (*iter).second * (*iter).first;
        ++iter;
    } 
 
    // last sector : average in [0, upperAngles[lastIdx]-180)
    fTestAvrg= (fASumWA - 360.*fCSumW)/fASumW; // average for sector, that minimizes SumDiffSqr
 
    if ((fTestAvrg >= 0.) && (fTestAvrg < fUpperBound))                          // if fTestAvrg is within sector
        TestSum(fTestAvrg, SumSqrC(fTestAvrg, fCSumW, fCSumWC));                 // check if fTestAvrg generates lower SumSqr
 
    // ----------------------------------------------
    return MinAvrgVals;
}
Note that the complexity of this function is the complexity of the sort operation, hence O(a log a). Given sorted input, the complexity can be O(a). Here is a sample use:
C++
vector<pair<CircVal<UnsignedDegRange>,double>> angles;

angles.push_back(make_pair( 90., 0.3));
angles.push_back(make_pair(180., 0.5));
angles.push_back(make_pair(270., 0.7));

auto avg= WeightedCircAverage(angles);

28. Median of circular values

Definition 7 

Given the multiset Circular-Values/148.png:
Let sequence Circular-Values/149.png be the non-decreasingly sorted permutation of  Circular-Values/150.png

 

For linear values: 

If n is odd:
Image 154 
If n is even: 
Image 155

 

For circular values:  

Since there is no "first" nor "last" element, we'll evaluate all rotations of S:

For each k∈{1…n}:

If n is odd: 
Image 156 
If n is even:
Image 157  (The avrg set may contain more than one element) 

The set of all candidate medians:
Image 158

Since an optimality property of the median is the minimization of the sum of absolute distances: 

Image 159 (The median set may contain more than one element)

 

Here are some examples for circular values, given the range [0,360): 

Image 160

Image 161

Image 162

Image 163

 

The function CircMedian calculates the median set of a given vector of circular values.
The return value is a set of all median values: 

C++
template<typename T>
set<CircVal<T>> CircMedian(vector<CircVal<T>> const& A)
{
    set <CircVal<T>> X;           // results set
 
    // ----------------------------------------------
    set<CircVal<T>> B;
    if (A.size() % 2 == 0)        // even number of values
    {
        auto S= A;
        sort(S.begin(), S.end()); // A, sorted
 
        for (size_t m= 0; m < S.size(); ++m)
        {
            size_t n= m+1; if (n==S.size()) n= 0;
            double d= CircVal<T>::Sdist(S[m],S[n]);
 
            // insert average set of each two circular-consecutive values
            B.insert(((double)S[m] + d/2.));
            if (d == -CircVal<T>::GetR()/2.)
                B.insert(((double)S[n] + d/2.));
        }
    }
    else                          // odd number of values
        for (size_t m= 0; m < A.size(); ++m)
            B.insert(A[m]);       // convert vector to set - remove duplicates
 
    // ----------------------------------------------
    double fMinSum= numeric_limits<double>::max();
 
    for (const auto& b : B)
    {
        size_t nL  = 0 ; // number of elements in A for which Sdist(Ai,b)<0
        size_t nH  = 0 ; // number of elements in A for which Sdist(Ai,b)>0
        double fSum= 0.; // sum(Sdist(a,b))
 
        for (const auto& a : A)
        {
            double d= CircVal<T>::Sdist(b,a);
            if (d == 0.) continue       ;
            if (d <  0.) ++nL; else ++nH;
            fSum+= abs(d);
        }
 
        if ((nL<= A.size()/2.) && (nH<= A.size()/2.))
                 if (fSum==fMinSum)              X.insert(b);
            else if (fSum< fMinSum) { X.clear(); X.insert(b); fMinSum= fSum; }
    }
 
    // ----------------------------------------------
    return X;
}

Here is a sample use:

C++
vector<CircVal<UnsignedDegRange>> angles;
 
angles.push_back( 90.);
angles.push_back(180.);
angles.push_back(270.);
 
auto mdn= CircMedian(angles);

29. Circular parameter estimation based on noisy measurements

For linear values:

For many observations/measurements techniques/instruments, normal (Gaussian) distribution is often used as a first approximation to model observations/measurements error. As Jonas Ferdinand Gabriel Lippmann said: “Everyone believes in the [normal] law of errors: the mathematicians, because they think it is an experimental fact; and the experimenters, because they suppose it is a theorem of mathematics.” Hence, for our reference problem:

  • Given a multiset of independent distance measurements from a stationary transmitter to a stationary receiver, using a measurement technique with a normal-distributed error – Estimate the distance.

It can be shown that given a multiset of independent measurements of a constant parameter, were the measurement technique/instrument has a normal distributed error, the maximum likelihood estimation (MLE) of the measured parameters equals to the arithmetic average of the measurements.

Hence, for our reference problem, the maximum likelihood estimation of the distance is the arithmetic average of the measurements.

For circular value:

Here, wrapped normal distribution, or wrapped truncated normal distribution may be used to model a measurement error (according to the physical phenomena we want to model). The maximal measurement error is bounded to half circle.

Given a multiset of measurements of a constant circular parameter, were the measurement technique/instrument has a wrapped/wrapped-truncated normal distributed error.

Here is our reference problem:

  • Given a multiset of independent direction measurements from a stationary transmitter to a stationary receiver, using a measurement technique with a wrapped/wrapped-truncated normal-distributed error – Estimate the direction.
We will consider two methods to estimate the circular parameter:
  1. The circular average of the samples - as described above
  2. Consider each sample as a 2D unit vector, and calculate the weighted vector Image 164 - The suggested method in most textbooks

Simulation: As a test case, we will generate 50,000 trails. Each trail will contain 1000 samples of measurements with wrapped/wrapped-truncated normal distributed error (90° truncation span). We will use the two methods to calculate the estimate for each trail, and then calculate the RMS error for each method.

We will repeat the test for different standard deviations – between 1° and 100°.

The following code will collect the data for the two methods:

C++
std::default_random_engine rand_engine;
std::random_device         rnd_device ;
rand_engine.seed(rnd_device()); // reseed engine
 
ofstream f1("log1.txt");
 
Concurrency::parallel_for(1, 101, [&](int nStdDev) // for each value of standard-deviation
{
    uniform_real<double> ud(0., 360.);
 
    double fSumSqrErr1= 0.;
    double fSumSqrErr2= 0.;
 
    const size_t nTrails = 50000;            // number of trails
    const size_t nSamples=  1000;            // number of observations per trail
            
    vector<CircVal<UnsignedDegRange>> vInput(nSamples);
 
    const double fAvrg= ud(rand_engine);     // our const parameter for this trail
    wrapped_normal_distribution          <double> r_wnd1(fAvrg, nStdDev,                       0., 360.);
 // wrapped_truncated_normal_distribution<double> r_wnd1(fAvrg, nStdDev, fAvrg-45., fAvrg+45., 0., 360.);
 
    for (size_t t= 0; t<nTrails; ++t)
    {
        for (auto& Sample : vInput)
            Sample= r_wnd1(rand_engine); // generate "noisy" observation
 
        set<CircVal<UnsignedDegRange>> sAvrg1= CircAverage(vInput);                  // avrg - method 1 (new method)
 
        double fSigSin= 0.;
        double fSigCos= 0.;
 
        for (const auto& Sample : vInput)
        {
            fSigSin+= sin(Sample);
            fSigCos+= cos(Sample);
        }
 
        CircVal<UnsignedDegRange> Avrg2= atan2<UnsignedDegRange>(fSigSin, fSigCos); // avrg - method 2 (conventional method)
 
        const double fErr1= CircVal<UnsignedDegRange>::Sdist(*sAvrg1.begin(), fAvrg); // error of estimate - method 1
        const double fErr2= CircVal<UnsignedDegRange>::Sdist(Avrg2          , fAvrg); // error of estimate - method 2
 
        fSumSqrErr1+= Sqr(fErr1);
        fSumSqrErr2+= Sqr(fErr2);
    }
 
    const double fRMS1= sqrt(fSumSqrErr1/ (nTrails-1)); // root mean square error - method 1
    const double fRMS2= sqrt(fSumSqrErr2/ (nTrails-1)); // root mean square error - method 2
 
    f1 << nStdDev << "\t" << fRMS1 << "\t" << fRMS2 << endl; // save RMS results to file
} );

Note that parallel_for is used for parallelizing tests for different standard deviations.

So which method is better?

Circular-Values/158.png

The graph above shows the RMS error of the average estimation for 50000 trails. Each trail averages 1000 measurements with wrapped normal distributed error with standard deviation according to the X axis.

Method 1 turns out to give more accurate estimations when the standard deviation is less than 65°-80° (this value grows when there are more samples per trail: 65° for 10 samples per trail, and 78° for 1000 samples per trail). If the data is noisier - method 2 becomes better. Since the standard deviation in most real scenarios is much less than 65°, method 1 is better for almost any practical purpose.

Circular-Values/159.png

The graph above shows the RMS error of the average estimation for 50000 trails. Each trail averages 1000 measurements with wrapped truncated normal distributed error with standard deviation according to the X axis.
The truncation span is 90° around the parameter value, and the wrapping range is [0,360).

Method 1 proves to be better for wrapped truncated normal standard deviation error, for any value of standard deviation.

30. Interpolation and average estimation of sampled continuous-time circular signal

Given a sequence of samples Circular-Values/160.png of a continuous-time circular signal, sampled at times Circular-Values/161.png respectively, where the measurement technique/instrument is accurate - Estimate the mean value of the signal at period Circular-Values/162.png.

Linear example:

  • Given airplane velocity measured each second [MPH] over a period of one hour – Estimate the mean velocity
Circular example:
  • Given airplane heading measured each second [deg] over a period of one hour - Estimate the mean heading

For circular values, it is not obvious in which “direction” the signal changed between two consecutive samples. For example, if one sample was 10°, and the next sample was 300°– did the signal changed in the increasing direction, or in the decreasing direction?

Since we can’t tell, our only way is to sample fast enough to ensure that the difference between consecutive measurements will always be less than 180°, so we can always be sure about the direction in which the signal changed. If the measure-rate is not high enough, a 330° change in the measured property would be incorrectly interpreted as a 30° change in the other direction, and will insert an error into the calculation.

Basically, in order to estimate the average, we should reconstruct the signal (interpolation), calculate the area bounded between the signal and the t-axis (integration), and divide it by Circular-Values/163.png.

Linear Interpolation

The simplest reconstruction technique is linear interpolation (connecting every two consecutive samples by a straight line). We will discuss the circular values equivalent to linear values linear interpolation.

For linear values, the average value of a signal, based on linear interpolation, is the weighted average of

Image 171
where, as defined previouslyImage 172,

For circular values, it is the weighted circular average of

Circular-Values/165.png

For the above example, assume that the airplane heading measurements vector is {200, 240, 320, 40, 100, 20, 300}.

Circular-Values/166.png

Here is how to calculate the average:

  • Connect the samples such that the difference between consecutive samples will always be less than 180° (see the red line intervals in the figure above).
  • Calculate the average circular value, and the weight of each interval k in the range [2,n]:Image 175
    In our examples the intervals average values are {220, 280, 0, 70, 60, 340}
  • Calculate the weighted circular average of all intervals average values
    (non-weighted circular average may be used when the samples are equally spaced).
    In our examples, assuming equally spaced sampling, it is 332.857.
Note that this method is more accurate than the Mitsuta method, described in the U.S. Meteorological Monitoring Guidance for Regulatory Modeling Applications (http://www.epa.gov/scram001/guidance/met/mmgrma.pdf), which gives equal weight to all measurements (For interpolation – the 1st and last measurements should have only half-weight) The Mitsuta method is also limited to equally-spaced samples.

The following code calculates the average estimation based on samples:

C++
// estimate the average of a sampled continuous-time circular signal, using circular linear interpolation
// T is a circular value type defined with the CircValTypeDef macro
template<typename T>
class CAvrgSampledCircSignal
{
    size_t                           m_nSamples ;
    CircVal<T>                       m_fPrevVal ; // previous value
    double                           m_fPrevTime; // previous time
    vector<pair<CircVal<T>, double>> m_Intervals; // vector of (avrg,weight) for each interval
 
public:
    CAvrgSampledCircSignal()
    {
        m_nSamples= 0;
    }
 
    void AddMeasurement(CircVal<T> fVal, double fTime)
    {
        if (m_nSamples)
        {
            assert(fTime > m_fPrevTime);
 
            double fIntervalAvrg  = CircVal<T>::Wrap((double)m_fPrevVal     + 
                                    CircVal<T>::Sdist(m_fPrevVal, fVal)/2.)   ;
            double fIntervalWeight= fTime-m_fPrevTime                         ;
            m_Intervals.push_back(make_pair(fIntervalAvrg, fIntervalWeight));
        }
 
        m_fPrevVal = fVal ;
        m_fPrevTime= fTime;
        ++m_nSamples;
    }
 
    // calculate the weighted average for all intervals
    bool GetAvrg(CircVal<T>& fAvrg)
    {
        switch (m_nSamples)
        {
        case 0:
            fAvrg= CircVal<T>::GetZ();
            return false;
 
        case 1:
            fAvrg= m_fPrevVal;
            return true;
 
        default:
            fAvrg= *WeightedCircAverage(m_Intervals).begin();
            return true;
        }
    }
};

Here is a sample use:

C++
CAvrgSampledCircSignal<UnsignedDegRange> A1;
A1.AddMeasurement(CircVal<UnsignedDegRange>(200.), 1);
A1.AddMeasurement(CircVal<UnsignedDegRange>(300.), 2);
A1.AddMeasurement(CircVal<UnsignedDegRange>( 20.), 6);
 
CircVal<UnsignedDegRange> ad1;
A1.GetAvrg(ad1);   

License

This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)


Written By
Israel Israel
I consider software development a craft that mixes engineering, science, and art, and is an utmost form of self-expression.

Comparing software languages to musical instruments, C is the harpsichord, C++ is the piano, and C# and java are electronic synthesizers. I would choose the piano!

Enjoying a grand piano (C++ 11/14/17), I am playing Classical music (that's mathematics and algorithms). Sometimes Jazz. When needed, I'll play Big band or even Acid. I'll leave Rave for others.

I can be reached at koganlior1 at gmail dot com

Comments and Discussions

 
GeneralMy vote of 5 Pin
Md. Marufuzzaman15-Jun-11 19:58
professionalMd. Marufuzzaman15-Jun-11 19:58 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Praise Praise    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.