数据源为.mat文件,点击这里免费下载
% 基于等效旋转矢量姿态更新算法 解算三轴陀螺仪的固连载体位姿,画出得到的三个姿态角
delta0 = [0, 0, 0]; %上一周期角增量
delta = [0, 0, 0]; %当前周期角增量
Q = [1, 0, 0, 0]; % 姿态四元数初始值
q = [0, 0, 0, 0]; % 姿态变化四元数初始值
eular = [0, 0, 0]; %用来保存Q对应的姿态角
load('gyro_data.mat');
% 求出欧拉角随时间的变化(固定轴欧拉角)
for i = 1:6950 % 数组索引必须为正整数或逻辑值,所以i从1开始
t(i) = i/1000; % 采样频率1kHz, 画图用
% 当前时刻角增量delta
delta_roll = gyro_data(i,1) * pi /180; % x翻滚角 弧度
delta_yaw = gyro_data(i,2) * pi /180; % y偏航角 弧度
delta_pitch = gyro_data(i,3) * pi /180; % z俯仰角 弧度
delta = [delta_roll, delta_yaw, delta_pitch];
% 旋转矢量fai (单子样+前一周期字样)
Fai = delta + cross(delta0, delta) / 12; % 3维旋转向量fai
% 旋转矢量fai的模,也就是旋转角度theta大小
theta = sqrt( sum(Fai.*Fai) ); %也可以直接用二范数 norm(Fai,2)
% 由旋转矢量 计算姿态变化四元数
% q_w = cos(theta/2); //对计算公式理解不对
% q_i = cos(theta/2) + Fai(1) / theta * sin(theta/2);
% q_j = cos(theta/2) + Fai(2) / theta * sin(theta/2);
% q_k = cos(theta/2) + Fai(3) / theta * sin(theta/2);
q_w = cos(theta/2);
q_i = Fai(1) / theta * sin(theta/2);
q_j = Fai(2) / theta * sin(theta/2);
q_k = Fai(3) / theta * sin(theta/2);
q = [q_w, q_i, q_j, q_k]; %构建旋转四元数
% 四元数乘法实现方法1
Q = quatmultiply(Q, q);
% 四元数乘法实现方法2 -自定义函数
% Q = q1oq2(Q, q);
% 四元数乘法实现方法3 -反对称矩阵+矩阵乘法
% AntisA = quat2Antisymmetry(Q(1), Q(2), Q(3), Q(4)) % Q的反对称矩阵
% Q = transpose (AntisA * transpose(q));
%四元数转为欧拉角,绘制姿态角
eular = quat2eul(Q);
% 记录载体历史姿态
eul(i,1) = eular(1) *180 /pi; % x-roll 累积角度随时间变化曲线
eul(i,2) = eular(2) *180 /pi; % y-yaw 累积角度随时间变化曲线
eul(i,3) = eular(3) *180 /pi; % z-pitch累积角度随时间变化曲线
angle(i,1) = Fai(1) *180 /pi; % x-roll增量随时间变化曲线
angle(i,2) = Fai(2) *180 /pi; % y-yaw 增量随时间变化曲线
angle(i,3) = Fai(3) *180 /pi; % z-pitch增量随时间变化曲线
% 上一周期变量
delta_roll0 = delta_roll;
delta_yaw0 = delta_yaw;
delta_pitch0= delta_pitch;
end
t = t - 0.001; % i从1开始,初始时刻是从0,要减1
figure(1);
subplot(6,1,1); plot(t, eul(:,1), 'r'); title('x-roll /°'); grid on;
subplot(6,1,2); plot(t, eul(:,2), 'g'); title('y-yaw /°'); grid on;
subplot(6,1,3); plot(t, eul(:,3), 'b'); title('z-pitch /°'); grid on;
subplot(6,1,4); plot(t, angle(:,1), 'r'); title('x-roll 增量随时间变化曲线 /°'); grid on;
subplot(6,1,5); plot(t, angle(:,2), 'g'); title('y-yaw 增量随时间变化曲线 /°'); grid on;
subplot(6,1,6); plot(t, angle(:,3), 'b'); title('z-pitch 增量随时间变化曲线 /°'); grid on;
% 四元数乘法实现 方法2-直接返回结果
function Q_ = q1oq2(q1, q2)
a1 = q1(1)*q2(1) - q1(2)*q2(2) - q1(3)*q2(3) - q1(4)*q2(4);
b1 = q1(1)*q2(2) + q1(2)*q2(1) + q1(3)*q2(4) - q1(4)*q2(3);
c1 = q1(1)*q2(3) + q1(3)*q2(1) + q1(4)*q2(2) - q1(2)*q2(4);
d1 = q1(1)*q2(4) + q1(4)*q2(1) + q1(2)*q2(3) - q1(3)*q2(2);
Q_ = [a1, b1, c1, d1];
end
% 四元数乘法实现 方法3-求向量q的反对称矩阵,再借助矩阵乘法计算,调用见主函数
function Aq_ = quat2A(qw, qi, qj, qk)
Aq_ = [ qw -qi -qj -qk;
qi qw -qk qj;
qj qk qw -qi;
qk -qj qi qw ];
end