diff --git a/utils/matlab/libs/f_error_sim.m b/utils/matlab/libs/f_error_sim.m new file mode 100644 index 000000000..a0f6881f7 --- /dev/null +++ b/utils/matlab/libs/f_error_sim.m @@ -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 diff --git a/utils/matlab/libs/find_f_error_scans.m b/utils/matlab/libs/find_f_error_scans.m new file mode 100644 index 000000000..c0b6cc80e --- /dev/null +++ b/utils/matlab/libs/find_f_error_scans.m @@ -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 diff --git a/utils/matlab/libs/load_f_error_dump.m b/utils/matlab/libs/load_f_error_dump.m new file mode 100644 index 000000000..268fe1c65 --- /dev/null +++ b/utils/matlab/libs/load_f_error_dump.m @@ -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 diff --git a/utils/matlab/libs/plot_f_error_scan.m b/utils/matlab/libs/plot_f_error_scan.m new file mode 100644 index 000000000..72f058a2f --- /dev/null +++ b/utils/matlab/libs/plot_f_error_scan.m @@ -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 diff --git a/utils/matlab/plot_all_f_error_scans.m b/utils/matlab/plot_all_f_error_scans.m new file mode 100644 index 000000000..ed59fec54 --- /dev/null +++ b/utils/matlab/plot_all_f_error_scans.m @@ -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