sgolay
Расчет фильтра Савицкого — Голея.
| Библиотека |
|
Аргументы
Выходные аргументы
#
b —
коэффициенты КИХ-фильтра, изменяющиеся во времени
матрица
Details
Коэффициенты КИХ-фильтра, изменяющиеся во времени, возвращаемые в виде матрицы fl×fl. В реализации сглаживающего фильтра, последние (fl−1)/2 строк (каждая — КИХ-фильтр) применяются к сигналу во время начального переходного процесса, а первые (fl−1)/2 строк — во время конечного переходного процесса. Центральная строка применяется к сигналу в устойчивом состоянии.
#
g —
матрица дифференциальных фильтров
матрица
Примеры
Сглаживание зашумленной синусоиды по методу Савицкого — Голея
Details
Сгенерируем сигнал, представляющий собой синусоиду с частотой 0.2 Гц, к которой добавлен белый гауссовский шум, с частотой дискретизации пять раз в секунду в течение 200 секунд.
import EngeeDSP.Functions: randn, sgolay,conv
dt = 1/5
t = (0:dt:200-dt)'
x = 5*sin.(2*pi*0.2*t) + randn(1, length(t))Используем sgolay для сглаживания сигнала. Используем кадры с 21 отсчетом и полиномы четвертой степени.
order = 4
framelen = 21
b,g = sgolay(order,framelen)Вычислим установившуюся часть сигнала, свернув его со средним рядом вектора b.
ycenter = conv(x,b[Int((framelen+1)/2), :],"valid")'Вычислим переходные процессы. Используем последние строки матрицы b для запуска и первые строки матрицы b для завершения.
ybegin = b[end:-1:Int((framelen+3)/2), :] * x[framelen:-1:1]
yend = b[Int((framelen-1)/2):-1:1, :] * x[end:-1:end-(framelen-1)]Объединим переходные процессы и установившуюся часть сигнала, чтобы получить полный сглаженный сигнал. Построим график исходного сигнала и сигнала, сглаженного по методу Савицкого — Голея.
y = [ybegin; ycenter; yend]
plot(x', label = "Noisy Sinusoid")
plot!(y, label = "S-G smoothed sinusoid")
Фильтр Савицкого — Голея с весами
Details
Спроектируем сглаживающий фильтр Савицкого — Голея с полиномами 5-го порядка для фильтрации зашумленного сигнала, дискретизированного с частотой 8192 Гц. Используем окно Хэмминга. Длина кадра составляет 101 отсчет.
import EngeeDSP.Functions: hamming, freqz, sgolay
Fs = 8192
m = 5
fl = 101
weights = hamming(fl)
b,g = sgolay(m,fl,weights)Построим частотную характеристику для центральной строки b. Используем 1024 точки БПФ.
NFFT = 1024
b=b[Int((fl+1)/2), :]
freqz(b',1,NFFT,Fs,out=:plot)
Дифференцирование по методу Савицкого — Голея
Details
Сгенерируем сигнал, представляющий собой синусоиду с частотой 0.2 Гц, к которой добавлен в белый гауссовский шум, с частотой дискретизации четыре раза в секунду в течение 20 секунд.
import EngeeDSP.Functions: randn, sgolay,conv
dt = 0.25
t = (0:dt:20-1)'
x = 5*sin.(2*pi*0.2*t) + 0.5*randn(1, length(t))Оценим первые три производные синусоиды, используя метод Савицкого — Голея. Используем 25-точечные окна и полиномы пятого порядка. Разделим столбцы на степени dt, чтобы правильно масштабировать производные.
b,g = sgolay(5,25)
dx = zeros(length(x), 4)
for p in 0:3
kernel = factorial(p) / (-dt)^p * g[:, p+1]
dx[:, p+1] = conv(x, kernel, "same")
endПостроим графики исходного сигнала, сглаженной последовательности и производных.
plot(x', marker = :., label = "x")
plot!(dx, label = ["x (smoothed)" "x'" "x''" "x'''"])
title!("Savitzky-Golay Derivative Estimates")
Алгоритмы
Фильтры Савицкого — Голея используются для сглаживания зашумленных сигналов с большим диапазоном частот. Сглаживающие фильтры Савицкого — Голея, как правило, отфильтровывают меньше высокочастотного содержимого сигнала, чем стандартные усредняющие КИХ-фильтры. Однако они менее эффективны в подавлении шума, когда его уровень особенно высок.
В общем случае фильтрация заключается в замене каждой точки сигнала некоторой комбинацией значений сигнала, содержащихся в скользящем окне с центром в этой точке, исходя из предположения, что близлежащие точки измеряют почти одно и то же базовое значение. Например, фильтры скользящего среднего заменяют каждую точку данных локальным средним значением окружающих точек данных. Если точка данных имеет точек слева и точек справа, при общей длине окна , то фильтр скользящего среднего производит замену
Фильтры Савицкого — Голея обобщают эту идею, аппроксимируя по методу наименьших квадратов полином порядка к значениям сигнала в окне и принимая вычисленную центральную точку кривой полинома за новую сглаженную точку данных. Для данной точки
или, в терминах матриц,
Чтобы найти оценки Савицкого — Голея, используем псевдоинверсию для вычисления , а затем умножим слева на :
Чтобы избежать плохой обусловленности, функция sgolay использует QR-разложение для вычисления экономного разложения как , в терминах которого .
Вычислять необходимо только один раз. Оценки Савицкого — Голея для большинства точек сигнала получаются в результате свертки сигнала с центральной строкой . В результате получается установившаяся часть отфильтрованного сигнала. Первые строк в дают начальный переходный процесс, а последние строк в — конечный переходный процесс. Можно улучшить подавление шума, увеличив длину окна, но это приводит к появлению звона, аналогичного явлению Гиббса, вблизи любых переходных процессов.
Литература
-
Orfanidis, Sophocles J. Introduction to Signal Processing. Englewood Cliffs, NJ: Prentice Hall, 1996.
-
Press, William H., Saul A. Teukolsky, William T. Vetterling, and Brian P. Flannery. Numerical Recipes in C: The Art of Scientific Computing. New York: Cambridge University Press, 1992.
-
Schafer, Ronald W. «What Is a Savitzky-Golay Filter? [Lecture Notes].» IEEE Signal Processing Magazine Vol. 28, Number 4, July 2011, pp. 111-117. https://doi.org/10.1109/MSP.2011.941097.