Forem: William Pederzoli

The Navigation Paradox: Why Perfect Sensors Can't Find Their Way Home

William Pederzoli — Tue, 11 Nov 2025 10:24:58 +0000

Introduction

There's a fascinating paradox in autonomous navigation: the most precise motion sensors become useless for positioning within minutes. This isn't a failure of technology, but a fundamental law of physics that every navigation engineer must confront.

I've been exploring this paradox through IMU mechanization. The art of transforming raw sensor data into meaningful motion estimates. What I've learned challenges conventional wisdom about what's possible with inertial navigation.

The Drift Dilemma: When Precision Meets Uncertainty

Imagine having a sensor that can detect movements finer than a human hair's width, yet can't tell you where you are after a short walk. This is the reality of inertial measurement units (IMUs). The reason being Error accumulation.

Every microscopic imperfection in sensor readings, whether from temperature variations, manufacturing tolerances, or electrical noise, compounds over time. What starts as nanometer-scale errors grows into meter-scale uncertainties in minutes. It's mathematical certainty, not engineering failure.

Beyond Simple Math: The Complexity of Motion Tracking

Mechanization isn't simply integrating acceleration to get position. The process involves:

Orientation estimation that accounts for both vehicle rotation and Earth's spin
Gravity compensation that separates actual motion from gravitational effects
Coordinate transformations between body-fixed and global reference frames
Error propagation modeling that anticipates how uncertainties evolve

Each layer introduces complexity that can't be solved with textbook equations alone.

The Strategic Insight: Embracing Sensor Limitations

Our biggest breakthrough came from recognizing that fighting drift is a losing battle. Instead, we've learned to work within sensors' inherent limitations:

Short-term excellence: IMUs provide unparalleled motion tracking between GPS updates
Complementary strengths: They excel where other sensors fail—in tunnels, urban canyons, and signal-denied environments
Fusion-friendly design: By understanding IMU behavior intimately, we can design better fusion algorithms

The Visualization Advantage: Seeing the Invisible

One of our most powerful tools has been advanced visualization. By creating intuitive representations of complex error patterns, we can:

Identify characteristic "signatures" of different error sources
Validate system performance against physical ground truth
Communicate complex concepts to non-technical stakeholders
Accelerate debugging and optimization cycles

Looking Forward: The Fusion Frontier

Mechanization isn't the end goal but the essential foundation for what comes next. By thoroughly understanding how inertial sensors behave in isolation, we're building the expertise needed to create robust multi-sensor systems that navigate reliably in challenging environments.

The real magic happens when we combine complementary technologies, creating navigation solutions that are greater than the sum of their parts.

To fellow innovators in autonomous systems: How do you balance the trade-offs between different sensing modalities? What strategies have you found for working within fundamental physical limitations?

Building an IMU Simulator from the Ground Up: A Journey Through Inertial Navigation

William Pederzoli — Wed, 05 Nov 2025 16:07:55 +0000

Introduction

In my previous articles, I've built a complete GNSS positioning engine. But for robust autonomous navigation, we need more than satellites... We need an Inertial Measurement Unit (IMU) that works when GNSS signals drop. However, testing fusion algorithms with real hardware is expensive and slow. So the best solution is to build a high-fidelity IMU simulator.

This article is the process of how to implement an IMU simulator from first principles, breaking down the problem into four logical levels of increasing complexity.

🎯 LEVEL 1: IMU Fundamentals - More Than Just Raw Data

An IMU doesn't measure what you might intuitively think. Understanding this distinction is crucial for anyone working with inertial systems.

1.1 What's Really in an IMU?

At its core, an IMU contains two fundamental sensors:

Accelerometers: Measure specific force (all reaction forces except gravity)
Gyroscopes: Measure angular rate (how fast the body is rotating in 3D space)

Think of it this way: if you're in a elevator accelerating upward, the accelerometer measures both the elevator's acceleration AND the gravity you're pushing against.

1.2 The Critical Insight: Specific Force ≠ Acceleration

This is one of the most fundamental concepts:

What accelerometers measure = (All real forces except gravity) / mass

Or mathematically:

s p ec i f i c f orce = t r u e a cce l er a t i o n - g r a v i t y v ec t or

This matter tremendously because when our rover is perfectly stationary on Earth, the accelerometers don't read zero, they measure approximately 9.8 m/s² upward (the normal force counteracting gravity). This is why we must "remove gravity" during inertial navigation.

1.3 Why Simulate Errors? The Realism vs. Theory Gap

Real IMUs are far from perfect. Simulation lets us bridge this gap:

Test algorithms under controlled, repeatable error conditions
Understand how different error types accumulate over time
Develop robust filters before investing in expensive hardware
Isolate specific error effects that would be impossible to separate in real data

🎯 LEVEL 2: Calculating the "Ground Truth" IMU Data

Before we can add realistic errors, we need perfect measurements derived from the rover's known trajectory. This reverse-engineering is surprisingly complex.

2.1 Rotation Matrices: The Language of 3D Orientation

While Euler angles (roll, pitch, yaw) are intuitive for humans, rotation matrices are computationally superior for navigation. Euler angles can suffer from "gimbal lock" (losing a degree of freedom when two axes align), while rotation matrices avoid this entirely.

A rotation matrix describes how to transform vectors from one coordinate frame to another. For the aerospace sequence (yaw → pitch → roll), we build the rotation from three elemental rotations:

The rotation from body to navigation frame (NED) is constructed as:

R_b^{n} = R_x(\phi) \cdot R_y(\theta) \cdot R_z(\psi)

Where the elemental rotation matrices are:

R_z(\psi) = \begin{bmatrix} \cos\psi & -\sin\psi & 0 \\ \sin\psi & \cos\psi & 0 \\ 0 & 0 & 1 \end{bmatrix}

R_y(\theta) = \begin{bmatrix} \cos\theta & 0 & \sin\theta \\ 0 & 1 & 0 \\ -\sin\theta & 0 & \cos\theta \end{bmatrix}

R_x(\phi) = \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos\phi & -\sin\phi \\ 0 & \sin\phi & \cos\phi \end{bmatrix}

When expanded, the complete rotation matrix becomes:

R_b^{n} = \begin{bmatrix} c_\theta c_\psi & c_\theta s_\psi & -s_\theta \\ s_\phi s_\theta c_\psi - c_\phi s_\psi & s_\phi s_\theta s_\psi + c_\phi c_\psi & s_\phi c_\theta \\ c_\phi s_\theta c_\psi + s_\phi s_\psi & c_\phi s_\theta s_\psi - s_\phi c_\psi & c_\phi c_\theta \end{bmatrix}

where
c=cos, s=sin.

Geometric Interpretation:

Column 1: Direction of body X-axis in NED coordinates
Column 2: Direction of body Y-axis in NED coordinates
Column 3: Direction of body Z-axis in NED coordinates

Key Properties:

The matrix is orthonormal: $R^-1 = R^T$
Determinant is +1 (preserves orientation)
Composition of rotations is simple matrix multiplication

Pseudo-code Implementation

FUNCTION create_rotation_matrix(roll_angle, pitch_angle, yaw_angle)

    // Precompute trigonometric functions
    cos_roll = COS(roll_angle)
    sin_roll = SIN(roll_angle)
    cos_pitch = COS(pitch_angle) 
    sin_pitch = SIN(pitch_angle)
    cos_yaw = COS(yaw_angle)
    sin_yaw = SIN(yaw_angle)

    // Build the combined rotation matrix (roll × pitch × yaw)
    rotation_matrix = [
        // Row 1
        [cos_pitch * cos_yaw, 
         cos_pitch * sin_yaw, 
         -sin_pitch],

        // Row 2  
        [sin_roll * sin_pitch * cos_yaw - cos_roll * sin_yaw,
         sin_roll * sin_pitch * sin_yaw + cos_roll * cos_yaw, 
         sin_roll * cos_pitch],

        // Row 3
        [cos_roll * sin_pitch * cos_yaw + sin_roll * sin_yaw,
         cos_roll * sin_pitch * sin_yaw - sin_roll * cos_yaw,
         cos_roll * cos_pitch]
    ]

    RETURN rotation_matrix

END FUNCTION

// Usage example:
body_to_ned_rotation = create_rotation_matrix(roll_rad, pitch_rad, yaw_rad)

// Transform a vector from body to NED frame:
vector_ned = MATRIX_MULTIPLY(body_to_ned_rotation, vector_body)

2.2 The Complete Rotation Model: It's Not Just About the Vehicle

The orientation change between time steps comes from two sources:

Body rotation: The rover's own movement (turning, tilting)
Earth rotation: The Earth spinning under us at ~15 degrees per hour

// Both rotations must be combined
body_rotation = calculate_from_angular_velocity(gyro_measurements)
earth_rotation = calculate_earth_rotation_at_current_latitude()

// The total rotation is the combination of both
total_rotation = combine_rotations(body_rotation, earth_rotation)

2.3 Rodrigues' Rotation Formula: Efficient 3D Updates

For discrete time steps, we need an efficient way to update rotation matrices. Rodrigues' rotation formula gives us this:

\sin(\theta)K + (1 - \cos(\theta))K^2

Where:

\begin{bmatrix} 0 & -k_z & k_y \ k_z & 0 & -k_x \ -k_y & k_x & 0 \end{bmatrix}

// Given an angular velocity vector [wx, wy, wz] and time step dt
rotation_vector = angular_velocity × time_step
rotation_matrix = rodrigues_formula(rotation_vector)

// This efficiently computes the rotation matrix
// It's really useful for real-time applications

2.4 The Specific Force Transformation

Now we can compute what ideal accelerometers would measure:

// Step 1: Get the true acceleration from position changes
true_acceleration_ecef = calculate_acceleration_from_trajectory()

// Step 2: Get gravity at current position  
gravity_vector_ecef = calculate_gravity_model(current_latitude, current_height)

// Step 3: Compute specific force in ECEF frame
specific_force_ecef = true_acceleration_ecef - gravity_vector_ecef

// Step 4: Transform to body frame (what the physical sensors would measure)
specific_force_body = transpose(rotation_matrix) × specific_force_ecef

🎯 LEVEL 3: Simulating Realistic IMU Errors

Now we transform our ideal measurements into the noisy, imperfect readings that real sensors produce. This is where simulation meets reality.

3.1 The Zoo of IMU Errors

Real IMUs suffer from multiple error sources that interact in complex ways:

Biases: Constant offsets that slowly drift over time
- Example: Always reading 0.1 m/s² even when stationary
Scale Factors: Gain errors in the measurement chain
- Example: 1% error means 1 m/s² becomes 1.01 m/s²
Cross-Coupling: Axis misalignment and interference
- Example: X-axis acceleration affecting Y-axis measurements
Noise: Random high-frequency variations
- Example: Small random fluctuations around the true value
Quantization: Digital resolution limits
- Example: Can only measure in 0.01 m/s² increments

3.2 Mathematical Error Model

Each error type has a specific mathematical representation:

// The complete error model for one axis
measured_value = (
    (true_value + constant_bias + time_varying_bias) 
    × (1 + scale_factor_error)
    + cross_coupling_from_other_axes
    + random_noise
    + quantization_error
)

3.3 Implementing Realistic Error Simulation

Here's how we apply errors in practice:

FUNCTION simulate_imu_errors(ideal_measurements, error_parameters)

    // Initialize output measurements
    realistic_measurements = copy(ideal_measurements)

    FOR each time_step IN realistic_measurements:

        // Apply biases (constant + slow drift)
        realistic_measurements.gyro_x[time_step] = (
            ideal_measurements.gyro_x[time_step] 
            + error_parameters.gyro_bias_x
            + calculate_bias_drift(time_step)
        )

        // Apply scale factors
        realistic_measurements.gyro_x[time_step] = (
            realistic_measurements.gyro_x[time_step]
            × (1 + error_parameters.gyro_scale_x)
        )

        // Apply cross-coupling between axes
        realistic_measurements.gyro_x[time_step] = (
            realistic_measurements.gyro_x[time_step]
            + error_parameters.cross_coupling_xy × ideal_measurements.gyro_y[time_step]
            + error_parameters.cross_coupling_xz × ideal_measurements.gyro_z[time_step]
        )

        // Add random noise
        realistic_measurements.gyro_x[time_step] = (
            realistic_measurements.gyro_x[time_step]
            + generate_gaussian_noise(error_parameters.noise_std)
        )

        // Apply quantization based on sensor resolution
        realistic_measurements.gyro_x[time_step] = (
            round(realistic_measurements.gyro_x[time_step] / resolution) × resolution
        )

    END FOR

    RETURN realistic_measurements

3.4 The Complete Measurement Pipeline

Putting it all together:

// Input: Perfect trajectory data
// Output: Realistic IMU measurements

ideal_angular_rates = calculate_ideal_gyro_measurements(true_trajectory)
ideal_specific_forces = calculate_ideal_accel_measurements(true_trajectory)

realistic_angular_rates = simulate_imu_errors(ideal_angular_rates, gyro_error_params)
realistic_specific_forces = simulate_imu_errors(ideal_specific_forces, accel_error_params)

The Complete Simulation Pipeline

Putting all three levels together gives us a powerful IMU data generator:

// Input: Reference trajectory (ground truth)
// Output: Realistic IMU measurements

true_trajectory = load_reference_data()

// Level 2: Work backward from trajectory to ideal measurements
ideal_imu = calculate_ideal_imu_measurements(true_trajectory)

// Level 3: Add realistic errors
realistic_imu = simulate_imu_errors(ideal_imu, imu_error_parameters)

// Output: Ready-to-use IMU data for testing fusion algorithms
RETURN realistic_imu

Validation and Critical Insights

The validation process revealed several key insights:

Visual Analysis: By plotting the simulated IMU data against the ground truth, I could immediately spot if the error models were behaving realistically. The noise should be random, biases should be constant, and scale factors should proportionally affect all measurements.
Error Correlation: Different error types create distinctive patterns in the data. White noise appears as high-frequency jitter, while bias drift creates low-frequency wander in the measurements.
The Gravity Test: A quick sanity check—when the rover is stationary and level, the accelerometers in the body frame should measure the gravity vector. If the rover is perfectly aligned with the NED frame, this would be [0, 0, 9.8] m/s², but in practice, the measured specific force depends on the rover's orientation relative to the local vertical.

Conclusion: Foundation for Sensor Fusion

Building this IMU simulator has been an incredible journey through the physics of inertial navigation. From understanding the fundamental concept of specific force to implementing sophisticated error models, each level built upon the previous one.

The most valuable lesson? Simulation forces you to understand the system at a deeper level. You can't simulate what you don't understand.

This IMU simulator now provides the perfect foundation for the next exciting phase: IMU Mechanization. In the next article, I'll show how we can take these simulated measurements and reconstruct the rover's trajectory through numerical integration and discover why IMUs alone can never provide long-term accurate navigation.

But that's a story for next time...

For those working with inertial systems: What aspects of IMU error modeling have you found most challenging? Are there particular error sources that surprised you with their impact on navigation performance?

If you've made it this long, thank you and hope this is useful content for your endeavors. See you in the next one!

From Measurements to Map: A Deep Dive into the Least Squares Algorithm for GNSS

William Pederzoli — Mon, 27 Oct 2025 20:21:48 +0000

Introduction

In my previous article, I explained how we generate corrected GNSS measurements. Today we take the crucial step: transforming those measurements into estimated position and velocity using the Least Squares algorithm. This isn't just any implementation, it's the mathematical heart of any satellite navigation system.

The Fundamental Problem: Too Many Unknowns

Each pseudorange measurement gives us a nonlinear equation:

ρᵢ = ||p_satᵢ - p_rcv|| + c·δt + εᵢ

Where:

ρᵢ = Measured pseudorange for satellite i

p_satᵢ = Known satellite position

p_rcv = Unknown receiver position (x, y, z)

c = Speed of light

δt = Receiver clock bias

εᵢ = Residual errors

With 4 satellites, we have 4 equations and 4 unknowns. But in practice, we use more satellites, creating an overdetermined system that we solve with** Least Squares**.

Linearization :The Key to Making It Solvable

The problem is that our equations are nonlinear (due to the ||·|| norm). The solution is to linearize around an initial estimate:

Δρ = G·Δx + ε

Where:

Δρ = Vector of residuals (measured - calculated)

G = Geometry matrix (Jacobian)

Δx = Correction to our state estimate

ε = Vector of errors

Detailed Anatomy of the LSQ Algorithm

1. System State
Our state vector has 8 components:

state = [
    position_x,           // ECEF X coordinate [meters]
    position_y,           // ECEF Y coordinate [meters] 
    position_z,           // ECEF Z coordinate [meters]
    clock_bias,           // Receiver clock bias converted to meters [m]
    velocity_x,           // ECEF X velocity [m/s]
    velocity_y,           // ECEF Y velocity [m/s]
    velocity_z,           // ECEF Z velocity [m/s]
    clock_drift           // Receiver clock drift [m/s]
]

2. The Iterative Process

// Initialize state (often with zeros or approximate position)
current_position_x = 0.0
current_position_y = 0.0  
current_position_z = 0.0
current_clock_bias = 0.0
current_velocity_x = 0.0
current_velocity_y = 0.0
current_velocity_z = 0.0
current_clock_drift = 0.0

FOR iteration = 1 TO maximum_iterations DO

    // Step 1: Calculate predictions and residuals for all satellites
    pseudorange_residuals = []
    range_rate_residuals = []
    geometry_matrix_pseudorange_rows = []
    geometry_matrix_range_rate_rows = []

    FOR EACH satellite IN visible_satellites DO

        // Calculate satellite-to-receiver geometry
        delta_x = satellite_position_x - current_position_x
        delta_y = satellite_position_y - current_position_y  
        delta_z = satellite_position_z - current_position_z

        geometric_distance = SQUARE_ROOT(delta_x² + delta_y² + delta_z²)

        // Unit vector pointing from receiver to satellite
        unit_vector_x = delta_x / geometric_distance
        unit_vector_y = delta_y / geometric_distance
        unit_vector_z = delta_z / geometric_distance

        // Predicted pseudorange = geometric distance + clock bias
        predicted_pseudorange = geometric_distance + current_clock_bias

        // Residual = measured - predicted
        pseudorange_residual = measured_pseudorange - predicted_pseudorange

        // Geometry matrix row for pseudorange measurement
        pseudorange_geometry_row = [
            -unit_vector_x, -unit_vector_y, -unit_vector_z,  // Position partials
            1.0,                                            // Clock bias partial  
            0.0, 0.0, 0.0, 0.0                             // Velocity and clock drift don't affect pseudorange
        ]

        // For range rate (Doppler) measurements
        delta_velocity_x = satellite_velocity_x - current_velocity_x
        delta_velocity_y = satellite_velocity_y - current_velocity_y
        delta_velocity_z = satellite_velocity_z - current_velocity_z

        // Predicted range rate = relative velocity projection + clock drift
        predicted_range_rate = (delta_velocity_x * unit_vector_x + 
                              delta_velocity_y * unit_vector_y + 
                              delta_velocity_z * unit_vector_z) + current_clock_drift

        range_rate_residual = measured_range_rate - predicted_range_rate

        // Geometry matrix row for range rate measurement  
        range_rate_geometry_row = [
            0.0, 0.0, 0.0, 0.0,                            // Position and clock bias don't affect range rate
            -unit_vector_x, -unit_vector_y, -unit_vector_z, // Velocity partials
            1.0                                            // Clock drift partial
        ]

        // Collect all residuals and geometry rows
        ADD pseudorange_residual TO pseudorange_residuals
        ADD range_rate_residual TO range_rate_residuals
        ADD pseudorange_geometry_row TO geometry_matrix_pseudorange_rows
        ADD range_rate_geometry_row TO geometry_matrix_range_rate_rows
    END FOR

    // Step 2: Build complete least squares system
    all_residuals = COMBINE(pseudorange_residuals, range_rate_residuals)
    full_geometry_matrix = STACK(geometry_matrix_pseudorange_rows, geometry_matrix_range_rate_rows)

    // Step 3: Solve for state corrections
    information_matrix = TRANSPOSE(full_geometry_matrix) × full_geometry_matrix
    covariance_matrix = INVERSE(information_matrix)
    correlation_vector = TRANSPOSE(full_geometry_matrix) × all_residuals
    state_corrections = covariance_matrix × correlation_vector

    // Step 4: Update state estimate
    current_position_x = current_position_x + state_corrections[0]
    current_position_y = current_position_y + state_corrections[1] 
    current_position_z = current_position_z + state_corrections[2]
    current_clock_bias = current_clock_bias + state_corrections[3]
    current_velocity_x = current_velocity_x + state_corrections[4]
    current_velocity_y = current_velocity_y + state_corrections[5]
    current_velocity_z = current_velocity_z + state_corrections[6]
    current_clock_drift = current_clock_drift + state_corrections[7]

    // Step 5: Check convergence
    correction_magnitude = NORM(state_corrections)
    IF correction_magnitude < convergence_tolerance THEN
        BREAK  // Solution has converged
    END IF
END FOR

Mathematical Explanation

Why This Formula Works
We minimize the cost function:

J(Δx) = ‖Δρ - G·Δx‖²

The optimal solution occurs when the gradient is zero:

∇J(Δx) = -2Gᵀ(Δρ - G·Δx) = 0

This leads to the normal equations:

GᵀG·Δx = GᵀΔρ

And finally:

Δx = (GᵀG)⁻¹GᵀΔρ

Geometric Interpretation:

G projects state corrections from state space to measurement space
(GᵀG)⁻¹Gᵀ is the pseudo-inverse that finds the minimum-norm solution
Each iteration brings us closer to the point where residuals are orthogonal to G's column space

Quality Metrics: Dilution of Precision (DOP)

From the covariance matrix we extract essential metrics:

covariance_matrix = (GᵀG)⁻¹
GDOP = SQRT(SUM(diagonal(covariance_matrix)))
PDOP = SQRT(SUM(diagonal(covariance_matrix[0:3,0:3])))
TDOP = SQRT(covariance_matrix[3,3])

DOP Interpretation:

GDOP < 3: Excellent geometry
GDOP 3-6: Good geometry
GDOP > 6: Poor geometry, degraded accuracy

Validation in Local NED Frame

We transform to geodetic coordinates and calculate errors:

estimated_latitude, estimated_longitude, estimated_altitude = CONVERT_ECEF_TO_GEODETIC(
    current_position_x, current_position_y, current_position_z
)

estimated_velocity_north, estimated_velocity_east, estimated_velocity_down = 
    CONVERT_VELOCITY_ECEF_TO_NED(
        current_velocity_x, current_velocity_y, current_velocity_z,
        estimated_latitude, estimated_longitude
    )

// Calculate errors in NED frame
radius_north, radius_east = CALCULATE_EARTH_RADII_OF_CURVATURE(true_latitude)

north_position_error = (estimated_latitude - true_latitude) × (radius_north + true_altitude)
east_position_error = (estimated_longitude - true_longitude) × (radius_east + true_altitude) × COS(true_latitude)
down_position_error = -(estimated_altitude - true_altitude)

horizontal_position_error = SQUARE_ROOT(north_position_error² + east_position_error²)
vertical_position_error = ABSOLUTE_VALUE(down_position_error)

Key Lessons Learned

Initialization: Starting near the true solution significantly speeds up convergence.
Satellite Geometry: The distribution of satellites in the sky (not the number) determines final accuracy.
Numerical Conditioning: When GᵀG is nearly singular (aligned satellites), small measurement errors produce large position errors.
Convergence: Linearization is only valid near the true solution, hence the need for iterations.

Conclusion: Where Theory Meets Practice

Implementing Least Squares from scratch reveals the elegant connection between geometry, linear algebra, and statistical estimation. Every component of the algorithm has deep physical meaning:

Unit vectors capture directional sensitivity
G matrix encodes satellite-receiver geometry
Residuals measure model consistency
Covariance matrix quantifies inherent uncertainty

This deep understanding is essential before advancing to more sophisticated techniques like the Kalman Filter for sensor fusion.

For those who have worked in navigation: Which theoretical aspects of LSQ did you find most fundamental for robust implementation? How do you handle cases with poor satellite geometry?

From Data to Physics: Building the GNSS Measurement Engine

William Pederzoli — Mon, 20 Oct 2025 19:04:26 +0000

Introduction

In my last article, I discussed the data files that power my sensor fusion project. But raw data alone isn't enough. The real challenge begins when you must transform those satellite ephemerides and reference trajectories into physically meaningful measurements that a rover can use to navigate.

Today I want to share how I built the core of my system: the GNSS Measurement Engine. This is the component that takes static data and converts it into a dynamic simulation of how satellite signals interact with my rover.

The Purpose: Where Data Meets Physics

The engine's main goal is simple in concept but complex in execution: generate corrected GNSS measurements that are realistic enough to feed into my future end to end system.

In practical terms, this means taking known satellite positions (SP3), antenna corrections (APO), and the rover's true trajectory, then producing:

Corrected pseudoranges: The apparent satellite-to-receiver distance, including simulated errors
Pseudorange rates: For Doppler measurements and velocity estimation
Line-of-sight information: Crucial metadata about each visible satellite

The Architecture: Three Stages of Data Transformation

My implementation follows a three-stage flow that mirrors real GNSS signal processing:

Start with raw satellite data and rover trajectory
↓
Filter out unusable satellites
    - Wrong constellation? Remove
    - Too low elevation? Remove  
    - Missing corrections? Remove
↓
Propagate satellite orbits to exact times
    - Take sparse position points (every 5 minutes)
    - Interpolate to get positions for every second
    - Achieve centimeter-level precision between known points
↓
Generate realistic measurements for each satellite
    - Calculate geometric distance rover-satellite
    - Compute elevation and azimuth angles
    - Apply physical corrections and error models
    - Output final pseudoranges with error breakdowns

Stage 1: Feasibility Filtering

Not all satellites in the data are useful. My filtering process eliminates satellites that don't meet practical operational criteria.

Stage 2: Bridging the Time Gaps with Lagrangian Interpolation

Here's where I faced a fundamental challenge: precise satellite orbits are only provided every 300 seconds in the SP3 files, but my simulation requires positions every single second. Five-minute gaps are eternity in GNSS time.

This is where Lagrangian interpolation becomes essential. The process works like this:

Take 10 known position points around my target time
    - Five points before, four points after the gap
    - Use polynomial fitting to create smooth trajectory
    - Evaluate polynomial at each required second
    - Repeat for X, Y, Z coordinates independently
    - Do the same for satellite velocities

The beauty of this approach is that it doesn't just guess positions, it reconstructs the physical continuity of orbital motion. A 10th-order polynomial captures the subtle accelerations and curvatures that simple linear interpolation would miss. The result is Satellite positions with centimeter-level accuracy at any moment, even between the sparse reference points.

Stage 3: The Measurement Kitchen - Where Physics Comes Alive

The measurement generation is where the true data-to-physics transformation occurs. Most of these values aren't in the raw data. We derive them through computation:

For each satellite that passes our filters:

    Calculate the straight-line distance to rover
    Convert positions to local navigation frame
    Derive elevation angle from vertical/horizontal components
    Compute azimuth for directional relationships

    Apply relativistic timing corrections
    Model atmospheric delays that bend and slow signals
    Account for satellite clock imperfections

    Combine everything into final distance measurements
    Package comprehensive metadata for analysis

Important: We compute elevation and azimuth ourselves because these angles drive many physical models. Atmospheric errors, signal strength, and even satellite selection all depend on these derived values.

Visual Validation: From Numbers to Insight

The real test isn't just producing numbers but producing numbers that make physical sense. I use visualization to transform raw output into intuitive understanding by plotting whatever result I'm interested in. For example:

Create sky plots showing satellite positions
    - Plot each satellite as point on circular chart
    - Position shows direction (azimuth)
    - Distance from center shows elevation
    - Color indicates signal quality

Analyze error patterns vs elevation
    - Plot different error types on separate charts
    - Show how errors decrease at higher elevations
    - Verify models match expected physical behavior

Track measurements over time
    - Graph pseudorange changes for each satellite
    - Highlight when satellites appear/disappear
    - Compare with geometric truth for validation

...

These visualizations serve as my "sanity check". They instantly reveal patterns that would be invisible in spreadsheets of numbers. Seeing atmospheric errors follow expected patterns, or watching satellites move predictably across the sky, confirms the physics is working correctly.

The Crucial Lesson: Understanding Orders of Magnitude

The biggest challenge wasn't implementing algorithms, but developing physical intuition about what the numbers should represent. I had to internalize the expected scales:

Satellite distances: Around 20,000 kilometers
Relativistic effects: Dozens of meters
Atmospheric errors: Meters to tens of meters
Clock errors: Typically 1-2 meters

Without this "numerical sense", debugging would be nearly impossible. The visualizations became my truth-teller, instantly flagging when results drifted from physical reality.

Current Limitations and Next Steps

Currently, hardware/local errors are simulated within realistic bounds. This works for prototyping, but eventually I'll need more sophisticated models based on real estimated data.

The combination of corrected measurements and comprehensive metadata provides both inputs for future filtering and diagnostic tools to understand system behavior.

Conclusion: From Theory to Practice

Building this measurement engine taught me that in precision navigation, the devil is in the physical details. It's not enough to know the equations; you must develop deep intuition about expected magnitudes and interactions. The visualizations bridge the gap between mathematical correctness and physical plausibility.

For those working on positioning systems: What's been your biggest challenge when modeling realistic GNSS measurements? Any counter-intuitive physical effects that surprised you? What visualization techniques have you found most valuable?

As I continue refining this measurement engine, I'm reminded that engineering is rarely a solitary pursuit. The insights and shared experiences from mentors, colleagues and community are invaluable. I'll be diving into the IMU engine implementation next, and I'd love to hear your stories and advice as I take that next step.

The Data Trinity: The Three Files Powering My GNSS Engine

William Pederzoli — Sun, 12 Oct 2025 19:20:54 +0000

Introduction

In my last post, I talked about the "overwhelm" of starting a sensor fusion project. A big part of that initial confusion was simply understanding the data I was supposed to be fusing. It felt like being handed three different maps, each in a unique language and coordinate system, and being told to navigate with them simultaneously.

Before I can even think about fancy algorithms, I need to get the basics right. And that starts with the inputs. So, let's break down the three key file types that are the lifeblood of my project's current phase: the trajectory file, the SP3 file, and the APC file.

Think of this as a prelude to the "GNSS Measurement Engine" I'm building. You can't build an engine without knowing your fuel.

The Rover's Truthful Diary: The Trajectory File

This file (traj_file) is my ground truth, the rover's honest log of where it actually was at every moment. In a real-world scenario, this would come from a highly precise, survey-grade GNSS system. In my simulation-driven stage, it's the absolute reference.

You'll typically find it as a simple text-based file (like CSV) with columns for a time stamp, position (X, Y, Z coordinates or Lat/Lon/Height), and often attitude (roll, pitch, yaw).

Its core purpose is to be the answer key. When my GNSS engine starts estimating positions, I will compare them against this diary to calculate error and validate my system's accuracy. Without this truth, I'd be navigating completely blind.

The Satellite's Precise Dance Card: Navigating SP3 Files

To figure out where you are, you first need to know where the satellites you're listening to actually are. That's where the SP3 file comes in. This file, whose name stands for Standard Product 3, contains the precise orbits of GPS and Galileo satellites.

It’s a standardized text file that acts as a dance card for each satellite, listing its exact position (X, Y, Z coordinates in kilometers in an Earth-Centered frame) and clock correction for specific times.

Using rough, broadcasted orbits introduces significant error. The SP3 file provides "precise orbit" data, which is non-negotiable for high-precision positioning. My engine uses this to get the true source locations for all its calculations.

The Millimeter Adjustment: Correcting with APC Data

Precision is in the details. The position in the SP3 file is for the satellite's center of mass, but the radio signal is generated from its transmitting antenna. The physical offset between these two points is the Antenna Phase Center (APC).

The apc_data file, often from the International GNSS Service (IGS), contains a table that maps satellite types to their specific antenna phase center offsets. It gives us a vector (dX, dY, dZ) to correct the mass-center position and find the true point of signal transmission.

For a basic system, you might ignore this. But for the precision I'm aiming for, these millimeter-to-centimeter corrections are a critical step in reducing systematic errors and respecting the underlying physics of the problem.

Synchronizing the Flow

So, here's the initial data symphony for my GNSS measurement engine: The traj_file tells me where my rover was, the sp3 file tells me where the satellites were, and the apc_data refines those satellite positions to their exact transmission points.

My engine's first job is to consume these three inputs and generate the fundamental "pseudorange" measurements (the estimated distances that form the core of all subsequent positioning math).

Understanding these files wasn't just a bureaucratic step; it was the key to unlocking the entire problem. It transformed the "black box" of GNSS into a structured data processing pipeline.

Now I'm curious for those in the field: Are there other foundational data sources or correction files you found indispensable when building your first positioning engine? What's the one "gotcha" in these standard files that tripped you up? Let me know in the comments.

From Code to Constellations: My first project as a space engineer

William Pederzoli — Mon, 06 Oct 2025 15:03:04 +0000

Introduction

I've always seen myself as a problem-solver. For a decade, those problems lived in electronics and code. But a spark ignited during my master's in GNSS, where I had to bridge two worlds: my background in software and robotics, and the complexity of satellite constellations.

This project is that bridge. It's my first attempt to build a system that takes a rover's trajectory data and satellite data (from GPS and Galileo) and estimates its positions to see how close it gets to "reality." It might sound abstract, but this is the very heart of reliable autonomous navigation.

And, like any first attempt, I hit a wall quickly. It wasn't a specific bug; it was the sheer overwhelm of everything new. Suddenly, I was:

Deciphering the sparse, esoteric structure of RINEX files.
Navigating different coordinate reference frames, trying to understand how to convert one into another.
Wrestling with how the rover's own attitude influences GNSS estimates.

It was like learning to read all over again, but with a dictionary written in a language that was only vaguely familiar. The temptation to give up whispered in my ear, but the problem-solver inside me said, "This is just another system to understand."

The Internal Debate: Choosing the First Tool

My software engineer instinct screamed: "Research, compare, choose the optimal tool!" I knew the end goal was to implement a Kalman Filter, the gold-standard tool for this job. But I also knew you don't start by running a marathon; you learn to walk first.

So, what was my decision? To start with a simple Lagrangian interpolation.

It's not the most glamorous solution, but it's a solid foundation. It allows me to understand the data flow, validate the inputs and outputs, and, most importantly, establish a baseline. This baseline will be crucial for measuring the improvement once I implement the Kalman Filter later. Sometimes, the best choice is the one that lets you make progress today, not the perfect one you might build tomorrow.

The Current State: The "It Works!" Moment

Right now, the project is in what I'd call a "promising alpha" stage. The code runs. It can generate GNSS measurement estimates, produce Line-of-Sight (LOS) and Position (POS) files, and plot the initial results.

The most rewarding part? Seeing the results fall within the expected orders of magnitude. That first plot, the one that actually resembles reality, is the biggest motivational boost you can get. It’s a tangible sign that you're on the right path.

The Open Question: Looking Ahead

The foundation is laid. Now, the path forks into multiple optimization branches. The Kalman Filter is the obvious next step, but I know there's more to robustness than just one algorithm.

So, I'd like to open the floor to those who have walked this path before. If I could ask one thing to an expert:

"Beyond implementing the Kalman Filter, what single improvement had the most significant impact on the performance or robustness of your sensor fusion system when you moved from a working prototype to a robust solution?"

Was it a specific outlier detection method? A particular way to model sensor errors? A data fusion architecture that proved exceptionally resilient?

I'm all ears for your war stories and wisdom. Let me know in the comments.