Фильтр Савицкого – Голея против линейного фильтра IIR или FIR

11
  • Традиционный БИХ / КИХ-фильтр (низкочастотный для удаления колебаний высокой частоты), например, скользящее среднее,

  • или фильтр Савицкого-Голея

все это может быть полезно для сглаживания сигнала, такого как сигнал огибающей:

введите описание изображения здесь

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

Чем он отличается от стандартного фильтра и что он добавляет по сравнению со стандартными фильтрами?

Адаптируется ли он к входным данным?

Это лучше для временного сохранения?


Вы когда-нибудь оказывались в инженерной ситуации, когда решили: «Давайте использовать фильтр SG вместо скользящего среднего или другого FIR lowpass! Лучше, потому что это, это и это ...» ? Тогда этот вопрос для вас!

g6kxjv1ozn
источник

Ответы:

4

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

Прежде всего, я хотел бы повторить то, что выяснилось в дискуссии, порожденной другими ответами: категоризация сглаживающих фильтров в вопросе на линейные и не зависящие от времени (LTI) FIR / IIR фильтры, с одной стороны, и Фильтры Савицкого-Голея, с другой стороны, вводят в заблуждение. Фильтр Савиткзи-Голея - это просто стандартный КИХ-фильтр, разработанный в соответствии с определенным критерием (приближение локального полинома). Таким образом, все фильтры, упомянутые в вопросе, являются фильтрами LTI.

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

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

Если спектр шума существенно перекрывается со спектром сигнала, требуется более осторожный подход, и затухание методом "грубой силы" не будет работать должным образом, потому что либо вы оставляете слишком много шума (выбирая слишком высокую частоту среза), либо искажаете желаемое Слишком сильный сигнал. В этом случае фильтры Savitzky-Golay (SG) могут быть хорошим выбором. Ценой, которую нужно заплатить, является посредственное затухание в полосе пропускания, но одно преимущество состоит в том, что некоторые свойства сигнала сохраняются очень хорошо. Это связано с тем, что фильтры SG имеют плоскую полосу пропускания, т.е.

(1)dКЧАС(еJω)dωК|ωзнак равно0знак равно0Кзнак равно1,2,...,р

где р - порядок аппроксимирующего полинома, а ЧАС(еJω) - частотная характеристика фильтра. Свойство (1) гарантирует, что первые р моментов входного сигнала сохраняются на выходе, что означает, что ширина и высота пиков в требуемом сигнале хорошо сохраняются.

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

Для тех, кто интересуется теорией фильтров SG, наиболее подходящие ссылки, которые я могу порекомендовать, следующие:

Мэтт Л.
источник
2

Как и во всем, иногда одни инструменты лучше других.

Фильтры скользящего среднего (MA) могут использоваться для сглаживания данных и являются FIR. Это в значительной степени самый простой фильтр, который вы можете придумать, и они хорошо работают для многих задач, если вы не пытаетесь смоделировать какие-либо внезапные скачки или полиномиальные тренды. Имейте в виду, что, по сути, это всего лишь фильтр нижних частот, поэтому они лучше всего работают, когда данные, которые вам нужны в сигнале, являются низкочастотными и довольно узкополосными.

Фильтры Савицкого-Голея (SG) представляют собой особую группу КИХ-фильтров, которые по существу соответствуют полиному для вашего временного ряда, когда свертка скользит вдоль сигнала. Фильтры SG полезны для сигналов, в которых вы заботитесь не обязательно о низкой частоте и достаточно узкой полосе.

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

РЕДАКТИРОВАТЬ :

Поскольку кажется, что существует некоторая путаница, что фильтры SG в некотором роде являются «адаптивными» на базовом уровне, я включил простой пример MATLAB. Как отметил Дэн, их можно сделать адаптивными, но их базовая реализация часто таковой не является. Изучив код, вы увидите, что это просто матричный поиск с некоторой специальной обработкой. Ничто в этом фильтре не является «адаптивным» в традиционном смысле: вы просто выбираете полиномиальную подгонку и длину, на которую этот полином будет вписываться в сигнал посредством свертки; СГ имеет важное значение FIR. Сценарий, который я имею ниже, производит этот сюжет: Сравнение фильтра SG с MA

Глядя на этот рисунок, вы можете видеть, что MA и SG по существу выполняют одно и то же, но с некоторыми важными различиями:

  1. МА отлично справляется с подавлением шума, но плохо справляется с захватом переходных скачков в сигнале. Мы можем подавить еще больший шум, увеличив длину фильтра, но тогда он будет работать еще хуже на переходных процессах; этот эффект будет выглядеть как «размазывание» вблизи любых переходных процессов, что вы должны увидеть на рисунке.
  2. SG отлично справляется с захватом переходных процессов в сигнале, но не так хорошо справляется с подавлением шума (по крайней мере, по сравнению с MA такого же размера). Мы можем улучшить подавление шума вблизи переходных процессов, увеличив длину кадра, но это приведет к появлению звонка, аналогичного феномену Гибба, вблизи любых переходных процессов.

Чтобы вы лучше поняли, как работают эти фильтры, я бы посоветовал вам взять здесь код, манипулировать им и посмотреть, как работают все отдельные части фильтра SG.

Код для примера MATLAB:

% Generate a signal "s" that has square waves, and scale it with a
% polynomial of order 5
up = 1*ones(1,100);
down = zeros(1,150);
s = [down down up up down up down up down up up up down down down down down];
n = numel(s);
nn = linspace(0,4,numel(s));
s = s .* (nn .^5);
sn = (s + 4*randn(size(s))).';

% smooth it with SMA of length 15
sz = 15;
h = 1/sz * ones(1,sz);
sn_sma = conv(sn,h,'same');

% smooth it with sgolay of frame length 15
m = (sz-1)/2;
% look up the SG matrix for this order and size
B = sgolay(5, sz);

% compute the steady state response for the signal, i.e. everywhere that
% isnt the first or last "frame" for the SG filter
steady = conv(sn, B(m+1,:), 'same');
% handle the transiet portion at the start of the signal
y_st   = B(1:m,:)*sn(1:sz);
% handle the transient portion at the end of the signal
y_en   = B(sz -m+1 : sz, :) * sn(n - sz+1:n);

% combine our results
sn_sg  = steady;
sn_sg(1:m) = y_st;
sn_sg(n-m+1:n) = y_en;

% plots
figure(1);
hold off;
plot(sn,'Color',[0.75 0.75 0.75]);
hold on;
plot(sn_sma,'b');
plot(sn_sg,'r');

legend('Noisy Signal','MA Smoothing','SG Smoothing, order 5','Location','NorthWest');
matthewjpollard
источник
1
Похоже, дело в том (увидев другой ответ), что SG-фильтр является «полностью зависящим от данных нелинейным изменяющимся во времени фильтром , коэффициенты которого пересчитываются для каждого короткого сегмента его входных данных».
g6kxjv1ozn
1
Фильтр SG, насколько я понимаю от его применения несколько раз, вовсе не является адаптивным фильтром, особенно по сравнению с вашим средним адаптивным фильтром, таким как LMS или RLS. Я бы полностью не согласился с утверждением, что вес фильтра меняется во времени. Фильтры SG по сути являются поиском в таблице, вы фильтруете значения из таблицы для вычисления переходного ответа, а затем возвращаетесь и обрабатываете крайние случаи в начале / конце сигнала. Я отредактирую свой пост на примере MATLAB, чтобы показать это вам.
Matthewjpollard
2
@matthewjpollard Следует отметить, что лично у меня нет значительного опыта использования этого фильтра, но мне кажется, что наилучшим образом реализованный фильтр SG выглядит как реализация адаптивного фильтра с изменяющимися во времени коэффициентами. То, как вы применили фильтр в своем коде, не является (так как вы рассматривали всю последовательность как «подмножество данных»), а именно так, как описано в Wikipedia en.wikipedia.org/wiki/Savitzky%E2%80% 93Golay_filter и в самой статье Савицкого и Голея действительно адаптивны: pdfs.semanticscholar.org/4830/…
Дэн Бошен,
2
@matthewjpollard В ваших системах реального времени ваши данные постоянно передаются в потоковом режиме, так что вы пересчитываете коэффициенты через более короткие промежутки времени или вы всегда работаете с небольшими блоками данных?
Дэн Бошен
2
Спасибо Мэтт. Поэтому, возможно, мы могли бы связать ваши действия как адаптивные / изменяющиеся во времени в том смысле, что коэффициенты вычисляются для каждого сбора данных (однако одни и те же коэффициенты в пределах сбора данных при правильной обработке начала и конца, но варьируются от одного сбора к следующему, если я правильно понимаю). Спасибо, что поделились своим кодом и примером приложения.
Дэн Бошен
2

НОТА

мой предыдущий ответ (перед этим редактированием), в котором было указано, что фильтр Савицкого-Голея (SG) является нелинейной изменяющейся во времени зависимой от входных данных, был неверным из-за преждевременной неверной интерпретации того, как фильтр Савицкого-Голея (SG) вычисляет свои выходные данные по предоставленной вики-ссылке. Так что теперь я исправляю это для тех, кто также увидит, как фильтры SG реализуются фильтрацией FIR-LTI. Благодаря @MattL. за его исправление, за отличную связь, которую он обеспечил, и за терпение, которое он имел (что я никогда не мог проявить) во время моего исследования проблемы. Хотя я бы честно предпочел более подробные возражения, которые, тем не менее, явно не нужны. Также обратите внимание, что правильный ответ - другой, этот только для дополнительного разъяснения свойства LTI фильтров SG.

Теперь неудивительно, что когда кто-то (кто никогда раньше не использовал эти фильтры) сталкивается с определением фильтра SG как полиномиального соответствия LSE низкого порядка с данными, он / она сразу же приходит к выводу, что они зависят от данных, нелинейны и изменяющиеся во времени (смены), адаптивные фильтры.

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

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

2M+1Икс[-M],Икс[-M+1],,,,,Икс[0],Икс[1],,,,,Икс[M]Nзнак равно0п[N]NNзнак равно-M,-M+1,,,,,-1,0,1,,,,M

п[N]знак равноΣКзнак равно0NaКNКзнак равноa0+a1N+a2N2+,,,+aNNN

aКNTчасп[N]

Езнак равноΣ-MM(п[N]-Икс[N])2

Иксзнак равно[Икс[-M],Икс[-M+1],,,,,Икс[0],Икс[1],,,,,Икс[M]]T

aКЕ

(1)Еaязнак равно0   ,   за    язнак равно0,1,,,,N

Теперь для тех, кто знаком с процедурой полифизики LSE, я просто напишу полученное матричное уравнение (по ссылке), которое определяет оптимальный набор коэффициентов:

(2)aзнак равно(ATA)-1ATИксзнак равноЧАСИкс
Икс(2M+1)×1ЧАС2M+1NANAЧАСA

Aзнак равно[αN,я]знак равно[(-M)0(-M)1,,,(-M)N(-M+1)0(-M+1)1,,,(-M+1)N,,,(0)0(0)1,,,(0)N,,,(M)0(M)1,,,(M)N]

Теперь давайте на мгновение откинемся назад и обсудим вопрос здесь.

AЧАСNaКMNИкс[N]aК2Nd

... Это (полифит LSE) может повторяться для каждой выборки входных данных, каждый раз создавая новый многочлен и новое значение выходной последовательности y [n] ...

Так как же нам преодолеть этот удивительный сюрприз? Путем интерпретации и определения вывода фильтра SG следующим образом:

NNИкс[N]Y[N]п[N]Nзнак равно0

Y[N]знак равноY[0]знак равноΣмзнак равно0NaмNмзнак равноa0

2M+1Икс[N]Nзнак равноdY[N]a0п[N]Икс[N]Nзнак равноdY[d]Икс[d-M],Икс[d-M+1],,,,,Икс[d-1],Икс[d],Икс[d+1],,,,Икс[d+M]

a0Икс[N]Y[N]Икс[N]NИкс[N]час[N], Но тогда, каковы коэффициенты фильтра для этого фильтра SG? Посмотрим.

aК

aзнак равноЧАСИкс

[a0a1aN]знак равно[час(0,0)час(0,1),,,час(0,2M)час(1,0)час(1,1),,,час(1,2M),,,час(N,0)час(0,1),,,час(0,2M)][Икс[-M]Икс[-M+1],,,Икс[M]]

a0ЧАСИкс

a0знак равноЧАС(0,N)Иксзнак равноΣЧАС(0,К)Икс[К]знак равноЧАС(0,-N)Икс[N]

час[N]знак равноЧАС(0,-N)

N2M+1

Y[N]2M+1Икс[N]LчасN[N]

Y[N]знак равноИкс[N]часN[N]

КОММЕНТАРИЙ

aКчас[N]Y[N]Иксaзнак равноЧАСИксaКп[N]aКчас[N]

MATLAB / OCTVE CODE

час[N]час[N]

% Savitzky-Golay Filter
% 
clc; clear all; close all;

N = 3;                      % a0,a1,a2,a3 : 3rd order polynomial
M = 4;                      % x[-M],..x[M] . 2M + 1 data

A = zeros(2*M+1,N+1);
for n = -M:M
    A(n+M+1,:) = n.^[0:N];
end

H = (A'*A)^(-1)* A';        % LSE fit matrix

h = H(1,:);                 % S-G filter impulse response (nancausal symmetric FIR)

figure,subplot(2,1,1)
stem([-M:M],h);
title(['Impulse response h[n] of Savitzky-Golay filter of order N = ' num2str(N), ' and window size 2M+1 =  ' , num2str(2*M+1)]);

subplot(2,1,2)
plot(linspace(-1,1,1024), abs(fftshift(fft(h,1024))));
title('Frequency response magnitude of h[n]');

Выход:

введите описание изображения здесь

Надеюсь, что это проясняет проблему.

Fat32
источник
2
@ Fat32 Я думаю, это потому, что это был длинный список комментариев, поэтому, чтобы держать доску в чистоте, они обычно перемещают ее «в чат». Это все еще там, только не загромождая главную страницу. Вот почему система предлагает переместить его в чат, когда вперед и назад становится длинным. Не волнуйтесь, все еще любят вас.
Дэн Бошен,
1
@ g6kxjv1ozn Я согласен с этим ... пожалуйста, подождите ...
Fat32
2
@ Fat32 Отличная работа! Я прочитал это, но мне нужно будет прочитать это, и оно написано очень четко, мне просто нужно будет пройтись карандашом и бумагой шаг за шагом, чтобы полностью увидеть это, как вы сейчас делаете. Спасибо, что выложили все это здесь.
Дэн Бошен
4
1ωзнак равно0
2
aКчас[N]