Как стать автором
Обновить

Комментарии 38

Левая картинка смотрится лучше

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

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

оригинальная картинка
оригинальная картинка

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

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

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

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

Фильтр Калмана
Еще один фильтр, принцип работы которого не столь очевиден для каждого. Из всего кода вам необходимо знать, что r — это приблизительная амплитуда шума, а q следует подбирать самостоятельно.


Хотелось бы понять, что такое амплитуда шума (или шире — амплитуда случайного процесса), к тому же еще и приблизительная? И по каким принципам следует подбирать q?

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

Спасибо, я учту это и доработаю статью

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

  • из 5 значений находим максимальное и минимальное;

  • вычитаем их из общей суммы и результат делим на 3.

Интересно, у этого фильтра есть название? Не могу нагуглить ничего похожего.

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

Интересно, у этого фильтра есть название? Не могу нагуглить ничего похожего.


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

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

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

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

Кстати, кроме прямоугольного и экспоненциального окна, рассмотренных в статье, существует огромное множество других сглаживающих ядер. Лично мне больше нравится гауссиан - он довольно эффективно эмулирует частотную фильтрацию для сигналов с пропусками (а у нас они есть всегда). Дело в том, что обычные частотные фильтры имеют знакопеременные фильтрующие (сверточные) последовательности, и при наличии пропусков это часто приводит к проблемам. Если же пересчитывать фильтр для каждого нового набора пропущенных наблюдений в окне, то это ужасно долго. (Заранее рассчитать все варианты фильтра можно только для очень маленьких окон - число комбинаций растет, как 2^N).

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

Зачем же писать свой moving average? В Матлабе уже есть movmean().

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

Спасибо за ваш комментарий, благодаря нему прочитал вашу статью про сглаживание. Подскажите, смотреть на АЧХ фильтра есть смысл ведь только если известна АЧХ шума? Но как правильно усреднять(уточнять) результаты измерения, если измеряемая величина постоянна, а измерений произведено несколько в случайные моменты времени?

Видимо да, т.е. если известна СПМ шума. Последний вопрос не понял: если величина константа, то зачем ее несколько раз измерять?

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

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

Тогда именно СПМ шума важна, когда мы измеряем в пространственной точке. Это фильтр Винера.

Квадрат АЧХ фильтра, G(f), должен стремиться к обратной СПМ шума: G(f) ~ 1 / N(f)

и, скорее всего, калмановский фильтр сойдется к винеровскому при правильной его настройке, т.е. при учете именно СПМ шума (т.е. динамической модели системы, генерирующей шум).

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

Фильтр Калмана, в свою очередь, допускает нестационарность процессов.

Еще зависит от цели: отфильтровать шум, при этом не исказив сильно сигнал, т.е. зависит от баланса между степенью фильтрации шума и степенью искажения сигнала. По факту, АЧХ фильтра всё равно важна, равно как и ФЧХ.

Как то странно реализован фильтр "Калмана" у Вас и тот что по ссылке. Совсем не похож на тот что приведен в вики.

Довольно очевидная оптимизация фильтра Калмана:

def kalman(f, q=0.25, r=0.7):
    if not hasattr(kalman, "Accumulated_Error"):
        kalman.Accumulated_Error_Sq = 1
        kalman.kalman_adc_old = 0

    if abs(f-kalman.kalman_adc_old) / 50 > 0.25:  # здесь необходимо пояснение констант нот
        Old_Input = f*0.382 + kalman.kalman_adc_old*0.618  # возможно, здесь тоже
    else:
        Old_Input = kalman.kalman_adc_old

    Old_Error_Sq = kalman.Accumulated_Error_Sq + q**2
    H = Old_Error_Sq / (Old_Error_Sq + r**2)
    kalman_adc = Old_Input + H * (f - Old_Input)
    kalman.Accumulated_Error_Sq = (1 - H) * Old_Error_Sq
    kalman.kalman_adc_old = kalman_adc

    return kalman_adc

Можно ещё оптимизировать, вынеся квадраты q и r в переменные доступные извне и вызывая kalman(f) с одним параметром.

И тогда тяжёлых вычислений не останется

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

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


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

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

Если я правильно понял, фильтр Калмана использует только статистическую модель сигнала, без учета динамики (формы кривой)?
Вообще, откуда он берет модель сигнала? Сам считает постепенно, или ему её надо в готовом виде подать?

Вообще, откуда он берет модель сигнала? Сам считает постепенно, или ему её надо в готовом виде подать?

Модель системы задается матрицами F, B, Q. Дальше основываясь на этой модели фильтр на каждом шаге рассчитывает ковариационную матрицу P и прогноз текущего состояния объекта x.

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

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

Получается, модель в «состав фильтра» не входит, её пользователь сам должен как-то определить?

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

И тогда не совсем понятно, если у нас уже есть модель (например «сигнал — это синусоида»)

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

Вобще для чистого наблюдения с течением времени ФК вырождается в ФНЧ 1-го порядка:

y = y + (x-y)*k,

где: y - выход ФК, x - вход, k = Q/(R+Q). Q,R - константы (параметры модели).

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

Наилучшая по какому критерию? Для каких моделей?

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

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

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

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

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

Также надо знать, допустима ли задержка, или требуется отклик в реальном времени. Если задержка допустима, то лучше использовать симметричные фильтры, чем "скользящее среднее", которое вообще пора запретить в реальных приложениях :)

Симметричные фильтры для задач подавления широкополосного шума мало подходят. Намного лучше подходит использовать в качестве ядра свёртки непосредственно оконные функции, как классические, так и новые. А ещё лучше использовать МНК с тригонометрическим базисом и перекрытием в 2 или 4 раза.

А про запрет скользящего среднего полностью поддерживаю.

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

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

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

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

Суть статьи: есть такие фильтры, их можно иногда использовать, а ещё бывают помехи разные. можно попробовать каждый фильтр по отдельности или совместно... где результат будет похож на "правду", тот и используйте.

Треугольный сигнал? Может всё же пилообразный?

Пиловидный сигнал по сути подвид треугольного :)

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

Сигнал треугольной формы, почему нет

Зарегистрируйтесь на Хабре, чтобы оставить комментарий

Публикации

Истории