Н.В. Клиначев

Программы генераторов псевдослучайных последовательностей для микроконтроллеров с целочисленным АЛУ

Рабочие файлы: [IQsin] [IQatan2] [IQsqrt] [IQlog] [IQgauss] [IQ-Математика]

Одна из задач, которую решает разработчик цифровых систем управления – это проверка качества работы технической системы в условиях воздействия помех. С этой целью, включают программный генератор шума и оптимизируют настройки системы управления. Требования к параметрам генераторов псевдослучайных последовательностей (ГПСЧ), используемых для решения таких задач, не высокие. Они могут иметь малый период, в пределах которого возможны циклические последовательности. Их распределение может нестрого совпадать с типовым. Но их программный код должен исполняться микроконтроллером без математического сопроцессора. И все промежуточные преобразования не должны существенно менять динамический диапазон чисел. Можно рекомендовать три генератора: IQnoise – ГПСЧ с равномерным распределением, IQgauss и IQgaussian – генераторы шума Гаусса с распределением близким к нормальному.

В библиотеке IQ-математики все генераторы шума используют регистр сдвига с линейной обратной связью (LFSR) [1]. Формально, LFSR – это псевдослучайный генератор последовательности битов. Но числа в таком регистре не повторяются на периоде – распределение равномерное. Программная реализация состоит из одной строки кода:

static unsigned long lfsr = 0x1;  // регистр сдвига != 0
lfsr = (lfsr & 0x1) ? (lfsr >> 1) ^ 0x829 : (lfsr >> 1);

где 0x829 – любое число из набора для 12-битного регистра:
829 834 83D 83E 84C 868 875 883 88F 891 89D 8A7 8AB 8B0 8B5
8C2 8D9 8EC 8EF 906 91B 91E 933 939 93F 95C 960 965 987 98E
990 99C 99F 9A6 9B8 9CC 9D1 9D4 A03 A18 A1B A27 A2E A33 A3A
A53 A56 A69 A87 A8E AA6 AC9 AE2 AEB AEE AF5 B04 B23 B2A B2C
B52 B5E B8A B8C BA1 BA2 BBA BC4 BD6 BD9 BDF BE0 C2B C2E C48
C4B C5C C77 C8D C9A CA0 CB2 CBD CC5 CD8 CDE CE4 CE7 CF3 D0D
D15 D19 D34 D45 D68 D70 D7A D85 D89 D8F
Источник: http://users.ece.cmu.edu/~koopman/lfsr/index.html

Для изменения разрядности регистра сдвига достаточно поменять число [2], кодирующее линейную обратную связь. Ограничение длины LFSR-регистра 12-тью битами определено, во-первых, разрядностью ЦАП-а микроконтроллера, во-вторых, требованием сокращения динамического диапазона при преобразовании чисел. В функции IQgaussian используется преобразование Бокса – Мюллера [3], в котором вычисляется натуральный логарифм от веса младшего разряда LFSR-регистра:

2 * _IQlog(_IQ(1/4096)) == 2 * _IQ(-8.315).

Т.о., при выполнении преобразований, число 2 * 8.315 – это предел для преобразуемых числовых значений. Для кодирования таких чисел, в 32-х битном целом под дробную часть можно отвести до 26 бит. Но поскольку используется функция IQsqrt, константа GL_Q должна быть чётной. Т.е. её максимальное значение равно 26. Минимальное – 16. Этот диапазон отвечает требованиям точности любых систем управления.

Чертёж 1

Листинг 1. Программный код dll-блока (javascript)

Ниже по тексту представлен листинг 2 с фрагментом кода h-файла библиотеки IQMath. В нём отсутствуют определения нескольких макросов (_IQ, _IQmpy, _IQcos), но на сайте имеются описания проектов [4], в которых используются эта библиотека. Фрагмент можно подключить к любой её версии.

/*******************************************************************************
* Function Name : _IQgauss
* Description :   Функция генерирует нормально распределённую псевдослучайную
*     последовательность чисел (шум Гаусса) со стандартными параметрами:
*     математическое ожидание == _IQ(0.0), дисперсия == _IQ(1.0).
*     Генератор использует метод, предполагающий суммирование нескольких
*     псевдослучайных равномерно распределенных чисел. Для генерации которых
*     используется 12-ти битный LFSR-регистр сдвига с линейной обратной связью
*     (см. числа: http://users.ece.cmu.edu/~koopman/lfsr/index.html).
*     Примечание. У хорошего ГПСЧ, дисперсия суммы 12 чисел == 1.
*                 У используемого ГПСЧ на LFSR, дисперсия суммы 8 чисел == 1.057
*     Рекомендуемое значение константы GL_Q от 16 до 26.
*     Время исполнения: код 3u125, функция 3u487, inline 3u070, без опт. 13u6
*******************************************************************************/
inline _iq _IQgauss(void);
#pragma optimize=speed       // 4 IAR
inline _iq _IQgauss(void) {
    static _iq rg0 = 0x1;
    unsigned char i = 0;
    _iq rnd8 = 0;
    for (i = 0; i < 16; i += 1) {
        rg0 = (rg0 & 1) ? (rg0 >> 1) ^ 0x829 : (rg0 >> 1); // LFSR 12bit
        rnd8 += rg0;                                       // 2 * rnd (0.0, 8.0)
    }
    rg0 = (rg0 & 1) ? (rg0 >> 1) ^ 0x829 : (rg0 >> 1);     // Увеличим период
    rnd8 = (rnd8 << (GL_Q - 12 - 1)) - _IQ(4.0);           // Вычтем среднее
    rnd8 = _IQmpy(rnd8, _IQ(1.057));                       // Приведем дисперсию
    return rnd8;
}

/*******************************************************************************
* Function Name : _IQsqrt
* Description :   Функция вычисления значения квадратного корня _IQsqrt(X).
*     Аргумент и возвращаемое значение - 32-х битные целые числа без знака.
*     Тип данных аргумента и возвращаемого значения - GL_Q. GL_Q - чётная.
*     Абсолютная погрешность вычисления результата - зависит от GL_Q.
*     http://model.exponenta.ru/k2/Jigrein/JS/fwlink.htm#B2C7
*     Время исполнения от значений аргументов не зависит - 1.4 uS.
*******************************************************************************/
inline _iq _IQsqrt(unsigned long val);
#pragma optimize=speed       // 4 IAR
inline _iq _IQsqrt(unsigned long val) {
    unsigned long temp = 0, g = 0;

    if (val >= (1 << 30)) { g = 1 << 15; val -= 1 << 30; }

    #define ITERATION_STEP(s) {                                                 \
        temp = (g << s) + (1 << ((s << 1) - 2));                                \
        if (val >= temp) { g += 1 << (s - 1); val -= temp; }                    \
    }

    ITERATION_STEP(15);
    ITERATION_STEP(14);
    ITERATION_STEP(13);
    ITERATION_STEP(12);
    ITERATION_STEP(11);
    ITERATION_STEP(10);
    ITERATION_STEP( 9);
    ITERATION_STEP( 8);
    ITERATION_STEP( 7);
    ITERATION_STEP( 6);
    ITERATION_STEP( 5);
    ITERATION_STEP( 4);
    ITERATION_STEP( 3);
    ITERATION_STEP( 2);

    #undef ITERATION_STEP

    temp = (g << 1) + 1;
    if (val >= temp) { g += 1; }
    return (_iq)(g << (GL_Q >> 1));
}

/*******************************************************************************
* Function Name : _IQlog
* Description :   Функция вычисляет натуральный логарифм числа.
*     Используется полином наилучшего приближения на интервале от 0.5 до 1.
*     Время исполнения: код , функция , inline , без опт.  us
*******************************************************************************/
// Вычисление натурального логарифма полином наилучшего приближения
inline _iq _IQlog(_iq x);
#pragma optimize=speed       // 4 IAR
inline _iq _IQlog(_iq x) {
    _iq n = 0;
    while (x >= _IQ(1.0)) { x = x >> 1; n += _IQ(0.6931471805599453); } // M_LN2
    while (x <  _IQ(0.5)) { x = x << 1; n -= _IQ(0.6931471805599453); } // M_LN2
    return x -= _IQ(1.0039), _IQmpy(x, (_IQ(0.891) - _IQmpy(x, _IQ(0.9472)))) + n;
}

/*******************************************************************************
* Function Name : _IQgaussian
* Description :   Функция генерирует нормально распределённую псевдослучайную
*     последовательность чисел (шум Гаусса) со стандартными параметрами:
*     математическое ожидание == _IQ(0.0), дисперсия == _IQ(1.0).
*     Генератор использует метод, предполагающий изменение распределения
*     псевдослучайных чисел (из равномерного в нормальное) преобразованием
*     Бокса – Мюллера:
*     -------- gauss1 = sin(2 * pi * rnd1) * sqrt(-2 * ln(rnd2)) --------
*     -------- gauss2 = cos(2 * pi * rnd1) * sqrt(-2 * ln(rnd2)) --------
*     Для генерации псевдослучайных чисел используется 12-ти битный
*     LFSR-регистр сдвига с линейной обратной связью
*     (см. числа: http://users.ece.cmu.edu/~koopman/lfsr/index.html).
*     Примечание. У хорошего ГПСЧ, дисперсия суммы 12 чисел == 1.
*                 У используемого ГПСЧ на LFSR, дисперсия суммы 8 чисел == 1.057
*     Константа GL_Q должны быть четной (см. _IQsqrt) от 16 до 26.
*     //2 * Math.log((1 << (GL_Q - 12)) / (1 << GL_Q)) == 2 * -8.314095497131347
*     Время исполнения: код , функция 3u6..4u79, inline 3u7..4u68, без опт. 7u84
*******************************************************************************/
inline _iq _IQgaussian(void);
#pragma optimize=speed       // 4 IAR
inline _iq _IQgaussian(void) {
    static _iq rg1 = 0x1, rg2 = 0x1;
    _iq rnd1 = 0, rnd2 = 0, hlp = 0;
    rg1 = (rg1 & 1) ? (rg1 >> 1) ^ 0x100D : (rg1 >> 1);     // LFSR 13bit
    rg1 = (rg1 & 1) ? (rg1 >> 1) ^ 0x100D : (rg1 >> 1);     // LFSR 13bit
    rg2 = (rg2 & 1) ? (rg2 >> 1) ^ 0x834  : (rg2 >> 1);     // LFSR 12bit
    rnd1 = (rg1 << (GL_Q - 12)) - _IQ(1.0);                 // rnd (-1.0, 1.0)
    rnd2 = (rg2 << (GL_Q - 12));                            // rnd ( 0.0, 1.0)
    return _IQmpy( _IQcos(_IQmpy(_IQ(M_PI), rnd1), hlp),
                   _IQsqrt(-_IQlog(rnd2) << 1) );
}

Листинг 2. Фрагмент файла IQMath.h (язык Си)

Для оценки времени исполнения функций была составлена программа для демонстрационной платы STM32F3DISCOVERY с 32-х разрядным микроконтроллером и тактовой частотой генератора 72 МГц. В программе, по прерыванию таймера, вызывалась функция IQgaussian. Генерируемые псевдослучайные числа выводились в ЦАП, см. осциллограмму на рис. 1. Код транслировался с оптимизацией и без. В листинге 2, в комментариях, указано время исполнения соответствующих функций (IQgaussian – 4.68 мкс).

Осциллограмма шума Гаусса сгенерированная функцией IQgaussian
Рис. 1. Шум Гаусса, сгенерированный функцией IQgaussian и выведенный
в ЦАП демонстрационной платы STM32F3DISCOVERY

Литература

  1. Регистр сдвига с линейной обратной связью // Википедия – свободная энциклопедия. (дата обращения: 23.07.2018).
  2. Philip Koopman. Maximal Length LFSR Feedback Terms // Philip Koopman's personal website. – URL: http://users.ece.cmu.edu/~koopman/lfsr/index.html (дата обращения: 23.07.2018).
  3. Преобразование Бокса – Мюллера // Википедия – свободная энциклопедия. (дата обращения: 23.07.2018).
  4. Клиначёв Н.В. Векторное управление в автоматизированных комплексах: Методические указания к лабораторным работам // Моделирующая программа Jigrein: Теория, программа, руководство, модели. – 2006-2018 гг. – URL: http://model.exponenta.ru/k2/Jigrein/foc_cntnt.htm (дата обращения: 23.07.2018).
  5. Набатчиков А.М. "Дедовский" способ генерации стандартного шума из равномерного // Личный сайт Александра Набатчикова. – URL: http://nabatchikov.com/blog/view/randn (дата обращения: 23.07.2018).

20.07.2018