Measured IMU position inside the iPhone 17

Code:
% Estimation of the IMU position (iPhone 17, phyphox) from several turntable runs
%
% Rigid-body model for a rotation about the device z-axis (phone frame, phyphox
% convention: specific force, +x right, +y up, +z out of the display). With
% l = (lx, ly) being the vector from the IMU to the rotation axis:
%   ax = wz²*lx + dwz/dt*ly + bx
%   ay = wz²*ly - dwz/dt*lx + by
% (bx, by: sensor offset + gravity leakage due to tilt)
%
% Combination of the runs: the phone sits in the same pocket in every run, only
% its orientation changes. With c = vector housing center -> rotation axis in the
% phone frame of the reference run (display up, see turntable.png) and p = IMU
% position relative to the housing center:
%   l_k = S_k*c - p
% S_k maps the reference placement into the phone frame of run k:
%   display up   (reference)            S = [ 1 0; 0  1]
%   display down, flipped about y-axis  S = [-1 0; 0  1]
%   display down, flipped about x-axis  S = [ 1 0; 0 -1]
%   display up,   rotated by 180°       S = [-1 0; 0 -1]
% Components of c that the runs cannot resolve are taken from the setup prior
% (R, phi). With the runs 'up' + 'down about y' this applies to the y-direction.

% clean up
clear all;
close all;
clc;

% settings
runs = struct( ...
  'file', {'Translation und Rotation 2026-09-29 16-10-10.xlsx', ...
           'Translation und Rotation 2026-10-26 15-35-40.xlsx', ...
           'Translation und Rotation 2026-10-27 13-36-20.xlsx'}, ...
  'S',    {[1 0; 0 1], [-1 0; 0 1], [1 0; 0 -1]}, ...
  'up',   {true, false, false});  % display up? (checked against az)
R      = 120;    % prior: distance housing center -> rotation axis (mm)
sigR   = 1;      % uncertainty of R (mm, 1 sigma)
phi    = 0;      % prior: angle of the rotation axis seen from the housing center (deg)
sigphi = 2;      % uncertainty of phi (deg, 1 sigma)
sigpos = 1;      % repositioning tolerance of the phone in the pocket (mm, 1 sigma)
Tsm    = 0.10;   % smoothing window (s)
taus   = -0.05 : 0.001 : 0.05;  % search grid for the gyro time shift (s)
Lblk   = 1.0;    % block length for the moving-block bootstrap (s)
B      = 2000;   % number of bootstrap samples
conf   = 0.95;   % confidence level
m2mm   = 1000;
rand('seed', 42);
randn('seed', 42);

% load a phyphox export (cached as .mat next to the .xlsx)
function [ACC, GYRO] = load_run(file)
  [~, name] = fileparts(file);
  cache = [name, '.mat'];
  if exist(cache, 'file')
    load(cache);
  else
    pkg load io;
    pkg load windows;
    ACC  = xlsread(file, 'Accelerometer');  % t (s), ax, ay, az (m/s²)
    GYRO = xlsread(file, 'Gyroscope');      % t (s), wx, wy, wz (rad/s)
    save('-v7', cache, 'ACC', 'GYRO');
  end
end

% build the (smoothed) linear model for a given gyro time shift tau
function [Y, X] = build_model(ta, ax, ay, tw, wz, tau, nsm, keep)
  w  = interp1(tw + tau, wz, ta, 'linear', 'extrap');
  dw = gradient(w, ta);
  % the same linear filter is applied to data and regressors -> model stays exact
  f  = @(v) movmean(v, nsm);
  W2 = f(w.^2)(keep);
  A  = f(dw)(keep);
  ax = f(ax)(keep);
  ay = f(ay)(keep);
  o  = ones(size(W2));
  z  = zeros(size(W2));
  Y  = [ax; ay];
  X  = [W2,  A, o, z;
        -A, W2, z, o];
end

% fit a single run: time shift, least squares, moving-block bootstrap
function r = fit_run(ACC, GYRO, Tsm, taus, Lblk, B)
  ta  = ACC(:, 1) - ACC(1, 1);
  tw  = GYRO(:, 1) - ACC(1, 1);
  dt  = median(diff(ta));
  nsm = max(1, round(Tsm/dt));
  % only use samples covered by the gyro, minus the filter edges
  keep = ta > tw(1) + max(abs(taus)) & ta < tw(end) - max(abs(taus));
  keep(1 : nsm) = false;
  keep(end - nsm + 1 : end) = false;
  n = sum(keep);

  ssr = zeros(size(taus));
  for k = 1 : numel(taus)
    [Y, X] = build_model(ta, ACC(:, 2), ACC(:, 3), tw, GYRO(:, 4), taus(k), nsm, keep);
    ssr(k) = sum((Y - X*(X\Y)).^2);
  end
  [~, kbest] = min(ssr);
  r.tau = taus(kbest);
  r.ssr = ssr;

  [Y, X] = build_model(ta, ACC(:, 2), ACC(:, 3), tw, GYRO(:, 4), r.tau, nsm, keep);
  th = X\Y;
  r.th  = th;
  r.rms = sqrt(mean((Y - X*th).^2));
  r.lx_rad = X(1:n, 1)\(Y(1:n) - th(3) - X(1:n, 2)*th(2));              % ax ~ wz²*lx
  r.lx_tan = X(n+1:end, 1)\(Y(n+1:end) - th(4) - X(n+1:end, 2)*th(2));  % ay ~ -dwz/dt*lx
  r.t = ta(keep);
  r.Y = Y;
  r.Ym = X*th;
  r.n = n;
  r.az = mean(ACC(:, 4));

  % residuals are strongly autocorrelated -> resample blocks
  L   = max(1, round(Lblk/dt));
  nb  = ceil(n/L);
  thb = zeros(B, numel(th));
  for b = 1 : B
    s   = randi(n - L + 1, nb, 1);
    idx = (s + (0 : L - 1))';
    idx = idx(1 : n)';
    thb(b, :) = ([X(idx, :); X(n + idx, :)]\[Y(idx); Y(n + idx)])';
  end
  r.thb = thb;
end

% combine runs: weighted least squares for u = [c; p] incl. setup prior on c
function u = combine(l, S, c0, Wc, sigpos)
  M = [];
  y = [];
  for k = 1 : numel(l)
    M = [M; [S{k}, -eye(2)]/sigpos];
    y = [y; l{k}/sigpos];
  end
  M = [M; [Wc, zeros(2)]];
  y = [y; Wc*c0];
  u = M\y;
end

% 1) fit every run
nr = numel(runs);
for k = 1 : nr
  [ACC, GYRO] = load_run(runs(k).file);
  res(k) = fit_run(ACC, GYRO, Tsm, taus, Lblk, B);
  if (res(k).az > 0) != runs(k).up
    warning('run %d: az = %.2f m/s² does not match up = %d', k, res(k).az, runs(k).up);
  end
end

% 2) combine the runs (nominal + Monte Carlo over bootstrap, prior and repositioning)
c0  = R*[cosd(phi); sind(phi)];
Ec  = [cosd(phi), -sind(phi); sind(phi), cosd(phi)];   % radial/tangential directions
Sc  = Ec*diag([sigR, R*sigphi*pi/180].^2)*Ec';          % prior covariance of c
Wc  = chol(inv(Sc));
S   = {runs.S};
l   = arrayfun(@(r) r.th(1:2)*m2mm, res, 'UniformOutput', false);
u   = combine(l, S, c0, Wc, sigpos);
ub  = zeros(B, 4);
for b = 1 : B
  lb = arrayfun(@(r) r.thb(b, 1:2)'*m2mm + sigpos*randn(2, 1), res, 'UniformOutput', false);
  ub(b, :) = combine(lb, S, c0 + chol(Sc)'*randn(2, 1), Wc, sigpos)';
end
c  = u(1:2);
p  = u(3:4);
pb = ub(:, 3:4);
% how much does each coordinate of p depend on the prior? (1 = only prior, 0 = only data)
J  = zeros(2);
for j = 1 : 2
  dc = zeros(2, 1);
  dc(j) = 1;
  du = combine(l, S, c0 + dc, Wc, sigpos) - u;
  J(:, j) = du(3:4);
end

% report
q  = [(1 - conf)/2, (1 + conf)/2];
ci = @(v) quantile(v, q);
printf('Single runs (vector IMU -> rotation axis in device coordinates, %2.0f%% CI, block bootstrap):\n', conf*100);
for k = 1 : nr
  r = res(k);
  printf('\n  Run %d: %s\n', k, runs(k).file);
  printf('    display %s (az = %+.2f m/s²), tau = %+.0f ms, RMS = %.3f m/s²\n', ...
         {'down', 'up  '}{runs(k).up + 1}, r.az, r.tau*1000, r.rms);
  printf('    lx = %8.2f mm  [%8.2f, %8.2f]   (radial %.2f / tangential %.2f)\n', ...
         r.th(1)*m2mm, ci(r.thb(:, 1)*m2mm), r.lx_rad*m2mm, r.lx_tan*m2mm);
  printf('    ly = %8.2f mm  [%8.2f, %8.2f]\n', r.th(2)*m2mm, ci(r.thb(:, 2)*m2mm));
  printf('    bx = %6.3f m/s², by = %6.3f m/s²\n', r.th(3), r.th(4));
end
printf('\nConsistency: residuals of the combined equations l_k = S_k*c - p (mm):\n');
for k = 1 : nr
  printf('  Run %d: [%+6.2f, %+6.2f]\n', k, S{k}*c - p - l{k});
end
printf('\nRotation axis seen from the housing center (reference placement): c = [%.2f, %.2f] mm\n', c);
printf('  => R = %.2f mm (prior %.1f +- %.1f), phi = %.2f deg (prior %.1f +- %.1f)\n', ...
       norm(c), R, sigR, atan2d(c(2), c(1)), phi, sigphi);
printf('\nIMU position relative to the housing center (x right, y up, looking at the display):\n');
printf('  x = %7.2f mm   [%7.2f, %7.2f]   share of prior R/phi: %3.0f %%\n', p(1), ci(pb(:, 1)), 100*norm(J(1, :)));
printf('  y = %7.2f mm   [%7.2f, %7.2f]   share of prior R/phi: %3.0f %%\n', p(2), ci(pb(:, 2)), 100*norm(J(2, :)));

% visualization
figure(1);
for k = 1 : nr
  r = res(k);
  subplot(nr, 2, 2*k - 1);
  plot(r.t, r.Y(1:r.n), '-', r.t, r.Ym(1:r.n), '--');
  grid on;
  ylabel(sprintf('run %d: a_x (m/s²)', k));
  subplot(nr, 2, 2*k);
  plot(r.t, r.Y(r.n+1:end), '-', r.t, r.Ym(r.n+1:end), '--');
  grid on;
  ylabel(sprintf('run %d: a_y (m/s²)', k));
end
xlabel('t (s)');
legend('measurement (smoothed)', 'model');

figure(2, 'position', [100 100 800 700]);
hold on;
rectangle('Position', [-71.5/2, -149.6/2, 71.5, 149.6], 'Curvature', 0.15);
plot(pb(:, 1), pb(:, 2), '.', 'Color', [0.7 0.7 0.9]);
plot(p(1), p(2), 'r+', 'MarkerSize', 14, 'LineWidth', 2);
plot(0, 0, 'k+');
axis equal;
axis([-60 60 -80 80]);
grid on;
xlabel('x (mm)');
ylabel('y (mm)');
title('iPhone 17: IMU position (looking at the display)');
legend('housing 71.5 x 149.6 mm', 'bootstrap / Monte Carlo', 'estimate', 'housing center', ...
       'location', 'eastoutside');
text(p(1) + 4, p(2) + 22, sprintf('x = %.1f mm [%.1f, %.1f]\ny = %.1f mm [%.1f, %.1f]\n(%2.0f%% CI)', ...
     p(1), ci(pb(:, 1)), p(2), ci(pb(:, 2)), conf*100));
print(2, 'imu_position_estimate.png', '-dpng', '-r120');