What is causing a negative bias in my super-sampling simulation?

Question

I'm developing a little test program to see if it is feasible to add noise to an ADC to gain oversampled bits. A little theory, before we begin. Nyquist's sampling theorem suggests that an increase in one bit of resolution requires four additional samples, and in general, n more bits requires 2^(n+1) more samples. I'm simulating a perfect 10 bit ADC which returns a value from 0..1023 monotonically and with no noise for an input from 0-2V.

To gain more bits, it is necessary to add randomly distributed noise (it doesn't have to be actually random, but it does have to appear random, like white noise.) The problem I'm having is although the resolution is increasing, the actual reading is offset by some small negative amount. Here's one sample output for an input of 1 volt (reference is 2 volts, so the count should be exactly half for an monotonic ADC):

10 bits:  512         volts:  1.0
11 bits:  1024        volts:  1.0
12 bits:  2046        volts:  0.9990234375
13 bits:  4093        volts:  0.999267578125
14 bits:  8189        volts:  0.999633789062
15 bits:  16375       volts:  0.999450683594
16 bits:  32753       volts:  0.999542236328
17 bits:  65509       volts:  0.999588012695
18 bits:  131013      volts:  0.999549865723
24 bits:  8384565     volts:  0.999518036842
28 bits:  134152551   volts:  0.999514393508

In fact, no matter how many times I run the simulation I'm always ending up with around ~0.9995, instead of 1; and the last value should be 134,217,728, not 134,152,551, which is about 65,771 out - or around 1/4 of the extra 18 bits of resolution (coincidence? I am diving by 4.) I suspect my PRNG is biased in some way, but I am using the default Mersenne Twister that comes with Python.

#!/usr/bin/python
# 
#  Demonstrates how oversampling/supersampling with noise can be used
#  to improve the resolution of an ADC reading.
#
#  Public domain.
#

import random, sys

volts = 1.000
reference = 2.000
noise = 0.01
adc_res = 10

def get_rand_bit():
    return random.choice([-1, 1])

def volts_with_noise():
    if get_rand_bit() == 1:
        return volts + (noise * random.random() * get_rand_bit())
    else:
        return volts

def sample_adc(v):
    # Sample ADC with adc_res bits on given voltage.
    frac = v / reference
    frac = max(min(frac, 1.0), 0.0) # clip voltage
    return int(frac * (2 ** adc_res))

def adc_do_no_noise_sample():
    return sample_adc(volts)

def adc_do_noise_sample(extra_bits_wanted):
    # The number of extra samples required to gain n bits (according to 
    # Nyquist) is 2^(n+1). So for 1 extra bit, we need to sample 4 times.
    samples = 2 ** (extra_bits_wanted + 1)
    print "Sampling ", samples, " times for ", extra_bits_wanted, " extra bits."
    # Sample the number of times and add the totals.
    total = 0
    for i in range(samples):
        if i % 100000 == 99999:
            print float(i * 100) / samples
            sys.stdout.flush()
        total += sample_adc(volts_with_noise())
    # Divide by two (to cancel out the +1 in 2^(n+1)) and return the integer part.
    return int(total / 2)

def convert_integer_to_volts(val, num_bits):
    # Get a fraction.
    frac = float(val) / (2 ** num_bits)
    # Multiply by the reference.
    return frac * reference

if __name__ == '__main__':
    # First test: we want a 10 bit sample.
    _10_bits = adc_do_no_noise_sample()
    # Next, create additional samples.
    _11_bits = adc_do_noise_sample(1)
    _12_bits = adc_do_noise_sample(2)
    _13_bits = adc_do_noise_sample(3)
    _14_bits = adc_do_noise_sample(4)
    _15_bits = adc_do_noise_sample(5)
    _16_bits = adc_do_noise_sample(6)
    _17_bits = adc_do_noise_sample(7)
    _18_bits = adc_do_noise_sample(8)
    _24_bits = adc_do_noise_sample(14)
    _28_bits = adc_do_noise_sample(18)
    # Print results both as integers and voltages.
    print "10 bits: ", _10_bits, "  volts: ", convert_integer_to_volts(_10_bits, 10)
    print "11 bits: ", _11_bits, "  volts: ", convert_integer_to_volts(_11_bits, 11)
    print "12 bits: ", _12_bits, "  volts: ", convert_integer_to_volts(_12_bits, 12)
    print "13 bits: ", _13_bits, "  volts: ", convert_integer_to_volts(_13_bits, 13)
    print "14 bits: ", _14_bits, "  volts: ", convert_integer_to_volts(_14_bits, 14)
    print "15 bits: ", _15_bits, "  volts: ", convert_integer_to_volts(_15_bits, 15)
    print "16 bits: ", _16_bits, "  volts: ", convert_integer_to_volts(_16_bits, 16)
    print "17 bits: ", _17_bits, "  volts: ", convert_integer_to_volts(_17_bits, 17)
    print "18 bits: ", _18_bits, "  volts: ", convert_integer_to_volts(_18_bits, 18)
    print "24 bits: ", _24_bits, "  volts: ", convert_integer_to_volts(_24_bits, 24)
    print "28 bits: ", _28_bits, "  volts: ", convert_integer_to_volts(_28_bits, 28)

I'd appreciate any suggestions on this. My plan is eventually to take this to a low cost microcontroller to implement a high resolution ADC. Speed will be fairly important, so I will probably be using an LFSR to generate PRNG bits, which won't be half as good as a Mersenne twister but should be good enough for most uses, and hopefully good enough for this.

Answer 1

In sample_adc(..) you probably want to round instead of truncating (systematically round towards negative infinity), i.e. do:

return int(frac * (2 ** adc_res) + 0.5)

instead of

return int(frac * (2 ** adc_res))

Then the deviations from one are not always on the same side:

10 bits:  512   volts:  1.0
11 bits:  1025   volts:  1.0009765625
12 bits:  2046   volts:  0.9990234375
13 bits:  4100   volts:  1.0009765625
14 bits:  8196   volts:  1.00048828125
15 bits:  16391   volts:  1.00042724609
16 bits:  32784   volts:  1.00048828125
17 bits:  65528   volts:  0.999877929688
18 bits:  131111   volts:  1.00029754639
24 bits:  8388594   volts:  0.99999833107
28 bits:  134216558   volts:  0.999991282821

Although to check the bias one would e.g. call adc_do_noise_sample(..) e.g. 10'000 times (for each resolution) and calculate the mean bias and the uncertainty on this mean (and check whether how compatible it is with zero).