8-900-374-94-44
[email protected]
Slide Image
Меню

Сглаживание сигнала: Сглаживание цифровых сигналов / Хабр

Содержание

Сглаживание цифровых сигналов / Хабр

Введение

Данную статью меня заставил написать пост habrahabr.ru/post/183986, где не совсем правильно используется некоторый алгоритм сглаживания изображения.

Сразу перейдём к сути дела.

Математические модели цифровых сигналов — вектора и матрицы, элементами которых являются числа. Числа могут быть двоичными (бинарный сигнал), десятичными («обычный» сигнал) и так далее. Любой звук, любое изображение и видео могут быть преобразованы в цифровой сигнал1: звук — в вектор, изображение — в матрицу, а видео — в последовательный набор матриц. Поэтому цифровой сигнал — это, можно сказать, универсальный объект для представления информации.

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

Задача сглаживания может использоваться при прореживании сигналов, то есть когда, например, необходимо отобразить большую картинку на небольшой экран. Или когда частота дискретизации звука снижается, например, с 48000 Гц до 44100 Гц. Понижение частоты выборок — коварная операция, требующая предварительной обработки сигнала (низкочастотной фильтрации), но это — тема отдельного разговора…

Приведём пример «плохого» сглаживания

Казалось бы, обычное усреднение и сигнал на выходе должен быть «гладким». Но как определить, насколько он стал «глаже»? Не переборщили ли мы? А может быть некоторые коэффициенты выбрать не по 1/3? А может быть усреднить по пяти точкам? Как определить насколько ослабляются частотные составляющие в сигнале? Как найти свой (то есть для конкретной задачи) оптимум?
На эти и некоторые другие вопросы я постараюсь ответить так, чтобы «обычный» программист смог обосновать свой алгоритм, — надеюсь, не только алгоритм на тему «Сглаживание», так как идеи будут излагаться весьма общие, заставляющие думать самому


Сигналы

В этом параграфе под сигналом понимается вектор, то есть массив из определённого количества чисел. Например, вектор из четырёх элементов s = (2,5; 5; 0; -5).

Для простоты будут рассматриваться только десятичные вещественные числа.

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

Как уже было сказано, один из самых простых способов сглаживания, это
(1)
где s — исходный сигнал, v — сглаженный сигнал.

Способ (1) основан на сглаживающем свойстве суммирования, ведь каждому ясно, что средняя величина, вычисляемая как сумма многих случайных чисел, с ростом количества суммируемых чисел становится всё менее и менее похожей на случайную величину 2, которая, попросту говоря, и есть шум 3.

Но на чём основан выбор коэффициентов в уравнении (1)? На том, что так вычисляется среднее? Вроде бы да, но… А если взять не три слагаемых, а шестнадцать? А тридцать два?.. Почему всё более отстоящие от центрального элемента s[i] отсчёты должны браться с одинаковым весом? Ведь может оказаться так, что связь между отсчётами будет постепенно ослабевать с ростом расстояния

4 между ними?

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

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

Как оценить насколько сильно связаны соседние отсчёты? Вычислить автокорреляционную функцию (АКФ). Желающим можно предложить провести эксперимент по записи разных слов, фраз, повторов фраз и последующему построению АКФ (благо, например, программа Matlab позволяет это сделать, особо не задумываясь над кодом и формулами).

Так как всё-таки выбрать коэффициенты фильтра в (1)?

В данном случае удобно рассмотреть реакцию фильтра на единичное воздействие, то есть на сигнал вида

Например, фильтр (1) даст следующий отклик (импульсную характеристику)

откуда мы можем сделать вывод, что после сглаживания длительность сигнала стала равна трём элементам. А если взять фильтр из пяти элементов?.. Правильно, длительность выходного сигнала будет равна пяти элементам. Насколько это полезно, определяется конкретной ситуацией (задачей).

Кстати, а долгожданный артефакт уже налицо! Импульсная характеристика (1) — это, по сути, прямоугольный импульс, нисколько не являющийся гладким!.. Странно, да? А если взять пятиточечный фильтр? Тогда на выходе получим более длинный прямоугольный импульс, но с меньшей амплитудой. Не очень хорошо выходит… Простейший тест говорит о непригодности простого усреднения для сглаживания.

Посмотрим на фильтр (1) с частотной стороны (с временной мы уже посмотрели).

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

амплитудно-частотная характеристика фильтра в рабочем диапазоне частот должна быть достаточно равномерной.

Пропустим через рассматриваемый фильтр один однотональный сигнал определённой частоты, естественно, не выходя за пределы теоремы отсчётов. Пусть шаг дискретизации по времени Td равен единице, то есть отсчёты идут с шагом одна секунда. Возьмём сигнал с частотой f = 1/(3Td) = 1/3 Гц, то есть

Ограничимся двумя периодами

Отклик фильтра (1) будет равен

Как ни странно, получили почти ноль… Получается, можно потерять некоторые составляющие полезного сигнала.

Проверим отклик на несколько более высокочастотный сигнал

Как видим, форму сигнала не потеряли. В чём же дело?..

Дело в том, что амплитудно-частотная характеристика фильтра (1) не является монотонной в рабочей полосе частот (в полосе от нуля до половины частоты дискретизации) и имеет один нуль на частоте, в три раза меньшей частоты дискретизации. Как это показать?

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

Обозначим спектр сигнала s[i] на входе как S( f ), тогда спектр задержанного на один такт Td сигнала s[i-1] выразится через спектр исходного сигнала как
(2)
где j — мнимая единица.

Спектр опережённого сигнала s[i+1] выразится как
(3)

Что означает мнимая единица? И как можно обосновать (2) и (3)?
Если записать известные [1, 2] ряды для синуса, косинуса и экспоненты
(4)
то подобрав число j, такое, что j2 = –1, можно последний ряд выразить через два первых
(5)
что означает то, что любое комплексное число можно записать через экспоненту с мнимым показателем. Модуль числа (5) равен единице (корень квадратный из суммы квадратов мнимой и реальной частей), поэтому для записи любого комплексного числа в форме (5) его необходимо разделить и умножить на свой модуль
(6)

Из (5) и (6) следует, что если в показателе экспоненты можно выделить мнимую единицу, умноженную на некоторое действительное число, то это число является аргументом комплексного числа.

В данном случае рассматриваются сигналы, поэтому модулю комплексного числа соответствует амплитуда гармонического сигнала, а аргументу — фаза. Если, например, взять сигнал вида
(7)
то можно выделить амплитуду A и фазу Ф. Множитель — это также фаза, и в некоторых случаях её выносят за скобки. Например, при прохождении сигнала (7) через некоторый фильтр важно знать разность фаз сигналов на входе и выходе, которую вносит фильтр на заданной частоте f.

Если сигнал (7) задержать на величину t0, то получится тот же самый сигнал, но сдвинутый по фазе на постоянную величину
(8)
то есть при задержке произвольного сигнала все его частотные составляющие сдвигаются по фазе на величину, зависящую от текущей частоты и величины задержки. Этим можно объяснить формулы (2) и (3). Поэтому при анализе какого-либо алгоритма важна и фазо-частотная характеристика, которая показывает на какое время задерживает фильтр (алгоритм) каждую частотную составляющую входного сигнала. Низкие частоты и высокие в общем случае будут иметь разную задержку в фильтре.

Из (1), (2) и (3) следует, что частотная характеристика фильтра (передаточная функция) будет иметь вид
(9)

Так как спектр выходного сигнала линейно выражается через спектр входного, то при выводе (9) спектр входного сигнала успешно сокращается. Далее замечаем, что частотная характеристика фильтра (1) — действительная, то есть данный фильтр не вносит фазовых искажений 6. Этого мы достигли (скорее всего, неосознанно) за счёт симметричности формулы (1): каждый отсчёт на выходе фильтра равен сумме текущего и двух соседних отсчётов.

Физически такой алгоритм реализуем только при наличии запоминающих устройств, так как в нём используются опережающие отсчёты (для вычисления отсчёта s[i] требуется отсчёт s[i+1]). В настоящее время это не является большой проблемой и, как правило, используют симметричные алгоритмы. Если фазовые искажения окажутся полезными, то — пожалуйста, главное осознанно применять какой-либо алгоритм и смотреть на него с разных сторон: с частотной и временной.

Построить график зависимости (9) от частоты не составляет труда. Для упрощения вводят нормированную частоту

f0  = f Td , полезный диапазон изменения которой [0…0,5]. Составляющие сигнала с частотами выше половины частоты дискретизации по теореме отсчётов должны отсутствовать (перед оцифровкой сигнал пропускают через соответствующие фильтры нижних частот). Частота дискретизации показывает количество выдаваемых цифровым устройством отсчётов в секунду. Если, например, один такт Td равен одной миллисекунде, то за одну секунду должна быть выдана тысяча отсчётов.

Анализируя (9) можно заметить, что на некотором промежутке коэффициент передачи меньше нуля, а амплитуда — это по определению положительная величина… Выход из ситуации — построить модуль передаточной функции, убрав знак минус в фазовую характеристику, которая, всё-таки, не является константой (нулём). Если взять число «минус единица», то его по формуле (5) можно представить как
(10)
то есть комплексным числом, имеющим модуль «единица» и фазу «180 градусов» («пи»).

Таким образом, трёхточечный симметричный алгоритм (1) для некоторых «высоких» частот вносит сдвиг по фазе на 180 градусов, то есть попросту инвертирует входной сигнал. Этот эффект можно заметить, анализируя рассмотренный выше отклик на нормированную частоту 2/5 Гц.


Рис. 1. Амплитудно-частотная и фазо-частотная характеристики трёхточечного симметричного алгоритма сглаживания (1)

Из рис. 1 сразу следует, что сигнал с частотой 1/3 будет данным алгоритмом подавлен, а сигналы, имеющие частоту выше 1/3 будут инвертированы. Таким образом, полезный рабочий диапазон частот можно смело сократить с [0…0,5] до [0…1/3]. Если нас не устраивает быстрое спадание коэффициента передачи, придётся определять другой алгоритм, имеющий более прямоугольную амплитудно-частотную характеристику и при этом — ещё неизвестно какую фазовую…

По сути, полученная немонотонная частотная характеристика — следствие прямоугольной формы импульсной характеристики …, 0, 0, …, 1/3, 1/3, 1/3, 0, 0, …. Поэтому импульсная характеристика — простой способ заметить грубые изъяны в алгоритмах. Частотная характеристика, несмотря на некоторую сложность её вычисления, удобна тем, что позволяет заметить и устранить более тонкие изъяны, чем мы и займёмся.

Если теперь записать алгоритм (1) в более общем виде
(11)
то, основываясь на знании частотной характеристики, можно попытаться подобрать коэффициенты так,
чтобы амплитудно-частотная характеристика стала монотонной в рабочем диапазоне частот (0…0,5). Для этого, как минимум, необходимо отсутствие нулей внутри рабочего диапазона.

Так как у нас нет оснований выделять запаздывающий отсчёт s[i–1] по отношению к опережающему s[i+1], то приравняем коэффициенты a1 и a3. После запишем коэффициент передачи
(12)

Попытаемся переместить первый нуль на частоту f0 = 0,5. Для этого должно выполняться равенство a2 = 2a1 , то есть вес у боковых отсчётов должен быть в два раза меньше веса центрального. Тогда более оптимальный алгоритм будет выглядеть так
(13)

Как в (13) выбрать единственный коэффициент a1?

Взглянем на алгоритм (13) с точки зрения импульсной переходной характеристики. Для этого найдём отклик на единичный скачок :
(14)

Как видим, чтобы сохранить выходную амплитуду в установившемся режиме, равную единице, требуется выбрать коэффициент a1=1/4. То есть сумма всех коэффициентов должна быть равна единице.

Наконец, для готового алгоритма
(16)
можно построить (рис. 2) частотные характеристики: амплитудную и фазовую.


Рис. 2. Амплитудно-частотная и фазо-частотная характеристики трёхточечного симметричного алгоритма сглаживания (16)

Анализ рис. 2 показывает, что фазовые искажения исчезли, а амплитудная характеристика стала монотонной в рабочей полосе частот [0…0,5]. В каком-то смысле мы выжали из трёхточечного фильтра «всё».

Теперь становится очевидно, что простое усреднение далеко не всегда является оптимальным, особенно когда усредняемых отсчётов много.

Изображения

Что касается изображений, то здесь в какой-то мере всё то же самое, что и для сигналов: амплитудные и фазовые искажения, импульсные характеристики, энергия сигнала, — основное отличие заключается в том, что используются матрицы вместо векторов. Также следует учесть, что для изображений нет координаты времени, а есть координата пространства, поэтому есть пространственные частоты дискретизации. Данное отличие скорее формальное, так как на математику алгоритмов оно не влияет вообще.

Рассмотрим теперь «самый простой» алгоритм сглаживания изображения по соседним точкам (рис. 3). Отсчёт v00 на выходе фильтра
(17)


Рис. 3. Схема сглаживания изображения по соседним точкам

В формуле (17) специально выделены три слагаемых A, B и C, так как четыре соответствующих внутренних слагаемых у B и C имеют свои расстояния от центрального отсчёта s00. Здесь естественно предположить, что максимальный вес будет у слагаемого A, затем в порядке убывания — у B и C.

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

Фильтр (17) имеет следующий коэффициент передачи
(18)

Если взять, например, около тридцати дискретных частот и построить (рис. 4) контурный график модуля частотной характеристики (18), то можно увидеть искажения, аналогичные искажениям на рис. 1. Левому нижнему углу соответствуют самые низкие частоты. Наблюдается провал частотной характеристики на частотах, составляющих 1/3 от частоты дискретизации (частоты дискретизации в данном случае равны 64 по обеим координатам).


Рис. 4. Амплитудно-частотная характеристика девятиточечного симметричного алгоритма сглаживания (17)

Перетасовывая коэффициенты для отсчётов на рис. 3, можно получить следующий алгоритм сглаживания
(19)

Причём если единичное воздействие поместить в центре изображения, то отклик фильтра (19) будет иметь вид
(20)

По сути, (20) является импульсной характеристикой фильтра (это девять главных отсчётов, так как остальные равны нулю). Сумма всех отсчётов импульсной характеристики равна единице. Частотная характеристика фильтра (19) имеет вид, показанный на рис. 5.


Рис. 5. Амплитудно-частотная характеристика девятиточечного симметричного алгоритма сглаживания (19)

Опять мы видим (сравни рис. 4 и рис. 5) улучшение алгоритма обычного усреднения.

Резюме

При разработке алгоритмов обработки цифровых сигналов (сглаживания и не только) не следует доверять интуитивным алгоритмам вроде простого усреднения.

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

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

Большую роль играют амплитудно-частотные характеристики фильтров, которые определяются коэффициентами весового суммирования. Они позволяют доказать, что заданный фильтр пропускает определённый диапазон частот и заграждает другой. Причём можно определить неравномерность в полосе пропускания, в полосе заграждения и так далее. Чем больше порядок фильтра, тем больше степеней свободы и тем лучше можно подобрать форму амплитудно-частотной характеристики.

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

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

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

Сноски

1. Всегда с потерями, так как частота работы и разрядность цифровых устройств ограничены
2. Стабилизирующее свойство средней величины справедливо при неизменности характеристик генератора случайных чисел, то есть в идеале генератор должен выдавать случайную величину с заданным законом распределения
3. Если, конечно, наблюдаемый шум — это не зашифрованный полезный сигнал, который становится случайным для тех, кто не имеет ключа…
4. Для звукового сигнала расстояние между отсчётами — это некоторый промежуток времени
5. Можно считать, что любой сигнал можно представить в виде суммы синусоид кратных частот, но звуковой по природе обязан «хорошо» раскладываться в ряд по частотам
6. Упрощённо можно сказать и так, подробности смотри ниже по тексту
Источники

1. Тригонометрические функции [Электронный ресурс], режим доступа: свободный, ru.wikipedia.org/wiki/%D2%F0%E8%E3%EE%ED%EE%EC%E5%F2%F0%E8%F7%E5%F1%EA%E8%E5_%F4%F3%ED%EA%F6%E8%E8
2. Экспонента [Электронный ресурс], режим доступа: свободный, ru.wikipedia.org/wiki/%DD%EA%F1%EF%EE%ED%E5%ED%F2%E0

зашумляем сигнал, чтобы улучшить его / Хабр


Введение


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

Блокнот Mathematica для воспроизведения результатов можно найти здесь, а pdf-версия находится здесь.

Что такое дизеринг?


Дизеринг (Dithering) можно описать как намеренное/осознанное внесение в сигнал шума для предотвращения ошибок большого масштаба/низкого разрешения, возникающих вследствие дискретизации или субдискретизации.

Если вы когда-нибудь работали с:

  • Аудиосигналами,
  • Палитровыми форматами изображений 90-х

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

Однако я обнаружил в Википедии довольно удивительный факт о том, как впервые был определён и использован дизеринг:

…Одна из первых областей применения дизеринга возникла во время Второй мировой войны. Самолёты-бомбардировщики для навигации и вычислений траектории бомб использовали механические компьютеры. Любопытно, что эти компьютеры (ящики, заполненные сотнями шестерёнок и зубчатых передач) работали точнее при полёте на борту самолёта, чем на земле. Инженеры поняли, что вызываемая самолётом вибрация снижает ошибку, вызываемую липкими подвижными деталями. На земле они двигались короткими отрывистыми движениями, а в воздухе их движение было более непрерывным. В компьютеры встроили небольшие вибромоторчики, и их вибрацию назвали «дизером» (dither), от среднеанглийского слова «didderen», означающего «дрожать». Сегодня, когда вы стучите по механическому измерителю, чтобы повысить его точность, вы применяете dither. При использовании в небольших количествах dither успешно превращает систему оцифровки в более аналоговую, в хорошем смысле этого слова.

— Кен Полманн, Principles of Digital Audio


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

Но хватит истории, давайте для начала рассмотрим процесс дизеринга в 1D-сигналах, например, в аудио.

Дискретизация дизерингом постоянного сигнала


Мы начнём с анализа самого скучного в мире сигнала — постоянного сигнала. Если вы знаете немного о цифровой обработке сигналов, связанных со звуком, то можете сказать: но ты же обещал рассмотреть аудио, а в звуке по определению нет постоянной составляющей! (Более того, и в ПО, и в оборудовании обработки звука намеренно устраняется так называемый сдвиг постоянной составляющей (DC offset).)

Это правда, и вскоре мы рассмотрим более сложные функции, но начнём мы сначала.

Представьте, что мы выполняем 1-битную дискретизацию нормализованного сигнала с плавающей запятой. Это значит, что мы имеем дело только с конечными двоичными значениями, 0 или 1.

Если сигнал равен 0,3, то простое округление без дизеринга будет самой скучной функцией — просто нулём!

Погрешность тоже постоянна и равна 0,3, следовательно, и средняя погрешность равна 0,3. Это означает, что мы внесли довольно большое отклонение в сигнал и полностью потеряли информацию исходного сигнала.

Мы можем попробовать выполнить дизеринг этого сигнала и посмотреть на результаты.

Дизеринг в этом случае (при использовании функции округления) просто применяет обычный случайный белый шум (случайное значение для каждого элемента, что создаёт равномерный спектр шума) и прибавляет в сигнал перед дискретизацией случайное в интервале (-0.5, 0.5).

quantizedDitheredSignal =
Round[constantSignalValue + RandomReal[] – 0.5] & /@ Range[sampleCount];



Здесь сложно что-то увидеть, теперь результат дискретизации — это просто набор случайных единиц и нулей. С (ожидаемо) большим количеством нулей. Сам по себе этот сигнал не особо интересен, однако довольно интересен график погрешностей и средняя погрешность.
Итак, как мы и ожидали, погрешность тоже варьируется, но пугает то, что погрешность иногда стала больше (абсолютное значение 0,7)! То есть максимальная погрешность к сожалению стала хуже, однако средний шум имеет значение:
Mean[ditheredSignalError]
0.013

Намного лучше, чем первоначальная погрешность в 0,3. При значительно большом количестве сэмплов эта погрешность будет стремиться к нулю (к пределу). Итак, погрешность постоянной составляющей стала намного меньше, но давайте взглянем на частотный график всех погрешностей.
Красный график/всплеск = частотный спектр погрешности при отсутствии дизеринга (постоянный сигнал без частот). Чёрный — с дизерингом при помощи белого шума.

Всё становится интереснее! Это демонстрирует первый вывод из этого поста — дизеринг распределяет погрешность/отклонение дискретизации среди множества частот.

В следующем разделе мы узнаем, как это нам поможет.

Частотная чувствительность и низкочастотная фильтрация


Выше мы наблюдали за дизерингом дискретизированного постоянного сигнала:
  • Он увеличил максимальную погрешность.
  • Почти обнулил среднюю погрешность.
  • Добавил к спектру частот погрешностей постоянный белый шум (с полным покрытием спектра), снизив низкочастотную погрешность.

Сам по себе он нам не очень полезен… Однако мы рассматриваем не дискретизацию любой произвольной математической функции/сигнала. Мы рассматриваем сигналы, которые будут восприниматься человеком. Очевидно, что восприятие человека ограничено. Вот несколько примеров этого:
  • Наше зрение имеет предел остроты. У многих людей есть близорукость и без очков они видят размытые изображения далёких объектов.
  • Мы воспринимаем средний масштаб деталей гораздо лучше, чем очень высокие или очень низкие частоты (мелкие детали очень плавных градиентов могут быть незаметными).
  • Наш слух работает в определённом диапазоне частот (20 Гц — 20 кГц, но со временем ухудшается) и наиболее чувствительны мы к среднему диапазону — 2 кГц — 5 кГц.

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

Более того, наши медиаустройства становятся всё лучше и лучше, обеспечивая большую избыточную дискретизацию (oversampling). Например, в случае телевизоров и мониторов у нас есть технология «retina» и 4K-дисплеи (на которых невозможно разглядеть отдельный пиксель), в области звука мы используем форматы файлов с дискретизацией не менее 44 кГц даже для дешёвых динамиков, которые часто не могут воспроизводить больше, чем 5-10 кГц.

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


Красный — желаемый недискретизированный сигнал. Зелёный — дискретизированный сигнал с дизерингом. Синий — низкочастотный фильтр этого сигнала.

Сигнал начинает выглядеть гораздо более близким к исходной, недискретизированной функции!

К сожалению, мы начинаем видеть низкие частоты, которые очень заметны и отсутствуют в исходном сигнале. В третьей части серии мы попробуем исправить при помощи синего шума. А пока вот как график может выглядеть с функцией псевдо-шума, имеющей содержимое с гораздо меньшей частотой:

Это возможно, потому что наша псевдослучайная последовательность имеет следующий спектр частот:
Но давайте закончим рассматривать простые, постоянные функции. Взглянем на синусоиду (если вы знакомы с теоремой Фурье, то знаете, что она является строительным блоком любого периодического сигнала!).

Дискретизация синусоиды


Если мы дискретизируем синусоиду 1-битной дискретизацией, то получим простой прямоугольный сигнал.
Прямоугольный сигнал довольно интересен, потому что включает в себя и базовую частоту, и нечётные гармоники.

Это интересное свойство, которое активно используется в аналоговых субтрактивных синтезаторах для создания звучания полых/медных инструментов. Субтрактивный синтез берёт сложный, гармонически богатый звук и фильтрует его, устраняя некоторые частоты (параметры фильтра варьируются со временем), чтобы придать звукам нужную форму.

Спектр частот прямоугольного сигнала:


Но в этом посте нас больше интересую погрешности дискретизации! Давайте создадим график погрешности, а также спектр частот погрешности:
В этом случае ситуация гораздо лучше — средняя погрешность близка к нулю! К сожалению, у нас по-прежнему присутствует множество нежелательных низких частот, очень близких к нашей основной частоте (нечётных множителей с уменьшающейся величиной). Это явление называется алиасингом или шумом дизеринта — возникают частоты, отсутствовавшие в исходном сигнале, и они имеют довольно большие величины.

Даже низкочастотная фильтрация не сможет значительно помочь сигналу. Погрешность имеет очень много низких частот:


Дискретизированная синусоида с низкочастотной фильтрацией
Погрешность дискретизированной синусоиды с низкочастотной фильтрацией

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


Однако если мы рассмотрим это как изображение, то оно начинает выглядеть лучше:
Заметьте, что погрешности дискретизации снова распределены среди различных частот:
Выглядит очень многообещающе! Особенно учитывая то, что теперь мы можем попробовать выполнить фильтрацию:
Это немного искажённая синусоида, но она выглядит намного ближе к исходной, чем версия без дизеринга, за исключением фазового сдвига, внесённого асимметричным фильтром (я не буду объяснять этого здесь; скажу только, что проблему можно устранить, применив симметричные фильтры):
Красный — исходная синусоида. Зелёный — подвергнутый низкочастотной фильтрации сигнал без дизеринга. Синий — подвергнутый низкочастотной фильтрации сигнал с дизерингом.

Графики обеих погрешностей численно подтверждают, что погрешность намного меньше:


Красный — погрешность подвергнутого низкочастотной фильтрации сигнала без дизеринга. Синий — погрешность подвергнутого низкочастотной фильтрации сигнала с дизерингом.

Наконец, давайте вкратце рассмотрим сигнал с более качественной функцией дизеринга, содержащей только высокие частоты:


Верхнее изображение — функция белого шума. Нижнее изображение — функция, содержащая более высокие частоты.
Версия с низкочастотной фильтрацией, дизерингом и улучшенной функцией — почти идеальные результаты, если не учитывать фазовый сдвиг, вызванный фильтром!

И наконец, сравнение всех трёх спектров погрешностей:


Красный — спектр погрешности дискретизации без дизеринга. Чёрный — спектр погрешности дискретизации с дизерингом белым шумом. Синий — спектр погрешности дискретизации с дизерингом с более высокими частотами.

Итог


На этом первая часть серии заканчивается. Основные выводы:
  • Дизеринг распределяет погрешность/отклонение дискретизации среди множества частот, и это зависит от функции дизеринга; они не сосредотачиваются в области низких частот.
  • Восприятие человеком любого сигнала (звука, зрения) лучше всего работает в очень конкретных диапазонах частот. Сигналы часто избыточно дискретизируются на границах спектра восприятия, где восприятие малочувствительно. Например, стандартные частоты дискретизации звука позволяют воспроизводить сигналы, которые вообще не может услышать большинство взрослых. Из-за предыдущего пункта это делает очень привлекательной идею использования дизеринга и сдвига погрешностей в этот диапазон частот.
  • Различные функции шума создают разные спектры погрешностей, которые можно использовать, зная, какой спектр погрешностей наиболее желателен.

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

Курс лекций «Основы цифровой обработки сигналов» / Хабр

Всем привет!

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

Большая часть обучающего материала для наглядного и интерактивного представления реализована с использованием Jupyter Notebook. Предполагается, что читатель имеет базовые знания из области высшей математики, а также немного владеет языком программирования Python.



Список лекций


Этот курс содержит материалы в виде законченных лекций по разным тематикам из области цифровой обработки сигналов. Материалы представлены с использованием библиотек на языке Python (пакеты numpy, scipy, matplotlib, и т.д.). Основная информация для этого курса взята из моих лекций, которые я, будучи аспирантом, читал студентам Московского Энергетического Института (НИУ МЭИ). Частично информация из этих лекций была использована на обучающих семинарах в Центре Современной Электроники, где я выступал в качестве лектора. Кроме того, в этот материал входит перевод различных научных статей, компиляция информации из достоверных источников и литературы по тематике цифровой обработки сигналов, а также официальная документация по прикладным пакетам и встроенным функциям библиотек scipy и numpy языка Python.

Для пользователей MATLAB (GNU Octave) освоение материала с точки зрения программного кода не составит труда, поскольку основные функции и их атрибуты во многом идентичны и схожи с методами из Python-библиотек.

Все материалы сгруппированы по основным тематикам цифровой обработки сигналов:

  1. Сигналы: аналоговые, дискретные, цифровые. Z-преобразование,
  2. Преобразование Фурье: амплитудный и фазовый сигнала, ДПФ и БПФ,
  3. Свертка и корреляция. Линейная и циклическая свертка. Быстрая свёртка,
  4. Случайные процессы. Белый шум. Функция плотности вероятностей,
  5. Детерминированные сигналы. Модуляция: АМ, ЧМ, ФМ, ЛЧМ. Манипуляция,
  6. Фильтрация сигналов: БИХ, КИХ фильтры,
  7. Оконные функции в задачах фильтрации. Детектирование слабых сигналов,
  8. Ресемплинг: децимация и интерполяция. CIC-фильтры, фильтры скользящего среднего,
  9. Непараметрические методы спектрального анализа,
  10. Усреднение по частоте и по времени. Полифазный БПФ.

Список лекций — достаточный но, разумеется, неполный для вводного знакомства с областью ЦОС. При наличии свободного времени я планирую поддерживать и развивать этот проект.

Где найти?


Все материалы — абсолютно бесплатны и доступны в виде открытого репозитория на моем гитхабе как opensource проект. Материалы представлены в двух форматах — в виде тетрадок Jupyter Notebook для интерактивной работы, изучения и редактирования, и в виде скомпилированных из этих тетрадок HTML-файлов (после скачивания с гитхаба имеют вполне пригодный формат для чтения и для печати).

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

Сигналы. Z-преобразование


Вводный раздел, в котором содержится основная информация по типам сигналов. Вводится понятие дискретной последовательности, дельта-функции и функции Хевисайда (единичный скачок).

Все сигналы по способу представления на множестве можно разделить на четыре группы:

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

Для правильного восстановления аналогового сигнала из цифрового без искажений и потерь используется теорема отсчетов, известная как Теорема Котельникова (Найквиста-Шеннона).

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

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

Также в этом разделе описывается Z-преобразование и его свойства, показывается представление дискретных последовательностей в Z-форме.

Пример конечной дискретной последовательности:

x(nT) = {2, 1, -2, 0, 2, 3, 1, 0}
.
Пример этой же последовательности в Z-форме:

X(z) = 2 + z-1 — 2z-2 + 2z-4 + 3z-5 + 1z-6

Преобразование Фурье. Свойства. ДПФ и БПФ


В этом разделе описывается понятие временной и частотной области сигнала. Вводится определение дискретного преобразования Фурье (ДПФ). Рассмотрены прямое и обратное ДПФ, их основные свойства. Показан переход от ДПФ к алгоритму быстрого преобразования Фурье (БПФ) по основанию 2 (алгоритмы децимации по частоте и по времени). Отражена эффективность БПФ в сравнении с ДПФ.

В частности, в этом разделе описывается Python пакет scipy.ffpack для вычисления различных преобразований Фурье (синусное, косинусное, прямое, обратное, многомерное, вещественное).

Преобразование Фурье позволяет представить любую функцию в виде набора гармонических сигналов! Преобразование Фурье лежит в основе методов свертки и проектировании цифровых корреляторов, активно применяется при спектральном анализе, используется при работе с длинными числами.

Особенности спектров дискретных сигналов:
1. Спектральная плотность дискретного сигнала – периодическая функция с периодом, равным частоте дискретизации.
2. Если дискретная последовательность вещественная, то модуль спектральной плотности такой последовательности есть четная функция, а аргумент – нечетная функция частоты.

Спектр гармонического сигнала:

Сравнение эффективности ДПФ и БПФ

Эффективность алгоритма БПФ и количество выполняемых операций линейно зависит от длины последовательности N:
Как видно, чем больше длина преобразования, тем больше экономия вычислительных ресурсов (по скорости обработки или количеству аппаратных блоков)!

Любой сигнал произвольной формы можно представить в виде набора гармонических сигналов разных частот. Иными словами, сигнал сложной формы во временной области имеет набор комплексных отсчетов в частотной области, которые называются *гармоники*. Эти отсчеты выражают амплитуду и фазу гармонического воздействия на определенной частоте. Чем больше набор гармоник в частотной области, тем точнее представляется сигнал сложной формы.

Свертка и корреляция


В этом разделе вводится понятие корреляции и свертки для дискретных случайных и детерминированных последовательностей. Показана связь автокорреляционной и взаимнокорреляционной функций со сверткой. Описываются свойства свертки, в частности, рассмотрены методы линейной и циклической свертки дискретного сигнала с подробным разбором на примере дискретной последовательности. Кроме того, показан метод вычисления «быстрой» свертки с помощью алгоритмов БПФ.

В реальных задачах часто ставится вопрос о степени похожести одного процесса на другой или же о независимости одного процесса от другого. Иными словами, требуется определить взаимосвязь между сигналами, то есть найти корреляцию. Методы корреляции используются в широком диапазоне задач: поиск сигналов, компьютерное зрение и обработка изображений, в задачах радиолокации для определения характеристик целей и определения расстояния до объекта. Кроме того, с помощью корреляции производится поиск слабых сигналов в шумах.

Свертка описывает взаимодействие сигналов между собой. Если один из сигналов — импульсная характеристика фильтра, то свертка входной последовательности с импульсной характеристикой есть ни что иное, как реакция цепи на входное воздействие. Иными словами, результирующий сигнал отражает прохождение сигнала через фильтр.

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

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

Как видно, для длин БПФ до 64, быстрая свёртка проигрывает у прямого метода. Однако, при увеличении длины БПФ результаты меняются в обратную сторону — быстрая свертка начинает выигрывать у прямого метода. Очевидно, чем больше длина БПФ, тем лучше выигрыш частотного метода.


Случайные сигналы и шум


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

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

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

В задачах ЦОС случайные сигналы делятся на два класса:

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

С помощью случайных величин можно моделировать воздействие реальной среды на прохождение сигнала от источника к приёмнику данных. При прохождении сигнала через какое-то шумящее звено, к сигналу добавляется так называемый белый шум. Как правило, спектральная плотность такого шума равномерно (одинаково) распределена на всех частотах, а значения шума во временной области распределены нормально (Гауссовский закон распределения). Поскольку белый шум физически добавляется к амплитудам сигнала в выбранные отсчеты времени, он называется аддитивный белый гауссовский шум (AWGN — Additive white Gaussian noise).

Сигналы, модуляция и манипуляция


В этом разделе показаны основные способы изменения одного или нескольких параметров гармонического сигнала. Вводятся понятия амплитудной, частотной и фазовой модуляции. В частности, выделяется линейная частотная модуляция, применяемая в задачах радиолокации. Показаны основные характеристики сигналов, спектры модулированных сигналов в зависимости от параметров модуляции.

Для удобства на языке Python создан набор функций, осуществляющих перечисленные виды модуляции. Пример реализации ЛЧМ-сигнала:

def signal_chirp(amp=1.0, freq=0.0, beta=0.25, period=100, **kwargs):
    """
    Create Chirp signal

    Parameters
    ----------
    amp : float
        Signal magnitude
    beta : float
        Modulation bandwidth: beta < N for complex, beta < 0.5N for real
    freq : float or int
        Linear frequency of signal
    period : integer
        Number of points for signal (same as period)
    kwargs : bool
        Complex signal if is_complex = True
        Modulated by half-sine wave if is_modsine = True
    """
    is_complex = kwargs.get('is_complex', False)
    is_modsine = kwargs.get('is_modsine', False)

    t = np.linspace(0, 1, period)
    tt = np.pi * (freq * t + beta * t ** 2)
    
    if is_complex is True:
        res = amp * (np.cos(tt) + 1j * np.sin(tt))
    else:
        res = amp * np.cos(tt)

    if is_modsine is True:
        return res * np.sin(np.pi * t)
    return res

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

Цифровые фильтры — БИХ и КИХ


Достаточно большой раздел, посвященный вопросам цифровой фильтрации дискретных последовательностей. В задачах цифровой обработки сигналов данные проходят через цепи, которые называются фильтрами. Цифровые фильтры, как и аналоговые, обладают различными характеристиками — частотные: АЧХ, ФЧХ, временная: импульсная характеристика, а также передаточная характеристика фильтра. Цифровые фильтры используются в основном для улучшения качества сигнала — для выделения сигнала из последовательности данных, либо для ухудшения нежелательных сигналов — для подавления определенных сигналов в приходящих последовательностях отсчетов.

В разделе перечислены основные преимущества и недостатки цифровых фильтров (в сравнении с аналоговыми). Вводится понятие импульсной и передаточной характеристик фильтра. Рассматривается два класса фильтров — с бесконечной импульсной характеристикой (БИХ) и конечной импульсной характеристикой (КИХ). Показан способ проектирования фильтров по канонической и прямой форме. Для КИХ фильтров рассматривается вопрос о способе перехода к рекурсивной форме.

Для КИХ фильтров показан процесс проектирования фильтра от стадии разработки технического задания (с указанием основных параметров), до программной и аппаратной реализации — поиска коэффициентов фильтра (с учетом формы представления числа, разрядности данных и т.д.). Вводятся определения симметричных КИХ фильтров, линейной ФЧХ и её связи с понятием групповой задержки.

Оконные функции в задачах фильтрации


В задачах цифровой обработки сигналов используются оконные функции различной формы, которые при наложении на сигнал во временной области, позволяют качественно улучшить его спектральные характеристики. Большое количество всевозможных окон обусловлено в первую очередь одной из главных особенностей любого оконного наложения. Эта особенность выражается во взаимосвязи уровня боковых лепестков и ширины центрального лепестка. Правило:
Чем сильнее подавление боковых лепестков спектра, тем шире главный лепесток спектра и наоборот.

Одно из применений оконных функций: обнаружение слабых сигналов на фоне более сильных путём подавления уровня боковых лепестков. Основные оконные функции в задачах ЦОС — **треугольное, синусоидальное, окно Ланцоша, Ханна, Хэмминга, Блэкмана, Харриса, Блэкмана-Харриса, окно с плоской вершиной, окно Наталла, Гаусса, Кайзера** и множество других. Большая часть из них выражена через конечный ряд путём суммирования гармонических сигналов с определенными весовыми коэффициентами. Такие сигналы отлично реализуются на практике на любых аппаратных устройствах (программируемые логические схемы или сигнальные процессоры).

Ресемплинг. Децимация и интерполяция


В этом разделе рассматриваются вопросы многоскоростной обработки сигналов — изменения частоты дискретизации. Многоскоростная обработка сигналов (multirate processing) предполагает, что в процессе линейного преобразования цифровых сигналов возможно изменение частоты дискретизации в сторону уменьшения или увеличения, либо в дробное число раз. Это приводит к более эффективной обработке сигналов, так как открывается возможность использования минимально допустимых частот дискретизации и, как следствие, значительного уменьшения требуемой вычислительной производительности проектируемой цифровой системы.

Децимация (прореживание) – понижение частоты дискретизации. Интерполяция – повышение частоты дискретизации.

Также в разделе рассматривается класс однородных КИХ фильтров, которые называются интегрально-гребенчатыми фильтрами (CIC, Cascaded integrator–comb). Показана реализация, основные свойства и особенности CIC фильтров. В силу линейности математических операций, происходящих в CIC фильтре возможно каскадное соединение нескольких фильтров подряд, что дает пропорциональное уменьшение уровня боковых лепестков, но также увеличивает «завал» главного лепестка амплитудно-частотной характеристики.

График АЧХ фильтра в зависимости от коэффициента децимации:

Также в этом разделе обсуждается вопрос увеличения разрядности данных на выходе CIC фильтра в зависимости от его параметров. Это особенно важно в задачах программной реализации, в частности на ПЛИС.

Для практической реализации CIC фильтров на Python разработан отдельный класс CicFilter, реализующий методы децимации и интерполяции. Также показаны примеры изменения частоты дискретизации с помощью встроенных методов из scipy пакета Python.

Python CicFilter Class for Digital Signal Processing
class CicFilter:
    """
    Cascaded Integrator-Comb (CIC) filter is an optimized class of
    finite impulse response (FIR) filter.
    CIC filter combines an interpolator or decimator, so it has some
    parameters:

    R - decimation or interpolation ratio,
    N - number of stages in filter (or filter order)
    M - number of samples per stage (1 or 2)*

    * for this realisation of CIC filter just leave M = 1.

    CIC filter is used in multi-rate processing. In hardware
    applications CIC filter doesn't need multipliers, just only
    adders / subtractors and delay lines.

    Equation for 1st order CIC filter:
    y[n] = x[n] - x[n-RM] + y[n-1].


    Parameters
    ----------
    x : np.array
        input signal
    """

    def __init__(self, x):
        self.x = x

    def decimator(self, r, n):
        """
        CIC decimator: Integrator + Decimator + Comb

        Parameters
        ----------
        r : int
            decimation rate
        n : int
            filter order
        """

        # integrator
        y = self.x[:]
        for i in range(n):
            y = np.cumsum(y)

        # decimator

        y = y[::r]
        # comb stage
        return np.diff(y, n=n, prepend=np.zeros(n))

    def interpolator(self, r, n, mode=False):
        """
        CIC inteprolator: Comb + Decimator + Integrator

        Parameters
        ----------
        r : int
            interpolation rate
        n : int
            filter order
        mode : bool
            False - zero padding, True - value padding.
        """

        # comb stage
        y = np.diff(self.x, n=n,
                    prepend=np.zeros(n), append=np.zeros(n))

        # interpolation
        if mode:
            y = np.repeat(y, r)
        else:
            y = np.array([i if j == 0 else 0 for i in y for j in range(r)])

        # integrator
        for i in range(n):
            y = np.cumsum(y)

        if mode:
            return y[1:1 - n * r]
        else:
            return y[r - 1:-n * r + r - 1]

Наконец, в этом разделе приведен особый класс фильтров — скользящего среднего. Показано три способа реализации: через свертку сигналов, с помощью КИХ-фильтра и БИХ-фильтра.

Заключение


Надеюсь, этот курс лекций в совокупности с моими предыдущими статьями по цифровой обработке сигналов на ПЛИС принесет практическую пользу и поможет читателю лучше понять основы цифровой обработки сигналов. Этот проект будет улучшаться и дополняться новым полезным и не менее интересным материалом. Следите за развитием!

Дополнительно к этому материалу я поддерживаю и развиваю свой проект по основным модулям ЦОС (на языке Python). Он содержит пакет генерации различных сигналов, класс CIC фильтров для задач децимации и интерполяции, алгоритм расчета коэффициентов корректирующего КИХ-фильтра, фильтр скользящего среднего, алгоритм вычисления сверх-длинного БПФ через методы двумерного преобразования (последнее очень пригодилось в работе при аппаратной реализации на ПЛИС).

UPD: 20.04.2020


В курс добавлено две лекции:
  1. Непараметрические методы спектрального анализа (Владимир Фадеев)
  2. Усреднение по частоте и по времени. Полифазный БПФ.

TODO:
  1. Вейвлет анализ
  2. STFT, мел-спектрограммы, преобразование Гриффина-Лима

Спасибо за внимание!

Скользящее среднее / Блог компании Нерепетитор.ру / Хабр

Принято считать, что две базовые операции «машинного обучения» — это регрессия и классификация. Регрессия — это не только инструмент для выявления параметров зависимости y(x) между рядами данных x и y (чему я уже посвятил несколько статей), но и частный случай техники их сглаживания. В этом примере мы пойдем чуть дальше и рассмотрим, как можно проводить сглаживание, когда вид зависимости y(x) заранее неизвестен, а также, как можно отфильтровать данные, которые контролируются разными эффектами с существенно разными временными характеристиками.

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


Приведенный график демонстрирует биржевое значение курса доллара к рублю с интервалом в 1 час. Исходные данные представлены синей кривой, а сглаженные — красной. Даже невооруженным глазом видно, что колебания курса имеют несколько характерных частот, что и является предметом одного из направлений технического анализа рынков.

Сглаживание при помощи «скользящего среднего»
Принцип сглаживания на основе «скользящего среднего» (МА — от англ. «moving average») состоит в расчете для каждого значения аргумента yi среднего значения по соседним w данным.

Число w называют окном скользящего усреднения: чем оно больше, тем больше данных участвуют в расчете среднего, соответственно, тем более гладкая кривая получается. На верхнем рисунке окно w=50, а вот как будет выглядеть скользящее усреднение при w=200.

Как меняется Фурье-спектр данных при скользящем усреднении?
Очевидно, что при малых w сглаженные кривые практически повторяют ход изменения данных, а при больших w — отражают лишь закономерность их медленных вариаций. Это типичный пример фильтрации данных, т.е. устранения одной из составляющих зависимости y(xi). Наиболее часто целью фильтрации является подавление быстрых вариаций y(xi), которые обычно обусловлены шумом. В результате из быстроосциллирующей зависимости y(xi) получается другая, сглаженная зависимость, в которой доминирует более низкочастотная составляющая.

Эти рассуждения плавно перевели нас к терминологии спектров. Давайте нарисуем график Фурье-преобразования («Фурье-спектр») исходных данных:

и убедимся в том, что спектр скользящего среднего вырезает из него высокие частоты (начиная примерно с частоты 0.005 Гц):

Пояснение: Фурье-спектр суммы синусов и ее МА
Для того чтобы пояснить принцип расчета Фурье-спектра, рассмотрим вместо (случайных) исходных данных простую модель суммы нескольких детерминированных сигналов (синусоид с разной частотой и амплитудой) и псевдослучайного шума:

Приведем графики этой суммы и ее МА (с тем же окном w=200):

а также их Фурье-спектры:

Из них видно, что скользящее усреднение вырезает из сигнала высокие частоты, начиная с частоты 0.005 Гц. Лучше это видно на крупном плане низкочастотной области спектра:

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

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

Также интерес представляют смешанные задачи выделения среднемасштабных вариаций путем подавления как более быстрых, так и более медленных вариаций. Одна из возможностей решения связана с применением «полосовой фильтрации», которая реализуется так:
1. Устранение из сигнала y высокочастотной составляющей, имеющее целью получить сглаженный сигнал middle, например, с помощью скользящего усреднения с малым окном (например, w=200).
2. Выделение из сигнала middle низкочастотной составляющей slow, например, путем скользящего усреднения с большим окном w.
3. Вычитание из сигнала middle тренд slow, тем самым выделяя среднемасштабную составляющую исходного сигнала y.

Оставляю заинтересовавшемуся читателю возможность самому реализовать полосовую фильтрацию в Mathcad Express (с небольшой оговоркой, что в бесплатной версии Маткада алгоритм БПФ расчета Фурье-спектров отключен, и с ним поиграться не получится). Сами расчеты находятся здесь.

Литература:
1. Кирьянов Д.В., Кирьянова Е.Н. Вычислительная физика (PDF, гл.1, п.6 и 7). М.: Полибук Мультимедиа, 2006.
2. Бат М. Спектральный анализ в геофизике. М., Наука, 1980.

Об одном забавном подходе к фильтрации унимодальных сигналов

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

Задача оффлайновой фильтрации сигналов в случае, когда ожидаемая форма сигнала известна с точностью до нескольких неизвестных параметров, сводится к задаче аппроксимации. Например, если известно, что сигнал линейно растет на рассматриваемом промежутке, задача сведётся к линейной регрессии, а если можно предположить, что шум — нормален, то правильным методом будет МНК. Но однажды мы столкнулись с задачей оценки формы профиля рентгеновского микрозонда (пучка), про которую априори было достоверно известно только одно: профиль унимодален, а именно имеет ровно один максимум. Оказывается, и в этом случае можно наилучшим (в смысле, например, L2 метрики) образом приблизить экспериментальный сигнал функцией, принадлежащей известному множеству (множеству унимодальных функций). Причём — с приемлемой ассимптотикой вычислительной сложности.

===> ===>

Итак, речь пойдёт о задаче аппроксимации некоторого зашумлённого сигнала как функции дискретного аргумента , причем известно, что исходный сигнал — унимодальный. В рассматриваемом случае дополнительно известно, что на краях сигнал ассимптотически приближается к нулю. Следует отметить, что рассматриваемый ниже метод сам по себе не требует такого допущения.

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

Если предположить, что шум нормально распределен, то задача сведётся к минимизации квадратичного функционала в присутствии неравенств вида . Это — задача квадратичного программирования, и для нее на просторах интернета можно найти некоторое количество готовых решателей. Но, во-первых, задача квадратичного программирования в общем случае обсчитывается довольно медленно, а во-вторых, задача аппроксимации унимодального сигнала не выписывается в подобном виде, если положение экстремума не известно заранее.

Существует альтернативное решение, не обладающее указанными недостатками, но при его использовании шкалу сигнала придётся считать дискретной. Тут стоит отметить, что задача не интересна в случае слабого шума, поскольку обычная линейная фильтрация, скорее всего, снимет все практические проблемы. Что же такое интересный, сильный шум? На рис. 1 можно увидеть модельный сигнал до и после такого зашумления. Шум здесь задан нормальным распределением . В подобных условиях дискретизация аппроксимирующего сигнала не выглядит существенной потерей. Кроме того, как будет видно далее, при увеличении числа отсчетов сложность вычислений будет возрастать линейно.

Рис. 1: a) Модельный монотонный сигнал до налолжения шума. b) Модельный монотонный сигнал после наложения шума.

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

О том, что такое динамическое программирование, можно почитать здесь, а как и к каким задачам его применяют — здесь.

Фильтрация монотонного сигнала методом динамического программирования


В рассматриваемой задаче таблица динамического программирования — двумерная, её столбцам i соответствуют возрастающие слева направо значения дискретного аргумента экспериментальных точек, а строкам j — возрастающие снизу вверх возможные значения аппроксимирующей функции . Минимизируемый функционал — сумма квадратов отклонений от вдоль путей с допустимыми переходами. Очевидно, что допустимыми будут переходы на 1 вправо с одновременным неотрицательным смещением вверх.

Внешний цикл ведётся по координате i слева направо. Для каждого значения выполняется внутренний цикл по координате j снизу вверх. Минимальное значение функционала определяется по формуле


причём за пределами таблицы функционал считается равным нулю.

При этом, чтобы не пересчитывать каждый раз минимум на меньших значениях j, его полезно сохранять в процессе заполнения вместе со значением функционала, а поскольку в конце требуется восстановить путь, сохранять удобнее значение m, на котором этот минимум был достигнут:


Тогда вычисление самого функционала можно выполнять за O(1) для каждой ячейки:
После завершения прямого прохода остаётся получить аппроксимацию в явном виде. Для этого ищется минимум среди значений функционала для крайнего правого столбца (столбца с наибольшей координатой i), а затем выполняется обратный проход по i, с восстановлением индекса значения отфильтрованного сигнала j* в каждом столбце:
Окончательные значения определяются как
На рис. 2 показан пример результата описанного алгоритма для монотонного сигнала.

Рис. 2: Фильтрация монотонно возрастающего сигнала динамическим программированием. Синим цветом показан входной зашумлённый сигнал, зелёным — результат работы алгоритма, красным — идеальный незашумлённый сигнал.

На этом рассмотрение вспомогательного монотонного случая закончено.

Унимодальный сигнал


Используя полученный результат, можно решить задачу об унимодальном сигнале, заполнив 2 таблицы ДП: одна заполняется как для монотонной функции, вторая — таким же образом, но в обратную сторону. В результате вычисляется 2 функционала: и , а также значения и для восстановления путей по таблицам. После этого среди всех значений i и j ищется такоая пара , для которой сумма функционалов минимальна — это будут индексы максимума аппроксимированной функции. Значения функции с каждой стороны от него восстанавливаются по соответствующей таблице.


На рисунках 3 и 4 показан пример исходного унимодального сигнала и результата его обработки предложенным алгоритмом.

Рис. 3: Зашумлённый унимодальный сигнал.

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

Ограничение на производную


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

Правда, чтобы не проиграть в трудоёмкости, нужно иметь возможность получать минимальное значение функционала на отрезке [j — dj, j] за константное время.
Это можно сделать с помощью алгоритма Ван Херка — Гила — Вермана, описанного, например, здесь.

Результат применения алгоритма с ограничением на производную показан на рис. 5.

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

Разоблачение чёрной магии


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

Чтобы рассеять ощущение возможного мошенничества, на рис. 6 и 7 приведены примеры обработки сигнала гауссовским сглаживанием и медианной фильтрацией с зашкаливанием в нуле. Параметры фильтров подбирались с целью получения оптимального результата, но без фанатизма.

Рис. 6: Фильтрация унимодального сигнала гауссовским сглаживанием. Синим цветом показан входной зашумлённый сигнал, зелёным — результат фильтрации, красным — идеальный незашумлённый сигнал.

Рис. 7: Медианная фильтрация унимодального сигнала. Синим цветом показан входной зашумлённый сигнал, зелёным — результат фильтрации, красным — идеальный незашумлённый сигнал.

Как видно, на достаточно произвольных данных метод динамического программирования дает результат получше, чем простейшие общеизвестные фильтры. При этом сравнивать их на модельных примерах подробно — абсолютно бессмысленно. Изложенный выше метод ценен ровно тем, что является примером “костюма, сшитого по мерке”: формализованная постановка задачи была положена в саму основу метода. В рассматриваемом случае это было условие унимодальности. Аналогично, линейные методы “по мерке” сигналам, частотно разделимым с шумом. Еще один инструмент в наборе — что может быть лучше для инженера?

Алгоритмы антиалиасинга в реальном времени / Хабр


Алиасинг (aliasing) — это, возможно, наиболее фундаментальный и самый широко обсуждаемый артефакт 3D-рендеринга всех времён. Однако в игровом сообществе его часто недопонимают. В этой статье я подробно расскажу о теме сглаживания (антиалиасинга, anti-aliasing, AA) в реальном времени, особенно о том, что касается игр, и в то же время буду излагать всё достаточно простым языком.

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

Эту программу можно скачать здесь.

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


«Если ты знаешь себя и знаешь врага, то не подвергнешься опасности и в сотне битв»

Как учит нас Сунь Цзы, чтобы победить врага, нам нужно сначала понять его. Врагом — простите меня за излишнюю драматичность — методов сглаживания являются артефакты алиасинга. Поэтому нам первым делом нужно понять, как и откуда появляется алиасинг.

Термин алиасинг был впервые введён в области обработки сигналов, в которой он изначально описывал эффект, возникающий, когда разные непрерывные сигналы становятся неразличимыми (или начинают искажать друг друга) при дискретизации. В 3D-рендеринге этот термин обычно имеет более конкретное значение: он относится ко множеству нежелательных артефактов, которые могут возникать, когда 3D-сцена рендерится для отображения на экране, состоящем из фиксированной сетки пикселей.

В этом случае 3D-сцена является непрерывным сигналом, а процесс генерирования значений цветов для каждого пикселя дискретизирует этот сигнал для создания выходных данных рендеринга. Цель методов антиалиасинга заключается в том, чтобы выходные данные как можно точнее походили на сцену на заданной сетке пикселей, при этом минимизируя визуально искажающие артефакты.

На Рисунке 1 показан алиасинг в простой сцене, состоящей из единственного белого треугольника на чёрном фоне. На этапе растеризации стандартного рендеринга сэмплируется центральная позиция каждого пикселя: если он находится в треугольнике, то пиксель будет закрашен белым, в противном случае он закрашивается чёрным. В результате получается хорошо заметный эффект «лесенки», один из самых узнаваемых артефактов алиасинга.

При идеальном сглаживании для каждого пикселя определяется, какая часть его площади закрыта треугольником. Если пиксель закрыт на 50%, то он должен быть заполнен цветом на 50% между белым и чёрным (средним серым). Если он закрыт меньше, то должен быть пропорционально темнее, если больше — то пропорционально светлее. Полностью закрытый пиксель является белым, полностью незакрытый — чёрным. Результат этого процесса показан на четвёртом рисунке. Однако выполнение этого вычисления в реальном времени в общем случае является невыполнимой задачей.

Рисунок 1. Простейший алиасинг.


1-1. Сетка 8x8 с помеченными центрами
1-2. Сетка 8x8 с треугольником
1-3. Сетка 8x8 с растеризированным треугольником
1-4. Сетка 8x8 с идеально сглаженными выходными даннымиТо же самое в виде GIF

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

Поэтому чтобы в полной мере обсудить относительные сильные и слабые стороны техник сглаживания, мы сгруппировали артефакты алиасинга, возникающие при 3D-рендеринге, в пять отдельных категорий. Это группирование зависит от точных условий генерирования артефактов. На Рисунке 2 показаны эти типы алиасинга на реальном примере, отрендеренном с помощью OpenGL.


Рисунок 2: Различные типы алиасинга. Слева направо, сверху вниз:
• Единственный выровненный относительно экрана прямоугольник с частично прозрачной текстурой.
• «Мельница», состоящая из выровненных относительно экрана переменных белых и чёрных треугольников.
• Несколько чёрных линий различной ширины, начиная с 1 пикселя сверху до 0,4 пикселя снизу, и белая линия толщиной 0,5, отображающая синусоиду.
• Куб, состоящий из шести плоских закрашенных прямоугольников
• Наклонная плоскость, текстурированная высокочастотной текстурой травы.
• Выровненный относительно экрана прямоугольник с пиксельным шейдером, определяющим цвет каждого пикселя на основе функции синуса.

Самым распространённым типом алиасинга, о котором мы уже говорили, является геометрический алиасинг. Этот артефакт возникает, когда какой-то примитив сцены (обычно треугольник) частично пересекается с пикселем, но это частичное перекрытие не учитывается в процессе рендеринга.

Алиасинг прозрачности возникает в текстурированных примитивах с частичной прозрачностью. Верхнее левое изображение на Рисунке 2 отрендерено с использованием одного прямоугольника, заполненного частично прозрачной текстурой сетчатого забора. Поскольку сама текстура — это просто фиксированная сетка пикселей, её нужно сэмплировать в точках, на которые накладывается каждый пиксель отрендеренного изображения, и для каждой такой точки должно приниматься решение, нужна ли в нём прозрачность. В результате возникает та же проблема сэмплирования, которую мы уже встречали на сплошной геометрии.

Несмотря на то, что фактически он является типом геометрического алиасинга, подпиксельный алиасинг требует особого рассмотрения, так как он ставит уникальные задачи для аналитических методов сглаживания, которые недавно получили большую популярность в рендеринге игр. Мы подробно рассмотрим их в статье. Подпиксельный алиасинг возникает тогда, когда растеризируемая структура накладывается менее чем на один пиксель в сетке буфера кадров. Такое чаще всего происходит в случае узких объектов — шпилей, телефонных или электрических линий, или даже мечей, когда они находятся достаточно далеко от камеры.

Рисунок 3. Иллюстрация подпиксельного алиасинга.


3-1. Сетка 8x8 с отмеченными центрами
3-2. Сетка 8x8 с двумя отрезками прямых
3-3. Сетка 8x8 с растеризированными отрезками, без AA
3-4. Сетка 8x8 с идеально сглаженным треугольникомТо же самое в виде GIF
На Рисунке 3 показан подпиксельный алиасинг в простой сцене, состоящей из двух отрезков прямых. Верхний имеет ширину в один пиксель, и хотя при растеризации он демонстрирует знакомый артефакт-«лесенку» геометрического алиасинга, результат всё равно в целом соответствует по форме входным данным. Нижний отрезок имеет ширину полпикселя. При растеризации часть пересекаемых им столбцов пикселей не имеет одного центра пикселя в пределах отрезка. В результате он разделяется на несколько несвязанных фрагментов. То же самое можно заметить на прямых линиях и кривой синусоиды на Рисунке 2.

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

Рисунок 4: Анимированная высокочастотная текстура с артефактами мерцания
Текстурный алиасинг обычно можно предотвратить использованием mip-текстурирования и фильтрацией высококачественных текстур, но он всё равно иногда остаётся проблемой, особенно с некоторыми версиями драйверов популярных видеопроцессоров, субдискретизирующих высокоанизотропные текстуры. На него также влияют различные методы антиалиасинга, поэтому он тоже включён в демонстрационную программу.

И, наконец, шейдерный алиасинг возникает, когда программа пиксельного шейдера, выполняемая для каждого пикселя и определяющая его цвет, генерирует результат с алиасингом. Такое часто случается в играх с шейдерами, создающими контрастное освещение, например, зеркальные засветы на основании карты нормалей, или с техниками контрастного освещения типа cel shading или задней подсветкой. В демонстрационной программе это аппроксимируется одним шейдером, вычисляющим функцию синуса для координат текстуры и закрашивающего все отрицательные результаты чёрным, а все положительные — белым.


Вооружившись пониманием артефактов алиасинга и всех типов алиасинга, которые могут возникнуть при рендеринге 3D-сцены, мы можем начать исследование техник антиалиасинга. Эти техники можно разбить на две категории: техники, пытающиеся снизить алиасинг увеличением количества генерируемых при рендеринге сэмплов и техники, пытающиеся смягчить артефакты алиасинга анализом и постобработкой сгенерированных изображений. Категория техник сглаживания на основе сэмплирования в более проста, поэтому стоит начать с неё.

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

Интуитивным решением будет простое увеличение количества сэмплов, взятых на пиксель. Эта концепция показана на Рисунке 5.


Рисунок 5: треугольник, растеризированный с четырьмя упорядоченными сэмплами на пиксель

Центры пикселей снова помечены красными точками. Однако в каждом пикселе сэмплируется на самом деле четыре отдельных места (они помечены бирюзовыми точками). Если треугольник не закрывает ни один из этих сэмплов, то пиксель считается чёрным, а если закрывает их всех, то белым. Здесь интересна ситуация, когда закрыта только часть пикселей: если закрыт один из четырёх, то пиксель будет на 25% белым и на 75% чёрным. В случае двух из четырёх соотношение 50/50, а при трёх закрытых сэмплах результатом будет более светлый оттенок в 75% белого.

Эта простая идея является фундаментом всех методов антиалиасинга на основе сэмплирования. В этом контексте также стоит заметить, что когда количество сэмплов на пиксель стремится к бесконечности, то результат этого процесса будет стремиться к «идеальному» сглаженному примеру, показанному ранее. Очевидно, что качество результата сильно зависит от количества использованных сэмплов — но и производительность тоже. Обычно в играх используется 2 или 4 сэмплов на пиксель, а 8 и более обычно применяются только в мощных PC.

Существуют и другие важные параметры, изменение которых может влиять на качество получаемых результатов методов антиалиасинга на основе сэмплирования. В основном это расположение сэмплов, тип сэмплов и группирование сэмплов.

Расположение сэмплов


Расположение сэмплов внутри пикселя сильно влияет на конечный результат, особенно в случае небольшого количества сэмплов (2 или 4), которое чаще всего используется в графике реального времени. В предыдущем примере сэмплы располагаются так, как будто они являются центрами отрендеренного изображения в четыре раза больше исходного (16×16 пикселей). Это интуитивно понятно и легко достигается простым рендерингом изображения бОльших размеров. Этот метод известен как антиалиасинг на упорядоченной сетке (ordered grid anti-aliasing, OGAA), также его иногда называют субдискретизацией (downsampling). В частности, его реализуют принудительным увеличением разрешения рендеринга по сравнению с разрешением монитора.

Однако упорядоченная сетка часто неоптимальна, особенно для почти вертикальных и почти горизонтальных линий, в которых как раз наиболее очевидны артефакты алиасинга. На Рисунке 6 показано, почему так происходит, и как повёрнутая или разреженная сетка сэмплирования обеспечивает гораздо лучшие результаты:


6-1. Сцена с почти вертикальной линией
6-2. Идеально сглаженная растеризация
6-3. Растеризация с четырьмя упорядоченными сэмплами
6-4. Сглаживание с четырьмя разреженными сэмплами

В этом почти вертикальном случае идеальный результат с четырьмя сэмплами должен иметь пять различных оттенков серого: черный при полностью незакрытых сэмплах, 25% белого при одном закрытом сэмпле, 50% при двух и так далее. Однако растеризация с упорядоченной сеткой даёт нам всего три оттенка: чёрный, белый и 50/50. Так происходит, потому что упорядоченные сэмплы расположены в два столбца, а потому, когда один из них закрывается почти вертикальным примитивом, другой тоже скорее всего будет закрыт.

Как показано на изображении с разреженным сэмплингом, эту проблему можно обойти, изменив положение сэмплов внутри каждого пикселя. Идеальным расположением сэмплов для сглаживания является разреженное. Это означает, что при N сэмплов никакие два сэмпла не имеют одного общего столбца, строки или диагонали в сетке NxN. Такие паттерны соответствуют решениям задачи об N ферзях. Методы антиалиасинга, в которых используются такие сетки, называют выполняющими антиалиасинг на разреженных сетках (sparse grid anti-aliasing, SGAA).

Типы сэмплов


Самый простой подход к антиалиасингу изображения на основе сэмплирования заключается в том, что все вычисления выполняются для «реального» пикселя каждого сэмпла. Хотя этот подход высокоэффективен для удаления всех типов артефактов алиасинга, он также является очень вычислительно затратным, потому что при N сэмплах увеличивает в N раз затраты на затенение, растеризацию, занимаемую полосу пропускания и память. Техники, при которых все вычисления выполняются для каждого отдельного сэмпла, называются сглаживанием суперсэмплингом (super-sampling anti-aliasing, SSAA).

Примерно в начале этого века в графическое оборудование была встроена поддержка антиалиасинга мультисэмплингом (multi-sample anti-aliasing, MSAA), являющегося оптимизацией суперсэмплинга. В отличие от случая SSAA, в MSAA каждый пиксель затеняется только один раз. Однако для каждого сэмпла вычисляются значения глубины и стенсила, что обеспечивает то же качество сглаживания на рёбрах геометрии, что и в SSAA, при значительно меньшем снижении производительности. Кроме того, возможны дальнейшие улучшения производительности, особенно занятой полосы пропускания, если поддерживается сжатие Z-буфера и буфера цвета. Они поддерживаются во всех современных архитектурах видеопроцессоров. Из-за способа оптимизации сэмплирования MSAA, с алиасингом прозрачности, текстур и шейдеров таким образом напрямую справляться невозможно.

Третий тип сэмплирования был представлен компанией NVIDIA в 2006 году в технологии антиалиасинг покрытия сэмплирования (coverage sampling anti-aliasing, CSAA). MSAA отделяет затенение от попиксельного вычисления глубины и стенсила, а CSAA добавляет сэмплы покрытия, которые не содержат значений цвета, глубины или стенсила — в них хранится только двоичное значение покрытия. Такие двоичные сэмплы используются для помощи в смешивании готовых сэмплов MSAA. То есть режимы CSAA добавляют сэмплы покрытия к режимам MSAA, но не имеет смысла выполнять сэмплирование покрытия без создания множества сэмплов MSAA. В современном оборудовании NVIDIA используется три режима CSAA: 8xCSAA (4xMSAA / 8 сэмплов покрытия), 16xCSAA (4x/16), 16xQCSAA (8x/16) и 32xCSAA (8x/32). У AMD есть похожая реализация с 4x EQAA (2x/4), 8xEQAA (4x/8) и 16xEQAA (8x/16). Дополнительные сэмплы покрытия обычно только незначительно влияют на производительность.

Группировка сэмплов


Последним ингредиентом методов AA на основе сэмплирования является способ группировки сэмплов, то есть то, как отдельные сэмплы, сгенерированные при рендеринге, собираются в конечный цвет каждого пикселя. Как показано на Рисунке 7, для этой цели используются различные фильтры группировки. На рисунке показаны пиксели 3×3 — бирюзовые точки обозначают позиции сэмплов, а жёлтый оттенок обозначает фильтр группировки сэмплов.
7-1. Фильтр Box
7-2. Фильтр Quincunx
7-3. Фильтр Tent

Очевидный и самый распространённый метод группировки просто накапливает каждый сэмпл в квадратной области, представляющей пиксель с равными весами. Это называется фильтром «box», и используется во всех обычных режимах MSAA.

Одним из первых подходов, пытавшихся улучшить эффект сглаживания с помощью малого количества сэмплов, является антиалиасинг «quincunx». В нём на пиксель вычисляется всего два сэмпла: один в центре, и один смещённый на полпикселя вверх и влево. Однако вместо этих двух сэмплов накапливается пять сэмплов, составляющих паттерн, показанный на Рисунке 7. Это приводит к значительному снижению алиасинга, но в то же время размывает всё изображение, потому что значения цветов окружающих пикселей группируются в каждый пиксель.

Более гибкий подход был представлен в 2007 году компанией AMD в серии видеопроцессоров HD 2900. В них используется программируемая группировка сэмплов, что позволяет реализовать режимы группировки «narrow tent» и «wide tent». Как показано выше, каждый сэмпл не имеет одинаковый вес. Вместо этого используется функция взвешивания, зависящая от расстояния до центра пикселя. Узкий (narrow) и широкий (wide) варианты используют разный размер ядра фильтра. Эти способы группировки можно сочетать с различным количеством сэмплов, и некоторые из полученных результатов показаны на общем сравнении. Что касается quincunx AA, то эти методы представляют собой компромисс между резкостью изображения и снижением алиасинга.

Сравнение AA сэмплирования


На Рисунке 8 показано сравнение всех рассмотренных нами методов AA на основе сэмплирования с различным количеством сэмплов. На изображении «ground truth» показано ближайшее к «реальному», идеальное представление сцены. Оно создано сочетанием 8xSGSSAA и 4×4 OGSSAA.

Стоит заметить аналогичное качество SGMSAA и SGSSAA с одинаковым количеством сэмплов при геометрическом алиасинге, и нехватку антиалиасинга прозрачности, текстур и шейдеров в случае MSAA. Недостатки упорядоченных паттернов сэмплинга, особенно для почти горизонтальных и почти вертикальных линий сразу заметны при сравнении 4x SGSSAA и 2×2 OGSSAA. При всего двух сэмплах на пиксель OGSSAA ограничен только горизонтальным (2×1) или только вертикальным (1×2) AA, а разреженный паттерн в какой-то мере может покрывать оба типа рёбер.

Методы AA с фильтрами группировки сэмплов, отличающиеся от обычного фильтра box, обычно обеспечивают более качественное снижение алиасинга на сэмпл, но страдают от эффекта размытия всего изображения.

Нужно заметить ещё один важный пункт — особенно в свете последующего обсуждения аналитических методов AA – все эти методы на основе сэмплирования одинаково хорошо применяются и к подпиксельному алиасингу, и обычному геометрическому алиасингу.

Рисунок 8: Обработка различных типов алиасинга различными методами AA на основе сэмплирования.


«Истинное» изображение
Без AA
2x MSAA
2x SGSSAA
4x MSAA
4x SGSSAA
8x MSAA
8x SGSSAA
8x MSAA + alpha-to-coverage
2x1 OGSSAA
1x2 OGSSAA
2x2 OGSSAA
4x Narrow Tent
6x Narrow Tent
6x Wide Tent
8x Wide TentТо же самое в виде GIF

Техники на основе сэмплинга интуитивно понятны и довольно хорошо работают с достаточно большим количеством сэмплов, но в то же время они затратны с точки зрения вычислений. Эта проблема усугубляется, когда используются способы рендеринга (например, отложенное затенение), которые могут усложнить применение эффективных типов сэмплирования с аппаратным ускорением. Поэтому исследуются другие способы уменьшения визуальных артефактов, создаваемых алиасингом при 3D-рендеринге. Такие методы рендерят обычное изображение с одним сэмплом на пиксель, а затем пытаются идентифицировать и устранить алиасинг, анализируя изображение.

Краткое введение и знакомство


Хотя идея сглаживания сгенерированных компьютером изображений популяризирована благодаря статье Решетова 2009 года о морфологическом антиалиасинге (который часто называют MLAA) [1], она ни в коем случае не является новой. Жуль Блументаль привёл сжатое описание этой техники в своей статье 1983 года для SIGGRAPH «Edge Inference with Applications to Antialiasing», которая активно применяется в современных методах [2]:
«Ребро, сэмплируемое по точкам для отображения на растровом устройстве, и непараллельное оси дисплея, выглядит как лесенка. Этот артефакт алиасинга часто возникает в компьютерных изображениях, сгенерированных двухмерными и трёхмерными алгоритмами. Точная информация о ребре часто больше недоступна, но из множества вертикальных и горизонтальных сегментов, формирующих эту лесенку, можно извлечь аппроксимацию исходного ребра с точностью, превосходящей точность растра. Таким образом можно обеспечить сглаживание ребра лесенки.

Такие извлечённые рёбра можно использовать для повторного затенения пикселей, которые они пересекают, таким образом сглаживая извлечённые рёбра. Сглаженные извлечённые рёбра оказываются более привлекательной аппроксимацией реальных рёбер, чем из аналоги с алиасингом».


В 1999 году Ишики и Куниэда представили первый вариант этой техники, предназначенной для использования в реальном времени, которая выполнялась сканированием пар строк и столбцов изображения, и могла быть реализована аппаратно [3].

В общем случае, все чисто аналитические методы антиалиасинга выполняются в три этапа:

  1. Распознавание разрывов в изображении.
  2. Воссоздание геометрических рёбер из паттерна разрывов.
  3. Сглаживание пикселей, пересекающих эти рёбра, смешиванием цветов каждой их стороны.

Отдельные реализации аналитического антиалиасинга различаются тем, как реализованы эти шаги.

Распознавание разрывов


Простейший и самый распространённый вариант распознавания разрывов — это простое изучение конечного отрендеренного цветового буфера. Если разность цветов двух соседних пикселей (их расстояние) больше какого-то порогового значения, то присутствует разрыв, в противном случае его нет. Эти показатели расстояний часто вычисляются в цветовом пространстве, которое лучше моделирует человеческое зрение, чем RGB, например, в HSL.

На Рисунке 9 показан пример отрендеренного изображения, а также вычисленных из них горизонтальных и вертикальных разрывов.


Рисунок 9: Распознавание разрывов в цветовом буфере. Слева: изображение без AA. В центре: горизонтальные разрывы. Справа: вертикальные разрывы.

Чтобы ускорить процесс распознавания разрывов или снизить количество ложноположительных распознаваний (например, в текстуре, или вокруг текста на Рисунке 9), можно использовать другие буферы, генерируемые в процессе рендеринга. Обычно для прямого и обратного рендерера доступен Z-буфер (буфер глубин). В нём хранится значение глубины для каждого пикселя, и его можно использовать для распознавания рёбер. Однако это только работает для распознавания рёбер-силуэтов, то есть внешних рёбер 3D-объекта. Чтобы рассматривать и рёбра внутри объекта, то нужно использовать ещё один буфер вместо или в дополнение к Z-буферу. Для отложенных рендереров часто генерируется буфер, хранящий направление поверхностей нормалей каждого пикселя. В таком случае для распознавания рёбер подходящей метрикой будет угол между соседними нормалями.

Воссоздание рёбер и смешивание


Способ воссоздания геометрических рёбер из разрывов немного различается в разных методах аналитического AA, но все они выполняют похожие действия по сопоставлению паттернов на горизонтальных и вертикальных разрывах для распознавания типичного паттерна «лесенки» артефактов алиасинга. На Рисунке 10 показаны паттерны, используемые по описанию Решетова в MLAA и способ воссоздания из них рёбер.

Рисунок 10: Паттерны MLAA и их воссозданные рёбра


Распознанные паттерны
Паттерны разрывов, использованные в MLAA

После воссоздания геометрических рёбер достаточно просто вычислить, насколько верхний/нижний или правый/левый пиксель вдоль ребра должен вносить вклад в смешанный цвет пикселя, чтобы генерировать сглаженный внешний вид.

Достоинства и недостатки аналитического сглаживания


По сравнению с методами антиалиасинга на основе сэмплирования, аналитические решения имеют несколько важных преимуществ. При правильной работе (на правильно распрознанных геометрических рёбрах) они могут обеспечить качество, равное качеству методов на основе сэмплирования с очень высоким количеством сэмплов, при этом тратя меньше вычислительных ресурсов. Более того, они легко применимы во множестве случаев, в которых AA на основе сэмплирования реализуется сложнее, например, в случае отложенного затенения.

Однако аналитический AA — это не панацея. Неотъемлемой проблемой чисто аналитических методов на основании единственного сэмпла является то, что они не могут справляться с подпиксельным алиасингом.

Если передать аналитическому алгоритму сглаживания структуру пикселей, показанную в правом верхнем углу Рисунка 2, то он никак не сможет понять, что разделённые группы пикселей на самом деле составляют линию. На этом этапе есть два равно неприятных способа решения этой проблемы: или размытие пикселей, снижающее видимый алиасинг, но и уничтожающее детали, или консервативная обработка только чётко определённых артефактов-«лесенок», которые совершенно точно вызваны алиасингом; при этом подпиксельный алиасинг сохранится и будет искажён.

Ещё одна проблема аналитических методов заключается в ложных срабатываниях. Когда часть изображения распознаётся как ребро с алиасингом, но на самом деле им не является, оно будет искажено смешиванием. Это особенно очевидно на тексте, и тоже требует идти на компромиссы: более консервативное распознавание рёбер приведёт к меньшему количеству ложных срабатываний, но при нём также будут упущены некоторые рёбра, которые на самом деле имеют алиасинг. С другой стороны, при расширении порога распознавания рёбер эти рёбра тоже будут включены, но это приведёт к ложноположительным срабатываниям. Поскольку анатилический антиалиасинг в своей основе пытается экстраполировать больше информации из растеризированного изображения, чем на самом деле в нём есть, от этих проблем невозможно избавиться полностью.

Наконец, интерпретация рёбер этими методами может сильно меняться в зависимости от разницы в единственном пикселе. Поэтому при использовании однопиксельных аналитических методов антиалиасинга может увеличиться или даже добавиться мерцание и временная нестабильность изображения: единственный изменившийся пиксель в исходном изображении может превратиться в сглаженных выходных данных в целую мерцающую линию.

На Рисунке 11 показаны некоторые из успешных и неудачных результатов использования аналитических методов AA на примере стандартных алгоритмов FXAA и SMAA1x. Последний обычно считается лучшим чисто аналитическим однопиксельным алгоритмом, который можно использовать в реальном времени.

Рисунок 11: Аналитические методы AA


Без AA
FXAA
1x SMAA

Сравнение аналитических методов антиалиасинга


На Рисунке 12 показано сравнение между результатами работы FXAA, SMAA1x и «идеальным» изображением, и изображениями без AA, с 4xMSAA и 4xSGSAAA из предыдущего сравнения.

Рисунок 12: Обработка различных видов алиасинга разными аналитическими и сэмплирующими методами AA


«Идеальное» изображение
Без AA
4x MSAA
4x SGSSAA
FXAA
SMAAТо же самое в виде GIF
Заметьте, что в отличие от MSAA, эти аналические методы не интересует, были ли причинами артефактов алиасинга геометрия, прозрачность или даже вычисление шейдеров. Все рёбра обрабатываются одинаково. К сожалению, то же самое относится и к рёбрам экранного текста, хотя искажение при SMAA1x и меньше, чем при FXAA.

Оба метода не справляются с антиалиасингом в случае подпиксельного алиасинга, однако они обрабатывают эту неудачу по-разному: SMAA1x просто решает вообще не влиять на отдельные белые пиксели синусоиды, а FXAA смешивает их с их окружением. Более желательная обработка зависит от контекста и личных предпочтений.

Более объективно можно сказать, что SMAA1x обрабатывает некоторые углы линий в тесте 2D-геометрии и кривые в примере с шейдерным алиасингом точно лучше, чем FXAA, обеспечивая более гладкий результат, который приятнее эстетически и ближе к «идеальному» изображению. Так получилось благодаря более сложному этапу воссоздания рёбер и смешивания, подробности которого объяснены в статье об SMAA 2012 года Хименеса et alia [4].


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

Ну, будущее уже частично здесь: комбинированные методы на основе сэмплирования и аналитического антиалиасинга. Эти алгоритмы создают несколько сэмплов сцены — посредством традиционного мульти- или суперсэмплирования или с помощью временного накопления между кадрами — и сочетают их с анализом, генерируя конечное изображение с антиалиасингом. Это позволяет им снизить проблемы подпиксельного алиасинга и временной нестабильности односэмпловых чисто аналитических методов, однако всё равно даёт гораздо лучшие результаты на геометрических рёбрах, чем чистые методы сэмплирования со схожими характеристиками производительности. Очень простую комбинацию дополнительного сэмплирования и аналитического AA можно получить сочетанием односэмпловой аналитической техники наподобие FXAA с субдискретизацией из буфера с более высоким разрешением. Более сложными примерами таких методов являются SMAA T2x, SMAA S2x и SMAA 4x, а также TXAA. Методы SMAA подробно объясняются и сравниваются в этой статье, в то время как NVIDIA реализовала собственный подход к TXAA здесь. Высока вероятность, что такие комбинированные методы будут более широко использоваться в будущем.

Ещё одним вариантом, пока не получившим широкого распространения, но имеющим большой потенциал на будущее, является кодирование в процессе рендеринга дополнительной геометрической информации, которая позже будет использоваться на этапе аналитического антиалиасинга. Примеры такого подхода — это антиалиасинг с геометрической постобработкой (Geometric Post-process Anti-Aliasing, GPAA) и антиалиасинг с использованием буфера геометрии (Geometry Buffer Anti-Aliasing, GBAA), демо которых выложены здесь.

Наконец, общий пул памяти ЦП и видеопроцессора новых консольных платформ и будущих архитектур PC могут позволить использовать техники, предназначенные для эксплуатации таких общих ресурсов. В недавней статье «Asynchronous Adaptive Anti-aliasing Using Shared Memory» Бэрринджер и Мёллер описывают технику, выполняющую традиционный односэмпловый рендеринг, в то же время распознавая важные пиксели (например, находящиеся на ребре) и растеризируя для них в ЦП дополнительные разреженные сэмплы [5]. Хотя это требует серьёзной реструктуризации процесса выполнения рендеринга, результаты выглядят многообещающе.

Справочные материалы


[1] A. Reshetov, «Morphological antialiasing», в Proceedings of the Conference on High Performance Graphics 2009 HPG ’09, New York, NY, USA, 2009, pp. 109–116.

[2] J. Bloomenthal, ‘Edge Inference with Applications to Antialiasing’, ACM SIGGRAPH Comput. Graph., vol. 17, no. 3, pp. 157–162, Jul. 1983.

[3] T. Isshiki and H. Kunieda, ‘Efficient anti-aliasing algorithm for computer generated images’, в Proceedings of the 1999 IEEE International Symposium on Circuits and Systems ISCAS ’99, Orlando, FL, 1999, vol. 4, pp. 532–535.

[4] J. Jimenez, J. I. Echevarria, T. Sousa, and D. Gutierrez, ‘SMAA: Enhanced Subpixel Morphological Antialiasing’, Comput. Graph. Forum, vol. 31, no. 2pt1, pp. 355–364, May 2012.

[5] R. Barringer and T. Akenine-Möller, ‘A 4: asynchronous adaptive anti-aliasing using shared memory’, ACM Trans. Graph., vol. 32, no. 4, pp. 100:1–100:10, Jul. 2013.

Цифровая фильтрация на ПЛИС – Часть 1 / Хабр

Всем привет!

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

В этих статьях я постараюсь минимизировать теоретическое описание тех или иных алгоритмов и большую часть материала посвятить практическим тонкостям, с которыми столкнулся лично я и мои коллеги, и знакомые, так или иначе связанные с разработкой на ПЛИС. Надеюсь, данный цикл статей принесет пользу, как начинающим инженерам, так и матёрым разработчикам.

Часть 1: CIC фильтр

В первой части рассмотрим простейший CIC фильтр. CIC – «cascaded integral-comb», по-русски – каскадный интегрально-гребенчатый фильтр типа БИХ (с бесконечной импульсной характеристикой). Класс таких фильтров широко используется в задачах, где требуется работа на нескольких скоростях передачи данных. CIC фильтры активно применяются для децимации и интерполяции, то есть для понижения и повышения частоты дискретизации. CIC фильтр сам по себе есть не что иное, как фильтр нижних частот (ФНЧ). То есть такой фильтр пропускает нижние частоты спектра, обрубая верхние за частотой среза. АЧХ фильтра строится по закону ~sin(x)/x. Главное преимущество CIC фильтров состоит в том, что они совсем не требуют операций умножения (в отличие от другого типа фильтров, например, КИХ).
Введение

Из названия можно догадаться, что в основе CIC фильтра лежит два базовых блока: интегратор и гребенчатый фильтр (дифференциатор). Интегрирующее звено (int) представляет собой обычный БИХ-фильтр первого порядка, выполненный как самый простой аккумулятор. Гребенчатый фильтр (comb) является КИХ-фильтром первого порядка.

Между интегратором и гребенчатым фильтром часто ставится узел повышения или понижения частоты дискретизации в целое число раз — R.

  • В случае понижения частоты дискретизации из входной последовательности выбирается каждый R-отсчет, образуя прореженную выходную последовательность.
  • В случае повышения частоты дискретизации между отсчетами входной последовательности просто вставляются нули, которые затем сглаживаются в интегрирующей секции, образуя последовательность на увеличенной частоте дискретизации.

Формулы для передаточной и амплитудно-частотной характеристик приведены ниже:

Более подробно со всеми математическими выкладками обо всех аспектах децимации и интерполяции можно почитать в других источниках, на которые в конце статьи я приведу ссылки.
Дециматор

Если CIC-фильтр используется для понижения частоты дискретизации, то он называется дециматором. В таком случае первым звеном идет интегратор, затем происходит понижение частоты дискретизации и, наконец, идет звено дифференцирующего фильтра.
Интерполятор

Если CIC-фильтр используется для повышения частоты дискретизации, то он называется интерполятором. В таком случае дифференцирующее звено стоит на первом месте, затем происходит повышение частоты дискретизации и, наконец, идет звено интегрирующего фильтра.

В зависимости от задержки входного сигнала в дифференцирующем звене, можно получать различные частотные характеристики фильтра. Известно, что при увеличении параметра задержки D, увеличивается количество «нулей» амплитудно-частотной характеристики (АЧХ) фильтра.

Заметим, что для связки интегратора и гребенчатого фильтра (CIC фильтра) при увеличении параметра D в дифференцирующей секции нули АЧХ смещаются к центру – изменяется частота среза фильтра Fc = 2 pi / D.

Каскадное соединение интегратора и гребенчатого фильтра без операций децимации и интерполяции называется фильтром «скользящего среднего». Уровень первого бокового лепестка такого фильтра составляет всего -13 дБ, что достаточно мало для серьезных задач ЦОС.
В силу линейности математических операций, происходящих в CIC фильтре возможно каскадное соединение нескольких фильтров подряд. Это дает пропорциональное уменьшение уровня боковых лепестков, но также увеличивает завал главного лепестка спектра (под спектром я часто буду понимать АЧХ фильтра). Таким образом, при N-каскадном соединении однотипных CIC фильтров происходит перемножение идентичных передаточных характеристик. Как правило, секции интеграторов и гребенчатых фильтров объединяются вместе по типу. Например, сначала последовательно ставится N секций однотипных интеграторов, затем N секций однотипных дифференцирующих фильтров.

На следующем рисунке приведена АЧХ фильтра при различных параметрах коэффициента дискретизации R (расчет сделан в MathCAD 14).

АЧХ CIC фильтра полностью эквивалентна частотной характеристике FIR фильтра с прямоугольной импульсной характеристикой (ИХ). Общая ИХ фильтра определяется как свертка всех импульсных характеристик каскадов связки интегратора и гребенчатого фильтра. С ростом порядка CIC фильтра, его ИХ интегрируется соответствующее число раз. Таким образом, для CIC фильтра первого порядка ИХ – прямоугольник, для фильтра второго порядка ИХ – равнобедренный треугольник, для третьего порядка ИХ – парабола и т.д.

Рост разрядности данных

К несчастью, увеличение величины задержки D в гребенчатой структуре и увеличение порядка фильтра N приводят к росту коэффициента передачи. Это в свою очередь приводит к увеличению разрядности на выходе фильтра. В задачах ЦОС, где применяются CIC фильтры нужно всегда об этом помнить и следить, чтобы передаваемые сигналы не выходили за используемую разрядную сетку. К примеру, негативный эффект роста разрядности проявляется в значительном увеличении используемых ресурсов кристалла ПЛИС.

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

Дециматор: CIC фильтр-дециматор очень чувствителен к параметрам D, R и N, от которых зависит разрядность промежуточных и выходных данных. И дифференцирующее звено, и интегратор влияют на конечную разрядность выходного сигнала.

В этих формулах: B — разрядность входных данных, Bmax — разрядность выходных данных, R — коэффициент дискретизации, D — параметр задержки, N — порядок фильтра (количество каскадов).

Замечание! В статье Хогенауэра описаны принципы выбора разрядности для каждого каскада дециматора. Xilinx и Altera при реализации своих фильтров учитывают негативный эффект роста разрядности фильтра и борятся с этим явлением методами, описанными в статье.

Xilinx CIC Filter

Так как я моя работа на 99% связана с микросхемами фирмы Xilinx, я приведу описание IP-ядра фильтра для этого вендора. Но смею вас заверить, что для Altera все практически аналогично.
Для того, чтобы создать CIC фильтр, необходимо зайти в приложение CORE Generator и создать новый проект, в котором указать тип используемого кристалла ПЛИС и различные другие несущественные в данном случае настройки.

CIC Compiler — Вкладка 1:

Component name имя компонента (используются латинские буквы a-z, цифры 0-9 и символ "_").

Filter Specification:

  • Filter type — тип фильтра: интерполирующий / децимирующий.
  • Number of stages — количество каскадов интеграторов и гребенчатых фильтров: 3-6.
  • Differential delay — задержка в дифференциальных ячейках фильтра: 1-2.
  • Number of channels — количество независимых каналов: 1-16.

Sample Rate Change Specification:
  • Fixed / Programmable — тип коэффициента дискретизации R: постоянный / программируемый.
  • Fixed or Initial Rate — значение коэффициента дискретизации R: 4..8192.
  • Minimum Rate — минимальное значение коэффициента дискретизации R: 4..8.
  • Maximum Rate — максимальное значение коэффициента дискретизации R: 8..8192.

Hardware Oversampling Specification: эти параметры влияют на выходную частоту дискретизации, количество тактов, требуемых для обработки данных. От этих параметров также зависит уровень параллелизма внутри ядра и количество занимаемых ресурсов.
  • Select format — выбор частотных соотношений фильтра: Frequency Specification / Sample period.
  • Frequency Specification — Частотная спецификация: пользователь задает частоту дискретизации и частоту обработки данных.
  • Sample period — Тактовая спецификация: пользователь задает отношение частоты обработки к тактовой частоте данных.
  • Input Sampling Frequency — входная частота дискретизации: *.
  • Clock frequency — частота обработки фильтра: *.
  • Input Sampling period — отношение частоты обработки к частоте входного тактового сигнала: *.

* — диапазон зависит от общих настроек и коэффициента дискретизации R.

CIC Compiler — Вкладка 2:

Numerical Precision:

  • Input Data Width — разрядность входных данных: 2..20.
  • Quantization — округление выходных данных: полная точность / округление разрядной сетки.
  • Output Data Width — разрядность выходных данных, диапазон зависит от коэффициентов N, D и R (максимальное значение — 48 битов).

Optional:
  • Use Xtreme DSP Slice — использовать встроенные DSP-блоки для реализации фильтра.
  • Use Streaming Interface — использовать потоковый интерфейс для многоканальной реализации фильтра.

Control Options:
  • ND — «New data», входной сигнал, определяющий поступление данных на вход фильтра.
  • SCLR — синхронный сброс фильтра (логическая единица на этом входе производит сброс).
  • CE — «Clock Enable», сигнал разрешения тактирования фильтра.

CIC Compiler — Вкладка 3:

Summary — эта вкладка в виде списка отражает конечные настройки фильтра (количество каскадов, параметры частот, разрядность входных, выходных и промежуточных данных, задержка в фильтре и т.д.).

В левой части окна CIC Compiler есть три полезные дополнительные вкладки:

  • IP-symbol — схематичный вид IP-блока с активными портами ввода/вывода.
  • Freq. response — передаточная характеристика CIC-фильтра.
  • Resource estimate — оценка занимаемых ресурсов.

После установки всех настроек необходимо нажать на кнопку Generate. В результате приложение CORE Generator через какое-то время выдаст целый набор файлов, из которых нам нужны самые основные:
  • *.VHD (или *.V) — файл исходных кодов для моделирования на VHDL или Verilog.
  • *.VHO — бесполезный файл, но из него можно взять описание компонента и портирование для вставки в проект.
  • *.NGC — файл списка соединений. Содержит описание архитектуры IP-ядра (используемые компоненты и связи сигналов между ними) для выбранного кристалла ПЛИС.
  • *.XCO — лог-файл, в котором хранятся все параметры и настройки IP-ядра. Полезный файл при работе в среде Xilinx ISE Design Suite.

Если вы работаете в среде ISE Design Suite, то CORE Generator автоматически создаст нужные файлы в рабочем каталоге. Для других средств разработки (типа Modelsim или Aldec Active-HDL) необходимо перенести нужные файлы в соответствующий рабочий каталог.
CIC Filter in MATLAB

Пример 1: Для моделирования очень удобным средством является программа MATLAB. Для примера возьмем модель CIC-фильтра 4 порядка, сделанного на логических элементах из System Generator Toolbox от Xilinx. Децимация и интерполяция не используется (CIC вырождается в фильтр скользящего среднего с окном 16). Параметры фильтра: R = 1, N = 4, D = 16. На следующем рисунке приведена модель одного каскада в среде MATLAB.

Посмотрим, как выглядит импульсная характеристика после каждого каскада фильтра, для этого подадим на вход системы периодический единичный импульс.

Видно, что сигнал на выходе первого звена образует прямоугольный импульс длительностью = D, на выходе второго звена — треугольный сигнал длительностью 2D, на выходе третьего звена — параболический импульс, на выходе третьего — кубическая парабола. Результат полностью согласуется с теорией.

Пример 2: непосредственно IP-ядро CIC фильтра. Параметры: N = 3, R = 4, D = 1. На следующем рисунке представлена модель фильтра.

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

Резюме

На этом хотелось бы подвести итог. CIC фильтры используются во многих задачах, где требуется изменить частоту дискретизации. CIC фильтры применяются в системах, работающих на нескольких частотах дискретизации (multirate processing), например, в аудио-технике для изменения бит-рейта (из 44.1кГц в 48кГц и обратно). CIC фильтры применяются в системах связи для реализации DDC (digital down converter) и DUC (digital up converters). Пример использования CIC-фильтров: микросхема цифрового приема AD6620 от Analog Devices.

Реализация собственного фильтра на ПЛИС на HDL языках часто не требуется, и можно смело пользоваться готовыми ядрами от вендоров, либо готовыми opensource-проектами. Если все же возникла необходимость в реализации собственного CIC фильтра для прикладной задачи, то нужно помнить следующие принципы.

CIC фильтры обладают рядом особенностей:

  1. Простые в реализации и не требуют операций умножения.
  2. Децимация и интерполяция на CIC фильтрах используется повсеместно для быстрого изменения частоты дискретизации, как в целое, так и в дробное число раз.
  3. С ростом порядка фильтра N и величины задержки D увеличивается разрядность промежуточных и выходных данных.
  4. С ростом порядка фильтра N увеличивается подавление боковых лепестков и увеличивается неравномерность главного лепестка АЧХ.
  5. Рекомендуется использовать фильтры порядка не более 6-8, т.к. с увеличением порядка усложняется реализация, увеличивается объем занимаемых ресурсов, а также происходят искажения АЧХ фильтра в пределах полосы пропускания.
  6. С ростом параметра задержки D гребенчатого фильтра изменяется частота среза фильтра, но в практических целях при каскадном соединении параметр D < 3.
  7. При децимации в R раз существенно увеличивается разрядность на выходе фильтра.
  8. При интерполяции основной вклад в разрядность промежуточных и выходных данных вносят только интегрирующие звенья.
  9. АЧХ CIC фильтра эквивалентна АЧХ FIR фильтра с прямоугольной импульсной характеристикой. Общая ИХ фильтра определяется как свертка всех импульсных характеристик каскадов связки интегратора и гребенчатого фильтра.
  10. При изменении частоты на выходе фильтра в ПЛИС используют сигнал разрешения «clock enable», а частоту обработки не изменяют.
  11. Если отношение «частота обработки / частота дискретизации» >> 1, в ПЛИС возможно повторно использовать ресурсы фильтра, тем самым для многоканальной системы реализовать обработку с минимальной затратой ресурсов кристалла.
  12. В современных ПЛИС CIC фильтры реализуются на блоках DSP (Xilinx, Altera), но при отсутствии свободных ресурсов возможна реализация на логических ячейках (SLICEs).
  13. После CIC фильтра рекомендуется ставить умножитель с программируемым коэффициентом усиления (gain multiplier), который будет подстраивать уровень сигнала до нужного динамического диапазона
  14. CIC фильтры вносят искажение в спектр выходного сигнала, поэтому после CIC фильтра необходимо ставить компенсирующий FIR фильтр (методика расчета представлена в даташите Altera, для расчета необходим MATLAB).

Литература

Продолжение следует...

Обработка сигналов (scipy.signal) - SciPy v1.5.2 Справочное руководство

билинейный (b, a [, fs])

Возврат цифрового БИХ-фильтра из аналогового с помощью билинейного преобразования.

билинейный_зпк (z, p, k, fs)

Возврат цифрового БИХ-фильтра из аналогового с помощью билинейного преобразования.

findfreqs (число, ден, N [, вид])

Найдите массив частот для вычисления отклика аналогового фильтра.

firls (numtaps, band, желаемый [, вес, nyq, fs])

Разработка КИХ-фильтра с использованием минимизации ошибок методом наименьших квадратов.

firwin (numtaps, cutoff [, width, window,…])

Построение КИХ-фильтра оконным методом.

firwin2 (numtaps, freq, gain [, nfreqs,…])

Построение КИХ-фильтра оконным методом.

частот (b, a [, работа, график])

Вычислить частотную характеристику аналогового фильтра.

freqs_zpk (z, p, k [, worN])

Вычислить частотную характеристику аналогового фильтра.

частота (b [, a, worN, целое, график, fs,…])

Вычислить частотную характеристику цифрового фильтра.

freqz_zpk (z, p, k [, worN, whole, fs])

Вычислить частотную характеристику цифрового фильтра в форме ZPK.

СОСФРЕКЗ (СОС [, ТОРГ, ЦЕЛЫЙ, ФС])

Вычислить частотную характеристику цифрового фильтра в формате SOS.

group_delay (system [, w, whole, fs])

Вычислить групповую задержку цифрового фильтра.

iirdesign (wp, ws, gpass, gstop [, аналог,…])

Полная конструкция цифрового и аналогового фильтра БИХ.

iirfilter (N, Wn [, rp, rs, btype, аналог,…])

Конструкция цифрового и аналогового БИХ-фильтров для заданного порядка и критических точек.

kaiser_atten (число, ширина)

Вычислить ослабление КИХ-фильтра Кайзера.

kaiser_beta (а)

Вычислить параметр Кайзера beta , учитывая затухание a .

kaiserord (волнистость, ширина)

Определите параметры окна фильтра для метода окна Кайзера.

minimum_phase (h [, method, n_fft])

Преобразование линейно-фазового КИХ-фильтра в минимальную фазу

savgol_coeffs (длина_окна, polyorder [,…])

Вычислить коэффициенты для одномерного КИХ-фильтра Савицки-Голея.

Remez (числа, диапазоны, желаемый [, вес, Гц,…])

Рассчитайте минимаксный оптимальный фильтр, используя алгоритм обмена Ремеза.

unique_roots (p [, tol, rtype])

Определите уникальные корни и их кратности из списка корней.

остаток (b, a [, tol, rtype])

Вычислить частичное дробное расширение b (s) / a (s).

остатки (b, a [, tol, rtype])

Вычислить частичное дробное разложение b (z) / a (z).

invres (r, p, k [, tol, rtype])

Вычислить b (s) и a (s) по частичному расширению фракции.

invresz (r, p, k [, tol, rtype])

Вычислить b (z) и a (z) по частичному расширению фракции.

Плохие коэффициенты

Предупреждение о плохо согласованных коэффициентах фильтра

.

Справка в Интернете - Учебники - Сглаживание

Сглаживание

Сводка

Сглаживание - распространенный метод удаления шума из сигналов. Origin предоставляет несколько методов сглаживания, включая усреднение по соседству, Савицки-Голея, процентильный фильтр и фильтр БПФ. Кроме того, доступен инструмент на основе вейвлетов.

Для матричных данных сглаживание Origin выполняется двумя способами. Если количество столбцов или строк меньше 32, матрица сначала расширяется, а затем возвращается к исходному размеру.Если больше 31, матрица сначала сжималась, а затем расширялась.

Что вы узнаете

Из этого туториала Вы узнаете, как:

  • Сглаживание сигналов несколькими способами.
  • Сглаживание сигналов с использованием вейвлетов.
  • Сгладьте 3D-поверхность.

Сглаживание несколькими методами

  1. Начать с пустой книги.
  2. Выберите пункт меню Данные: Импорт из файла: Отдельный ASCII ... , чтобы импортировать данные Сигнал с высокочастотным шумом.dat в папке \ Samples \ Signal Processing \ .
  3. Выделите столбец B на листе. Затем выберите пункт меню Analysis: Signal Processing: Smooth ... (если инструмент использовался ранее, выберите Open Dialog ), чтобы открыть диалоговое окно Smooth .

Усреднение по соседнему

Метод усреднения по соседству выполняет широкое сглаживание данных.

  1. Выберите Усреднение по соседству для метода . Установите Точки окна на 100 и Граничное условие на Периодическое . Чтобы просмотреть результат на правой панели, установите флажок Auto Preview .
  2. Нажмите ОК , чтобы сгенерировать результат.

Савицкий-Голай

Метод Савицкого-Голея хорош для сохранения формы пиков сигнала.

  1. Снова выберите столбец B. В меню Analysis щелкните Signal Processing: Smooth: Open Dialog ....
  2. В диалоговом окне сглаживание установите Method на Savizky-Golay . Установите Точки окна на 100 , Граничное условие на Периодический и полиномиальный порядок с на 3 .
  3. Щелкните ОК .

Процентильный фильтр

  1. Снова выберите столбец B.В меню Анализ щелкните Обработка сигнала: Сглаживание: открыть диалог .... .
  2. Для Method выберите Percentile Filter . Установите Points of Window на 100 , Boundary Condition to Percentile , а для Percentile примите значение по умолчанию 50 .
  3. Нажмите ОК кнопку.

Фильтр БПФ

  1. Снова выберите столбец B.В меню Анализ щелкните Обработка сигнала: Сглаживание: открыть диалог .... .
  2. Установите Method на FFT Filter . Установите точек окна на 100 .
  3. Нажмите ОК , чтобы закрыть диалоговое окно.
  4. Обратите внимание, что на листе теперь есть четыре добавленных столбца, содержащих результаты ваших четырех операций сглаживания. Выделите эти четыре столбца (C, D, E, F), затем в меню выберите Plot> 2D: Line: Line , чтобы построить линейный график с этими четырьмя наборами данных.
  5. Сравнивая результаты четырех методов, вы можете видеть, что метод Савицкого-Голея лучше всего справлялся с сохранением пиков в данных, в то время как метод фильтра БПФ справлялся хуже всего с сохранением пиков.

Сглаживание с помощью вейвлета

  1. Начать новую книгу.
  2. Выберите пункт меню Data: Import from File: Single ASCII ... , чтобы импортировать данные Signal with High Frequency Noise.dat в папку \ Samples \ Signal Processing \ .
  3. Выделите столбец B и в меню выберите Analysis: Signal Processing: Wavelet: Smooth ... , чтобы открыть диалоговое окно Smooth: wtsmooth .
  4. В диалоговом окне установите Wavelet Type на DB6 и Cutoff (%) на 98 . Установите флажок Auto Preview , чтобы просмотреть результат на правой панели.
  5. Нажмите ОК , чтобы закрыть диалоговое окно и сгенерировать результат.
  6. Чтобы увидеть разницу между исходным сигналом и сглаженным сигналом, выделите все столбцы на листе и в меню выберите График> 2D: Линия: Линия .

Матрица сглаживания

  1. Начать с новой матричной книги.
  2. Выберите меню Data: Import from File: Image to Matrix ... , чтобы импортировать изображение scale.jpg в папке \ Samples \ Image Processing and Analysis \ (если диалоговое окно impImage открывается, просто примите значения по умолчанию и нажмите OK ).
  3. Origin не может применяться для выполнения сглаживания изображения, поэтому мы должны сначала преобразовать изображение в матричные данные. В меню выберите Изображение: Преобразование: Преобразовать в данные ... . Примите настройки диалогового окна по умолчанию и нажмите ОК .
  4. При активированной преобразованной матрице выберите меню View: Image Mode , чтобы просмотреть эту матрицу как изображение.
  5. Чтобы выполнить сглаживание, выберите Анализ: Обработка сигнала: Сглаживание... . Откроется диалоговое окно Smooth: msmooth .
  6. Примите настройки по умолчанию и нажмите ОК . Чтобы просмотреть результат, выберите в меню View: Image Mode .
.Сглаживание сигнала

- это ... Что такое сглаживание сигнала?

  • Обработка сигналов - это область системной инженерии, электротехники и прикладной математики, которая занимается операциями или анализом сигналов в дискретном или непрерывном времени. Представляющие интерес сигналы могут включать звук, изображения, переменные во времени измерения… Wikipedia

  • Сглаживание - В статистике и обработке изображений сглаживание набора данных означает создание аппроксимирующей функции, которая пытается уловить важные закономерности в данных, исключая шум или другие мелкомасштабные структуры / быстрые явления.Много разных…… Википедия

  • Численное сглаживание и дифференцирование - Экспериментальное значение нулевой точки можно концептуально описать как сумму сигнала и некоторого шума, но на практике эти два вклада нельзя разделить. Целью сглаживания является увеличение отношения сигнал / шум без значительного…… Wikipedia

  • Экспоненциальное сглаживание - это метод, который может применяться к данным временных рядов, либо для получения сглаженных данных для представления, либо для составления прогнозов.Сами данные временного ряда представляют собой последовательность наблюдений. Наблюдаемое явление может быть по сути случайным…… Wikipedia

  • Ложный сигнал - В техническом анализе ложный сигнал относится к указанию будущих ценовых движений, что дает неточное представление об экономической реальности. Ложные сигналы могут возникать из-за ряда факторов, в том числе временных задержек, неточностей в данных…… Инвестиционный словарь

  • Усреднение обобщенного сигнала - Во многих случаях доступно только одно изображение с шумом, и затем выполняется усреднение в локальной окрестности.Результаты приемлемы, если шум меньше по размеру, чем самые маленькие интересующие объекты на изображении, но размытие краев - это…… Wikipedia

  • Реализация пространства масштабирования - Масштабное пространство Аксиомы масштабируемого пространства Реализация масштабного пространства Обнаружение функций Обнаружение краев Обнаружение больших двоичных объектов Обнаружение углов… Wikipedia

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

  • Масштабное пространство - теория представляет собой основу для многомасштабного представления сигналов, разработанную сообществами компьютерного зрения, обработки изображений и сигналов с дополнительными мотивами из физики и биологического видения. Это формальная теория обработки… Wikipedia

  • Хендрик Уэйд Боде - Инфобокс Имя ученого = Хендрик Уэйд Боде ширина изображения = 162 caption = Хендрик Уэйд Боде дата рождения = дата рождения | 1905 | 12 | 24 | df = y дата смерти = дата и возраст смерти | 1982 | 6 | 21 | 1905 | 12 | 24 | df = y место рождения = Мэдисон, штат Висконсин место смерти = Кембридж,…… Википедия

  • Скользящее среднее - Для использования в других целях см. Скользящее среднее (значения).В статистике скользящее среднее, также называемое скользящим средним, скользящим средним или скользящим средним, представляет собой тип фильтра с конечной импульсной характеристикой, используемый для анализа набора точек данных путем создания…… Wikipedia

  • .

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

    Ваш адрес email не будет опубликован. Обязательные поля помечены *