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

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

Введение

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

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

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

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

Поэтому для некоторого алгоритма сглаживания в первую очередь необходимо определить как сильно ослабляются разные частотные составляющие. Другими словами, необходимо построить амплитудно-частотную характеристику соответствующего фильтра, иначе велика вероятность «нарваться» на артефакты.

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

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

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

свой (то есть для конкретной задачи) оптимум?
На эти и некоторые другие вопросы я постараюсь ответить так, чтобы «обычный» программист смог обосновать свой алгоритм, — надеюсь, не только алгоритм на тему «Сглаживание», так как идеи будут излагаться весьма общие, заставляющие думать самому…

Сигналы

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

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

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

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

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

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

4 между ними?

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

f.

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

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

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

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

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

Анализируя (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], то приравняем коэффициенты a₁ и a₃. После запишем коэффициент передачи
(12)

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

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

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

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

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

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

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

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

Изображения

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

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

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

В формуле (17) специально выделены три слагаемых A, B и C, так как четыре соответствующих внутренних слагаемых у B и C имеют свои расстояния от центрального отсчёта s₀₀. Здесь естественно предположить, что максимальный вес будет у слагаемого 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

Сглаживание

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

1) Цифровая фильтрация — чаще всего первого порядка или оптимальная фильтрация;

2) Сглаживание измеренных значений за счет формирования средних -метод скользящего среднего или прогрессивная интерполяция;

3) Фильтрация сигналов с помощью фильтра нижних частот (ФНЧ)).

Цифровая фильтрация первого порядка

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

где — фильтруемая входная величина,

— отфильтрованная выходная величина.

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

где — значение выборки сигнала в момент времени .

После решения относительно отфильтрованной выходной величины:

где — интервал дискретизации сигнала.

Такую рекурсивную форму можно использовать в качестве алгоритма работы ЭВМ:

где — выбираемый фактор сглаживания, 0<А<1.

Таким образом, работа фильтра заключается в получении отфильтрованных значений сигнала из текущего зашумленного с учетом предыдущего значения. А — постоянный коэффициент, принимающий значения 1>A³0,5 для более динамичных сигналов, передаваемых в условиях слабых помех, 0,5>A>1 -наоборот, для медленно изменяющихся сигналов на фоне сильных помех.

Оптимальная фильтрация

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

где — весовая показательная функция, n — весовой показатель.

Такую рекурсивную формулу удобно применять в качестве алгоритма на ЭВМ. В качестве фактора сглаживания выбирается весовой показатель n. Очевидно, что при n=1 получаем среднее арифметическое значение предыдущего отфильтрованного значения сигнала и текущего зашумленного. Именно это значение берется в качестве текущего отфильтрованного значения сигнала. Таким образом, при n>1 будем получать соответствующие моменты более высокого порядка, причем при увеличении n влияние предыдущего отфильтрованного значения сигнала каждый раз значительно увеличивается.

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

Прогрессивная интерполяция (метод усреднения)

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

Отдельные значения S суммируются. После выполнения суммирования вычисляется среднее значение путем деления суммы на число N. Шум, в данном случае, рассматривается в виде стационарного некоррелируемого процесса и характеризуется значением математического ожидания E(n), которое при достаточном числе проб стремится к нулю E(n)=0. Среднее арифметическое из суммы выборок сигнала с наложенным на него шумом будет равно:,

где N — количество выборок в группе,

i — номер отсчета сигнала,

k — текущая выборка.

N должно выбираться не более необходимого для нахождения отсчета сигнала. Можно записать: .

Если частота выборок сигнала превышает высшую частоту спектра сигнала , то , откуда .

Поскольку среднее значение шума равно нулю, то значит в пределе при достаточно большом N мы должны получить в результате суммирования значение сигнала без примеси шума. Так как E(n)=0, доля примеси шума зависит от величины его дисперсии. Учитывая, что дисперсия случайной величины x равна запишем:

Если считать, что выборки не коррелируют друг с другом, можно окончательно записать, т.к. :

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

.

Отсюда следует, что отношение сигнал/шум должно расти при усреднении пропорционально значению .

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

Применение данного метода сглаживания помех имеет и ряд существенных ограничений:

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

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

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

Фильтрация сигналов с помощью фильтра нижних частот (ФНЧ)

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

Оптимальное значение можно найти из соотношения: .

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

Учитывая что ; ;

,тогда

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

python — Как правильно сгладить кривую?

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

 | время работы в с | время работы в с
метод | список питонов | пустой массив
------|---------------|------------
регрессия ядра | 23.93405 | 22.75967
низкий | 0,61351 | 0,61524
наивный средний | 0,02485 | 0,02326
другие* | 0,00150 | 0,00150
фт | 0,00021 | 0,00021
свернуть | 0,00017 | 0,00015
* savgol с различными функциями подгонки и некоторыми методами numpy
 

Регрессия ядра плохо масштабируется, Lowess немного быстрее, но оба дают плавные кривые. Savgol — это золотая середина по скорости и может давать как скачкообразные, так и плавные выходные данные, в зависимости от класса полинома. БПФ очень быстрый, но работает только с периодическими данными.

Методы скользящего среднего с numpy быстрее, но, очевидно, создают график с шагами.

Я создал 1000 точек данных в форме кривой синусоиды:

 размер = 1000
х = np.linspace (0, 4 * np.pi, размер)
y = np.sin(x) + np.random.random(размер) * 0,2
данные = {"х": х, "у": у}
 

Я передаю их в функцию для измерения времени выполнения и построения графика соответствия:

 def test_func(f, label): # f: дескриптор функции для одного из методов сглаживания
    начало = время ()
    для я в диапазоне (5):
        обр = f (данные ["у"], 20)
    print(f"{метка:26с} - время: {время() - начало:8.5f} ")
    plt.plot (данные ["x"], обр, "-", метка = метка)
 

Я протестировал множество различных функций сглаживания. arr — это массив значений y для сглаживания, а span — параметр сглаживания. Чем ниже, тем лучше аппроксимация будет приближаться к исходным данным, чем выше, тем более гладкой будет результирующая кривая.

 def smooth_data_convolve_my_average (пример, диапазон):
    re = np.convolve(arr, np.ones(span * 2 + 1) / (span * 2 + 1), mode="same")
    # Часть "my_average": сужает окно усреднения со стороны
    # выходит за пределы данных, сохраняя другую сторону того же размера, что и заданный
    # по "промежутку"
    re[0] = np.average(обр[:диапазон])
    для i в диапазоне (1, интервал + 1):
        re[i] = np.average(arr[:i + span])
        re[-i] = np.average(arr[-i - span:])
    вернуться повторно
def smooth_data_np_average(arr, span): # мой оригинальный, наивный подход
    вернуть [np.average (arr [val - span: val + span + 1]) для val в диапазоне (len (arr))]
def smooth_data_np_convolve (пример, диапазон):
    вернуть np.convolve(arr, np.ones(span * 2 + 1)/(span * 2 + 1), mode="same")
def smooth_data_np_cumsum_my_average (пример, диапазон):
    cumsum_vec = np. cumsum (обр)
    moving_average = (cumsum_vec[2 * span:] - cumsum_vec[:-2 * span]) / (2 * span)
    # Снова часть "my_average". Немного отличается от предыдущего, потому что
    # скользящее среднее из cumsum короче входного и должно быть дополнено
    спереди, сзади = [np.average(arr[:span])], []
    для i в диапазоне (1, диапазон):
        front.append(np.average(arr[:i + span]))
        back.insert(0, np.average(arr[-i - span:]))
    back.insert (0, np.average (обр [-2 * диапазон:]))
    вернуть np.concatenate ((спереди, moving_average, сзади))
def smooth_data_lowess (обр, диапазон):
    х = np.linspace (0, 1, длина (обр))
    return sm.nonparametric.lowess(arr, x, frac=(5*span / len(arr)), return_sorted=False)
def smooth_data_kernel_regression (обозначение, диапазон):
    # Параметр сглаживания "span" игнорируется. Если вы знаете, как
    # включите это в регрессию ядра, пожалуйста, прокомментируйте ниже.
    kr = KernelReg (arr, np.linspace (0, 1, len (arr)), 'c')
    вернуть kr.fit () [0]
def smooth_data_savgol_0 (обр, диапазон):
    вернуть savgol_filter (обр, диапазон * 2 + 1, 0)
def smooth_data_savgol_1 (обр, диапазон):
    вернуть savgol_filter (обр, диапазон * 2 + 1, 1)
def smooth_data_savgol_2 (обр, диапазон):
    вернуть savgol_filter (обр, диапазон * 2 + 1, 2)
def smooth_data_fft(arr, span): # масштабирование "span" открыто для предложений
    w = fftpack. rfft(обр)
    спектр = ш ** 2
    cutoff_idx = спектр < (spectrum.max() * (1 - np.exp(-span/2000)))
    w[cutoff_idx] = 0
    вернуть fftpack.irfft(w)
 

Скорость

Время выполнения более 1000 элементов, протестировано в списке Python, а также в массиве numpy для хранения значений.

 метод | список питонов | пустой массив
------|--------------|------------
регрессия ядра | 23,93405 с | 22,75967 с
низкий | 0,61351 с | 0,61524 с
среднее значение | 0,02485 с | 0,02326 с
савгол 2 | 0,00186 с | 0,00196 с
савгол 1 | 0,00157 с | 0,00161 с
савгол 0 | 0,00155 с | 0,00151 с
numpy свернуться + я | 0,00121 с | 0,00115 с
numpy cumsum + я | 0,00114 с | 0,00105 с
фт | 0,00021 с | 0,00021 с
свернуть | 0,00017 с | 0,00015 с
 

Особенно регрессия ядра очень медленно вычисляет более 1k элементов, lowess также дает сбой, когда набор данных становится намного больше. numpy convolve и fft особенно быстры. Я не исследовал поведение во время выполнения (O(n)) при увеличении или уменьшении размера выборки.

Поведение края

Я разделю эту часть на две части, чтобы изображение было понятным.

Методы на основе Numpy + savgol 0 :

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

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

• savgol 0 : продолжается константой до края данных ( savgol 1 и savgol 2 заканчиваются линией и параболой соответственно)
• numpy Average : останавливается, когда окно достигает левой части данных и заполняет эти места в массиве Nan , то же поведение, что и метод my_average с правой стороны
• numpy convolve : довольно точно следует данным. Я подозреваю, что размер окна уменьшается симметрично, когда одна сторона окна достигает края данных
• my_average / me : мой собственный метод, который я реализовал, потому что другие меня не устраивали. Просто сжимает часть окна, выходящую за пределы данных, к краю данных, но сохраняет исходный размер окна с другой стороны, заданный с помощью span
.

Сложные методы:

Все эти методы хорошо подходят к данным. савгол 1 заканчивается линией, савгол 2 параболой.

Поведение кривой

Для демонстрации поведения различных методов в середине данных.

Различные фильтры savgol и medium дают грубую линию, lowess , fft и регрессия ядра дают гладкую аппроксимацию. lowess срезает углы при изменении данных.

У меня есть данные регистрации Raspberry Pi для развлечения, и визуализация оказалась небольшой проблемой. Все точки данных, за исключением использования ОЗУ и трафика Ethernet, записываются только дискретными шагами и/или зашумлены по своей природе. Например, датчик температуры выдает только целые градусы, но отличается до двух градусов между последовательными измерениями. Из такого графика рассеяния нельзя получить никакой полезной информации. Поэтому для визуализации данных мне понадобился какой-то метод, не слишком затратный в вычислительном отношении и производящий скользящее среднее. Я также хотел получить хорошее поведение на краях данных, так как это особенно влияет на последнюю информацию при просмотре оперативных данных. Я остановился на 9Метод 0019 numpy convolve с my_average для улучшения поведения края.

фильтр нижних частот - Как сгладить сигнал?

спросил 10 лет, 11 месяцев назад

Изменено 10 лет, 11 месяцев назад

Просмотрено 5к раз

$\begingroup$

Самый простой способ сгладить сигнал — это скользящее среднее значение окна.

Более продвинутый способ — использовать фильтр Савицкого-Голея. Из википедии:

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

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

Сглаживание и спектральная утечка — одно и то же?

Когда следует использовать Савицкого-Голея и когда следует использовать Гаусса, Ханна, Хэмминга и т.д.?

Заранее спасибо за любые ответы!

• Фильтр нижних частот

$\endgroup$

3

$\begingroup$

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

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

Адаптивная полиномиальная регрессия удалит различные спектральные характеристики сигнала в разные моменты времени, что может либо выявить, либо скрыть «важные» элементы сигнала, в зависимости от сигнала и того, что вы ищете.

$\endgroup$

2

$\begingroup$

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