continuous averaging routine

[Yves]

Hello all,

Is anyone has a simple routine to generate a continuous averaging on 5 points. I'm having a noisy signal generated from an electro chemical reaction. The AD is generated a mV value every seconds and I will like to display on a OLED display a a smoother line using the previous 5 points and display the average. your help will be most appreciated.

Cheers,

Yves 

[trastikata]

Better using a power of two averaging, this way you can avoid floating point math. Something of the sort:

Code Select Expand

Dim wBuffer[8] as Word
Dim wTemp as word
dim wResult as word
Dim i as Byte

'Initialize buffer
wTemp = 0
For i = 0 to 7
wBuffer[i] = AdIn 0
wTemp = wTemp + wBuffer[i]
Next

While 1 = 1
'Shift buffer and adjust accumulator
wTemp = wTemp - wBuffer[0]
For i = 0 to 6
wBuffer[i] = wBuffer[i + 1]
next
wBuffer[7] = AdIn 0
wTemp = wTemp + wBuffer[7]
'Average the reaqdings
wResult = wTemp >> 3
Wend


[david]

I've used the classic LPF based on -
newavg = ((2^n-1)*oldavg + newsample)/2^n    (the higher the value of n, the longer the time constant is)
Example 1  (7*oldavg + newsample)/8
Example 2  (127*oldavg + newsample)/128

Ideally it needs floating point for the final value and if you're using it with a display then it should be rounded (up or down) to the nearest whole number.   
Here's a piece of working code used to provide a long averaging time constant for battery voltage sensing in a model boat speed controller.
 'Low Batt routine
        sample=ADIn 1                 '
        batt=255*oldbatt      
        batt=batt+sample         'add new value from ADC
        batt=batt/256            'divide total by 256 (bit shift)
        If batt<642 Then batflag=1  'change to suit battery used 
        oldbatt=batt            'update oldavg value with integer result

You can also test any of these functions with a simple spreadsheet by loading in some typical values and applying the formula to them and plot input and output curves.

Cheers,
David

[TimB]

For something like that I would make a ring buffer.

Every time you get an ADC value save it to an array and move the pointer. If the pointer goes past the end make the pointer 0

Then when needed just add up all the elements in the array in say a large enough variable.

Finally / by no elements and convert the raw to the the required unit

[Yves]

@TimB,

I like that it is a very elegant way of doing it. Thank you all for your valuable support.

Regards,

Yves

[top204]

Continuous averaging works well and I do sometimes use it as long as the samples never go too far apart in value. For unpredictable ADC samples I prefer a method using a sort mechanism of several ADC samples that are placed into the final array after they are sorted and the middle element's value used. I use a simple bubble sort procedure so the middle value can be taken or the lowest value or the highest value. It can be used as a low pass, band pass or high pass filter for values depending on which part of the sorted array is extracted for the returned value.

For example, the procedure below is a 16-bit median filter for an 18F or enhanced 14-bit core device that returns the middle (median) value from an array containing ADC samples. It is a bubble sort mechanism:

Code Select Expand

'------------------------------------------------------------------------------------------------
' An unsigned 16-bit bubble sort median filter procedure that uses indirect access to the array passed to it
' Input     : pArrayAddr holds the address of the array to sort
'           : pArraySize holds the size of the array to sort (3 to 255)
' Output    : Sorts the array passed to it
'           : Returns the median value from the array (middle value)
' Notes     : None
'
Proc Median16(ByRef pArrayAddr As Word, pArraySize As Byte), Word
    Dim tSwap_Occured As Bit                                    ' Indicates if the sort is complete
    Dim bElementCount As Byte                                   ' Holds a count of the array's elements
    Dim wElement_0    As Word                                   ' Holds the value of an element within the array
    Dim wElement_1    As Word                                   ' Holds the value of an element + 1 within the array
    Dim wArrayAddress As Word                                   ' Holds the start address of the array

    Repeat
        wArrayAddress = pArrayAddr                              ' Place the start address of the array into wArrayAddress
        tSwap_Occured = 0                                       ' Clear the flag that indicates a swap occured
        bElementCount = 1                                       ' Reset the array element counter
        Repeat                                                  ' Loop for each cell of the array...
            wElement_0 = Ptr16(wArrayAddress++)                 ' Extract an element from the array then increment the address
            wElement_1 = Ptr16(wArrayAddress--)                 ' Extract an element from the array then decrement the address
            If wElement_0 > wElement_1 Then                     ' Is the contents of element + 0 greater than the contents of element + 1?
                Ptr16(wArrayAddress++) = wElement_1             ' Yes. So store the element + 1 into element + 0 then increment the address
                Ptr16(wArrayAddress) = wElement_0               ' Store element + 0 into element + 1
                tSwap_Occured = 1                               ' Set the flag if a swap occurred
            Else                                                ' Otherwise....
                wArrayAddress = wArrayAddress + 2               ' Increment the array address to examine the next element pair
            EndIf
            Inc bElementCount                                   ' Check the next element of the array
        Until bElementCount = pArraySize                        ' Until the end of the array is reached
    Until tSwap_Occured = 0                                     ' Keep sorting until no more swaps occur
    pArraySize = pArraySize / 2                                 ' \
    pArraySize = pArraySize * 2                                 ' | Calculate the address of the middle element within the array
    wArrayAddress = pArrayAddr + pArraySize                     ' /
    Result = Ptr16(wArrayAddress)                               ' Return the value from the middle address of the array
EndProc


[top204]

The demonstration code listing below shows a way of using the above Median16 procedure to extract the median value from an array of samples:

Code Select Expand

'
'   /\\\\\\\\\
'  /\\\///////\\\
'  \/\\\     \/\\\                                                 /\\\          /\\\
'   \/\\\\\\\\\\\/        /\\\\\     /\\\\\\\\\\     /\\\\\\\\   /\\\\\\\\\\\  /\\\\\\\\\\\  /\\\\\\\\\
'    \/\\\//////\\\      /\\\///\\\  \/\\\//////    /\\\/////\\\ \////\\\////  \////\\\////  \////////\\\
'     \/\\\    \//\\\    /\\\  \//\\\ \/\\\\\\\\\\  /\\\\\\\\\\\     \/\\\         \/\\\        /\\\\\\\\\\
'      \/\\\     \//\\\  \//\\\  /\\\  \////////\\\ \//\\///////      \/\\\ /\\     \/\\\ /\\   /\\\/////\\\
'       \/\\\      \//\\\  \///\\\\\/    /\\\\\\\\\\  \//\\\\\\\\\\    \//\\\\\      \//\\\\\   \//\\\\\\\\/\\
'        \///        \///     \/////     \//////////    \//////////      \/////        \/////     \////////\//
'                                  Let's find out together what makes a PIC Tick!
'
' Demonstrate a bubble sort median filter for an unsigned 16-bit array
'
' Written for the Positron8 compiler by Les Johnson.
'
    Device = 18F25K20                                           ' Tell the compiler what device to compile for
    Declare Xtal = 16                                           ' Tell the compiler what frequency the device will be operating at (in MHz)
'
' Setup USART1
'
    Declare Hserial1_Baud = 9600                                ' Set the Baud rate to 9600
    Declare HRSOut1_Pin = PORTC.6                               ' Set the TX pin
'
' Create any variables for the demo here
'
    Dim wArray[7] As Word                                       ' Holds the 16-bit array to sort
    Dim wMedianValue As Word                                    ' Holds the median value from the array
    Dim bIndex As Byte                                          ' Holds the index to the array for displaying purposes

'-------------------------------------------------------------------------
' The main program starts here
' Sort a 16-bit array and extract its median (middle) value and display it on a serial terminal
'
Main:
    wArray = 800, 900, 700, 0, 120, 230, 10                     ' Load the array with some values to sort

    HRsoutLn "Unsorted:"                                        ' \
    For bIndex = 0 To Bound(wArray)                             ' | Display the unsorted array's contents on a serial terminal
        HRsout1 Dec wArray[bIndex], ","                         ' |
    Next                                                        ' /

    wMedianValue = Median16(wArray, SizeOf(wArray))             ' Sort the array passed to it and return its median value

    HRsoutLn "\rSorted:"                                        ' \
    For bIndex = 0 To Bound(wArray)                             ' | Display the sorted array's contents on a serial terminal
        HRsout1 Dec wArray[bIndex], ","                         ' |
    Next                                                        ' /

    HRsoutLn "\rMedian value is ", Dec wMedianValue             ' Display the median (middle) value extracted from the array


'------------------------------------------------------------------------------------------------
' An unsigned 16-bit bubble sort median filter procedure that uses indirect access to the array passed to it
' Input     : pArrayAddr holds the address of the array to sort
'           : pArraySize holds the size of the array to sort (3 to 255)
' Output    : Sorts the array passed to it
'           : Returns the median value from the array (middle value)
' Notes     : None
'
Proc Median16(ByRef pArrayAddr As Word, pArraySize As Byte), Word
    Dim tSwap_Occured As Bit                                    ' Indicates if the sort is complete
    Dim bElementCount As Byte                                   ' Holds a count of the array's elements
    Dim wElement_0    As Word                                   ' Holds the value of an element within the array
    Dim wElement_1    As Word                                   ' Holds the value of an element + 1 within the array
    Dim wArrayAddress As Word                                   ' Holds the start address of the array

    Repeat                                                      ' Create a loop for the sorting
        wArrayAddress = pArrayAddr                              ' Place the start address of the array into wArrayAddress
        tSwap_Occured = 0                                       ' Clear the flag that indicates a swap occured
        bElementCount = 1                                       ' Reset the array element counter
        Repeat                                                  ' Loop for each cell of the array...
            wElement_0 = Ptr16(wArrayAddress++)                 ' Extract an element from the array then increment the address
            wElement_1 = Ptr16(wArrayAddress--)                 ' Extract an element from the array then decrement the address
            If wElement_0 > wElement_1 Then                     ' Is the contents of element + 0 greater than the contents of element + 1?
                Ptr16(wArrayAddress++) = wElement_1             ' Yes. So store the element + 1 into element + 0 then increment the address
                Ptr16(wArrayAddress) = wElement_0               ' Store element + 0 into element + 1
                tSwap_Occured = 1                               ' Set the flag if a swap occurred
            Else                                                ' Otherwise....
                wArrayAddress = wArrayAddress + 2               ' Increment the array address to examine the next element pair
            EndIf
            Inc bElementCount                                   ' Check the next element of the array
        Until bElementCount = pArraySize                        ' Until the end of the array is reached
    Until tSwap_Occured = 0                                     ' Keep sorting until no more swaps occur
    pArraySize = pArraySize / 2                                 ' \
    pArraySize = pArraySize * 2                                 ' | Calculate the address of the middle element within the array
    wArrayAddress = pArrayAddr + pArraySize                     ' /
    Result = Ptr16(wArrayAddress)                               ' Return the value from the middle address of the array
EndProc

On the serial terminal, it will display:

Unsorted:
800,900,700,0,120,230,10,
Sorted:
0,10,120,230,700,800,900,
Median value is 230

It shows the unsorted and sorted values contained in the array, and it can be seen that the median (middle) value of the array's contents is 230.

If required, place the median values of a few samples into a continuous averager and it will produce very stable values because any extreme values are sorted out and the middle values are in the middle of the sorted array.

Or just use the middle value from the Median filter because the more of the same sample values there are, the closer they will be to the centre of the array, and after each ADC sample, load it into the same sorted array and run the median sample procedure again to sort it again. This way, an average will be gotten because the extremes will be removed.

[david]

What is the advantage of using a buffer or array over the simple averaging as shown in the attached?    The routine used uses no buffer and only requires one value stored.  It uses only a multiply, an add and a divide (bit shift) so can be made very fast.  By changing the value of n the degree of averaging can be simply changed from close to sampled data to a very flat, slow moving average as indicated by the filtering of a square wave by various degrees.
I know buffers are often used in averaging but I can't understand what it offers over the simple approach so hopefully someone can explain this to me.

Cheers,
David

[TimB]

David

Would you upload your sheet?

I have tried a few  averaging systems but not yours and would like to see how it performs. I'm thinking perhaps I try add a ring to the list and look at what the result is on the graph using your data

[trastikata]

I know buffers are often used in averaging but I can't understand what it offers over the simple approach so hopefully someone can explain this to me.

Hello David,

in your Excel table you are working with double floating point values and in the PIC it is supposed to be integer math.

Now try the same spreadsheet but to the result of each division add INT() function to extract the resultant integer. Note ROUND() is not the same as it will round to the nearest integer.

Besides what you are using is a IIR (infinite impulse response) whereas the buffered filtering is a FIR (finite impulse response) filter.

There are advantages and disadvantages.



Here's a demonstration of both filters with averaging of 4. Which one will be used is up to designer depending on the wave form. Say in this example the IIR filter takes 4 values of 250 and returns something else ... which might not be acceptable. With more complex wave forms the situation is no longer so obvious as to what filter to use.

[david]

Hi Tim,
I've attached a file with a fudged extension - just change the .txt to .xlsx and hopefully it should work. If anyone wants it in .ODS format let me know.

Cheers,
David

[david]

Trastikata - I will make the changes and get back to you shortly.
I'm not sure I follow the subtle difference in needing to use Int rather than Round but I will investigate both.
I was aware it was an IIR filter but was unaware of any limitations it may have.  Many thanks for your post.

Cheers,
David

[david]

Hi All,
Well it seems that using an integral function does cause problems on the more aggressive filtering but then this is equivalent to continuously rounding down and hence applies a bias to the values.
The multiplication and addition process are both integer functions so if you used a float for the division and rounded it to the nearest integer then the bias is removed and results look much better.   The 255:256 filter still has a difference but it doesn't look terminal.  I will extend the sample source out a bit more to determine where it is going.
Where are all the smart college kids when you need them?  Trastikata is probably asleep now but there must be some in the UK.
I've attached screen shots of the two different methods with an example formula near the top.

Cheers,
David

[david]

Just to round off this saga here is a plot of more typical data (constrained random values) rather than looking at the LPF effects on a square wave.
I think it's unlikely that you would need to go more than 15:16 and generally no more than 3:4 given that these samples are 33% noise.

Cheers,
David

[top204]

I often use different moving average filters, and it is a preference of mine to use sorting of samples for an average, not because it is better. But in my findings, I have found the sorting method a bit faster for slow moving samples because a lump of them can be placed in the single array and an average taken from them.

I have taken your pseudo code snippet from above and created a procedure around it, with a demo showing it in operation. It is listed below:

Code Select Expand

'
'   /\\\\\\\\\
'  /\\\///////\\\
'  \/\\\     \/\\\                                                 /\\\          /\\\
'   \/\\\\\\\\\\\/        /\\\\\     /\\\\\\\\\\     /\\\\\\\\   /\\\\\\\\\\\  /\\\\\\\\\\\  /\\\\\\\\\
'    \/\\\//////\\\      /\\\///\\\  \/\\\//////    /\\\/////\\\ \////\\\////  \////\\\////  \////////\\\
'     \/\\\    \//\\\    /\\\  \//\\\ \/\\\\\\\\\\  /\\\\\\\\\\\     \/\\\         \/\\\        /\\\\\\\\\\
'      \/\\\     \//\\\  \//\\\  /\\\  \////////\\\ \//\\///////      \/\\\ /\\     \/\\\ /\\   /\\\/////\\\
'       \/\\\      \//\\\  \///\\\\\/    /\\\\\\\\\\  \//\\\\\\\\\\    \//\\\\\      \//\\\\\   \//\\\\\\\\/\\
'        \///        \///     \/////     \//////////    \//////////      \/////        \/////     \////////\//
'                                  Let's find out together what makes a PIC Tick!
'
' A continuous moving average filter demo, based on David's pseudo code from the compiler's forum
'
' Written for the Positron8 compiler by Les Johnson.
'
    Device = 18F25K20                                               ' Tell the compiler what device to compile for
    Declare Xtal = 16                                               ' Tell the compiler what frequency the device will be operating at (in MHz)
'
' Setup USART1
'
    Declare Hserial1_Baud = 9600                                    ' Set the Baud rate to 9600
    Declare HRSOut1_Pin   = PORTC.6                                 ' Set the TX pin

    Include "Amicus18_ADC.inc"                                      ' Load the 18F25K20 ADC routines into the program
'
' Create the variables for the demonstration
'
    Dim wSampleCount  As Word = 1                                   ' Holds the amount of samples taken
    Dim wRaw_Sample   As Word                                       ' Holds the raw ADC sample value
    Dim wAverageValue As Word                                       ' Holds the averaged value
    Dim wPrevAverage  As Word = 0                                   ' Holds the previous averaged value for the display

'-------------------------------------------------------------------------
' The main program starts here
' Take 10-bit ADC readings and average them and display them on a serial terminal
'
Main:
'
' Open the ADC for a 10-bit reading and make pin AN0 analogue
'
    OpenADC(ADC_FOSC_RC & ADC_RIGHT_JUST & ADC_2_TAD, ADC_REF_VDD_VSS, ADC_1ANA)

    Do                                                              ' Create a loop
        wRaw_Sample = Adin 0                                        ' Take an ADC reading for pin AN0
        wAverageValue = Get_Average16(wRaw_Sample)                  ' Get an average of it
        '
        ' Display the values on a serial terminal, but only if they are different to the last values
        '
        If wAverageValue <> wPrevAverage Then                       ' Is the current average value the same as the previous value?
            HRsoutLn "Smpl ", Dec wSampleCount, " : ADC value is ", Dec wRaw_Sample, " : Average Is ", Dec wAverageValue ' No. So display the raw and averages ADC values
            wPrevAverage = wAverageValue                            ' Load the previous value with the current value
        EndIf
        Inc wSampleCount                                            ' Increment the amount of samples taken
    Loop                                                            ' Do it forever

'-------------------------------------------------------------------------
' A continuous moving average filter for unsigned 16-bit samples
' Input     : pSample holds the 16-bit sample value
' Output    : Returns the averaged 16-bit value
' Notes     : Uses floating point for the internal calculations
'
Proc Get_Average16(pSample As Word), Word
Static Dim fPrevValue As Float = 0                                  ' Holds the previous averaged value
    Dim fTemp As Float                                              ' A temporary variable for calculations
    Symbol cDivMulsReq = 32                                         ' The larger the number, the more samples required for the final average value
    Symbol cDivMulsReqPlus1 = (cDivMulsReq + 1)

    fTemp = (cDivMulsReq * fPrevValue) + pSample                    ' Add the sample value
    fTemp = fTemp / cDivMulsReqPlus1                                ' Divide the total by the amount of multiplications plus one
    fPrevValue = fTemp + 0.001                                      ' Update fPrevValue value with rounded up result
    Result = fPrevValue                                             ' Return the integer result
EndProc

The averaging calculations are using floating point within the procedure, but I think they could be changed to integers with some multiplications and divisions of 10 or 100 for the values so they do not underflow. With the original multiplication of 255 value in the filter, it takes a lot of samples before the average is calculated, as the screenshots of the code simulated show, and the smaller the multiplications/divisions the less samples it requires, but 'I would imagine' the less averaging is performed as well.

The "Samples_255.jpg" image shows that it requires 1774 samples before the averaging calculations give the averaged value back with the multiplication of 255, and the "Samples_32.jpg" image show it takes 293 samples before it gives the correct average back with a multiplication of 32. Is this correct David? The maths seems to suggest it is, and the code demo.

If the procedure's Static fPrevValue variable is loaded with an initial value closer to the final average value, or greater than 0, the averaging takes a lot less samples. Is there a preferred initial value for fPrevValue david?

[top204]

To counteract the initial previous average value being low and taking too many samples for the final average value, I have altered the procedure so that the fPrevValue variable is pre-loaded with the initial pSample value, but only once. This speeds up the averaging a heck of a lot:

Code Select Expand

'-------------------------------------------------------------------------
' A continuous moving average filter for unsigned 16-bit samples
' Input     : pSample holds the 16-bit sample value
' Output    : Returns the averaged 16-bit value
' Notes     : Uses floating point for the internal calculations
'
Proc Get_Average16(pSample As Word), Word
Static Dim tPrevLoaded As Bit = 0                                   ' Set to 1 when the fPrevValue variable is pre-loaded with the initial pSample value
    Dim fPrevValue As Float                                         ' Holds the previous averaged value
    Dim fTemp As Float                                              ' A temporary variable for calculations
    Symbol cDivMulsReq = 16                                         ' The larger the number, the more samples required for the final average value
    Symbol cDivMulsReqPlus1 = (cDivMulsReq + 1)

    If tPrevLoaded = 0 Then                                         ' Is tPrevLoaded 0?
        tPrevLoaded = 1                                             ' Yes. So set it to 1
        fPrevValue = pSample                                        ' Load fPrevValue with the contents of pSample, but only once
    EndIf
    fTemp = (cDivMulsReq * fPrevValue) + pSample                    ' Add the sample value
    fTemp = fTemp / cDivMulsReqPlus1                                ' Divide the total by the amount of multiplications plus one
    fPrevValue = fTemp + 0.001                                      ' Update fPrevValue value with rounded up result
    Result = fPrevValue                                             ' Return the integer result
EndProc


[david]

Interesting stuff.
Using a 255:256 filter will take a very long time and represents a massive amount of filtering.  As I mentioned in last night's post, I would doubt you would ever need to go past 7:8 or 15:16 at the most.  As you found, allocating the first value to be somewhere near the expected average is critical for the quicker startup.  I guess the procedure could be adapted to have a dynamic filter by changing the one value.
What we need now is someone who has some logged data we can experiment on....

Cheers,
David

[TimB]

I found this thread really really interesting. I'm forever trying to average readings.

My conclusion is that none of the routines I have tried do much more than using a ring buffer and averaging the buffer when you need the data out.

The routine is really fast as to add to the buffer is just put in an array, When you need to get the average of the past x values its all addition and one division. If you stick to array sizes that are 2's complement eg 2,4,8,16,32 the division (if not a float) evne this division is just a few shifts

I preload the buffer in my initialisation routine and I'm out the gate flying

[david]

Hello Tim,
We need a maths expert with a background in digital signal processing. 
The buffer approach is widely used but would get a bit weighty if trying to do 256 samples whereas the IIR approach is very scalable. 
The standard, basic IIR formula-
newavg=(2^n-1*oldavg + newsample)/2^n or as an example newavg=(15*oldavg+newsample)/16  (multiply, add, bit shift divide)
could be moified (I think) to this-
newavg=(2^n*oldavg -oldavg + newsample)/2^n or as an example newavg=(16*oldavg-oldavg+newsample)/16  (bitshift, subtract, add, bit shift)
Perhaps someone can comment if this would be a somewhat faster routine.

In both methods the division can introduce fractional values which should I think be rounded up or down to the nearest integer but at least in the buffered approach they will exit the function whereas in the IIR approach fractional values could tend to accumulate and rounding would be essential (as opposed to Int down) to avoid the possibility of bias.
As I mentioned before, I think it's unlikely that you would need to use more than 16 values to average data samples although I have used 255 where I simply wanted a very long time constant.
I'm sure we will all carry on using what we're familiar with but because of my limited understanding of the IIR filters I am still a little apprehensive about it even though floating point maths indicates it works well.

Cheers,
David

[david]

Some of you whose eyes haven't yet glazed over on this subject may have noticed that the simplified formula above can be further simplified-

newavg=(2^n*oldavg -oldavg + newsample)/2^n    is the same as-

newavg=oldavg +(newsample-oldavg)/2^n

This reduces the IIR filter down to an add, a subtract and a bit shift divide.

As for the pros and cons of IIR vs FIR filters-   (not my words)
▪ A single-pole IIR filter is easy to implement and requires a storage space of just one memory element to store the past output. However, the settling time increases as the attenuation factor increases and the filter takes infinite samples to settle to the input value.
▪ A moving average FIR filter is relatively complex to implement and requires a storage space of N elements, where N is the number of samples averaged.  The filter takes just N samples to settle to the input value.
▪ A moving average filter with N elements has a frequency response roll-off comparable to an IIR filter with attenuation factor, a = N/2. Or, in other words, the roll-off of a 32-tap moving average filter can be matched by an IIR filter with an attenuation factor of just 16.
However, if you compare time responses, it takes 108 samples for the IIR filter to settle to 99.9% of the input while the moving average FIR takes just 32 samples to settle to 100% of the input signal.

Cheers,
David