Indoor Walking Trajectory Estimation Using Mobile Device Sensors for Hand-Held and Hand-Swinging Modes ^†

Abstract

We propose an indoor location estimation method using sensors of mobile devices. First, we perform attitude estimation using each sensor. This estimation is used to estimate the attitude of the mobile device with respect to the earth. Based on the acceleration and other information obtained from the attitude estimation, we then estimate the step detection, step length, and direction of the step. Finally, the location is calculated using all the estimation results. To eliminate the need to hold the mobile device in place during the estimation process, the method is configured so that estimates may be performed while walking, while looking at the screen, and while walking and holding the device in one hand. As the proposed method does not use indoor location fingerprinting or machine learning, real-time estimation can be performed. Although the accuracy could be higher, our experimental results show that the proposed method is able to effectively estimate the location and walking trajectory.

1. Introduction

With the spread of mobile devices such as smartphones and tablets in recent years, applications that provide users with location information obtained by using the devices’ sensors have become popular, and the demand for such applications has been increasing. In general applications, Global Positioning System (GPS) and Light Detection And Ranging (LiDAR) are mainly used for position estimation. However, because GPS uses the transmission time and reception of radio waves, accurate location estimation requires an outdoor location. On the other hand, LiDAR requires nonmoving objects and the presence of obstructions or objects in the vicinity. GPS is unreliable in underground and indoor facilities due to the presence of obstacles that block or attenuate radio waves. It can also be difficult to estimate position using LiDAR due to the presence of humans around the facility. Although there are methods that use machine learning, it is difficult to solve the problem of non-periodic geomagnetic interference from other devices, and the need for pretraining makes high introduction costs a challenge for general users.

In indoor location estimation, a method called Pedestrian Dead Reckoning (PDR) is particularly applicable to pedestrians. PDR is a method of calculating relative position from a certain point based on information such as acceleration, angular velocity, geomagnetism, and atmospheric pressure. There are various approaches to position estimation. Essentially, there are three steps: step number estimation, step length estimation, and direction estimation. In step count estimation, methods include record acceleration and determining the number of steps by the magnitude of the acceleration value [1,2,3]. In stride estimation, different methods use fixed values [4,5,6,7] or acceleration [8,9,10]. In direction estimation, methods use the acceleration, angular velocity, and geomagnetism [3]. A number of the above methods have actually already been implemented in factories and commercial facilities; however, these often require a dedicated terminal to be attached to the pedestrian. Previous studies have obtained measurements with an Inertial Measurement Unit (IMU) attached to the foot [11] or a PDR with a dedicated sensor unit attached to the belt [12,13]. Most of these methods require fixation to the body or the use of a dedicated terminal. As mentioned at the beginning of this paper, applications that provide location information to users have become popular because they are easy to operate and use, with no cost and no complicated operations. Therefore, these currently introduced applications are unlikely to be accepted by general users. In order for the system to be usable by anyone, it is necessary to eliminate restrictions such as the need for it to be fixed to the body or reliance on a dedicated terminal. Although PDR using mobile terminals has been studied in the past, it is difficult to perform unconstrained self-position estimation because it is necessary to perform estimation in hand-held walking and hand-swinging walking conditions. In addition, although geomagnetic sensors are suitable for estimating absolute direction, they are difficult to use in practical situations because other electronic devices and structures affect the geomagnetic field and cause errors. Device-based positioning estimation uses the geomagnetic values [10,11,14]. However, previous studies have often failed to consider the influence of structures and other electronic devices.

In a recent trend, indoor positioning fingerprints such as RSSI using BLE or WiFi are sometimes used to correct estimated positions [15,16,17,18]. However, when using BLE it is necessary to install beacons in the indoor environment where estimation is performed, and even WiFi-based methods require prior setup. Additionally, estimation methods that use markers or indoor images [13,19] require directing the camera to a specified position, and also require pretraining when utilizing machine learning. In all cases, these methods result in significant setup costs for general users and lead to loss of real-time performance capability due to increased computation time.

In this study, we propose an indoor location estimation approach using only mobile devices and perform trajectory estimation based on multiple data types obtained from the posture and sensors of the mobile device while the user is walking. The proposed approach has low deployment cost for general users. It does not use indoor location fingerprinting or machine learning, allowing for real-time estimation. In addition, it implements an algorithm that is insensitive to geomagnetic effects from other devices, meaning that it can be used in any environment. Based on the positional relationship between a pedestrian and a mobile terminal, we propose a walking trajectory and location estimation system that divides the walking method into walking while looking at a screen (hand-held walking) and walking while holding a mobile device by hand (hand-swinging walking). In particular, we implement the proposed algorithm for estimating the walking trajectory during walking, then evaluate the trajectory and errors by conducting actual experiments.

2. System

2.1. System Configuration

In this study, trajectory estimation is performed using a group of sensors: a triaxial acceleration sensor, triaxial angular rate sensor, and triaxial geomagnetic sensor. These sensors are commonly installed in mobile devices. This work is an extended and refined version of research initially presented at the IEEE International Conference on Consumer Electronics—Taiwan [20]. The findings presented here include additional analysis and expanded discussion that goes beyond the scope of the original conference presentation.

In this study, PDR is divided into five steps: attitude estimation, step detection, step length estimation, step direction estimation, and position estimation. The estimation method is varied based on whether the subject is walking with a mobile device held in front of the user’s body (hand-held walking) or walking while swinging their hands to the side (hand-swinging walking). Figure 1 shows an overview of the proposed system configuration. The alphabets in Figure 1 indicate the order of processing.

The process flow in Figure 1 is shown below.

A

Acquire acceleration, angular velocity, and geomagnetic data from sensors

B

Attitude estimation using acceleration, angular velocity, and geomagnetic data

C

Perform step detection based on acceleration, angular velocity, and attitude estimation results

D

Estimate step length based on acceleration, angular velocity, and step detection results

E

Estimate step direction based on step detection and attitude estimation results

F

Perform position estimation based on the results of step detection, step length estimation, and step direction estimation

Development Environment

The device used in this study was a Google Pixel 4 with Android version 12. Each estimation was performed by an application developed using Unity. The sampling frequency of each sensor was 60 Hz. While some terminals can be set to a sampling frequency higher than 60 Hz, due to the specification that Unity processes at each drawing timing, we set the sampling frequency to 60 Hz, which is a common frame rate. Noise and offsets are problematic during the actual estimation process using sensors; thus, it is necessary to adjust the obtained data before performing each estimation.

2.2. Adjustment of Data Obtained from Sensors

2.2.1. Acceleration Sensor

The values obtained from the acceleration sensor require offset and coefficient adjustments.

Offset Adjustment

Let $g=9.798$ be the acceleration of gravity in Japan. When the sensor is placed on a plane, the value of acceleration should be $(x,y,z)=(0,0,g)$. In reality, the values obtained from the sensor are not $(x,y,z)=(0,0,g)$ due to the presence of errors in either the positive or negative direction. Therefore, it is necessary to set the value of the error and the opposite sign as the offset value. For example, if the measured value is $(x,y,z)=(0.1,-0.2,9.95)$ when placed on a horizontal plane, then the offset value $(x_{offset},y_{offset},z_{offset})$ is as follows:

$$\begin{array}{cc}\hfill (x_{offset},y_{offset},z_{offset})& =(0,0,g)-(0.1,-0.2,9.95)\hfill \\ \hfill & =(-0.1,0.2,-0.152).\hfill \end{array}$$

Coefficient Adjustment

The value after the offset adjustment should be g for the axis pointing in the direction of gravitational acceleration; in fact, however, it may be larger or smaller than g. Therefore, it is necessary to obtain the coefficients by making measurements with each axis of the sensor pointing vertically upward and downward. For example, if the x-axis is placed vertically upward and downward and the measured values are $x=-9.514,x=-10.501$, then the coefficient values in the positive and negative directions of the x-axis $x_{P_{coeff}},x_{N_{coeff}}$ are as follows:

$$\begin{array}{cc}\hfill x_{P_{coeff}}& =\left|{\displaystyle \frac{g}{-9.514}}\right|\hfill \\ \hfill & \fallingdotseq 1.030\hfill \\ \hfill x_{N_{coeff}}& =\left|{\displaystyle \frac{g}{10.501}}\right|\hfill \\ \hfill & \fallingdotseq 0.933.\hfill \end{array}$$

2.2.2. Angular Velocity Sensor

The values obtained from the angular velocity sensor require offset adjustment.

Offset Adjustment

When the sensor is placed and kept stationary, all the values of the angular rate sensors should be 0. In practice, however, a deviation exists in either the positive or negative direction. Therefore, the offset must be set to a value with the opposite sign to the deviated value. The method for determining the offset is the same as in the case of acceleration (Section 2.2.1).

2.2.3. Geomagnetic Sensor

The values obtained from the geomagnetic sensor require offset adjustment. When the values obtained from a three-axis geomagnetic sensor are plotted in three-dimensional space, they should ideally be distributed on the surface of a sphere or an ellipsoid. The center of this point cloud should correspond to the origin at $(x,y,z)=(0,0,0)$. However, if there are magnetic objects or iron products in the surrounding environment, the center may deviate significantly from the origin. Therefore, it is necessary to calculate the center of the point cloud obtained from the geomagnetic sensor and set it as the offset. To determine the offset, the geomagnetic sensor should be rotated in all directions while collecting and storing N sets of 3D point cloud data. In the following, the calculation method is explained, with the n-th 3D point cloud data for offset calculation represented as $(x_n,y_n,z_n)$.

Offset Adjustment

As mentioned earlier, the values obtained from the geomagnetic sensor are distributed on the surface of a sphere or an ellipsoid, making it necessary to calculate the center. The equation of an ellipsoid is expressed as follows:

$${\displaystyle \frac{x^2}{a^2}}+{\displaystyle \frac{y^2}{b^2}}+{\displaystyle \frac{z^2}{c^2}}=1.$$

(1)

In the Equation (1), when $a=b=c$, the shape becomes a sphere; in other words, if the center can be determined based on the equation of an ellipsoid, then the center of a sphere can also be obtained using the same approach.

In this study, the point cloud obtained from the geomagnetic sensor is fitted to an ellipsoid, then the center coordinates are calculated by determining the equation of the ellipsoid using the least squares method. While the equation of the ellipsoid is expressed as shown in Equation (1), this assumes that the x, y, and z axes are the principal axes. However, because the principal axes of the point cloud are often rotated, the following equation must be used instead:

$$A{x}^2+B{y}^2+C{z}^2+Dxy+Eyz+Fzx+Gx+Hy+Iz-K=0.$$

(2)

Here, because K does not affect any of the variables x, y, or z, both sides of Equation (2) are divided by K. The parameters obtained after division are re-expressed using the same symbols:

$$A{x}^2+B{y}^2+C{z}^2+Dxy+Eyz+Fzx+Gx+Hy+Iz-1=0$$

(3)

where J is defined as the least squares evaluation function. If we define the left-hand side of Equation (3) as $f(x,y,z)$, then $f(x_n,y_n,z_n)=0$ should hold for a point $(x_n,y_n,z_n)$ on the ellipsoid surface. In reality, there will be some deviation, meaning that it will not be exactly zero; nonetheless, we seek a constant that approaches zero as closely as possible for N points. Thus, the least squares evaluation function J is defined as in Equation (4).

$$J=\sum _{n=0}^N{f(x_n,y_n,z_n)}^2$$

(4)

To calculate this in such a way that the evaluation function J is minimized, it is necessary to solve the system of nine simultaneous equations formed by setting the partial derivatives concerning the nine parameters from A to I equal to zero in order to find the parameters. The equations obtained by partially differentiating the evaluation function J with respect to each parameter are shown below.

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial J}{\partial A}}& =\sum _{n=0}^{N}f(x,y,z)·x_n^2\hfill \\ \hfill & =\sum _{n=0}^N\left(A{x}_n^4+B{y}_n^2{x}_n^2+C{z}_n^2{x}_n^2+D{x}_n^3{y}_n+E{y}_n{z}_n{x}_n^2+F{z}_n{x}_n^3+G{x}_n^3+H{y}_n{x}_n^2+I{z}_n{x}_n^2-x_n^2\right)\hfill \\ \hfill & =0\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial J}{\partial B}}& =\sum _{n=0}^{N}f(x,y,z)·y_n^2\hfill \\ \hfill & =\sum _{n=0}^N\left(A{x}_n^2{y}_n^2+B{y}_n^4+C{z}_n^2{y}_n^2+D{x}_n{y}_n^3+E{y}_n^3{z}_n+F{z}_n{x}_n{y}_n^2+G{x}_n{y}_n^2+H{y}_n^3+I{z}_n{y}_n^2-y_n^2\right)\hfill \\ \hfill & =0\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial J}{\partial C}}& =\sum _{n=0}^{N}f(x,y,z)·z_n^2\hfill \\ \hfill & =\sum _{n=0}^N\left(A{x}_n^2{z}_n^2+B{y}_n^2{z}_n^2+C{z}_n^4+D{x}_n{y}_n{z}_n^2+E{y}_n{z}_n^3+F{z}_n^3{x}_n+G{x}_n{z}_n^2+H{y}_n{z}_n^2+I{z}_n^3-z_n^2\right)\hfill \\ \hfill & =0\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial J}{\partial D}}& =\sum _{n=0}^{N}f(x,y,z)·x_n{y}_n\hfill \\ \hfill & =\sum _{n=0}^N\left(A{x}_n^3{y}_n+B{y}_n^3{x}_n+C{z}_n^2{x}_n{y}_n+D{x}_n^2{y}_n^2+E{y}_n^2{z}_n{x}_n+F{z}_n{x}_n^2{y}_n+G{x}_n^2{y}_n+H{y}_n^2{x}_n+I{z}_n{x}_n{y}_n-x_n{y}_n\right)\hfill \\ \hfill & =0\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial J}{\partial E}}& =\sum _{n=0}^{N}f(x,y,z)·y_n{z}_n\hfill \\ \hfill & =\sum _{n=0}^N\left(A{x}_n^2{y}_n{z}_n+B{y}_n^3{z}_n+C{z}_n^3{y}_n+D{x}_n{y}_n^2{z}_n+E{y}_n^2{z}_n^2+F{z}_n^2{x}_n{y}_n+G{x}_n{y}_n{z}_n+H{y}_n^2{z}_n+I{z}_n^2{y}_n-y_n{z}_n\right)\hfill \\ \hfill & =0\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial J}{\partial F}}& =\sum _{n=0}^{N}f(x,y,z)·z_n{x}_n\hfill \\ \hfill & =\sum _{n=0}^N\left(A{x}_n^3{z}_n+B{y}_n^2{z}_n{x}_n+C{z}_n^3{x}_n+D{x}_n^2{y}_n{z}_n+E{y}_n{z}_n^2{x}_n+F{z}_n^2{x}_n^2+G{x}_n^2{z}_n+H{y}_n{z}_n{x}_n+I{z}_n^2{x}_n-z_n{x}_n\right)\hfill \\ \hfill & =0\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial J}{\partial G}}& =\sum _{n=0}^{N}f(x,y,z)·x_n\hfill \\ \hfill & =\sum _{n=0}^N\left(A{x}_n^3+B{y}_n^2{x}_n+C{z}_n^2{x}_n+D{x}_n^2{y}_n+E{y}_n{z}_n{x}_n+F{z}_n{x}_n^2+G{x}_n^2+H{y}_n{x}_n+I{z}_n{x}_n-x_n\right)\hfill \\ \hfill & =0\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial J}{\partial H}}& =\sum _{n=0}^{N}f(x,y,z)·y_n\hfill \\ \hfill & =\sum _{n=0}^N\left(A{x}_n^2{y}_n+B{y}_n^3+C{z}_n^2{y}_n+D{x}_n{y}_n^2+E{y}_n^2{z}_n+F{z}_n{x}_n{y}_n+G{x}_n{y}_n+H{y}_n^2+I{z}_n{y}_n-y_n\right)\hfill \\ \hfill & =0\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial J}{\partial I}}& =\sum _{n=0}^{N}f(x,y,z)·z_n\hfill \\ \hfill & =\sum _{n=0}^N\left(A{x}_n^2{z}_n+B{y}_n^2{z}_n+C{z}_n^3+D{x}_n{y}_n{z}_n+E{y}_n{z}_n^2+F{z}_n^2{x}_n+G{x}_n{z}_n+H{y}_n{z}_n+I{z}_n^2-z_n\right)\hfill \\ \hfill & =0\hfill \end{array}$$

(13)

From Equations (5) to (13), by moving the last term with a negative sign to the right-hand side, we obtain the system of simultaneous equations shown below.

$$\left\{\begin{array}{c}{\displaystyle \sum _{n=0}^N\left(A{x}_n^4+B{y}_n^2{x}_n^2+C{z}_n^2{x}_n^2+D{x}_n^3{y}_n+E{y}_n{z}_n{x}_n^2+F{z}_n{x}_n^3+G{x}_n^3+H{y}_n{x}_n^2+I{z}_n{x}_n^2\right)=\sum _{n=0}^N{x}_n^2}\hfill \\ {\displaystyle \sum _{n=0}^N\left(A{x}_n^2{y}_n^2+B{y}_n^4+C{z}_n^2{y}_n^2+D{x}_n{y}_n^3+E{y}_n^3{z}_n+F{z}_n{x}_n{y}_n^2+G{x}_n{y}_n^2+H{y}_n^3+I{z}_n{y}_n^2\right)=\sum _{n=0}^N{y}_n^2}\hfill \\ {\displaystyle \sum _{n=0}^N\left(A{x}_n^2{z}_n^2+B{y}_n^2{z}_n^2+C{z}_n^4+D{x}_n{y}_n{z}_n^2+E{y}_n{z}_n^3+F{z}_n^3{x}_n+G{x}_n{z}_n^2+H{y}_n{z}_n^2+I{z}_n^3\right)=\sum _{n=0}^N{z}_n^2}\hfill \\ {\displaystyle \sum _{n=0}^N\left(A{x}_n^3{y}_n+B{y}_n^3{x}_n+C{z}_n^2{x}_n{y}_n+D{x}_n^2{y}_n^2+E{y}_n^2{z}_n{x}_n+F{z}_n{x}_n^2{y}_n+G{x}_n^2{y}_n+H{y}_n^2{x}_n+I{z}_n{x}_n{y}_n\right)=\sum _{n=0}^N{x}_n{y}_n}\hfill \\ {\displaystyle \sum _{n=0}^N\left(A{x}_n^2{y}_n{z}_n+B{y}_n^3{z}_n+C{z}_n^3{y}_n+D{x}_n{y}_n^2{z}_n+E{y}_n^2{z}_n^2+F{z}_n^2{x}_n{y}_n+G{x}_n{y}_n{z}_n+H{y}_n^2{z}_n+I{z}_n^2{y}_n\right)=\sum _{n=0}^N{y}_n{z}_n}\hfill \\ {\displaystyle \sum _{n=0}^N\left(A{x}_n^3{z}_n+B{y}_n^2{z}_n{x}_n+C{z}_n^3{x}_n+D{x}_n^2{y}_n{z}_n+E{y}_n{z}_n^2{x}_n+F{z}_n^2{x}_n^2+G{x}_n^2{z}_n+H{y}_n{z}_n{x}_n+I{z}_n^2{x}_n\right)=\sum _{n=0}^N{z}_n{x}_n}\hfill \\ {\displaystyle \sum _{n=0}^N\left(A{x}_n^3+B{y}_n^2{x}_n+C{z}_n^2{x}_n+D{x}_n^2{y}_n+E{y}_n{z}_n{x}_n+F{z}_n{x}_n^2+G{x}_n^2+H{y}_n{x}_n+I{z}_n{x}_n\right)=\sum _{n=0}^N{x}_n}\hfill \\ {\displaystyle \sum _{n=0}^N\left(A{x}_n^2{y}_n+B{y}_n^3+C{z}_n^2{y}_n+D{x}_n{y}_n^2+E{y}_n^2{z}_n+F{z}_n{x}_n{y}_n+G{x}_n{y}_n+H{y}_n^2+I{z}_n{y}_n\right)=\sum _{n=0}^N{y}_n}\hfill \\ {\displaystyle \sum _{n=0}^N\left(A{x}_n^2{z}_n+B{y}_n^2{z}_n+C{z}_n^3+D{x}_n{y}_n{z}_n+E{y}_n{z}_n^2+F{z}_n^2{x}_n+G{x}_n{z}_n+H{y}_n{z}_n+I{z}_n^2\right)=\sum _{n=0}^N{z}_n}\hfill \end{array}\right.$$

(14)

In Equation (14), by expanding the left-hand side and factoring out the nine parameters A to I from the summation, it can be expressed as the following matrix using A to I. Note that the summation ${\sum }_{n=0}^N$ applies to all elements in the matrix shown below; however, it is omitted here for simplicity.

$$\left[\begin{array}{ccccccccc}{x}_n^4& y_n^2{x}_n^2& z_n^2{x}_n^2& x_n^3{y}_n& y_n{z}_n{x}_n^2& z_n{x}_n^3& x_n^3& y_n{x}_n^2& z_n{x}_n^2\\ x_n^2{y}_n^2& y_n^4& z_n^2{y}_n^2& x_n{y}_n^3& y_n^3{z}_n& z_n{x}_n{y}_n^2& x_n{y}_n^2& y_n^3& z_n{y}_n^2\\ x_n^2{z}_n^2& y_n^2{z}_n^2& z_n^4& x_n{y}_n{z}_n^2& y_n{z}_n^3& z_n^3{x}_n& x_n{z}_n^2& y_n{z}_n^2& z_n^3\\ x_n^3{y}_n& y_n^3{x}_n& z_n^2{x}_n{y}_n& x_n^2{y}_n^2& y_n^2{z}_n{x}_n& z_n{x}_n^2{y}_n& x_n^2{y}_n& y_n^2{x}_n& z_n{x}_n{y}_n\\ x_n^2{y}_n{z}_n& y_n^3{z}_n& z_n^3{y}_n& x_n{y}_n^2{z}_n& y_n^2{z}_n^2& z_n^2{x}_n{y}_n& x_n{y}_n{z}_n& y_n^2{z}_n& z_n^2{y}_n\\ x_n^3{z}_n& y_n^2{z}_n{x}_n& z_n^3{x}_n& x_n^2{y}_n{z}_n& y_n{z}_n^2{x}_n& z_n^2{x}_n^2& x_n^2{z}_n& y_n{z}_n{x}_n& z_n^2{x}_n\\ x_n^3& y_n^2{x}_n& z_n^2{x}_n& x_n^2{y}_n& y_n{z}_n{x}_n& z_n{x}_n^2& x_n^2& y_n{x}_n& z_n{x}_n\\ x_n^2{y}_n& y_n^3& z_n^2{y}_n& x_n{y}_n^2& y_n^2{z}_n& z_n{x}_n{y}_n& x_n{y}_n& y_n^2& z_n{y}_n\\ x_n^2{z}_n& y_n^2{z}_n& z_n^3& x_n{y}_n{z}_n& y_n{z}_n^2& z_n^2{x}_n& x_n{z}_n& y_n{z}_n& z_n^2\end{array}\right]·\left[\begin{array}{c}A\\ B\\ C\\ D\\ E\\ F\\ G\\ H\\ I\end{array}\right]=\left[\begin{array}{c}{x}_n^2\\ y_n^2\\ z_n^2\\ x_n{y}_n\\ y_n{z}_n\\ z_n{x}_n\\ x_n\\ y_n\\ z_n\end{array}\right]$$

(15)

In Equation (15), let the coefficient matrix be $\mathit{A}$, the matrix of parameters from A to I (the solution) be $\mathit{x}$, and the right-hand side matrix be $\mathit{b}$. Then, we have

$$\mathit{A}\mathit{x}=\mathit{b}.$$

(16)

By using the inverse matrix of $\mathit{A}$, we can solve for the parameters from A to I as follows:

$$\mathit{x}={\mathit{A}}^{-1}\mathit{b}.$$

(17)

If the center coordinates are $(x_0,y_0,z_0)$, then the parameters from A to I can be expressed as follows:

$$-\left[\begin{array}{ccc}2A& D& F\\ D& 2B& E\\ F& E& 2C\end{array}\right]·\left[\begin{array}{c}{x}_0\\ y_0\\ z_0\end{array}\right]=\left[\begin{array}{c}G\\ H\\ I\end{array}\right].$$

(18)

Thus, the center coordinates of the point cloud obtained from the geomagnetic sensor can be calculated as follows:

$$\left[\begin{array}{c}{x}_0\\ y_0\\ z_0\end{array}\right]=-{\left[\begin{array}{ccc}2A& D& F\\ D& 2B& E\\ F& E& 2C\end{array}\right]}^{-1}·\left[\begin{array}{c}G\\ H\\ I\end{array}\right].$$

(19)

2.3. Attitude Estimation

In conventional studies, only the $Yaw$ axis value necessary to calculate the number of steps is obtained; however, we always perform estimation in all axes of $Roll$, $Pitch$, and $Yaw$ using all of the acceleration, angular velocity, and geomagnetism sensors, and always obtain the direction in which the sensor is facing in the world coordinate system. This has the advantage that the acceleration in the $x,y$ direction can be separated from the acceleration in the z direction no matter what state the sensor is in, which facilitates the estimation procedure.

There are several methods for attitude estimation, including methods using the extended Kalman filter [21] and the Madgwick filter [22]. However, the extended Kalman filter is computationally expensive, while the Madgwick filter is slow in convergence. Both of these filters require time to obtain an accurate attitude. Therefore, we use the complementary filter in this study. In addition, because the use of Euler angles in calculating angles results in the appearance of many trigonometric functions and requires a great deal of computation time, calculations are performed using quaternions.

2.3.1. Complementary Filter Using Quaternions [23]

The complementary filter used in this study is based on the values estimated by the angular rate sensor and modified by the accelerometer and the geomagnetic sensor.

We define $\mathit{a}=(a_x,a_y,a_z)$ [${\mathrm{m/s}}^2$] as the value obtained from the acceleration sensor, $\mathit{\omega }=({\omega }_x,{\omega }_y,{\omega }_z)$ [rad/s] as the value obtained from the geomagnetic sensor, $\mathit{m}=(m_x,m_y,m_z)$ [T] as the value obtained from the geomagnetic sensor, $\mathit{R}\left(\mathit{q}\right)$ as the rotation matrix that rotates a vector by a quaternion q, $\mathrm{\Delta }T\left[s\right]$ as the difference between the previous processing time $t-1$ and the current processing time t, $\mathit{g}=9.798{\mathrm{m/s}}^2$ as the gravitational acceleration, ${\mathit{q}}_{\mathit{t}}={[q_{t,0},q_{t,1},q_{t,2},q_{t,3}]}^T$ as the quaternion representing the attitude at time t, and ${\mathit{q}}_{\mathit{t}-\mathbf{1}}={[q_{t-1,0},q_{t-1,1},q_{t-1,2},q_{t-1,3}]}^T$ as the quaternion representing the attitude at time $t-1$. The initial value for a quaternion is $\mathit{q}={[1,0,0,0]}^T$.

2.3.2. Attitude Estimation Using Angular Velocity Sensor

Let $\dot{\mathit{q}}={[{\dot{q}}_0,{\dot{q}}_1,{\dot{q}}_2,{\dot{q}}_3]}^T$ be the quaternion that changed during $\mathrm{\Delta }T$, which can be calculated as in Equation (20).

$$\left[\begin{array}{c}{\dot{q}}_0\\ {\dot{q}}_1\\ {\dot{q}}_2\\ {\dot{q}}_3\end{array}\right]={\displaystyle \frac{1}{2}}\left[\begin{array}{ccc}-q_{t-1,1}& -q_{t-1,2}& -q_{t-1,3}\\ q_{t-1,0}& -q_{t-1,3}& q_{t-1,2}\\ q_{t-1,3}& q_{t-1,0}& -q_{t-1,1}\\ -q_{t-1,2}& q_{t-1,1}& q_{t-1,0}\end{array}\right]·\left[\begin{array}{c}{\omega }_x\\ {\omega }_y\\ {\omega }_z\end{array}\right]·\mathrm{\Delta }T$$

(20)

The attitude ${\mathit{q}}_{\mathit{t}}$ at time t can be calculated as in Equation (21).

$$\left[\begin{array}{c}{q}_{t,0}\\ q_{t,1}\\ q_{t,2}\\ q_{t,3}\end{array}\right]=\left[\begin{array}{c}{q}_{t-1,0}\\ q_{t-1,1}\\ q_{t-1,2}\\ q_{t-1,3}\end{array}\right]+\left[\begin{array}{c}{\dot{q}}_0\\ {\dot{q}}_1\\ {\dot{q}}_2\\ {\dot{q}}_3\end{array}\right]$$

(21)

2.3.3. Correction of Attitude Estimates Using an Acceleration Sensor

Using the values obtained from the acceleration sensor, the absolute angles of the $Pitch$ and $Roll$ axes with respect to the world coordinates can be obtained. This absolute angle is used to correct the attitude estimation.

Rotating the acceleration sensor value based on the attitude calculated from the angular velocity sensor should always result in $(a_x,a_y,a_z)=(0,0,g)$. For example, if the sensor is rotated 90 degrees around the x-axis in world coordinates, then the sensor value is $(a_x,a_y,a_z)=(0,g,0)$. If we know that the current attitude is rotated 90 degrees around the x-axis, then rotating the sensor value by $-90$ degrees around the x-axis in the opposite direction will result in $(a_x,a_y,a_z)=(0,0,g)$. However, in reality the estimated attitude by the angular velocity sensor deviates slightly from $(0,0,g)$ because of a small error in the attitude estimated by the angular rate sensor. Letting $\mathit{G}=(G_x,G_y,G_z)$ be the vector on which the error exists, it can be calculated as in Equation (22):

$$\mathit{G}=\mathit{R}\left({\mathit{q}}_{\mathit{t}}\right)\mathit{a}.$$

(22)

Then, by finding the quaternion to rotate from $\mathit{G}$ to $(0,0,g)$ and multiplying it by ${\mathit{q}}_{\mathit{t}}$ and rotating it, we can calculate the attitude corrected by the acceleration sensor. The quaternions that rotate from $\mathit{G}$ to $(0,0,g)$ can be calculated using inner and outer products. We define the reference $(0,0,g)$ vector as ${\mathit{a}}_{\mathbf{0}}=(0,0,1)$, the rotation axis unit vector representing the direction of rotation in the quaternion as ${\mathit{a}}_{\mathit{n}}=(an,x,an,y,an,z)$, and the scalar representing the rotation angle as ${\mathit{\theta }}_{\mathit{a}}$. Then, ${\mathit{a}}_{\mathit{n}}$ and ${\mathit{\theta }}_{\mathit{a}}$ are as follows:

$${\mathit{a}}_{\mathit{n}}={\displaystyle \frac{\mathit{G}\times {\mathit{a}}_{\mathbf{0}}}{|\mathit{G}\times {\mathit{a}}_{\mathbf{0}}|}},$$

(23)

$${\theta }_a=arccos\left(\mathit{G}·{\mathit{a}}_{\mathbf{0}}\right).$$

(24)

The corrected quaternion ${\mathit{q}}_{\mathit{a}}$ due to the acceleration can be calculated as in Equation (25).

$$\begin{array}{cc}\hfill \left[\begin{array}{c}{q}_{a,0}\\ q_{a,1}\\ q_{a,2}\\ q_{a,3}\end{array}\right]& =\left[\begin{array}{c}cos{\displaystyle \frac{{\theta }_a}{2}}\\ a_{n,x}sin{\displaystyle \frac{{\theta }_a}{2}}\\ a_{n,y}sin{\displaystyle \frac{{\theta }_a}{2}}\\ a_{n,z}sin{\displaystyle \frac{{\theta }_a}{2}}\end{array}\right]\hfill \\ \hfill & =\left[\begin{array}{c}\sqrt{{\displaystyle \frac{G_z+1}{2}}}\\ {\displaystyle \frac{G_y}{\sqrt{2(G_z+1)}}}\\ -{\displaystyle \frac{G_x}{\sqrt{2(G_z+1)}}}\\ 0\end{array}\right]\hfill \end{array}$$

(25)

As shown in Equation (25), the z-axis component is 0, i.e., only the corrections to the $Pitch$ and $Roll$ axes are calculated. If the obtained correction quaternion ${\mathit{q}}_{\mathit{a}}$ is used as it is, it will be affected by all the effects of the noise from the acceleration sensor. Thus, it is necessary to use linear interpolation with the non-rotating quaternion ${\mathit{q}}_{\mathit{I}}=(1,0,0,0)$. We define $\alpha$ as the value of the correction strength of the correction quaternion and set a new correction quaternion ${\mathit{q}}_{\mathit{a}}^{\prime }$.

$${\mathit{q}}_{\mathit{a}}^{\prime }=\left[\begin{array}{c}\alpha q_{a,0}+(1-\alpha )\\ \alpha q_{a,1}\\ \alpha q_{a,2}\\ \alpha q_{a,3}\end{array}\right]$$

(26)

Although $\alpha$ can be a fixed value, in this study the correction intensity is adjusted according to the magnitude of the vectors obtained from the acceleration sensor. If the sensor is not accelerating or decelerating, the magnitude $\sqrt{a_x^2+a_y^2+a_z^2}$ of the acceleration vector obtained from the sensor should be equal to g. Therefore, we define the correction strength $\alpha$ as

$$\alpha ={\alpha }_0·f\left(e_g\right),$$

(27)

where ${\alpha }_0$ is the fixed correction intensity and $f\left(e_g\right)$ is a function of the error with respect to the acceleration of gravity. We define $e_g$ and $f\left(x\right)$ as follows.

$$e_g={\displaystyle \frac{\left|\left|\left|\mathit{a}\right|\right|-g\right|}{g}}$$

(28)

$$f\left(x\right)=\left\{\begin{array}{cc}1\hfill & (x\le e_1)\hfill \\ {\displaystyle \frac{e_2-x}{e_2-e_1}}\hfill & (e_1<x\le e_2)\hfill \\ 0\hfill & (x\ge e_2)\hfill \end{array}\right.$$

(29)

It is necessary to conduct experiments to obtain appropriate values for ${\alpha }_0$ and $e_1,e_2$.

From the above, the quaternion representing the attitude corrected by the acceleration sensor can be calculated as in Equation (30):

$${\mathit{q}}_{\mathit{t}}^{\prime }={\mathit{q}}_{\mathit{a}}^{\prime }{\mathit{q}}_{\mathit{t}}.$$

(30)

2.3.4. Correction of Attitude Estimates Using a Geomagnetic Sensor

Using the values obtained from the geomagnetic sensor, the absolute angles of the $Yaw$ axes with respect to the world coordinates can be obtained.

By rotating the geomagnetic sensor values based on the attitude calculated from the acceleration and angular velocity, we can obtain the same values $\mathit{l}=(l_x,l_y,l_z)$ when the sensor is horizontal.

$$\mathit{l}=\mathit{R}\left({\mathit{q}}_{\mathit{t}}^{\prime }\right)\mathit{m}$$

(31)

Considering $l_x,l_y$, we can calculate the error concerning the direction at the beginning of the process. If the sensor was facing north at the beginning of the process, then the error concerning north can be obtained. As in the case of the accelerometer, the corrected quaternion ${\mathit{q}}_{\mathit{m}}$ from the geomagnetic sensor can be obtained as in Equation (32).

$${\mathit{q}}_{\mathit{m}}=\left[\begin{array}{c}\sqrt{{\displaystyle \frac{(l_x^2+l_y^2)+l_x\sqrt{l_x^2+l_y^2}}{2(l_x^2+l_y^2)}}}\\ 0\\ 0\\ {\displaystyle \frac{l_y}{\sqrt{2\left\{(l_x^2+l_y^2)+l_x\sqrt{l_x^2+l_y^2}\right\}}}}\end{array}\right]$$

(32)

If the obtained correction quaternion ${\mathit{q}}_{\mathit{m}}$ is used as it is, it will be affected by all the effects of the noise from the geomagnetic sensor. Thus, it is necessary to use linear interpolation with the non-rotating quaternion ${\mathit{q}}_{\mathit{I}}=(1,0,0,0)$. We define $\beta$ as the value of the correction strength of the correction quaternion and set a new correction quaternion ${\mathit{q}}_{\mathit{m}}^{\prime }$.

$${\mathit{q}}_{\mathit{m}}^{\prime }=\left[\begin{array}{c}\beta q_{m,0}+(1-\beta )\\ \beta q_{m,1}\\ \beta q_{m,2}\\ \beta q_{m,3}\end{array}\right]$$

(33)

From the above, the quaternion representing the attitude corrected by the geomagnetic sensor can be calculated as in Equation (34):

$${\mathit{q}}_{\mathit{t}}^{″}={\mathit{q}}_{\mathit{m}}^{\prime }{\mathit{q}}_{\mathit{t}}^{\prime }$$

(34)

where ${\mathit{q}}_{\mathit{t}}^{″}$ is the quaternion representing the corrected attitude in all axes, i.e., $Pitch$, $Roll$, and $Yaw$.

2.3.5. Dealing with Geomagnetic Sensor Disturbances

The value of the geomagnetic sensor is affected by magnets or other magnetizing objects contained in other devices or worn on the body. While disturbances are often not considered in conventional studies, they must be taken into account in practical applications. Because the geomagnetic sensor is the only sensor that can obtain the direction in world coordinates, it is important to deal with disturbances to estimate the direction of motion.

In this study, the value of $\beta$ in Equation (33) is changed according to the magnitude of the error of the value obtained from the geomagnetic sensor compared with the normal condition. In addition, to accommodate changes in offset values caused by variations in surrounding structures, updates to the offset adjustment described in the geomagnetic sensor section (Section 2.2.3) are also performed. However, this update is not calculated when the value of $\beta$ in Equation (35) is 0, as it is then affected by other devices.

We redefine the correction strength $\beta$ as follows:

$$\beta ={\beta }_0·f\left(e_m\right)$$

(35)

where ${\beta }_0$ is the fixed correction intensity and $f\left(e_m\right)$ is a function of the error concerning the geomagnetic sensor value $M_a$ under normal conditions. We define $e_m$ and $f\left(x\right)$ as follows.

$$e_m=\sqrt{l_x^2+l_y^2+l_z^2}-M_a$$

(36)

$$f\left(x\right)=e^{-c_1{x}^2}$$

(37)

It is necessary to conduct experiments to obtain appropriate values for ${\beta }_0$ and $c_1$.

2.4. Step Detection

Step detection is performed based on the acceleration, angular velocity, and attitude estimation results. Previous studies have commonly used only the z-axis acceleration for estimation. However, because the magnitude of acceleration changes frequently depending on the way of walking, this introduces the possibility of a large error in the number of steps. Therefore, in this study estimation is performed using not only the acceleration along the z-axis but also the acceleration along the $x,y$-axis, walking interval, and magnitude of angular velocity. Based on the results of attitude estimation, we determine whether the subject is engaged in hand-held walking or hand-swinging walking. If the $Roll$ axis (y-axis) is tilted by more than a certain angle, the user is judged to be engaged in hand-swinging walking; otherwise, they are judged to be engaged in hand-held walking. If the angular velocity around the $Roll$ axis is larger than a certain value, this is considered to indicate a switch between hand-held walking and hand-swinging walking, and the number of steps is not counted.

2.5. Step Length Estimation

Step length estimation is calculated based on how much larger or smaller the stride length is compared to the reference stride length. When estimating the step length, it is necessary to determine the magnitude of the acceleration compared to that of normal walking. Let $A_z$ be the acceleration in the direction of gravitational acceleration and let $A_x$ and $A_y$ be the acceleration in the horizontal direction. These accelerations are a rotated value of the acceleration sensor value $(a_x,a_y,a_z)$ based on the result of attitude estimation. We define the maximum value of $A_z$ during step counting as $A_{z,max}$ and the minimum value as $A_{z,min}$. In addition, we define the norm of the maximum absolute value $A_{x,max},A_{y,max}$ of $A_x,A_y$ as $A_{xy,max}$. The magnitude of acceleration $A_{ave}$ is obtained using the average, defined as follows:

$$A_{ave}={\displaystyle \frac{(A_{z,max}-A_{z,min})+A_{xy,max}}{2}}.$$

(38)

In addition, because we need to compare the acceleration to that of normal walking, we redefine it so that the magnitude of the acceleration during normal walking is 0:

$$A_{adjust}=A_{ave}-C_{normal}$$

(39)

where $C_{normal}$ is a constant that must be obtained by calculating the magnitude of acceleration during normal walking by experiment.

Let l be the estimated step length and $l_s$ the reference step length; then, l can be calculated by the following Equation (40):

$$l=l_s+K_l·g\left(A_{adjust}\right).$$

(40)

The reference step length $l_s$ is calculated by the following Equation (41) [4] when the height is h [cm]:

$$l_s=h-100$$

(41)

where $K_l$ is the maximum value of the step length adjustment, e.g., $K_l=30$ cm, and the function $g\left(x\right)$ changes the value from $-30$ cm to $+30$ cm, then is added to the standard step length $l_s$ to compute l. The function $g\left(x\right)$ is defined as follows.

$$g\left(x\right)=\left\{\begin{array}{cc}1\hfill & (x\ge 1)\hfill \\ x^3\hfill & (-1<x<1)\hfill \\ -1\hfill & (x\le -1)\hfill \end{array}\right.$$

(42)

By cubing, even if there is some error in $A_{adjust}$, it will be small if the error is not large.

From Equations (40)–(42), l can be calculated as in Equation (43):

$$l=(h-1)+K_l·g\left({\displaystyle \frac{(A_{z,max}-A_{z,min})+A_{xy,max}}{2}}-C_{normal}\right).$$

(43)

Therefore, it is possible to estimate the step length by experimentally obtaining the constants $K_l$ and $C_{normal}$ for the cases of hand-held walking and hand-swinging walking, respectively.

2.6. Step Direction Estimation

Step direction estimation is performed based on the step detection and attitude estimation results. For both hand-held walking and hand-swinging walking, the direction in which the sensor was facing at the moment of the maximum acceleration during the step count is taken as the direction of the step. In the case of hand-held walking, the direction of the step is the direction when $A_{diff,max}$ is measured, while in the case of hand-swinging walking the direction of the step is the direction when $A_{z,max}$ is measured. In this case, the direction is the angle of rotation around the $Yaw$ axis in world coordinates, i.e., the direction with respect to north.

2.7. Position Estimation

Position estimation is performed based on the results of step detection, step length estimation, and step direction estimation. Processing is performed when the number of steps is counted, and the position is computed by adding up the results based on the step length obtained from the step length estimation and the direction of travel obtained from the step direction estimation. We define the position at the count one step ago as ${\mathit{P}}_{\mathit{t}-\mathbf{1}}={(P_{t,x},P_{t,y})}^T$, the current position as ${\mathit{P}}_{\mathit{t}}={(P_{t-1,x},P_{t-1,y})}^T$, the estimated step length as l, and the step direction as $\theta$. Then, ${\mathit{P}}_{\mathit{t}}$ can be obtained by the following Equation (44).

$$\left[\begin{array}{c}{P}_{t,x}\\ P_{t,y}\end{array}\right]=\left[\begin{array}{c}{P}_{t-1,x}\\ P_{t-1,y}\end{array}\right]+\left[\begin{array}{c}lcos\theta \\ lsin\theta \end{array}\right]$$

(44)

Note that ${\mathit{P}}_{\mathit{t}}$ is the position where y-axis positive is north and x-axis positive is east. This can be changed in any direction by adding offsets.

3. Results

3.1. Summary

3.1.1. Device

The device used in this study was a Google Pixel 4 with Android version 12. Each estimation was performed by an application developed using Unity.

3.1.2. Experiment Goals

The goals of this experiment were threefold: to estimate walking trajectories using only sensors embedded in mobile devices, to handle geomagnetic disturbances, and to achieve accurate estimation in real-world environments. These goals correspond to Section 3.4.1, Section 3.4.2, and Section 3.4.3, respectively.

3.1.3. Implementation Method

Section 3.4.1 described the experiments to confirm the feasibility of trajectory estimation. These experiments were conducted in a relatively disturbance-free open space. The estimation was performed for both hand-held walking and hand-swinging walking scenarios. Hand-held walking refers to a state where the device is held in one hand while the user is looking at the screen, whereas hand-swinging walking refers to a state where the user holds the device in one hand with the arm relaxed and swinging back and forth while walking. Section 3.4.2 addresses estimation in scenarios where magnetic disturbances from other devices occur, simulating a crowded environment. Section 3.4.3 focuses on indoor positioning scenarios and performing estimation along complex trajectories in environments with numerous household appliances and other devices. In each experiment, walking is performed continuously without stopping at each step.

3.2. Result of Adjustment

3.2.1. Acceleration Sensor

The offset values were determined using the average of 500 data points obtained when the sensor was placed on a horizontal plane.

Table 1. The average of 500 data points from the accelerometer when placed on a horizontal plane.

	${\mathit{ax}}_{\mathit{ave}}$	${\mathit{ay}}_{\mathit{ave}}$	${\mathit{az}}_{\mathit{ave}}$
Average [${\mathrm{m/s}}^2$]	0.0468	0.0295	9.9853

shows the average acceleration of the 500 data points obtained when the sensor was placed on a horizontal plane.

From Table 1, as calculated and explained in Section 2.2.1, $(a{x}_{offset},a{y}_{offset},a{z}_{offset})=(-0.0468,-0.0295,-0.1873)$ was determined.

The coefficient values were determined using the average of 500 data points obtained when each x, y, and z axis of the sensor was oriented vertically upward and vertically downward relative to the horizontal plane.

Table 2. The average of 500 data points from the accelerometer for each x, y, and z axis when oriented vertically upward and downward.

	${\mathit{ax}}_{\mathit{ave}}$		${\mathit{ay}}_{\mathit{ave}}$		${\mathit{az}}_{\mathit{ave}}$
	Upward	Downward	Upward	Downward	Upward	Downward
Average [${\mathrm{m/s}}^2$]	9.7850	−9.7914	9.7922	−9.7807	9.7352	−9.7143

shows the average acceleration of 500 data points for each x, y, and z axis when oriented vertically upward and vertically downward.

From Table 2, as calculated and explained in Section 2.2.1, $(a{x}_{P_{coeff}},a{x}_{N_{coeff}})=(1.0013,1.0007),(a{y}_{P_{coeff}},a{y}_{N_{coeff}})=(1.0006,1.0018),(a{z}_{P_{coeff}},a{z}_{N_{coeff}})=(1.0065,1.0086)$ was determined.

3.2.2. Angular Velocity Sensor

The offset values were determined using the average of 500 data points obtained when the sensor was placed on a flat surface and kept stationary.

Table 3. The average of 500 data points from the angular velocity sensor when kept stationary.

	${\mathit{gx}}_{\mathit{ave}}$	${\mathit{gy}}_{\mathit{ave}}$	${\mathit{gz}}_{\mathit{ave}}$
Average  ${10}^{-3}$ [rad/s]	−2.672	−3.499	8.724

shows the average angular velocity of 500 data points obtained when the sensor was placed on a flat surface and kept stationary.

From Table 3, as calculated and explained in Section 2.2.2, $(g{x}_{offset},g{y}_{offset},g{z}_{offset})=(0.002672,0.003499,-0.008724)$ was determined.

3.2.3. Geomagnetic Sensor

The offset values were determined by calculating the center from 3D point cloud data using the ellipsoid center calculation method based on the least squares method, as described in Section 2.2.3. As explained in Section 2.3.5, it should be noted that the offset is updated over time; therefore, the following data represent test data used to verify whether the calculations were correctly implemented. The following figure shows the 3D point cloud data obtained from the geomagnetic sensor before adjustment.

Because it is difficult to see in a 3D graph, the 3D point cloud data viewed from the x-y plane, x-z plane, and y-z plane are shown below.

As shown in Figure 2 and Figure 3, the values form a sphere displaced from the origin. As a result of using the least squares method to calculate the center, the coordinates $(x_0,y_0,z_0)=(39.6975,-69.0510,105.0062)$ were obtained.

As explained in Section 2.2.3, the offset values were calculated and determined to be $(m{x}_{offset},m{y}_{offset},m{z}_{offset})=(-39.6975,69.0510,-105.0062)$.

Using the obtained offset values, when the values in Figure 2 and Figure 3 are corrected, the result is as shown in Figure 4 and Figure 5, confirming that the offset values have been correctly determined.

3.3. Results and Evaluation of Attitude Estimation

In order to use the complementary filter, it is necessary to find appropriate values for four parameters: three parameters ${\alpha }_0,e_1,e_2$, and a parameter $\beta$ for finding parameters ${\beta }_0,c_1$ for magnetic disturbance. The ${\alpha }_0,e_1,e_2$ parameters are only related to $Roll,Pitch$ correction, while $\beta$ is only related to $Yaw$ correction. Therefore, we evaluated ${\alpha }_0,e_1,e_2$ using the time until $Roll,Pitch$ becomes stable and evaluated $\beta$ using the time until $Yaw$ becomes stable. In the complementary filter, increasing ${\alpha }_0$ means increasing the considered value of the current acceleration. Because the acceleration should be 0 at the moment of stationary, larger ${\alpha }_0$ indicates shorter time to stability. However, a larger considered value of the current acceleration value means that the estimation error during acceleration and deceleration becomes larger. This point can be somewhat mitigated by setting appropriate values of $e_1,e_2$; however, if ${\alpha }_0$ is set too large, then the error is considered to be unacceptable. Hence, we set a new criterion for evaluating the change in attitude during acceleration and deceleration based on the fact that no rotation around the $Roll,Pitch$ axis occurs when the horizontal plane is moved. Figure 6 shows the results of $Roll,Pitch$ estimation when the sensor is placed on the horizontal plane and moved back and forth.

From Figure 6, it can be seen that rotation around the $Roll,Pitch$ axis should not occur when the sensor is moved back and forth in the horizontal plane; however, the error is caused by increasing ${\alpha }_0$, which is the considered value of the current acceleration value. In finding the values of $e_1,e_2$, we calculate the norm $R{P}_{norm}$ of the maximum value of the rotation around $Roll,Pitch$ axis in the round-trip movement and find $e_1,e_2$ such that $R{P}_{norm}$ becomes as small as possible.

Table 4. Average of $R{P}_{norm}$ when $e_1$ of the complementary filter is varied.

${\mathit{e}}_1$	${\mathit{RP}}_{\mathit{norm}}$ [rad]
0.01	0.282
0.001	0.222
0.0001	0.212

shows the values of $R{P}_{norm}$ when $e_1$ is varied with ${\alpha }_0=0.4,e_2=0.1$. Note that $R{P}_{norm}$ is the average of ten times.

From Table 4, it can be seen that the error $R{P}_{norm}$ becomes smaller by decreasing $e_1$.

Table 5. Average of $R{P}_{norm}$ when $e_2$ of the complementary filter is varied.

${\mathit{e}}_1$	${\mathit{RP}}_{\mathit{norm}}$ [rad]
0.1	0.212
0.01	0.179
0.001	0.073

shows the values of $R{P}_{norm}$ when $e_2$ is varied with ${\alpha }_0=0.4,e_1=0.0001$. Note that $R{P}_{norm}$ is the average of ten times.

From Table 5, it is found that the error $R{P}_{norm}$ becomes smaller by decreasing $e_2$. The smaller $e_1,e_2$ are, the smaller $R{P}_{norm}$ is; however, if $e_2$ is set too small from Equation (29), the value of acceleration is hardly taken into account during acceleration, and a large error is considered to occur in the estimation results during acceleration. Therefore, we use $e_2=0.01$ and set $(e_1,e_2)=(0.0001,0.01)$.

Table 6. Average of $R{P}_{norm}$ when ${\alpha }_0$ of the complementary filter is varied.

${\mathit{\alpha }}_0$	${\mathit{RP}}_{\mathit{norm}}$ [s]
0.1	0.072
0.2	0.103
0.3	0.166
0.4	0.179

shows the ten averages of $R{P}_{norm}$ when $(e_1,e_2)=(0.0001,0.01)$ is fixed and ${\alpha }_0$ is varied.

From Table 6, it can be seen that the error $R{P}_{norm}$ becomes smaller by decreasing ${\alpha }_0$. However, the larger ${\alpha }_0$ is, the shorter the time to stability becomes; thus, we set ${\alpha }_0=0.2$ to improve the estimation results for both $R{P}_{norm}$ and $t_{stab}$. From the above, the parameters for the correction of $Roll,Pitch$ are determined to be ${\alpha }_0=0.2,(e_1,e_2)=(0.0001,0.01)$.

Next, we evaluated how the change in $\beta$ affects the stability and obtained the parameters. Figure 7 shows the results of attitude estimation when $\beta =0.01,0.2$ and the sensor is kept stationary.

As can be seen from Figure 7 (Right), the estimated value is not stable even after the sensor becomes stationary when $\beta$ was increased.

Therefore, the parameter for the correction of $Yaw$ is set to $\beta =0.01$.

Next, the parameters ${\beta }_0,c_1$ for coping with magnetic disturbance are obtained. The value of Equation (37) takes a value from 0 to 1, and the final maximum value of Equation (35) must be $\beta =0.01$; thus, we set ${\beta }_0=0.01$. The value of $c_1$ should be determined by the magnitude of the difference between the measured geomagnetic field under normal conditions and that under disturbance. In the following

Table 7. Geomagnetic readings when other devices are close by.

Distance [cm]	Devices
	Device Affecting the Geomagnetic Field (iPhone)		Device Affecting the Geomagnetic Field (AppleWatch)
	Measured Value [μT]	Absolute Error [μT]	Measured Value [μT]	Absolute Error [μT]
40	45.3	0.05	45.2	0.11
30	45.3	0.04	45.3	0.03
20	45.8	0.45	45.2	0.15
10	55.2	9.85	48.6	3.22
0	2395.3	2349.95	215.2	169.85

, we show the measured values when an iPhone and an AppleWatch are placed close to the measured device. The measured values are averaged over 10 s. The measured geomagnetic field under normal conditions was averaged for 100 s, and the result of the measurement was $45.349\mathrm{\mu }$ T.

From Table 7, the error is found to occur at a distance as close as 20 cm, and becomes very large when the distance is closer than 10 cm. In practical applications, smartwatches and similar devices may be closer than 10 cm; thus, errors that occur in such cases must be dealt with. In order to avoid large errors in the estimates, we do not consider the geomagnetic values when the absolute error is about $3\mathrm{\mu }$ T. From the shape of the graph of $e^{-a{x}^2}$, the value of $c_1$ is set to about $0.8$.

From the above results, the parameters of the correction in the complementary filter are determined as follows: ${\alpha }_0=0.2,(e_1,e_2)=(0.0001,0.01),{\beta }_0=0.01,c_1=0.8$.

To evaluate how the geomagnetic sensor copes with disturbances, we compare the results between the cases with and without a geomagnetic sensor value as well as the cases with a fixed value of $\beta$ and a variable value of $\beta$. Because the only error that can occur in each case is the $Yaw$ value, only the $Yaw$ value is considered for the error assessment. Each estimated value is the error from the initial state when the sensor is rotated in various directions for 10 s and then stopped.

The following

Table 8. Error in $Yaw$ axis estimation with and without geomagnetic sensor.

	Absolute Error When Not Used [deg]	Absolute Error When Used [deg]
1	10.51	0.47
2	15.21	0.65
3	11.42	0.13
4	10.15	0.35
5	13.32	0.08
6	14.23	0.12
7	10.57	0.88
8	10.91	0.43
9	11.25	0.12
10	12.89	0.78
Ave.	12.04	0.40

shows the error values of $Yaw$ axis estimation without and with the geomagnetic sensor in addition to the angular velocity and acceleration sensors.

From Table 8, the errors are within 1 degree in each case when the geomagnetic sensor is used, whereas the errors are more than about 10 degrees in the case when the geomagnetic sensor is not used. This indicates that the geomagnetic sensor is important for attitude estimation.

The following

Table 9. Error in $Yaw$ axis estimation with and without change in $\beta$.

	Absolute Error at Fixed $\mathit{\beta }$ [deg]	Absolute Error When $\mathit{\beta }$ Changes Dynamically [deg]
1	96.26	1.78
2	73.12	0.63
3	102.66	0.82
4	80.29	0.91
5	161.81	1.21
6	89.11	1.98
7	133.12	1.03
8	148.99	0.91
9	103.31	1.87
10	129.73	0.22
Average	111.84	1.13

shows the error of the $Yaw$ axis estimation with the AppleWatch in close contact both when the value of $\beta$ is fixed and when it is changed dynamically to cope with the disturbance.

Table 9 shows that the error is within 2 degrees when $\beta$ is changed dynamically, while the error is more than about 100 degrees when $\beta$ is fixed. Thus, the dynamic change of $\beta$ is very important to cope with the disturbance. Although the error is a little larger than in the case without disturbance (Table 8), if $c_1$ is set such that the geomagnetic field is not taken into account when there is even a small deviation from the normal geomagnetic field value, then even a small change in value due to changes in the surrounding environment will cause the geomagnetic field to not be considered. Therefore, we believe that the value set in this study is appropriate. If the error between the geomagnetic field and the normal geomagnetic field value becomes small, then the geomagnetic field value is referred to again to correct the correct attitude; thus, the error in Table 9 is not considered to be a problem.

Various estimations were performed using the complementary filters obtained from the above results.

3.4. Results and Evaluation of Trajectory Estimation

3.4.1. Trajectory Estimation in Environments with Minimal Disturbances

The evaluation trajectories consisted of a rectangular trajectory with $(x,y)=(5,10)$ [m] and a circular trajectory with a radius of 3 m. These are shown in the figures below. This environment is on a ground floor with good visibility, providing an environment with minimal magnetic interference. The building is constructed with reinforced concrete.

The evaluation is conducted using the estimated trajectory during walking along the same rectangular and circular trajectories, as shown in Figure 8. The evaluation is based on the difference between the starting point and the endpoint. This difference is calculated as the distance between the final point of the estimated trajectory and the origin, as the origin is used as the starting point for the evaluation. Figure 9 and Figure 10 show the estimation results for hand-held walking and hand-swinging walking in rectangular trajectory, and Figure 11 and Figure 12 show the estimation results for hand-held walking and hand-swinging walking in circular trajectory.

The error between the starting point and the endpoint is shown in

Table 10. Absolute error between the starting point and endpoint during position estimation.

	Rectangular Trajectory		Circular Trajectory
	First Time (L)	Second Time (R)	First Time (L)	Second Time (R)
Absolute error of hand-held [cm]	70.1	110.5	101.7	115.3
Absolute error of hand-swinging [cm]	46.0	138.7	135.0	57.9

.

3.4.2. Trajectory Estimation in Environments with Significant Geomagnetic Disturbances

To evaluate whether the influence of geomagnetic fields from other devices is suppressed, an estimation was performed when other devices were located near the walking trajectory. Figure 13 below shows the actual walking trajectory and the positions of the other devices. Note that the other devices were installed at the same height as the measurement device and that the distance was varied from the walking trajectory at 0.5 m, 0.3 m, and 0.2 m. Walking started from the origin, passing beside Other Device 1 on the right and Other Device 2 on the left. To avoid collisions with the other devices, the minimum distance was set to 0.2 m. Other Device 1 was an iPhone 12 and Other Device 2 was an Apple Watch (Series 5). The estimation results shown in Figure 14 were obtained by simultaneously estimating with (proposed method) and without (conventional method) dynamically changing $\beta$ while using exactly the same sensor values.

From Figure 14, it can be seen that there is no significant error when the distance is about $0.5$ m. However, when the distance is reduced to below $0.3$ m, the estimated trajectory deviates. In the case of $0.2$ m, when $\beta$ is fixed, there is a deviation of $0.258$ m in the x-direction; when $\beta$ is dynamically adjusted, the deviation is reduced to $0.083$ m. Although there were only two other devices within a 4 m range in this experiment, in a real-world indoor environment there are likely to be even more devices passing by. Thus, the effectiveness of the proposed method in improving the estimation results is expected to be significant.

3.4.3. Trajectory Estimation in Real-World Environments

To evaluate whether estimation is feasible in an actual indoor environment, walking trajectory estimation was conducted during hand-held walking indoors. The environment used for estimation included desks, chairs, PCs, TVs, refrigerators, etc., inside a wooden structure. The following Figure 15 shows the actual walking trajectory and the environment. The arrows indicate the order in which the trajectory is walked. By following the arrows that are not opposite to the direction of movement, return to the starting point. The thin black lines in the figure are spaced at 50 cm intervals, the thick black lines represent walls, and the gray areas indicate places where desks or shelves are located so as to prevent walking.

Figure 16 shows the estimated results obtained by walking along the trajectory shown in Figure 15. In the figure, “S” represents the starting point and “G” represents the estimated position when returning. The estimation results shown in Figure 16 were obtained by simultaneously estimating with (proposed method) and without (conventional method) dynamically changing $\beta$ while using exactly the same sensor values.

From the results shown in Figure 16, it is evident that preventing geomagnetic disturbances caused by other devices and equipment improves the accuracy of walking trajectory estimation. This is probably because the geomagnetic values vary slightly from room to room depending on the other devices, leading to the direction gradually deviating if this variation is not properly addressed.

4. Conclusions

In this study, we used an acceleration sensor, angular velocity sensor, and geomagnetic sensor to estimate attitude, step count, step length, step direction, and position. As a result of our experiments, we conclude the following:

• For attitude estimation, our proposed approach is able to cope with disturbance of the geomagnetic sensor and counteract the influence of other devices when estimating the attitude in a real environment.
• For step detection, although some errors may occur at the beginning and the end of walking, the estimation procedure was sufficiently accurate during walking.
• For step length estimation, the error was limited to about $\pm 10$ cm at the maximum, and the accuracy was about $\pm 5$ cm under normal conditions.
• For step direction estimation, the error during hand-swinging walking was about ±6 degrees, while the error in hand-held walking was estimated with an accuracy of about ±2 degrees.
• For position estimation, although it was not possible to estimate the position precisely, it was possible to estimate the position with an accuracy of about $1.5$ m for the trajectory evaluated in this study. In addition, the shape of the trajectory could be discriminated sufficiently.
• For the position estimation with disturbances from other devices, although the influence could not be completely eliminated, the error in the estimation results was significantly reduced.

In conclusion, we found that trajectory estimation is approximately feasible for hand-held walking and hand-swinging walking, and we believe that our results are a sufficient, if not perfect, solution for the problem of unconstrained self-position estimation.

The following are among of the issues to be addressed in future research:

• Algorithm for determining the number of steps taken at the beginning and end of a walk.
• Modification of the algorithm for estimating step length during hand-held walking.
• Relationship between acceleration and step direction during step direction estimation and related modification of algorithm.
• Investigation of estimation algorithms for walking methods other than hand-held walking and hand-swinging walking.
• Estimation of errors due to individual differences and adjustment of parameters (experiments with multiple participants).