Что вызывает отрицательный уклон в моей симуляции с суперсэмплированием?

#python #algorithm #statistics

#python #алгоритм #Статистика

Вопрос:

Я разрабатываю небольшую тестовую программу, чтобы посмотреть, возможно ли добавить шум в АЦП для получения битов с избыточной выборкой. Небольшая теория, прежде чем мы начнем. Теорема Найквиста о дискретизации предполагает, что увеличение разрешения на один бит требует четырех дополнительных выборок, и в общем случае для увеличения количества n бит требуется на 2 ^ (n 1) больше выборок. Я моделирую идеальный 10-разрядный АЦП, который монотонно и без шума возвращает значение от 0 .. 1023 для ввода от 0 до 2V.

Чтобы получить больше битов, необходимо добавить случайно распределенный шум (он не обязательно должен быть на самом деле случайным, но он должен казаться случайным, как белый шум.) Проблема, с которой я сталкиваюсь, заключается в том, что, хотя разрешение увеличивается, фактическое значение смещается на некоторую небольшую отрицательную величину. Вот один пример вывода для входа в 1 вольт (эталонное значение равно 2 вольтам, поэтому для монотонного АЦП количество должно быть ровно вдвое меньше):

 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
  

На самом деле, независимо от того, сколько раз я запускаю симуляцию, я всегда получаю около ~ 0.9995 вместо 1; и последнее значение должно быть 134 217 728, а не 134 152 551, что составляет около 65 771 из — или около 1/4 дополнительных 18 бит разрешения (совпадение? Я погружаюсь на 4.) Я подозреваю, что мой PRNG каким-то образом предвзят, но я использую стандартный Mersenne Twister, который поставляется с 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)
  

Я был бы признателен за любые предложения по этому поводу. В конечном итоге я планирую перенести это на недорогой микроконтроллер для реализации АЦП с высоким разрешением. Скорость будет довольно важна, поэтому я, вероятно, буду использовать LFSR для генерации битов PRNG, которые будут и вполовину не так хороши, как Mersenne twister, но должны быть достаточно хороши для большинства применений, и, надеюсь, достаточно хороши для этого.

Комментарии:

1. Не могли бы вы добавить подсказку о значении «ADC»?

2. Я также немного не уверен в ваших намерениях. «Разрешение» для меня означает «мельчайшие различимые детали». Если у вас есть подробности ниже разрешения, вам не хватает этой информации. Как добавление шума может добавить эту информацию обратно? Вы используете какое-то другое значение «разрешения»?

3. @Svante — Аналого-цифровой преобразователь.

4. @Svante, источник напряжения можно считать бесконечным по разрешению. Добавление шума приводит к тому, что значения в АЦП немного скачут, что увеличивает разрешение. Это трудно объяснить без знаний в области теории информации (которой у меня на самом деле нет, но она почерпнута в основном из книг и веб-сайтов).

5. не то, чтобы это, похоже, устраняло / добавляло предвзятость, но ваша реализация volts_with_noise() выглядит немного забавно для меня: в 50% выборок будет добавлен нулевой шум. Это специально?

Ответ №1:

В sample_adc(..) вы, вероятно, хотите округлять вместо усечения (систематически округлять до отрицательной бесконечности), т. Е. делать:

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

вместо

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

Тогда отклонения от единицы не всегда находятся на одной стороне:

 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
  

Хотя для проверки смещения можно, например, вызвать adc_do_noise_sample(..) , например, 10 000 раз (для каждого разрешения) и вычислить среднее смещение и неопределенность в этом среднем (и проверить, насколько оно совместимо с нулем).

Ответ №2:

Похоже, что ваш источник шума смещен. Я не знаком с Mersenne Twister, но я точно знаю, что генераторы псевдослучайного шума LFSR всегда имеют небольшое смещение. Вы можете сделать смещение сколь угодно малым, увеличив длину LFSR, но оно всегда будет на месте.

Комментарии:

1. Я провел быстрый тест, взяв среднее значение из 100 000 чисел из моего генератора случайных чисел, и оно составило от -0.0001 до 0.0001, так что смещение совсем небольшое (примерно то, что можно было бы предсказать). Проблема, похоже, заключается в очень постоянном уклоне в сторону ~ 0.9995.

Ответ №3:

Я не знаком с этой областью, но если вы используете python2, вы можете страдать от непреднамеренного целочисленного деления. Есть несколько подходов для решения этой проблемы.

 >>> 10 / 3
3
>>> 10 * 1.0 / 3
3.3333333333333335
>>> from __future__ import division
>>> 10 / 3
3.3333333333333335
>>> 10 // 3
3