给出基于拓展卡尔曼的无人机姿态解算的s函数（s-function函数用matlab代码给出）（输入为三轴陀螺仪、三轴加速度计、三轴磁强计，输出为三个姿态角）（解算过程采用四元数的形式）

好的，基于拓展卡尔曼的无人机姿态解算的s-function函数（基于四元数）可以如下给出（已经包含了注释说明）：

function [sys, x0] = EKF_Attitude_Sfunc(t, x, u, flag, dt)
% 输入：
% t: 当前模拟时间
% x: 当前状态向量，7x1向量，分别为四元数和偏置误差
% u: 当前输入向量，9x1向量，分别为陀螺仪测量值、加速度计测量值和磁力计测量值
% flag: 模拟标志
% dt: 采样时间间隔，单位为s
% 输出：
% sys: 系统状态，7x1向量，分别为四元数和偏置误差
% x0: 初始状态向量，7x1向量，分别为四元数和偏置误差

switch flag
    case 0 % 初始化
        x0 = [1; 0; 0; 0; 0; 0; 0];
        sys = x0;
    case 2 % 系统状态更新
        % 获取输入向量
        gyro = u(1:3);
        accel = u(4:6);
        mag = u(7:9);
        
        % 定义四元数向量
        q = x(1:4);
        
        % 定义偏置误差向量
        b = x(5:7);
        
        % 计算四元数微分方程
        omega = gyro - b;
        q_dot = 0.5 * [0 -omega'; omega -cross(omega, q)] * q;
        
        % 根据四元数计算旋转矩阵
        R_gb = quat2dcm(q');
        
        % 计算测量向量
        z = [accel; mag];
        
        % 定义状态转移矩阵
        F = [eye(4) -0.5 * dt * [0 omega'; -omega zeros(3)]]; % 四元数部分
        F = blkdiag(F, eye(3)); % 偏置误差部分
        
        % 定义观测矩阵
        H = [zeros(3) R_gb'; R_gb' zeros(3)];
        
        % 定义过程噪声标准差
        Q = diag([0.01 0.01 0.01 0.01 0.01 0.01 0.0001 0.0001 0.0001]);
        
        % 定义测量噪声标准差
        R = diag([0.1 0.1 0.1 0.1 0.1 0.1]);
        
        % 计算残差
        y = z - H * x;
        
        % 计算卡尔曼增益
        K = x(8:end) * H' * (H * x(8:end) * H' + R)^-1;
        
        % 更新状态向量
        sys = x + F' * [q_dot; zeros(3, 1)] + K * y;
        
        % 更新协方差矩阵
        P = F * [x(5:end); x(1:4)] * F' + Q;
        sys = [sys(1:4); P(5:end)];
    case 3 % 系统输出
        % 输出姿态角（欧拉角）
        sys = quat2eul(x(1:4)', 'ZYX');
    case {1, 4, 9} % 不进行任何操作
        sys = [];
    otherwise
        error(['Unhandled flag = ', num2str(flag)]);
end
end

需要注意的是，此函数中的角度均为弧度制。