A lightweight amateur satellite Doppler shift correction algorithm based on Chebyshev polynomials
500-byte Flash memory, SGP4 comparable accuracy
Abstract
Multi-point frequency shift correction in low-Earth orbit satellite communication typically requires running an SGP4 orbit propagator, which consumes over 30 KB of Flash memory for its code and ephemeris data. This is difficult to accommodate on resource-constrained Cortex-M0 microcontrollers (commonly used amateur radio firmware). This paper proposes a lightweight alternative: using a Chebyshev polynomial family to fit the normalized radial velocity curve, with transit duration, satellite altitude, and maximum elevation angle as parameters. By using Q1.16 fixed-point format, the algorithm + coefficients require only approximately 500 bytes of Flash storage, eliminating reliance on floating-point libraries. In a one-satellite cross-validation experiment over 10 FM amateur satellites, 1420 transits, and a 30-day window, the average maximum error in the VHF band was 397 Hz, and in the UHF band it was 868 Hz, fully meeting the ±5 kHz frequency deviation tolerance for narrowband FM. Further generalization experiments showed:（1）The error shows no significant degradation within 90 days, and after one year, the deviation from the normalized curve shape only increases by 6%.（2）The models trained in Taian were directly applied to locations across six different latitudes, with UHF error deviations remaining within 13%.（3）The estimated error for Δf_max increases by 75% after one year (primarily due to orbital decay, requiring regular updates to the ephemeris), but the total error remains within the usable range of FM. The algorithm only requires the user to input three parameters: satellite selection, AOS time, and LOS time.
1. Introduction
1.1 Background of the Problem
Amateur radio enthusiasts use FM satellites (such as SO-50, AO-91, and ISS inter-band repeaters) for voice communication. These satellites operate in near-circular orbits at altitudes of 400–1400 km, with a speed of approximately 7.5 km/s. The resulting Doppler shift can reach ±3.5 kHz at 145 MHz (2m band) and ±10.5 kHz at 435 MHz (70 cm band). Without real-time frequency correction, the received signal will drift outside the receiver's passband within a few seconds.
Recently, while developing the firmware for the Baofeng UV-K61, we discovered that running the SGP4 orbit predictor on a low-end amateur radio microcontroller requires storing the complete TLE (Two-Line Ephemeris) and implementing the SGP4 algorithm. This resulted in Flash memory usage exceeding 30 KB (considering the overhead of soft-float code). For Cortex-M0 radios with only 64-128 KB of Flash, this presented a difficult choice for developers: either abandon other radio functions to run a dedicated SGP4 firmware, or use a computer to pre-download the necessary data.
1.2 Existing solutions and their shortcomings
Currently, open-source firmware primarily employs two strategies, taking Quansheng UV-K5 as an example (excluding high-performance DMR radios):
- Full-featured SGP4 firmware: Dedicate all Flash functionality to SGP4, and disable other radio functions.
- Pre-calculated lookup table（LUT）: Before each passage, upload the correction table from the PC via serial port. A computer is required, and the device needs to be re-programmed each time.
1.3 Contributions of this paper
We have demonstrated that the Doppler curves observed during satellite transit can be obtained fromA family of low-order Chebyshev polynomialsPrecise approximation, requiring only three parameters known to the user – satellite identity, AOS time, and LOS time. This algorithm:
- Using Q1.16 Fixed-point number, eliminate floating-point library dependencies, requiring only Approximately 500 bytes Flash (saves more than 100 times compared to SGP4)
- Achieve 397 Hz（VHF）868 Hz（UHF） Mean error (vs. true values from SGP4)
- Each frequency update requires only 30 CPU cycles(Cortex-M0 integer instructions)
- No need to upload data each time you pass through.— Users simply need to select a satellite and enter the transit time.
- Shape of the normalized curveTime-invariant— After one year, the shape error only decreased by 6%.
- Global generalization— A set of coefficients applies to all frequencies (UHF offset)<13 %）
2. Data and Methods
2.1 Data Sources
TLE Star Chart: August 3, 2026, from CelestraK amateur radio group Obtain the latest TLE data for 10 commonly used amateur FM satellites.
Observatory: Located in Tai'an City, Shandong Province (36.18° N, 117.13° E), with a minimum observed angle of 5°.
Analyze satellite data：
| Satellite | NORAD | Downward frequency | Height | Number of crossings in 30 days |
|---|
| SO-50 (SaudiSat-1C) | 27607 | 436.795 MHz | 627 km | 163 |
| AO-91 (RadFxSat/Fox-1B) | 43017 | 145.960 MHz | 533 km | 109 |
| AO-85 (Fox-1A) | 40967 | 145.980 MHz | 601 km | 138 |
| AO-95 (Fox-1Cliff) | 43770 | 145.940 MHz | 490 km | 112 |
| LilacSat-2 (CAS-3H) | 40908 | 437.200 MHz | 415 km | 103 |
| ISS (Zarya) | 25544 | 437.800 MHz | 426 km | 208 |
| PO-101 (Diwata-2) | 43678 | 437.500 MHz | 572 km | 150 |
| AO-73 (FUNcube-1) | 39444 | 145.950 MHz | 543 km | 118 |
| NO-44 (PCSAT) | 26931 | 145.825 MHz | 792 km | 157 |
| RS-44 | 44909 | 145.935 MHz | 1350 km | 162 |
Total: 1420 transits (796 VHF + 624 UHF), within a 30-day window (August 3, 2026 to September 2, 2026).
2.2 Truth Value Generation
Pass through each orbit using the SGP4 model.skyfield Python library (encapsulating standard)sgp4(Package) Propagate TLE, calculate true Doppler shift at 150 evenly spaced time sampling points:
Δftrue(t) = -fCarrier $\cdot \frac{v_r(t)}{c}$ (1)
Specifically, $v_r(t)$ represents the radial velocity (positive when moving away from the observer), which is obtained by calculating the difference in relative position between the satellite and the observer in the geocentric coordinate system. Also, $c = 299\;792\;458$ m/s is the speed of light.
2.3 Key Insights: Unrelatedness to Normalization and Frequency
The key physical insights of this method are:The shape of the Doppler curve is independent of the carrier frequency.. Unifying Doppler shift measurements from AOS time series:
$$\text{norm}(t_{\text{norm}}) = \frac{\Delta f(t)}{\Delta f(\text{AOS})} = \frac{v_r(t)}{v_r(\text{AOS})} \tag{2}$$
The carrier frequency $f_{\text{carrier}}$ is precisely compensated for in the normalization. The normalization curve depends only on the satellite's orbital geometry – specifically:
- Satellite altitude $h$: Determines the orbital angular velocity $\omega = \sqrt{GM/(R+h)}$3}$
- Maximum angle of elevation $\theta_{\max}$: Determines the symmetry of the crossing and the steepness of the S-curve near the TCA.
Inference: Any error $\epsilon$ in modeling the normalized radial velocity, converted to a frequency error as $\epsilon \cdot f$.{\text{carrier}}/c$. Due to $f$The ratio of UHF to VHF error, i.e., UHF/f<sub>VHF</sub> ≈ 3, implies that the UHF error is inherently three times greater than the VHF error.
2.4 Approximate Model
We evaluated six models with increasing complexity:
Model A: Cosine Baseline
$$\Delta \hat{f}(t) = \Delta f$${max} ⋅ cos(π ⋅ t)(norm), t{\text{norm}} = \frac{t - t}{AOS} tLOS - t{AOS} ∈ [0, 1] (3)
Model B: Polynomial with a single term (degree $d$)
$$\Delta \hat{f}(t) = \Delta \hat{f}$${max} ⋅ Σ∑_{i=0}^{d} c\_i * t\_normi (4)
The coefficients $c_i$ are obtained by fitting a normalized SGP4 curve using least squares. The values are calculated using the Horner method.
Model C: Chebyshev polynomials (degree $d$)
$$\Delta \hat{f}(t) = \Delta \hat{f}$${max} ⋅ Σ∑_{i=0}^{d} a_i * T_i(x) , where x = 2 t_{\text{norm}} - 1 ∈ [-1, 1] (5)
Among them, $T_i(x) = \cos(i \cdot \arccos(x))$ is the first-order Chebyshev polynomial on $[-1, 1]$. The coefficients $a_i$ are obtained through least squares fitting. The evaluation method is adopted.Clenshaw's backward recursion：
$$\begin{cases} b{d + 2} = b{d+1} = 0 \ b\_k = a\_k + 2x * b{k+1} - b{k+2}, where k = d, d-1, ..., 1 \ \text{norm} = a_0 + x \cdot b_1 - b_2
Clenshaw's backward difference method, when used with float32 precision on the Cortex-M0, provides more stable numerical results compared to the Horner method of the same order.
Model D: 3-bin oblique Chebyshev-7
The passing curves were divided into three groups based on their maximum angle of elevation (low/medium/high, using the 33rd and 67th percentiles of the training set as split points), and each group was fitted independently with a Chebyshev-7 polynomial.
$$\Delta \hat{f}(t) = \Delta \hat{f}$${max} ⋅ Σ{i=0 to 7} a_i^{(\text{bin})} \cdot T_i(2t)(norm - 1), bin = {low} & θ{max} ≤ θ{33} \ middle & θ{33} < θ{max} ≤ θ{67} \ high & θMaximum > θ{67} \end{cases} (7)$$
Maximum Δf estimation model
All polynomial models require an estimate of the maximum Doppler shift, $\Delta \hat{f}_{\max}$. We useHigh-order perceptual model：
$$\Delta \hat{f}_{max} = \alpha_1 T$$2 + α₂T + α₃Th + α₄h + α₅ (8)
In this case, $T = t$LOS - tAOS is the transit duration (in seconds), and h is the satellite altitude (in km). This model captures the following physical facts: both the transit duration and the peak Doppler shift depend on the satellite altitude and the maximum elevation angle. Knowing the satellite's identity (and therefore its altitude) eliminates the ambiguity of "different satellites having the same transit duration corresponding to different Δf_max."
2. Verification Method
Main experiment: UsingLeave-one-satellite Cross-Validation (LOSO CV): For each satellite $s$ in the dataset:
- Use all transits of satellite $s$ as the test set.
- Train all models using the data from the remaining 9 satellites.
- Evaluate along each transit path of satellite $s$.
This ensures that the data from the test satellite is accurate.Never encountered during trainingIt provides a realistic assessment of the generalization capabilities of unknown satellites.
Generalization experiment：
- 90-day generalizationTrained on data from August 2026 (days 1-30), tested on data from September-October 2026 (days 31-90)
- One year generalization: Trained on data from August 2026, tested on data from August-September 2027 (days 365-395), to assess the robustness of polynomial coefficients against TLE aging and orbital decay.
Evaluation metricsCalculate the prediction error at each test path using 150 sampling points and compare it to the true value obtained from SGP4.Maximum absolute errorReport the mean, median, P95, and worst-case value for each crossing boundary error.
3. Experiment and Results
3.1 Why the Cosine Model Failed
Figure 1: Comparison of Cosine Model and SGP4
Figure 1a–b The cosine models for the high-angle passes of SO-50 (436.8 MHz) and AO-91 (146.0 MHz) were compared with the true values from SGP4. The cosine curves showed a fundamental mismatch with the shape of the SGP4 model: it was not steep enough near TCA (when the satellite is closest to the horizon and the sky motion speed is fastest), and it was too steep near the horizon (when the satellite's apparent motion is slow). The error curve (green) shows typical "double peak" residuals. The worst case for SO-50 reached 3.3 kHz, which is 66% of the FM frequency deviation.
Figure 1c The 141 normalized Doppler curves were overlaid. The gray curves representSpreadRepresents any irreducible shape variation that cannot be captured by a single polynomial. The black solid line represents the Chebyshev-7 mean.
Figure 1d DemonstratedMulti-degree polynomial platform effect: For VHF and UHF bands, the CV error stops decreasing after the 7th order. The error floor (approximately 500 Hz for VHF and 1440 Hz for UHF) represents the spread in shape between curves – rather than a lack of flexibility in the polynomial. The CV curve closely matches the in-sample curve, confirming that overfitting does not occur even at the 15th order.Chebyshev and Monomial have identical accuracy when they are of the same degree.— Both use the same polynomial space, but Chebyshev coefficients decay faster, resulting in greater numerical stability under float32.
3.2 The Physical Basis of Normalizing Doppler Curves
Figure 2: Physical Mechanism and Angle Dependence
Figure 2a Demonstrates the geometry of transit: the radial velocity $v_r(t)$ is the projection of the satellite's orbital velocity vector onto the line of sight of the observer. At AOS and LOS times (low elevation angles), the velocity is almost perpendicular to the line of sight, resulting in slow Doppler shifts. Near TCA (high elevation angle), the velocity aligns with the line of sight, leading to rapid Doppler shifts.
Figure 2b VerifiedThe 3x zoom rule: Multiply the 145 MHz Doppler curve by 3 to precisely match the 435 MHz curve under the same geometric conditions. This confirms the frequency independence of $\text{norm}(t) = v_r(t)/v_r(\text{AOS})$.
Figure 2c Revealed a key finding:For a single satellite（SO-50）The shape of the normalized curve changes dramatically with the maximum angle.The crossing height asymmetry at a 9° angle of attack (TCA significantly off-center), while the crossing height was nearly symmetrical at 87°. The shape difference reached.47% at full scale。
Figure 2d Color all normalized curves according to the maximum angle of incidence, which visually demonstrates that the angle of incidence is the primary driver of changes in curve shape. High-angle transitions (yellow) are clustered around the symmetrical S-shaped curves, while low-angle transitions (purple) are more dispersed.
3.3 Model Comparison
Figure 3: Model Comparison
Figure 3a Demonstrated the mapping relationship between transit time and the Doppler peak. The high-resolution secondary model achieved an R value of approximately 70 cm in the 70 cm band.2 = 0.982 (RMS residual at 397 Hz), while a band-level quadratic model without high-frequency information had an R of2 = 0.889$. The estimated maximum residual values (Δf_max) for each satellite are shown in Table 1.
Table 1: Estimated residual values of Δf_max (70 cm bandwidth, height-aware model)
| Satellite | Mean residual | Standard deviation | Maximum residual |
|---|
| ISS | -318 Hz | 185 Hz | 635 Hz |
| LilacSat-2 | +482 Hz | 133 Hz | 687 Hz |
| PO-101 | +340 Hz | 260 Hz | 782 Hz |
| SO-50 | -212 Hz | 246 Hz | 614 Hz |
Figure 3b Compared the modelsShape error(CV, assuming perfect Δf_max):
| Model | Flash | 2m average | Average of 70 cm |
|---|
| Cosine | 8 B | 1013 Hz | 3305 Hz |
| Poly-3 | 16 B | 600 Hz | 1931 Hz |
| Poly-7 (Baseline) | 32 B | 499 Hz | 1440 Hz |
| 3-bin Cheb-7 | 96 B | 313 Hz | 714 Hz |
| SGP4 (full) | >10 KB | 0 Hz | 0 Hz |
From Cosine to Poly-7, the shape error was reduced by 50% (2m) / 56% (70cm). From Poly-7 to 3-bin, the reduction further decreased by 37% (2m) / 50% (70cm).
Figure 3c DemonstratedTotal error CDF(Shape + Estimated Δf_max). Key statistics are shown in Table 2.
Table 2: Total Error Budget (3-bin Cheb-7 + Height-Aware Δf_max)
| Metrics | 2m（VHF） | 70 cm（UHF） |
|---|
| Average maximum error | 397 Hz | 868 Hz |
| Maximum median error | 373 Hz | 905 Hz |
| Worst-case scenario | 883 Hz | 1753 Hz |
| P95 | 707 Hz | 1464 Hz |
| < 500 Hz | 79 % | 18 % |
| < 1000 Hz | 100 % | 67 % |
| < 2000 Hz | 100 % | 100 % |
Figure 3d Place each scheme at the forefront of Pareto with Flash precision. The 3-bin Chebyshev-7 (96 bytes) is located at the "knee point" on the front end.
3.4 Generalization Experiment
3.4.1 90-day generalization
The coefficients, trained using data from August 2026 (days 1-30), were tested on data from September-October (days 31-90).
Table 3: Generalization results over 90 days (3-bin Cheb-7)
| Dataset | Average maximum error of 2m | 2m - worst | Average maximum error of 70 cm | 70 cm is the worst |
|---|
| TRAIN (August) | 397 Hz | 883 Hz | 868 Hz | 1753 Hz |
| TEST-1 (September) | 390 Hz (-2%) | 872 Hz | 866 Hz (-0%) | 1820 Hz |
| TEST-2 (October) | 381 Hz (-4%) | 925 Hz | 922 Hz (increased by +6%) | 1697 Hz |
The errors observed in each satellite over the two test months were almost identical to those observed during the training month. Within 90 days.No apparent generalization degradation.— This demonstrates that the physical nature of the normalized radial velocity curve is time-invariant: the geometric relationship between satellites passing through does not change with TLE aging.
3.4.2 Generalization over one year
Test on data from August-September 2027 (TLE age 365-395 days) using the coefficients trained in August 2026.
Table 4: Generalization results over one year (3-bin Cheb-7)
| Dataset | Total error of 2m | 2m shape error | Estimated error in 2m Δf | Total error of 70 cm | 70 cm shape error | Estimated error in Δf (70 cm) |
|---|
| TRAIN（2026.08） | 397 Hz | 287 Hz | 234 Hz | 868 Hz | 691 Hz | 416 Hz |
| +1 year（2027.08） | 413 Hz (increased by +4%) | 308 Hz (increased by +7%) | 244 Hz (±5%) | 1165 Hz (+34%) | 729 Hz (increased by +6%) | 730 Hz (increased by +75%) |
| January 2023（2027.09） | 440 Hz (increased by +11%) | 298 Hz (increased by +4%) | 277 Hz (increased by +18 %) | 1220 Hz (+41%) | 725 Hz (with a +5% increase) | 806 Hz (increased by +94%) |
Table 5: Degradation details for each satellite over one year (70 cm band)
| Satellite | Height | AVERAGE TRAIN | +1 year, average | Degradation rate |
|---|
| ISS | 426 km | 800 Hz | 761 Hz | −5 %* |
| LilacSat-2 | 415 km | 1049 Hz | 2004 Hz | +91 % |
| PO-101 | 572 km | 767 Hz | 977 Hz | +27 % |
| SO-50 | 627 km | 804 Hz | 907 Hz | +13 % |
*The ISS has regular orbit maintenance (reboosting), so its altitude remains almost constant.
The results indicate that,
The shape of the normalized curve is almost time-invariant.: The shape errors of 2m and 70cm increased by only 7% and 6% respectively after one year. This means that the Chebyshev polynomial coefficientsNo update requiredBecause they encode purely physical relationships between objects passing through a space.
The degradation is mainly caused by the drift in the estimated value of Δf_max.: The estimated error in the Δf for the 70 cm bandwidth increased from 416 Hz to 730 Hz (+75%), accounting for most of the total degradation. The root cause isTrack degradation, low-Earth orbit satellites (particularly<The LilacSat-2, with a 500 km altitude and a degradation rate of +91%, experiences significant atmospheric drag, which causes its orbital height to slowly decrease. This alters the relationship between its transit time and peak Doppler shift. In contrast, the ISS (426 km), due to its regular orbit maintenance, exhibits a degradation rate of -5%.
Instructions for practical deployment: To maintain accuracy, re-estimate the Δf_max coefficient (5 floats, 20 bytes) every 6–12 months using the latest TLE. The shape polynomial coefficients can be used permanently.
3. 5 Final Algorithm
Figure 4: The final algorithm
Figure 4a Demonstrates the effects of three sets of Chebyshev-7 fitting segments overlaid on a normalized curve collection, each at different angles. Each set of polynomials captures the characteristic shape within its angle range: high angles result in flatter S-curves, while low angles produce more asymmetric curves.
Figure 4b Demonstrates the overall error distribution for each satellite under the LOSO CV (Constellation Volumetric) model. The differences between satellites are minimal, validating the robustness of the model.
Figure 4c Taking SO-50 (which was not part of the training data in CV) as an example, we demonstrate the excellent consistency between the 3-bin model and the true values of SGP4 over multiple elevation angles.
Figure 4d Place all errors within the context of ±5 kHz FM frequency deviation: Even the worst-case scenario (1753 Hz) represents only 35% of the FM frequency deviation, and is entirely within the typical AFC capture range of a standard FM receiver (±3 kHz).
Algorithm pseudocode
// One-time initialization (written to Flash during firmware burn)
// Chebyshev-7 coefficients (6 groups × 8 floats = 192 bytes)
// dur→Δf_max estimation coefficients (40 bytes)
// Satellite altitude table (10 satellites × 1 float = 40 bytes)
// Initialization on each pass
function doppler_init(aos_unix, los_unix, sat_id):
g_aos = aos_unix
g_dur = los_unix - aos_unix // Transit duration
h = sat_altitude_table[sat_id] // Lookup satellite altitude
d = g_dur
g_df_max = DUR_COEFF[0]*d² + DUR_COEFF[1]*d // Formula(8)
+ DUR_COEFF[2]*d*h + DUR_COEFF[3]*h + DUR_COEFF[4]
θ_est = estimate_elevation(g_dur, h) // Estimate the maximum elevation angle
if θ_est ≤ θ_LOW: g_cheb = cheb_coeffs_low
elif θ_est ≤ θ_HIGH: g_cheb = cheb_coeffs_mid
else: g_cheb = cheb_coeffs_high
// Call every second (during transit)
function doppler_update(now_unix):
t_norm = (now_unix - g_aos) / g_dur
if t_norm < 0 or t_norm > 1: return 0 // The transition has not started or has ended.
x = 2 * t_norm - 1 // [0,1] → [-1,1]
// Clenshaw's recurrence formula(6)
bk2 = 0; bk1 = 0
for k = 7 down to 1:
bk = g_cheb[k] + 2 * x * bk1 - bk2
bk2 = bk1; bk1 = bk
norm = g_cheb[0] + x * bk1 - bk2 // ∈ [-1, 1]
return norm * g_df_max // Hz
The computational cost for each update.: 1 division (t_norm) + Clenshaw's iterative algorithm (7 multiplications + 7 FMA operations) ≈ approximately 30 cycles on a Cortex-M0.
3.6 Fixed-point number implementation
Cortex-M0 does not have a hardware floating-point unit.（FPU）。 Any float The operations are linked to a software simulation library, which significantly increases Flash usage (Table 5).
Table 5: Size of Arm GCC Cortex-M0 software floating-point library
| Function | Applications | Volume |
|---|
__addsf3 | Floating-point addition | 504 B |
__subsf3 | Floating-point subtraction | 428 B |
__mulsf3 | Floating-point multiplication | 416 B |
__divsf3 | Floating-point division | 524 B |
__fixsfsi | float → int | 112 B |
__floatsisf | int → float | 168 B |
CMPSF2 | Floating-point comparison | 180 B |
| Total | | 2,332 B |
This 2.3 KB floating-point library is a true Flash killer – it's orders of magnitude larger than the data and code associated with the algorithm itself.
3.6.1 Q1.16 Fixed-point format
We use Q1.16 Fixed-point format eliminates reliance on floating-point numbers. All values are stored as signed 32-bit integers, with the lower 16 bits representing the fractional part:
$x$Q16 = round(x * 2^16), where x = (x)Q16 = 2^16
- Possible values: [-32768, 32768 - 2^(-16)] (Covers normalized Doppler [-1, 1] and Chebyshev mean)
- Resolution: $2^{-16} \approx 1.53 \times 10^{-5}$ (which is equivalent to 0.15 Hz at a maximum frequency difference of 10 kHz)
- Multiplication：
(int64_t) a * b >> 16(The Cortex-M0 has) SMULL/UMULL Hardware instructions (no library required)
- Clenshaw's midpoint value：
2 * x * bk1 The three-factor product needs... (int64_t) 2 * x_q * bk1 >> 32, 64-bit intermediate results do not overflow
3.6.2 Bit-width scanning
Figure 5: Point-wise quantization analysis
Figure 5a Demonstrates the relationship between quantization error and the number of decimal places. Within the Q8-Q16 range, quantization noise decreased from a steady 46 Hz to 0.2 Hz (70 cm).Q18 and above, Δf_max value (10 kHz x 2^18 > The use of $2^{31}$ as the exponent leads to an integer overflow, resulting in a rapid increase in error. Q16 represents the optimal balance between precision and avoiding overflow.
Figure 5b The comparison showed the difference in Flash usage. The floating-point scheme used approximately 2.7 KB (of which 2.3 KB was for the floating-point library), while the fixed-point Q16 scheme required only about 500 bytes (without relying on a floating-point library).Save 91%.
Figure 5c Demonstrates the quantized noise distribution of Q16 (compared to float32 reference). VHF standard deviation is 0.15 Hz, UHF standard deviation is 0.50 Hz – compared to the inherent 400–900 Hz error of the model itself.Completely negligible。
Figure 5d The fixed-point scheme is placed at the forefront of the Flash-precision Pareto front, with Q16 3-bin Cheb-7 located at the optimal knee point.
3.6.3 Q16 Fixed-point coefficient table
Table 6: Chebyshev-7 coefficients (70 cm band) for Q1.16 format
| bin | c0 | c1 | c2 | c3 | c4 | c5 | c6 | c7 |
|---|
| low (5–22 degrees) | 15 | −69374 | −61 | 5090 | 22 | −797 | −7 | 178 |
| middle (22–55 degrees) | −34 | −77378 | 17 | 16074 | −6 | −5089 | 5 | 2407 |
| High (55–90 degrees) | −7 | −79773 | 16 | 20900 | −10 | −8055 | 11 | 5000 |
Multiply by $2^{-16}$ to restore the floating-point value. The complete Q16 coefficients and the "C" header file are shown in Appendix C.
3.6.4 Fixed-point Clenshaw Algorithm (C implementation)
// Q1.16 Clenshaw - No floating-point dependency, native Cortex-M0 hardware instructions
static inline int32_t clenshaw_q16(const int32_t c[8], int32_t x_q) {
int32_t bk2 = 0, bk1 = 0, bk;
int32_t two_q = 2 << 16;
for (int k = 7; k < 16; k++) { > 0; k--) {
// (int64)2 * x_q * bk1 >> 32 SMULL hardware instruction, no library calls
int32_t term = (int32_t)(((int64_t)two_q * x_q * bk1) >> 32);
bk = c[k] + term - bk2;
bk2 = bk1; bk1 = bk;
}
int32_t term_f = (int32_t)(((int64_t)x_q * bk1) >> 16);
return c[0] + term_f - bk2;
}
// Doppler update: Called every second during transit, returns an integer Hz
int32_t doppler_update(uint32_t now_unix) {
int32_t elapsed = (int32_t)(now_unix - (uint32_t)g_aos);
if (elapsed <= 0 || elapsed >= g_dur_sec); return 0;
int32_t t_q = ((int64_t)elapsed << 16) / g_dur_sec);
int32_t x_q = (t_q << 1) - (1 << 16); // x = 2t - 1
int32_t norm_q = clenshaw_q16(g_cheb, x_q);
return (((int64_t)norm_q * g_df_max_q); >> 16); // Hz
Each update requires only 7 times. SMULL(64-bit multiplication) + approximately 20 ALU instructions, without any library calls.
3.7 Generalization Across Locations
To assess the model's dependence on observer location, the model was trained using data from Taian (36°N), and then tested at six different locations around the world with varying latitudes (see Figure 6).
Figure 6: Cross-location generalization experiment
Table 7: Generalization results across locations (3-bin Cheb-7 Q16)
| Location | Latitude | Average of 2 meters | 2 meters (worst) | Average of 70 cm | 70 cm of worsted yarn | vs Tai'an (70cm) |
|---|
| Taian (Baseline) | 36° North | 377 Hz | 575 Hz | 815 Hz | 1428 Hz | — |
| Singapore | 1° North | 438 Hz | 820 Hz | 811 Hz | 1507 Hz | −0.4 % |
| Helsinki | 60 degrees North | 378 Hz | 625 Hz | 830 Hz | 1307 Hz | +1.9 % |
| San Diego | 33° South | 469 Hz | 872 Hz | 881 Hz | 1450 Hz | +8.1 % |
| Ancrechy | 61° North | 410 Hz | 698 Hz | 873 Hz | 1476 Hz | +7.1 % |
| Sydney | 34° South | 465 Hz | 836 Hz | 916 Hz | 1287 Hz | +12.4 % |
Figure 6a–b The error distributions for each location in the UHF and VHF frequency bands were presented. All locations in the 70 cm band were within 13% of the Tai'an baseline, while all locations in the 2 m band were within 25%.Figure 6c Display deviation percentage.Figure 6d With latitude as the horizontal axis, the error in the southern hemisphere is slightly higher than that in the northern hemisphere (approximately +100 Hz), primarily due to the slight asymmetry caused by the geometry of the southern hemisphere passing through. However, this is still far within the acceptable range for FM.
The excellent physical reasons, as explained in section 2.3, have been clarified: Normalization eliminates frequency dependence, and the main differences between different latitudes (transit zenith angle distribution) are covered by a 3-bin zenith angle segmentation model. The difference in linear velocity due to Earth's rotation between different latitudes (460 m/s at the equator vs 0 m/s at the poles) only results in approximately 6% correction of the observed velocity, and this difference is essentially eliminated through normalization.
4. Conclusion
4.1 Summary
This article demonstrates that3-bin oblique-angle segmented Chebyshev-7 polynomialCooperationHighly perceptive secondary Δf_max estimator, enabling sub-kHz precision FM satellite Doppler correction on resource-constrained microcontrollers. By using Q1.16 fixed-point format to eliminate reliance on floating-point libraries, the algorithm and coefficients require only approximately 500 bytes of Flash memory (more than 100 times less than SGP4), without requiring data uploads each time a pass occurs. Each frequency update requires only approximately 30 CPU cycles. The model can be globally generalized after training in Tai'an, with UHF error deviations for 6 locations at different latitudes.<13%. A one-year generalization experiment showed that the normalized curve shape is almost time-invariant. The maximum Δf value should be estimated every 6–12 months using the latest TLE update.
4.2 Key Findings
- The shape of the Doppler curve represents the normalized radial velocity.: The carrier term is normalized, and a set of polynomial coefficients applies to all frequencies.
- UHF error = 3 × VHF errorThis is a necessary consequence of the Doppler physical law, and any model (including a complete SGP4) cannot avoid it.
- The Chebyshev and Monomial precision are completely consistent.: Chebyshev is recommended because Clenshaw's iterative method provides more stable numerical results in float32, and the coefficients decay quickly, allowing for safe truncation.
- The maximum angle of elevation is the strongest predictor factor for the curve shape.: Segmenting by aspect angle (3-bin) reduces shape error by 50% compared to a single polynomial.
- Add satellite altitude informationEstimate Δf_max and then use it to calculate R2The ambiguity caused by different satellite altitudes has been eliminated, increasing from 0.89 to 0.98 (70 cm wavelength).
- The TCA information gain is limited.: TCA offset is related to pitch angle, but pitch angle is a direct driver of shape; TCA is merely a byproduct.
- The shape of the normalization curve is time-invariant.: After one year, the shape error only degraded by 6%, and the Chebyshev polynomial coefficients could be used permanently.
- The estimated maximum Δf (frequency difference) decreases with orbital decay and degradation.: The estimated error in Δf is expected to increase by 75% (70 cm) within one year. It is recommended to update the latest TLE every 6-12 months (using 5 Q16 coefficients and a 20-byte data).
- Point-by-point quantification with minimal loss.Q1.16 format: Quantized noise is only 0.5 Hz.（UHF）The difference between this and the model error is negligible. After removing the floating-point library, Flash was reduced from 2.7 KB to 500 bytes.
- Global generalization does not require retraining.: The coefficients trained in Tai'an (36°N) are directly used for calculations from the equator to 60°N/S latitude, with an error of ±70 cm.<13 %.
4.3 Practical Deployment Recommendations
The algorithm is suitable for integration into amateur radio firmware.
- One-time programming: Chebyshev Q16 coefficients (6 groups x 8 x 4 = 192 bytes), dur2Δf Q16 coefficients (2 groups x 5 x 4 = 40 bytes), satellite height table (10 satellites)16=160 bytes), code ((180 bytes). Total approximately 572 bytes (including satellite table), without reliance on floating-point libraries.
- Setting for each crossingUser selects satellite from menu and enters AOS and LOS times. The MCU estimates Δf_max and the maximum angle of elevation based on duration and altitude, and then chooses the corresponding Chebyshev array.
- During transit: Calculates $t_{\text{norm}}$ (Q16 integer division) per second, evaluates the Chebyshev polynomial using Clenshaw's algorithm (7x64-bit multiply-add), multiplies by Δf_max (1x64-bit multiply+shift), and applies frequency correction. All operations use Cortex-M0 hardware integer instructions.
4.4 Limitations
- This model is designed forNear-circular LEO orbitVerification (altitude 400–1400 km, eccentricity)<(0.05). Highly eccentric orbits (such as AO-10, AO-7) require separate treatment.
- A 3x UHF/VHF error ratio is an inherent constraint of Doppler physics, and cannot be reduced by any model. It can only be improved by using a lower frequency band.
- Almost reached the top>The rapid Doppler shift (up to 85 degrees) poses a significant challenge for any polynomial model. The 3-bin scheme mitigates this problem through the use of specialized high-angle polynomial components.
- The TLE needs to be updated every 6–12 months to maintain the accuracy of the estimated Δf_max value.
- The error in the Southern Hemisphere is slightly higher than that in the Northern Hemisphere (+12%), but the absolute difference (100 The frequency of 1 Hz has no significant impact on FM communication.
Appendix A: Chebyshev coefficient table
A.1 70 cm band, 3-bin Chebyshev-7
// Low angle (5–22°), CV average maximum error 569 Hz
static const float cheb7_70cm_low[8] = {
+0.0002335556, -1.0585669335, -0.0009278245, +0.0776628764,
+0.0003356745, -0.0121636282, -0.0001125546, +0.0027210660
};
// Medium angle (22–55°), CV average maximum error 752 Hz
static const float cheb7_70cm_mid[8] = {
-0.0005193342, -1.1806953138, +0.0002602060, +0.2452703434,
-0.0000972926, -0.0776532436, +0.0000699958, +0.0367319887
};
// High angle (55–90°), CV average maximum error 820 Hz
static const float cheb7_70cm_high[8] = {
-0.0001019328, -1.2172468235, +0.0002460633, +0.3189162096,
-0.0001456953, -0.1229115857, +0.0001671119, +0.0762920805
};
A.2 70 cm bandwidth – dur → Δf_max estimation coefficient
// Δf_max = c[0]·dur² + c[1]·dur + c[2]·dur·alt + c[3]·alt + c[4]
static const float dur2df_70cm[5] = {
0.00114662f, 35.69112240f, -0.03579820f, 1.84041431f, -978.11635070f
};
A.3 Satellite Altitude Reference Table
// Calculate a from the average movement of TLE: a = (GM/(n/60))²(1/3), h = a - 6371 km
static const float sat_altitudes[10] = {
// SO-50 AO-91 AO-85 AO-95 LilacSat-2 ISS PO-101 AO-73 NO-44 RS-44
627, 533, 601, 490, 415, 426, 572, 543, 792, 1350
};
A.4 Angle of elevation threshold
2 MHz band: $\theta${low} ≤ 26.1•$, $\theta{high} > 56.7•$
70 cm bandwidth: $\theta${low} ≤ 22.1•$, $\theta{high} > 54.7•$
Appendix B: Clenshaw's Backward Difference C implementation
// Evaluate the Chebyshev series for t_norm ∈ [0, 1]
// coeffs[0..7] correspond to T_0 to T_7 on x ∈ [-1, 1]
float chebyshev_eval(const float coeffs[8], float t_norm) {
float x = 2.0f * t_norm - 1.0f; // [0,1] → [-1,1]
float bk2 = 0.0f, bk1 = 0.0f, bk;
for (int k = 7; k >= 0; k--) {
bk = coeffs[k];
bk2 = bk * bk1 + bk * coeffs[k+1];
bk1 = bk;
x = bk2 + x;
}
return x;
} > 0; k--) {
bk = coeffs[k] + 2.0f * x * bk1 - bk2;
bk2 = bk1;
bk1 = bk;
}
return coeffs[0] + x * bk1 - bk2; // ∈ [-1, 1]
}
// Complete Doppler update — call every second during transition
float doppler_update(uint32_t now_unix) {
float t_norm = (float)(now_unix - g_aos_time) / g_duration_sec;
if (t_norm < 0. 0f || t_norm > 1.0f) return 0.0f; // Inactive
// Select Chebyshev basis array based on the elevation angle bin
const float *cheb;
if (g_max_el <= EL_THRESH_HIGH)
cheb = cheb7_high;
else if (g_max_el <= EL_THRESH_HIGH)
cheb = cheb7_mid;
else
cheb = cheb7_high;
float norm = chebyshev_eval(cheb, t_norm);
return norm * g_df_max; // Hz
Appendix C: Q1.16 Fixed-Point Coefficient Table
C.1 2 MHz band – 3-bin Chebyshev-7（Q1.16）
// Low angle (5–26°), CV average maximum error 215 Hz
static const int32_t cheb7_2m_low[8] = {
-869, -70013, 376, 5967, -58, -944, 10, 216
};
// Medium angle (26–57°), CV average maximum error 328 Hz
static const int32_t cheb7_2m_mid[8] = {
-204, -77140, 280, 15652, -88, -4802, 24, 2187
};
// High angle (57–90°), CV average maximum error 406 Hz
static const int32_t cheb7_2m_high[8] = {
-224, -79212, 467, 19698, -226, -7196, 159, 4131
};
C.2 70cm band – 3-bin Chebyshev-7（Q1.16）
// Low angle (5–22°), average maximum error of CV is 569 Hz
static const int32_t cheb7_70cm_low[8] = {
15, -69374, -61, 5090, 22, -797, -7, 178
};
// Medium angle (22–55°), average maximum error of CV is 752 Hz
static const int32_t cheb7_70cm_mid[8] = {
-34, -77378, 17, 16074, -6, -5089, 5, 2407
};
// High angle (55–90°), average maximum error of CV is 820 Hz
static const int32_t cheb7_70cm_high[8] = {
-7, -79773, 16, 20900, -10, -8055, 11, 5000
};
C.3 dur → Δf_max Estimation coefficient（Q1.16）
// 2m: Δf_max = c0·dur² + c1·dur + c2·dur·alt + c3·alt + c4
static const int32_t dur2df_2m[5] = {
-156, 461355, -27, -140350, 81041129
};
// 70cm: Same formula
static const int32_t dur2df_70cm[5] = {
75, 2339053, -2346, 120613, -64101833
};
C.4 Angle of Elevation Threshold
| Band | θlow | θ_high |
|---|
| 2m | ≤ 26.1° | > 56.7° |
| 70cm | ≤ 22.1° | > 54.7° |
Data: 1420 passes, 10 satellites, 30-day training window, with one satellite for cross-validation.
TLE Source: Celestrak Amateur Radio Group（2026-08-03）。
General experiment: 90 days (until October 2026) + 1 year (until September 2027). The maximum duration of TLE propagation is 395 days.
Cross-location experiment: 6 latitudes (1°N – 61°N, 33°S – 34°S), 14-day testing window.
P.S. Recently, we've been working on the firmware for the Baofeng UV-K61 (yes, the UV-K61!), and it has all the features you would expect. Stay tuned!
P.S. Why doesn't HamCQ support LaTeX?
闽ICP备2021006864号-7
闽公网安备35020602002794号