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 measured on a 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 
is used, sometimes 
. Sometimes it is 
and sometimes 
- even in the same application. For time-of-day computations, the range may be, or something else. We need a scheme to perform operations between values in different representations.
2. Features
In this article I am going to present:
- Definition of circular value type and its associated operators
- A C++0x infrastructure for doing mathematics with circular values
- C++0x 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++ 2010.
The given code relies heavily of C++0x features. It also uses some Microsoft-specific syntax.
The code may be easily converted to any other C++0x 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:
template <typename T>
T Sqr(const T& x)
{
return x*x;
}
- Modulo function:
and 
are linear values, and the result is also a linear value.
Floating-point modulo operator can be defined in several ways: The sign of the result may be equal to the sign of the dividend, equal to the sign of the divisor, always non-negative or the result may be closest to zero.
We will use the following definition.
Definition 1
The quotient rounds towards negative infinity and the remainder has the same sign as the divisor.
When 
, the function may be undefined, or defined, or defined as (It is not relevant for our usage).
Note, that the given definition is not similar to fmod: In C, the result has the same sign as the dividend; In C++, according to the 2003 standard, if one of the operands is negative - the sign of the remainder is implementation-defined.
This definition was given by Knuth in The Art of Computer Programming, Vol.1
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);
}
Theorem 1
Proof
- 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:
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);
}
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 
is defined by three constants:
- The range
- The zero-value
.
Note the use of a right-open interval range.
For example, for angles in the range 
, it is natural to define the zero-value is in the middle, while for angles in the range 
it is natural to define zero-value on the edge. 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.
We will use the following macro to define a circular value type:
#define CircValTypeDef(_Name, _L, _H, _Z) \
struct _Name \
{ \
static const double L ; \
static const double H ; \
static const double Z ; \
static const double R ; \
static const double R_2; \
\
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:
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:
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.
template <typename Type>
class CircVal
{
…
}
Here is an example of defining a circular value variable with a given type:
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.
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 , 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
The CircVal::Wrap static function implements it:
static double Wrap(double r)
{
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;
return Mod(r - Type::L, Type::R) + Type::L;
}
8. A Walk
Definition 4
For linear values, a walk 
is defined by:
- A start point
- The walk length
A walk is a movement along the linear axis, from the start point, in the given direction, with a given length.
When the length is positive, we call it ‘an increasing walk’ (toward positive infinity).
When the length is negative, we call it ‘a decreasing walk’ (toward negative infinity).
The end point 
is the point reached at the end of the walk. We’ll call a walk with a start point and an end point “A walk from to”.
There is one possible walk from s to e: <s,,e — s >
For any circular value type, a walk is defined by:
- A start point
- The walk length
A walk is a movement along the circular axis, from the start point, in the given direction, with a given length.
When the length is positive, we call it ‘an increasing walk’ (clockwise in the picture above).
When the length is negative, we call it ‘a decreasing walk’ (counterclockwise in the picture above).
The end point 
is the point reached at the end of the walk.
We’ll call a walk with a start point and an end point “A walk from to”.
Since the walk is wrapped around the circle:
There is infinite number of walks from 
:
9. Distance between two values
Definition 5
For every two linear values 
we can define a signed distance function from 
:
Which is the walk length from to, such that:
For example, 
For every two circular values 
we can define a signed distance function from 
to:
Which is the length of the walk from to, with the lowest absolute-value length (The walk may be increasing or decreasing).
It can easily be deduced that:
Since
And since
We can write:
which we will use for faster implementation.
The CircVal::Sdist static function calculates sdist():
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 ;
}
We'll define another distance function for circular values:
- The length of the increasing-walk from c1 to c2 with the lowest length. It is always in the range 
.
There is no similar distance function for linear values.
Definition 6
which may also be expressed as:
It can easily be deduced that:
The CircVal::Pdist static function calculates:
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
Proof
Rephrasing the definitions:
Theorem 3
(Very intuitive)
Proof
According to the definition above:
And also
, which are equal in all cases.
Theorem 4
For any two circular values: 
Proof
Theorem 5
For any two circular values: 
Proof
Theorem 6
For any two circular values: 
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:
CircVal<UnsignedDegRange> d1= 10.;
CircVal<SignedDegRange > d2= -10.;
CircVal<SignedRadRange > d3 ;
d3= d1+d2;
To convert a circular value 
from circular type 
to circular type 
, we’ll retain the arc-length between and 
. For that, we can use either of these formulae:
Theorem 7
I and II are equivalent
Proof
Therefore, in both cases I =II
The conversion can be accomplished easily be implementing appropriate constructor and assignment operator:
template<typename CircVal2>
CircVal(const CircVal2& c2)
{
double val2= c2.Pdist(c2.GetZ(), c2);
val= Wrap(val2*Type::R/c2.GetR() + Type::Z);
}
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 circular object is implemented with operator double() const overloading.
CircVal(double r)
{
val= Wrap(r);
}
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 is defined for any circular value. The result is a circular value.
Negation of a circular value is the symmetric value relative to (the zero-value).
- Opposite operator is defined for any circular value. The result is a circular value of the class.
Opposite of a circular value is the circular value pointing to the opposite direction
- Addition and subtraction operators
are defined for any two circular values.
The result is a circular value of the class - 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).
- Multiplication operator (c·r) is defined for any circular value c and any linear value r.
The result is a circular value of the class - the point reached and the end of the walk:
Starting from z, walking distance r·pdist(c,z).
- Division operator (c/r) is defined for any circular value c and any non-zero linear value r.
The result is a circular value of the class - the point reached and the end of the walk:
Starting from z, walking 
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 < or = 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:
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:
15. Comparison operators
Defining operators = and
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:
All these requirements, however, tell us nothing about the meaning of the > operator.
For linear values we can define:
Based on this, we can define the circular values operator in different ways:
<>Hence, for two circular values
c1,
c2, for each comparison operator
, we can use:
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 default implementation is (i), but it may be different for different applications.
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; }
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
Theorem 8
The unary operator can be equivalently defined by each of these two equivalent formulae:
The and the binary operators can be equivalently defined by each of these two equivalent formulae:
Proof
Here is the implementation of these operators:
template <typename Type>
class CircVal
{
double val;
public:
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)); } const CircVal operator~ ( ) const { return Wrap(val+Type::R_2 ); }
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:
template <typename Type>
class CircValTester
{
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);
}
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++0x random generation, which is a very useful addition to the language:
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. ); std::uniform_real_distribution<double> t_uni_dist(-1. , 1. );
std::random_device rnd_device;
rand_engine.seed(rnd_device());
for (unsigned i= 10000; i--;)
{
CircVal<Type> c1(c_uni_dist(rand_engine)); CircVal<Type> c2(c_uni_dist(rand_engine)); CircVal<Type> c3(c_uni_dist(rand_engine)); double r (r_uni_dist(rand_engine)); double a1(t_uni_dist(rand_engine)); double a2(t_uni_dist(rand_engine));
assert (c1 == CircVal<Type>((double)c1) );
AssertCircAlmostEq(+c1 , c1 ); AssertCircAlmostEq(-(-c1) , c1 ); AssertCircAlmostEq(c1 + c2 , c2 + c1 ); AssertCircAlmostEq(c1 + (c2 +c3) , (c1 + c2) + c3 ); AssertCircAlmostEq(c1 + -c1 , ZeroVal ); AssertCircAlmostEq(c1 + ZeroVal , c1 );
AssertCircAlmostEq(c1 - c1 , ZeroVal ); AssertCircAlmostEq(c1 - ZeroVal , c1 ); AssertCircAlmostEq(ZeroVal - c1 , -c1 ); AssertCircAlmostEq(c1 - c2 , -(c2 - c1) );
AssertCircAlmostEq(c1 * 0. , ZeroVal ); AssertCircAlmostEq(c1 * 1. , c1 ); AssertCircAlmostEq(c1 / 1. , c1 );
AssertCircAlmostEq((c1 * (1./(r+1.))) / (1./(r+1.)) , c1 ); AssertCircAlmostEq((c1 / ( r+1.) ) * ( r+1. ) , c1 );
AssertCircAlmostEq(~(~c1) , c1 ); AssertCircAlmostEq(c1 - (~c1) , ToC<Type>(Type::R/2.) );
AssertAlmostEq (sin(ToR(CircVal<SignedRadRange>(c1))), sin(c1) ); AssertAlmostEq (cos(ToR(CircVal<SignedRadRange>(c1))), cos(c1) ); AssertAlmostEq (tan(ToR(CircVal<SignedRadRange>(c1))), tan(c1) );
AssertAlmostEq (sin(-c1) , -sin(c1) ); AssertAlmostEq (cos(-c1) , cos(c1) ); AssertAlmostEq (tan(-c1) , -tan(c1) );
AssertAlmostEq (sin(c1+ToC<Type>(Type::R/4.)) , cos(c1) ); AssertAlmostEq (cos(c1+ToC<Type>(Type::R/4.)) , -sin(c1) ); AssertAlmostEq (sin(c1+ToC<Type>(Type::R/2.)) , -sin(c1) ); AssertAlmostEq (cos(c1+ToC<Type>(Type::R/2.)) , -cos(c1) );
AssertAlmostEq (Sqr(sin(c1))+Sqr(cos(c1)) , 1. );
AssertAlmostEq (sin(c1)/cos(c1) , tan(c1) );
AssertCircAlmostEq(asin<Type>(a1) , CircVal<SignedRadRange>(asin(a1))); AssertCircAlmostEq(acos<Type>(a1) , CircVal<SignedRadRange>(acos(a1))); AssertCircAlmostEq(atan<Type>(a2) , CircVal<SignedRadRange>(atan(a2)));
AssertCircAlmostEq(asin<Type>(a1) + asin<Type>(-a1) , ZeroVal ); AssertCircAlmostEq(acos<Type>(a1) + acos<Type>(-a1) , ToC<Type>(Type::R/2.) ); AssertCircAlmostEq(asin<Type>(a1) + acos<Type>( a1) , ToC<Type>(Type::R/4.) ); AssertCircAlmostEq(atan<Type>(a2) + atan<Type>(-a2) , ZeroVal );
assert (c1 > c2 == (c2 < c1) ); assert (c1 >= c2 == (c2 <= c1) ); assert (c1 >= c2 == ( (c1 > c2) || (c1 == c2)) ); assert (c1 > c2 == (!(c1 == c2) && !(c1 < c2)) ); assert (c1 == c2 == (!(c1 > c2) && !(c1 < c2)) ); assert (c1 < c2 == (!(c1 == c2) && !(c1 > c2)) ); assert (c1 < c2 == (!(c1 == c2) && !(c1 > c2)) ); assert (!(c1>c2) || !(c2>c3) || (c1>c3) );
AssertCircAlmostEq(c1 , ToC<Type>(ToR( c1) ) ); AssertCircAlmostEq(-c1 , ToC<Type>(ToR(-c1) ) ); AssertCircAlmostEq(c1 + c2 , ToC<Type>(ToR(c1)+ToR(c2)) ); AssertCircAlmostEq(c1 - c2 , ToC<Type>(ToR(c1)-ToR(c2)) ); AssertCircAlmostEq(c1 * r , ToC<Type>(ToR(c1)*r ) ); AssertCircAlmostEq(c1 / r , ToC<Type>(ToR(c1)/r ) );
}
}
public:
CircValTester()
{
Test();
}
};
The code below executes these tests for all 8 circular-value defined types:
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:
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++0x 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 
.
Therefore, the probability density function 
of the wrapped normal distribution is:
The PDF is described in the following image, for same 
, but different values of 
:
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:
#include "WrappedNormalDist.h"
std::default_random_engine rand_engine;
std::random_device rnd_device ;
rand_engine.seed(rnd_device());
double fAvrg = 0;
double fSigma= 45;
double fL = -180; double fH = 180;
wrapped_normal_distribution<double> r_wrp(fAvrg, fSigma, fL, fH);
double r1= r_wrp(rand_engine);
The implementation is based on VC 2010 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 
. The PDF is described in the following image, for different values of 
:
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:
#include "TruncNormalDist.h"
std::default_random_engine rand_engine;
std::random_device rnd_device ;
rand_engine.seed(rnd_device());
double fAvrg = 0;
double fSigma= 45;
double fA = -40; double fB = 40;
truncated_normal_distribution<double> r_trn(fAvrg, fSigma, fA, fB);
double r2= r_trn(rand_engine);
The implementation is based on VC 2010 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 , and then wrapped around the range
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, 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 hours.
After wrapping, it is possible that - See picture below, which may depict human estimation of an azimuth between two points, where the real azimuth is 0°.
Wrapped truncated normal distributed values can be generated as follows:
#include "WrappedTruncNormalDist.h"
std::default_random_engine rand_engine;
std::random_device rnd_device ;
rand_engine.seed(rnd_device());
double fAvrg = 0;
double fSigma= 100;
double fA = -500; double fB = 500; double fL = 0; double fH = 360;
wrapped_truncated_normal_distribution<double> r_ wrp_trn(fAvrg, fSigma, fA, fB, fL, fH);
double d= r_wrp_trn(rand_engine);
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 type of problem, I’ll give a description, and a 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 of a continuous-time [circular] signal, sampled at times respectively, where the measurement technique/instrument is accurate - Estimate the mean value of the signal at period]
Usage: Estimating the average of a continuous-signal, in a given period, based on constant-interval 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, , we can start at and take a walk with length 
. The end-point of the walk is the average. We can, by symmetry, start at and take a walk with length .
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:
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 the 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 
for any mathematical field, is the value of x that minimizes 
. Formally: 
, where dist(v1,v2) is a functions that measures the minimal distance between v1 and v2.
For linear values, we can calculate the average as follow: 
Let y denote the expression we want to minimize.
To find the value of x that minimizes y, we need the derivative 
to be 0:
for 
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 
, 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
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.
The values of y that minimizes is our average.
Here are some examples:
Example 1: For the input {0,30,60,90}, we’ll plot 
as a function of x:
We can see that the average is {45} which matches out intuition.
Example 2: For the input {0,90,180,270}, we’ll plot 
as a function of x:
In this case the average is {45.135, 225, 315}
(For all these values, 
is equal).
Example 3: For the input, we’ll plot as a function of:
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:
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 each (auto a in Angles2)
fSum+= Sqr(__min(abs(x-a), 360.-abs(x)));
f0 << x << "\t" << fSum << endl;
}
To find the values of that minimizes, we need the derivative 
to be 0. However, due to the min and abs operators, the expression 
is not derivable.
Let us try to rephrase the expression.
We’ll divide the multiset into 3 distinct multi-subsets:
We’ll denote with a, b, c, d – the number of elements in A, B, C, D respectively.
Of course 
The minimum of y may be where 
, or where 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 
. 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:
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
(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, calculateas follows:
- Assuming x=180, all values are in multiset B:
- Assuming is in the range, multiset range is and multiset is empty:
- Assuming is in the range), multiset range is and multiset is empty:
- Find the sectors(s) that has the lowest ; return their
The Function CircAverage calculates the average set of a given vector of circular values.
The return value is a set of all average values:
template<typename T>
set<CircVal<T>> CircAverage(vector<CircVal<T>> const& A)
{
set<CircVal<T>> MinAvrgVals ;
double fSum = 0.; double fSumSqr = 0.; double fMinSumSqrDiff ; vector<double> LowerAngles ; vector<double> UpperAngles ; double fTestAvrg ;
auto SumSqr = [&]() -> double
{
return 32400.*A.size() - 360.*fSum + fSumSqr;
};
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.);
};
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.);
};
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 each (auto a in A)
{
double v= CircVal<UnsignedDegRange>(a); 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() ); sort(UpperAngles.begin(), UpperAngles.end(), greater<double>());
MinAvrgVals.clear();
MinAvrgVals.insert(CircVal<UnsignedDegRange>(180.));
fMinSumSqrDiff= SumSqr();
double fLowerBound= 0.; double fSumD = 0.;
auto iter= LowerAngles.begin();
for (unsigned d= 0; d < LowerAngles.size(); ++d)
{
fTestAvrg= (fSum + 360.*d)/A.size();
if ((fTestAvrg > fLowerBound+180.) && (fTestAvrg <= *iter+180.)) TestSum(fTestAvrg, SumSqrD(fTestAvrg, d, fSumD));
fLowerBound= *iter ;
fSumD += fLowerBound;
++iter;
}
fTestAvrg= (fSum + 360.*LowerAngles.size())/A.size();
if ((fTestAvrg < 360.) && (fTestAvrg > fLowerBound)) TestSum(fTestAvrg, SumSqrD(fTestAvrg, LowerAngles.size(), fSumD));
double fUpperBound= 360.; double fSumC = 0.;
iter= UpperAngles.begin();
for (unsigned c= 0; c < UpperAngles.size(); ++c)
{
fTestAvrg= (fSum - 360.*c)/A.size();
if ((fTestAvrg >= *iter-180.) && (fTestAvrg < fUpperBound-180.)) TestSum(fTestAvrg, SumSqrC(fTestAvrg, c, fSumC));
fUpperBound= *iter ;
fSumC += fUpperBound;
++iter;
}
fTestAvrg= (fSum - 360.*UpperAngles.size())/A.size();
if ((fTestAvrg >= 0.) && (fTestAvrg < fUpperBound)) TestSum(fTestAvrg, SumSqrC(fTestAvrg, UpperAngles.size(), fSumC));
return MinAvrgVals;
}
Note that the complexity of this function is the complexity of the sort operation, hence. Given sorted input, the complexity can be.
Here is a sample use:
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 
with weights 
respectively, for any mathematical field is the value x that minimizes the expression 
. Formally: 
.
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:
Let y denote the expression we want to minimize:
To find the value of x that minimizes y, we need the derivative 
to be 0:
We got the well-known weighted arithmetic mean formula for linear values.
For circular values, we showed that
Let y denote the expression we want to minimize:
Same as for non-weighted average, we will break the input into 3 distinct multi-subsets: B, C and D.
And the algorithm for weighted average is:
- 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 we want to average.
- For each sector: Assume that
(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, calculateas follows:
- Assuming x=180, all values are in multiset B:
- Assuming is in the range, multiset range is and multiset is empty:
- Assuming is in the range), multiset range is and multiset is empty:
- Find the sectors(s) that has the lowest ; return their
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:
template<typename T>
set<CircVal<T>> WeightedCircAverage(vector<pair<CircVal<T>,double>> const& A) {
set<CircVal<T>> MinAvrgVals ;
double fASumW = 0.; double fASumWA = 0.; double fASumWA2 = 0.; double fMinSumSqrDiff ; vector<pair<double, double>> LowerAngles ; vector<pair<double, double>> UpperAngles ; double fTestAvrg ;
auto SumSqr = [&]() -> double
{
return 32400.*fASumW - 360.*fASumWA + fASumWA2;
};
auto SumSqrC= [&](double x ,
double fCSumW , double fCSumWC ) -> double {
return fASumWA2 + x*x*fASumW -2*x*fASumWA - 720*fCSumWC + (129600+720*x)*fCSumW;
};
auto SumSqrD= [&](double x ,
double fDSumW , double fDSumWD ) -> double {
return fASumWA2 + x*x*fASumW -2*x*fASumWA + 720*fDSumWD + (129600-720*x)*fDSumW;
};
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 each (auto a in A)
{
double v= CircVal<UnsignedDegRange>(a.first); double w= a.second; 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() ); sort(UpperAngles.begin(), UpperAngles.end(), greater<pair<double,double>>());
MinAvrgVals.clear();
MinAvrgVals.insert(CircVal<UnsignedDegRange>(180.));
fMinSumSqrDiff= SumSqr();
double fLowerBound= 0.; double fDSumW = 0.; double fDSumWD = 0.;
auto iter= LowerAngles.begin();
for (unsigned d= 0; d < LowerAngles.size(); ++d)
{
fTestAvrg= (fASumWA + 360.*fDSumW)/fASumW;
if ((fTestAvrg > fLowerBound+180.) && (fTestAvrg <= (*iter).first+180.)) TestSum(fTestAvrg, SumSqrD(fTestAvrg, fDSumW, fDSumWD));
fLowerBound= (*iter).first ;
fDSumW += (*iter).second ;
fDSumWD += (*iter).second * (*iter).first;
++iter;
}
fTestAvrg= (fASumWA + 360.*fDSumW)/fASumW;
if ((fTestAvrg < 360.) && (fTestAvrg > fLowerBound)) TestSum(fTestAvrg, SumSqrD(fTestAvrg, fDSumW, fDSumWD));
double fUpperBound= 360.; double fCSumW = 0.; double fCSumWC = 0.;
iter= UpperAngles.begin();
for (unsigned c= 0; c < UpperAngles.size(); ++c)
{
fTestAvrg= (fASumWA - 360.*fCSumW)/fASumW;
if ((fTestAvrg >= (*iter).first-180.) && (fTestAvrg < fUpperBound-180.)) TestSum(fTestAvrg, SumSqrC(fTestAvrg, fCSumW, fCSumWC));
fUpperBound= (*iter).first ;
fCSumW += (*iter).second ;
fCSumWC += (*iter).second * (*iter).first;
++iter;
}
fTestAvrg= (fASumWA - 360.*fCSumW)/fASumW;
if ((fTestAvrg >= 0.) && (fTestAvrg < fUpperBound)) TestSum(fTestAvrg, SumSqrC(fTestAvrg, fCSumW, fCSumWC));
return MinAvrgVals;
}
Note that the complexity of this function is the complexity of the sort operation, hence. Given sorted input, the complexity can be.
Here is a sample use:
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
We’ll use the following notation: Given multiset M,
Definition 7
For both linear and circular values: Given the multiset 
:
If |A| is odd:
We’ll define B as the set of all elements in A. (Note that B is set and A is multiset)
If |A| is even:
Let sequence 
be the non-decreasingly sorted permutation of 
For linear values: Let set 
For circular values: Let set 
(note that the average set of two elements may contain two elements)
For linear values, the median set will always contains one element; therefore the median is commonly denoted as a value, and not as a set of values. Here are some examples for linear values:
Therefore, Median(A) = {3.5}
Here are some examples for circular values, given the range [0,360):
Therefore, Median(A) = {200}
Therefore, Median(A) = {70}
Therefore, Median(A) = {45,135,225,315}
Therefore, Median(A) = {270,270.5}
The function CircMedian calculates the median set of a given vector of circular values.
The return value is a set of all median values:
template<typename T>
set<CircVal<T>> CircMedian(vector<CircVal<T>> const& A)
{
set <CircVal<T>> X;
set<CircVal<T>> B;
if (A.size() % 2 == 0) {
auto S= A;
sort(S.begin(), S.end());
for (unsigned m= 0; m < S.size(); ++m)
{
int n= m+1; if (n==S.size()) n= 0;
double d= CircVal<T>::Sdist(S[m],S[n]);
B.insert(((double)S[m] + d/2.));
if (d == -CircVal<T>::GetR()/2.)
B.insert(((double)S[n] + d/2.));
}
}
else for (unsigned m= 0; m < A.size(); ++m)
B.insert(A[m]);
double fMinSum= numeric_limits<double>::max();
for each (auto b in B)
{
size_t nL = 0 ; size_t nH = 0 ; double fSum= 0.;
for each (auto a in 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:
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:
- The circular average of the samples - as described above
- Consider each sample as a 2D unit vector, and calculate the weighted vector . 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:
std::default_random_engine rand_engine;
std::random_device rnd_device ;
rand_engine.seed(rnd_device());
ofstream f1("log1.txt");
Concurrency::parallel_for(1, 101, [&](int nStdDev) {
uniform_real<double> ud(0., 360.);
double fSumSqrErr1= 0.;
double fSumSqrErr2= 0.;
const int nTrails = 50000; const int nSamples= 1000;
vector<CircVal<UnsignedDegRange>> vInput(nSamples);
double fAvrg= ud(rand_engine); wrapped_normal_distribution <double> r_wnd1(fAvrg, nStdDev, 0., 360.);
for (int t= 0; t<nTrails; t++)
{
for (int i= 0; i<nSamples; i++)
vInput[i]= r_wnd1(rand_engine);
set<CircVal<UnsignedDegRange>> sAvrg1= CircAverage(vInput);
double fSigSin= 0.;
double fSigCos= 0.;
for (size_t i= 0; i<nSamples; i++)
{
fSigSin+= sin(vInput[i]);
fSigCos+= cos(vInput[i]);
}
CircVal<UnsignedDegRange> Avrg2= atan2<UnsignedDegRange>(fSigSin, fSigCos);
double fErr1= CircVal<UnsignedDegRange>::Sdist(*sAvrg1.begin(), fAvrg); double fErr2= CircVal<UnsignedDegRange>::Sdist(Avrg2 , fAvrg);
fSumSqrErr1+= Sqr(fErr1);
fSumSqrErr2+= Sqr(fErr2);
}
double fRMS1= sqrt(fSumSqrErr1/ (nTrails-1)); double fRMS2= sqrt(fSumSqrErr2/ (nTrails-1));
f1 << nStdDev << "\t" << fRMS1 << "\t" << fRMS2 << endl; } );
Note that parallel_for is used for parallelizing tests for different standard deviations.
So which method is better?
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.
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 
of a continuous-time circular signal, sampled at times 
respectively, where the measurement technique/instrument is accurate - Estimate the mean value of the signal at period 
.
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 
(average).
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
For the above example, assume that the airplane heading measurements vector is {200, 240, 320, 40, 100, 20, 300}.
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
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 1
st 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:
template<typename T>
class CAvrgSampledCircSignal
{
size_t m_nSamples ;
CircVal<T> m_fPrevVal ; double m_fPrevTime; vector<pair<CircVal<T>, double>> m_Intervals;
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++;
}
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:
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);
Prize winner in Competition
"Best C++/MFC article of May 2011" 
General 
News 
Suggestion 
Question 
Bug 
Answer 
Joke 
Rant 
Admin