How to average large numbers on a PIC?

[Mr RB]

Hi, does anyone have an easy way to keep a long term running average of large numbers, signed long 32bit numbers which are limited by definition to a max size of about 2.14 billion.

The goal is to average a high number of samples that are very similar over a long time.

The samples are based on the PIC xtal freq so typical samples for a 4MHz xtal would be;
4000010
4000013
4000009 (etc)

and may even be more similar samples;
4000010
4000011
4000010 (etc)

and may be as high as a 20MHz xtal;
20000010
20000011 (etc)

As the PIC 16F has limited amounts of RAM this rules out having hundreds of 32bit variables, and an accumulator system can only be 100* samples or so without risking overflowing the signed long variable.

I have had some success with a 100* accumulator;
accum -= (accum / 100);
accum += new_sample;

But this is very poor at handling the average of such fine differences.

An ideal result would be after a long time to have the accumulator with an accurate total based on the average;
4000010.37 (where the decimal point is imaginary)

Is there a simple system that will fit within these constraints? Thanks!

 

[misterT]

If your compiler doesn't support 64bit integers, you could try multiple accumulators.
Or you could try subtracting some big value from each number.. to get rid of the "DC" value

EDIT:
You could use the very first sample as an estimate of the average. Call that number M and save it in a variable.
Subtract that number M from every sample that follows before adding it to the accumulator.

SUM = SUM + NewValue - M;
N = N+1;

'N' is the number of samples. The average is then:

Average = (SUM/N) + M;

 

[3v0]

I recall reading that keeping the average has the same result as keeping the samples it was obtained from.

What I use is similar to what misterT suggested but retains the sum of the last 8 readings. If you do a get prior to the 8th sample it will be bogus.

I have attempted to make this generic enough so that non c programmers can understand it

// Each time you obtain a new sample do a put
void put(int newReading)
{
  static unsigned char samples = 0;  // this line executes once during program initialization
  if(samples < 8)  // add newReading's till adcAbs8 is sum of 8 readings
  {
    samples = samples + 1;
    adcAbs8 = adcAbs8 + newReading;
  }
  else
  {
    ave = adcAbs8 / 8;  
    adcAbs8 = adcAbs8  -ave + newReading;
  }
}

// to use the average value of the parameter do a get
int get(void)
{
   val = adcAbs8 / 8;  
   return val;
}

The code is a simplification of existing working code. Could have botched the simplification. Always a risk to post untested code, even simple stuff.

 

[dougy83]

A simple method of getting a running average is using a LPF, e.g. the simplest would be:

A = ß.A + (1-ß).x

, where A is the average, ß is the weighting of the old average vs the new value and x is the new sample value. 'A' should be initialised to the expected average (e.g. the first sample reading will do) as the convergence can be a little slow. A useful value for ß might be ~0.99 - 0.999; but use power of 2^(-1) to allow simple use of fractional math.

misterT's suggestion re removing the 'DC' to obtain the error value can also be used here to again reduce the size of the accumulator.

 

[crutschow]

Another variation for doing a running average is to add a percentage of the difference between the average and the new sample to the previous average. This somewhat emulates a LP RC filter.

A = A + (x-A)*β

where A is the running average, x is the new sample, and β is the factional part of the difference to save.

The smaller the β, the slower the rate of change of the average. β is rather like the inverse time-constant of the running average.

 

[Diver300]

A = A + (x-A)*β

A = ß.A + (1-ß).x

The only difference between those is that β in one is (1-ß) in the other.

I have used that technique quite often to filter the result from A2D converters.

 

[crutschow]

The only difference between those is that β in one is (1-ß) in the other.

True.

My form of the equation has only one multiplication which may make it easier to implement.

 

[Mr RB]

3v0- Thanks for the suggestion to average X number of samples, I have used that system before but in this case had already ruled it out as I need to average hundreds or thousands of 32bit samples over a long time.

Thanks MisterT, that's correct my compiler does not have variables >32bit for PIC16F. Double and long double are still handled as 32bit vars.

Dougy83 and Crutschow, thanks for the low pass filter suggestions, I used a similar one in the orig post; A = 0.99*A + X but the trouble of course is that it does not have the resolution needed for fine averaging.

Thanks also to everyone for the suggestion to separate the error value from the large static value (or "DC" ) I had a think last night and came to a similar conclusion. The real problem is in the implementation, I like the low pass filter accumulator but when averaging thousands of samples that differ only in 1 integer count it must have some system for handling extra resolution.

So based on everyone's input it looks like the best system is going to be something along these lines of using Crutchow's filter;
A = A + ((x-A)*β)
1. accumulator holds the average as 32bit integer (same scale as the samples)
2. another 32bit (16bit?) value holds all the sub-integer resolution of the average (makes a 64bit or 48bit accumulator?)
3. an "error" value is calculated for every new sample; error = new - average
4. error value is scaled down, ie; error = (error / 100)
5. error value is added to the 64bit average in the accumulator

I think that looks great, it should work well. Now the problem is in handling that sub-integer variable... It's going to get ugly.

Here are some inital thoughts (please comment if you have improvements);
accum_int = signed 32bit
accum_sub = unsigned 16bit
new_sample = signed 32bit

Since new_sample must be an integer we can get an integer error to start with; (signed)error_int = new_sample - accum_int

Then that error value can be adjusted by subtracting the (unsigned) value from the accum sub-integer variable;
error = (error_int << 16); // shift error to make it the same scale as the sub-integer data
error -= accum_sub; // include the sub-integer data

at this point error holds the full (48bit resolution) error value? Then it can be scaled; error = (error / 100)
Now how to break this signed error value back up into int-and sub-int data and modify the 2 accum variables? And then I have to display it in decimal as; int.sub
There has to be a neater way...

 

[dougy83]

Thanks also to everyone for the suggestion to separate the error value from the large static value (or "DC" ) I had a think last night and came to a similar conclusion. The real problem is in the implementation, I like the low pass filter accumulator but when averaging thousands of samples that differ only in 1 integer count it must have some system for handling extra resolution.

That's why it was suggested to use fractional math.

Re the implementation you have outlined, perhaps something _like_ the following would work:

union Accumulator
{
    struct
    {
        uint16_t fractional;
        uint32_t whole;
    };

    struct
    {
        uint32_t low32;
        uint16_t high16;
    };
} acc;

uint32_t updateAcc(uint32_t sample)
{
    int32_t error = acc.whole - sample;

    // set beta scale factor here by shifting left...
    error <<= 0;    // beta = 2^(shift value - 16)

    uint32_t old32 = acc.low32;
    acc.low32 += error;

    // this needs optimising
    if(error < 0 && old32 < acc.low32)
        acc.high16--;       // borrow
    else if (error > 0 && old32 > acc.low32)
        acc.high16++;       // carry

    return acc.whole;
}

It might be better to implement the function in assembler so you can have a 'native' 48-bit variable to work on

 

[3v0]

My method implements A = A + ((x-A)*β) . My values are scaled by 8 which retains 3 bits more information and β=1/8.

 

[Mr RB]

My apologies 3v0, your code does indeed do a running low pass filter! I only glanced through your code and didn't notice it was a running filter, funny because it was an almost identical accumulator filter to the one I posted in #1.

Dougy83 thanks for the code, I would never have thought of using Unions as a way of accessing a 48bit variable in 2 different ways. I think that definitely solves some of the problems with accessing the sub-integer data and it looks close to a pretty good working solution.

However before we put it to rest, I had an idea which on reflection should be credited to MisterT as on re-reading his post #2 seems to be pretty much what he was saying;

If an arbitrary large integer is chosen and subtracted from the new_sample, the number to be averaged becomes a very small figure and can easily be averaged in a single 32bit variable. More importantly, this "error" value is now an integer for every sample. Now it doesn't matter too much what that arbitrary number is.

As an example consider these incoming new samples;
4000010
4000011
4000008
4000004

if the arbitrary number is 4000005, then the "error" values look like this;
4000010 +5
4000011 +6
4000008 +3
4000004 -1

So all we need to do seems to be to average the small "error" integers; 5+6+3-1 = 13 /4 so for these samples the error_average is 13/4 or 3.25

Now the running average is found by adding arbitrary number with error_average;
4000005 + 3.25 = 4000008.25

This has the big advantage of speed and simplicity as it only has to average small integers now, and the closest thing to a difficult task is the final display of the large average number. That should be easy enough to separate the integer from the sub-integer and just do separate display systems on either side of the decimal point?

Maybe I'm missing something, but this does look a little "too good to be true"...?

 

[MrAl]

Hi there MrRB,

It looks sound. Here's why...

If we take four original numbers a,b,c,d, and add them together and divide by 4 we get the average:
avg=(a+b+c+d)/4

but if we first subtract some number n from each sample first we have:
y=((a-n)+(b-n)+(c-n)+(d-n))/4;

The question then becomes what can we do to y to get the average?

Simplifying y we get:
y=(-4*n+d+c+b+a)/4

and expanding that a little we get:
y=-n+d/4+c/4+b/4+a/4

and now expanding the average we get:
avg=a/4+b/4+c/4+d/4

Now it becomes clear to see that the only difference between avg and y is the negative of the number, so if we add n back to y we'll again get the average:
y+n = (-n+d/4+c/4+b/4+a/4)+n = d/4+c/4+b/4+a/4 = avg

and that's the original average so this proves that's the same way of doing it at least numerically.

 

[Mr RB]

Thank you MrAL for the math proof! As I don't have your math skill I tried a "nuts and bolts" approach this morning, I coded it up, injected known data into the front end and displayed the result coming out the back end.

My results back up yours (of course), it doesn't matter what the arbitrary number is as all the variance of the data is always correctly averaged in the error values.

Next I took the logical step of "trimming" the error value to force the "arbitrary" average integer to be the FULL integer value, and ensure the sub-integer data (the error average) was always stored as decimal in the sub-integer range 0-0.999999. That value was chosen for display simplicity, so the sub-integer data is basically scaled so 1000000 = 1 integer (ie; 0-999999 range). This gives a lot of nice decimal places for the averaging and greatly simplifies display.

I tested it farly comprehensively and here is the relevant source code used for the averaging and display;

#define FILT1 40        // faster filter constant used for first 100 seconds
#define FILT2 600       // slow filter constant used after that
unsigned int filter;    // holds filter constant
signed long ierror;            // for calcs
signed long real_xtal_freq;    // holds new sample
signed long real_xtal_av;      // average, as integer
signed long real_xtal_avsub;   // average extra resolution data (sub-integer)
unsigned char cycles;          // how many samples were tested
unsigned char txt[14];         // for text string display

    // (loop for each sample) the new sample is in; real_xtal_freq
    // calculate a running average in a Low Pass Filter accumulator.
    // only do averaging after it has settled down (more than 4 captures)
    if(cycles < 255)  cycles++;        // will stop counting at 255
    if(cycles == 4)  real_xtal_av = real_xtal_freq;   // starting value for average
    if(cycles > 4)   
    {
      // Now do the averaging. The process is to average the small integer error 
      // between the captured value and the integer portion of the average.
      // these small integer errors are then averaged in real_xtal_avsub,
      // and this becomes extra resoluton data (sub-integer data) of the average.
      // the sub-integer data is kept as decimal, scaled so; 1000000 = 1 integer
      if(cycles == 104)
      {
        filter = FILT2;             // use finer filtering after first 100 samples
        RomanLCD_Out(0,15,"f");     // show "fine" tag on LCD
      }
      ierror = (real_xtal_freq - real_xtal_av); // get small integer error
      real_xtal_avsub -= ((real_xtal_avsub+(filter/2))/filter);   // remove a rounded % of average
      ierror = (ierror * (1000000/filter));     // scale new error to % of filter value
      real_xtal_avsub += ierror;                // add into sub-integer average

      // if at any point the sub-integer average becomes >1 or <0 we remove
      // integer counts from it and put them in the integer average.
      while(real_xtal_avsub >= 1000000)    // if sub is > 1 integer
      {
        real_xtal_avsub -= 1000000;
        real_xtal_av += 1;
      }
      while(real_xtal_avsub < 0)           // if sub is < 0 integer
      {
        real_xtal_avsub += 1000000;
        real_xtal_av -= 1;
      }

      // at this point we must have an integer average in real_xtal_av
      // and the 0-0.999999 sub-integer resolution in real_xtal_avsub 
            
      // display the 0-0.999999 data that will appear after the decimal point
      // this is currently in a 6 digit decimal format; 0-999999
      ierror = (real_xtal_avsub / 1000);  // make 6 dec digits into 3 decimal digits
      LongToStr(ierror,txt);              // format into text string
      txt[7] = '.';                       // manually insert dec point
      if(txt[8] == ' ') txt[8] = '0';     // add leading zeros where needed
      if(txt[9] == ' ') txt[9] = '0'; 
      RomanLCD_Out(1,4,txt);              // display it
    
      // and lastly display the integer average
      LongToStr(real_xtal_av,txt);      // format into text string
      RomanLCD_Out(1,0,txt);            // display it
      RomanLCD_Out(1,0,"av");
      // the result displays like this; 4000000.000
    }

Here is the LCD readout from my 4MHz PIC xtal, this project was for the PIC to sample its own xtal against a GPS 1pps pulse, with zero-accumulated-error systems. The averaging is also done with zero error systems (if I got the math right) and the only error in the entire system is that the last decimal digit .--x is rounded down which I can't be bothered fixing as it fluctuates a couple of counts anyway.

**broken link removed**

This is my PIC xtal running at exactly 4.000148341 MHz (give or take a couple of fluctuating counts).

I was able to determine my 4MHz xtal drops speed by about 1.5 PPM (Parts Per Million) when heated 10'C above ambient. The measurement resolution and accuracy are so good I could easily see the xtal speed increasing by 0.03 PPM or so when I stood up and paced around a little, I believe this was from drafts created by my body movement cooling the xtal by about 0.2'C.

Likewise when leaning over the xtal to take a photo with my body heat maybe 18 to 20 inches from the xtal I could see the freq rising by 0.05 PPM or more. I'm sure this was temperature related not capacitive as putting a metal object near the xtal had much less effect.

The GPS module was a LEA5S which according to it's datasheet has a 1pps pulse accurate to 50nS timing. However that timing spec does not matter that much as the zero-error measurement systems retain all edge error and timer error and average it into the final result.

I will write up the whole project soon and put it on my web page with credit given to everyone who helped in this thread.

 

[MrAl]

Hi again Mr RB,

That's an interesting idea really. I wouldnt mind trying it myself. Are you saying the module has a pulse output then?
Where did you get the module and how much did it cost approximately?

 

[Mr RB]

The module I am using is this one;
**broken link removed**
at $59.95 and has 1pps output and standard logic-level serial data output;
MikroElektronika - SmartGPS Board - GPS Development Tool (LEA6S)

Actually mine is an older model with LEA-5S, the new one seems to have the LEA-6S module.

They also have a cheaper version at $49.95;
MikroElektronika - EasyGPS Board - GPS Prototype Development Tool (LEA6S)

You can buy the LEA6S and make your own PCB but these modules have the RF PCB design and SMD soldering taken care of and have inbuilt SMPS supply and 5v/3v regulators etc etc.

Since I have a web page full of algorithms to generate very high resolution frequencies etc It's been on my list of things to do to calibrate a PIC with GPS and allow generation or measurement of frequencies with extremely high accuracy. I've never been fortunate enough to own a $8000 HP frequency meter...

(edit) Wow I just downloaded the LEA-6S datasheet and there are some nice changes. My LEA-5S has a 1 Hz timebase pulse output with 50nS timing accuracy.

The new LEA-6S has a configurable timepulse of 0.1Hz to 1kHz with 30nS accuracy and the special unit LEA-6T has a second timepulse configurable from 1/60Hz to 10MHz ! Looks like the one module would be a GPS discliplined 10MHz freq reference! Very nice.

Howver i don't know if anyone sells a finished PCB with LEA-6T on it yet.

 

[MrAl]

Hi again Mr RB,

Do any of these versions output the pulse without having to upload some boot code to the device or board?
In other words, hook it up to +5vdc and see a 1 second pulse output, no microcontroller to help it boot up with config code.

 

[Mr RB]

Mr AL, based on my experience with the SmartGPS yes they do. All you need to do is connect 5v (or any reasonable DC supply) and put the antenna where it can find some satellites. I just have mine inside near a window, it was finding 5 to 7 satellites ok there.

Once power is connected it sends out serial GPS data (if you care) saying it has no satellites, then as sats are detected the GPS serial data updates the details. And the clock and location data are output and get more accurate.

After it gets a solid lock on enough sats (which can take 5 minutes) it starts sending out the 1pps pulse. This is a standard feature of the uBlox LEA series GPS modules.

The LEA-6S is a little unknown as it is "configurable" for different timepulse output frequencies. Probably under serial control. However as far as I know the 1pps output is an industry standrd so I would expect the LEA-6S to be like the LEA-5S and output the 1pps signal as default once it gets a sat lock.

I would suggest you look at the EasyGPS board above as it is more suited for being installed in equipment, with the SmartGPS you pay for more config options for general development purposes.

-------------------------

The "PIC measures its own xtal using GPS" project is now up on my web page;
High accuracy PIC timing systems

 

[Mr RB]

Mr Al, I had a read through the UBX protocol that is used to config the LEA modules (by serial). Both my LEA-5S and the new LEA-6S can be configured for a timepulse of 1ms to 60 seconds (1kHz to 1/60 Hz). The modules always boot up in default 1Hz mode. If they have a 3v battery connected like the SmartGPS they will retain config settings so you could set one to 1kHz output if needed and it would then boot up that way until the battery failed.

The UBX protocol is a little messy to change the timepulse frequency you need to send a 28 byte serial packet to the LEA module.

A minor update on my web page; I added a 4 sample rotating average to the new-sample. This gives 4 times less jitter now as it sums the latest 4 timer samples to get xtal freq instead of getting each new timer sample and *4 to convert to xtal freq.

I also found a fault, my XTAL_Calib2.c project always showed the xtal at 4 Parts Per Million slower than the first XTAL_Calib.c project! I pulled my hair out for 20 mins changing code until I realised it was the osc config. One project was set to XT osc, the other to HS osc. So now I know that with a 4MHz xtal it runs 4 PPM slower on HS osc than it does on XT osc...

And I recalibrated my trusty old freq meter to within about 0.25 PPM accuracy now by measuring the PIC xtal. The scope lead loaded the PIC osc down from 4000167 Mhz to 4000078 Mhz, and it was nice and stable so I calibrated my freq meter to 4000078 Mhz.

Now I need to build some nice little xtal ovens...