Введение
Данную статью меня заставил написать пост habrahabr.ru/post/183986, где не совсем правильно используется некоторый алгоритм сглаживания изображения.
Сразу перейдём к сути дела.
Математические модели цифровых сигналов — вектора и матрицы, элементами которых являются числа. Числа могут быть двоичными (бинарный сигнал), десятичными («обычный» сигнал) и так далее. Любой звук, любое изображение и видео могут быть преобразованы в цифровой сигнал1: звук — в вектор, изображение — в матрицу, а видео — в последовательный набор матриц. Поэтому цифровой сигнал — это, можно сказать, универсальный объект для представления информации.
Задача сглаживания — это, по сути, задача фильтрации сигнала от скачкообразных (ступенчатых) изменений. Считается, что полезный сигнал их не содержит. Ступенчатый сигнал за счёт множества резких, но небольших по амплитуде, перепадов уровня содержит высокочастотные составляющие, которых нет в сглаженном сигнале. Поэтому для некоторого алгоритма сглаживания в первую очередь необходимо определить как сильно ослабляются разные частотные составляющие. Другими словами, необходимо построить амплитудно-частотную характеристику соответствующего фильтра, иначе велика вероятность «нарваться» на артефакты.
Задача сглаживания может использоваться при прореживании сигналов, то есть когда, например, необходимо отобразить большую картинку на небольшой экран. Или когда частота дискретизации звука снижается, например, с 48000 Гц до 44100 Гц. Понижение частоты выборок — коварная операция, требующая предварительной обработки сигнала (низкочастотной фильтрации), но это — тема отдельного разговора…
Приведём пример «плохого» сглаживания
![](https://habrastorage.org/storage2/caa/4a9/045/caa4a90451517ce0cbe7f528645e5911.gif)
Казалось бы, обычное усреднение и сигнал на выходе должен быть «гладким». Но как определить, насколько он стал «глаже»? Не переборщили ли мы? А может быть некоторые коэффициенты выбрать не по 1/3? А может быть усреднить по пяти точкам? Как определить насколько ослабляются частотные составляющие в сигнале? Как найти свой (то есть для конкретной задачи) оптимум?
На эти и некоторые другие вопросы я постараюсь ответить так, чтобы «обычный» программист смог обосновать свой алгоритм, — надеюсь, не только алгоритм на тему «Сглаживание», так как идеи будут излагаться весьма общие, заставляющие думать самому…
Сигналы
В этом параграфе под сигналом понимается вектор, то есть массив из определённого количества чисел. Например, вектор из четырёх элементов s = (2,5; 5; 0; -5).
Для простоты будут рассматриваться только десятичные вещественные числа.
Одним из наиболее распространённых и понятных сигналов является оцифрованный звук. Размер сигнала зависит от длительности звука и от частоты, с которой делают выборки (от частоты дискретизации). Элементы-числа сигнала зависят от текущей амплитуды звука, измеряемой устройством выборки и хранения.
Как уже было сказано, один из самых простых способов сглаживания, это
![](https://habrastorage.org/storage2/c4b/402/3ca/c4b4023caa728489dc36a9d34b708c7f.png)
где s — исходный сигнал, v — сглаженный сигнал.
Способ (1) основан на сглаживающем свойстве суммирования, ведь каждому ясно, что средняя величина, вычисляемая как сумма многих случайных чисел, с ростом количества суммируемых чисел становится всё менее и менее похожей на случайную величину 2, которая, попросту говоря, и есть шум 3.
Но на чём основан выбор коэффициентов в уравнении (1)? На том, что так вычисляется среднее? Вроде бы да, но… А если взять не три слагаемых, а шестнадцать? А тридцать два?.. Почему всё более отстоящие от центрального элемента s[i] отсчёты должны браться с одинаковым весом? Ведь может оказаться так, что связь между отсчётами будет постепенно ослабевать с ростом расстояния 4 между ними?
Если рассмотреть пример произношения слова «арбуз» десять раз подряд и попытаться отследить связь между отсчётами записанного сигнала, то можно обнаружить ослабление зависимости между всё более отстоящими друг от друга отсчётами. Естественно, что если рассмотреть «большие расстояния», то звуки будут повторяться за счёт повтора одного и того же слова и зависимость будет нарастать и снова спадать, и так далее. Но, как правило, «большие расстояния» при сглаживании не рассматривают, так как шумы проявляются в окрестности отдельных звуков, а не слов, фраз и предложений. Шум, действующий на уровне слов или даже фраз — это явно искусственный (звуковые эффекты) или экзотический естественный (эхо). Это — уже «неслучайный» шум, требующий отдельного исследования. Здесь рассматривается «чистый» шум, который, говоря простым языком, раздражающе шумит и нисколько не похож на какой-либо полезный сигнал.
На основании простых рассуждений становится очевидным, что количество слагаемых в (1) (порядок фильтра) должно зависеть от того, насколько сильно зависят друг от друга соседние отсчёты. Например, нет смысла брать фильтр тридцатого порядка, если наблюдается зависимость только лишь десяти отсчётов, следующих друг за другом. На самом деле даже не то, что «нет смысла», а — нельзя, так как если отсчёты практически не связаны, то начнётся чрезмерное сглаживание полезного сигнала («съедение» слогов). Но и фильтр третьего порядка здесь не будет оптимальным по степени использования информации о полезном сигнале, так как, как уже было сказано, наблюдается зависимость порядка десяти соседних отсчётов. Поэтому можно «попытать счастья» с помощью фильтра девятого порядка, естественно, увеличив нагрузку на процессор-вычислитель. Здесь уже требуется определить, скорее всего экспериментально, а стоит ли данная игра свеч?..
Как оценить насколько сильно связаны соседние отсчёты? Вычислить автокорреляционную функцию (АКФ). Желающим можно предложить провести эксперимент по записи разных слов, фраз, повторов фраз и последующему построению АКФ (благо, например, программа Matlab позволяет это сделать, особо не задумываясь над кодом и формулами).
Так как всё-таки выбрать коэффициенты фильтра в (1)?
В данном случае удобно рассмотреть реакцию фильтра на единичное воздействие, то есть на сигнал вида
![](https://habrastorage.org/storage2/88e/f51/126/88ef51126cc78cab8ba4727d69b3505f.png)
Например, фильтр (1) даст следующий отклик (импульсную характеристику)
![](https://habrastorage.org/storage2/824/438/79e/82443879ed0ac6de39a0fe33deeadf18.png)
откуда мы можем сделать вывод, что после сглаживания длительность сигнала стала равна трём элементам. А если взять фильтр из пяти элементов?.. Правильно, длительность выходного сигнала будет равна пяти элементам. Насколько это полезно, определяется конкретной ситуацией (задачей).
Кстати, а долгожданный артефакт уже налицо! Импульсная характеристика (1) — это, по сути, прямоугольный импульс, нисколько не являющийся гладким!.. Странно, да? А если взять пятиточечный фильтр? Тогда на выходе получим более длинный прямоугольный импульс, но с меньшей амплитудой. Не очень хорошо выходит… Простейший тест говорит о непригодности простого усреднения для сглаживания.
Посмотрим на фильтр (1) с частотной стороны (с временной мы уже посмотрели).
Если сигнал звуковой, то он достаточно хорошо описывается набором гармонических сигналов 5 («синусоид») и степень ослабления конкретной синусоиды зависит от её частоты. Опять же, при грамотном выборе сглаживающего фильтра никакая из полезных синусоид не должна пропадать полностью, то есть амплитудно-частотная характеристика фильтра в рабочем диапазоне частот должна быть достаточно равномерной.
Пропустим через рассматриваемый фильтр один однотональный сигнал определённой частоты, естественно, не выходя за пределы теоремы отсчётов. Пусть шаг дискретизации по времени Td равен единице, то есть отсчёты идут с шагом одна секунда. Возьмём сигнал с частотой f = 1/(3Td) = 1/3 Гц, то есть
![](https://habrastorage.org/storage2/94e/524/0bd/94e5240bd1b6ad6b6fcbb2daa977b00c.gif)
Ограничимся двумя периодами
![](https://habrastorage.org/storage2/92f/471/cd0/92f471cd0801ded8adade0cb55494e57.gif)
Отклик фильтра (1) будет равен
![](https://habrastorage.org/storage2/cf5/c38/9dd/cf5c389dd6fed39b0a924a3d9a6113a2.png)
Как ни странно, получили почти ноль… Получается, можно потерять некоторые составляющие полезного сигнала.
Проверим отклик на несколько более высокочастотный сигнал
![](https://habrastorage.org/storage2/fca/f90/43f/fcaf9043fec33d502371eaf5ee360c69.png)
![](https://habrastorage.org/storage2/296/53f/8d8/29653f8d89c61969e0abac25bf5a125a.png)
Как видим, форму сигнала не потеряли. В чём же дело?..
Дело в том, что амплитудно-частотная характеристика фильтра (1) не является монотонной в рабочей полосе частот (в полосе от нуля до половины частоты дискретизации) и имеет один нуль на частоте, в три раза меньшей частоты дискретизации. Как это показать?
Попросту говоря, чтобы определить частотную характеристику фильтра, необходимо найти отношение спектра на выходе фильтра к спектру на входе.
Обозначим спектр сигнала s[i] на входе как S( f ), тогда спектр задержанного на один такт Td сигнала s[i-1] выразится через спектр исходного сигнала как
![](https://habrastorage.org/storage2/c72/900/f0f/c72900f0f984caa58aee8cf7a5ab9fdd.gif)
где j — мнимая единица.
Спектр опережённого сигнала s[i+1] выразится как
![](https://habrastorage.org/storage2/b26/c87/fe4/b26c87fe45058d633a8611dad7f3d7b5.gif)
Что означает мнимая единица? И как можно обосновать (2) и (3)?
Если записать известные [1, 2] ряды для синуса, косинуса и экспоненты
![](https://habrastorage.org/storage2/5ce/a38/fd7/5cea38fd765c42afc282a90fb818c2ec.gif)
то подобрав число j, такое, что j2 = –1, можно последний ряд выразить через два первых
![](https://habrastorage.org/storage2/588/871/d9a/588871d9a2f37485307174bc013143ea.gif)
что означает то, что любое комплексное число можно записать через экспоненту с мнимым показателем. Модуль числа (5) равен единице (корень квадратный из суммы квадратов мнимой и реальной частей), поэтому для записи любого комплексного числа в форме (5) его необходимо разделить и умножить на свой модуль
![](https://habrastorage.org/storage2/ec3/e80/f78/ec3e80f78affdde6e4832e73f5cf9a99.gif)
Из (5) и (6) следует, что если в показателе экспоненты можно выделить мнимую единицу, умноженную на некоторое действительное число, то это число является аргументом комплексного числа.
В данном случае рассматриваются сигналы, поэтому модулю комплексного числа соответствует амплитуда гармонического сигнала, а аргументу — фаза. Если, например, взять сигнал вида
![](https://habrastorage.org/storage2/8a5/c8b/13a/8a5c8b13a5241de6c0a727c3982673b8.gif)
то можно выделить амплитуду A и фазу Ф. Множитель
![](https://habrastorage.org/storage2/7f1/6a6/8f3/7f16a68f3e3df43cb22149dfd87e80c8.gif)
Если сигнал (7) задержать на величину t0, то получится тот же самый сигнал, но сдвинутый по фазе на постоянную величину
![](https://habrastorage.org/storage2/34a/6ee/ffa/34a6eeffa7b41dfac96cceffe9d3433c.gif)
то есть при задержке произвольного сигнала все его частотные составляющие сдвигаются по фазе на величину, зависящую от текущей частоты и величины задержки. Этим можно объяснить формулы (2) и (3). Поэтому при анализе какого-либо алгоритма важна и фазо-частотная характеристика, которая показывает на какое время задерживает фильтр (алгоритм) каждую частотную составляющую входного сигнала. Низкие частоты и высокие в общем случае будут иметь разную задержку в фильтре.
Из (1), (2) и (3) следует, что частотная характеристика фильтра (передаточная функция) будет иметь вид
![](https://habrastorage.org/storage2/448/8ea/daf/4488eadaf2b1f4195f89515a0de28359.gif)
Так как спектр выходного сигнала линейно выражается через спектр входного, то при выводе (9) спектр входного сигнала успешно сокращается. Далее замечаем, что частотная характеристика фильтра (1) — действительная, то есть данный фильтр не вносит фазовых искажений 6. Этого мы достигли (скорее всего, неосознанно) за счёт симметричности формулы (1): каждый отсчёт на выходе фильтра равен сумме текущего и двух соседних отсчётов.
Физически такой алгоритм реализуем только при наличии запоминающих устройств, так как в нём используются опережающие отсчёты (для вычисления отсчёта s[i] требуется отсчёт s[i+1]). В настоящее время это не является большой проблемой и, как правило, используют симметричные алгоритмы. Если фазовые искажения окажутся полезными, то — пожалуйста, главное осознанно применять какой-либо алгоритм и смотреть на него с разных сторон: с частотной и временной.
Построить график зависимости (9) от частоты не составляет труда. Для упрощения вводят нормированную частоту f0 = f Td , полезный диапазон изменения которой [0…0,5]. Составляющие сигнала с частотами выше половины частоты дискретизации по теореме отсчётов должны отсутствовать (перед оцифровкой сигнал пропускают через соответствующие фильтры нижних частот). Частота дискретизации показывает количество выдаваемых цифровым устройством отсчётов в секунду. Если, например, один такт Td равен одной миллисекунде, то за одну секунду должна быть выдана тысяча отсчётов.
Анализируя (9) можно заметить, что на некотором промежутке коэффициент передачи меньше нуля, а амплитуда — это по определению положительная величина… Выход из ситуации — построить модуль передаточной функции, убрав знак минус в фазовую характеристику, которая, всё-таки, не является константой (нулём). Если взять число «минус единица», то его по формуле (5) можно представить как
![](https://habrastorage.org/storage2/2ac/051/2f4/2ac0512f4004e122f4f2595cf82ca70e.gif)
то есть комплексным числом, имеющим модуль «единица» и фазу «180 градусов» («пи»).
Таким образом, трёхточечный симметричный алгоритм (1) для некоторых «высоких» частот вносит сдвиг по фазе на 180 градусов, то есть попросту инвертирует входной сигнал. Этот эффект можно заметить, анализируя рассмотренный выше отклик на нормированную частоту 2/5 Гц.
![](https://habrastorage.org/storage2/778/7c0/323/7787c0323449fe552509dd7d0d0c3cf3.png)
Рис. 1. Амплитудно-частотная и фазо-частотная характеристики трёхточечного симметричного алгоритма сглаживания (1)
Из рис. 1 сразу следует, что сигнал с частотой 1/3 будет данным алгоритмом подавлен, а сигналы, имеющие частоту выше 1/3 будут инвертированы. Таким образом, полезный рабочий диапазон частот можно смело сократить с [0…0,5] до [0…1/3]. Если нас не устраивает быстрое спадание коэффициента передачи, придётся определять другой алгоритм, имеющий более прямоугольную амплитудно-частотную характеристику и при этом — ещё неизвестно какую фазовую…
По сути, полученная немонотонная частотная характеристика — следствие прямоугольной формы импульсной характеристики …, 0, 0, …, 1/3, 1/3, 1/3, 0, 0, …. Поэтому импульсная характеристика — простой способ заметить грубые изъяны в алгоритмах. Частотная характеристика, несмотря на некоторую сложность её вычисления, удобна тем, что позволяет заметить и устранить более тонкие изъяны, чем мы и займёмся.
Если теперь записать алгоритм (1) в более общем виде
![](https://habrastorage.org/storage2/e5a/507/26a/e5a50726a975e497a0c3e4b637bd696e.gif)
то, основываясь на знании частотной характеристики, можно попытаться подобрать коэффициенты так,
чтобы амплитудно-частотная характеристика стала монотонной в рабочем диапазоне частот (0…0,5). Для этого, как минимум, необходимо отсутствие нулей внутри рабочего диапазона.
Так как у нас нет оснований выделять запаздывающий отсчёт s[i–1] по отношению к опережающему s[i+1], то приравняем коэффициенты a1 и a3. После запишем коэффициент передачи
![](https://habrastorage.org/storage2/805/06f/29f/80506f29fe900602d81d0fe7a07e442b.gif)
Попытаемся переместить первый нуль на частоту f0 = 0,5. Для этого должно выполняться равенство a2 = 2a1 , то есть вес у боковых отсчётов должен быть в два раза меньше веса центрального. Тогда более оптимальный алгоритм будет выглядеть так
![](https://habrastorage.org/storage2/160/3c6/627/1603c6627f416752148e8acfd92da0b7.gif)
Как в (13) выбрать единственный коэффициент a1?
Взглянем на алгоритм (13) с точки зрения импульсной переходной характеристики. Для этого найдём отклик на единичный скачок
![](https://habrastorage.org/storage2/7d6/4b1/21a/7d64b121a46913d243582478509976ad.gif)
![](https://habrastorage.org/storage2/3ab/e06/6ac/3abe066ac89ce37c2a775399a1607b92.gif)
Как видим, чтобы сохранить выходную амплитуду в установившемся режиме, равную единице, требуется выбрать коэффициент a1=1/4. То есть сумма всех коэффициентов должна быть равна единице.
Наконец, для готового алгоритма
![](https://habrastorage.org/storage2/827/0ce/f3b/8270cef3b56b2e5b2e76ce6df03146f9.gif)
можно построить (рис. 2) частотные характеристики: амплитудную и фазовую.
![](https://habrastorage.org/storage2/50f/4cf/7d9/50f4cf7d938b34b2b710063067c97a0e.png)
Рис. 2. Амплитудно-частотная и фазо-частотная характеристики трёхточечного симметричного алгоритма сглаживания (16)
Анализ рис. 2 показывает, что фазовые искажения исчезли, а амплитудная характеристика стала монотонной в рабочей полосе частот [0…0,5]. В каком-то смысле мы выжали из трёхточечного фильтра «всё».
Теперь становится очевидно, что простое усреднение далеко не всегда является оптимальным, особенно когда усредняемых отсчётов много.
Изображения
Что касается изображений, то здесь в какой-то мере всё то же самое, что и для сигналов: амплитудные и фазовые искажения, импульсные характеристики, энергия сигнала, — основное отличие заключается в том, что используются матрицы вместо векторов. Также следует учесть, что для изображений нет координаты времени, а есть координата пространства, поэтому есть пространственные частоты дискретизации. Данное отличие скорее формальное, так как на математику алгоритмов оно не влияет вообще.
Рассмотрим теперь «самый простой» алгоритм сглаживания изображения по соседним точкам (рис. 3). Отсчёт v00 на выходе фильтра
![](https://habrastorage.org/storage2/941/e54/e65/941e54e65d514db52432308ac19b69c6.gif)
![](https://habrastorage.org/storage2/1d0/072/176/1d007217671352f6f974eb404cf1ece0.png)
Рис. 3. Схема сглаживания изображения по соседним точкам
В формуле (17) специально выделены три слагаемых A, B и C, так как четыре соответствующих внутренних слагаемых у B и C имеют свои расстояния от центрального отсчёта s00. Здесь естественно предположить, что максимальный вес будет у слагаемого A, затем в порядке убывания — у B и C.
Для изображений спектр и коэффициент передачи будут двумерными, то есть будут зависеть от двух частот: первая частота — по горизонтали, вторая — по вертикали (здесь всё условно, для определённости).
Фильтр (17) имеет следующий коэффициент передачи
![](https://habrastorage.org/storage2/433/1b6/c0c/4331b6c0c8161149f10f48e7dff5c98a.gif)
Если взять, например, около тридцати дискретных частот и построить (рис. 4) контурный график модуля частотной характеристики (18), то можно увидеть искажения, аналогичные искажениям на рис. 1. Левому нижнему углу соответствуют самые низкие частоты. Наблюдается провал частотной характеристики на частотах, составляющих 1/3 от частоты дискретизации (частоты дискретизации в данном случае равны 64 по обеим координатам).
![](https://habrastorage.org/storage2/2d6/438/a6f/2d6438a6f6c578919e66bdfc60b2e88a.png)
Рис. 4. Амплитудно-частотная характеристика девятиточечного симметричного алгоритма сглаживания (17)
Перетасовывая коэффициенты для отсчётов на рис. 3, можно получить следующий алгоритм сглаживания
![](https://habrastorage.org/storage2/e76/55e/56a/e7655e56a766e11ffd7c8cd7ed14d45d.gif)
Причём если единичное воздействие поместить в центре изображения, то отклик фильтра (19) будет иметь вид
![](https://habrastorage.org/storage2/8eb/bd9/bc4/8ebbd9bc4f7196286a449c30751d0d81.gif)
По сути, (20) является импульсной характеристикой фильтра (это девять главных отсчётов, так как остальные равны нулю). Сумма всех отсчётов импульсной характеристики равна единице. Частотная характеристика фильтра (19) имеет вид, показанный на рис. 5.
![](https://habrastorage.org/storage2/270/5e9/e07/2705e9e07bed10bf92403c5f0f63429c.png)
Рис. 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