# Filtering Angular Measurements

I am trying to apply a filter to reduce noise in angular measurements. Since the data comes with noise, the filter needs to be able to handle measurements that are not smoothly progressing; for example, the filter should be able to work with a series like below:

2, 358, 3, 5, 10, 18 ... 304, 315, 330, 348, 9, 15


As you can see, the second measurement of 358 is an erroneous measurement. How can I reliably put this into a linear scale so that I can apply a regular filter like a moving average filter to the data.

I was initially thinking of setting a maximum change parameter, like 90 degrees, and if such a change is observed to correct that.

I am very new to signal processing. Maybe this is a very common problem that is easily handled. Hope you can put me in the right direction.

Median filtering is non-linear, and pretty awesome about removing outliers. You just need to adjust the length of the filter based on the estimate of the frequency of the errored samples.

• That's true, but the data doesn't actually has any outliers. The OP is just looking at it the wrong way. – Hilmar Dec 12 '20 at 16:36
• I didn't look close enough to realize that the measurements were in degrees. – IanJ Dec 12 '20 at 16:41
• Thanks, I can use mean of circular quantities for this right? en.wikipedia.org/wiki/Mean_of_circular_quantities – c00der Dec 12 '20 at 16:59
• If you did have errored samples you could. Its a little trickier to do median filtering with circular quantities, because the first thing it to sort them in order. One way of doing this could be to first check if all the values are within $90^\circ$ of each other. Or do they all lie within 2 of the 8 $45\circ$ wedges of the circle. If you go with the $45\circ$ method then you would find the 2 wedges with the most samples (others are outliers). Median filter or average the results. Make sure that you don't so a simple numeric average when the circle wraps around. – IanJ Dec 12 '20 at 17:33

As you can see, the second measurement of 358 is an erroneous measurement.

Why would that be erroneous? Phase is periodic with $$2\pi$$ or 360 degrees. That means 358 degrees is the same angle as -2 degrees. Your data looks perfectly fine if you look at it as

 2, -2, 3, 5, 10, 18 ,  ...


What you probably need is phase "unwrapping". If you see a phase difference of more than 180 degrees than simply add or subtract integer multiples of 360 until it's less than 180 degrees.

A good way to deal with circular (directional) averages is to turn it into a vector average.

To find the average angle $$\bar{\theta}$$ of several angles $$\theta_n$$ (in radians) then: $$\bar{\theta} = \arg \left ( \sum_{n=0}^{N-1} e^{i \theta_n} \right )$$ where $$\arg$$ is the argument (angle) of the resulting complex sum.

• Thanks. How I finally handled it was taking the sine and cos of each angle, then averaging them separately and converting that back to an angle (atan). Is that correct? en.wikipedia.org/wiki/Mean_of_circular_quantities – c00der Dec 13 '20 at 12:06
• @c00der Yes, that has precisely the same effect. $e^{i \theta} = \cos(\theta) + i \sin(\theta)$ so summing the $e^{i \theta}$s will have the same effect. Doing it my way saves the "divide by $N$" to average it... in fact, you probably don't need to average it, just sum them and take the arctan of the sums. The factor of $N$ will cancel. – Peter K. Dec 13 '20 at 21:22

If a first-order IIR will do, modify that slightly, and you're done.

So the usual first-order low-pass filter can be defined as $$y_n = h(\theta_n)$$ such that $$y_n = y_{n-1} + a(\theta_n - y_{n-1})$$. This works great for $$\theta_n \in \mathbb{R}$$.

You want a low-pass filter that's defined on an interval that spans $$360^\circ$$. For reasons that will become apparent in a moment, you want $$\theta_n \in \left [-180^\circ, 180^\circ \right)$$, where addition wraps around* sorta-kinda modulo $$360^\circ$$, i.e., $$170^\circ + 30^\circ = -160^\circ$$.

For want of a better notation, let $$\mathrm{wrap_{360}}(\phi)$$ be a sorta-kinda modulo function** that maps the real number line onto the interval $$\left[-180^\circ, 180^\circ \right)$$ modulo $$360^\circ$$, where $$\mathrm{wrap_{360}}(\phi) = \phi - 360^\circ \left \lfloor \frac{\phi + 180^\circ}{360^\circ} \right \rfloor$$

Then, the desired filter becomes $$y_n = h(\theta_n)$$ such that $$y_n = \mathrm{wrap_{360}}\left ( y_{n-1} + a\, \mathrm{wrap_{360}}(\theta_n - y_{n-1}) \right)$$.

The action of this filter is that for angles that stay close enough to $$0^\circ$$, the $$\mathrm{wrap_{360}}$$ function just acts as a passthrough. But for angles that get close to wrap-over (or for angles like yours, that start by being expressed in the range $$\theta \in [0^\circ, 360^\circ)$$) the actual angular difference is computed into something that expresses what it is physically -- that is, the difference between $$1^\circ$$ and $$359^\circ$$ is just $$2^\circ$$, not $$358^\circ$$.

Notes

• If your filter starts $$180^\circ$$ off, it will be at a false, but thankfully unstable, equilibrium. It may take an extra-long while to settle from startup.
• If your input angle is in some finite range around $$0^\circ$$, like $$350^\circ$$ to $$10^\circ$$, just recast it to $$-10^\circ$$ to $$10^\circ$$ and have fun.
• I can see a couple of ways that this notion could be extended to higher-order IIR filters, but it makes my brain hurt. If you need to do that, ask a separate question and someone will go to the effort to help you.
• If you want to use a FIR filter, you could first take your time sequence and "unwrap" it (Scilab and numpy both have "unwrap" functions that work on vectors, so I'm sure Matlab does to), then apply your FIR to the unwrapped version, then do whatever modulo operations you wanted to get it back to your desired interval.
• Alternately, you could do some sort of histogramming to determine the rough center of your average, shift everything, do the FIR, and shift back -- I think that would be computationally more expensive than unwrapping, though.
• Back when we actually had to worry about processor clock ticks, it was really convenient to map a 360 degree rotation into a 16- or 32-bit word, depending on the processor. On nearly all processors these days, in C, with a native $$n$$-bit word, if a and b are int, c = a + b; implements $$c = \mathrm{wrap_{2^n}}(a + b)$$. I always experienced an underlying thread of glee when I wrote code that took advantage of this, because usually integer overflow is the bane of your existence when you're doing fixed-bit signal procassing. Actually making it work for me made me feel like Tom Sawyer in the whitewashed fence story.

* I should be enough of a mathematician to tell you what such an entity is. I'm not. I don't think it counts as a field -- so we're stuck with lots of extra words.

** I just had a conversation with a mathematician who tells me that it's not a sorta-kinda modulo function: it's a plain old modulo function with a choice of equivalence classes that leads to the result being in $$[-\pi, \pi)$$ (I didn't want to offend him with degrees, so I talked in radians). Or more formally, that I'm just working with the additive quotient group, $$\mathbb R/2\pi\mathbb Z$$, with representatives chosen from $$[-\pi, \pi)$$, and I should just use the usual $$\mod 2\pi$$ notation. I'm not going to edit my text -- but there is terminology for this.