Думай, когда используешь готовые формулы

Подходит ко мне коллега и говорит: "У меня почему-то шум в данных, вычисленных мной, раз в 10 больше чем в оригинальных данных от прибора. Наверное, в приборе используется какое-то хитрое сглаживание, так как обычное усреднение не помогает. Не знаешь, что здесь можно применить?"

Пошли, посмотрели. Выглядело это примерно так:


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

Шум действительно большой. На какое-то хитрое сглаживание совершенно не похоже. Что бы такое сгладить надо сглаживать по несколько десяткам точкам. Надо искать в чем дело. На рисунке простейшее скользящее усреднение по 7 точкам.

Ну что же надо разбираться.

А что при других значениях угла α? Выяснялось, что проблема существует только при углах близких к 0 или 180. Ну вот и зацепка.

По оси X - установленный угол прибора, по оси Y - результат измерений. Вверху (синий) - сигнал с прибора, внизу   (красный) - расчетный сигнал,

 

Введение.

Предупреждение. Реальная ситуация была сильно упрощена, чтобы сделать это сообщение не очень объемным и понятным. Реальный прибор был полноценный 3-х мерный инклинометр и магнитометр. Здесь мы ограничимся 2-х мерным случаем и сделаем многие другие упрощения. Все иллюстрации были искусственно сгенерированы и не являются реальными данными с прибора. Но суть явления не затронута и все выводы соответствуют действительно произошедшему.

Существует старый прибор, который в числе прочего измерял угол своего отклонения от вертикали α. В качестве выходных данных он выдает сырые данные - проекции вектора ускорения свободного падения \(\vec g\) на оси прибора - (X, Y) и конечный результат, вычисляемый по этим данным - угол  α. Нам нужно было модернизировать  работу прибора, учесть некоторые неточности в механической конструкции, и улучшить точность показаний прибора. Причем одним из условий была совместимость выдаваемых результатов со старым прибором. Поэтому нужно было взять сырые данные с прибора, сделать вычисления по старой схеме и сравнить с результатами выдаваемым прибором. Только потом мы могли вносить свои улучшения.

Начинаем разбираться

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

 \(\cos(\alpha) = Y\)  и соответственно для вычисления он использовал  \(\alpha = \arccos(Y)\)

Коллега встретился с некоторыми трудностями при использовании этой формулы - например при некоторых данных величина Y становилась немного больше 1. Стандартная функция arccos естественно отказывалась работать в "таких условиях. Он вышел из положения довольно просто - он написал специальную функцию, которая при таких "неправильных" значениях Y выдавал ответ 0, а в "нормальных" вычисляла стандартный arccos.

Конечно более правильная формула - \(\cos(\alpha) =Y/\sqrt(X^2 + Y^2)\), но как известно ускорение свободного падения для практических целей у нас не меняется со временем, да и в различных местах на Земле величина g тоже не сильно меняется. Поэтому авторы документации для упрощения формул положили \(g=1\). Да и прибор тоже следовал этому соглашению: его показания величины \(X\) и \(Y\) были нормированы на величину g в какой-то точки Земли.

Объяснение

Глядя на формулу, которая использовалась при вычислениях, становится понятно, почему неприятности были при малых углах. Пусть при каком-то малом угле \(\alpha_0\) выполняется условие:

 \(\cos(\alpha_0) = Y_0,\ при\ этом\ Y_0\ близко\ к\ 1\)

используя приближенную формулу \(\cos(x) \approx 1 - x^2/2\) получаем

\(1 - \alpha_0^2/2  \approx Y_0\)

Пусть у нас есть измерение Y с небольшим шумом - \(Y_s = Y_0 + \delta Y\), тогда для вычисленного по этим данным угла мы имеем:

 \(\cos(\alpha_s) = Y_s\)

 и положив  \(\alpha_s = \alpha_0 + \delta \alpha\) ( \(\alpha_0 \ll 1,  \delta \alpha \ll 1\)) получаем:

\(1 - \alpha_s^2/2 \approx 1 - \alpha_0^2/2 - \alpha_0*\delta\alpha = Y_s = Y_0 + \delta Y\)

\(-\alpha_0*\delta\alpha \approx  \delta Y\)

\(\delta\alpha \approx -\delta Y/\alpha_0\)

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

Happy end

Исправить ситуацию оказалось очень просто - надо использовать правильную формулу:

\(\cos(\alpha) =Y/\sqrt(X^2 + Y^2)\)

Ну а еще более правильно при всяком вычислении углов по обратным тригонометрическим формулам использовать специальную функцию - atan2. Она выдает угол в полном диапазоне \(2\pi\) в отличии от \(\arccos и \arcsin\) и не имеет особых точек, где в принципе может теряется точность вычислений.

Как только формула вычисления угла была заменена на    \(\alpha = atan2(X, Y)\) все стало на свои места. Шумы при вычислений углов сравнялись с шумами углов, которые выдавал прибор.

Мораль

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