Фильтр второго порядка Опции

[Yevdokimenko]

Добрый день.
Я не знаю в этом ли разделе должна быть эта тема, но всё же...
В ходе разработки устройства, одной из функций которого является определение ускорения, столкнулся с проблемой вычисления оного.
Итак, на автомобиле (на ступице колеса) установлен зубчатый венец, который вращается "вместе" с колесом.
Этот венец своими зубьям создает "наводки" на датчик, который "генерирует" сигнал с частотой, пропорциональной частоте вращения колеса (скорости движения).
Если венец и возможно изготовить в промышленных масштабах с хорошей точностью (но все же неидеально), то колесо как таковое не является абсолютно упругим телом, и имеет переменный радиус обката, что приводит к переменной частоте сигнала даже при условии равномерного движения (с постоянной скоростью).
Судя по проведенным тестам разброс частоты относительно некой средней подчиняется нормальному закону распределения случайной величины.
В случае же неравномерного (и даже не равноускоренного) движения разброс частоты становится еще большим.
Если для вычисления скорости (первой производной) такая ситуация вполне допустима, то при вычислении ускорения (второй производной) расчетный параметр ОЧЕНЬ зашумлен.
Посоветуйте, пожалуйста, с какой стороны подойти к решению данной задачи?
Пока что в голову пришла лишь мысль об аппроксимации МНК (методом наименьших квадратов) для нахождения линейной функции по последним N наблюдениям... Наклон данной функции и будет искомым ускорением.
Достаточно ли функции линейного вида, либо нужны функции больших порядков (квадратичные и т.п.)?

[Xenia]

Предположим, что мы пока остановились на N=7 (± 3 точки от центральной) и сглаживающем полиноме 2-ой степени (параболе). Тогда искомое уравнение для МНК в матричном виде будет иметь вид:
Y = X * C
где:
Y - семь сырых отсчетов в окрестности точки, в которой ищется производная (Y[i-3],Y[i-2],Y[i-1],Y[i],Y[i+1],Y[i+2],Y[i+3]).
С - коэффициенты квадратичного полинома y = с2*x^2 + с1*x + с0, напихнутые в массив/вектор C[3].
X - матрица степеней индексов.
Здесь вектор C - искомое, вектор Y - заданное, а матрицу X придется построить:
1-ый столбец X - все 7 индексов в нулевой степени, т.е. все единички:
1 1 1 1 1 1 1
2-ый столбец X - все 7 индексов в 1-ой степени, т.е. такие, как они есть:
-3 -2 -1 0 +1 +2 +3
(здесь можно нумеровать и в порядке "1 2 3 4 5 6 7" или "0 1 2 3 4 5 6", но я предпочитаю именно этот вариант, надеясь на симметрию в будущем, позволяющую уменьшить число умножений).
3-ый столбец X - все 7 индексов в 2-ой степени, т.е. квадраты предыдущего столбца:
9 4 1 0 1 4 9
Если у вас квадратичный полином, то это всё, а если был бы кубический, то пришлось бы еще добавить столбец из кубов.
Итого получаем матрицу X в виде:

Код

1    -3     9
1    -2     4
1    -1     1
1     0     0
1     1     1
1     2     4
1     3     9

А дальше делаем этой матрице операцию pinv. Если в матричной алгебре волокёте, то и без меня знаете, что это такое, а если нет, то и разбираться в том не будем, т.к. большой необходимости в том нет:

Код

pinv(X)
-0.0952    0.1429    0.2857    0.3333    0.2857    0.1429   -0.0952
-0.1071   -0.0714   -0.0357    0.0000    0.0357    0.0714    0.1071
+0.0595   -0.0000   -0.0357   -0.0476   -0.0357    0.0000    0.0595

А дальше можно и на МК считать, если эти весовые коэффициентики в него зашить. Коэффициенты полинома получатся такими
c0 = -0.0952*Y[i-3]+0.1429*Y[i-2]+0.2857*Y[i-1]+0.3333*Y[i]+0.2857*Y[i+1]+0.1429*Y[i+2]-0.0952*Y[i+3]
с1 = -0.1071*Y[i-3]-0.0714*Y[i-2]-0.0357*Y[i-1]+0.0000*Y[i]+0.0357*Y[i+1]+0.0714*Y[i+2]+0.1071*Y[i+3]
с2 = +0.0595*Y[i-3]-0.0000*Y[i-2]-0.0357*Y[i-1]-0.0476*Y[i]-0.0357*Y[i+1]+0.0000*Y[i+2]+0.0595*Y[i+3]
Т.е. получили все 3 коэффициента квадратичного полинома y = с2*x^2 + с1*x + с0

Впрочем, считать еще рано! Т.к. не все коэффициенты нам понадобятся.
У полинома
y = с2*x^2 + с1*x + с0
1-ая производная такая:
y' = 2*с2*x^2 + с1
А 2-ая производная такая:
y" = 2*с2
Отсюда следует, что для вычисления 2-ой производной коэффициенты полинома c1 и c0 не нужны, а потому и вычислять их не надо. А чтобы не множить каждый раз на 2, можно сразу умножить на 2 весовые коэффициенты третьей строки:
y"[i] = 2*С2 = +0.1190*Y[i-3]-0.0000*Y[i-2]-0.0714*Y[i-1]-0.0952*Y[i]-0.0714*Y[i+1]+0.0000*Y[i+2]+0.1190*Y[i+3]
- это и есть искомая 2-ая производная, вычисляемая для средней точки в ряду из 7-ми последовательных.

Но это еще не всё, т.к. теперь можно извлечь профит из симметрии в коэффициентах, вынося за общую скобку!
y"[i] = +0.1190*(Y[i-3]+Y[i+3])-0.0000*(Y[i-2]+Y[i+2])-0.0714*(Y[i-1]+Y[i+1])-0.0952*Y[i]
от этого умножений стало почти в два раза меньше.
Теперь выкидываем члены с нулевым коэффициентом:
y"[i] = +0.1190*(Y[i-3]+Y[i+3])-0.0714*(Y[i-1]+Y[i+1])-0.0952*Y[i]
Видите, как хорошо стало?

Но..., можно сделать еще лучше для МК, если заменить вычисления с плавающей точкой на целочисленные. Например, умножить все коэффициенты на круглую степень числа 2 (скажем на 256), округлить, в у результата отрасывать младший байт, компенсируя умножение. Приблизительно получим:
y"[i]*256 = 15*(Y[i-3]+Y[i+3])-9*(Y[i-1]+Y[i+1])-12*Y[i]
МК от такой постановки задачи, несомненно, возрадуется . Хотя 256, конечно, маловато.