sgolay

Создание фильтра Savitzky-Golay

Описание

пример

b = sgolay(order,framelen) проектирует КИХ-фильтр сглаживания Savitzky-Golay с полиномиальным порядком order и структурируйте длину framelen.

b = sgolay(order,framelen,weights) задает вектор взвешивания, weights, который содержит действительные, веса с положительным знаком, которые будут использоваться во время минимизации наименьших квадратов.

пример

[b,g] = sgolay(___) возвращает матричный g из фильтров дифференцирования. Можно использовать эти выходные аргументы с любым из предыдущих входных синтаксисов.

Примеры

свернуть все

Сгенерируйте сигнал, который состоит из синусоиды на 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));

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

y = [ybegin; ycenter; yend];
plot([x y])
legend('Noisy Sinusoid','S-G smoothed sinusoid')

Figure contains an axes object. The axes object contains 2 objects of type line. These objects represent 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));

Оцените первые три производные синусоиды с помощью метода Savitzky-Golay. Используйте системы координат с 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')

Figure contains an axes object. The axes object with title Savitzky-Golay Derivative Estimates contains 5 objects of type line. These objects represent x, x (smoothed), x', x'', x'''.

Входные параметры

свернуть все

Полиномиальный порядок в виде положительного целого числа. Значение order должно быть меньшим, чем framelen. Если order = framelen – 1, затем спроектированный фильтр не производит сглаживания.

Структурируйте длину в виде положительного нечетного целого числа. Значение framelen должен быть больше order.

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

Выходные аргументы

свернуть все

Изменяющиеся во времени КИХ-коэффициенты фильтра в виде framelen- framelen матрица. В сглаживании фильтруют реализацию (например, sgolayfilt), последний (framelen-1)/2 строки (каждый КИХ-фильтр) применяются к сигналу во время переходного процесса запуска и первому (framelen-1)/2 строки применяются к сигналу во время терминального переходного процесса. Центральная строка применяется к сигналу в устойчивом состоянии.

Матрица дифференцирования фильтрует в виде матрицы. Каждый столбец g фильтр дифференцирования для производных порядка p-1, где p индекс столбца. Учитывая x сигнала из длины framelen, можно найти оценку pth закажите производную, xp, из его среднего значения от xp((framelen+1)/2) = (factorial(p)) * g(:,p+1)' * x.

Алгоритмы

Фильтры Savitzky-Golay используются, чтобы сгладить сигналы с шумом с большим промежутком частоты. Savitzky-Golay сглаживающие фильтры имеют тенденцию отфильтровывать меньше высокочастотного содержимого сигнала, чем стандартные КИХ-фильтры усреднения. Однако они менее успешны при отклонении шума, когда уровень шума особенно высок.

В общем случае фильтрация состоит из замены каждой точки сигнала некоторой комбинацией значений сигналов, содержавшихся в движущемся окне, сосредоточенном в точке, при условии, что соседние точки измеряют почти то же базовое значение. Например, фильтры скользящего среднего значения заменяют каждую точку данных на локальное среднее значение окружающих точек данных. Если точка определенных данных имеет точки k налево и точки k направо для общей длины окна L = 2k + 1, фильтр скользящего среднего значения делает замену

xsx^s=1Lr=kkxs+r.

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

[xskxs1xsxs+1xs+k]=[b0+b1(tskΔt)+b2(tskΔt)2++bn(tskΔt)nb0+b1(ts1Δt)+b2(ts1Δt)2++bn(ts1Δt)nb0+b1(ts0Δt)+b2(ts0Δt)2++bn(ts0Δt)nb0+b1(ts+1Δt)+b2(ts+1Δt)2++bn(ts+1Δt)nb0+b1(ts+kΔt)+b2(ts+kΔt)2++bn(ts+kΔt)n]=[a0+a1(k)+a2(k)2++an(k)na0+a1(1)+a2(1)2++an(1)na0+a1(0)+a2(0)2++an(0)na0+a1(1)+a2(1)2++an(1)na0+a1(k)+a2(k)2++an(k)n]

или, в терминах матриц,

x=[1k(k)2(k)n112(2)2(2)n11(1)2(1)n100011121n12222n11kk2kn][a0an]Ha.

Чтобы найти оценки Savitzky-Golay, используйте псевдоинверсию H, чтобы вычислить a и затем предварительно умножиться H:

x^=H(HTH)1HTx=Bx.

Постараться не плохо обусловливать, sgolay использует qr функция, чтобы вычислить разложение размера экономики H как H = QR, в терминах который B = QQT.

Необходимо вычислить B только однажды. Savitzky-Golay оценивает для большей части результата точек сигнала свертки к сигналу с центральной строкой B. Результатом является установившийся фрагмент отфильтрованного сигнала. Первые строки k B дают к начальному переходному процессу, и итоговые строки k B дают к итоговому переходному процессу. Смотрите sgolayfilt для примера. Возможно улучшить подавление шумов путем увеличения длины окна, но это вводит вызов, аналогичный Явлению Гиббса около любых переходных процессов.

Ссылки

[1] Orfanidis, Софокл Дж. Введение в обработку сигналов. Englewood Cliffs, NJ: Prentice Hall, 1996.

[2] Нажмите, Уильям Х., Саул А. Теукольский, Уильям Т. Веттерлинг и Брайан П. Флэннери. Числовые рецепты в C: Искусство научных вычислений. Нью-Йорк: Издательство Кембриджского университета, 1992.

[3] Schafer, Рональд В., “Что такое Фильтр Savitzky-Golay? [Читайте лекции Примечаниям]”. Издание 28 Журнала Обработки сигналов IEEE, Номер 4, июль 2011, стр 111–117. https://doi.org/10.1109/MSP.2011.941097.

Расширенные возможности

Генерация кода C/C++
Генерация кода C и C++ с помощью MATLAB® Coder™.

Смотрите также

| | |

Представлено до R2006a