% Accelerometer state vector
  Xa=[ba S1 d12 S2 d13 d23 S3 (FX1-FI1)]';

% State transition matrix
  Fa=I

% Unknown accelerometer biases, normalized to the magnitude of gravity
  ba=[ba1;
    ba2;
    ba3];

% Unknown acceleration-squared nonlinearity for acceleration along the acelerometer input axis
  FI1=[FI11;
     FI12;
     FI13];

% Unknown acceleration-squared nonlinearity for acceleration normal to the acelerometer input axis
  FX1=[FX11;
     FX12;
     FX13];

% Unknown accelerometer scale factor errors
  S1 S2 S3

% Unknown accelerometer axes nonorthogonalities
  d12 d13 d23

% Midpoint components of acceleration in platform coordinates
  bm^2=[b1^2 0  0;
    0  b2^2 0;
    0  0  b3^2];

  bm=[b1 b2 b3]';

Ha=[b1,b3,b1^2,b1,b2,b1,b3,b2^2,b2,b3,b3^2,(1-b1^2)*b1,(1-b2^2)*b2,(1-b3^2)*b3];

%%%%%%%%%%%%%%
% EKF module %
%%%%%%%%%%%%%%
% Input initial parameters Rk, xa, Za, Ha, FY
% xa=[0 0 0 0 0 0 0 0 0 0 0 0]
Id=[1 0 0 0 0 0 0 0 0 0 0 0;
    0 1 0 0 0 0 0 0 0 0 0 0;
    0 0 1 0 0 0 0 0 0 0 0 0;
    0 0 0 1 0 0 0 0 0 0 0 0;
    0 0 0 0 1 0 0 0 0 0 0 0;
    0 0 0 0 0 1 0 0 0 0 0 0;
    0 0 0 0 0 0 1 0 0 0 0 0;
    0 0 0 0 0 0 0 1 0 0 0 0;
    0 0 0 0 0 0 0 0 1 0 0 0;
    0 0 0 0 0 0 0 0 0 1 0 0;
    0 0 0 0 0 0 0 0 0 0 1 0;
    0 0 0 0 0 0 0 0 0 0 0 1];
% Computing Kalman gain
  Kg = Pk*Ha'*inv(Ha*Pk*Ha'+Rk); %Kg = (Pk*Hk')/(Hk*Pk*Hk'+Rk);

% Computing the predicted measurment
  za = Ha*xa;

% Conditioning the predicted estimate on the measurement
  Xa = xa + Kg*(Za - za);

% Computing a posteriori covariance matrix for update estimate
  P2k= (Id-Kg*Ha)*Pk;

% Computing the a priori covariance matrix ahead for next time step
  Pk = FY*P2k*FY' + Qk;

% Computing the predicted state estimate
  xaf = (FY*Xa)';

% Accelerometer state vector
  Xa=[ba' S1 d12 S2 d13 d23 S3 (FX1-FI1)']';

Ha=[b1,b2,b3,b1^2,b1*b2,b1*b3,b2^2,b2*b3,b3^2,(1-b1^2)*b1,(1-b2^2)*b2,(1-b3^2)*b3]';

%Difference equation for the accelerometers
Xa(i)=Fa*Xa(i-1)+wa(i-1)

%Accelerometer observation equation
Za(i)=Ha*Xa(i)+Va(i)

The midpoint acceleration error in terms of the unknown accelerometer parameters is given by:
dbm=bm*psi+ba+ha*bm+(bm^2)*(FI1-FX1)+db;
ba - (3x1) vector of unknown accelerometer biases, normalized to the magnitude of gravity
ha - (3x3) matrix
    |S1 d12 d13|
ha= |0  S2  d23|
    |0  0   S3 |
S1,S2,S3 - unknown accelerometer scale factor errors.
d12,d13,d23 - unknown accelerometer axes nonorthogonalities (misalignments)
db - represents other error terms, some of which are observable; for reason of practicality in our example they are not estimated, only compensated with factory calibrated values.
FI1 - (3x1) unknown acceleration-squared nonlinearity for acceleration along the accelerometer input axis
FX1 - (3x1) unknown acceleration-squared nonlinearity for acceleration normal to accelerometer input axis
bm - is a three vector (b1, b2, b3)' of midpoint components of acceleration in platform coordinates
      |b1^2 0    0   |
bm^2= |0    b2^2 0   |
      |0    0    b3^2|