Add Octave tooling to load, plot and simulate f_error dump scans

load_f_error_dump/find_f_error_scans/plot_f_error_scan (utils/matlab/libs)
read and visualize the per-bin CSV written by the tracking f_error_dump
diagnostic, and plot_all_f_error_scans (utils/matlab) batches
plot_f_error_scan over every scan in a dump, optionally saving PNGs.
f_error_sim (utils/matlab/libs) is a self-contained Monte Carlo
simulation of the same bin-scan/CN0-estimator logic, useful for sizing
f_error_step_num/f_error_accumulation/f_error_doppler_step without a
live capture.
This commit is contained in:
joebre
2026-09-30 07:42:39 +02:00
parent cd4f643a20
commit 5c24058459
5 changed files with 446 additions and 0 deletions
+212
View File
@@ -0,0 +1,212 @@
% f_error_sim.m
%
% Octave Monte Carlo simulation of the tracking "frequency error reduction"
% state (state 5) implemented in dll_pll_veml_tracking.cc: right after
% pull-in, a passive scan tests f_error_step_num Doppler bins centered on the
% pull-in Doppler estimate, spaced f_error_doppler_step Hz apart
% (bin 0 -> +0, then alternating +step, -step, +2*step, -2*step, ...,
% matching dll_pll_veml_tracking::f_error_bin_multiplier()). Each bin is
% incoherently accumulated for f_error_accumulation code periods (|Prompt|^2
% summed), and the bin with the highest accumulated power is kept as the
% corrected Doppler.
%
% This script reproduces that exact bin geometry and accumulation logic on
% simulated noisy correlator samples, for a given C/N0 and initial (pull-in)
% frequency error, fixed to the Galileo E1 code period (Td = 4 ms). It also
% reproduces the CN0 M2M4 estimator (lock_detectors.cc::cn0_m2m4_estimator)
% applied to the winning bin's raw samples, matching the receiver's own
% "Frequency error reduction: CN0 ..." diagnostic line.
%
% Usage:
% f_error_sim() % all defaults
% f_error_sim(CN0_dBHz, f_error_step_num, f_error_accumulation, ...
% f_error_doppler_step, f_error_init_Hz, n_trials, seed)
%
% f_error_sim(35, 9, 20, 62.5, 150) % typical EVK1029 fer.conf-like case
% f_error_sim(30, 5, 10, 125, -300, 5000, 2)
%
% What it prints:
% - the bin table (index, multiplier, Doppler tested, noiseless sinc loss)
% - Monte Carlo statistics: residual Doppler error after the scan, bin
% "capture" probability (did it pick the bin closest to the true
% frequency?), and CN0 estimate bias/spread on the winning bin
%
% What it plots:
% 1) One example trial's per-bin accumulated power vs. Doppler tested,
% with the true frequency marked -- the same data the receiver now logs
% in its "bin Doppler [Hz]" / "bin correlation power" diagnostic vectors.
% 2) Histogram of the residual Doppler error left after the scan, over all
% Monte Carlo trials.
% -------------------------------------------------------------------------
%
% GNSS-SDR is a Global Navigation Satellite System software-defined receiver.
% This file is part of GNSS-SDR.
%
% SPDX-FileCopyrightText: Joel Brenner, 2026. joel.brenner(at)saphyrion.ch
% SPDX-License-Identifier: GPL-3.0-or-later
%
% -------------------------------------------------------------------------
function f_error_sim(CN0_dBHz, f_error_step_num, f_error_accumulation, f_error_doppler_step, f_error_init_Hz, n_trials, seed)
%% ------------------------- PARAMETERS ----------------------------------
if nargin < 1 || isempty(CN0_dBHz), CN0_dBHz = 35; end % [dB-Hz]
if nargin < 2 || isempty(f_error_step_num), f_error_step_num = 9; end % number of Doppler bins
if nargin < 3 || isempty(f_error_accumulation), f_error_accumulation = 20; end % incoherent accum. per bin [code periods]
if nargin < 4 || isempty(f_error_doppler_step), f_error_doppler_step = 62.5; end % [Hz] bin spacing
if nargin < 5 || isempty(f_error_init_Hz), f_error_init_Hz = 150; end % [Hz] true residual Doppler at pull-in
if nargin < 6 || isempty(n_trials), n_trials = 2000; end % Monte Carlo trial count
if nargin < 7 || isempty(seed), seed = 1; end % RNG seed, [] = random
%% ------------------------------------------------------------------------
if ~isempty(seed)
rand('state', seed);
randn('state', seed);
end
Td = 4e-3; % Galileo E1 code period [s] -- fixed per this script's scope
% mirror Dll_Pll_Conf::SetFromConfiguration(): f_error_step_num must be odd
% (a center bin plus a symmetric number of +/- steps); round up if even.
f_error_step_num = max(1, round(f_error_step_num));
if mod(f_error_step_num, 2) == 0
f_error_step_num = f_error_step_num + 1;
printf('note: f_error_step_num rounded up to %d (must be odd)\n', f_error_step_num);
end
num_bins = f_error_step_num;
% bin 0 -> 0; then alternating outward: +1, -1, +2, -2, +3, -3, ...
% (matches dll_pll_veml_tracking::f_error_bin_multiplier(), 0-based bin_index)
bin_mult = zeros(1, num_bins);
for idx = 1:num_bins
bin_index = idx - 1;
if bin_index == 0
bin_mult(idx) = 0;
else
half_steps = floor((bin_index + 1) / 2);
if mod(bin_index, 2) == 1
bin_mult(idx) = half_steps;
else
bin_mult(idx) = -half_steps;
end
end
end
bin_doppler_offset = bin_mult * f_error_doppler_step; % [Hz], relative to pull-in estimate
bin_residual_error = f_error_init_Hz - bin_doppler_offset; % [Hz], true - tested, per bin
sinc0 = @(x) (x == 0) * 1 + (x ~= 0) .* (sin(x + (x == 0)) ./ (x + (x == 0)));
A = sqrt(2 * 10^(CN0_dBHz / 10) * Td); % correlator amplitude; SNR = 2*(C/N0)*Td, unit-variance complex noise per branch
printf('--- f_error_sim: CN0=%.1f dB-Hz, f_error_step_num=%d, f_error_accumulation=%d, f_error_doppler_step=%.2f Hz, f_error_init=%.1f Hz, Td=%.1f ms (Galileo E1) ---\n', ...
CN0_dBHz, num_bins, f_error_accumulation, f_error_doppler_step, f_error_init_Hz, Td * 1e3);
% noiseless per-bin loss, for the printed table
theta_bins = pi * bin_residual_error * Td;
bin_loss_dB = 20 * log10(abs(sinc0(theta_bins)) + eps);
printf('\n bin mult Doppler[Hz] residual[Hz] sinc loss[dB]\n');
[~, ideal_bin] = min(abs(bin_residual_error));
for idx = 1:num_bins
marker = '';
if idx == ideal_bin, marker = ' <- closest to truth'; end
printf(' %2d %+3d %8.2f %8.2f %6.2f%s\n', ...
idx - 1, bin_mult(idx), bin_doppler_offset(idx), bin_residual_error(idx), bin_loss_dB(idx), marker);
end
%% --------------------------- Monte Carlo ---------------------------------
winning_bin = zeros(1, n_trials);
residual_after = zeros(1, n_trials);
cn0_est_dBHz = zeros(1, n_trials);
example_power = []; % saved from trial 1, for plotting
for t = 1:n_trials
bin_power = zeros(1, num_bins);
bin_samples = cell(1, num_bins);
for b = 1:num_bins
theta = pi * bin_residual_error(b) * Td;
amp = A * sinc0(theta);
k = (1:f_error_accumulation).';
phase = (2 * k - 1) * theta; % midpoint-sampling phase ramp, one full 2*theta per epoch
samples = amp * exp(1j * phase) + (randn(f_error_accumulation, 1) + 1j * randn(f_error_accumulation, 1));
bin_power(b) = sum(abs(samples) .^ 2);
bin_samples{b} = samples;
end
[~, best_bin] = max(bin_power);
winning_bin(t) = best_bin;
residual_after(t) = f_error_init_Hz - bin_doppler_offset(best_bin);
cn0_est_dBHz(t) = cn0_m2m4(bin_samples{best_bin}, Td);
if t == 1
example_power = bin_power;
end
end
capture_prob = mean(winning_bin == ideal_bin);
printf('\nMonte Carlo results (%d trials):\n', n_trials);
printf(' Residual Doppler error after scan: mean=%.2f Hz, std=%.2f Hz, RMSE=%.2f Hz (initial error was %.1f Hz)\n', ...
mean(residual_after), std(residual_after), sqrt(mean(residual_after .^ 2)), f_error_init_Hz);
printf(' Bin-capture probability (picked bin closest to truth, bin %d): %.1f%%\n', ideal_bin - 1, 100 * capture_prob);
printf(' CN0 estimate on winning bin: mean=%.2f dB-Hz, std=%.2f dB-Hz (true CN0=%.1f dB-Hz)\n', ...
mean(cn0_est_dBHz), std(cn0_est_dBHz), CN0_dBHz);
%% ------------------------------ Plots -------------------------------------
figure('Name', 'Frequency error reduction: example bin scan');
bar(bin_doppler_offset, example_power);
hold on;
yl = ylim();
plot([f_error_init_Hz, f_error_init_Hz], yl, 'r--', 'LineWidth', 1.5);
hold off;
xlabel('Doppler tested [Hz]');
ylabel('Accumulated |Prompt|^2 power');
title(sprintf('Example bin scan (trial 1) -- true error = %.1f Hz', f_error_init_Hz));
legend('bin power', 'true frequency', 'Location', 'northeast');
grid on;
figure('Name', 'Frequency error reduction: residual error distribution');
hist(residual_after, min(40, max(10, round(n_trials / 50))));
xlabel('Residual Doppler error after scan [Hz]');
ylabel('Trial count');
title(sprintf('Residual error after %d-bin scan (initial error = %.1f Hz)', num_bins, f_error_init_Hz));
grid on;
end
function cn0_dBHz = cn0_m2m4(prompt_samples, coh_integration_time_s)
% Octave port of lock_detectors.cc::cn0_m2m4_estimator(), applied to a
% column vector of complex Prompt samples, for exact parity with the
% receiver's own state-5 CN0 diagnostic.
n = numel(prompt_samples);
if n == 0 || coh_integration_time_s == 0
cn0_dBHz = -100;
return;
end
Psig = mean(abs(real(prompt_samples)));
Psig = Psig^2;
p2 = abs(prompt_samples) .^ 2;
m2 = mean(p2);
m4 = mean(p2 .^ 2);
aux = sqrt(2 * m2^2 - m4);
if isnan(aux)
denominator = m2 - Psig;
if denominator == 0
cn0_dBHz = -100;
return;
end
SNR_aux = Psig / denominator;
else
denominator = m2 - aux;
if denominator == 0
cn0_dBHz = -100;
return;
end
SNR_aux = aux / denominator;
end
if SNR_aux <= 0
cn0_dBHz = -100;
return;
end
cn0_dBHz = 10 * log10(SNR_aux) - 10 * log10(coh_integration_time_s);
end
+29
View File
@@ -0,0 +1,29 @@
% find_f_error_scans.m
%
% Returns the indices into a scans struct array (as returned by
% load_f_error_dump) matching a given PRN, optionally restricted to a
% specific channel.
%
% Usage:
% idxs = find_f_error_scans(scans, 13); % all scans for PRN 13, any channel
% idxs = find_f_error_scans(scans, 13, 2); % PRN 13 on channel 2 only
%
% See also: load_f_error_dump, plot_f_error_scan
% -------------------------------------------------------------------------
%
% GNSS-SDR is a Global Navigation Satellite System software-defined receiver.
% This file is part of GNSS-SDR.
%
% SPDX-FileCopyrightText: Joel Brenner, 2026. joel.brenner(at)saphyrion.ch
% SPDX-License-Identifier: GPL-3.0-or-later
%
% -------------------------------------------------------------------------
function idxs = find_f_error_scans(scans, prn, channel)
mask = [scans.prn] == prn;
if nargin >= 3 && ~isempty(channel)
mask = mask & ([scans.channel] == channel);
end
idxs = find(mask);
end
+78
View File
@@ -0,0 +1,78 @@
% load_f_error_dump.m
%
% Loads the CSV written by dll_pll_veml_tracking's state-5 (frequency error
% reduction) diagnostic -- see Dll_Pll_Conf::f_error_dump_filename -- into an
% Octave struct array, one element per scan (one satellite's pass through
% state 5). Each row of the CSV is one Doppler bin; this groups rows by
% scan_id.
%
% CSV columns (header row is skipped):
% scan_id, prn, system_char, channel, cn0_dBHz, selected_doppler_hz,
% selected_bin, bin_index, doppler_hz, power
%
% Usage:
% scans = load_f_error_dump(); % default: ./f_error_dump.csv
% scans = load_f_error_dump('/path/to/dump.csv');
% load_f_error_dump(); % no output arg: also prints a summary table
%
% Each scans(i) has fields:
% scan_id, prn, system_char, channel, cn0_dBHz, selected_doppler_hz, selected_bin,
% doppler_hz (row vector, increasing Doppler order), power (row vector, same order)
%
% See also: find_f_error_scans, plot_f_error_scan
% -------------------------------------------------------------------------
%
% GNSS-SDR is a Global Navigation Satellite System software-defined receiver.
% This file is part of GNSS-SDR.
%
% SPDX-FileCopyrightText: Joel Brenner, 2026. joel.brenner(at)saphyrion.ch
% SPDX-License-Identifier: GPL-3.0-or-later
%
% -------------------------------------------------------------------------
function scans = load_f_error_dump(filename)
if nargin < 1 || isempty(filename)
filename = 'f_error_dump.csv';
end
raw = csvread(filename, 1, 0); % skip header row
if isempty(raw)
error('load_f_error_dump: no data rows found in %s', filename);
end
scan_id_col = raw(:, 1);
unique_ids = unique(scan_id_col);
n = numel(unique_ids);
scans = struct('scan_id', {}, 'prn', {}, 'system_char', {}, 'channel', {}, 'cn0_dBHz', {}, ...
'selected_doppler_hz', {}, 'selected_bin', {}, ...
'doppler_hz', {}, 'power', {});
for i = 1:n
rows = raw(scan_id_col == unique_ids(i), :);
[~, order] = sort(rows(:, 8)); % bin_index, defensive (writer already emits increasing Doppler order)
rows = rows(order, :);
scans(i).scan_id = rows(1, 1);
scans(i).prn = rows(1, 2);
scans(i).system_char = rows(1, 3);
scans(i).channel = rows(1, 4);
scans(i).cn0_dBHz = rows(1, 5);
scans(i).selected_doppler_hz = rows(1, 6);
scans(i).selected_bin = rows(1, 7);
scans(i).doppler_hz = rows(:, 9)';
scans(i).power = rows(:, 10)';
end
if nargout == 0
printf('%6s %5s %4s %9s %14s %8s %6s\n', 'scan', 'PRN', 'ch', 'CN0[dBHz]', 'selDoppler[Hz]', 'nBins', '');
for i = 1:n
printf('%6d %5d %4d %9.2f %14.2f %8d\n', scans(i).scan_id, scans(i).prn, scans(i).channel, ...
scans(i).cn0_dBHz, scans(i).selected_doppler_hz, numel(scans(i).doppler_hz));
end
printf('\n%d scans loaded from %s\n', n, filename);
printf('PRNs present: %s\n', mat2str(unique([scans.prn])));
clear scans; % avoid echoing the struct array to the console when called as a script
end
end
+59
View File
@@ -0,0 +1,59 @@
% plot_f_error_scan.m
%
% Bar-plots one state-5 frequency-error-reduction scan (as loaded by
% load_f_error_dump), with the selected/winning Doppler marked. Optionally
% overlays the theoretical sinc(pi*Delta_f*Td)^2 envelope (normalized to the
% scan's own peak) for comparison, if the coherent integration time Td [s] is
% given -- e.g. Td=0.004 for Galileo E1, Td=0.001 for GPS L1 C/A.
%
% Usage:
% scans = load_f_error_dump();
% plot_f_error_scan(scans, 5); % plot scans(5)
% plot_f_error_scan(scans, 5, 0.004); % + theoretical Galileo E1 sinc^2 overlay
%
% See also: load_f_error_dump, find_f_error_scans
% -------------------------------------------------------------------------
%
% GNSS-SDR is a Global Navigation Satellite System software-defined receiver.
% This file is part of GNSS-SDR.
%
% SPDX-FileCopyrightText: Joel Brenner, 2026. joel.brenner(at)saphyrion.ch
% SPDX-License-Identifier: GPL-3.0-or-later
%
% -------------------------------------------------------------------------
function plot_f_error_scan(scans, idx, Td)
if idx < 1 || idx > numel(scans)
error('plot_f_error_scan: idx %d out of range (1..%d)', idx, numel(scans));
end
s = scans(idx);
figure();
bar(s.doppler_hz, s.power);
hold on;
if nargin >= 3 && ~isempty(Td)
% theoretical envelope centered on the selected (winning) Doppler, since that's
% this scan's best estimate of the true residual frequency
theta = pi * (s.doppler_hz - s.selected_doppler_hz) * Td;
sinc0 = @(x) (x == 0) * 1 + (x ~= 0) .* (sin(x + (x == 0)) ./ (x + (x == 0)));
envelope = (sinc0(theta) .^ 2) * max(s.power);
plot(s.doppler_hz, envelope, 'k--', 'LineWidth', 1.3);
end
yl = ylim();
plot([s.selected_doppler_hz, s.selected_doppler_hz], yl, 'r--', 'LineWidth', 1.5);
hold off;
xlabel('Doppler tested [Hz]');
ylabel('Accumulated correlation power');
if nargin >= 3 && ~isempty(Td)
legend('bin power', 'theoretical sinc^2 (centered on selected)', 'selected Doppler', 'Location', 'northeast');
else
legend('bin power', 'selected Doppler', 'Location', 'northeast');
end
title(sprintf('PRN%d ch%d scan %d: CN0=%.1f dB-Hz, selected=%.1f Hz (bin %d/%d)', ...
s.prn, s.channel, s.scan_id, s.cn0_dBHz, s.selected_doppler_hz, s.selected_bin, numel(s.doppler_hz) - 1));
grid on;
end
+68
View File
@@ -0,0 +1,68 @@
% plot_all_f_error_scans.m
%
% Bar-plots every scan in a scans struct array (as returned by
% load_f_error_dump), one figure per scan -- a batch wrapper around
% plot_f_error_scan. Optionally restrict to a subset of scans (e.g. from
% find_f_error_scans) and/or save each figure to disk instead of leaving a
% window open per scan (recommended once there are more than a handful).
%
% Usage:
% scans = load_f_error_dump();
% plot_all_f_error_scans(scans); % every scan, one window each
% plot_all_f_error_scans(scans, 0.004); % + theoretical sinc^2 overlay (Galileo E1 Td)
% plot_all_f_error_scans(scans, 0.004, find_f_error_scans(scans, 13)); % only PRN 13's scans
% plot_all_f_error_scans(scans, 0.004, [], 'plots'); % save PNGs into ./plots/ instead of opening windows
%
% See also: load_f_error_dump, find_f_error_scans, plot_f_error_scan
% -------------------------------------------------------------------------
%
% GNSS-SDR is a Global Navigation Satellite System software-defined receiver.
% This file is part of GNSS-SDR.
%
% SPDX-FileCopyrightText: Joel Brenner, 2026. joel.brenner(at)saphyrion.ch
% SPDX-License-Identifier: GPL-3.0-or-later
%
% -------------------------------------------------------------------------
function plot_all_f_error_scans(scans, Td, idxs, save_dir)
if ~exist('load_f_error_dump.m', 'file')
addpath('./libs')
end
if nargin < 2
Td = [];
end
if nargin < 3 || isempty(idxs)
idxs = 1:numel(scans);
end
if nargin < 4
save_dir = '';
end
if ~isempty(save_dir)
% note: do NOT set DefaultFigureVisible='off' here -- combined with an
% invisible figure, saveas()/print() falls back to the fltk renderer,
% which then fails without a real X display. Leaving figures nominally
% visible (but headlessly rendered via gnuplot) is what actually works.
graphics_toolkit('gnuplot');
if ~exist(save_dir, 'dir')
mkdir(save_dir);
end
end
for idx = idxs
s = scans(idx);
plot_f_error_scan(scans, idx, Td);
if ~isempty(save_dir)
fname = fullfile(save_dir, sprintf('scan%04d_prn%d_ch%d.png', s.scan_id, s.prn, s.channel));
saveas(gcf(), fname);
close(gcf());
end
end
if ~isempty(save_dir)
set(0, 'DefaultFigureVisible', 'on');
printf('Saved %d plots to %s\n', numel(idxs), save_dir);
end
end