Алгоритм долговременной генерации розового шума во временной области

Тема данной статьи была навеяна публикацией Автокорреляционная функция фликкер-шума / Хабр, в которой даны теоретические оценки автокорреляционной функции фликкер-шума. Однако, практическое сравнение было дано лишь для броуновского шума, потому что он тривиально получается из белого шума интегрированием. Проверка розового шума осталась под вопросом, потому что для его генерации во временной области требуется «расщепить» интегратор на два одинаковых каскада, что, как будет показано далее, нетривиально, а практическая реализация может быть лишь приближенной.

Аналоговые способы (фактически, в частотной области) генерации розового шума оставляем за скобкой, потому что это отдельная тема: там используется свойство, что мощность шума одинакова в каждой октаве, то есть диапазоны частот равномощны. Также популярны алгоритмы, основанные на быстром преобразовании Фурье (это тоже в частотной области). Однако, они требуют специальных вычислений (хотя сейчас это не проблема) и выдают результат в виде вектора отсчетов — это их главная особенность. Также Фурье алгоритмы требуют специальных ухищрений: наложения оконной функции, сдвига максимума спектра, обработки полюса на нулевой частоте.

Сфокусируемся на цифровом способе долговременной генерации розового шума во временной области, то есть отсчет за отсчетом. Обзор показал, что такой алгоритм, как ни странно, почти отсутствует, что и побудило к созданию такового. Исключение составляет рекомендация [1], в которой без объяснений предлагается делить каждый отсчет белого шума на корень квадратный из номера отсчета. Однако, это не подходит для долговременной генерации шума, потому что для каждого отсчета шума придется вычислять новое значение корня, что затратно и рано или поздно приведет к исчерпанию разрядности, и, вообще, выглядит подозрительно, хотя небольшие реализации таким способом получать можно. Попытка пояснить находится в спойлере.

Скрытый текст

Пусть белый шум масштабируется корнем квадратным от времени, тогда его преобразование Фурье будет выглядеть так

$P(f) = 2 \int_{0}^{\infty} {n(t) \over \sqrt{t}} {e}^{-2 \pi j f t} dt.$

Внесем этот корень под знак дифференциала

$P(f) = 4 \int_{0}^{\infty} n(t) {e}^{-2 \pi j f t} d \sqrt{t}.$

Сделаем замену переменной ft = x

$P(f) = {4 \over {\sqrt{f}}} \int_{0}^{\infty} n\left({x \over f}\right) {e}^{-2 \pi j x} d \sqrt{x}.$

Так как шум — белый, то он инвариантен к масштабированию времени (и его нельзя нарисовать во временной области): квадрат модуля интеграла в среднем не будет зависеть от коэффициента масштабирования , то есть от частоты. Заметим, что подынтегральная экспонента не зависит от частоты, в отличие от формулы до замены переменной. Значит

$\overline{{\left|P(f) \right|}^{2}} = {{C} \over {f}},$

где — константа. Это подтверждает то, что масштабирование отсчетов белого шума корнем квадратным от номера отсчета дает розовый шум, как бы это ни казалось странным… Но белость исходного шума здесь очень важна.

В качестве исходной посылки будем использовать то, что z-образ интегратора равен

${K}_{\int}(z) = {{1} \over {1 - {z}^{-1}}}.$

Фактически, эта функция является обратной к передаточной функции дифференциатора, который считает «дельты» некоторого процесса. Если разложить в ряд Тейлора передаточную функцию интегратора, то коэффициентами при степенях ${z}^{-1}$ будут единицы. Откликом интегратора на единицу (цифровую дельта-функцию) будет ступенчатая функция — интегратор имеет бесконечную память. Если на интегратор подавать случайные числа, то откликом будет сумма всех ступенек, уровни которых равны этим числам — интегратор является линейным звеном.

Если сделать стандартную подстановку

$z = {e}^{~j 2 \pi f T_d},$

то можно построить частотную характеристику (ЧХ) звена, амплитудную и фазовую, где T_d — интервал дискретизации. ЧХ будет периодической, потому что рассматривается дискретный по времени процесс. Для малых частот ЧХ будет пропорциональной 1/f , что согласуется с ЧХ непрерывного интегратора. Для частот возле половинной частоты дискретизации ЧХ будет пропорциональна малой по сравнению с уровнем плотности в окрестности нулевой частоты константе. Эта константа является результатом наложения идеальных ЧХ, периодически продолженных кратно частоте дискретизации. Таким образом, высокочастотные компоненты, подаваемые на вход дискретного интегратора, будут несколько искаженно интегрироваться. Это можно увидеть явно, если подать знакочередующуюся последовательность единиц

которая соответствует синусоидальному сигналу с частотой, равной половине частоты дискретизации. Откликом интегратора (с нулевым начальным состоянием) будет следующая последовательность

Фактически, чередование как было, так и осталось: основная гармоника осталась равной половине частоты дискретизации, фаза сигнала не поменялась, сместился лишь уровень.

Чтобы увидеть контраст, когда интегрирование действительно работает, можно подать на вход прямоугольный импульс, например, из единиц

Откликом интегратора будет следующая последовательность

Как говорится, разница в уровнях налицо. Здесь уровень подрос до против при интегрировании знакочередующегося процесса. Чем более низкочастотен процесс, тем больше у него условных единичек и прямо пропорционально выше отклик интегратора — это и отражает амплитудно-частотная характеристика 1/f интегратора.

Каскадному соединению двух звеньев соответствует перемножение передаточных функций, поэтому расщепление интегратора в частотной области тривиально

$K(z) = {{1} \over {\sqrt{1 - {z}^{-1}}}}.$

Разложение в ряд Тейлора данной иррациональности известно

$K(z) = \sum_{n=0}^{\infty}{{\binom{-1/2}{n}}{(-1)}^{n}{z}^{-n}},$

где

$\binom{-1/2}{n} \approx {{1} \over {\sqrt{\pi n}}} + O\left({n}^{-3/2}\right), ~~~~n \gg 1$

биномиальный коэффициент, который для больших убывает обратно пропорционально корню квадратному, то есть довольно-таки медленно. Коэффициенты при степенях ${z}^{-1}$ являются импульсной характеристикой (ИХ) фильтра, генерирующего в данном случае розовый шум, если на его вход подать белый шум

$K(z) = \sum_{n \geq 0}^{} {h_n {z}^{-n}} = 1 + {{1}\over{2}} {z}^{-1} + {{3}\over{8}} {z}^{-2} + {{5}\over{16}} {z}^{-3} + ....$

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

${{{h}_{n}} \over {{h}_{n-1}}} = 1 - {{1} \over {2n}} \neq \text{const},$

что не позволяет найти рекуррентное уравнение с постоянными коэффициентами относительно h_n , то есть построить БИХ-фильтр некоторого небольшого порядка, который бы точно отрабатывал исходную ИХ генератора розового шума. Это прямо следует уже из того, что передаточная функция содержит радикал, то есть в знаменателе фактически находится полином бесконечной степени, что говорит о бесконечных порядках как БИХ-фильтра, так и КИХ. Значит, точная реализация предполагает бесконечную память, а на практике остается надеяться только на разные аппроксимации.

Опишем трюк, который позволяет ускорить затухание исходной ИХ.

Откликом линейного фильтра с постоянными коэффициентами h_n на входную последовательность x_n является свертка

$y_n = \sum_{j \geq 0}^{} {{h}_{n-j} {x}_{j}}.$

Здесь предполагается, что процессы развиваются начиная с нулевого индекса, то есть верно

Можно показать, что следующий отсчет отклика равен

${y}_{n+1} = y_n + h_0 {x}_{n+1} - 1/2 \sum_{j \geq 0}^{} {{{h}_{n-j} x_j}\over{n+1-j}}.$

Не считая слагаемого $h_0 x_{n+1}$ , следующий отсчет отклика содержит предыдущий отсчет за вычетом масштабированного коэффициентом 1/2 отклика фильтра с модифицированной ИХ

$\widetilde{h_n} = {{h_n} \over {n+1}}, ~~~~n \geq 0,$ ${y}_{n+1} = y_n + h_0 {x}_{n+1} - {1\over2} \widetilde{y_n}.$

Коэффициенты модифицированного фильтра связаны уравнением

${{\widetilde{h_n}}\over{\widetilde{{h}_{n-1}}}} = {1 - {{3}\over{2(n+1)}}}.$

Данные коэффициенты убывают быстрее с ростом индекса

$\widetilde{h_n} = 2{(-1)}^{n}\binom{1/2}{n+1} = {{1}\over{\sqrt{\pi}}} \left[ { {\left({1}\over{n}\right)}^{3/2} - {{9}\over{8}} {\left({1}\over{n}\right)}^{5/2} + {{145}\over{128}} {\left({1}\over{n}\right)}^{7/2} - ... } \right].$

Назовем такую характеристику как «ИХ 3/2». Для аппроксимации возьмем три первых слагаемых. Однако, остается задача аппроксимировать степенные зависимости , 5/2 и 7/2 некоторым БИХ-фильтром.

К счастью, есть такая вещь как аппроксимация степенных зависимостей суммой кратных экспонент (навеяно методом Прони из [2], хотя вытащить оттуда нижеприведенную формулу нетривиально, но спасибо моему научному руководителю, который когда-то помог понять суть и написать практичную формулу [3])

${\left({1}\over{n}\right)^p} \approx {V{\sum_{k=-R}^{R} {\tau}^{kp} {e}^{-{\tau}^{k}n}}}, ~~~~ n > 0.$

Здесь — степень аппроксимируемого степенного закона, 2R+1 — количество экспонент, $\tau$ — параметр, типовые значения которого зависят от и который подбирается так, чтобы минимизировать ошибку аппроксимации для заданного диапазона по индексу, — нормировочный множитель, который вычисляется при начальном индексе n=1 .

Физика данной аппроксимации основана на суммировании экспонент с разной постоянной времени и разными амплитудами так, что более растянутые экспоненты имеют малые амплитуды, менее растянутые — большие, и таким образом покрывается некоторый интервал степенной зависимости: чем больше экспонент, тем больше покрытие.

Каждая экспонента может быть тривиально реализована БИХ-фильтром первого порядка (как геометрическая прогрессия), поэтому, в итоге, общий порядок фильтра будет равен 2R+1 .

Так как мы взяли три слагаемых «ИХ 3/2», то после такой аппроксимации остается некоторая ошибка, наибольшее значение которой сосредоточено возле малых индексов. Ее можно учесть КИХ-фильтром небольшого порядка: в данном случае оптимален восьмой порядок (подбирается графически).

Для реализации генератора розового шума был выбран -й порядок БИХ-части для степени и порядки и для степеней 5/2 и 7/2 , соответственно. Это дает хорошее соответствие итоговой ИХ до индексов $n \approx 500000$ . При этом оптимальный вектор параметра «тау» равен $\tau = (3.933, 1.005, 1)$ .

Таким образом, найдены выражения, позволяющие реализовать генератор розового шума с помощью линейного фильтра с постоянными параметрами небольшого порядка (-го порядка в данном примере).

Предложенная схема на очень больших индексах — если построить ИХ данного фильтра — выдает некоторый уровень постоянной составляющей, своего рода DC offset, который является артефактом аппроксимации: вспомогательный фильтр $\widetilde{y_n}$ при аппроксимации его экспонентами имеет быстро затухающую до нуля ИХ, что не доводит до конца интегрирование ${y}_{n+1} = y_n + ...$ , и итоговый фильтр как бы замирает на некотором пусть малом, но ненулевом уровне. Это не позволяет использовать данный генератор на большие времена — уровень DC offset при подаче на вход белого шума постепенно выйдет за пределы разрядности вычислителя, поэтому необходима некоторая коррекция.

Отрезок импульсной характеристики генератора розового шума: синяя кривая - эталон, красная - результат аппроксимации фильтром 41-го порядка: видна постоянная составляющая минус 67 дБ. — Отрезок импульсной характеристики генератора розового шума: синяя кривая — эталон, красная — результат аппроксимации фильтром 41-го порядка: видна постоянная составляющая минус 67 дБ.

Если бы мы продолжили выше описанный трюк с повышением скорости затухания ИХ, то могли бы аппроксимировать суммами экспонент, начиная со степени 5/2 , что привлекает, однако, артефактом аппроксимации при этом был бы не DC offset, а линейный тренд (за счет двойного интегрирования), и т. д. Поэтому остановимся на аппроксимации с итоговым артефактом типа DC offset и попытаемся его уменьшить.

Коррекция DC offset основана на введении БИХ-фильтра второго порядка дополнительно к выходу вспомогательного фильтра $\widetilde{y_n}$ . Фильтр-корректор как бы активируется к заданному моменту времени и компенсирует недостающую сумму. Передаточная функция корректора DC offset предлагается следующей

${K}_{\text{DC}}(z) = {{\epsilon} \over {1 - {\alpha}_{1} z^{-1}}} - {{\epsilon} \over {1 - {\alpha}_{2} z^{-1}}},$

где $\epsilon$ — некоторый уровень (подбирается графически), $\alpha_i, ~i={1,2}$ — два параметра, которые показывают интервал времени, когда отклик фильтра наибольший. Были выбраны следующие значения

$\alpha_1 = 1 - {10}^{-6},$ $\alpha_2 = 1 - 2 \cdot {10}^{-6},$ $\epsilon = 1,685 \cdot 10^{-9}.$

Индекс, при котором достигается половина амплитуды $\epsilon$ , равен

$n_i \approx {{2} \over {|1 - \alpha_i |}}.$

Данная аппроксимация справедлива для параметра «альфа» близкого к единице.

ИХ фильтра-корректора нарастает от до примерно $\epsilon$ , и затем убывает до нуля — такая структура фильтра была выбрана намеренно, чтобы не портить аппроксимацию для малых индексов и скорректировать аппроксимацию для больших индексов ИХ. Остаточный уровень DC offset при этом около минус дБ.

Отрезок импульсной характеристики генератора розового шума: синяя кривая - эталон, красная - результат аппроксимации фильтром 43-го порядка: виден результат коррекции постоянной составляющей с минус 67 дБ до минус 132 дБ. — Отрезок импульсной характеристики генератора розового шума: синяя кривая — эталон, красная — результат аппроксимации фильтром 43-го порядка: виден результат коррекции постоянной составляющей с минус 67 дБ до минус 132 дБ.

В данной работе для работы генератора розового шума использовался формат float 32 бит, а для расчета точной ИХ использовался double 64 бит.

Неустранимый остаток DC offset предлагается корректировать динамически: нелинейным образом, накапливая текущий DC offset (как среднее арифметическое) и вычитая его один раз через каждые отсчетов, после чего сбрасывая как счетчик, так и сам DC offset для последующего накопления, и т. д. Выбрано $N = {2}^{22}$ , то есть около $4 \cdot {10}^{6}$ , что соответствует половине интервала работы корректора DC offset. Конечно, это примерные прикидки, чтобы убедиться в качестве работы генератора розового шума на очень больших временах: отсутствии смещений, которые могли бы привести к перегрузке (насыщению). Моделирование показало работоспособность предлагаемой коррекции.

Реализация розового шума: 8192 отсчета. Соответствует подаче на вход разработанного генератора отсчетов белого шума со стандартным нормальным распределением N(0, 1) амплитуд.

Типовой DC offset при коррекции составляет не более-х единиц, чаще менее, при этом в процессе генерации розового шума интервал «гуляния» составляет (-10, +10) , что является нормальным низкочастотным поведением шума, ведь плотность мощности сильно возрастает в окрестности нулевой частоты.

Оценка спектра плотности мощности генерируемого розового шума методом БПФ с усреднением квадратов амплитуд экспоненциальным сглаживанием.

Попробовать генератор в действии и убедиться в его долговременной работе можно здесь nawww83/pink_noise_generator: Генератор розового шума

Источники

Как сгенерировать шум в python и применить его в своих проектах?
Марпл С. Л. Цифровой спектральный анализ и его приложения, Мир, 1990.
Новиков А. В. Эффективная реализация оператора дискретного прозрачного граничного условия для двумерного параболического уравнения, 2011 International Siberian Conference on Control and Communications (SIBCON), https://conf.sfu-kras.ru/uploads/Novikov-AV-sibconIEEE.doc

ссылка на оригинал статьи https://habr.com/ru/articles/895628/

Алгоритм долговременной генерации розового шума во временной области

Комментарии

Добавить комментарий Отменить ответ