sgolay

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

Синтаксис

b = sgolay(order,framelen)

b = sgolay(order,framelen,weights)

[b,g] = sgolay(___)

Описание

b = sgolay(order,framelen) проектирует фильтр сглаживания конечной импульсной характеристики Савицкого-Голея с полиномиальным порядком order и длина системы координат framelen.

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

пример

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

Примеры

свернуть все

Савицкий-Голай сглаживание шумного синусоида

Открыть Live Script

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

Figure contains an axes. The axes contains 2 objects of type line. These objects represent Noisy Sinusoid, S-G smoothed sinusoid.

Дифференциация Савицкого-Голая

Открыть Live Script

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

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

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

свернуть все

`order` - Полиномиальный порядок
положительное целое число

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

`framelen` - Длина системы координат
положительное нечетное целое число

Длина системы координат, заданная как положительное нечетное целое число. Значение framelen должно быть больше order.

`weights` - Вектор взвешивания
действительный положительный вектор

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

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

свернуть все

`b` - Изменяющиеся во времени коэффициенты конечной импульсной характеристики
матрица

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

`g` - Матрица дифференцировочных фильтров
матрица

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

Алгоритмы

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

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

$x_{s} \to {\hat{x}}_{s} = \frac{1}{L} \sum_{r = - k}^{k} x_{s + r} .$

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

$[\begin{matrix} x_{s - k} \\ ⋮ \\ x_{s - 1} \\ x_{s} \\ x_{s + 1} \\ ⋮ \\ x_{s + k} \end{matrix}] = [\begin{matrix} b_{0} + b_{1} (t_{s} - k Δ t) + b_{2} {(t_{s} - k Δ t)}^{2} + \dots + b_{n} {(t_{s} - k Δ t)}^{n} \\ ⋮ \\ b_{0} + b_{1} (t_{s} - 1 Δ t) + b_{2} {(t_{s} - 1 Δ t)}^{2} + \dots + b_{n} {(t_{s} - 1 Δ t)}^{n} \\ b_{0} + b_{1} (t_{s} - 0 Δ t) + b_{2} {(t_{s} - 0 Δ t)}^{2} + \dots + b_{n} {(t_{s} - 0 Δ t)}^{n} \\ b_{0} + b_{1} (t_{s} + 1 Δ t) + b_{2} {(t_{s} + 1 Δ t)}^{2} + \dots + b_{n} {(t_{s} + 1 Δ t)}^{n} \\ ⋮ \\ b_{0} + b_{1} (t_{s} + k Δ t) + b_{2} {(t_{s} + k Δ t)}^{2} + \dots + b_{n} {(t_{s} + k Δ t)}^{n} \end{matrix}] = [\begin{matrix} a_{0} + a_{1} (- k) + a_{2} {(- k)}^{2} + \dots + a_{n} {(- k)}^{n} \\ ⋮ \\ a_{0} + a_{1} (- 1) + a_{2} {(- 1)}^{2} + \dots + a_{n} {(- 1)}^{n} \\ a_{0} + a_{1} (0) + a_{2} {(0)}^{2} + \dots + a_{n} {(0)}^{n} \\ a_{0} + a_{1} (1) + a_{2} {(1)}^{2} + \dots + a_{n} {(1)}^{n} \\ ⋮ \\ a_{0} + a_{1} (k) + a_{2} {(k)}^{2} + \dots + a_{n} {(k)}^{n} \end{matrix}]$

или, с точки зрения матриц,

$x = [\begin{matrix} 1 & - k & {(- k)}^{2} & \dots & {(- k)}^{n} \\ 1 & ⋮ & ⋮ & ⋰ & ⋮ \\ 1 & - 2 & {(- 2)}^{2} & \dots & {(- 2)}^{n} \\ 1 & - 1 & {(- 1)}^{2} & \dots & {(- 1)}^{n} \\ 1 & 0 & 0 & \dots & 0 \\ 1 & 1 & 1^{2} & \dots & 1^{n} \\ 1 & 2 & 2^{2} & \dots & 2^{n} \\ 1 & ⋮ & ⋮ & ⋱ & ⋮ \\ 1 & k & k^{2} & \dots & k^{n} \end{matrix}] [\begin{matrix} a_{0} \\ ⋮ \\ a_{n} \end{matrix}] \equiv H a .$

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

$\hat{x} = H {(H^{T} H)}^{- 1} H^{T} x = B x .$

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

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

Ссылки

[1] Orfanidis, Sophocles J. Введение в обработку сигналов. Englewood Cliffs, Нью-Джерси: Prentice Hall, 1996.

[2] Press, William H., Saul A. Teukolsky, William T. Vetterling, and Brian P. Flannery. Численные рецепты в C: Искусство научных вычислений. New York: Cambridge University Press, 1992.

[3] Schafer, Ronald W. "What is a Savitzky-Golay Filter? [Примечания к лекциям] ". Журнал обработки сигналов IEEE, том 28, номер 4, июль 2011, стр. 111-117. https://doi.org/10.1109/MSP.2011.941097.

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

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

См. также

filter | fir1 | firls | sgolayfilt

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

Документация

sgolay

Синтаксис

Описание

Примеры

Савицкий-Голай сглаживание шумного синусоида

Дифференциация Савицкого-Голая

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

`order` - Полиномиальный порядок
положительное целое число

`framelen` - Длина системы координат
положительное нечетное целое число

`weights` - Вектор взвешивания
действительный положительный вектор

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

`b` - Изменяющиеся во времени коэффициенты конечной импульсной характеристики
матрица

`g` - Матрица дифференцировочных фильтров
матрица

Алгоритмы

Ссылки

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

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

См. также

Signal Processing Toolbox документация

Поддержка

Документация

sgolay

Синтаксис

Описание

Примеры

Савицкий-Голай сглаживание шумного синусоида

Дифференциация Савицкого-Голая

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

order - Полиномиальный порядок положительное целое число

framelen - Длина системы координат положительное нечетное целое число

weights - Вектор взвешивания действительный положительный вектор

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

b - Изменяющиеся во времени коэффициенты конечной импульсной характеристики матрица

g - Матрица дифференцировочных фильтров матрица

Алгоритмы

Ссылки

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

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

См. также

Signal Processing Toolbox документация

Поддержка

`order` - Полиномиальный порядок
положительное целое число

`framelen` - Длина системы координат
положительное нечетное целое число

`weights` - Вектор взвешивания
действительный положительный вектор

`b` - Изменяющиеся во времени коэффициенты конечной импульсной характеристики
матрица

`g` - Матрица дифференцировочных фильтров
матрица

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