Python
June 2009 17

Случайные числа из звуковой карты

Многие когда-либо интересовались случайными числами. Хочу поделиться моими экспериментами по получению истинно случайных чисел с помощью «аппаратного генератора» встроенного в практически любой компьютер — звуковой карты.

При подготовке материала, я переписал свой старый Си код на Питоне, поэтому данный опус также является примером по использованию Windows DLL из Питона с использованием стандартной библиотеки ctypes.

В конце статьи сравниваются данные полученные от двух звуковых карт Realtek и Audigy 2, приведены результаты статистических тестов на случайность.

UPD Исправил пропавшие в коде нули, которые съело НЛО.

Скучное теоретическое введение


Практически все языки программирования предоставляют несколько функций генерации так называемых ПСЧ Псевдослучайных Чисел. ПСЧ не являются случайными в математическом смысле слова, они получаются по некоторому известному алгоритму и, зная исходные параметры алгоритма, полученную последовательность можно всегда повторить снова. Для многих задач качественный генератор ПСЧ вполне может заменить истинно случайные числа. Компьютерные игры, моделирование процессов, интегрирование Монте-Карло, генетические алгоритмы… список задач, для которых достаточно хорошего ПСЧ, можно продолжать долго.

С другой стороны, есть ограниченный круг проблем, для которых важна истинная случайность полученной последовательности. Главный пример — криптография. Именно в контексте шифрования люди всегда интересовались вопросом получения истинно случайных чисел. Простейший пример — единственный шифр, для которого доказана абсолютная криптографическая стойкость, Шифр Вернама (англ. One-time pad) требует для ключа последовательность истинно случайных чисел равную по длине секретному сообщению. Для обеспечения криптостойкости случайные данные используемые для генерации ключей (будь то шифр Вернама, AES, RSA или что-то еще) никогда не должны использоваться повторно. Что приводит нас к вопросу о поиске надежного источника случайных чисел.

Звуковая карта, в отличие от большинства остальных компонентов компьютера, имеет в себе не только цифровую часть, но и аналоговую.

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

Нас интересует пункт 2. Как известно, любой аналоговый электрический сигнал неизбежно содержит шум, компоненты шума можно примерно разделить на несколько категорий:
  • радиопомехи и наводки от соседних устройств и радиоэфира (слышали когда-нибудь как плохо экранированный усилитель ловит радио?)
  • помехи электропитания (принципиально тоже можно отнести к радиопомехам)
  • тепловой шум случайного движения электронов в компонентах электрической схемы

Если случайность помех и питания вопрос спорный, то третий тип шума чисто квантовый и истинно случайный.

Собственно вопрос: насколько этот шум проникает в цифровую часть звуковой карты? Опыт показывает, что да — проникает, его можно оцифровать и сохранить.

В режиме записи 16 бит случайную информацию несет только самый младший бит из каждого сэмпла, в 24 битном режиме — несколько младших, но надежней всего всегда брать только один самый младший бит. Как это сделать дальше и пойдет речь на примере Python программы для Windows.

Кому не интересен Питон: анализ результатов и выводы в конце, после описания программы.


Запись звука в Windows


Простейшим способом записать звук в Windows является использование интерфейсов Waveform Audio из библиотеки winmm.dll. Стандартной библиотеки для работы со звуком в Питоне нет, поэтому мы будет пользоваться библиотекой ctypes, которая предоставляет интерфейс к обычным DLL.

Импортируем библиотеки sys (для доступа к стандартному выводу) и time (для доступа к функции sleep). Также импортируем все имена из ctypes.

import sys
import time
from ctypes import *


Для С-функций возвращающих целое число (обычно код ошибки) Питон позволяет задать функцию, которая будет проверять код ошибки. Мы воспользуемся этим, чтобы добавить минимальный контроль ошибок в программе: если С-функция возвратит ошибку, MMSYSERR_NOERROR вызовет Питон-исключение и в консоли уже будет видно, где именно проблема.

Затем идет цикл, который для каждой функции библиотеки winmm.dll (Питон объект windll.winmm импортирован из ctypes) из списка, мы создаем переменную в текущем контексте vars(), это позволит позднее обращаться к функции просто по имени (waveInOpen вместо windll.winmm.waveInOpen). Также мы присваиваем возвращаемому типу нашу «контролирующую» функцию MMSYSERR_NOERROR.

def MMSYSERR_NOERROR(value):
    if value ! 0 :
        raise Exception("Error while running winmm.dll function", value)
    return value
 
for funcname in ["waveInOpen""waveInPrepareHeader",
                 "waveInAddBuffer""waveInStart",
                 "waveInStop""waveInReset",
                 "waveInUnprepareHeader""waveInClose"]:
    vars()[funcname] = windll.winmm[funcname]
    vars()[funcname].restype = MMSYSERR_NOERROR


Определим необходимые для работы с Windows Audio Си-структуры. Класс структуры должен наследоваться от класса Structure импортированного из ctypes, и должен содержать поле _fields_, перечисляющее имена и Си-типы элементов структуры. Классы для Си-типов мы импортировали из ctypes, их названия говорят сами за себя: c_int, c_uint и тд

Первая структура WAVEFORMATEX содержит информацию о формате звуковых данных. Без параметров конструктор создаст структуру с типичными значениями для большинства звуковых карт: 16бит 48кГц моно.

Вторая WAVEHDR описывает буфер для звуковых данных. В качестве параметра конструктор требует объект типа WAVEFORMATEX и выделяет буфер способный хранить 1 секунду аудио данного формата. Массив С-символов создается функцией create_string_buffer.

class WAVEFORMATEX(Structure):
    WAVE_FORMAT_PCM = 1
    _fields_ = [("wFormatTag", c_ushort),
                ("nChannels", c_ushort),
                ("nSamplesPerSec", c_uint),
                ("nAvgBytesPerSec", c_uint),
                ("nBlockAlign", c_ushort),
                ("wBitsPerSample", c_ushort),
                ("cbSize", c_ushort)]
 
    def __init__(self, samples=48000, bits=16, channels=1):
        self.wFormatTag = WAVEFORMATEX.WAVE_FORMAT_PCM
        self.nSamplesPerSec = samples
        self.wBitsPerSample = bits
        self.nChannels = channels
        self.nBlockAlign = self.nChannels*self.wBitsPerSample/8
        self.nAvgBytesPerSec = self.nBlockAlign*self.nSamplesPerSec
        self.cbSize =  0
 
class WAVEHDR(Structure):
    _fields_ = [("lpData", POINTER(c_char)),
                ("dwBufferLength", c_uint),
                ("dwBytesRecorded", c_uint),
                ("dwUser", c_uint)# User data dword or pointer
                ("dwFlags", c_uint),
                ("dwLoops", c_uint),
                ("lpNext", c_uint)# pointer reserved
                ("reserved", c_uint)] # pointer reserved
    def __init__(self, waveformat):
        self.dwBufferLength = waveformat.nAvgBytesPerSec
        self.lpData = create_string_buffer('\000' * self.dwBufferLength)
        self.dwFlags =  0


Далее мы создаем объект waveFormat. И три буфера для звуковых данных.

К сожалению с большинством драйверов winmm.dll (тк это достаточно старый интерфейс) не позволяет оцифровывать точнее 16 бит даже если аудиокарта это поддерживает. Мне известна только одна карта: SB Live24bit, которая это могла. Сейчас её под рукой нет, а есть Audigy 2 Notebook, которая пишет 24 бит только в DirectX или ASIO. Поэтому сегодняшний пример рассчитан на 16 бит (место которое надо изменить для 24 бит помечено комментарием, на случай если ваша карта это поддерживает).

waveFormat = WAVEFORMATEX(samples=48000,bits=16)
waveBufferArray = [WAVEHDR(waveFormat) for i in range(3)]


Следующаяя функция основная в программе, она будет вызываться Windows, когда winmm.dll заполнит один из наших буферов.

Тк это Си-callback, сначала необходимо создать класс. Это делается функцией WINFUNCTYPE из ctypes, аргументы: возвращаемый тип, типы аргументов. Тк наша функция ничего не должна возвращать, первый аргумент None, остальные согласно MSDN. Обратите внимание на аргумент POINTER(c_uint) — это указатель на пользовательские данные, в примере это просто число, но там может быть что угодно, например класс указывающий куда писать данные или сколько данных необходимо и тд. POINTER(WAVEHDR) это указатель на буфер с данными.

Параметр uMsg указывает на причину вызова, нас интересует MM_WIM_DATA — доступны аудиоданные.

WRITECALLBACK = WINFUNCTYPE(None, c_uint, c_uint, POINTER(c_uint), POINTER(WAVEHDR), c_uint)
def pythonWriteCallBack(HandleWaveIn, uMsg, dwInstance, dwParam1, dwParam2):
    MM_WIM_CLOSE = 0x3BF
    MM_WIM_DATA = 0x3C0
    MM_WIM_OPEN = 0x3BE
    if   uMsg == MM_WIM_OPEN:
        print "Open handle =", HandleWaveIn
    elif uMsg == MM_WIM_CLOSE:
        print "Close handle =", HandleWaveIn
    elif uMsg == MM_WIM_DATA:
        #print "Data handle =", HandleWaveIn
        wavBuf = dwParam1.contents
        if wavBuf.dwBytesRecorded >  0 :
            bits = [ord(wavBuf.lpData[i]) & 1 for i in range( 0 ,wavBuf.dwBytesRecorded,2)]
            # для 24 бит: заменить в конце 2 на 3 в предыдущей строке
            bias = [bits[i] for i in range( 0 ,len(bits),2) if bits[i] != bits[i+1]]
            bytes =  [chr(reduce(lambda v,b:v<<1|b,bias[i-8:i], 0 )) for i in range(8,len(bias),8)]
            rndstr = ''.join(bytes)
            #print bytes,
            sys.stdout.write(rndstr)
        if wavBuf.dwBytesRecorded == wavBuf.dwBufferLength:
            waveInAddBuffer(HandleWaveIn, dwParam1, sizeof(waveBuf))
        else:
            print "Releasing one buffer from", dwInstance[ 0 ]
            dwInstance[ 0 ]-=1
    else:
        raise "Unknown message"


Ключевые моменты:
bits = [ord(wavBuf.lpData[i]) & 1 for i in range( 0 ,wavBuf.dwBytesRecorded,2)]

Тк данные упакованы по 2 байта (3 байта в случае 24 бит) мы проходим по массиву с шагом 2: range(0,wavBuf.dwBytesRecorded,2) и отбираем младший бит младшего байта в паре: ord(wavBuf.lpData[i])&1. Результатом будет список bits младших битов каждого сэмпла.

Небольшое метематическое отступление:

Случайные данные могут иметь различное распределение, т.е. частоту выпадения того или иного числа. Например, если у нас последовательность бит 0000100000000, и положение единицы меняется случайно, это тоже случайная последовательность, но очевидно, что вероятность выпадения нуля гораздо выше чем единицы. Наиболее удобно т.н. «равномерное распределение» в котором вероятность появления нуля и единицы равна. Процедура приведения к равномерному распределению называется unbiasing. Наиболее простой метод это замена: 10 -> 1, 01 -> 0, 11 -> отбрасываем, 00 -> отбрасываем.

Unbiasing выполняет следующая строка
bias = [bits[i] for i in range( 0 ,len(bits),2) if bits[i] != bits[i+1]]

Проходим с шагом 2: range(0,len(bits),2), выкидывая одинаковые пары: bits[i] != bits[i+1], из оставшихся берем первый: bits[i].

Наконец, это выражение собирает наши биты в байты и сливает все в строку
 bytes =  [chr(reduce(lambda v,b:v<<1|b,bias[i-8:i], 0 )) for i in range(8,len(bias),8)]

Проходим с шагом 8, для каждых 8ми битов вызывается функция reduce(lambda v,b:v<<1|b,bias[i-8:i],0). Здесь использован элемент функционального программирования, анонимная функция(назовем её временно F) lambda v,b:v<<1|b, вызывается функцией reduce так: F(...F(F(0,bias[i-8]),bias[i-7]),...,bias[i-1]) — получается байт, который преобразуется в символ функцией chr.

Список байт преобразуется в строку, которая пишется в stdout, тк это двоичные данные лучше вывод перенаправить в файл. Если хотите писать в консоль, то лучше это делать через "print bytes,", тк при бинарном выводе в консоль Винда плохо ловит Ctrl+C.

В конце функции стоит проверка, полностью ли был заполнен буфер. Если да, то мы возвращаем его обратно системе для заполнения функцией waveInAddBuffer. Если нет, это значит что вывод данных остановлен (устройство закрыто), уменьшаем счетчик занятых буферов на единицу (счетчик хранится в пользовательских данных).

Создаем экземпляр С-функции класса WRITECALLBACK из нашей функции pythonWriteCallBack.

Далее, определив несколько полезных констант, открываем устройство WAVE_MAPPER (можно напрямую задать номер звковой карты, они нумеруются, начиная с нуля: 0, 1, 2) сначала с параметром WAVE_FORMAT_QUERY, чтобы проверить что формат поддерживается, затем с параметром CALLBACK_FUNCTION, указав нашу функцию и пользовательские данные (в нашем случае число ExitFlag).

Обратите внимание как передаются указатели из Питона в Си-функции, с помощью функции byref. Также мы передали указатель на наши «пользовательские данные» в byref(ExitFlag), который Windows будет передавать нашему callback в виде dwInstance при каждом событии (заполнения буфера например).

Затем мы вызываем waveInPrepareHeader для каждого созданного буфера и передаем их winmm.dll с помощью waveInAddBuffer. Наконец, вызов waveInStart(HandleWaveIn) дает команду начать ввод аудио. Конца мы ждем в цикле time.sleep(1).

Выход из программы осуществляется по перехвату комбинации клавиш Ctrl+C (Питон KeyboardInterrupt).
writeCallBack=WRITECALLBACK(pythonWriteCallBack)
try:
    ExitFlag = c_uint(3)
    HandleWaveIn = c_uint( 0 )
    WAVE_MAPPER = c_int(-1)
    WAVE_FORMAT_QUERY = c_int(1)
    CALLBACK_FUNCTION = c_int(0x30000)
 
    waveInOpen( 0 , WAVE_MAPPER, byref(waveFormat) 0 0 , WAVE_FORMAT_QUERY)
    waveInOpen(byref(HandleWaveIn), WAVE_MAPPER, byref(waveFormat), writeCallBack, byref(ExitFlag), CALLBACK_FUNCTION)
 
    for waveBuf in waveBufferArray:
        waveInPrepareHeader(HandleWaveIn, byref(waveBuf), sizeof(waveBuf))
        waveInAddBuffer(HandleWaveIn, byref(waveBuf), sizeof(waveBuf))
 
    waveInStart(HandleWaveIn)
 
    while 1:
        time.sleep(1)
 
except KeyboardInterrupt:
    waveInReset(HandleWaveIn)
 
    while 1:
        time.sleep(1)
        if ExitFlag.value ==  0 :
            break
 
    for waveBuf in waveBufferArray:
        waveInUnprepareHeader(HandleWaveIn, byref(waveBuf), sizeof(waveBuf))
 
    waveInClose(HandleWaveIn)


Чтобы отключить автоматическое преобразование конца строки, вызов скрипта следует осуществлять с параметром "-u":
c:\...\python.exe -u Скрипт.py > файл.rnd

Выводы



Записывать можно без внешнего сигнала, просто «тишину», младшие биты всё равно получают случайные значения. Делая запись без внешнего сигнала (источника), мы получаем характеристики шума естественно присутствующего в цепях аудиокарты. Другим вариантом может быть запись какого либо сигнала, тогда шум будет проявляться в ошибках оцифровки младших битов.

Ошибки оцифровки есть всегда, с другой стороны «нулевой» шум сильно зависит от физического устройства аудиокарты и может не проявляться на некоторых моделях. Далее обсуждаются результаты по анализу «нулевого» шума.

Оптимальные настройки микшера для записи, полученные экспериментальным путем: выбираем канал записи (не воспроизведения!) Line-In, громкость канала в максимум, все остальные каналы отключаем.

Канал Microphone показал худшее качество случайных чисел, объясняется это тем, что многие аудиокарты пытаются «улучшить» сигнал с микрофона и применяют различные цифровые фильтры. На одной из протестированных аудиокарт (а именно Realtek) микрофонный канал выдавал негодный для получения СЧ вывод. На Audigy 2 включение усиления микрофонного сигнала «Mic Boost +20DBA», также приводило к удалению случайной компоненты.

Протестировать качество полученных случайных чисел можно при помощи нескольких программ. Простейшая в использовании ent (скачать с http://www.fourmilab.ch/random/random.zip). Запускается из консоли
> type data.rnd | ent.exe
Где data.rnd двоичный файл со случайными данными (не забудьте закомментировать лишние print'ы в программе, они могут немного подпортить статистику). Оптимальный размер файла для теста ок. 500КБ. Программа подсчитывает несколько параметров:
  • Энтропию (для случайных данных: стремится к 8 — количество бит в байте)
  • Арифметическое среднее (для случайных данных: стремится к 127.5)
  • Число Пи методом Монте-Карло (для 500кБ данных ошибка порядка 0.1%)
  • Хи-Квадрат (в идеале лежит в пределах 10-90%)
  • Коэффициент (для случайных данных близок к 0)


Следует помнить, что критерии случайности — статистические по своей природе и применяются не к последовательности самой по себе, а к источнику. Поэтому для получения надежных результатов надо провести серию тестов на выборках одного источника. Наример 20 кусков по 500КБ. Если большинство прошло тест (ок. 90%), то источник случаен с определенной вероятностью.

К сожалению, 100%-ых критериев в статистике нет, т.к. с определенной вероятностью ваша звуковая карта может сгенерировать эту статью (собственно так я её и получил — шучу).

Также, я протестировал случайность «нулевого» вывода обеих звуковых карт с помощью более навороченного пакета программ от NIST (Национального Института Стандартов США). Обе карты показали хорошее качество случайных чисел.

Занимательно, что встроенное аудио от Realtek, дает чуть лучшую равномерность распределения, видимо из-за более низкого качества АЦП и высокой шумности. Audigy 2 была бы вне конкуренции в 24 битном режиме, который Реалтеком не поддерживается (в любом случае для него нужен DirectX, а DirectX в Питоне это отдельная история).

Ссылки на другие тестовые пакеты:
http://stat.fsu.edu/pub/diehard/
http://www.phy.duke.edu/~rgb/General/dieharder.php
http://csrc.nist.gov/groups/ST/toolkit/rng/documentation_software.html

Пример вывода программы ent на файле полученном записью тишины с линейного входа Audigy 2:
Entropy = 7.999609 bits per byte.

Optimum compression would reduce the size
of this 500000 byte file by 0 percent.

Chi square distribution for 500000 samples is 270.78, and randomly
would exceed this value 23.75 percent of the times.

Arithmetic mean value of data bytes is 127.5890 (127.5 = random).
Monte Carlo value for Pi is 3.139644559 (error 0.06 percent).
Serial correlation coefficient is 0.001109 (totally uncorrelated = 0.0).
+131
10.3k 73
Comments 72
Top of the day