Фильтр Савицкого – Голея

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

Фильтр Савицкого-Голея — это цифровой фильтр , который можно применять к набору точек цифровых данных с целью сглаживания данных, то есть для повышения точности данных без искажения тенденции сигнала. Это достигается в процессе, известном как свертка , путем аппроксимации последовательных подмножеств соседних точек данных полиномом низкой степени с помощью метода линейных наименьших квадратов . Когда точки данных расположены на одинаковом расстоянии друг от друга, можно найти аналитическое решение уравнений наименьших квадратов в форме одного набора «коэффициентов свертки», который можно применять ко всем подмножествам данных, чтобы дать оценки сглаженного сигнал (или производные сглаженного сигнала) в центральной точке каждого подмножества. Метод, основанный на установленных математических процедурах, был популяризирован Авраамом Савицким и Марселем Дж. Голеем , которые опубликовали таблицы коэффициентов свертки для различных полиномов и размеров подмножеств в 1964 году. Исправлены некоторые ошибки в таблицах. Метод был расширен для обработки 2- и 3-мерных данных.

Статья Савицкого и Голея является одной из наиболее широко цитируемых статей в журнале Analytical Chemistry. и классифицируется этим журналом как одна из «10 плодотворных статей», в которой говорится, что «можно утверждать, что зарождение аналитического прибора с компьютерным управлением можно проследить до этой статьи».

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

Приложения

Данные состоят из набора точек $\{x_{j}$, $y_{j}\};j=1,...,n$, где $x_{j}$ является независимой переменной и $y_{j}$ является наблюдаемым значением. Их лечат комплексом $m$ коэффициенты свертки, $C_{i}$, согласно выражению

$Y_{j}=\sum _{i={\tfrac {1-m}{2}}}^{\tfrac {m-1}{2}}C_{i}\,y_{j+i},\qquad {\frac {m+1}{2}}\leq j\leq n-{\frac {m-1}{2}}$

Выбранные коэффициенты свертки показаны в таблицах ниже . Например, для сглаживания 5-точечным квадратичным полиномом: $m=5,i=-2,-1,0,1,2$ и $j^{th}$ сглаженная точка данных, $Y_{j}$, определяется

$Y_{j}={\frac {1}{35}}(-3y_{j-2}+12y_{j-1}+17y_{j}+12y_{j+1}-3y_{j+2})$,

где, $C_{-2}=-3/35,C_{-1}=12/35$и т. д. Существует множество применений сглаживания, например, для предотвращения распространения шума по цепочке алгоритмов или иногда просто для того, чтобы данные выглядели менее зашумленными, чем они есть на самом деле.

Ниже приведены приложения численного дифференцирования данных. Примечание. При вычислении n- й производной применяется дополнительный коэффициент масштабирования ${\frac {n!}{h^{n}}}$ может применяться ко всем рассчитанным точкам данных для получения абсолютных значений (см. выражения для ${\frac {d^{n}Y}{dx^{n}}}$ниже, подробности).

(1) Синтетический лоренциан + шум (синий) и 1-я производная (зеленый)

(2) Кривая титрования (синяя) малоновой кислоты и 2-го производного (зеленая). Часть в голубой рамке увеличена в 10 раз.

(3) Лоренциан на экспоненциальной базовой линии (синий) и 2-й производной (зеленый)

(4) Сумма двух лоренцианов (синий) и 2-й производной (зеленый)

(5) 4-я производная суммы двух лоренцианов

1. Расположение максимумов и минимумов на кривых экспериментальных данных. Это была заявка, которая впервые вдохновила Савицкого. Первая производная функции равна нулю в максимуме или минимуме. На диаграмме показаны точки данных, принадлежащие синтетической кривой Лоренца , с добавленным шумом (синие ромбы). Данные нанесены в масштабе половинной ширины относительно максимума пика при нуле. Сглаженная кривая (красная линия) и первая производная (зеленая) были рассчитаны с помощью 7-точечных кубических фильтров Савицкого – Голея. Линейная интерполяция значений первой производной в положениях по обе стороны от перехода через нуль дает положение пикового максимума. Для этой цели также можно использовать третьи производные.
2. Расположение конечной точки на кривой титрования . Конечная точка — это точка перегиба , где вторая производная функции равна нулю. Кривая титрования малоновой кислоты иллюстрирует эффективность метода. Первая конечная точка при объеме 4 мл едва видна, но вторая производная позволяет легко определить ее значение путем линейной интерполяции для нахождения перехода через нуль.
3. Выравнивание базовой линии. В аналитической химии иногда необходимо измерить высоту полосы поглощения относительно изогнутой базовой линии. Поскольку кривизна базовой линии намного меньше кривизны полосы поглощения, вторая производная эффективно выравнивает базовую линию. Тремя мерами высоты производной, которая пропорциональна высоте полосы поглощения, являются расстояния «от пика до впадины» h1 и h2 и высота от базовой линии h3.
4. Повышение разрешения в спектроскопии. Полосы второй производной спектроскопической кривой уже, чем полосы в спектре: они имеют уменьшенную полуширину . Это позволяет частично перекрывающимся полосам «разлагаться» на отдельные (отрицательные) пики. Диаграмма показывает, как это можно использовать и для химического анализа , используя измерение расстояний «от пика до впадины». В этом случае долины являются свойством 2-й производной лоренциана. ( Положение оси X указано относительно положения максимума пика на шкале половины ширины и половины высоты ).
5. Повышение разрешения с помощью 4-й производной (положительные пики). Минимумы являются свойством 4-й производной лоренциана.

Скользящее среднее

Основная статья: Скользящее среднее

«Фильтр скользящего среднего» — это тривиальный пример фильтра Савицкого-Голея, который обычно используется с данными временных рядов для сглаживания краткосрочных колебаний и выделения долгосрочных тенденций или циклов. Каждое подмножество набора данных соответствует прямой горизонтальной линии, а не полиному более высокого порядка. Невзвешенный фильтр скользящего среднего — это самый простой фильтр свертки.

Скользящее среднее часто используется для быстрого технического анализа финансовых данных, таких как цены на акции, доходность или объемы торгов. Он также используется в экономике для изучения валового внутреннего продукта, занятости или других макроэкономических временных рядов.

Он не был включен в некоторые таблицы коэффициентов свертки Савицкого-Голея, поскольку все значения коэффициентов идентичны, причем значение ${\frac {1}{m}}$.

Вывод коэффициентов свертки

Когда точки данных расположены на одинаковом расстоянии друг от друга, можно найти аналитическое решение уравнений наименьших квадратов. Это решение лежит в основе сверточного метода численного сглаживания и дифференцирования. Предположим, что данные состоят из набора из n точек ( x _j , y _j ) ( j = 1, ..., n ), где x _j — независимая переменная, а y _j — исходное значение. Полином будет аппроксимироваться линейным методом наименьших квадратов набору из m (нечетное число) соседних точек данных, каждая из которых разделена интервалом h . Сначала производится замена переменной

$z={{x-{\bar {x}}} \over h}$

где ${\bar {x}}$ значение центральной точки. z принимает значения ${\tfrac {1-m}{2}},\cdots ,0,\cdots ,{\tfrac {m-1}{2}}$ (например, m = 5 → z = −2, −1, 0, 1, 2). Полином степени k определяется как

$Y=a_{0}+a_{1}z+a_{2}z^{2}\cdots +a_{k}z^{k}.$

Коэффициенты a ₀ , a ₁ и т. д. получаются путем решения нормальных уравнений (жирным шрифтом a обозначен вектор , жирным шрифтом J обозначена матрица ).

${\mathbf {a} }=\left({{\mathbf {J} }^{\mathbf {T} }{\mathbf {J} }}\right)^{-{\mathbf {1} }}{\mathbf {J} }^{\mathbf {T} }{\mathbf {y} },$

где $\mathbf {J}$ является матрицей Вандермонда , то есть $i$-й ряд $\mathbf {J}$ имеет ценности $1,z_{i},z_{i}^{2},\dots$.

Например, для кубического полинома, сопоставленного с 5 точками, z = −2, −1, 0, 1, 2, нормальные уравнения решаются следующим образом.

$\mathbf {J} ={\begin{pmatrix}1&-2&4&-8\\1&-1&1&-1\\1&0&0&0\\1&1&1&1\\1&2&4&8\end{pmatrix}}$

$\mathbf {J^{T}J} ={\begin{pmatrix}m&\sum z&\sum z^{2}&\sum z^{3}\\\sum z&\sum z^{2}&\sum z^{3}&\sum z^{4}\\\sum z^{2}&\sum z^{3}&\sum z^{4}&\sum z^{5}\\\sum z^{3}&\sum z^{4}&\sum z^{5}&\sum z^{6}\\\end{pmatrix}}={\begin{pmatrix}m&0&\sum z^{2}&0\\0&\sum z^{2}&0&\sum z^{4}\\\sum z^{2}&0&\sum z^{4}&0\\0&\sum z^{4}&0&\sum z^{6}\\\end{pmatrix}}={\begin{pmatrix}5&0&10&0\\0&10&0&34\\10&0&34&0\\0&34&0&130\\\end{pmatrix}}$

Теперь нормальные уравнения можно разделить на две отдельные группы уравнений, переставив строки и столбцы, при этом

$\mathbf {J^{T}J} _{\text{even}}={\begin{pmatrix}5&10\\10&34\\\end{pmatrix}}\quad \mathrm {and} \quad \mathbf {J^{T}J} _{\text{odd}}={\begin{pmatrix}10&34\\34&130\\\end{pmatrix}}$

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

$(\mathbf {J^{T}J} )_{\text{even}}^{-1}={1 \over 70}{\begin{pmatrix}34&-10\\-10&5\\\end{pmatrix}}\quad \mathrm {and} \quad (\mathbf {J^{T}J} )_{\text{odd}}^{-1}={1 \over 144}{\begin{pmatrix}130&-34\\-34&10\\\end{pmatrix}}$

Обычные уравнения становятся

${\begin{pmatrix}{a_{0}}\\{a_{2}}\\\end{pmatrix}}_{j}={1 \over 70}{\begin{pmatrix}34&-10\\-10&5\end{pmatrix}}{\begin{pmatrix}1&1&1&1&1\\4&1&0&1&4\\\end{pmatrix}}{\begin{pmatrix}y_{j-2}\\y_{j-1}\\y_{j}\\y_{j+1}\\y_{j+2}\end{pmatrix}}$

и

${\begin{pmatrix}a_{1}\\a_{3}\\\end{pmatrix}}_{j}={1 \over 144}{\begin{pmatrix}130&-34\\-34&10\\\end{pmatrix}}{\begin{pmatrix}-2&-1&0&1&2\\-8&-1&0&1&8\\\end{pmatrix}}{\begin{pmatrix}y_{j-2}\\y_{j-1}\\y_{j}\\y_{j+1}\\y_{j+2}\\\end{pmatrix}}$

Умножая и удаляя общие факторы,

${\begin{aligned}a_{0,j}&={1 \over 35}(-3y_{j-2}+12y_{j-1}+17y_{j}+12y_{j+1}-3y_{j+2})\\a_{1,j}&={1 \over 12}(y_{j-2}-8y_{j-1}+8y_{j+1}-y_{j+2})\\a_{2,j}&={1 \over 14}(2y_{j-2}-y_{j-1}-2y_{j}-y_{j+1}+2y_{j+2})\\a_{3,j}&={1 \over 12}(-y_{j-2}+2y_{j-1}-2y_{j+1}+y_{j+2})\end{aligned}}$

Коэффициенты при y в этих выражениях известны как коэффициенты свертки . Они являются элементами матрицы

$\mathbf {C=(J^{T}J)^{-1}J^{T}}$

В общем,

$(C\otimes y)_{j}\ =Y_{j}=\sum _{i=-{\tfrac {m-1}{2}}}^{\tfrac {m-1}{2}}C_{i}\,y_{j+i},\qquad {\frac {m+1}{2}}\leq j\leq n-{\frac {m-1}{2}}$

В матричной записи этот пример записывается как

${\begin{pmatrix}Y_{3}\\Y_{4}\\Y_{5}\\\vdots \end{pmatrix}}={1 \over 35}{\begin{pmatrix}-3&12&17&12&-3&0&0&\cdots \\0&-3&12&17&12&-3&0&\cdots \\0&0&-3&12&17&12&-3&\cdots \\\vdots &\vdots &\vdots &\vdots &\vdots &\vdots &\vdots &\ddots \end{pmatrix}}{\begin{pmatrix}y_{1}\\y_{2}\\y_{3}\\y_{4}\\y_{5}\\y_{6}\\y_{7}\\\vdots \end{pmatrix}}$

Таблицы коэффициентов свертки, рассчитанные таким же образом для m до 25, были опубликованы для сглаживающего фильтра Савицкого – Голея в 1964 г. Значение центральной точки z = 0 получается из одного набора коэффициентов: для _{й производной и т. д .} 1 - _{0 для сглаживания, 1} Численные производные получаются путем дифференцирования Y . Это означает, что производные рассчитываются для сглаженной кривой данных . Для кубического многочлена

${\begin{aligned}Y&=a_{0}+a_{1}z+a_{2}z^{2}+a_{3}z^{3}=a_{0}&{\text{ at }}z=0,x={\bar {x}}\\{\frac {dY}{dx}}&={\frac {1}{h}}\left({a_{1}+2a_{2}z+3a_{3}z^{2}}\right)={\frac {1}{h}}a_{1}&{\text{ at }}z=0,x={\bar {x}}\\{\frac {d^{2}Y}{dx^{2}}}&={\frac {1}{h^{2}}}\left({2a_{2}+6a_{3}z}\right)={\frac {2}{h^{2}}}a_{2}&{\text{ at }}z=0,x={\bar {x}}\\{\frac {d^{3}Y}{dx^{3}}}&={\frac {6}{h^{3}}}a_{3}\end{aligned}}$

В общем случае полиномы степени (0 и 1), (2 и 3), (4 и 5) и т. д. дают одинаковые коэффициенты сглаживания и четные производные. Полиномы степени (1 и 2), (3 и 4) и т. д. дают одинаковые коэффициенты при нечетных производных.

Алгебраические выражения

Не обязательно всегда использовать таблицы Савицкого–Голея. Суммы в матрице J ^Т J можно оценить в закрытой форме ,

${\begin{aligned}\sum _{z=-{\frac {m-1}{2}}}^{\frac {m-1}{2}}z^{2}&={m(m^{2}-1) \over 12}\\\sum z^{4}&={m(m^{2}-1)(3m^{2}-7) \over 240}\\\sum z^{6}&={m(m^{2}-1)(3m^{4}-18m^{2}+31) \over 1344}\end{aligned}}$

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

Сглаживание, полиномиальная степень 2,3: $C_{0i}={\frac {\left({3m^{2}-7-20i^{2}}\right)/4}{m\left({m^{2}-4}\right)/3}};\quad {\frac {1-m}{2}}\leq i\leq {\frac {m-1}{2}}$ (диапазон значений i также применим к выражениям ниже)

1-я производная: полиномиальная степень 3,4 $C_{1i}={\frac {5\left({3m^{4}-18m^{2}+31}\right)i-28\left({3m^{2}-7}\right)i^{3}}{m\left({m^{2}-1}\right)\left({3m^{4}-39m^{2}+108}\right)/15}}$

2-я производная: полиномиальная степень 2,3 $C_{2i}={\frac {12mi^{2}-m\left({m^{2}-1}\right)}{m^{2}\left({m^{2}-1}\right)\left({m^{2}-4}\right)/30}}$

3-я производная: полиномиальная степень 3,4 $C_{3i}={\frac {-\left({3m^{2}-7}\right)i+20i^{3}}{m\left({m^{2}-1}\right)\left({3m^{4}-39m^{2}+108}\right)/420}}$

Более простые выражения, которые можно использовать с кривыми, не имеющими точки перегиба:

Сглаживание, степень полинома 0,1 (скользящее среднее): $C_{0i}={\frac {1}{m}}$

1-я производная, полиномиальная степень 1,2: $C_{1i}={\frac {i}{m(m^{2}-1)/12}}$

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

Использование ортогональных полиномов

Альтернативой аппроксимации m точек данных простым полиномом во вспомогательной переменной z является использование ортогональных полиномов .

$Y=b_{0}P^{0}(z)+b_{1}P^{1}(z)\cdots +b_{k}P^{k}(z).$

где Р ⁰ , ..., П ^к представляет собой набор взаимно ортогональных многочленов степени 0,..., k . Полную информацию о том, как получить выражения для ортогональных многочленов и взаимосвязь между коэффициентами b и a, дает Гость. Выражения для коэффициентов свертки легко получить, поскольку матрица нормальных уравнений J ^Т J , является диагональной матрицей , поскольку произведение любых двух ортогональных многочленов равно нулю в силу их взаимной ортогональности. Следовательно, каждый ненулевой элемент обратного ему элемента является просто обратной величиной соответствующего элемента в матрице нормального уравнения. Вычисление еще больше упрощается за счет использования рекурсии для построения ортогональных полиномов Грама . Весь расчет можно закодировать в нескольких строках PASCAL — компьютерного языка, хорошо адаптированного для вычислений, включающих рекурсию.

Лечение первой и последней точки

Фильтры Савицкого – Голея чаще всего используются для получения сглаженного или производного значения в центральной точке z = 0 с использованием одного набора коэффициентов свертки. ( m − 1)/2 точки в начале и конце серии не могут быть вычислены с помощью этого процесса. Чтобы избежать этого неудобства, можно использовать различные стратегии.

• Данные можно было искусственно расширить, добавив в обратном порядке копии первых ( m − 1)/2 точек в начале и копии последних ( m − 1)/2 точек в конце. Например, при m = 5 две точки добавляются в начале и конце данных y ₁ , ..., y _n .

y ₃ , y ₂ , y ₁ , ... , y _n , y _{n -1} , y _{n -2} .

• Снова взглянув на подгоночный полином, становится очевидным, что данные можно вычислить для всех значений z , используя все наборы коэффициентов свертки для одного полинома a ₀ .. a _k .

Для кубического многочлена

${\begin{aligned}Y&=a_{0}+a_{1}z+a_{2}z^{2}+a_{3}z^{3}\\{\frac {dY}{dx}}&={\frac {1}{h}}(a_{1}+2a_{2}z+3a_{3}z^{2})\\{\frac {d^{2}Y}{dx^{2}}}&={\frac {1}{h^{2}}}(2a_{2}+6a_{3}z)\\{\frac {d^{3}Y}{dx^{3}}}&={\frac {6}{h^{3}}}a_{3}\end{aligned}}$

• Коэффициенты свертки для отсутствующих первой и последней точек также можно легко получить. Это также эквивалентно оснащению первых ( m + 1)/2 точек одним и тем же полиномом, и аналогично для последних точек.

Взвешивание данных

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

$U=\sum _{i}w_{i}(Y_{i}-y_{i})^{2}$

минимизация в процессе наименьших квадратов имеет единичные веса, w _i = 1. Когда веса не одинаковы, нормальные уравнения становятся

$\mathbf {a} =\left(\mathbf {J^{T}W} \mathbf {J} \right)^{-1}\mathbf {J^{T}W} \mathbf {y} \qquad W_{i,i}\neq 1$,

Если для всех подмножеств данных используется один и тот же набор диагональных весов, $W={\text{diag}}(w_{1},w_{2},...,w_{m})$, можно записать аналитическое решение нормальных уравнений. Например, для квадратичного многочлена

$\mathbf {J^{T}WJ} ={\begin{pmatrix}\sum w_{i}~~~&\sum w_{i}z_{i}&\sum w_{i}z_{i}^{2}\\\sum w_{i}z_{i}&\sum w_{i}z_{i}^{2}&\sum w_{i}z_{i}^{3}\\\sum w_{i}z_{i}^{2}&\sum w_{i}z_{i}^{3}&\sum w_{i}z_{i}^{4}\end{pmatrix}}$

Явное выражение для обратной этой матрицы можно получить с помощью правила Крамера . Затем набор коэффициентов свертки можно получить как

$\mathbf {C} =\left(\mathbf {J^{T}W} \mathbf {J} \right)^{-1}\mathbf {J^{T}W} .$

В качестве альтернативы коэффициенты C можно рассчитать в электронной таблице, используя встроенную процедуру обращения матрицы для получения обратной матрицы нормальных уравнений. Этот набор коэффициентов, однажды вычисленный и сохраненный, может использоваться во всех расчетах, в которых применяется одна и та же схема взвешивания. Для каждой схемы взвешивания необходим свой набор коэффициентов.

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

Коэффициенты двумерной свертки

Двумерное сглаживание и дифференцирование также можно применять к таблицам значений данных, например, значений интенсивности в фотографическом изображении, состоящем из прямоугольной сетки пикселей. Такая сетка называется ядром, а точки данных, составляющие ядро, называются узлами. Хитрость заключается в том, чтобы преобразовать прямоугольное ядро ​​в одну строку, просто упорядочив индексы узлов. В то время как одномерные коэффициенты фильтра находятся путем подгонки полинома во вспомогательной переменной z к набору из m точек данных, двумерные коэффициенты находятся путем подгонки полинома во вспомогательных переменных v и w к набору значений в узлы ядра m × n . Следующий пример для двумерного полинома общей степени 3, m = 7 и n = 5 иллюстрирует процесс, который аналогичен процессу для одномерного случая, описанному выше.

$v={\frac {x-{\bar {x}}}{h(x)}};w={\frac {y-{\bar {y}}}{h(y)}}$

$Y=a_{00}+a_{10}v+a_{01}w+a_{20}v^{2}+a_{11}vw+a_{02}w^{2}+a_{30}v^{3}+a_{21}v^{2}w+a_{12}vw^{2}+a_{03}w^{3}$

Прямоугольное ядро ​​из 35 значений данных, d ₁ − d _35.

v В	−3	−2	−1	0	1	2	3
−2	д 1	д 2	д 3	д 4	д 5	д 6	d 7
−1	д 8	д 9	д 10	д 11	д 12	д 13	д 14
0	д 15	д 16	д 17	д 18	d19	д 20	д 21
1	д 22	д 23	д 24	д 25	д 26	д 27	д 28
2	д 29	д 30	д 31	д 32	д 33	д 34	д 35

становится вектором, когда строки располагаются одна за другой.

d = ( d ₁ ... d ₃₅ ) ^Т

Якобиан имеет 10 столбцов, по одному для каждого из параметров a ₀₀ − a ₀₃ , и 35 строк, по одной для каждой пары значений v и w . Каждая строка имеет вид

$J_{\text{row}}=1\quad v\quad w\quad v^{2}\quad vw\quad w^{2}\quad v^{3}\quad v^{2}w\quad vw^{2}\quad w^{3}$

Коэффициенты свертки рассчитываются как

$\mathbf {C} =\left(\mathbf {J} ^{T}\mathbf {J} \right)^{-1}\mathbf {J} ^{T}$

Первая строка C содержит 35 коэффициентов свертки, которые можно умножить на 35 значений данных соответственно, чтобы получить полиномиальный коэффициент. $a_{00}$, которое представляет собой сглаженное значение в центральном узле ядра (т. е. в 18-м узле приведенной выше таблицы). Аналогично, другие строки C можно умножить на 35 значений, чтобы получить другие полиномиальные коэффициенты, которые, в свою очередь, можно использовать для получения сглаженных значений и различных сглаженных частных производных в разных узлах.

Никитас и Паппа-Луизи показали, что в зависимости от формата используемого полинома качество сглаживания может существенно различаться. Они рекомендуют использовать полином вида

$Y=\sum _{i=0}^{p}\sum _{j=0}^{q}a_{ij}v^{i}w^{j}$

поскольку такие полиномы позволяют добиться хорошего сглаживания как в центральной, так и в приграничных областях ядра, и поэтому их можно уверенно использовать при сглаживании как во внутренних, так и в приграничных точках данных выборочной области. Чтобы избежать плохой обусловленности при решении задачи наименьших квадратов, p < m и q < n . Информацию о программном обеспечении, которое рассчитывает двумерные коэффициенты, и о базе данных таких C см. в разделе, посвященном многомерным коэффициентам свертки, ниже.

Многомерные коэффициенты свертки

Идея двумерных коэффициентов свертки может быть просто распространена и на более высокие пространственные измерения: путем организации многомерного распределения узлов ядра в один ряд. После вышеупомянутого открытия Никитаса и Паппа-Луизи в двумерных случаях в многомерных случаях рекомендуется использовать следующий вид полинома:

$Y=\sum _{i_{1}=0}^{p_{1}}\sum _{i_{2}=0}^{p_{2}}\cdots \sum _{i_{D}=0}^{p_{D}}(a_{i_{1}i_{2}\cdots i_{D}}\times u_{1}^{i_{1}}u_{2}^{i_{2}}\cdots u_{D}^{i_{D}})$

где D — размерность пространства, $a$'s - полиномиальные коэффициенты, а u 's - координаты в различных пространственных направлениях. Алгебраические выражения для частных производных любого порядка, смешанного или другого, можно легко вывести из приведенного выше выражения. Обратите внимание, что C зависит от способа расположения узлов ядра в ряду и от способа расположения различных членов расширенной формы вышеуказанного многочлена при построении якобиана.

Точное вычисление C в многомерных случаях становится сложной задачей, поскольку точность стандартных чисел с плавающей запятой, доступных в языках компьютерного программирования, больше не остается достаточной. Недостаточная точность приводит к тому, что ошибки усечения с плавающей запятой становятся сравнимыми с величинами некоторых элементов C , что, в свою очередь, серьезно ухудшает его точность и делает его бесполезным. Чандра Шекхар представил два программного обеспечения с открытым исходным кодом : Advanced Convolution Coefficient Calculator (ACCC) и Precision Convolution Coefficient Calculator (PCCC) , которые адекватно справляются с этими проблемами точности. ACCC выполняет вычисления, используя числа с плавающей запятой, итеративным образом. Точность чисел с плавающей запятой постепенно увеличивается на каждой итерации с помощью GNU MPFR . Как только полученные C в двух последовательных итерациях начинают иметь одинаковые значащие цифры до достижения заранее заданного расстояния, считается, что сходимость достигнута. Если расстояние достаточно велико, вычисление дает высокоточное значение C . PCCC использует вычисления рациональных чисел с использованием арифметической библиотеки множественной точности GNU и выдает полностью точный C в формате рациональных чисел . В конце концов, эти рациональные числа преобразуются в числа с плавающей запятой, до достижения заранее заданного количества значащих цифр.

Доступна база данных C , рассчитанных с использованием ACCC для симметричных ядер и как симметричных, так и асимметричных полиномов на узлах ядра с единичным интервалом в 1-, 2-, 3- и 4-мерном пространствах. Чандра Шекхар также изложил математическую основу, которая описывает использование C, рассчитанного на узлах ядра с единичным интервалом, для выполнения фильтрации и частичного дифференцирования (различного порядка) на узлах ядра с неравномерным расположением. позволяющее использовать C, представленный в вышеупомянутой базе данных. Хотя этот метод дает лишь приблизительные результаты, они приемлемы в большинстве инженерных приложений при условии, что неравномерность узлов ядра является слабой.

Некоторые свойства свертки

1. Сумма коэффициентов свертки для сглаживания равна единице. Сумма коэффициентов при нечетных производных равна нулю.
2. Сумма квадратов коэффициентов свертки для сглаживания равна значению центрального коэффициента.
3. Сглаживание функции оставляет область под функцией неизменной.
4. Свертка симметричной функции с коэффициентами четной производной сохраняет центр симметрии.
5. Свойства производных фильтров.

Искажение сигнала и снижение шума

Неизбежно, что сигнал будет искажен в процессе свертки. Согласно свойству 3, приведенному выше, при сглаживании данных, имеющих пик, высота пика будет уменьшена, а полуширина увеличена . Как степень искажения, так и улучшение отношения сигнал/шум ( отношение сигнал/шум ):

• уменьшаются по мере увеличения степени полинома
• увеличивается по мере увеличения ширины m функции свертки

Влияние сглаживания на точки данных с некоррелированным шумом единичного стандартного отклонения

Например, если шум во всех точках данных некоррелирован и имеет постоянное отклонение стандартное σ , стандартное отклонение шума будет уменьшено путем свертки с функцией сглаживания m -точки до

полиномиальная степень 0 или 1: ${\sqrt {1 \over m}}\sigma$ ( скользящее среднее )

полином 2 или 3 степени: ${\sqrt {\frac {3(3m^{2}-7)}{4m(m^{2}-4)}}}\sigma$.

Эти функции показаны на графике справа. Например, с помощью 9-точечной линейной функции (скользящее среднее) удаляется две трети шума, а с помощью 9-точечной квадратичной/кубической функции сглаживания удаляется только около половины шума. Большая часть оставшегося шума — это низкочастотный шум (см. Частотные характеристики сверточных фильтров ниже).

Хотя функция скользящего среднего обеспечивает лучшее подавление шума, она не подходит для сглаживания данных, кривизна которых превышает m точек. Функция квадратичного фильтра не подходит для получения производной кривой данных с точкой перегиба , поскольку у квадратичного многочлена ее нет. Оптимальный выбор порядка полинома и количества коэффициентов свертки будет компромиссом между шумоподавлением и искажениями.

Многопроходные фильтры

Один из способов смягчить искажения и улучшить удаление шума — использовать фильтр меньшей ширины и выполнить с ним более одной свертки. Для двух проходов одного и того же фильтра это эквивалентно одному проходу фильтра, полученного сверткой исходного фильтра с самим собой. Например, 2 прохода фильтра с коэффициентами (1/3, 1/3, 1/3) эквивалентны 1 проходу фильтра с коэффициентами (1/9, 2/9, 3/9, 2/9, 1/9).

Недостатком многопроходного режима является то, что эквивалентная ширина фильтра для $n$ пропуски $m$–точечная функция $n(m-1)+1$ поэтому многопроходная передача приводит к более серьезным конечным эффектам. Тем не менее, многопроходная передача используется с большим преимуществом. Например, около 40–80 проходов данных с соотношением сигнал/шум всего 5 дали полезные результаты. Приведенные выше формулы снижения шума не применимы, поскольку корреляция между рассчитанными точками данных увеличивается с каждым проходом.

Частотные характеристики сверточных фильтров

Преобразование Фурье 9-точечной квадратичной/кубической функции сглаживания

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

$FT(\theta )=\sum _{j={\tfrac {1-m}{2}}}^{\tfrac {m-1}{2}}C_{j}\cos(j\theta )$

θ изменяется от 0 до 180 градусов , после чего функция просто повторяется. График 9-точечной квадратичной/кубической функции сглаживания является типичным. Под очень малым углом график почти плоский, а это означает, что низкочастотные компоненты данных практически не изменятся в результате операции сглаживания. По мере увеличения угла значение уменьшается, так что более высокочастотные компоненты все больше и больше ослабляются. Это показывает, что сверточный фильтр можно описать как фильтр нижних частот : удаляемый шум представляет собой в основном высокочастотный шум, а низкочастотный шум проходит через фильтр. Некоторые высокочастотные компоненты шума ослабляются сильнее, чем другие, о чем свидетельствуют волнистости преобразования Фурье под большими углами. Это может привести к небольшим колебаниям в сглаженных данных. и обращение фазы, т. е. высокочастотные колебания в данных инвертируются с помощью фильтрации Савицкого – Голея.

Свертка и корреляция

Свертка влияет на корреляцию между ошибками в данных. Эффект свертки можно выразить как линейное преобразование.

$Y_{j}\ =\sum _{i=-{\frac {m-1}{2}}}^{\frac {m-1}{2}}C_{i}\,y_{j+i}$

По закону распространения ошибок дисперсионно -ковариационная матрица данных A будет преобразована в B согласно

$\mathbf {B} =\mathbf {C} \mathbf {A} \mathbf {C} ^{T}$

Чтобы увидеть, как это применяется на практике, рассмотрим влияние 3-точечного скользящего среднего на первые три расчетные точки, Y ₂ − Y ₄ , предполагая, что точки данных имеют одинаковую дисперсию и что между ними нет корреляции. A будет единичной матрицей , умноженной на константу σ ² , дисперсия в каждой точке.

$\mathbf {B} ={\sigma ^{2} \over 9}{\begin{pmatrix}1&1&1&0&0\\0&1&1&1&0\\0&0&1&1&1\end{pmatrix}}{\begin{pmatrix}1&0&0&0&0\\0&1&0&0&0\\0&0&1&0&0\\0&0&0&1&0\\0&0&0&0&1\end{pmatrix}}{\begin{pmatrix}1&0&0\\1&1&0\\1&1&1\\0&1&1\\0&0&1\end{pmatrix}}={\sigma ^{2} \over 9}{\begin{pmatrix}3&2&1\\2&3&2\\1&2&3\end{pmatrix}}$

В этом случае коэффициенты корреляции ,

$\rho _{ij}={\frac {B_{ij}}{\sqrt {B_{ii}B_{jj}}}},\ (i\neq j)$

между расчетными точками i и j будет

$\rho _{i,i+1}={2 \over 3}=0.66,\,\rho _{i,i+2}={1 \over 3}=0.33$

В общем, расчетные значения коррелируют, даже если наблюдаемые значения не коррелируют. Корреляция распространяется на m − 1 расчетную точку одновременно.

Многопроходные фильтры

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

$\mathbf {CBC^{T}} ={\sigma ^{2} \over 81}{\begin{pmatrix}1&1&1&0&0\\0&1&1&1&0\\0&0&1&1&1\end{pmatrix}}{\begin{pmatrix}3&2&1&0&0\\2&3&2&0&0\\1&2&3&0&0\\0&0&0&0&0\\0&0&0&0&0\end{pmatrix}}{\begin{pmatrix}1&0&0\\1&1&0\\1&1&1\\0&1&1\\0&0&1\end{pmatrix}}={\sigma ^{2} \over 81}{\begin{pmatrix}19&16&10&4&1\\16&19&16&10&4\\10&16&19&16&10\\4&10&16&19&16\\1&4&10&16&19\end{pmatrix}}$

После двух проходов стандартное отклонение центральной точки уменьшилось до ${\sqrt {\tfrac {19}{81}}}\sigma =0.48\sigma$, по сравнению с 0,58 σ за один проход. Снижение шума немного меньше, чем было бы получено при одном проходе 5-точечного скользящего среднего, что при тех же условиях привело бы к тому, что сглаженные точки имели бы меньшее стандартное отклонение 0,45 σ .

Корреляция теперь распространяется на диапазон из 4 последовательных точек с коэффициентами корреляции.

$\rho _{i,i+1}={16 \over 19}=0.84,\rho _{i,i+2}={10 \over 19}=0.53,\rho _{i,i+3}={4 \over 19}=0.21,\rho _{i,i+4}={1 \over 19}=0.05$

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

Сравнение с другими фильтрами и альтернативами

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

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

См. также

• Сглаживание ядра . Разная терминология для многих одних и тех же процессов, используемых в статистике.
• Локальная регрессия — методы LOESS и LOWESS.
• Числовое дифференцирование - Применение к дифференцированию функций.
• Сглаживающий сплайн
• Трафарет (численный анализ) – Применение к решению дифференциальных уравнений
• Ходрик-Прескотт_filter
• Калман_фильтр

Приложение

Таблицы выбранных коэффициентов свертки

Рассмотрим набор точек данных ⁠ $(x_{j},y_{j})_{1\leq j\leq n}$⁠ . Таблицы Савицкого – Голея относятся к случаю, когда шаг ⁠ $x_{j}-x_{j-1}$⁠ постоянна, ч . примеры использования так называемых коэффициентов свертки с кубическим полиномом и размером окна m Ниже приведены , равным 5 точкам.

Сглаживание

$Y_{j}={\frac {1}{35}}(-3\times y_{j-2}+12\times y_{j-1}+17\times y_{j}+12\times y_{j+1}-3\times y_{j+2})$ ;

1-я производная

$Y'_{j}={\frac {1}{12h}}(1\times y_{j-2}-8\times y_{j-1}+0\times y_{j}+8\times y_{j+1}-1\times y_{j+2})$ ;

2-я производная

$Y''_{j}={\frac {1}{7h^{2}}}(2\times y_{j-2}-1\times y_{j-1}-2\times y_{j}-1\times y_{j+1}+2\times y_{j+2})$.

Выбранные значения коэффициентов свертки для полиномов степени 1, 2, 3, 4 и 5 приведены в следующих таблицах. Значения были рассчитаны с использованием кода PASCAL, предоставленного в Gorry.

Коэффициенты сглаживания Полиномиальный Степень квадратичный или кубический 2 или 3 квартика или квинтика 4 или 5 Размер окна 5 7 9 7 9 −4 −21 15 −3 −2 14 5 −55 −2 −3 3 39 −30 30 −1 12 6 54 75 135 0 17 7 59 131 179 1 12 6 54 75 135 2 −3 3 39 −30 30 3 −2 14 5 −55 4 −21 15 Нормализация 35 21 231 231 429

Коэффициенты для 1-й производной Полиномиальный Степень линейный или квадратичный 1 или 2 кубический или четвертичный 3 или 4 Размер окна 3 5 7 9 5 7 9 −4 −4 86 −3 −3 −3 22 −142 −2 −2 −2 −2 1 −67 −193 −1 -1 −1 −1 −1 −8 −58 −126 0 0 0 0 0 0 0 0 1 1 1 1 1 8 58 126 2 2 2 2 −1 67 193 3 3 3 −22 142 4 4 −86 Нормализация 2 10 28 60 12 252 1,188

Коэффициенты для 2-й производной Полиномиальный Степень квадратичный или кубический 2 или 3 квартика или квинтика 4 или 5 Размер окна 5 7 9 5 7 9 −4 28 −126 −3 5 7 −13 371 −2 2 0 −8 −1 67 151 −1 −1 −3 −17 16 −19 −211 0 −2 −4 −20 −30 −70 −370 1 −1 −3 −17 16 −19 −211 2 2 0 −8 −1 67 151 3 5 7 −13 371 4 28 −126 Нормализация 7 42 462 12 132 1716

Коэффициенты для 3-й производной Полиномиальный Степень кубический или четвертичный 3 или 4 квинтика или секстик 5 или 6 Размер окна 5 7 9 7 9 −4 −14 100 −3 −1 7 1 −457 −2 −1 1 13 −8 256 −1 2 1 9 13 459 0 0 0 0 0 0 1 −2 −1 −9 −13 −459 2 1 −1 −13 8 −256 3 1 −7 −1 457 4 14 −100 Нормализация 2 6 198 8 1144

Коэффициенты для 4-й производной Полиномиальный Степень квартика или квинтика 4 или 5 Размер окна 7 9 −4 14 −3 3 −21 −2 −7 −11 −1 1 9 0 6 18 1 1 9 2 -7 −11 3 3 −21 4 14 Нормализация 11 143

Ссылки

• Ганс, Питер (1992). Аппроксимация данных в химических науках: методом наименьших квадратов . ISBN  9780471934127 .

Внешние ссылки

• Расширенный калькулятор коэффициента свертки (ACCC) для многомерных фильтров наименьших квадратов
• Фильтр Савицкого – Голея в основах статистики
• Более широкий диапазон коэффициентов для различных размеров наборов данных, порядков соответствия и смещений от центральной точки.