Конструкция фильтра Савицкого-Голая
Создайте сигнал, который состоит из синусоиды 0,2 Гц, встроенной в белый гауссов шум и дискретизированный пять раз в секунду в течение 200 секунд.
dt = 1/5; t = (0:dt:200-dt)'; x = 5*sin(2*pi*0.2*t) + randn(size(t));
Использовать sgolay для сглаживания сигнала. Используйте 21-образные кадры и многочлены четвертого порядка.
order = 4; framelen = 21; b = sgolay(order,framelen);
Вычислите установившуюся часть сигнала, свернув ее с центральной строкой b.
ycenter = conv(x,b((framelen+1)/2,:),'valid');Вычислите переходные процессы. Использовать последние строки b для запуска и первых строк b для терминала.
ybegin = b(end:-1:(framelen+3)/2,:) * x(framelen:-1:1); yend = b((framelen-1)/2:-1:1,:) * x(end:-1:end-(framelen-1));
Соедините переходные процессы и установившуюся часть для формирования полного сглаженного сигнала. Постройте график исходного сигнала и оценки Савицкого-Голая.
y = [ybegin; ycenter; yend]; plot([x y]) legend('Noisy Sinusoid','S-G smoothed sinusoid')
Создайте сигнал, который состоит из синусоиды 0,2 Гц, встроенной в белый гауссов шум и дискретизированный четыре раза в секунду в течение 20 секунд.
dt = 0.25; t = (0:dt:20-1)'; x = 5*sin(2*pi*0.2*t)+0.5*randn(size(t));
Оцените первые три производные синусоиды по методу Савицкого - Голая. Используйте 25-образные кадры и многочлены пятого порядка. Разделить столбцы на степени dt для правильного масштабирования производных.
[b,g] = sgolay(5,25); dx = zeros(length(x),4); for p = 0:3 dx(:,p+1) = conv(x, factorial(p)/(-dt)^p * g(:,p+1), 'same'); end
Постройте график исходного сигнала, сглаженной последовательности и оценок производных.
plot(x,'.-') hold on plot(dx) hold off legend('x','x (smoothed)','x''','x''''', 'x''''''') title('Savitzky-Golay Derivative Estimates')
Фильтры Савицкого - Голая используются для сглаживания шумных сигналов с большим частотным диапазоном. Сглаживающие фильтры Савицки-Голая имеют тенденцию отфильтровывать меньше высокочастотного содержимого сигнала, чем стандартные усредняющие фильтры FIR. Однако они менее успешно отклоняют шум, когда уровни шума особенно высоки.
В общем, фильтрация состоит из замены каждой точки сигнала некоторой комбинацией значений сигнала, содержащихся в движущемся окне, центрированном в точке, исходя из предположения, что близлежащие точки измеряют почти одно и то же нижележащее значение. Например, фильтры скользящего среднего заменяют каждую точку данных локальным средним для окружающих точек данных. Если данная точка данных имеет k точек слева и k точек справа, для общей длины окна L = 2k + 1 фильтр скользящего среднего производит замену
Фильтры Савицки-Голая обобщают эту идею наименьшими квадратами, подгоняющими многочлен n-го порядка через значения сигнала в окне и принимающими вычисленную центральную точку аппроксимированной полиномиальной кривой в качестве новой сглаженной точки данных. Для данной точки xs,
na0+a1 (1) +a2 (1) 2 + ⋯ + (1) n⋮a0+a1 (k) +a2 (k) 2 + ⋯ + (k) n]
или, в терминах матриц,
n100⋯01112⋯1n1222⋯2n1⋮⋮⋱⋮1kk2⋯kn] [a0⋮an]≡Ha.
Чтобы найти оценки Савицки-Голая, используйте псевдоинверсию H, чтобы вычислить a и затем предварительно умножить на H:
1HTx = Bx.
Во избежание плохого кондиционирования, sgolay использует qr функция для вычисления разложения H по экономичному размеру как H = QR, в терминах которого B = QQT.
Вычислять B необходимо только один раз. Оценки Савицки-Голая для большинства точек сигнала являются результатом свертки сигнала с центральной строкой В. Результатом является установившаяся часть отфильтрованного сигнала. Первые k рядов В дают начальный переходный процесс, а конечные k рядов В дают конечный переходный процесс. Посмотрите sgolayfilt например. Можно улучшить подавление шума за счет увеличения длины окна, но это вводит звон, аналогичный явлению Гиббса вблизи любых переходных процессов.
[1] Орфанидис, Софокл Дж. Введение в обработку сигналов. Энглвуд Клиффс, Нью-Джерси: Прентис Холл, 1996.
[2] Пресса, Уильям Х., Сол А. Теукольский, Уильям Т. Веттерлинг и Брайан П. Фланнери. Численные рецепты в C: Искусство научных вычислений. Нью-Йорк: Cambridge University Press, 1992.
[3] Шефер, Рональд У. "Что такое фильтр Савицки-Голея? [Примечания к лекциям]. " IEEE Signal Processing Magazine Vol. 28, Number 4, July 2011, pp. 111-117. https://doi.org/10.1109/MSP.2011.941097.
filter | fir1 | firls | sgolayfilt