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

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

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

switch flag
    case 0 % 初始化
        x0 = zeros(6, 1);
        sys = x0;
    case 2 % 系统状态更新
        % 获取输入向量
        gyro = u(1:3);
        accel = u(4:6);
        mag = u(7:9);
        
        % 定义旋转矩阵
        R_roll = [1 0 0; 0 cos(x(1)) -sin(x(1)); 0 sin(x(1)) cos(x(1))];
        R_pitch = [cos(x(2)) 0 sin(x(2)); 0 1 0; -sin(x(2)) 0 cos(x(2))];
        R_yaw = [cos(x(3)) -sin(x(3)) 0; sin(x(3)) cos(x(3)) 0; 0 0 1];
        R_gb = R_yaw * R_pitch * R_roll;
        
        % 计算测量向量
        z = [accel; mag];
        
        % 定义状态转移矩阵
        F = [eye(3) -R_gb * dt; zeros(3) eye(3)];
        
        % 定义观测矩阵
        H = [zeros(3) eye(3)];
        
        % 定义过程噪声标准差
        Q = diag([0.01 0.01 0.01 0.01 0.01 0.01]);
        
        % 定义测量噪声标准差
        R = diag([0.1 0.1 0.1]);
        
        % 计算残差
        y = z - H * x;
        
        % 计算卡尔曼增益
        K = x(4:end) * H' * (H * x(4:end) * H' + R)^-1;
        
        % 更新状态向量
        sys = x + K * y;
        
        % 更新协方差矩阵
        P = F * [x(4:end); x(1:3)] * F' + Q;
        sys = [sys(1:3); P(4:end)];
    case 3 % 系统输出
        % 输出姿态角
        sys = x(1:3);
    case {1, 4, 9} % 不进行任何操作
        sys = [];
    otherwise
        error(['Unhandled flag = ', num2str(flag)]);
end
end

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