353

1 INTRODUCTION

The most elementary approach to estimate the azimuth

of a moving object is to employ a magnetic compass or

a magnetometer. In the context of aircraft and ships,

where the presence of magnetic interference is a

concern, the utilization of a gyrocompass becomes

imperative [1]. Gyrocompasses are substantial in size

and challenging to maintain, thus prompting recent

interest in methodologies employing MEMS gyro

sensors to estimate azimuth. The estimation of azimuth

through the utilization of small MEMS gyro sensors is

predicated on the initial azimuth estimation, which is

ordinarily derived from due north detection. As

demonstrated by the "North Finder” [2], high-precision

MEMS sensors facilitate the detection of due north with

an accuracy of within 1° although it is a portable

device. On the other hand, small mobile vehicles, such

as UUVs, possess constrained loading capacity and are

anticipated to be utilized in substantial groups.

Consequently, it is imperative that these devices are

constructed with cost-effective components.

Consequently, the present study focuses on a method

for detecting due north using low-cost MEMS gyro

sensors. In the study by Lucian et al. [3], a low-cost

MEMS gyro sensor was mounted on a horizontal

platform. By frequently changing the sensor’s angle to

avoid the effects of random walk, the Earth’s rotation

was detected, enabling due north detection with an

accuracy of 1°. Ariffin and Arsad [4] conducted

research using sensors that were even lower in cost

than those used in Lucian et al.’s study. In this study,

the device was placed on a rotating platform, and the

signals from the accelerometer and gyro sensor were

corrected using a complementary filter. The evaluation

of the error was conducted by orienting the device in

four directions: The angles of 0°, 90°, 180°, and 270° are

to be considered. The findings indicated that due north

could be determined with an accuracy of ± 1.5° even

when employing a low-cost gyro sensor (MPU6050).

Study of North Finding System with Low-cost MEMS

Gyro Sensor and 3-axis Small Turntable under

Non-horizontal and Magnetically Disturbed

Environments for UUVs

T. Hayashi & D. Terada

National Defense Academy, Yokosuka, Kanagawa, Japan

ABSTRACT: This study developed a north-finding system using low-cost sensors and a small 3-axis turntable,

even in non-horizontal and magnetically disturbed environments. First, the accuracy of the developed turntable

was examined based on image-based analysis. Second, the turntable was rotated to estimate an ellipsoid from the

resulting point cloud, enabling use non-horizontal environment. Finally, two experiments were conducted. The

first experiment is to identify the four directions. To enhance the precision of the results, it is necessary to increase

the amount of experimental data. In order to the time required, this study adopted a method of increasing the

amount of data through simulation by sampling without replacement. The second experiment is to detect the due

north under the condition that the due north is unknown. As a result, the proposed method achieved an accuracy

of ±10° in approximately 27 min, and ±0.59° in approximately 3.29 h, demonstrating enable to use these

environments.

http://www.transnav.eu

the International Journal

on Marine Navigation

and Safety of Sea Transportation

Volume 20

Number 2

June 2026

DOI: 10.12716/1001.20.02.10

354

Ariffin and Arsad [5] further demonstrated that by

using accelerometers, gyro sensors, and GNSS signals

to correct gyro sensors, due north could be detected

from multiple locations. Note that these studies have

the following issues.

− The utilization of these devices poses a significant

challenge due to their substantial physical

dimensions, which complicates their effective

transportation.

− The existence of a horizontal platform and a

mechanism capable of altering the azimuth angle is

imperative.

− Precise determination of the northward direction is

not feasible under an inclined condition due to the

utilization of an accelerometer in the calculation of

the angle.

− The verification process in instances where the

precise direction of due north remains unknown

has not been undertaken.

− The time required to detect due north is not

discussed in any detail.

The objective of this study is to develop a portable

due north detection system that can be utilized in non-

horizontal and magnetic disturbance environments,

which can be mounted on UUVs. The system is

composed of a low-cost MEMS sensor and a small 3-

axis turntable. The small 3-axis turntable has the

capacity to rotate 360° in each axis and has a function

that automatically searches for and maintains the

horizontal position. The signal measurement and

signal processing for the low-cost MEMS sensor is

performed by Raspberry Pi. Firstly, an investigation

was conducted into the horizontal holding

functionality of the proposed system. By measuring the

acceleration while rotating the proposed system, the

tilt of each axis can be obtained. The level can be

maintained by rotating the proposed system to

opposite direction to cancel the tilt. Based on these

studies, an experiment on due north detection was

conducted. In the experiment, the orientation of the

building at the measurement point was estimated from

a map, and then the proposed system was fixed. In the

initial experiment, the proposed system was rotated in

four distinct directions. The signals obtained were then

analyzed to identify the direction of each. The objective

of this experiment is to ascertain whether the proposed

system can accurately identify the four directions when

the device is rotated in a known direction.

Additionally, the coefficients to be employed in the

calculations were determined. Firstly, the proposed

system was rotated in 16 directions, and due north

detection was performed based on the results of fitting

a cosine curve to the obtained signals. In Section 1 the

background and objectives of this paper are described.

In Section 2, an overview of the proposed equipment is

provided. In Section 3 the signal processing methods

are explained in detail. In Section 4 the experimental

results and discussions are presented. Finally, in

Section 5 the obtained findings are summarized.

2 PROPOSED SYSTEM

The proposed system consists of a low-cost MEMS

sensor (acceleration: ADXL355 and angular velocity:

MPU9250) and a small 3-axis turntable. The small 3-

axis turntable utilizes a low-cost stepper motor (28BYJ-

48), capable of rotating 360° around the x-, y-, and z-

axes, respectively. The proposed system, coordinate

system, and rotation direction are illustrated in Figure

1. As illustrated in the figure, a stepper motor is

positioned around the MEMS sensor, thereby is not

possible of utilizing a magnetic sensor due to the

constant disruption of the magnetic field by the motor

drive. And the red arrows in the figure show the

coordinate system employed in the proposed system,

which is a right-handed system with the z-axis

pointing downward. The orange arrows in the figure

indicate the positive direction of rotation ϕ, θ, and ψ

indicate rotation about the x-, y-, and z- axes,

respectively. The stepper motor (28BYJ-48) utilized in

this study possesses a step angle of 5.63°/64 (output

shaft, gear reduction ratio 1/64), a mass of

approximately 35 g, and a rotational torque of 800 gf-

cm [6] In the proposed system, the angular velocities

regarding ϕ and θ, are configured to complete a 360°

rotation within 512 steps, while the angular velocity

regarding ψ, is configured to complete a 360° rotation

within 2048 steps. This is achieved through the

utilization of gears in the device. It is evident that ϕ and

θ rotate 7.03×10

-1

° in one step and ψ rotates 1.76×10

-1

in one step. The rotation is transmitted via Pulse Width

Modulation (PWM) control from the Raspberry Pi’s

four General Purpose Input/Output (GPIO) pins. The

PWM cycles were determined as 0.01 s through a

process of trial and error. The flow of signal processing

in the proposed system is illustrated in Figure 2.

Signals αn obtained from the accelerometer are

processed using the RKF to remove outliers. The

processed signals αn

are converted to rotation angles

using Equations 4 and 5. Furthermore, the data is

converted to rotation angles in the C-frame based on

Equations 6 and 7. These angles are then converted to

angular velocities, ωn

, through numerical

differentiation expressed by Equation 9. The angular

velocity, ωn

αr

, obtained from the final accelerometer, is

the consequence of the processing of ωn

with RKF.

Signals ωn

obtained from the gyro sensor are processed

using a RKF to remove outliers. The processed signals,

ωn

, are converted to angular velocities in the C-frame

using Equation 8. For due north estimation, the

weighted average value, ωn

, of the frequency

distribution of ωn

and ωn

αr

is used. This calculation is

performed using Equations 21 to 23.

Figure 1. Overview of proposed system, coordinate system,

and rotation direction.

355

Figure 2. Process flow of proposed system.

The MPU9250 gyro sensor incorporates the same

sensor core as the MPU6050, for which previous

studies have demonstrated the capability to detect

Earth’s rotation [5]. In this study, the gyro sensor scale

was configured to ± 250°/s, with a scale factor of 131

LSB/(°/s) [7]. The accelerometer was configured to a

scale of ± 2 "g" , with a scale factor of 256000 LSB/"g"

[8].

Allan variance was employed to evaluate the error

components of the gyro sensor and accelerometer. This

facilitates the identification of sensor characteristics

such as random walk, bias instability, and drift [9]. In

particular, a slope of -1/2 in the Allan variance plot

indicates random walk, a slope of 0 corresponds to bias

instability (i.e., the resolution limit of the sensor), and

a slope of +1 reflects drift. The measurement was

conducted over a period of 2 h. Before the

measurement, the proposed system was configured

such that the rotation angles ϕn

and θn

remained

below 1 step (7.03×10

-1

°), and it was oriented toward

the East (ψ=90°). Representative Allan variance plots

for the gyro sensor and accelerometer are shown in

Figure 3 and Figure 4, respectively.

From Figure 3, it can be seen that the ωx shows a

random walk of 2.42×10

-2

°/s, bias instability of 2.22×10

°/s, and drift at 2685 s. The ωy shows a random walk

of 2.92×10

-2

°/s, bias instability of 6.14×10

-3

°/s, and drift

at 73 s. In summary, the ωx demonstrates an error of

approximately 2.42×10

-2

°/s, a resolution of

approximately 2.22×10

-3

°/s, and reliable measurements

can be obtained up to 2685 s before drift occurs. In

contrast, the ωy should be limited to measurements

within 73 s, since drift onset is 73 s. Considering Earth’s

rotation of 360 ° per day (i.e., approximately 86400 s),

the corresponding angular velocity is about 4.17×10

-3

°/s. At the measurement site in Yokosuka, Kanagawa

Prefecture (latitude 35.26°N), the local rotation rate is

approximately 3.40×10

-3

°/s. Accordingly, the proposed

system configuration allows detection of Earth’s

rotation using the gyro sensor’s ωx, but not with the ωy.

From Figure 4, it can be seen that the αx shows a

random walk of 3.53×10

-4

m/s

, bias instability of

8.35×10

-5

m/s

, and drift at 3801 s. The αy shows a

random walk of 6.25×10

-4

m/s

, bias instability of

6.51×10

-5

m/s

, and drift at 3385 s. In summary, the αx

demonstrates an error of approximately 3.53×10

-4

m/s

a resolution of approximately 8.35×10

-5

m/s

, and

reliable measurement can be obtained up to 3801 s

before drift occurs. The αy demonstrates an error of

approximately 6.25×10

-4

m/s

, a resolution of

approximately 6.51×10

-5

m/s

, and reliable

measurement can be obtained up to 3385 s before drift

occurs.

According to Equations 4 and 5, the resulting

attitude estimation errors from the proposed system

are approximately ± 3.65×10

-3

° for ϕ and ± 2.06×10

-3

for θ, both of which are considered negligible.

Figure 3. Allan variance results for the gyro sensor: The

figure shows x- and y- axes from the top.

Figure 4. Allan variance results for the accelerometer: The

figure shows x- and y- axes from the top.

2.1 Verification of accuracy of azimuth angle

The performance test was conducted to verify whether

the proposed system can continuously rotate at a

specified angular increment during operation. In this

test, control signals corresponding to 22.5° increments

regarding ψ were repeatedly sent until the device

completed 100 full rotations. The initial position and

the position after 100 full rotations were examined.

Alignment markers, indicated by red circles in Figure

5, were used to visually assess rotational accuracy. By

aligning these markers before the test, any deviation at

the end of the test could be visually confirmed. Figures

6 and 7 show the state of the proposed system before

and after the test, respectively. A comparison of the

images reveals that the alignment markers remained

nearly unchanged in position, indicating that the

proposed system achieved rotation regarding ψ

consistent with the control signals.

356

Figure 5. Overview of alignment markers.

Figure 6. State of the proposed system before test.

Figure 7. State of the proposed system after test.

2.2 Verification of accuracy of roll angle and pitch angle

The accuracy of attitude estimation regarding ϕ and θ

in the proposed system is evaluated using image

processing of a spirit level attached to the rear of the

device, following the procedure illustrated in Figure 8.

First, in the “Estimate offset and scale error process”

the offset and scale errors are calculated using the

method described in Section 3.2. Next, in the “4 angle

orientation process” starting from ψ=90°, the “Surface

estimation process” and “Capture image” steps are

performed, followed by a 90° rotation regarding ψ.

This sequence is repeated until ψ=0°. In the “Surface

estimation process” the device is rotated regarding ϕ

and θ until both θn

and ϕn

fall below one step (i.e.,

7.03×10

-1

°), based on Equations 4 and 5. In the

“Capture image” step, the device is rotated 180°

regarding θ, and a photograph is taken. After

completing all steps shown in the figure, 100

photographs are taken at each of the following

orientations: East (ψ=90°), South (ψ=180°), West

(ψ=270°), and North (ψ=0°).

Figure 8. Flow of Verification process.

An example of a captured image is shown in Figure

9. The spirit level includes circular scale markings

indicating tilts of 1°, 3°, and 5°. The image is converted

into RGB values for each pixel, and the approximate

center of the circular bubble is estimated by averaging

the maximum and minimum values along the x- and y-

axes of the blue point cloud (defined as

[R,G,B]=[0,50±10,155±20]), using the top 10 extreme

values in each direction. Using this estimated center,

both the blue and black point clouds are translated to

bring their centers closer to the origin. Next, the black

point cloud (defined as [R,G,B]=[0±50,0±50,±150]) is

used to estimate the radii and centers of the 1°, 3°, and

5° circular scale markings using the least-squares

method. The centers of these circles are iteratively

adjusted so that their average position becomes

[x,y]=[0.01,0.01] pixels. This average center is then

subtracted from both the blue and black point clouds.

Finally, the circular shape of the bubble, represented by

the blue point cloud, is estimated to use the “minimize

function” from Python’s SciPy library. An example of

the estimated scale markings and bubble shape from a

photograph taken during the attitude estimation

process is shown in Figure 10.

357

Figure 9. Captured image during ϕ and θ verification.

Figure 10. Estimated scale markings and bubble shape from

a photograph during attitude estimation.

Based on the obtained results, the attitude angles ϕi

and θi

are estimated using the following equations.

Here, ϕi

is defined as the rotation about the vertical

axis of the image, and θi

as the rotation about the

horizontal axis of the image, both following a right-

handed coordinate system.



−

(1)



−

(2)

Here, b*i

(*: x, y) denotes the center of the bubble

circle, c*i

shows the average center of the 1°, 3°, and 5°

circles, ri

is the radius of the 1° circle, and ri

is the

radius of the bubble circle. The relationship between

ϕi

, θi

and ϕi

, -θi

is summarized in Table 1.

Table 1. Relationship between ϕi

, θi

and ϕi

, -θi

(°)

ϕi

θi

ϕi

-θi

180

θi

ϕi

270

-ϕi

θi

-θi

-ϕi

The root mean square (RMS) errors of ϕ

and θ

each direction are summarized in Table 2. It can be seen

that all values fall within a range of accuracy that is

considered practically acceptable.

Table 2. RMS errors of ϕ

and θ

in each direction.

(°)

RMS ϕ

(°)

RMS θ

(°)

2.41 ×10

-2

1.14 ×10

-2

180

5.49 ×10

-2

1.82 ×10

-2

270

3.88 ×10

-3

2.09×10

-2

8.60 ×10

-3

1.40 ×10

-2

3 METHOD

3.1 Rotation of the earth

The magnetic declination between magnetic north and

due north at the measurement location in Yokosuka,

Kanagawa Prefecture, is approximately 7° [10]. Since

the gyro sensor operates in a right-handed coordinate

system with the z-axis pointing downward, it detects

an angular velocity of approximately 3.40×10

-3

°/s when

oriented toward due north. Consequently, ωxn can be

calculated using the following equation.

( )

3.40 10 cos



−

=

(3)

3.2 Ellipsoid estimation to calculate the offset error and

the scale error in accelerometer

To level the proposed system, an accelerometer is

employed. Accelerometers are subject to offset and

scale errors, which can be corrected using

measurements taken in both horizontal and vertical

rotations on the x-y, y-z, and x-z planes. In other words,

when the accelerometer is rotated once around both the

x- and y-axes, the time-series signals of the x-, y-, and

z-axes form an ellipsoid in three-dimensional space.

Theoretically, these signals should form a sphere

centered at [x,y,z]=[0,0,0] with radii equal to "g" along

the x-, y-, and z-axes. Deviations from this ideal sphere

indicate the presence of offset and scale errors. By

estimating the center and radii of the resulting

ellipsoid, the offset and scale errors of the

accelerometer can be determined. The flow of ellipsoid

estimation is demonstrated in Figure 11. Following this

procedure, the ellipsoid is estimated from the x-, y-,

and z-axes signals obtained from the accelerometer

using the “minimize function” from Python’s SciPy

library. Based on the estimated radii and center of the

ellipsoid, the offset and scale errors of the

accelerometer are calculated.

358

Figure 11. Flow of rotation and measurement using the small

turntable for ellipsoid estimation.

As the step interval increases, the calibration time

decreases; however, it may also lead to reduced

accuracy. To verify the acceptable step interval,

measurements were taken 10 times each for a total of

31 different step intervals, ranging from 4 to 124 steps.

Figures 12 and 13 show the results obtained at 4-

step and 40-step intervals, respectively. In these

figures, the blue dots represent signals from the

accelerometer, and the green sphere indicates the

estimated ellipsoid.

Figure 14 shows the standard deviations of the

centers and radii of the resulting ellipsoids. In the

figure, the blue, orange, and green dots show the

standard deviations of the x-, y-, and z-axes of the

ellipsoid centers, respectively. Similarly, the red,

purple, and brown dots show the standard deviations

of the x-, y-, and z-axes of the ellipsoid radii. From the

figure, it can be observed that when the step interval is

set to 40, all standard deviations remain below 5.00×10

m/s

, which is comparable to the error levels in the x-

and y-axes obtained from the Allan variance. Using the

ellipsoid obtained from these measurements, the

proposed system is leveled by rotating it regarding ϕ

and θ such that the corresponding angles are reduced

to less than 7.03×10

-1

° (i.e., one step regarding ϕ and θ),

based on the accelerometer signals corrected for offset

and scale errors.

Figure 12. Results measured at 4-step intervals.

Figure 13. Results measured at 40-step intervals.

Figure 14. Standard deviations of the centers and radii of the

resulting ellipsoids.

3.3 Attitude estimation from accelerometer

Attitude estimation is carried out through the

utilization of accelerometer signals and the following

equations.

sin 1 1







−

   

   

= − −  

   

   

(4)

sin 1 1







−

   

   

= −  

   

   

(5)

3.4 Quaternion-based coordinate transformation

Since the ϕ and θ rotate by 7.03×10

-1

° per step, tilts

smaller than 7.03×10

-1

° can not corrected. Therefore, it

is necessary to consider the rotation of the other axes in

the C-frame. To address this, a coordinate

transformation to the C-frame must be performed, and

the resulting data should be used to compute ωxn

. The

rotation angles in the S-frame can be expressed in

quaternion form as follows [11].

359

( )

( ) ( )

cos / 2 cos / 2

sin / 2 cos / 2

sin / 2

sin / 2 cos / 2 sin / 2

sin / 2 sin / 2

n n n





  











  





   





   



   



   



   



   



   





   

−



(6)

Where, * is quaternion product. The conversion

from quaternions to Euler angles is performed by the

following equation.

( )

 

2 2 2 2

3 2 1 0

0 1 2 3

1 3 0 2

tan

sin 2

n n n n

q q q q





−



− − +











−





(7)

The angular velocity ωn

in the C-frame can be

calculated as follows:

cos

c r c

x x n

c r c

y y n

  

  

   

   

   

(8)

3.5 Angular velocity from acceleration

The angular velocity ωxn

based on the accelerometer

signals can be calculated by using the following

equation.







−



(9)

3.6 Outlier correction using RKF [12]

To remove outliers in the acceleration αn, angular

velocity ωn, and ωn

, a RKF is used. For outlier

correction of acceleration, as the notation T means the

transpose, let the state x=[αx, αy]

, and the observation

y=[αx, αy]

. For outlier correction of angular velocity, let

the state x=[ωx, ωy]

, and the observation y=[ωx, ωy]

Their state-space model is represented by the following

Equation 10.

1n n n

n n n n

−





= + +



x x v

y x w z

(10)

where, vn is the system noise vector and is the two-

dimensional Gaussian white noise sequence N(0,Q).

wn is the observation noise vector and is the two-

dimensional Gaussian white noise sequence N(0,R).

zn is the outlier vector, R is a variance-covariance

matrix of observation noise and are represented by

Equation 12.

( )

222

x y z

Q qdiag



(11)

x xy

xy y









(12)

Note that the coefficient q is an experimentally

determined parameter: it is set to 1.0 for the

accelerometer, 1.0×10

-4

for the angular velocity

obtained from the gyro sensor, and 1.0×10

-2

for the

angular velocity estimated from the accelerometer.

The state estimation of Equation 10 can be done by

using the RKF algorithm shown in Equations 13 to 20.

[One-step-ahead prediction]

| 1 1| 1

ˆˆ

n n n n− − −

=xx

(13)

| 1 1| 1

n n n n

− − −

=+PP

(14)

[Filtering]

( )

| | 1

ˆˆ

n n n n n n

−

= + −x x e z

(15)

( )

| | 1n n n n

−

=−PP

(16)

()

n n n

−

= P σ

(17)

n n n n−

=−e y x

(18)

( )

0 otherwise

n n n n

k k k k

k n k k k k

z e e k x y



−





= +  − =







(19)

|1n n n

−

=+σP

(20)

where, I is unit matrix of 2×2.

3.7 Frequency-weighted average

The angular velocity signal, from which outliers have

been removed using the RKF, is divided into 250 bins

sorted in ascending order of



. Let mi (1≤ i ≤ 250)

denote the number of samples in the i-th bin and let



represent the corresponding angular velocity value.

Let mmax be the maximum of all mi. A conditionally

weighted average is then computed, using only those

bins for which the number of samples exceeds a

threshold defined by γ

⋅

mmax, where γ is a time constant

in the range 0<γ<1. In other words, only bins with

sufficiently high sample counts are included in the

weighted average calculation. It is important to note

that, due to the presence of significant outliers in ωx

αr

an additional condition is represented as described in

Equations 21 and 22.

max

250











(21)

max

250

0.02

250

0.02























(22)

where γ is a calibration coefficient 0.4 in this paper.

360







and





are utilized to derive





according to the following equation.

( )

1 if 0.005

otherwise

rrc

x x x

  





  











+ − 









(23)

where β is the weight coefficient 0.9 in this paper.

4 EXPERIMENTAL RESULTS AND DISCUSSIONS

4.1 Experimental environment

The experiment was conducted at the "Education and

Research Building A" of the National Defense

Academy of Japan. Figure 15 shows the building where

the experiment took place. Using the two yellow pins

in the figure, the azimuth angle was calculated to be

ψ=83.802°, based on a map and an online azimuth

calculation tool provided by the Geospatial

Information Authority of Japan [13]. To align the x-axis

of the proposed system with ψ=83.802°, the device was

rotated regarding ψ. After the rotation, the position

was recorded using an alignment markers. Next, to

orient the proposed system toward East (ψ=90°), it

needed to be rotated by 6.198°, which corresponds to

35.260 steps regarding ψ. The device was then rotated

35 steps regarding ψ, and the position was recorded

using the alignment markers.

Figure 15. The "Education and Research Building A" of the

National Defense Academy of Japan [13].

4.2 Experimental results and discussions: Identification of

the four directions

The experiment was conducted by measuring in four

directions, starting from East (ψ=90°) and rotating 90°

at each step to complete a full revolution. Figures 16

and 17 show representative time-series data from the

gyro sensor and accelerometer, respectively. In Figure

16, the time series of ωx and its filtered version ωx

are

plotted. The blue and orange lines represent ωx and ωx

respectively. The figure shows that outliers in ωx are

effectively suppressed, and the signal stabilizes

approximately 5 s. Figure 17 shows the time series of

ωx

and its filtered counterpart ωx

αr

. The blue and

orange lines represent ωx

and ωx

αr

, respectively. The

figure shows that ωx

αr

shows a smoothing effect relative

to ωx

, is effectively suppressed, and the signal

stabilizes approximately 5 s. Figures 18 and 19 show

representative frequency distributions of angular

velocity obtained from the gyro sensor and

accelerometer, respectively. In these figures, the

vertical axis indicates the frequency of occurrence, and

the horizontal axis shows angular velocity. The blue

line, the yellow line, the orange line, and the green line

respectively show the observed distribution, the

threshold used in Equations 21 and 22, the weighted

average calculated using Equations 21 and 22, and the

mean of the entire distribution. These figures highlight

that the frequency-weighted average differs from the

simple mean. Empirically, the frequency-weighted

average was found to yield higher estimation accuracy

and was therefore adopted in subsequent analysis.

Figure 16. Time series of ωx and ωx

:The figure shows

direction ψ=90°, 180°, 270°, and 0° from the top.

361

Figure 17. Time series of ωx

and ωx

αr

:The figure shows

direction ψ=90°, 180°, 270°, and 0° from the top.

Figure 18. Frequency distribution of the ωx

: The figure shows

direction ψ=90°, 180°, 270°, and 0° from the top.

362

Figure 19. Frequency distribution of the ωx

αr

:The figure

shows direction ψ=90°, 180°, 270°, and 0° from the top.

If the signal obtained when the device is oriented

toward North (ψ=0°) shows the maximum value, or the

signal at South (ψ=180°) shows the minimum value, it

can be concluded that the four cardinal directions have

been correctly identified. A total of 214 experiments

were conducted. For each direction, the signal

measured over a 25 s interval was analyzed using the

values calculated from Equations 21 and 22. The

average values of ωx

, ωx

αr

, and ωx

across all 214 trials

were computed using data from these 5 to 10 s

intervals. The results are presented in Figure 20. In the

figure, the blue dots, the purple dots, the red dots, the

green dots, and the yellow dots respectively show

average values of ωx

, ωx

αr

, ωx

, the theoretical values,

and bias error range calculated based on Allan

variance, centered around the theoretical values. These

results confirm the reliable identification of the four

cardinal directions from the average outcomes of 214

trials.

Figure 20. Result of average for the four directions.

While the average across 214 trials indicates that the

four cardinal directions can be identified, the

probability of correctly identifying all four directions

in a single trial (i.e., one full rotation) was 58.4 %. This

implies that a single rotation does not always

guarantee successful identification of all directions. To

improve accuracy, combinations of the 214 trials were

used to simulate multiple rotations, and the probability

of successful direction identification was evaluated.

For 1 to 6 rotations, all possible combinations were

considered. However, exhaustively evaluating all

combinations beyond this range would impose a

significant computational burden. Therefore, for three

or more rotations, 100,000 combinations were obtained

by sampling without replacement to estimate the

identification probability. The combination was

executed by utilizing a calculation program that was

developed using Python. The results of these

simulations are shown in Figure 21. In the figure, red

circles indicate the probabilities calculated from all

combinations (for 1 to 6 rotations), while green dots

represent the probabilities estimated from the 100,000

randomly sampled combinations (for 3 to 10 rotations).

The close agreement between the red circles and green

dots confirms the reliability of the sampling-based

estimates. According to the figure, the probability of

correctly identifying the four cardinal directions

reaches 87.7 % with 10 simulated rotations. One full

rotation of the device takes 20.48 s, and an additional

10 s of measurement is required for each direction.

Therefore, identifying the directions takes about 1 min

per trial, resulting in a 58.4 % success rate for a single

trial and 87.7 % after approximately 10 min (i.e., 10

rotations).

Figure 21. Results of simulated multiple rotations and their

probabilities for identification of the four directions.

363

4.3 Experimental results and discussions: Detection of

due north

The experiment was conducted by rotating the

proposed system in increments of 22.5° regarding ψ,

orienting it in 16 directions to complete one full

rotation. The rotation was set to East (ψ=90°). The

obtained results were divided into two halves, and the

phase a in Equation 24 was estimated using Python’s

curve fitting function. The average of the estimated

phase values a shows the estimated initial angle at the

start of the measurement. Given the confirmed

accuracy of z-axis rotation, due north can be estimated

from the phase shift in the fitted cosine curve.

( )

cos

A a b



= + +

(24)

Note that A denotes the amplitude, which is set to

3.40×10

-3

° based on Equation 3, and b is set to zero

because the mean value has been subtracted from the

data. Figure 22 shows the average of 77 trials in which

the proposed system was rotated in 16 directions at

22.5° intervals, starting from ψ=90°. In the figure, the

blue dots, the purple dots, the red dots, the yellow

lines, the light blue dots, and the green dots

respectively show average values of ωx

, ωx

αr

, ωx

, the

estimated cosine wave, the estimated due north, and

the theoretical values. The results demonstrate that due

north was successfully estimated.

Figure 22. Average angular velocity results from 77 trials

with 16 orientations, starting at ψ=90°.

To simulate a scenario in which due north was not

known in advance, the measurement was started from

a position rotated 11.25° to the right of East (ψ=101.25°).

In this case, due north falls between the measurement

directions. Figure 23 shows the average of 77 trials in

which the proposed system was rotated in 16

directions at 22.5° intervals, starting from ψ=101.25°. In

the figure, the blue dots, the purple dots, the red dots,

the yellow lines, the light blue dots, and the green dots

respectively show average values of ωx

, ωx

αr

, ωx

, the

estimated cosine wave, the estimated due north, and

the theoretical values. The results demonstrate that due

north was successfully estimated, even when it was not

known beforehand. As in the four direction case,

simulated rotations were performed, and cosine waves

were fitted to the average values at each rotation.

Figure 24 shows the RMS error of due north

estimation via simulated multiple rotations. In the

figure, the RMS error was within ±10° for 10 rotations,

and within ±0.59° for 77 rotations. This corresponds to

an accuracy of ±10° in approximately 27 min, and

±0.59° in approximately 3.29 h. These results confirm

that the proposed method enables accurate estimation

of due north, even when its direction is not known in

advance. This performance is primarily attributed to

the following three factors:

− Dividing the signal into two halves for cosine fitting

− Estimating the phase component of the fitted curve

− Applying the frequency-weighted averaging

method using Equations 21 and 22

The achieved accuracy of ±0.59° surpasses the

previously reported value of less than ±1° using the

same sensor, demonstrating the effectiveness of the

proposed approach. Future work includes evaluating

the method’s performance without applying

corrections for rotations regarding ϕ and θ on the 3-

axis turntable, as well as improving computational

efficiency.

Figure 23. Average angular velocity results from 77 trials

with 16 orientations, starting at ψ=101.25°.

Figure 24. RMS error of due north estimation via simulated

multiple rotations.

5 CONCLUSIONS AND FUTURE WORKS

In this study, we developed a portable due north

detection system capable of operating under non-

horizontal and magnetically disturbed environments.

The system comprises a small 3-axis turntable, with

each axis capable of 360° rotation, and a device

equipped with a low-cost MEMS sensor mounted on a

Raspberry Pi. To evaluate the system, two types of

experiments were conducted. The key findings are

summarized as follows:

1. Experiment of identification of the four directions

Weight coefficients in Equation 23 were

determined, and multiple simulated rotations were

conducted for validation. As a result, the proposed

system was able to identify the four cardinal

directions with a probability of 58.4 % in

approximately 1 min, and 87.7 % in approximately

10 min.

2. Experiment of detection of due north

364

The device was rotated in increments of 22.5° to

complete a full circle. Two scenarios were tested:

one starting from ψ=90°, and another starting from

a position rotated 11.25° to the right of ψ=90°,

simulating a scenario in which due north was not

known in advance. The acquired data were divided

into two halves, and the phase parameter a in

Equation 24 was estimated using Python’s curve

fitting function. As a result, an accuracy of within

±10° was achieved in approximately 27 min, and

within ±0.59° in approximately 3.29 h. These results

confirm that the proposed method enables accurate

estimation of due north, even when its direction is

initially unknown.

This performance is attributed to the

implementation of the following three techniques:

1. Dividing the signal into two halves for cosine fitting

2. Estimating the phase component of the fitted curve

3. Applying the frequency-weighted averaging

method using Equations 21 and 22

The achieved accuracy of ±0.59° surpasses the

previously reported value of less than ±1° using the

same sensor [5], demonstrating the effectiveness of the

proposed approach.

Future work includes the following:

1. Since the data are already transformed into the

C-frame, corrections for the turntable’s ϕ and θ

angles are theoretically unnecessary. It is necessary

to verify whether due north can still be accurately

estimated under this condition.

2. Improving detection speed, for example by

increasing the number of sensors or optimizing the

processing algorithm.

NOMENCLATURE

S-frame

Sensor frame: The right-hand side system with z- axis

pointing downward.

C-frame

Computer frame: Rotation relative coordinate system.

RKF

Robust Kalman filter

αn

Acceleration measured in S-frame (m/s

)

αn

Acceleration filtered by RKF in S-frame (m/s

)

ωn

Angular velocity measured in S-frame (°/s)

ωn

Angular velocity filtered by RKF in S-frame (°/s)

ωn

Angular velocity in C-frame (°/s)

ωn

Angular velocity calculated from αn in C-frame (°/s)

ωn

αr

Angular velocity calculated from αn filtered by RKF in C-

frame (°/s)

ωψ

Combination angular velocity in C-frame (°/s)

Roll angle around the x- axis (°)

Pitch angle around y- axis (°)

Azimuth angle (°)

ϕn

Roll angle calculated from αn (°)

ϕn

Roll angle in C-frame (°)

θn

Pitch angle calculated from αn (°)

θn

Pitch angle in C-frame (°)

Δt

Sampling interval ("s")

"g"

Gravitational acceleration (m/s

)

ACKNOWLEDGMENTS

We would like to express our sincere gratitude to Professor

Takita for his valuable insights regarding the selection of the

sensor used in this study.

REFERENCES

[1] T. Mozai and M. Kobayashi, Theory and Practice of

Compass and Gyroscope. Kaibundo Publishing, 1971, pp.

69–73, (in Japanese).

[2] Sumitomo Precision Products Co., Ltd. “Northfinder

attitude & heading reference systems (AHRS) GCAH-

12C-02 (premium model) datasheet.” (2025), [Online].

Available: https://www.spp.co.

jp/mems/__assets/mems/en/product/prod04/GCAH-12C-

02_Brochure_20230913.pdf.

[3] L. I. Iozan, M. Kirkko-Jaakkola, J. Collin, J. Takala, and C.

Rusu, “North finding system using a mems gyroscope,”

in Proceedings of European Navigation Conference on

Global Navigation Satellite Systems, 2010.

[4] N. H. Ariffin and N. Arsad, “Mems gyroscope and

accelerometer based north finding system,” ARPN

Journal of Engineering and Applied Sciences, vol. 12, no.

22, pp. 6517–6522, 2017.

[5] N. H. Ariffin and N. Arsad, “Mems gyro and

accelerometer as north-finding system for bulk direction

marking,” Ieee Access, vol. 10, pp. 114 214–114 222, 2022.

[6] AKIZUKI DENSHI TSUSHO CO.,LTD. “28BYJ-48

datasheet.” (2019), [Online]. Available: https://

akizukidenshi.com/goodsaffix/28BYJ48-W01.pdf.

[7] TDK Corporation. “MPU9250 Product Specification

Revision 1.1.” (2016), [Online]. Available: https:

//invensense.tdk.com/wp-content/uploads/2015/02/PS-

MPU-9250A-01-v1.1.pdf.

[8] Analog Devices. “ADXL354/ADXL355: Low Noise, Low

Drift, Low Power, 3-Axis MEMS Accelerometers Data

Sheet (Rev.C).” (2024), [Online]. Available:

https://www.analog.com/media/en/ technical-

documentation/data-sheets/adxl354_adxl355.pdf.

[9] T. Imamura, “Introduction of GNSS/INS,” Journal of the

Robotics Society of Japan, vol. 37, no. 7, pp. 579–584, 2019,

(in Japanese).

[10] Geospatial Information Authority of Japan.

“Geomagnetic field values at epoch 2020.0 (as of 00:00 UT

on January 1, 2020).” (2022), [Online]. Available: https : /

/ www . gsi . go . jp / common / 000236992.pdf.

[11] Y. Shimada and S. Sasa, Flight Dynamics, expanded

edition. Morikita Publishing Co., Ltd., 2017, p. 271, (in

Japanese).

[12] Y. Kaneda and Y. Irizuki, “Development of an outlier

reduction filter,” Tokyo Metropolitan Industrial

Technology Research Center Research Report, no. 8, pp.

2–5, 2013, (in Japanese).

[13] Geospatial Information Authority of Japan. “Map of

geospatial information authority of japan.” (2023),

[Online]. Available: https://maps.gsi.go.jp/.