- Описание алгоритма
- Автоматическая калибровка магнитометра
- Акселерометр
- Калибровка исходных/необработанных данных (raw data)
- Код для pico
- Код для processing
- Магнитометр
- Общие сведения о датчике mpu9250
- Ориентация объекта в пространстве
- Получение “реальных” значений
- Преобразователь питания rt9193-18
- Преобразователь питания rt9193-33
- Работа с arduino и mpu6050 | alexgyver
- Расчет углов ориентации
- Считывание данных с датчика mpu9250 в плату arduino mega 2560
- Углы ориентации
- Фильтрация данных
Описание алгоритма
- Сформируем массивы случайных углов Эйлера roll, pitch, yaw. Они будут задавать наборы вариантов истинной ориентации объекта в модели.
- Из случайных углов roll, pitch, yaw формируется матрица преобразования из ССК XYZ в ЛСК ENU:
где
,
,
,
,
,
.
- Используя данную матрицу можно получить выражение для истинных ускорений в ССК:
— вектор, определяющий направление гравитационного ускорения, выраженный в единицах g,
— матрица преобразования координат из ЛСК в ССК (обратная матрице преобразования из ССК в ЛСК).
- Применяем модель измерения акселерометра:
- По измерениям акселерометра рассчитываются новые углы крена и тангажа (оценки) по формулам:
- Также необходимо сформировать матрицу пересчета в «горизонт» из этих углов, для этого воспользуемся функцией rpy2mat:
где углы roll’ и pitch’ — это углы, рассчитанные по измерениям акселерометра, а третий угол — нулевой.
- Возьмем вектор истинных магнитных плотностей в ЛСК ENU и пересчитаем его в ССК XYZ:
- Применяем модель измерений магнитометра:
- Осталось пересчитать измерения магнитометра из ССК в «горизонт»:
- По «горизонтированным» измерениям магнитометра расcчитывается угол магнитного азимута:
- Ошибки оценивания углов ориентации рассчитываются как разность между истинными углами roll, pitch, yaw и рассчитанными по измерениям датчиков — roll’, pitch’, yaw’.
Автоматическая калибровка магнитометра
Это одна из самых простых, но в то же время важных участков кода программы. Функция magcalMPU9250(float * dest1, float * dest2) производит калибровку магнитометра в то время когда вы перемещаете его по фигуре восьмерки (цифры “8”). Она сохраняет максимальные и минимальные значения, а затем на их основе вычисляет среднее.
Serial.println(“Mag Calibration: Wave device in a figure eight until done!”);
sample_count = 128;
for(ii = 0; ii < sample_count; ii ) {
readMagData(mag_temp); // Read the mag data
for (int jj = 0; jj < 3; jj ) {
if(mag_temp[jj] > mag_max[jj]) mag_max[jj] = mag_temp[jj];
if(mag_temp[jj] < mag_min[jj]) mag_min[jj] = mag_temp[jj];
}
delay(135); // at 8 Hz ODR, new mag data is available every 125 ms
}
// Get hard iron correction
mag_bias[0] = (mag_max[0] mag_min[0])/2; // get average x mag bias in counts
mag_bias[1] = (mag_max[1] mag_min[1])/2; // get average y mag bias in counts
mag_bias[2] = (mag_max[2] mag_min[2])/2; // get average z mag bias in counts
dest1[0] = (float) mag_bias[0]*mRes*magCalibration[0]; // save mag biases in G for main program
dest1[1] = (float) mag_bias[1]*mRes*magCalibration[1];
dest1[2] = (float) mag_bias[2]*mRes*magCalibration[2];
// Get soft iron correction estimate
mag_scale[0] = (mag_max[0] – mag_min[0])/2; // get average x axis max chord length in counts
mag_scale[1] = (mag_max[1] – mag_min[1])/2; // get average y axis max chord length in counts
mag_scale[2] = (mag_max[2] – mag_min[2])/2; // get average z axis max chord length in counts
float avg_rad = mag_scale[0] mag_scale[1] mag_scale[2];
avg_rad /= 3.0;
dest2[0] = avg_rad/((float)mag_scale[0]);
dest2[1] = avg_rad/((float)mag_scale[1]);
dest2[2] = avg_rad/((float)mag_scale[2]);
Serial.println(“Mag Calibration done!”);
}
voidmagcalMPU9250(float*dest1,float*dest2){ uint16_tii=,sample_count=; int32_tmag_bias[3]={,,},mag_scale[3]={,,}; int16_tmag_max[3]={0x8000,0x8000,0x8000},mag_min[3]={0x7FFF,0x7FFF,0x7FFF},mag_temp[3]={,,}; Serial.println(“MagCalibration:Wavedeviceinafigureeightuntildone!”); sample_count=128; for(ii=;ii<sample_count;ii ){ readMagData(mag_temp); // Read the mag data for(intjj=;jj<3;jj ){ if(mag_temp[jj]>mag_max[jj])mag_max[jj]=mag_temp[jj]; if(mag_temp[jj]<mag_min[jj])mag_min[jj]=mag_temp[jj]; } delay(135); // at 8 Hz ODR, new mag data is available every 125 ms } // Get hard iron correction mag_bias[] =(mag_max[] mag_min[])/2; // get average x mag bias in counts mag_bias[1] =(mag_max[1] mag_min[1])/2; // get average y mag bias in counts mag_bias[2] =(mag_max[2] mag_min[2])/2; // get average z mag bias in counts dest1[]=(float)mag_bias[]*mRes*magCalibration[]; // save mag biases in G for main program dest1[1]=(float)mag_bias[1]*mRes*magCalibration[1]; dest1[2]=(float)mag_bias[2]*mRes*magCalibration[2]; // Get soft iron correction estimate mag_scale[] =(mag_max[]–mag_min[])/2; // get average x axis max chord length in counts mag_scale[1] =(mag_max[1]–mag_min[1])/2; // get average y axis max chord length in counts mag_scale[2] =(mag_max[2]–mag_min[2])/2; // get average z axis max chord length in counts floatavg_rad=mag_scale[] mag_scale[1] mag_scale[2]; avg_rad/=3.0; dest2[]=avg_rad/((float)mag_scale[]); dest2[1]=avg_rad/((float)mag_scale[1]); dest2[2]=avg_rad/((float)mag_scale[2]); Serial.println(“MagCalibrationdone!”); } |
Акселерометр
Это датчик, измеряющий проекцию кажущегося ускорения на ось чувствительности. Кажущегося — потому что измеряет и силу тяжести тоже, даже в то время как акселерометр неподвижен. Проще всего представить акселерометр как грузик на пружинке, его выдаваемые измерения пропорциональны степени растяжения пружины. Если акселерометр покоится — пружина растянута лишь силой тяжести. Если ускоряется — то будет сумма сил: инерции грузика
и силы тяжести
Примем следующую модель измерений триады ортогональных (взаимно перпендикулярных) акселерометров:
где
– измеряемое ускорение в ССК (собственной системе координат) XYZ,
– матрица перекоса осей и масштабных коэффициентов акселерометра,
– вектор истинного ускорения в ССК XYZ,
– вектор смещения нуля акселерометра,
– шум измерений.
Матрица перекоса осей и масштабных коэффициентов выглядит следующим образом:
где элементы, расположенные по главной диагонали (
) — это масштабные коэффициенты и их погрешности по трём осям акселерометра, а остальные элементы матрицы — перекосы осей акселерометра.
Калибровка исходных/необработанных данных (raw data)
Получаемые с датчика MPU9250 данные должны быть откалиброваны в соответствии с оборудованием пользователя, в котором будет использоваться датчик. Калибровка магнитометра требуется для компенсации магнитного склонения (Magnetic Declination). Корректировочное значение в данном случае сильно зависит от местоположения датчика. Необходимо откалибровать две переменные: yaw и magbias.
Ниже приведенный фрагмент кода реализует калибровку переменной yaw для магнитного склонения, характерного для местности: Potheri, Chennai, India.
Данные о магнитном склонении в конкретной местности можно получить, к примеру, с помощью следующих сайтов:
В следующем фрагменте кода реализована калибровка переменной magbias, которая взята из функции magcalMPU9250(float * dest1, float * dest2) (из библиотеки, которую мы используем для работы с датчиком MPU9250).
Код для pico
- imu_visualization.py
# Импортируем необходимые модулиfrom amperka_pico_imu_10dof import icm20948 importtimeimportmath # Печатаем тестовую строку в консольprint("IMU-sensor Test Program")# Создаём объект для работы с IMU-сенсором imu=icm20948.ICM20948() whileTrue: # Считываем данных с акселерометра accelerometer = imu.read_accelerometer_g()# Считываем данных с гироскопа gyroscope = imu.read_gyroscope_rad()# Считываем данных с магнитометра magnetometer = imu.read_magnetometer_() ax = accelerometer[] ay = accelerometer[1] az = accelerometer[2] gx = gyroscope[] gy = gyroscope[1] gz = gyroscope[2] mx = magnetometer[] my = magnetometer[1] mz = magnetometer[2]# Обновляем входные данные в фильтр Маджвика imu.AHRS_update(gx, gy, gz, ax, ay, az, mx, my, mz)# Считываем кватернион quaternion = imu.read_quaternion()# Выводим кватернион в консоль для передачи в Processingprint(f"{quaternion[0]},{quaternion[1]},{quaternion[2]},{quaternion[3]}")time.sleep(0.1)
Код для processing
- imu_visualization_draw.pde
importprocessing.serial.*;importtoxi.geom.*;importtoxi.processing.*; // NOTE: requires ToxicLibs to be installed in order to run properly.// 1. Download from http://toxiclibs.org/downloads// 2. Extract into [userdir]/Processing/libraries// (location may be different on Mac/Linux)// 3. Run and bask in awesomeness // The serial port Serial port; String message; float[] q =newfloat[4]; Quaternion quat =new Quaternion(1, , , );// New line character in ASCIIfinalchar newLine ='n';String[] massQ =newString[4];float[] ypr =newfloat[3]; void setup(){// Size form 400x400 size(400, 400, P3D);// Open serial port// Replace "COM7" with the COM port on which your arduino is connected port =new Serial(this, "COM7", 9600);} void draw(){// Read and parse incoming serial message serialEvent();// Set background to black background(); printQuaternions(); printYawPitchRoll();// Set position to centre translate(width /2, height /2);// Begin object pushMatrix();float[] axis = quat.toAxisAngle(); rotate(axis[], axis[2], axis[3], axis[1]);// Draw main body in red drawBody();// Draw front-facing tip in blue drawCylinder();// Draw Triangles drawTriangles();// Draw Quads drawQuads();// End of object popMatrix();// Send character 's' to Arduino//port.write('s');} void serialEvent(){// Read from port until new line (ASCII code 13) message = port.readStringUntil(newLine);if(message !=null){// Split message by commas and store in String array massQ = split(message, ","); q[]=float(massQ[]); q[1]=float(massQ[1]); q[2]=float(massQ[2]); q[3]=float(massQ[3]);}// Print values to console print(q[]); print("t"); print(q[1]); print("t"); print(q[2]); print("t"); print(q[3]); println("t");// Set our toxilibs quaternion to new data quat.set(q[], q[1], q[2], q[3]); } void drawCylinder(){float topRadius =;float bottomRadius =20;float tall =20;int sides =8;// Begin object pushMatrix(); translate(, , -120); rotateX(PI/2); fill(, , 255, 200); float angle =;float angleIncrement = TWO_PI / sides; beginShape(QUAD_STRIP);for(int i =; i < sides 1; i){ vertex(topRadius * cos(angle), , topRadius * sin(angle)); vertex(bottomRadius * cos(angle), tall, bottomRadius * sin(angle)); angle = angleIncrement;} endShape(); // if it is not a cone, draw the circular top capif(topRadius !=){ angle =; beginShape(TRIANGLE_FAN);// Center point vertex(, , );for(int i =; i < sides 1; i ){ vertex(topRadius * cos(angle), , topRadius * sin(angle)); angle = angleIncrement;} endShape();} // If it is not a cone, draw the circular bottom capif(bottomRadius !=){ angle =; beginShape(TRIANGLE_FAN);// Center point vertex(, tall, );for(int i =; i < sides 1; i ){ vertex(bottomRadius * cos(angle), tall, bottomRadius * sin(angle)); angle = angleIncrement;} endShape();} popMatrix();} void drawBody(){ fill(255, , , 200); box(10, 10, 200);} void drawTriangles(){// Draw wings and tail fin in green fill(, 255, , 200); beginShape(TRIANGLES);// Wing top layer vertex(-100, 2, 30); vertex(, 2, -80); vertex(100, 2, 30);// Wing bottom layer vertex(-100, -2, 30); vertex(, -2, -80); vertex(100, -2, 30);// Tail left layer vertex(-2, , 98); vertex(-2, -30, 98); vertex(-2, , 70);// Tail right layer vertex(2, , 98); vertex(2, -30, 98); vertex(2, , 70); endShape();} void drawQuads(){ beginShape(QUADS); vertex(-100, 2, 30); vertex(-100, -2, 30); vertex(, -2, -80); vertex(, 2, -80); vertex(100, 2, 30); vertex(100, -2, 30); vertex(, -2, -80); vertex(, 2, -80); vertex(-100, 2, 30); vertex(-100, -2, 30); vertex(100, -2, 30); vertex(100, 2, 30); vertex(-2, , 98); vertex(2, , 98); vertex(2, -30, 98); vertex(-2, -30, 98); vertex(-2, , 98); vertex(2, , 98); vertex(2, , 70); vertex(-2, , 70); vertex(-2, -30, 98); vertex(2, -30, 98); vertex(2, , 70); vertex(-2, , 70); endShape();} void printQuaternions(){// Set text mode to shape textMode(SHAPE); textSize(13); fill(255, 255, 255); text("Quaternions:", 20, 20, 10); text(q[], 20, 40, 10); text(q[1], 20, 60, 10); text(q[2], 20, 80, 10); text(q[3], 20, 100, 10);} void printYawPitchRoll(){// Calculate yaw/pitch/roll angles ypr[]= atan2(2*q[1]*q[2]-2*q[]*q[3], 2*q[]*q[] 2*q[1]*q[1]-1)*57.2; ypr[1]= atan2(2* q[2]* q[3]-2* q[]* q[1], 2* q[]* q[] 2* q[3]* q[3]-1)*57.2; ypr[2]=-atan2(2*(q[]* q[2]- q[1]* q[3]), 1-2*(q[2]* q[2] q[1]*q[1]))*57.2; text("Yaw:", 150, 20, 10); text(ypr[], 150, 40, 10); text("Pitch:", 220, 20, 10); text(ypr[1], 220, 40, 10); text("Roll:", 290, 20, 10); text(ypr[2], 290, 40, 10);}
При запуске визуализации на «processinge» откроется окно с графическим отображением IMU-сенсора в виде самолёта. Самолёт на экране будет повторять перемещения IMU-сенсора в пространстве.
Магнитометр
Датчик, который измеряет проекцию индукции магнитного поля на ось чувствительности. Магнитометру свойственны искажения hard-iron и soft-iron. Hard-iron искажение — это аддитивный эффект, когда к измеряемому полю добавляется постоянная составляющая. Причиной может быть, например, действие постоянного магнита или собственное смещение нуля датчика.
Искажение soft-iron — мультипликативный эффект, отражающий изменение направления и/или ослабление вектора магнитной индукции. Этот эффект может быть вызван наличием металлического предмета в непосредственной близости от магнитометра или же собственными искажениями датчика — погрешностью масштабного коэффициента или перекосом его оси чувствительности.
Примем модель измерений триады магнитометров:
где
– измерения магнитометра в ССК XYZ,
– диагональная матрица перекоса осей и масштабных коэффициентов (которая описывает эффект soft–iron),
– вектор истинной магнитной индукции в ССК,
– смещение нулей магнитометра (описывает эффект hard–iron),
– шум измерений.
Матрица перекоса осей и масштабных коэффициентов магнитометра:
элементы, расположенные на главной диагонали (
) — масштабные коэффициенты и их погрешности по трём осям магнитометра, остальные элементы матрицы — перекосы осей магнитометра. Все элементы матрицы также учитывают эффект soft-iron.
Общие сведения о датчике mpu9250
Датчик MPU9250 реализует функции 3-осевого гироскопа, 3-осевого акселерометра и 3-осевого магнитометра, то есть он представляет собой 9-ти осевой IMU сенсор (Inertial Measurement Unit – инерционное измерительное устройство). Его внешний вид с распиновкой контактов показан на следующем рисунке.
Назначение контактов датчика MPU9250:
Обозначение контакта | Назначение контакта |
VCC | Внешнее питание 3.3 В |
GND | Общий |
SCL | Линия тактовых импульсов I2C и SPI |
SDA | Линия данных для I2C или SPI |
EDA | Линия данных при подключении внешних датчиков по шине I2C |
ECL | Линия тактов при подключении внешних датчиков по шине I2C |
AD0 | Для выставления адреса I2C в режиме I2C. В режиме SPI это линия данных от датчика |
INT | Линия прерываний. Срабатывание настраивается при конфигурировании датчика MPU-9250 |
NCS | В режиме SPI – линия выбора ведомого (chip select). В режиме I2C не соединяется ни с чем |
FSYNC | Зависит от конфигурации |
MPU9250 является самым миниатюрным в мире девятиосевым датчиком. Это говорит о высокой производительности его микросхемы. Состоит корпус датчика из двух мельчайших кристаллов, один из которых отвечает за гироскоп и акселерометр, а другой за магнитометр.
Рассмотрим более подробно работу составных частей датчика.
Гироскоп – это сенсор, реагирующий на изменение углов ориентации в пространстве.
Акселерометр сравнивает проекцию ускорения объекта с гравитационным ускорением и способен замерять линейную скорость объекта, а вкупе с гироскопом – и положение в пространстве.
Магнитометр представляет собой устройство для измерения интенсивности ближайшего магнитного поля, действующего на объект.
Датчик MPU9250 позволяет измерять такие параметры как угол yaw (рыскание), угол pitch (тангаж) и угол roll (крен). Графическое изображение этих параметров представлено на следующем рисунке.
Ориентация объекта в пространстве
Определим ориентацию объекта в пространстве с учётом направления на север. Для этого показания акселерометра, гироскопа и компаса используем как входные данные для фильтра Маджвика. На выходе получим кватернионы, которые для наглядности пересчитаем в самолётные углы Эйлера.
- imu_9DOF.py
# Импортируем необходимые модулиfrom amperka_pico_imu_10dof import icm20948 importtimeimportmath # Печатаем тестовую строку в консольprint("IMU-sensor Test Program")# Создаём объект для работы с IMU-сенсором imu=icm20948.ICM20948() whileTrue: # Считываем данных с акселерометра accelerometer = imu.read_accelerometer_g()# Считываем данных с гироскопа gyroscope = imu.read_gyroscope_rad()# Считываем данных с магнитометра magnetometer = imu.read_magnetometer_uT() ax = accelerometer[] ay = accelerometer[1] az = accelerometer[2] gx = gyroscope[] gy = gyroscope[1] gz = gyroscope[2] mx = magnetometer[] my = magnetometer[1] mz = magnetometer[2]# Обновляем входные данные в фильтр Маджвика imu.AHRS_update(gx, gy, gz, ax, ay, az, mx, my, mz)# Получаем из фильтра углы Эйлера: yaw, pitch и roll yaw = imu.get_yaw(); pitch = imu.get_pitch(); roll = imu.get_roll();print(f"Yaw:t{yaw}")print(f"Pitch:t{pitch }")print(f"Roll:t{roll}")print()time.sleep(0.1)
Получение “реальных” значений
Наконец, на завершающем этапе обработки данных, мы получим значения параметров yaw, pitch и roll из кватернионов.
yaw =atan2(2.0f*(q[1]*q[2] q[]*q[3]),q[]*q[] q[1]*q[1]–q[2]*q[2]–q[3]*q[3]); pitch=–asin(2.0f*(q[1]*q[3]–q[]*q[2])); roll =atan2(2.0f*(q[]*q[1] q[2]*q[3]),q[]*q[]–q[1]*q[1]–q[2]*q[2] q[3]*q[3]); pitch*=180.0f/PI; yaw *=180.0f/PI; yaw =1.34;/* Declination at Potheri, Chennail ,India Model Used: IGRF12 Help Latitude: 12.823640° N Longitude: 80.043518° E Date Declination 2022-04-09 1.34° W changing by 0.06° E per year ( ve for west )*/ roll *=180.0f/PI; Serial.print(“Yaw,Pitch,Roll:“); Serial.print(yaw 180,2); Serial.print(“,“); Serial.print(pitch,2); Serial.print(“,“); Serial.println(roll,2); |
Преобразователь питания rt9193-18
Понижающий DC-DC преобразователь RT9193-18 обеспечивает питание инерциального датчика и светодиодной индикации. Входное напряжение поступает от преобразователя RT9193-33. Выходное напряжение 1,8 В с максимальным выходным током 300 мА.
Преобразователь питания rt9193-33
Работа с arduino и mpu6050 | alexgyver
Для достижения максимальной точности измерений нужно откалибровать акселерометр и гироскоп. Калибровка акселерометра позволяет выставить “ноль” для вектора силы тяжести, а калибровка гироскопа уменьшает его “дрифт”, то есть статическое отклонение в режиме покоя. Идеально откалиброванный и лежащий горизонтально датчик должен показывать ускорение ~16384 по оси Z и нули по всем остальным осям ускорения и угловой скорости. Но это фантастика =)
Максимально правильно использовать калибровку в проекте нужно так: калибровка по запросу (кнопка, меню, и т.д.), затем запись калибровочных значений в EEPROM. При запуске – чтение и настройка оффсетов. Да, можно замерить значения по всем 6 осям в покое, сохранить их в переменные, а затем вычитать из свежих прочитанных в процессе работы. Такой способ работает для каких-то базовых операций с датчиком (определение перегрузок, тряски, наличия вращения, и т.д.). Для использования MPU6050 с целью максимально точного определения углов поворота платы такой вариант к сожалению не подходит: калибровать нужно
рекурсивно
.
Рассмотрим несколько примеров калибровки, первый – из библиотеки. Калибрует долго, но максимально точно.
При малых вибрациях и движениях датчика в процессе калибровки (даже от громкого звука) калибровка может не закончиться
. Второй вариант – мой упрощённый алгоритм калибровки, калибрует быстро, без возможности зависнуть при тряске, но даёт менее точный результат. Я делюсь примерами, в свой проект их нужно будет переносить вручную и аккуратно =)
Также в библиотеке есть готовые функции для калибровки акселя и гиро, они принимают количество итераций калибровки (до 15):
mpu.CalibrateAccel(6); mpu.CalibrateGyro(6);
Функции блокирующие, выполняются и автоматически выставляют оффсеты, т.е. датчик сразу калибруется. Для чтения оффсетов (например для записи в EEPROM) можно воспользоваться тем же способом, что и раньше:
mpu.getXAccelOffset(); mpu.getYAccelOffset(); mpu.getZAccelOffset(); mpu.getXGyroOffset(); mpu.getYGyroOffset(); mpu.getZGyroOffset();
Расчет углов ориентации
Благодаря наличию на Земле силы тяжести, акселерометры «чувствуют» направление вниз. Их измерения используются для расчета углов крена и тангажа. Формулы для расчёта можно найти
. Третий — угол рыскания (а в данном случае — магнитного азимута), может быть определен благодаря наличию у Земли магнитного поля. Вектор индукции магнитного поля измеряется магнитометрами и их измерения участвуют в расчете угла рыскания. Нужно отметить, что в расчёте магнитного азимута используются измерения магнитометра, пересчитанные в плоскость.
можно найти формулу для расчёта магнитного азимута.
где
— функция полного арктангенса,
— измерения акселерометра по трём осям в ССК,
— измерения магнитометра по осям X’, Y’ (измерения магнитометров пересчитаны в плоскость).
Считывание данных с датчика mpu9250 в плату arduino mega 2560
Для работы с датчиком MPU9250 существует достаточно много библиотек, но мы будем использовать одну из самых популярных среди них: MPU-9250 Arduino Library by Kriswiner.
Когда вы сохраните библиотеку в ваш каталог Arduino, можно приступать к работе с ней. Откройте пример MPU9250BasicAHRS.ino.
Также выполните следующие соединения в схеме:
MPU9250 Breakout ——— Arduino• VIN ———————- 5V• SDA ———————– SDA (Pin 20)• SCL ———————– SCL (Pin 21)• GND ———————- GND
Помните о том, что соединительные провода в нашем случае не должны быть очень длинными потому что протокол I2C плохо работает на больших расстояниях.
Углы ориентации
Будем понимать под углами ориентации объекта углы Эйлера — крен (roll), тангаж (pitch), рыскание (yaw), связывающие собственную систему координат XYZ объекта и локальную систему координат восток-север-верх (ENU — East North Up). Углы roll, pitch, yaw обозначают поворот, который нужно совершить осям XYZ чтобы перейти в оси ENU. Соответственно, нулевые углы означают, что ось X объекта смотрит на восток, ось Y объекта смотрит на север, ось Z — вверх.
Порядок поворота осей — начиная с последнего угла: сначала на yaw (вокруг оси Z), потом на pitch (вокруг оси Y), потом на roll (вокруг оси X).
Фильтрация данных
Поскольку исходные данные (raw data) содержат достаточно большое количество шума, мы будем использовать различные фильтры (Madgwick/Mahony/Kalman) чтобы преобразовать их в кватернионы (Quaternions).
voidMadgwickQuaternionUpdate(floatax,floatay,floataz,floatgx,floatgy,floatgz,floatmx,floatmy,floatmz) { floatq1=q[],q2=q[1],q3=q[2],q4=q[3]; // short name local variable for readability floatnorm; floathx,hy,_2bx,_2bz; floats1,s2,s3,s4; floatqDot1,qDot2,qDot3,qDot4; // Auxiliary variables to avoid repeated arithmetic float_2q1mx; float_2q1my; float_2q1mz; float_2q2mx; float_4bx; float_4bz; float_2q1=2.0f*q1; float_2q2=2.0f*q2; float_2q3=2.0f*q3; float_2q4=2.0f*q4; float_2q1q3=2.0f*q1*q3; float_2q3q4=2.0f*q3*q4; floatq1q1=q1*q1; floatq1q2=q1*q2; floatq1q3=q1*q3; floatq1q4=q1*q4; floatq2q2=q2*q2; floatq2q3=q2*q3; floatq2q4=q2*q4; floatq3q3=q3*q3; floatq3q4=q3*q4; floatq4q4=q4*q4; // Normalise accelerometer measurement norm=sqrt(ax*ax ay*ay az*az); if(norm==0.0f)return;// handle NaN norm=1.0f/norm; ax*=norm; ay*=norm; az*=norm; // Normalise magnetometer measurement norm=sqrt(mx*mx my*my mz*mz); if(norm==0.0f)return;// handle NaN norm=1.0f/norm; mx*=norm; my*=norm; mz*=norm; // Reference direction of Earth’s magnetic field _2q1mx=2.0f*q1*mx; _2q1my=2.0f*q1*my; _2q1mz=2.0f*q1*mz; _2q2mx=2.0f*q2*mx; hx=mx*q1q1–_2q1my*q4 _2q1mz*q3 mx*q2q2 _2q2*my*q3 _2q2*mz*q4–mx*q3q3–mx*q4q4; hy=_2q1mx*q4 my*q1q1–_2q1mz*q2 _2q2mx*q3–my*q2q2 my*q3q3 _2q3*mz*q4–my*q4q4; _2bx=sqrt(hx*hx hy*hy); _2bz=–_2q1mx*q3 _2q1my*q2 mz*q1q1 _2q2mx*q4–mz*q2q2 _2q3*my*q4–mz*q3q3 mz*q4q4; _4bx=2.0f*_2bx; _4bz=2.0f*_2bz; // Gradient decent algorithm corrective step s1=–_2q3*(2.0f*q2q4–_2q1q3–ax) _2q2*(2.0f*q1q2 _2q3q4–ay)–_2bz*q3*(_2bx*(0.5f–q3q3–q4q4) _2bz*(q2q4–q1q3)–mx) (–_2bx*q4 _2bz*q2)*(_2bx*(q2q3–q1q4) _2bz*(q1q2 q3q4)–my) _2bx*q3*(_2bx*(q1q3 q2q4) _2bz*(0.5f–q2q2–q3q3)–mz); s2=_2q4*(2.0f*q2q4–_2q1q3–ax) _2q1*(2.0f*q1q2 _2q3q4–ay)–4.0f*q2*(1.0f–2.0f*q2q2–2.0f*q3q3–az) _2bz*q4*(_2bx*(0.5f–q3q3–q4q4) _2bz*(q2q4–q1q3)–mx) (_2bx*q3 _2bz*q1)*(_2bx*(q2q3–q1q4) _2bz*(q1q2 q3q4)–my) (_2bx*q4–_4bz*q2)*(_2bx*(q1q3 q2q4) _2bz*(0.5f–q2q2–q3q3)–mz); s3=–_2q1*(2.0f*q2q4–_2q1q3–ax) _2q4*(2.0f*q1q2 _2q3q4–ay)–4.0f*q3*(1.0f–2.0f*q2q2–2.0f*q3q3–az) (–_4bx*q3–_2bz*q1)*(_2bx*(0.5f–q3q3–q4q4) _2bz*(q2q4–q1q3)–mx) (_2bx*q2 _2bz*q4)*(_2bx*(q2q3–q1q4) _2bz*(q1q2 q3q4)–my) (_2bx*q1–_4bz*q3)*(_2bx*(q1q3 q2q4) _2bz*(0.5f–q2q2–q3q3)–mz); s4=_2q2*(2.0f*q2q4–_2q1q3–ax) _2q3*(2.0f*q1q2 _2q3q4–ay) (–_4bx*q4 _2bz*q2)*(_2bx*(0.5f–q3q3–q4q4) _2bz*(q2q4–q1q3)–mx) (–_2bx*q1 _2bz*q3)*(_2bx*(q2q3–q1q4) _2bz*(q1q2 q3q4)–my) _2bx*q2*(_2bx*(q1q3 q2q4) _2bz*(0.5f–q2q2–q3q3)–mz); norm=sqrt(s1*s1 s2*s2 s3*s3 s4*s4); // normalise step magnitude norm=1.0f/norm; s1*=norm; s2*=norm; s3*=norm; s4*=norm; // Compute rate of change of quaternion qDot1=0.5f*(–q2*gx–q3*gy–q4*gz)–beta*s1; qDot2=0.5f*(q1*gx q3*gz–q4*gy)–beta*s2; qDot3=0.5f*(q1*gy–q2*gz q4*gx)–beta*s3; qDot4=0.5f*(q1*gz q2*gy–q3*gx)–beta*s4; // Integrate to yield quaternion q1 =qDot1*deltat; q2 =qDot2*deltat; q3 =qDot3*deltat; q4 =qDot4*deltat; norm=sqrt(q1*q1 q2*q2 q3*q3 q4*q4); // normalise quaternion norm=1.0f/norm; q[]=q1*norm; q[1]=q2*norm; q[2]=q3*norm; q[3]=q4*norm; } // Similar to Madgwick scheme but uses proportional and integral filtering on the error between estimated reference vectors and // measured ones. voidMahonyQuaternionUpdate(floatax,floatay,floataz,floatgx,floatgy,floatgz,floatmx,floatmy,floatmz) { floatq1=q[],q2=q[1],q3=q[2],q4=q[3]; // short name local variable for readability floatnorm; floathx,hy,bx,bz; floatvx,vy,vz,wx,wy,wz; floatex,ey,ez; floatpa,pb,pc; // Auxiliary variables to avoid repeated arithmetic floatq1q1=q1*q1; floatq1q2=q1*q2; floatq1q3=q1*q3; floatq1q4=q1*q4; floatq2q2=q2*q2; floatq2q3=q2*q3; floatq2q4=q2*q4; floatq3q3=q3*q3; floatq3q4=q3*q4; floatq4q4=q4*q4; // Normalise accelerometer measurement norm=sqrt(ax*ax ay*ay az*az); if(norm==0.0f)return;// handle NaN norm=1.0f/norm; // use reciprocal for division ax*=norm; ay*=norm; az*=norm; // Normalise magnetometer measurement norm=sqrt(mx*mx my*my mz*mz); if(norm==0.0f)return;// handle NaN norm=1.0f/norm; // use reciprocal for division mx*=norm; my*=norm; mz*=norm; // Reference direction of Earth’s magnetic field hx=2.0f*mx*(0.5f–q3q3–q4q4) 2.0f*my*(q2q3–q1q4) 2.0f*mz*(q2q4 q1q3); hy=2.0f*mx*(q2q3 q1q4) 2.0f*my*(0.5f–q2q2–q4q4) 2.0f*mz*(q3q4–q1q2); bx=sqrt((hx*hx) (hy*hy)); bz=2.0f*mx*(q2q4–q1q3) 2.0f*my*(q3q4 q1q2) 2.0f*mz*(0.5f–q2q2–q3q3); // Estimated direction of gravity and magnetic field vx=2.0f*(q2q4–q1q3); vy=2.0f*(q1q2 q3q4); vz=q1q1–q2q2–q3q3 q4q4; wx=2.0f*bx*(0.5f–q3q3–q4q4) 2.0f*bz*(q2q4–q1q3); wy=2.0f*bx*(q2q3–q1q4) 2.0f*bz*(q1q2 q3q4); wz=2.0f*bx*(q1q3 q2q4) 2.0f*bz*(0.5f–q2q2–q3q3); // Error is cross product between estimated direction and measured direction of gravity ex=(ay*vz–az*vy) (my*wz–mz*wy); ey=(az*vx–ax*vz) (mz*wx–mx*wz); ez=(ax*vy–ay*vx) (mx*wy–my*wx); if(Ki>0.0f) { eInt[] =ex; // accumulate integral error eInt[1] =ey; eInt[2] =ez; } else { eInt[]=0.0f; // prevent integral wind up eInt[1]=0.0f; eInt[2]=0.0f; } // Apply feedback terms gx=gx Kp*ex Ki*eInt[]; gy=gy Kp*ey Ki*eInt[1]; gz=gz Kp*ez Ki*eInt[2]; // Integrate rate of change of quaternion pa=q2; pb=q3; pc=q4; q1=q1 (–q2*gx–q3*gy–q4*gz)*(0.5f*deltat); q2=pa (q1*gx pb*gz–pc*gy)*(0.5f*deltat); q3=pb (q1*gy–pa*gz pc*gx)*(0.5f*deltat); q4=pc (q1*gz pa*gy–pb*gx)*(0.5f*deltat); // Normalise quaternion norm=sqrt(q1*q1 q2*q2 q3*q3 q4*q4); norm=1.0f/norm; q[]=q1*norm; q[1]=q2*norm; q[2]=q3*norm; q[3]=q4*norm; } |