B4R Question moving average - how for sensor data ?

[peacemaker]

HI, All

Magnetic compass info is very variative: the heading (of the fixed sensor) is say, 30, 31, 35, 28, 27, 32, 37, 30, 35... around the 30 degrees value. How correctly to filter the stream of float values ?

 

[emexes]

How correctly to filter the stream of float values ?

Lol correctly is subjective, but for real-time applications I usually use:

Dim FilteredValue As Float

Sub HandleNewReading(NewValue As Float)
    Dim N As Int = 10    'smoothing factor, 0 = none, higher = smoother
    FilteredValue = (FilteredValue * (N - 1) + NewValue) / N    'moves FilteredValue 1/Nth towards NewValue
End Sub

If that looks too simple but you need to impress, then refer to it as Exponential Moving Average

 

[hatzisn]

Moving average is easy. You just take the average of last N (let's say 10) values obviously. Obviously also it is valid after the 10 measurements this brief code. I would do something like that:

        Dim fValue(10) As Float
        Dim ii As Int =0
        Dim fAvg As Float

       Sub LogValue(fMeasuredValue As Float)
        ii = ii+1
        If ii>9 Then ii=0
        fValue(ii) = fMeasuredValue
        For jj = ii To ii+10
            fAvg = fAvg + fValue(jj Mod 10) 
        Next
        fAvg = fAvg/10
       End Sub

 

[peacemaker]

Thanks, all. For now chosen such variant:

Private Sub Timer1_Tick
    Ready_flag = False   
    Timer1.Enabled = False
    
    Private avg(5) As Double
    For i = 0 To avg.Length - 1
        RunNative("read_hmc5883",Null)
        avg(i) = heading
    Next
    SortArray(avg)
    
    azimuth = avg(0)   'minimal one among 5
    Ready_flag = True
    Timer1.Enabled = True
End Sub

'Bubble sort
Sub SortArray(Values() As Double)
    Dim temp As Double
    Dim sorted As Boolean = False
    Do While sorted = False
        sorted = True
        For i = 0 To Values.Length - 2
            If Values(i) > Values(i + 1) Then
                temp = Values(i)
                Values(i) = Values(i + 1)
                Values(i + 1) = temp
                sorted = False
            End If
        Next
    Loop
End Sub

 

[hatzisn]

Lol correctly is subjective, but for real-time applications I usually use:

Dim FilteredValue As Float

Sub HandleNewReading(NewValue As Float)
    Dim N As Int = 10    'smoothing factor, 0 = none, higher = smoother
    FilteredValue = (FilteredValue * (N - 1) + NewValue) / N    'moves FilteredValue 1/Nth towards NewValue
End Sub

If that looks too simple but you need to impress, then refer to it as Exponential Moving Average

Thank you, I learned something today...

 

[emexes]

Private Sub Timer1_Tick
...
SortArray(avg)

@peacemaker I write this because we're all friends on this forum

SortArray and Timer_Tick are natural enemies. I appreciate that it's probably fast enough at the moment, but as your program grows and more channels are added, or the averaging size is increased, or users run it on devices that don't have a double-precision FPU...

I'm wondering if you meant to take the median value rather than the minimal - certainly I've seen that approach used for filtering, where the need to avoid outlier values outweighs the O(n**2) nature of a basic sort.

But if you're happy with the minimal value, then maybe consider eliminating the computational and memory load of having a buffer being sorted:

Private Sub Timer1_Tick
    Ready_flag = False
    Timer1.Enabled = False

    '''Private avg(5) As Double
    Dim MinHeading As Double    'same type as heading would be even better
    For i = 0 To 4    '''avg.Length - 1
        RunNative("read_hmc5883",Null)
        '''avg(i) = heading
        If i = 0 Or heading < MinHeading Then MinHeading = heading
    Next
    '''SortArray(avg)
 
    '''azimuth = avg(0)   'minimal one among 5
    azimuth = MinHeading
    Ready_flag = True
    Timer1.Enabled = True
End Sub

 

[emexes]

Magnetic compass info is very variative: the heading (of the fixed sensor) is say, 30, 31, 35, 28, 27, 32, 37, 30, 35... around the 30 degrees value.

azimuth = avg(0) 'minimal one among 5

Hmm... what about when the heading is say, 352, 2, 1, 353... ie +/- 5 degrees around the 357 degrees value?

That's going to clobber all of the filters above. Your approach of using minimal value does better than most, but won't return azimuths 350 to 359 (given sensor variative +/- 5 degrees).

A first draft of an EMA that handles this is:

Dim FilteredHeading As Float = 0    'needle north at program start?
Dim Smoothing As Int = 10    '1 = none, higher = smoother

Sub HandleReading(NewHeading As Float)

    Dim HeadingChange As Float  = NewHeading - FilteredHeading
    if HeadingChange < -180 then
        HeadingChange = HeadingChange + 360
    ElseIf HeadingChange > 180 Then
        HeadingChange = HeadingChange - 360
    End If
 
    FilteredHeading = FilteredHeading + HeadingChange / Smoothing
 
    If FilteredHeading < 0 Then
        FilteredHeading = FilteredHeading + 360
    ElseIf FilteredHeading >= 360 Then
        FilteredHeading = FilteredHeading - 360
    End If
 
End Sub

 

[peacemaker]

first draft

Yes, it needs to handle. But the C-code output is already 0....359, positive range.
And it seems to me that the trouble root is the lib (Adafruit_HMC5883_Unified) that has not the calibration, and the output is so much variative (look like +/- 20%). I will try first another library.

 

[emexes]

the C-code output is already 0....359, positive range.

but if we are then filtering multiple outputs, that cross the 359-0-1 degree wraparound point, then we need to account for the wraparound.

eg let's say that the variativeness is half what you're seeing, ie +/- 10% = +/- 36 degrees,

and we have 5 readings to filter being 342 +/- 36 degrees ie 306, 324, 342, 0, 18 degrees, then:

1/ a simple average of those would be (306 + 324 + 342 + 0 + 18) / 5 = 990 / 5 = "average" bearing of 198 degrees which is obviously nowhere near correct
2/ the median is 306 which is off by 36 degrees
3/ the minimum is 0 which is off by 18 degrees

whereas the filter in this post should work:

first draft of an EMA that handles this

lol, you're going to make me actually test it, aren't you? not that I'm complaining, because you've highlighted an unexpectedly interesting puzzle

 

[emexes]

And it seems to me that the trouble root is the lib (Adafruit_HMC5883_Unified) that has not the calibration, and the output is so much variative (look like +/- 20%). I will try first another library.

I agree that +/- 20% does seem a bit high, if this is +/- 72 degrees (wtf?!?!) and is occurring whilst the sensor is stationary.

So definitely try other libraries.

But regardless of how low the variativeness is, the issue of how to filter/smooth/average the readings across the 359-0-1 wraparound point will remain.

 

[peacemaker]

No, no any test from you side, i'm trying myself
Thanks for help.
Yes +/-20% is maybe too much i told , but anyway no good.

And yes - 359-0-1-trouble is always actual, for any lib, for any compass... It seems, values has to be checked <360 and >0 directly... but later.

 

[peacemaker]

BTW, another Grove_3Axis_Digital_Compass_HMC5883L lib is .... making the looped restart of the ESP32
And the example sketch on Github is another, not same as downloaded together with the lib.

 

[peacemaker]

Ha! Found.
The Adafruit_HMC5883_Unified lib works OK for getting the azimuth to the North along the X-axle, but ... only if the board are horizontal very strictly, i.e. when Y is around 0. If even Y ~ +/-3 degrees - the X-azimuth calculation is crazy wobbling. Arctan2(y, x) makes it very dependable from Y.
If horizontal - the result is rather smooth 0....359 degrees.

But the raw X-word from the registers is 0....180 for real 0...180 azimuth - clockwize rotation, but if anti-clockwize - it's 65536 65535 and smaller from 359 ....till 270, and again higher (to 65535) from 270 ... to 180.

 

[emexes]

But the raw X-word from the registers is 0....180 for real 0...180 azimuth - clockwize rotation, but if anti-clockwize - it's 65365 and smaller

Is it 65535 and smaller?

65535 is -1 signed-16-bit interpreted as unsigned-16-bit or as 32-bit.

65365 is likewise -171.

The simplest way to convert it to a signed value (eg -180 to -1) is to subtract 65536, assuming that the reading is stored in a signed 32-bit variable.

Or even simpler: read the raw X-word into a signed 16-bit variable (if B4R has that type).

 

[peacemaker]

Linear substraction is wrong, it's changed sinusoidal between each 90 degrees. So, no idea how to make the compass to work well with non-horizontal board yet...
But it's offtopic

p.s.

HMC5883L digital compass module gives 3-axis magnetic field component values, which can be used to determine the Azimuth angle (Heading). However, this method will give accurate Heading only when the module is held horizontal to the gound.

That is, when the module is tilted, the calculated Heading value will be incorrect. To compensate for the tilt, an accelerometer sensor is used to calculate the tilt angles (Roll and Pitch). The tilt angles are applied to a correction formula to obtain the accurate Heading value. Here, accelerometer values are read from MPU6050 Accelerometer+Gyroscope module.

GitHub - visakhanc/eCompass: Tilt Compensated Digital Compass based on AVR using HMC5883L Magnetometer, MPU6050 Accelerometer and OLED display

Tilt Compensated Digital Compass based on AVR using HMC5883L Magnetometer, MPU6050 Accelerometer and OLED display - visakhanc/eCompass

github.com

 

[emexes]

The simplest way to convert it to a signed value (eg -180 to -1) is to subtract 65536, assuming that the reading is stored in a signed 32-bit variable.

Linear subtraction is wrong,

Admittedly I am flying blind without seeing the declarations for (1) the library function that returns the X reading and (2) the variable that the reading is stored in, and (3) the code that transfers the reading from (1) to (2).

I am 97% certain that somewhere in the above is the cause of the large values where negative values are expected.

 

[emexes]

So, no idea how to make the compass to work well with non-horizontal board yet...

It is to do with the "compass" sensor being 3d ie there should be Y and Z readings, as well as X.

I would expect that if you have the board aligned orthagonally along one of the other two axis, then you will again get plausible compass readings but from the Y or Z registers instead of the X register.

The math for handling intermediate alignments is too hard for me, but I'm sure it's on the web somewhere.

 

[emexes]

idea how to make the compass to work well with non-horizontal board

A poor man's solution might be to confirm the board *is* horizontal, if the board also has acceleration sensors.

When the board is horizontal, one acceleration axis will be +/- 1 g(ravity) ie 9.8 m/s/s or 32 feet/s/s, and the other two axis will be close to 0.

 

[peacemaker]

A bit higher i already added the link to such project
So, the trouble is known.

p.s. finally solved so: https://www.b4x.com/android/forum/threads/compass-magnet-sensor-group-xxx5883.157604/