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

Аналогичная задача была несколько лет назад. Только был датчик положения (расстояния).
В принципе, если задержка не важна, то приемлемые результаты дает фильтр кальмана.
Мы же остановились на smoothing spline. Но тут сложно подобрать параметры сглаживания.

И то и другое есть разновидность МНК.
И то и другое крайне неудобно реализовывать в дробной (целочисленной) арифметике.

Линейная аппроксимация некоторого набора данных - довольно просто реализуется.
Каждый раз суммировать заданное кол-во данных по точкам наблюдения необязательно, что сводит определение каждой суммы формулы расчета параметра a к простым матоперациям, где самое тяжелое - одно деление.
Если я где-то неправ - поправьте меня.

Если вам нужна только 1-ая производная, но простейшей функцией для сглаживания/усреднения будет линейный полином. Однако он не годится, если требуется оценка 2-ой производной, т.к. прямая линия повсюду имеет 2-ую производную равную нулю. Поэтому во втором случае требуется полином не ниже квадратичного (параболы).

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

Все коэффициенты полинома, вычисляемого с помощью МНК, можно свести к виду скалярного произведения вектора/массива коэффициентов (это другие коэффициенты - не полиномиальные, а весовые на каждую точку данных) на массив сырых данных в окрестности той точки, для которой вычисляется производная. Например, если среднее значение производной вам нужно вычислить в точке i, по 7-ми точкам (слева и справа по 3 точки), то это будут два массива по 7 элементов в каждом. В первом - предварительно вычисленные коэффициенты, а вот втором - сырые измерения от точки с индексом i-3, до точки с индксом i+3.

Вычисления 1-ой, 2-ой и старших производных отличаются ТОЛЬКО заменой весовых коффициентов! Причем их всегда можно посчитать заранее на персональном компьютере с помошью мат-пакетов. Перекладывать эту работу на МК не стоит. Нужно только определиться с тем, какую выбрать длину усредняемой области и степень полинома. Если у вас колесо крутится, и хочется усреднить до целого оборота, то длину отрезка вероятно лучше брать равной числу оцифрованных точек за полный оборот колеса. Но от этого вектора станут длиннее, и МК придется дольше считать сумму произведений.

Предположим, что мы пока остановились на 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 в виде:

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

А дальше можно и на МК считать, если эти весовые коэффициентики в него зашить. Коэффициенты полинома получатся такими
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, конечно, маловато.

Спасибо, суть расчета я если честно не понял (как и почему), со времен ВУЗа много времени прошло, но конечный результат ясен.
Мне нужна производная от скорости (ускорение), т.е. первая производная.
В Вашем описании нет упоминаний о наборе Х-координат (шаг непостоянный), он не нужен ни при каких условиях для данного метода?

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


У меня в качестве X-координаты используется индекс массива, т.к. предполагается, что измерения в нем следуют с постоянным временным интервалом.

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

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

Данная формула (вычисление a, b - ненужно), по моим подсчетам, будет эквивалентна 15-20 простым матоперациям:

Если не ошибаюсь, то это около 20мкс+-.
Судя по всему она вполне подходит для заданного набора точек (X;Y) - т.е. переменный шаг учтен.
Посмотрим что из этого выйдет если ускорение буду считать пусть раз 100 (или 200/500) в секунду...
Благодарю за участие!

Ваш оптимизм безоснователен - основные затраты времени происходят не там, где нужно вычесть/поделить друг на друга суммы, а чтобы те суммы набрать! Задайтесь реальным n и посмотрите, сколько у вас получится операций сложения+умножения при вычислении всех этих сумм.

Я указал количество операций с учетом накопления ЛЮБОГО количества наблюдений (N).
От N, при наличии некоторой области памяти, длительность расчетов НИКАК не зависит.

Откуда Вы взяли пары x,y?
У Вас измеряется лишь дельта времени между сигналами, а значит имеется одномерный поток данных dT[i].
Причем это уже формально есть угловая скорость вращения и для получения ускорения достаточно первой производной, ну и расстояние считается интегрированием ( тут я думаю, без проблем ).



Длительнось любых расчетов над элементами массива определяется длиной массива, если в расчет вовлечены все его элементы - это аксиома.

Например в фильтре скользящего среднего количество вычислений не зависит от длины фильтра

Не только скользящее среднее...

I32shr - побитовый сдвиг вправо на заданное кол-во разрядов с учетом знака.

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

>Не только скользящее среднее...
То, что Вы привели и есть скользящее среднее.

В общем же случае, взвешенная обработка массива ( нерекурсивный фильтр ) требует участия всех элементов.

P.S.
Вы не на то обратили внимание в моих словах
Я бы рекомендовал Вам поэкпериментировать с нерекурсивным фильтром выполняющим операцию дифференцирования.

Всегда думал что скользящее это:
Новое=Старое*(к-1)/к + Новое*к (с вышеизложенным это разные способы как по реализации, так и по итоговой точности).
Можно примером, если не затруднит.

http://ru.wikipedia.org/wiki/%D0%A1%D0%BA%...%BD%D1%8F%D1%8F

Вы правы. Тогда что за фильтр, который я назвал скользящим средним?

Один из вариантов рекурсивного. Это потому что он использует уже отфильтрованные (исправленные на прошлых шагах) значения.

Тест здесь довольно простой. Представьте себе, что ваше колесо остановилось , и с некоторого момента вы начинаете получать исключительно одни нулевые значения.
В этой ситуации нерекурсивный фильтр после n шагов (n - это предел суммирования тех самых сумм) полностью забудет предысторию. А на протяжении этих n шагов плавно опустит среднее до нуля. А ваш вариант рекурсивен, т.к. он не забудет "былое величие" никогда. Точнее говоря, оно будет рассасываться теоретически бесконечно долго (подобно тому как сходятся к нулю обратные величины натурального ряда).

Второй тест на ту же тему. Нерекурсивный фильтр дает одно и тоже значение среднего, независимо от того, двигается ли окно слева направо или справа налево. А для рекурсивных это не одно и тоже, т.е. слева и справа разная предыстория.

google: дифференцирующий цифровой фильтр

Например:
http://www.nbuv.gov.ua/portal/natural/Popu/2004_2/4-2.pdf

Или
Проектирование специализированных информационно-вычислительных систем
Высшая школа, 1984, под ред. проф. Ю.М. Смирнова

Например, коэф-ты дифф. фильтра 1-го порядка со сглаживанием

K[N,i] = 6/((N+1)*N) - 12*i/(N*(N-1)*(N+1))

В основе теории - синтез полиномиальных оптимальных по критерию СКО фильтров решением вариационной задачи минимизации целевой функции в условиях представления полезного сигнала непрерывной и дифференцируемой функцией (аппроксимация полиномом конечной степени), а помехи - некоррелированным белым шумом.

Всё верно, но мы подняли вопрос при оценке расчётоёмкости данным способом, где Вы его оценили как "не для МК"...
Рекурсивный? - Да (этого никто не отрицал). Тяжелый? - Нет. Кол-во затрат = всё тем же 15-20 операций на уровне +/-/* (одно деление так же оценено в виде эквивалента в простых операциях). Таких расчётов нужно сделать в 1000 раз меньше, чем способен МК!

Вам надо решение найти или поспорить тут?
Если спорить - Вы проспорите
Ну и решение заодно не найдете.

P.S.
Не, ну мне нравятся эти софтовые мальчики от сохи - вместо того, чтобы начать решать сверху, они заходят с заднего прохода. Вот и результат - в трех соснах.