Фильтр Савицкого-Голея и иже с ним

← →
MBo © (2009-12-02 17:29) [0]

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

← →
Jeer © (2009-12-02 18:25) [1]

В соответствии с Ван-Трисом все задачи фильтрации подразделяются на 3 группы:
1. известный сигнал и шумы
2. известный сигнал и неизв.шумы
3. неизвестные сигнал и шумы

Оптимальный фильтр для коррелированных шумов по гр.1 имеет ЧХ вида
H(f) = S(f) / [ S(f) + N(f) ]
S(f) - СП полезного случ.сигнала
N(f) - СП коррелированного шума
Если что-то неизвестно, то опт.фильтр невозможен.
При этом в некоррелированных шумах используют упомянутые полиномиальные нерекурсивные фильтры S-G. Теоретическая применимость их для коррелированных шумов неясна.

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

← →
Jeer © (2009-12-02 18:30) [2]

P.S.
Ну а вообще, мне трудно представить, чтобы Mbo что-то не знал :)

← →
Германн © (2009-12-03 00:44) [3]

> мне трудно представить, чтобы Mbo что-то не знал :)

+1
На моей памяти это первый его серьезный вопрос на ДМ.
:)

← →
MBo © (2009-12-03 07:31) [4]

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

← →
Jeer © (2009-12-03 09:28) [5]

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

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

Ну это так, в двух словах. Главное тут - степень свободы значительно выше.

← →
Jeer © (2009-12-04 16:20) [6]

2Mbo

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

Вот, для примера комплексированный БИХ и КИХ фильтр. Результаты лучше, чем по отдельности у каждого.
Вводя управление полосой среза для БИХ фильтра ( по примеру мной упомянутого ) сглаженные верхушки практически приближаются к исходной трапеции (на снапе - без управления полосой, потому верхушки сглажены).
Пояснения:
- исходный сигнал две синусоиды, ограниченные по амплитуде;
- шумы ( белый гауссовский некоррелированный с дисперсией 0.1 + пиковый с амплитудой 0.1 и вероятностью 0.25 на периоде )
- комплексный фильтр ( БИХ второго порядка и оконный ( Hamming ) 27 точек )
- формат чисел - плавающий.

Никаких предположений о полосе полезного сигнала не делалось, кроме того, что она Low-Pass.
Фазовый сдвиг будет, да, а иначе - к экстрасексам :)

Экран симулятора:
http://s61.radikal.ru/i172/0912/b2/3c4edfb4f9a0.jpg

← →
MBo © (2009-12-04 17:02) [7]

Да пока на модельной системе Savitzky-Golay хорошо себя показывает.

http://i081.radikal.ru/0912/28/48e5f8ea20ab.gif

средняя картинка - увеличенный излом.

Красный - исходный сигнал, СКЗ шума 0.95 единиц, синий - сглаживающий сплайн, не совсем удачный подбор параметров (выброс на изломе), черный - Савицкого-Голея 4-го порядка, белый - вполне достойная производная, получена Сав-Гол. второго порядка из исходного сигнала.

Внизу - сигнал с невысоким уровнем. Здесь видно, что результаты фильтра СГ нетрудно будет почистить дополнительным фильтром с небольшой апертурой, что даст приятность глазу.

← →
Jeer © (2009-12-04 17:17) [8]

Ok.

← →
MBo © (2009-12-04 17:24) [9]

вот данные (не те, что на рисунке)
http://www.onlinedisk.ru/file/284877/

← →
Jeer © (2009-12-04 17:31) [10]

Ok.

← →
Jeer © (2009-12-09 17:12) [11]

Извиняюсь за задержку с ответом - навалилось под конец года.

Итак, по существу.

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

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

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

3. Зона нарастания.
Вероятно реакция объекта на некоторые увеличивающиеся входные воздействия разной интенсивности.
Не "гладкий" характер реакции свидетельствует скорее о не "гладкости" задающего сигнала.

4. Зона "гребня".
Явно присутствует изменение задающего сигнала на уменьшение. Поскольку в момент смены задающего сигнала объект не достиг установившихся параметров, эта зона интересна только градиентами.

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

6. Зона выключения.
Объект, находясь в процессе снижения параметра, получил ступеньку задающего сигнала ( выключение ).
Эффективная постоянная времени примерно аналогична зоне 2.

Выводы:
Для идентификации параметров модели объекта допустимо использование зон 2 и 5.
Для идентификации параметров помех могут быть использованы зоны 1 (холодная зона) и 5 ( для выявления корреляции параметров
сигнала и помех).

На первый взгляд, корреляция отсутствует.

Параметры шумов:
Зона 1:
Mean: 5.33
Std: 0.94
Нормальность - да.
- skewness = 0.039
- kurtosis = -0.029
"Окрашенность" минимальна.

Двухпараметрическая регрессия: snap-02.png
(практически нулевой дрейф)
http://i052.radikal.ru/0912/2f/63c3fd90ca24.png

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

Для упрощения расчетной схемы фильтров, были использованы ( испытаны ) несколько FIR-фильтров на основе "оконного" метода синтеза.
Поскольку для данной модели объекта особых предпочтений между видами формирующих "окон" нет, то приводятся результаты для
фильтра Хемминга ( Hamming ), весовые коэффициенты которого рассчитываются как:

W[i] := a+b*Cos(2*pi*(i-1)/(N-1));

procedure winHamming(var value: double; arg: double);
begin
value := value * (0.54 + 0.46 * Cos(arg));
end;

На верхнем графике - импульсная переходная характеристика фильтра Хемминга.

Выходной сигнал объекта
01.png
http://s57.radikal.ru/i158/0912/03/a22d8cd2f30f.png

Переход от зоны 1 к зоне 2
02.png
http://i006.radikal.ru/0912/0f/71f01016fbb6.png

Зона 4 "гребень"
03.png
http://i048.radikal.ru/0912/d7/11403dcde0ea.png

Переход от зоны 5 "спад" к зоне 6 "выключение"
04.png
http://i056.radikal.ru/0912/25/3ac8acfa518c.png

Зона 5 "Выключение"
05.png
http://s54.radikal.ru/i143/0912/59/f8f2d0f3ff79.png

P.S.
Проверка на более примитивных (по расчету) фильтрах ( скользящее среднее, взвешенное среднее) не показала
особого ухудшения ситуации, впрочем как не показал себя и фильтр Савицкого-Голея в отношении явного улучшения.

P.P.S.
Поскольку за "кадром" остались вопросы использования результатов измерений ( измерение как цель, идентификация параметров, управление состоянием),
а также сложность реализации с учетом выч.ресурсов, то могут потребоваться и другие варианты решения задачи фильтрации.
На мой взгляд, для такого "удобного" объекта и с весьма небольшим уровнем шума простор для выбора оптимального режима фильтрации
весьма велик.
Это и синтез по Ремезу и билинейное преобразование и модели авторегрессии и нелинейная фильтрация, даже CIC-фильтры.

Если есть вопросы - welcome.

← →
Jeer © (2009-12-09 17:30) [12]

Отфильтрованный сигнал, для примера.
Начальный участок около 100 точек - переходной процесс.
Не стал выводить фильтр на среднее значение.
134 kb
http://slil.ru/28318360

← →
MBo © (2009-12-09 18:09) [13]

О, простой Хэмминг отлично сработал.
Я провел уже натурный эксперимент, и объяснили больше о задаче.
Про тепловой объект - действительно так.
Интересует характеристика охлаждения, а именно - зона быстрого спада, главным образом поведение производной.

http://s12.radikal.ru/i185/0912/83/c60996e69fb2.gif
сигнал и производная (dT/dt нарисована как функция от T), температура T по левой оси.

Использовал Савицкого-Голея + рассчитал IIR режекторный фильтр второго порядка для подавления 50 Гц.
Правда, в нужном масштабе отображения сетевой фон особо не мешает, а в начале, пока режекторный фильтр не установился, осцилляции изрядные вносит (длина этого периода зависит от добротности).

Завтра повожусь с Хэммингом. Спасибо за науку.

> Двухпараметрическая регрессия:

Где про регрессию лучше всего почитать?

> Для упрощения расчетной схемы фильтров, были использованы
> ( испытаны ) несколько FIR-фильтров на основе "оконного"
> метода синтеза.Поскольку для данной модели объекта особых
> предпочтений между видами формирующих "окон" нет, то приводятся
> результаты дляфильтра Хемминга ( Hamming ), весовые коэффициенты
> которого рассчитываются как:

Не понял. У тебя на картинке написанно Low-pass FIR фильтер, а тут ты говоришь что используешь в качестве коэффициентов FIR Хэминга используешь. А это же уже ближе к полосовому фильтру.

Или ты все таки используешь окна? Потом уже FIR фильтер. Тогда вопрос как сам FIR фильтр считал выбирал.

> билинейное преобразование

А откуда в одномерном сигнали "би" взялось? Это я к тому что про это я вообще слыхом не слыхивал.

> Pavia © (09.12.09 20:48) [14]

Все неверно тобой понято :)

>Где про регрессию лучше всего почитать?

В умных книжках :) типа "Регресиионный анализ".
Если конкретный вопрос, то спрашивай.

> Не понял. У тебя на картинке написанно Low-pass FIR фильтер,
> а тут ты говоришь что используешь в качестве коэффициентов
> FIR Хэминга используешь. А это же уже ближе к полосовому
> фильтру.
>

Есть разные методы проектирования цифровых фильтров:

- по аналоговому прототипу ( билинейное преобразование - переход от аналоговой передаточной функции к цифровой с помощью аппроксимации z-преобразования) - применяется для синтеза для БИХ-фильтров;
- аналитический синтез оптимальных фильтров через решение вариационной задачи минимизации целевой функции ( методы неопределенных множителей Лагранжа и т.п.)
- "чебышевская" аппрокцимация, например алгоритм Ремеза;
- синтез стохастических моделей типа авторегрессии и др. методами Берга, Бакса-Дженкинса, Прони и др;
- "прямой" синтез по желаемому виду АЧХ ( с использованием весовых функций);
- адаптивные фильтры ( Винера, Калмана и их клоны);

Для Mbo я использовал КИХ-фильтр ( FIR ) нижних частот с использованием "прямого" синтеза.
Для устранения ( снижение вляния) явления Гиббса использовалось конструирование весовых функций, также известное как метод "окон".
Т.е. мы переходим от желаемой АЧХ к импульсной характеристие и "накладываем" на нее одну из "оконных" функций.
В частости, было использовано окно Хемминга (а их много - Ханн, Ланцош, Гаусс, Казер, Барлетт и пр)

Это что касается методов синтеза.

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

Т.е. метод окон вполне может использоваться для проектирования фильтров любых типов.

Кроме того, "оконный" метод применяется для взвешивания исходного сигнала при выполнении преобразования Фурье для устранения того же явления Гиббса.

> > билинейное преобразование
>
> А откуда в одномерном сигнали "би" взялось?

Это один из ранних методов проектирования цифровых рекурсивных фильтров.
Основан на том, что имеются очень подробные справочники по аналоговым фильтрам с уже рассчитанными коэффициентами.
Чтобы не ставить заново задачу синтеза с нуля, используют метод "билинейного" преобразования передаточной функции аналогового фильтра-прототипа в цифровую передаточную функцию.
Еще существует обобщенное билинейное преобразование предложенное Константинидисом для перевода одного типа вида в другой: ФНЧ - ФВЧ, ФВЧ - ПФ и т.д.

> рассчитал IIR режекторный фильтр второго порядка для подавления
> 50 Гц.

Нецелесообразно смешивать IIR и FIR.
Лучше синтезировать решетчатый ( гребенчатый ) FIR-фильтр для устранения периодической помехи 50 Гц и ее гармоник и "добавить" его в основной фильтр, путем коррекции коэфф-тов весовой функции.
Ну или каскадная реализация - не суть.

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

Линейный знаю, интегрирующее звено, дифференцирующее, второго порядка периодическое, второго порядка аппериодическое. Билинейного не встречал.

> Кроме того, "оконный" метод применяется для взвешивания
> исходного сигнала при выполнении преобразования Фурье для
> устранения того же явления Гиббса.

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

> Для устранения ( снижение вляния) явления Гиббса использовалось
> конструирование весовых функций, также известное как метод
> "окон".

Я слышал что FIR фильтер в отличие от FFT не дает искажений Гибсса. Поэтому меня и удивило что вы используете окна. Прокомментируйте?

Еще один вопрос который меня мучает когда мы проектируем FIR фильтр НЧ или ВЧ, то он имеет в частотной области переход. Так вот как лучше проектировать так чтобы этот переход был как можно ближе к ступеньке или плавный. Т.е вопрос что лучше?

Понятно что это зависит от сигнала. И для АИМ лучше даже наличия лепестков. А если у нас вот такой вот сигнал как у MBO то как лучше?

Собственно меня интересует как вы строили фильтр. Как выбирали(определили) частоту среза или что вы там выбирали?

> Нецелесообразно смешивать IIR и FIR.Лучше синтезировать
> решетчатый ( гребенчатый ) FIR-фильтр для устранения периодической
> помехи 50 Гц и ее гармоник и "добавить" его в основной фильтр,
> путем коррекции коэфф-тов весовой функции.Ну или каскадная
> реализация - не суть.

А может лучше IIR так как FIR является подмножеством IIR фильтров?

> Pavia © (10.12.09 20:20) [18]
> Линейный знаю, интегрирующее звено, дифференцирующее, второго
> порядка периодическое, второго порядка аппериодическое.
> Билинейного не встречал.

Еще раз.
Билинейное преобразование, т.е. переход от аналогового прототипа к цифровому IIR-фильтру, а не билинейный фильтр.
Тип и порядок исходного аналогового фильтра ничем не ограничивается.

Выглядит это так:
В передаточной функции аналогового фильтра, записанной в операторном Лапласовском виде как H(s), делается замена ( подстановка) s -> k*(z-1)/(z+1), в результате чего получаем передаточную характеристику цифрового фильтра G(z) в терминах z-преобразования.

> Это я знаю. Ты не думай что я ничего не знаю. Просто вы
> знаете больше.

Хорошо, договоримся на этом. :)

>
> Я слышал что FIR фильтер в отличие от FFT не дает искажений
> Гибсса. Поэтому меня и удивило что вы используете окна.
> Прокомментируйте?

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

> Еще один вопрос который меня мучает когда мы проектируем
> FIR фильтр НЧ или ВЧ, то он имеет в частотной области переход.
> Так вот как лучше проектировать так чтобы этот переход
> был как можно ближе к ступеньке или плавный. Т.е вопрос
> что лучше?

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

Например, есть еще метод синтеза FIR-фильтров по точкам желаемой АЧХ.
Если не принимать дополнительных мер, мы получим АЧХ точно проходящую через желаемые точки, но поведение ее между точками будет малопредсказуемым, вплоть до появления пиков и провалов.

> Понятно что это зависит от сигнала. И для АИМ лучше даже
> наличия лепестков. А если у нас вот такой вот сигнал как
> у MBO то как лучше?

Поскольку у Mbo есть вполне реальный объект, который можно рассматривать как некую динамическую систему, можно сделать оценку параметров ее упрощенной модели по испытаниям, т.е. мы получаем определенную передаточную функцию. А это значит, что зная параметры входного детерминированного сигнала, мы можем определить выходной сигнал, например через интеграл Дюамеля.
Что это нам дает ? Это дает возможность определить спектр полезного сигнала, чтобы потом бороться с помехами наилучшим образом.

Для объекта Mbo это примерно два каскадно-соединенных апериодических звена ( звена первого порядка ). Есть методы, позволяющие делать оценку параметров моделей. Для простых моделей возможна аналитическая оценка, для более сложных - вычислительная.
Итак определив параметры ( постоянные времени ) модели, можно приступать к выбору типа фильтра и метода его синтезу фильтра.
Но предварительно надо разобраться с помехами, а для этого нужна оценка их спектральной плотности. По данным Mbo это можно сделать.
Именно для анализа параметров сигнала и помех я и сделал разбор "полета" выходного сигнала, разбив его на принципиально значимые зоны.

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

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

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

Ну а дальше - эксперименты :)

> А может лучше IIR так как FIR является подмножеством IIR
> фильтров?

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

Ну как-то вот так я ответил бы :)

Фильтр Савицкого-Голея и иже с ним Найти похожие ветки