Can boxcar averaging cause rounding error?

[Circlotron]

Got a setup where incoming pulses of variable frequency cause either a high or a low port pin output depending on whether the input frequency is above or below a certain frequency. Just sampling the counter, resetting it, then comparing it to a fixed value. No big deal. Without averaging it works dead nuts on. Perfect.

But if I boxcar average the last sixteen 16 bit counter samples by adding them up then right shifting this total four times, the frequency trip threshold is very slightly low. Instead of tripping at 120Hz it trips at something like 119.5Hz IIRC. I should run it again and get an exact figure. Are the four LS bits that are lost in the right shift the cause of this? It's no big deal, just want to know what's going on.

[Phoenix]

A right shift is a division with floor rounding. In order to correctly round you need to add half the divider to the start value.

Code: [Select]

value = (value+8)>>4

[Circlotron]

A right shift is a division with floor rounding. In order to correctly round you need to add half the divider to the start value.

Code: [Select]

value = (value+8)>>4

That is super interesting! I didn't know that. Thanks!

[newbrain]

A couple of notes:

1. I prefer to be explicit when coding (say what you mean, not what you think the compiler should be doing), so I'd use /16 rather than >>4 when I mean math and not bit manipulation. If the compiler is not a toy, it will use a shift even without optimizations (but see below). It's also easier to update, should the divisor change.

2. One of the improvement of C99 over C89 is the clear definition of integer division: a / will always round towards 0 (so -15 / 4 == -3) while a shift with a negative first operand is implementation defined. Other languages might differ: (ANSI) FORTH, e.g., rounds towards -∞.

3. Most implementations will perform an arithmetic shift on a signed operand, so the result will be rounded towards -∞ (so -15 >> 2 == -4), but the standard does not prevent using a logical shift (same as with an unsigned argument, so -15 >> 2 == 1'073'741'820, with surprising results).

Points 2 and 3 above are a bit moot in this case, as the OP is averaging a frequency - I don't expect it to take a negative value in this context.

[Circlotron]

Time for an update.

I should have mentioned, this thing is written in assembly, (NXP S08 micro) so the logical shift right isn't doing anything unexpected.

I am looking for 120Hz so 8,33333 mS.
 
1.25MHz counter so 8.33333 mS = 10,417 closest integer count for comparison.

I am polling the input pin with a 500nS loop so there is a small amount of latency. As a result, I am judging the threshold point as when the output pin dithers with approximately a 50% duty cycle.

No average = 119.99 Hz
Average 16 = 119.085 Hz
Average 16 with 8 added to sample = 119.175 Hz

So it is better, but not much. Hmmm...

[mtwieg]

Accumulating the 16 samples first and then dividing by 16 (with the +8 to ensure rounding) should also work, assuming it doesn't overflow your data type.

If it's only 16 samples, should be straightforward to look at the raw samples (as integers, not Hz) and check things by hand.

[Circlotron]

Accumulating the 16 samples first and then dividing by 16 (with the +8 to ensure rounding) should also work, assuming it doesn't overflow your data type.

I've done exactly that. I've got 24 bits to contain the total, and each sample can't exceed #46875 so there's plenty of room. Here's the code. I've cut some irrelevant parts out. Don't want to give things away. Besides, I'm no hotshot coder.

Code: [Select]

********************************************************************************
MAIN_LOOP LDHX #1250 ;3 1mS delay to prevent
NO_RETRIG AIX #-1 ;2 immediate retrigger
CPHX #0 ;3 8 cyc 0.8uS loop
BNE NO_RETRIG ;3

WAIT_TRIG BRSET TRIG,PORTD,WAIT_TRIG ;trigger on low pulse
LDHX TCNTH ;get counter
STA TCNTH ;reset counter
AIX #8 ;add 8 to sample
STHX TEMP_CTR_MSB

JSR PUSH_DOWN ;push old counter readings down one place
BRSET TOF,TSC,TOO_SLOW ;if under 37.5mS and set flag
BCLR SLOW,FLAGS

LDHX TEMP_CTR_MSB ;latest count
CONTINUE STHX MSB_1 ;update to latest
BCLR TOF,TSC ;just in case
JSR ADD_UP ;add up last 16 counter periods
JSR DIV_16 ;average = MSB_TOT:LSB_TOT
BRSET SLOW,FLAGS,NO_AVG ; > 37.5mS so don't use average
*********************************************
LDHX MSB_TOT ;averaged figure
CPHX #10417 ;120Hz, 8.333mS
BHS NO_BOOST
BSET 4,PORTB ;hi output
JMP MAIN_LOOP
*********************************************
TOO_SLOW LDHX #46875
BSET SLOW,FLAGS
JMP CONTINUE
*********************************************
NO_AVG MOV VOLTS_14,PORTC ;So first slow pulse won't have
JMP MAIN_LOOP ;to wait for 16 averages.
*********************************************
NO_BOOST BCLR 4,PORTB ;low output
JMP MAIN_LOOP
********************************************************************************
;push old counter readings down one place

PUSH_DOWN LDHX MSB_15
STHX MSB_16

LDHX MSB_14
STHX MSB_15

LDHX MSB_13
STHX MSB_14

LDHX MSB_12
STHX MSB_13

LDHX MSB_11
STHX MSB_12

LDHX MSB_10
STHX MSB_11

LDHX MSB_9
STHX MSB_10

LDHX MSB_8
STHX MSB_9

LDHX MSB_7
STHX MSB_8

LDHX MSB_6
STHX MSB_7

LDHX MSB_5
STHX MSB_6

LDHX MSB_4
STHX MSB_5

LDHX MSB_3
STHX MSB_4

LDHX MSB_2
STHX MSB_3

LDHX MSB_1
STHX MSB_2
RTS
*************************************
ADD_UP CLR     LSB_TOT
                CLR     MSB_TOT
                CLR     EXT_TOT

                LDA     LSB_TOT
                ADD     LSB_1
                STA     LSB_TOT
                LDA     MSB_TOT
                ADC     MSB_1
                STA     MSB_TOT
JSR EXT_ADD

                LDA     LSB_TOT
                ADD     LSB_2
                STA     LSB_TOT
                LDA     MSB_TOT
                ADC     MSB_2
                STA     MSB_TOT
JSR EXT_ADD

                LDA     LSB_TOT
                ADD     LSB_3
                STA     LSB_TOT
                LDA     MSB_TOT
                ADC     MSB_3
                STA     MSB_TOT
JSR EXT_ADD

                LDA     LSB_TOT
                ADD     LSB_4
                STA     LSB_TOT
                LDA     MSB_TOT
                ADC     MSB_4
                STA     MSB_TOT
JSR EXT_ADD

                LDA     LSB_TOT
                ADD     LSB_5
                STA     LSB_TOT
                LDA     MSB_TOT
                ADC     MSB_5
                STA     MSB_TOT
JSR EXT_ADD

                LDA     LSB_TOT
                ADD     LSB_6
                STA     LSB_TOT
                LDA     MSB_TOT
                ADC     MSB_6
                STA     MSB_TOT
JSR EXT_ADD

                LDA     LSB_TOT
                ADD     LSB_7
                STA     LSB_TOT
                LDA     MSB_TOT
                ADC     MSB_7
                STA     MSB_TOT
JSR EXT_ADD

                LDA     LSB_TOT
                ADD     LSB_8
                STA     LSB_TOT
                LDA     MSB_TOT
                ADC     MSB_8
                STA     MSB_TOT
JSR EXT_ADD

                LDA     LSB_TOT
                ADD     LSB_9
                STA     LSB_TOT
                LDA     MSB_TOT
                ADC     MSB_9
                STA     MSB_TOT
JSR EXT_ADD

                LDA     LSB_TOT
                ADD     LSB_10
                STA     LSB_TOT
                LDA     MSB_TOT
                ADC     MSB_10
                STA     MSB_TOT
JSR EXT_ADD

                LDA     LSB_TOT
                ADD     LSB_11
                STA     LSB_TOT
                LDA     MSB_TOT
                ADC     MSB_11
                STA     MSB_TOT
JSR EXT_ADD

                LDA     LSB_TOT
                ADD     LSB_12
                STA     LSB_TOT
                LDA     MSB_TOT
                ADC     MSB_12
                STA     MSB_TOT
JSR EXT_ADD

                LDA     LSB_TOT
                ADD     LSB_13
                STA     LSB_TOT
                LDA     MSB_TOT
                ADC     MSB_13
                STA     MSB_TOT
JSR EXT_ADD

                LDA     LSB_TOT
                ADD     LSB_14
                STA     LSB_TOT
                LDA     MSB_TOT
                ADC     MSB_14
                STA     MSB_TOT
JSR EXT_ADD

                LDA     LSB_TOT
                ADD     LSB_15
                STA     LSB_TOT
                LDA     MSB_TOT
                ADC     MSB_15
                STA     MSB_TOT
JSR EXT_ADD

                LDA     LSB_TOT
                ADD     LSB_16
                STA     LSB_TOT
                LDA     MSB_TOT
                ADC     MSB_16
                STA     MSB_TOT
JSR EXT_ADD
RTS
*************************************
EXT_ADD LDA     EXT_TOT
                ADC     #0
                STA     EXT_TOT
RTS
*************************************
DIV_16 MOV #4,LOOP_CTR
DIV_LOOP LSR     EXT_TOT
                ROR     MSB_TOT
                ROR     LSB_TOT
DBNZ LOOP_CTR,DIV_LOOP

RTS
********************************************************************************
********************************************************************************
                ORG     $FFFE
                FDB     START
********************************************************************************

















[Nominal Animal]

Shouldn't you set LSB_1 at some point, too?

[Circlotron]

Shouldn't you set LSB_1 at some point, too?

It is inherent. The upper and lower bytes of H:X are stored at MSB_x:LSB_x.  Provided that the latter are at adjacent memory locations, which they are.

[Circlotron]

Shouldn't you set LSB_1 at some point, too?

It gets set at the CONTINUE label either to the latest sampled counter value or to a fixed maximum value of #46875.  The latter is a consequence of the timer overflow flag being set so it goes through the TOO_SLOW label.

[Nominal Animal]

It is inherent.

What architecture and instruction set is this?

Wouldn't it make more sense to use a 24-bit accumulator you add the result to, instead of a 16 byte buffer?  Or do you really need the update rate of a sliding window?

[Circlotron]

NXP S08 architecture. Specifically 9S08FL micro.
Doing it in assembler so what is a different way to have a 24 bit accumulator other than what I have done?

The averaged value I am getting seems to be about 80 counts edit -> *less* than expected.

[KE5FX]

Normally rounding error just shows up as DC offset, not frequency error.  Could there be an aliasing effect between the clock frequency and the frequency being measured?  To rule that out, try some other frequencies and see how they influence the magnitude and sign of the error.

Hysteresis is another thing to consider.  A ZCD typically shouldn't actually be implemented as a "zero crossing detector," since that's where the noise is.  Instead look for falling and rising edges separately with nonzero threshold voltages, if you're not doing so already (I didn't read the code.)

[Circlotron]

 No ZCD involved. Nothing to do with power line frequency either. I’m looking at the repetition rate of a 20uS duration pulse.

[Circlotron]

Hey, I got it!

In the RAM equates list the MS byte for each 16 bit sample is listed first because whenever the H:X register is loaded from or saved to memory you just specify the MS byte location and the LS byte is the next one along. e.g. LDHX MSB_15 and it grabs MSB_15 and LSB_15.

Where I messed up was the 3 byte totals RAM. I had it listed as LS, MS, EXT(ra) so after the  samples were all added up and LSR'd four places to the right, H:X instead of being loaded with MS_TOT:LS_TOT was being loaded with MS_TOT:EXT_TOT, and seeing these three had been LSR'd that meant EXT_TOT (the wrongful LS byte being loaded into X) was empty so that would account for the approximately 80 lower count than expected.

Code: [Select]

MSB_15 EQU $5C
LSB_15 EQU $5D
MSB_16 EQU $5E
LSB_16 EQU $5F
LSB_TOT EQU $60
MSB_TOT EQU $61
EXT_TOT EQU $62
FLAGS EQU $63
Something else, it is actually more accurate without adding 8 to each sample as discussed earlier.

[mikerj]

Something else, it is actually more accurate without adding 8 to each sample as discussed earlier.

I'm not surprised, you aren't supposed to be adding 8 to each sample!  You add 8 to the total of 16 samples before dividing by 16.  This gives you rounding to the nearest integer rather than just truncating all fractional bits.

[jpanhalt]

@Circlotron
Can you post the 16 values used in post #5? 

John

[S. Petrukhin]

A right shift is a division with floor rounding. In order to correctly round you need to add half the divider to the start value.

Code: [Select]

value = (value+8)>>4

Is this true for negative values? 

[Circlotron]

@Circlotron
Can you post the 16 values used in post #5? 

John

They would be in the region of 10,417 if the measure pulse rate was 120 Hz. Lower frequency = longer time between pulses so greater count value. Counter rate = 1.25MHz.

[gf]

A right shift is a division with floor rounding. In order to correctly round you need to add half the divider to the start value.

Code: [Select]

value = (value+8)>>4

Is this true for negative values? 

It fails for very large positive values. For eample, if type is int, then values >= MAX_INT - 7 lead to a negative result.

[jpanhalt]

I was referring to this:

No average = 119.99 Hz
Average 16 = 119.085 Hz
Average 16 with 8 added to sample = 119.175 Hz

An alternative way to round, which I am sure you saw, is to check the carry bit after the last rotation (e.g., PIC equivalent is decrease counter and skip if zero) and add 1, if set.  My experience is limited to PIC Assembly, so the instructions and execution times are a bit different.

John

[S. Petrukhin]

Most often I use the moving average method according to the formula OvrValue:=(1-k)*OvrValue+k*CurValue where k is usually 0.2 or 0.1.

This smooths out fluctuations quite well, recovers quickly and does not require accumulation.
If using real mathematics is inconvenient, you can make a shift towards integer values.

[Circlotron]

I was referring to this:

No average = 119.99 Hz
Average 16 = 119.085 Hz
Average 16 with 8 added to sample = 119.175 Hz

John

If I understand you, the 16 values would be 1.25 MHz x the period of each of those frequencies.

[Circlotron]

An alternative way to round, which I am sure you saw, is to check the carry bit after the last rotation (e.g., PIC equivalent is decrease counter and skip if zero) and add 1, if set.  My experience is limited to PIC Assembly, so the instructions and execution times are a bit different.

John

Hey, That's interesting! Thanks for that. It took a moment to sink in (11:20pm here) but yeah, that's cool!

[S. Petrukhin]

No average = 119.99 Hz
Average 16 = 119.085 Hz
Average 16 with 8 added to sample = 119.175 Hz

Just add up the current values as they are, and then round up the result.
The sum of rounded values will never equal the rounding of the average value.