diff --git a/Classes/00_signals/Signal.m b/Classes/00_signals/Signal.m index 613b7a0..229bce9 100644 --- a/Classes/00_signals/Signal.m +++ b/Classes/00_signals/Signal.m @@ -24,9 +24,9 @@ classdef Signal obj.signal = signal; obj.signal = obj.signal; obj.fs = options.fs; - - [~,obj.gitSHA] = system('git rev-parse HEAD'); - [~,obj.gitStatus] = system('git status --porcelain'); + % + % [~,obj.gitSHA] = system('git rev-parse HEAD'); + % [~,obj.gitStatus] = system('git status --porcelain'); % [~,obj.gitPatch] = system('git diff'); %%% Stuff for Logbook %%% diff --git a/Classes/04_DSP/Equalizer/FFE.m b/Classes/04_DSP/Equalizer/FFE.m index e624137..0ae09ff 100644 --- a/Classes/04_DSP/Equalizer/FFE.m +++ b/Classes/04_DSP/Equalizer/FFE.m @@ -35,16 +35,19 @@ classdef FFE < handle % A1 moving-average input suppression dc_avg_bufferlength_a1 + dc_avg_update_blocklength_a1 dc_smoothing_a1 % A2 level-dependent residual suppression dc_level_avg_bufferlength_a2 + dc_level_update_blocklength_a2 dc_smoothing_a2 dc_level_weights_a2 % Adaptive DC-tracking loop dc_tracking_mu dc_tracking_adaptive_enabled + dc_tracking_persistence_gain dc_tracking_power_exponent dc_tracking_buffer_len @@ -110,14 +113,17 @@ classdef FFE < handle options.dd_len_fraction = 1; options.dc_avg_bufferlength_a1 = 0; + options.dc_avg_update_blocklength_a1 = 0; options.dc_smoothing_a1 = 0; options.dc_level_avg_bufferlength_a2 = 0; + options.dc_level_update_blocklength_a2 = 0; options.dc_smoothing_a2 = 0; options.dc_level_weights_a2 = 0; options.dc_tracking_mu = 0; options.dc_tracking_adaptive_enabled = false; + options.dc_tracking_persistence_gain = 0; options.dc_tracking_power_exponent = 2; options.dc_tracking_buffer_len = 1; @@ -150,15 +156,19 @@ classdef FFE < handle assert(obj.dc_tracking_buffer_len >= 0); assert(obj.ffe_update_buffer_len >= 0); assert(obj.dc_avg_bufferlength_a1 >= 0); + assert(obj.dc_avg_update_blocklength_a1 >= 0); assert(obj.dc_level_avg_bufferlength_a2 >= 0); + assert(obj.dc_level_update_blocklength_a2 >= 0); obj.e = zeros(obj.order,1); obj.e_dc = 0; obj.error = 0; obj.a2_level_weight_initial_stats = struct(); obj.dc_avg_bufferlength_a1 = floor(obj.dc_avg_bufferlength_a1); + obj.dc_avg_update_blocklength_a1 = floor(obj.dc_avg_update_blocklength_a1); obj.dc_smoothing_a1 = min(max(obj.dc_smoothing_a1,0),1); obj.dc_level_avg_bufferlength_a2 = floor(obj.dc_level_avg_bufferlength_a2); + obj.dc_level_update_blocklength_a2 = floor(obj.dc_level_update_blocklength_a2); obj.dc_smoothing_a2 = min(max(obj.dc_smoothing_a2,0),1); obj.dc_tracking_buffer_len = floor(obj.dc_tracking_buffer_len); obj.ffe_update_buffer_len = floor(obj.ffe_update_buffer_len); @@ -260,6 +270,24 @@ classdef FFE < handle x = [zeros(floor(obj.order/2),1); x; zeros(obj.order,1)]; lambda = mu; + n_symbols = ceil(N / obj.sps); + y = zeros(n_symbols,1); + d_hat = zeros(n_symbols,1); + err = zeros(n_symbols,1); + true_err = zeros(n_symbols,1); + constellation = obj.constellation; + + if obj.adaption_technique == adaption_method.nlms + adaption_code = 1; + elseif obj.adaption_technique == adaption_method.lms + adaption_code = 2; + elseif obj.adaption_technique == adaption_method.rls + adaption_code = 3; + else + builtin("error","FFE:InvalidAdaptionTechnique", ... + "Unsupported FFE adaption technique."); + end + adaption_is_rls = adaption_code == 3; if training mask = ones(obj.order,1); @@ -280,21 +308,31 @@ classdef FFE < handle epochs = 1; end - P_err = 0; - err_prev = 0; - dc_tracking_alpha = obj.dc_tracking_alpha; - dc_tracking_gamma = obj.dc_tracking_gamma; dc_tracking_mu_min = obj.dc_tracking_mu_min; dc_tracking_mu_max = obj.dc_tracking_mu_max; dc_tracking_mu_eff_min = obj.dc_tracking_mu_eff_min; dc_tracking_mu_eff_max = obj.dc_tracking_mu_eff_max; - dc_tracking_power_exponent = obj.dc_tracking_power_exponent; - dc_buffer_enabled = obj.dc_tracking_mu ~= 0 && obj.dc_tracking_buffer_len > 1; + dc_tracking_persistence_gain = 0; + if obj.dc_tracking_adaptive_enabled + dc_tracking_persistence_gain = obj.dc_tracking_persistence_gain; + end + dc_tracking_enabled = obj.dc_tracking_mu ~= 0; + if dc_tracking_enabled + obj.dc_tracking_mu = min(max(obj.dc_tracking_mu,dc_tracking_mu_min),dc_tracking_mu_max); + end + dc_tracking_mu = obj.dc_tracking_mu; + dc_tracking_use_persistence = dc_tracking_persistence_gain > 0; + dc_tracking_base_mu_eff = min(max(dc_tracking_mu,dc_tracking_mu_eff_min),dc_tracking_mu_eff_max); + dc_buffer_enabled = dc_tracking_enabled && obj.dc_tracking_buffer_len > 1; if dc_buffer_enabled - e_dc_buffer = NaN(obj.dc_tracking_buffer_len,1); + dc_tracking_err_buffer = NaN(obj.dc_tracking_buffer_len,1); + dc_tracking_err_buffer_pos = 0; + dc_tracking_err_sum = 0; + dc_tracking_abs_err_sum = 0; + dc_tracking_valid_count = 0; end - ffe_buffer_enabled = obj.ffe_update_buffer_len > 1 && obj.adaption_technique ~= adaption_method.rls; + ffe_buffer_enabled = obj.ffe_update_buffer_len > 1 && ~adaption_is_rls; if ffe_buffer_enabled ffe_update_buffer = NaN(obj.order,obj.ffe_update_buffer_len); end @@ -303,7 +341,10 @@ classdef FFE < handle if dc_avg_enabled dc_avg_buffer_a1 = NaN(obj.dc_avg_bufferlength_a1,1); dc_avg_est_a1 = 0; - dc_avg_update_blocklength_a1 = obj.dc_avg_bufferlength_a1; + dc_avg_update_blocklength_a1 = obj.dc_avg_update_blocklength_a1; + if dc_avg_update_blocklength_a1 <= 0 + dc_avg_update_blocklength_a1 = obj.dc_avg_bufferlength_a1; + end dc_avg_update_blocklength_a1 = max(1,floor(dc_avg_update_blocklength_a1)); dc_avg_input_offset_a1 = floor(obj.order/2); % Hardware-like A1 is causal; dc_smoothing_a1 is reserved for offline variants. @@ -311,12 +352,12 @@ classdef FFE < handle dc_level_enabled = obj.dc_level_avg_bufferlength_a2 > 1 && any(obj.dc_level_weights_a2(:) ~= 0); if dc_level_enabled - if isempty(obj.constellation) + if isempty(constellation) builtin("error","FFE:MissingConstellation", ... "A2 level-dependent MPI suppression requires obj.constellation to be set."); end - n_levels = numel(obj.constellation); + n_levels = numel(constellation); if isscalar(obj.dc_level_weights_a2) dc_level_weight_by_level = repmat(obj.dc_level_weights_a2,n_levels,1); elseif numel(obj.dc_level_weights_a2) == n_levels @@ -325,17 +366,24 @@ classdef FFE < handle builtin("error","FFE:InvalidDCLevelWeights", ... "dc_level_weights_a2 must be scalar or have one entry per constellation level."); end - dc_level_err_buffer = NaN(n_levels,obj.dc_level_avg_bufferlength_a2); + dc_level_buffer_len_a2 = obj.dc_level_avg_bufferlength_a2; + dc_level_err_buffer = NaN(n_levels,dc_level_buffer_len_a2); + dc_level_err_buffer_pos_by_level = zeros(n_levels,1); + dc_level_err_sum_by_level = zeros(n_levels,1); + dc_level_buffer_valid_count_by_level = zeros(n_levels,1); dc_level_mpi_est_by_level = zeros(n_levels,1); dc_level_valid_count_by_level = zeros(n_levels,1); - dc_level_update_blocklength_a2 = obj.dc_level_avg_bufferlength_a2; + dc_level_update_blocklength_a2 = obj.dc_level_update_blocklength_a2; + if dc_level_update_blocklength_a2 <= 0 + dc_level_update_blocklength_a2 = dc_level_buffer_len_a2; + end dc_level_update_blocklength_a2 = max(1,floor(dc_level_update_blocklength_a2)); dc_level_window_future_fraction = obj.dc_smoothing_a2; %#ok % reserved for delayed/offline A2 variants end debug_enabled = obj.save_debug; if debug_enabled - n_symbols_debug = ceil(N / obj.sps); + n_symbols_debug = n_symbols; obj.debug_struct.error = NaN(1,n_symbols_debug); obj.debug_struct.error_first_epoch = NaN(1,n_symbols_debug); obj.debug_struct.main_cursor = NaN(1,n_symbols_debug); @@ -392,15 +440,15 @@ classdef FFE < handle if training d_hat(symbol,1) = d(symbol); - if isempty(obj.constellation) + if isempty(constellation) symbol_idx = NaN; else - [~,symbol_idx] = min(abs(d_hat(symbol) - obj.constellation)); + [~,symbol_idx] = min(abs(d_hat(symbol) - constellation)); end else if ~always_ideal_decision - [~,symbol_idx] = min(abs(y(symbol) - obj.constellation)); % decision for closest constellation point - d_hat(symbol,1) = obj.constellation(symbol_idx); + [~,symbol_idx] = min(abs(y(symbol) - constellation)); % decision for closest constellation point + d_hat(symbol,1) = constellation(symbol_idx); else d_hat(symbol,1) = d(symbol); end @@ -412,26 +460,63 @@ classdef FFE < handle dc_level_mpi_est = dc_level_mpi_est_by_level(symbol_idx); dc_level_valid_count = dc_level_valid_count_by_level(symbol_idx); dc_level_weight = dc_level_weight_by_level(symbol_idx) * ... - min(dc_level_valid_count / obj.dc_level_avg_bufferlength_a2,1); + min(dc_level_valid_count / dc_level_buffer_len_a2,1); y(symbol,1) = y_raw - dc_level_weight * dc_level_mpi_est; if training d_hat(symbol,1) = d(symbol); else if ~always_ideal_decision - [~,symbol_idx] = min(abs(y(symbol) - obj.constellation)); - d_hat(symbol,1) = obj.constellation(symbol_idx); + [~,symbol_idx] = min(abs(y(symbol) - constellation)); + d_hat(symbol,1) = constellation(symbol_idx); else d_hat(symbol,1) = d(symbol); end end - dc_level_err_buffer(dc_level_symbol_idx,:) = circshift(dc_level_err_buffer(dc_level_symbol_idx,:),1,2); - dc_level_err_buffer(dc_level_symbol_idx,1) = mpi_err; + dc_level_err_buffer_pos = dc_level_err_buffer_pos_by_level(dc_level_symbol_idx) + 1; + if dc_level_err_buffer_pos > dc_level_buffer_len_a2 + dc_level_err_buffer_pos = 1; + end + dc_level_err_buffer_pos_by_level(dc_level_symbol_idx) = dc_level_err_buffer_pos; + + dc_level_old_err = dc_level_err_buffer(dc_level_symbol_idx,dc_level_err_buffer_pos); + if isfinite(dc_level_old_err) + dc_level_err_sum_by_level(dc_level_symbol_idx) = ... + dc_level_err_sum_by_level(dc_level_symbol_idx) - dc_level_old_err; + dc_level_buffer_valid_count_by_level(dc_level_symbol_idx) = ... + dc_level_buffer_valid_count_by_level(dc_level_symbol_idx) - 1; + end + + if isfinite(mpi_err) + dc_level_err_buffer(dc_level_symbol_idx,dc_level_err_buffer_pos) = mpi_err; + dc_level_err_sum_by_level(dc_level_symbol_idx) = ... + dc_level_err_sum_by_level(dc_level_symbol_idx) + mpi_err; + dc_level_buffer_valid_count_by_level(dc_level_symbol_idx) = ... + dc_level_buffer_valid_count_by_level(dc_level_symbol_idx) + 1; + else + dc_level_err_buffer(dc_level_symbol_idx,dc_level_err_buffer_pos) = NaN; + end + if mod(symbol,dc_level_update_blocklength_a2) == 0 - dc_level_valid_count_by_level = sum(isfinite(dc_level_err_buffer),2); - dc_level_mpi_est_by_level = mean(dc_level_err_buffer,2,"omitnan"); - dc_level_mpi_est_by_level(dc_level_valid_count_by_level == 0) = 0; + if dc_level_update_blocklength_a2 == 1 + dc_level_valid_count_by_level(dc_level_symbol_idx) = ... + dc_level_buffer_valid_count_by_level(dc_level_symbol_idx); + if dc_level_valid_count_by_level(dc_level_symbol_idx) == 0 + dc_level_mpi_est_by_level(dc_level_symbol_idx) = 0; + else + dc_level_mpi_est_by_level(dc_level_symbol_idx) = ... + dc_level_err_sum_by_level(dc_level_symbol_idx) / ... + dc_level_valid_count_by_level(dc_level_symbol_idx); + end + else + dc_level_valid_count_by_level = dc_level_buffer_valid_count_by_level; + dc_level_has_valid = dc_level_valid_count_by_level > 0; + dc_level_mpi_est_by_level(:) = 0; + dc_level_mpi_est_by_level(dc_level_has_valid) = ... + dc_level_err_sum_by_level(dc_level_has_valid) ./ ... + dc_level_valid_count_by_level(dc_level_has_valid); + end end end @@ -440,10 +525,10 @@ classdef FFE < handle true_err(symbol) = y(symbol) - d(symbol); % Instantaneous error - if training || obj.dd_mode - switch obj.adaption_technique + if 1 %training || obj.dd_mode + switch adaption_code - case adaption_method.nlms + case 1 % mu used as update weight (suggestion: 0.01-0.05; bit higher during tr) normU = ((U.'*U)) + eps; @@ -452,7 +537,7 @@ classdef FFE < handle update = grad * weight; - case adaption_method.lms + case 2 % mu used as update weight (suggestion: 0.001) weight = mu; @@ -461,7 +546,7 @@ classdef FFE < handle - case adaption_method.rls + case 3 % RLS‐Gain: @@ -486,30 +571,87 @@ classdef FFE < handle obj.e = obj.e + update; end - if obj.adaption_technique == adaption_method.rls + if adaption_is_rls obj.P = (1/lambda) * (obj.P - k * (U.' * obj.P)); end - if obj.dc_tracking_mu ~= 0 - if obj.dc_tracking_adaptive_enabled - delta_mu = dc_tracking_gamma * err(symbol) * err_prev * (U.'*U); - obj.dc_tracking_mu = min(max(obj.dc_tracking_mu + delta_mu,dc_tracking_mu_min),dc_tracking_mu_max); - err_prev = err(symbol); - P_err = dc_tracking_alpha*P_err + (1-dc_tracking_alpha)*err(symbol)^2; - dc_tracking_mu_eff = obj.dc_tracking_mu / ((P_err + eps)^dc_tracking_power_exponent); - else - dc_tracking_mu_eff = obj.dc_tracking_mu; - end - - dc_tracking_mu_eff = min(max(dc_tracking_mu_eff,dc_tracking_mu_eff_min),dc_tracking_mu_eff_max); + if dc_tracking_enabled if dc_buffer_enabled - e_dc_buffer(1) = obj.e_dc + dc_tracking_mu_eff * err(symbol); - e_dc_buffer = circshift(e_dc_buffer,1); + dc_tracking_err_buffer_pos = dc_tracking_err_buffer_pos + 1; + if dc_tracking_err_buffer_pos > obj.dc_tracking_buffer_len + dc_tracking_err_buffer_pos = 1; + end + + dc_tracking_old_err = dc_tracking_err_buffer(dc_tracking_err_buffer_pos); + if isfinite(dc_tracking_old_err) + dc_tracking_err_sum = dc_tracking_err_sum - dc_tracking_old_err; + dc_tracking_valid_count = dc_tracking_valid_count - 1; + if dc_tracking_use_persistence + dc_tracking_abs_err_sum = dc_tracking_abs_err_sum - abs(dc_tracking_old_err); + end + end + + dc_tracking_new_err = err(symbol); + if isfinite(dc_tracking_new_err) + dc_tracking_err_buffer(dc_tracking_err_buffer_pos) = dc_tracking_new_err; + dc_tracking_err_sum = dc_tracking_err_sum + dc_tracking_new_err; + dc_tracking_valid_count = dc_tracking_valid_count + 1; + if dc_tracking_use_persistence + dc_tracking_abs_err_sum = dc_tracking_abs_err_sum + abs(dc_tracking_new_err); + end + else + dc_tracking_err_buffer(dc_tracking_err_buffer_pos) = NaN; + end + if mod(symbol,obj.dc_tracking_buffer_len) == 0 - obj.e_dc = mean(e_dc_buffer,"omitnan"); + if dc_tracking_valid_count == 0 + dc_tracking_err_mean = 0; + if dc_tracking_use_persistence + dc_tracking_err_abs_mean = 0; + end + else + dc_tracking_err_mean = dc_tracking_err_sum / dc_tracking_valid_count; + if dc_tracking_use_persistence + dc_tracking_err_abs_mean = dc_tracking_abs_err_sum / dc_tracking_valid_count; + end + end + + if dc_tracking_use_persistence + dc_tracking_persistence_scale = abs(dc_tracking_err_mean) / (dc_tracking_err_abs_mean + eps); + dc_tracking_persistence_scale = min(max(dc_tracking_persistence_scale,0),1); + dc_tracking_mu_eff = dc_tracking_mu * ... + (1 + dc_tracking_persistence_gain * dc_tracking_persistence_scale); + dc_tracking_mu_eff = min(max(dc_tracking_mu_eff,dc_tracking_mu_eff_min),dc_tracking_mu_eff_max); + else + dc_tracking_mu_eff = dc_tracking_base_mu_eff; + end + obj.e_dc = obj.e_dc + dc_tracking_mu_eff * dc_tracking_err_mean; end else - obj.e_dc = obj.e_dc + dc_tracking_mu_eff * err(symbol); + dc_tracking_err_mean = err(symbol); + if isfinite(dc_tracking_err_mean) + if dc_tracking_use_persistence + dc_tracking_err_abs_mean = abs(dc_tracking_err_mean); + dc_tracking_persistence_scale = abs(dc_tracking_err_mean) / ... + (dc_tracking_err_abs_mean + eps); + else + dc_tracking_mu_eff = dc_tracking_base_mu_eff; + end + else + dc_tracking_err_mean = 0; + if dc_tracking_use_persistence + dc_tracking_persistence_scale = 0; + else + dc_tracking_mu_eff = dc_tracking_base_mu_eff; + end + end + if dc_tracking_use_persistence + dc_tracking_persistence_scale = min(max(dc_tracking_persistence_scale,0),1); + dc_tracking_mu_eff = dc_tracking_mu * ... + (1 + dc_tracking_persistence_gain * dc_tracking_persistence_scale); + dc_tracking_mu_eff = min(max(dc_tracking_mu_eff,dc_tracking_mu_eff_min),dc_tracking_mu_eff_max); + end + obj.e_dc = obj.e_dc + dc_tracking_mu_eff * dc_tracking_err_mean; end end end @@ -521,7 +663,7 @@ classdef FFE < handle if debug_enabled && epoch == epochs error_power = err(symbol) * err(symbol)'; - update_power = update.'*update ./ (rms(obj.e) + eps); + update_power = update.'*update ./ (sqrt((obj.e.'*obj.e) / obj.order) + eps); obj.debug_struct.error(1,symbol) = error_power; obj.debug_struct.main_cursor(1,symbol) = abs(obj.e(maincursor_pos)); obj.debug_struct.mu_nlms(1,symbol) = weight; @@ -611,9 +753,7 @@ classdef FFE < handle vars = optimizableVariable("dc_tracking_mu",[1e-5,1e-1],"Transform","log"); if obj.dc_tracking_adaptive_enabled vars = [vars, ... - optimizableVariable("dc_tracking_alpha",[0.85,0.995]), ... - optimizableVariable("dc_tracking_gamma",[1e-7,3e-5],"Transform","log"), ... - optimizableVariable("dc_tracking_power_exponent",[0,2]), ... + optimizableVariable("dc_tracking_persistence_gain",[0,2]), ... optimizableVariable("dc_tracking_mu_eff_max",[1e-3,3e-1],"Transform","log")]; end @@ -630,12 +770,10 @@ classdef FFE < handle best = obj.dc_tracking_optimization.XAtMinObjective; obj.dc_tracking_mu = best.dc_tracking_mu; if obj.dc_tracking_adaptive_enabled - obj.dc_tracking_alpha = best.dc_tracking_alpha; - obj.dc_tracking_gamma = best.dc_tracking_gamma; - obj.dc_tracking_power_exponent = best.dc_tracking_power_exponent; + obj.dc_tracking_persistence_gain = best.dc_tracking_persistence_gain; obj.dc_tracking_mu_eff_max = best.dc_tracking_mu_eff_max; - fprintf("\nFFE DC opt done: dc_tracking_mu=%9.3e, alpha=%6.3f, gamma=%9.3e, p=%5.2f, mu_eff_max=%9.3e, objective=%9.3e\n", ... - obj.dc_tracking_mu,obj.dc_tracking_alpha,obj.dc_tracking_gamma,obj.dc_tracking_power_exponent,obj.dc_tracking_mu_eff_max,obj.dc_tracking_optimization.MinObjective); + fprintf("\nFFE DC opt done: dc_tracking_mu=%9.3e, persistence_gain=%6.3f, mu_eff_max=%9.3e, objective=%9.3e\n", ... + obj.dc_tracking_mu,obj.dc_tracking_persistence_gain,obj.dc_tracking_mu_eff_max,obj.dc_tracking_optimization.MinObjective); else fprintf("\nFFE DC opt done: dc_tracking_mu=%9.3e, objective=%9.3e\n", ... obj.dc_tracking_mu,obj.dc_tracking_optimization.MinObjective); @@ -668,34 +806,52 @@ classdef FFE < handle [initial_weights,stats] = obj.a2LevelWeightInitialGuess(x_opt,d_opt); current_weights = obj.expandA2LevelWeights(n_levels); - initial_matrix = initial_weights(:).'; + initial_matrix = [zeros(1,n_levels); initial_weights(:).']; if any(current_weights ~= 0) initial_matrix = [initial_matrix; current_weights(:).']; end initial_matrix = min(max(initial_matrix,0),obj.a2_level_weight_max); initial_matrix = unique(initial_matrix,"rows","stable"); initial_x = array2table(initial_matrix,"VariableNames",cellstr(var_names)); + [baseline_ber,baseline_errors] = obj.a2LevelWeightBer(initial_x(1,:),x_opt,d_opt); + [x_val,d_val,N_val] = obj.a2ValidationSignals(x,d,N_opt); + stats.baseline_ber = baseline_ber; + stats.baseline_errors = baseline_errors; obj.a2_level_weight_initial_stats = stats; obj.a2_level_weight_optimization_iter = 0; max_evals = max(obj.a2_level_weight_optimization_max_evals,height(initial_x)); fprintf("FFE A2 opt uses fixed mu_tr=%9.3e, mu_dd=%9.3e on %d samples / %d symbols\n", ... obj.mu_tr,obj.mu_dd,N_opt,numel(d_opt)); - fprintf("FFE A2 opt init: var_slope=%9.3e, weights=%s\n", ... - stats.variance_slope,mat2str(initial_weights(:).',3)); + fprintf("FFE A2 opt validation uses %d samples / %d symbols\n",N_val,numel(d_val)); + fprintf("FFE A2 opt init: baseline BER=%9.3e (%d errors), var_slope=%9.3e, weights=%s\n", ... + baseline_ber,baseline_errors,stats.variance_slope,mat2str(initial_weights(:).',3)); - obj.a2_level_weight_optimization = bayesopt(@(p)obj.a2LevelWeightObjective(p,x_opt,d_opt),vars, ... + old_rng = rng; + cleanup_rng = onCleanup(@()rng(old_rng)); + rng(42,"twister"); + obj.a2_level_weight_optimization = bayesopt(@(p)obj.a2LevelWeightObjective(p,x_opt,d_opt,baseline_ber),vars, ... "MaxObjectiveEvaluations",max_evals, ... "InitialX",initial_x, ... "AcquisitionFunctionName","expected-improvement-plus", ... - "IsObjectiveDeterministic",false, ... + "IsObjectiveDeterministic",true, ... "Verbose",0, ... "PlotFcn",[]); + clear cleanup_rng - best = obj.a2_level_weight_optimization.XAtMinObjective; + [best,best_validation_ber,best_validation_errors] = obj.selectA2LevelWeightsByValidation( ... + obj.a2_level_weight_optimization.XTrace, ... + obj.a2_level_weight_optimization.ObjectiveTrace, ... + x_val,d_val,initial_x); + best_opt = obj.a2_level_weight_optimization.XAtMinObjective; + best_opt_weights = obj.a2LevelWeightsFromParams(best_opt); obj.dc_level_weights_a2 = obj.a2LevelWeightsFromParams(best); - fprintf("\nFFE A2 opt done: weights=%s, BER=%9.3e\n", ... - mat2str(obj.dc_level_weights_a2(:).',3),obj.a2_level_weight_optimization.MinObjective); + obj.a2_level_weight_initial_stats.validation_ber = best_validation_ber; + obj.a2_level_weight_initial_stats.validation_errors = best_validation_errors; + obj.a2_level_weight_initial_stats.validation_weights = obj.dc_level_weights_a2; + fprintf("\nFFE A2 opt done: opt_weights=%s, opt_obj=%9.3e, validation_weights=%s, validation_BER=%9.3e (%d errors)\n", ... + mat2str(best_opt_weights(:).',3),obj.a2_level_weight_optimization.MinObjective, ... + mat2str(obj.dc_level_weights_a2(:).',3),best_validation_ber,best_validation_errors); end function [x_opt,d_opt,N_opt] = optimizationSignals(obj,x,d,opt_len) @@ -719,6 +875,47 @@ classdef FFE < handle d_opt = d(1:n_symbols); end + function [x_val,d_val,N_val] = a2ValidationSignals(obj,x,d,N_opt) + N_available = min(numel(x),numel(d) * obj.sps); + N_val = min(N_opt,N_available); + if N_available <= N_opt + [x_val,d_val,N_val] = obj.optimizationSignals(x,d,N_opt); + return + end + + start_symbol = floor((N_available - N_val) / obj.sps) + 1; + start_sample = (start_symbol - 1) * obj.sps + 1; + N_val = obj.sps * floor((N_available - start_sample + 1) / obj.sps); + N_val = max(obj.sps,N_val); + n_symbols = N_val / obj.sps; + x_val = x(start_sample:start_sample+N_val-1); + d_val = d(start_symbol:start_symbol+n_symbols-1); + end + + function [best_params,best_ber,best_errors] = selectA2LevelWeightsByValidation(obj,x_trace,objective_trace,x_val,d_val,initial_x) + objective_trace = objective_trace(:); + objective_trace(~isfinite(objective_trace)) = inf; + [~,sort_idx] = sort(objective_trace,"ascend"); + n_trace_candidates = min(8,numel(sort_idx)); + candidate_x = x_trace(sort_idx(1:n_trace_candidates),:); + candidate_x = [initial_x; candidate_x]; + candidate_x = unique(candidate_x,"rows","stable"); + + n_candidates = height(candidate_x); + validation_ber = inf(n_candidates,1); + validation_errors = nan(n_candidates,1); + for candidate_idx = 1:n_candidates + [validation_ber(candidate_idx),validation_errors(candidate_idx)] = ... + obj.a2LevelWeightBer(candidate_x(candidate_idx,:),x_val,d_val); + end + + [best_ber,best_idx] = min(validation_ber); + best_errors = validation_errors(best_idx); + best_params = candidate_x(best_idx,:); + fprintf("FFE A2 validation: checked %d candidates, best weights=%s, BER=%9.3e, errors=%d\n", ... + n_candidates,mat2str(obj.a2LevelWeightsFromParams(best_params).',3),best_ber,best_errors); + end + function objective = muObjective(obj,params,x,d) old_debug = obj.save_debug; old_dc_tracking_mu = obj.dc_tracking_mu; @@ -762,7 +959,27 @@ classdef FFE < handle obj.dc_tracking_mu = old_dc_tracking_mu; end - function objective = a2LevelWeightObjective(obj,params,x,d) + function objective = a2LevelWeightObjective(obj,params,x,d,baseline_ber) + if nargin < 5 || ~isfinite(baseline_ber) + baseline_ber = inf; + end + + [ber,errors] = obj.a2LevelWeightBer(params,x,d); + objective = ber; + if isfinite(baseline_ber) + objective = objective + max(0,ber - baseline_ber); + end + if ~isfinite(objective) + objective = inf; + end + + obj.a2_level_weight_optimization_iter = obj.a2_level_weight_optimization_iter + 1; + weights = obj.a2LevelWeightsFromParams(params); + fprintf("\rFFE A2 opt %02d: weights=%s, BER=%9.3e, obj=%9.3e, errors=%d", ... + obj.a2_level_weight_optimization_iter,mat2str(weights(:).',3),ber,objective,errors); + end + + function [ber,errors] = a2LevelWeightBer(obj,params,x,d) state = obj.captureObjectiveState(); cleanup = onCleanup(@()obj.restoreObjectiveState(state)); @@ -782,14 +999,6 @@ classdef FFE < handle end [ber,errors] = obj.berObjective(signal,d); - objective = ber; - if ~isfinite(objective) - objective = inf; - end - - obj.a2_level_weight_optimization_iter = obj.a2_level_weight_optimization_iter + 1; - fprintf("\rFFE A2 opt %02d: weights=%s, BER=%9.3e, errors=%d", ... - obj.a2_level_weight_optimization_iter,mat2str(obj.dc_level_weights_a2(:).',3),ber,errors); end function objective = dcTrackingObjective(obj,params,x,d) @@ -821,8 +1030,8 @@ classdef FFE < handle obj.dc_tracking_optimization_iter = obj.dc_tracking_optimization_iter + 1; if obj.dc_tracking_adaptive_enabled - fprintf("\rFFE DC opt %02d: dc_tracking_mu=%9.3e, alpha=%6.3f, gamma=%9.3e, p=%5.2f, mu_eff_max=%9.3e, BER=%9.3e, delay=%7.0f, corr=%6.3f, obj=%9.3e, errors=%d", ... - obj.dc_tracking_optimization_iter,params.dc_tracking_mu,params.dc_tracking_alpha,params.dc_tracking_gamma,params.dc_tracking_power_exponent,params.dc_tracking_mu_eff_max, ... + fprintf("\rFFE DC opt %02d: dc_tracking_mu=%9.3e, persistence_gain=%6.3f, mu_eff_max=%9.3e, BER=%9.3e, delay=%7.0f, corr=%6.3f, obj=%9.3e, errors=%d", ... + obj.dc_tracking_optimization_iter,params.dc_tracking_mu,params.dc_tracking_persistence_gain,params.dc_tracking_mu_eff_max, ... ber,delay_symbols,delay_corr,objective,errors); else fprintf("\rFFE DC opt %02d: dc_tracking_mu=%9.3e, BER=%9.3e, delay=%7.0f, corr=%6.3f, obj=%9.3e, errors=%d", ... @@ -835,14 +1044,8 @@ classdef FFE < handle if any(var_names == "dc_tracking_mu") obj.dc_tracking_mu = params.dc_tracking_mu; end - if any(var_names == "dc_tracking_alpha") - obj.dc_tracking_alpha = params.dc_tracking_alpha; - end - if any(var_names == "dc_tracking_gamma") - obj.dc_tracking_gamma = params.dc_tracking_gamma; - end - if any(var_names == "dc_tracking_power_exponent") - obj.dc_tracking_power_exponent = params.dc_tracking_power_exponent; + if any(var_names == "dc_tracking_persistence_gain") + obj.dc_tracking_persistence_gain = params.dc_tracking_persistence_gain; end if any(var_names == "dc_tracking_mu_eff_max") obj.dc_tracking_mu_eff_max = params.dc_tracking_mu_eff_max; @@ -950,6 +1153,51 @@ classdef FFE < handle "initial_weights",weights); end + function [e_dc_next,stats] = dcTrackingBlockUpdate(obj,e_dc_current,err_block,options) + arguments + obj + e_dc_current (1,1) double + err_block (:,1) double + options.mu_dc (1,1) double = NaN + options.persistence_gain (1,1) double = 0 + options.mu_eff_min (1,1) double = NaN + options.mu_eff_max (1,1) double = NaN + end + + if isnan(options.mu_dc) + options.mu_dc = obj.dc_tracking_mu; + end + if isnan(options.mu_eff_min) + options.mu_eff_min = obj.dc_tracking_mu_eff_min; + end + if isnan(options.mu_eff_max) + options.mu_eff_max = obj.dc_tracking_mu_eff_max; + end + + valid_err = err_block(isfinite(err_block)); + if isempty(valid_err) + err_mean = 0; + err_abs_mean = 0; + else + err_mean = mean(valid_err,"omitnan"); + err_abs_mean = mean(abs(valid_err),"omitnan"); + end + + persistence_scale = abs(err_mean) / (err_abs_mean + eps); + persistence_scale = min(max(persistence_scale,0),1); + mu_eff = options.mu_dc * (1 + max(options.persistence_gain,0) * persistence_scale); + mu_eff = min(max(mu_eff,options.mu_eff_min),options.mu_eff_max); + update = mu_eff * err_mean; + e_dc_next = e_dc_current + update; + + stats = struct( ... + "err_mean",err_mean, ... + "err_abs_mean",err_abs_mean, ... + "persistence_scale",persistence_scale, ... + "mu_eff",mu_eff, ... + "update",update); + end + function state = captureObjectiveState(obj) state.e = obj.e; state.e_dc = obj.e_dc; @@ -959,6 +1207,7 @@ classdef FFE < handle state.dc_tracking_mu = obj.dc_tracking_mu; state.dc_tracking_alpha = obj.dc_tracking_alpha; state.dc_tracking_gamma = obj.dc_tracking_gamma; + state.dc_tracking_persistence_gain = obj.dc_tracking_persistence_gain; state.dc_tracking_power_exponent = obj.dc_tracking_power_exponent; state.dc_tracking_mu_eff_max = obj.dc_tracking_mu_eff_max; state.dc_level_weights_a2 = obj.dc_level_weights_a2; @@ -974,6 +1223,7 @@ classdef FFE < handle obj.dc_tracking_mu = state.dc_tracking_mu; obj.dc_tracking_alpha = state.dc_tracking_alpha; obj.dc_tracking_gamma = state.dc_tracking_gamma; + obj.dc_tracking_persistence_gain = state.dc_tracking_persistence_gain; obj.dc_tracking_power_exponent = state.dc_tracking_power_exponent; obj.dc_tracking_mu_eff_max = state.dc_tracking_mu_eff_max; obj.dc_level_weights_a2 = state.dc_level_weights_a2; @@ -989,8 +1239,8 @@ classdef FFE < handle rx_bits = mapper.demap(eq_signal_hd); tx_bits = mapper.demap(tx_symbols); [~,errors,ber,~] = calc_ber(rx_bits.signal,tx_bits.signal, ... - "skip_front",10, ... - "skip_end",10, ... + "skip_front",1000, ... + "skip_end",0, ... "returnErrorLocation",1); end diff --git a/Functions/EQ_recipes/mpi_recipe.m b/Functions/EQ_recipes/mpi_recipe.m new file mode 100644 index 0000000..cb542a2 --- /dev/null +++ b/Functions/EQ_recipes/mpi_recipe.m @@ -0,0 +1,404 @@ +function output = mpi_recipe_dev(Scpe_sig_raw, Symbols, Tx_bits, options) +%mpi_recipe_dev Minimal example recipe for the DSP job framework. +% This recipe intentionally performs only light preprocessing and records +% summary values. It demonstrates the recipe interface without running a +% full equalizer chain. + +arguments + Scpe_sig_raw + Symbols + Tx_bits + options.fsym + options.M + options.duob_mode + options.dataTable table + options.userParameters struct = struct() + options.debug_plots (1,1) logical = false +end + +Scpe_sig = preprocessSignal(Scpe_sig_raw, Symbols, options.fsym, ... + "mode", "auto", ... + "debug_plots", options.debug_plots); +output = struct(); + +%% conventional FFE +eq_core_settings = { ... + "sps", 2, ... + "order", 25, ... + "decide", 0, ... + "adaption_technique", "nlms"}; + +eq_training_settings = { ... + "len_tr", 4096, ... + "epochs_tr", 5, ... + "mu_tr", 0.04}; + +eq_dd_settings = { ... + "dd_mode", 1, ... + "epochs_dd", 3, ... + "mu_dd", 0.012}; + +eq_a1_settings = { ... + "dc_smoothing_a1", 0, ... + "dc_avg_bufferlength_a1", 0, ... + "dc_avg_update_blocklength_a1", 1}; + +eq_a2_settings = { ... + "dc_smoothing_a2", 0, ... + "dc_level_avg_bufferlength_a2", 0, ... + "dc_level_update_blocklength_a2", 0, ... + "dc_level_weights_a2", [0.6 0.6 0.6 0.6]}; + +eq_dc_tracking_settings = { ... + "dc_tracking_mu", 0, ... + "dc_tracking_adaptive_enabled", 0, ... + "dc_tracking_persistence_gain", 0, ... + "dc_tracking_buffer_len", 0}; + +eq_ffe_update_settings = { ... + "ffe_update_buffer_len", 1}; + +eq_optimizer_settings = { ... + "optmize_mus", 0, ... + "plot_mu_optimization", options.debug_plots, ... + "optimize_dc_tracking_params", 0, ... + "optimize_a2_level_weights", 0, ... + "a2_level_weight_optimization_len", 2^15, ... + "a2_level_weight_optimization_max_evals", 30, ... + "a2_level_weight_max", 1}; + +eq_debug_settings = { ... + "save_debug", false}; + +eq_settings = [ ... + eq_core_settings, ... + eq_training_settings, ... + eq_dd_settings, ... + eq_a1_settings, ... + eq_a2_settings, ... + eq_dc_tracking_settings, ... + eq_ffe_update_settings, ... + eq_optimizer_settings, ... + eq_debug_settings]; + +eq_ffe = FFE(eq_settings{:}); + +if 1 + + [ffe_results,equalized_signal] = ffe(eq_ffe, options.M, Scpe_sig, Symbols, Tx_bits, ... + "precode_mode", options.duob_mode, ... + 'showAnalysis', options.debug_plots, ... + "postFFE", [], ... + "eth_style_symbol_mapping", 0); + + ffe_results.metrics.print("description",sprintf('Immediate A2; SIR %d dB',options.dataTable.sir)); + output.conventional_ffe = ffe_results; + +end + +if options.debug_plots + showLevelScatter(Scpe_sig_raw, Symbols, ... + "fsym", options.fsym, ... + "fignum", 400, ... + "normalize", true); + + [~, avg_for_lvl] = showLevelScatter(equalized_signal, Symbols, ... + "fsym", options.fsym, ... + "fignum", 401, ... + "normalize", true); + +end + + +%% A2 +eq_core_settings = { ... + "sps", 2, ... + "order", 25, ... + "decide", 0, ... + "adaption_technique", "nlms"}; + +eq_training_settings = { ... + "len_tr", 4096, ... + "epochs_tr", 5, ... + "mu_tr", 0.04}; + +eq_dd_settings = { ... + "dd_mode", 1, ... + "epochs_dd", 3, ... + "mu_dd", 0.012}; + +eq_a1_settings = { ... + "dc_smoothing_a1", 0, ... + "dc_avg_bufferlength_a1", 0, ... + "dc_avg_update_blocklength_a1", 1}; + +eq_a2_settings = { ... + "dc_smoothing_a2", 1, ... + "dc_level_avg_bufferlength_a2", 256, ... + "dc_level_update_blocklength_a2", options.userParameters.block_update, ... + "dc_level_weights_a2", [0.6 0.6 0.6 0.6]}; + +eq_dc_tracking_settings = { ... + "dc_tracking_mu", 0, ... + "dc_tracking_adaptive_enabled", 0, ... + "dc_tracking_persistence_gain", 0, ... + "dc_tracking_buffer_len", 0}; + +eq_ffe_update_settings = { ... + "ffe_update_buffer_len", 1}; + +eq_optimizer_settings = { ... + "optmize_mus", 0, ... + "plot_mu_optimization", options.debug_plots, ... + "optimize_dc_tracking_params", 0, ... + "optimize_a2_level_weights", 1, ... + "a2_level_weight_optimization_len", 2^15, ... + "a2_level_weight_optimization_max_evals", 30, ... + "a2_level_weight_max", 1}; + +eq_debug_settings = { ... + "save_debug", false}; + +eq_settings = [ ... + eq_core_settings, ... + eq_training_settings, ... + eq_dd_settings, ... + eq_a1_settings, ... + eq_a2_settings, ... + eq_dc_tracking_settings, ... + eq_ffe_update_settings, ... + eq_optimizer_settings, ... + eq_debug_settings]; + +eq_ffe = FFE(eq_settings{:}); + +if 1 + + [ffe_results,equalized_signal] = ffe(eq_ffe, options.M, Scpe_sig, Symbols, Tx_bits, ... + "precode_mode", options.duob_mode, ... + 'showAnalysis', options.debug_plots, ... + "postFFE", [], ... + "eth_style_symbol_mapping", 0); + + ffe_results.metrics.print("description",sprintf('Immediate A2; SIR %d dB',options.dataTable.sir)); + output.a2_immediate_updates = ffe_results; + +end + +if options.debug_plots + showLevelScatter(Scpe_sig_raw, Symbols, ... + "fsym", options.fsym, ... + "fignum", 400, ... + "normalize", true); + + [~, avg_for_lvl] = showLevelScatter(equalized_signal, Symbols, ... + "fsym", options.fsym, ... + "fignum", 401, ... + "normalize", true); + +end + +%% A1 + +eq_a1_settings = { ... + "dc_smoothing_a1", 1 ... + "dc_avg_bufferlength_a1", 2048, ... + "dc_avg_update_blocklength_a1", options.userParameters.block_update}; + +eq_a2_settings = { ... + "dc_smoothing_a2", 0, ... + "dc_level_avg_bufferlength_a2", 0, ... + "dc_level_update_blocklength_a2", 0, ... + "dc_level_weights_a2", [0.6 0.6 0.6 0.6]}; + +eq_optimizer_settings = { ... + "optmize_mus", 0, ... + "plot_mu_optimization", options.debug_plots, ... + "optimize_dc_tracking_params", 0, ... + "optimize_a2_level_weights", 0, ... + "a2_level_weight_optimization_len", 2^15, ... + "a2_level_weight_optimization_max_evals", 30, ... + "a2_level_weight_max", 1}; + +eq_debug_settings = { ... + "save_debug", false}; + +eq_settings = [ ... + eq_core_settings, ... + eq_training_settings, ... + eq_dd_settings, ... + eq_a1_settings, ... + eq_a2_settings, ... + eq_dc_tracking_settings, ... + eq_ffe_update_settings, ... + eq_optimizer_settings, ... + eq_debug_settings]; + +eq_ffe = FFE(eq_settings{:}); + +if 1 + + [ffe_results,equalized_signal] = ffe(eq_ffe, options.M, Scpe_sig, Symbols, Tx_bits, ... + "precode_mode", options.duob_mode, ... + 'showAnalysis', options.debug_plots, ... + "postFFE", [], ... + "eth_style_symbol_mapping", 0); + + ffe_results.metrics.print("description",sprintf('Immediate A1; SIR %d dB',options.dataTable.sir)); + output.a1_immediate_updates = ffe_results; + +end + +%% Tracking + +eq_core_settings = { ... + "sps", 2, ... + "order", 25, ... + "decide", 0, ... + "adaption_technique", "nlms"}; + +eq_training_settings = { ... + "len_tr", 4096, ... + "epochs_tr", 5, ... + "mu_tr", 0.04}; + +eq_dd_settings = { ... + "dd_mode", 1, ... + "epochs_dd", 3, ... + "mu_dd", 0.012}; + +eq_a1_settings = { ... + "dc_smoothing_a1", 0, ... + "dc_avg_bufferlength_a1", 0, ... + "dc_avg_update_blocklength_a1", 1}; + +eq_a2_settings = { ... + "dc_smoothing_a2", 0, ... + "dc_level_avg_bufferlength_a2", 0, ... + "dc_level_update_blocklength_a2", 0, ... + "dc_level_weights_a2", [0.6 0.6 0.6 0.6]}; + +eq_dc_tracking_settings = { ... + "dc_tracking_mu", 0.002, ... + "dc_tracking_adaptive_enabled", 0, ... + "dc_tracking_persistence_gain", 0, ... + "dc_tracking_buffer_len", options.userParameters.block_update}; + +eq_ffe_update_settings = { ... + "ffe_update_buffer_len", 1}; + +eq_optimizer_settings = { ... + "optmize_mus", 0, ... + "plot_mu_optimization", options.debug_plots, ... + "optimize_dc_tracking_params", 1, ... + "optimize_a2_level_weights", 0, ... + "a2_level_weight_optimization_len", 2^15, ... + "a2_level_weight_optimization_max_evals", 30, ... + "a2_level_weight_max", 1}; + +eq_debug_settings = { ... + "save_debug", false}; + +eq_settings = [ ... + eq_core_settings, ... + eq_training_settings, ... + eq_dd_settings, ... + eq_a1_settings, ... + eq_a2_settings, ... + eq_dc_tracking_settings, ... + eq_ffe_update_settings, ... + eq_optimizer_settings, ... + eq_debug_settings]; + +eq_ffe = FFE(eq_settings{:}); + +if 1 + + [ffe_results,equalized_signal] = ffe(eq_ffe, options.M, Scpe_sig, Symbols, Tx_bits, ... + "precode_mode", options.duob_mode, ... + 'showAnalysis', options.debug_plots, ... + "postFFE", [], ... + "eth_style_symbol_mapping", 0); + + ffe_results.metrics.print("description",sprintf('Immediate A2; SIR %d dB',options.dataTable.sir)); + output.tracking_fixmu_immediate_updates = ffe_results; + +end + +%% Tracking Adaptive + +eq_core_settings = { ... + "sps", 2, ... + "order", 25, ... + "decide", 0, ... + "adaption_technique", "nlms"}; + +eq_training_settings = { ... + "len_tr", 4096, ... + "epochs_tr", 5, ... + "mu_tr", 0.04}; + +eq_dd_settings = { ... + "dd_mode", 1, ... + "epochs_dd", 3, ... + "mu_dd", 0.012}; + +eq_a1_settings = { ... + "dc_smoothing_a1", 0, ... + "dc_avg_bufferlength_a1", 0, ... + "dc_avg_update_blocklength_a1", 1}; + +eq_a2_settings = { ... + "dc_smoothing_a2", 0, ... + "dc_level_avg_bufferlength_a2", 0, ... + "dc_level_update_blocklength_a2", 0, ... + "dc_level_weights_a2", [0.6 0.6 0.6 0.6]}; + +eq_dc_tracking_settings = { ... + "dc_tracking_mu", 0.002, ... + "dc_tracking_adaptive_enabled", 1, ... + "dc_tracking_persistence_gain", 0, ... + "dc_tracking_buffer_len", options.userParameters.block_update}; + +eq_ffe_update_settings = { ... + "ffe_update_buffer_len", 1}; + +eq_optimizer_settings = { ... + "optmize_mus", 0, ... + "plot_mu_optimization", options.debug_plots, ... + "optimize_dc_tracking_params", 1, ... + "optimize_a2_level_weights", 0, ... + "a2_level_weight_optimization_len", 2^15, ... + "a2_level_weight_optimization_max_evals", 30, ... + "a2_level_weight_max", 1}; + +eq_debug_settings = { ... + "save_debug", false}; + +eq_settings = [ ... + eq_core_settings, ... + eq_training_settings, ... + eq_dd_settings, ... + eq_a1_settings, ... + eq_a2_settings, ... + eq_dc_tracking_settings, ... + eq_ffe_update_settings, ... + eq_optimizer_settings, ... + eq_debug_settings]; + +eq_ffe = FFE(eq_settings{:}); + +if 1 + + [ffe_results,equalized_signal] = ffe(eq_ffe, options.M, Scpe_sig, Symbols, Tx_bits, ... + "precode_mode", options.duob_mode, ... + 'showAnalysis', options.debug_plots, ... + "postFFE", [], ... + "eth_style_symbol_mapping", 0); + + ffe_results.metrics.print("description",sprintf('Immediate A2; SIR %d dB',options.dataTable.sir)); + output.tracking_adaptive_immediate_updates = ffe_results; + +end + +end diff --git a/Functions/EQ_recipes/mpi_recipe_dev.m b/Functions/EQ_recipes/mpi_recipe_dev.m index bb0fc13..bf50562 100644 --- a/Functions/EQ_recipes/mpi_recipe_dev.m +++ b/Functions/EQ_recipes/mpi_recipe_dev.m @@ -4,202 +4,469 @@ function output = mpi_recipe_dev(Scpe_sig_raw, Symbols, Tx_bits, options) % summary values. It demonstrates the recipe interface without running a % full equalizer chain. - arguments - Scpe_sig_raw - Symbols - Tx_bits - options.fsym - options.M - options.duob_mode - options.dataTable table - options.userParameters struct = struct() - options.debug_plots (1,1) logical = false - end +arguments + Scpe_sig_raw + Symbols + Tx_bits + options.fsym + options.M + options.duob_mode + options.dataTable table + options.userParameters struct = struct() + options.debug_plots (1,1) logical = false +end - Scpe_sig = preprocessSignal(Scpe_sig_raw, Symbols, options.fsym, ... - "mode", "auto", ... - "debug_plots", options.debug_plots); - output = struct(); +Scpe_sig = preprocessSignal(Scpe_sig_raw, Symbols, options.fsym, ... + "mode", "auto", ... + "debug_plots", options.debug_plots); +output = struct(); + +%% conventional FFE +eq_core_settings = { ... + "sps", 2, ... + "order", 50, ... + "decide", 0, ... + "adaption_technique", "nlms"}; + +eq_training_settings = { ... + "len_tr", 4096, ... + "epochs_tr", 5, ... + "mu_tr", 0.04}; + +eq_dd_settings = { ... + "dd_mode", 1, ... + "epochs_dd", 3, ... + "mu_dd", 0.012}; + +eq_a1_settings = { ... + "dc_smoothing_a1", 0, ... + "dc_avg_bufferlength_a1", 0, ... + "dc_avg_update_blocklength_a1", 1}; + +eq_a2_settings = { ... + "dc_smoothing_a2", 0, ... + "dc_level_avg_bufferlength_a2", 0, ... + "dc_level_update_blocklength_a2", 0, ... + "dc_level_weights_a2", [0.6 0.6 0.6 0.6]}; + +eq_dc_tracking_settings = { ... + "dc_tracking_mu", 0, ... + "dc_tracking_adaptive_enabled", 0, ... + "dc_tracking_persistence_gain", 0, ... + "dc_tracking_buffer_len", 1}; + +eq_ffe_update_settings = { ... + "ffe_update_buffer_len", 0}; + +eq_optimizer_settings = { ... + "optmize_mus", 0, ... + "plot_mu_optimization", options.debug_plots, ... + "optimize_dc_tracking_params", 0, ... + "optimize_a2_level_weights", 0, ... + "a2_level_weight_optimization_len", 2^15, ... + "a2_level_weight_optimization_max_evals", 30, ... + "a2_level_weight_max", 1}; + +eq_debug_settings = { ... + "save_debug", false}; + +eq_settings = [ ... + eq_core_settings, ... + eq_training_settings, ... + eq_dd_settings, ... + eq_a1_settings, ... + eq_a2_settings, ... + eq_dc_tracking_settings, ... + eq_ffe_update_settings, ... + eq_optimizer_settings, ... + eq_debug_settings]; + +eq_ffe = FFE(eq_settings{:}); + +% ffe_order_ffe = [50, 0, 0]; +% mu_ffe = [0.0001, 0.0008, 0.001]; +% mu_dfe = 0.0004; +% mu_ffe = [0, 0,0]; +% mu_dfe = 0.0004; +% mu_dc = 0.005; +% eq_ffe = EQ("Ne",ffe_order_ffe,"Nb",[0,0,0],"training_length",4096,"training_loops",5,"dd_loops",5,"K",2,"DCmu",0.00,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",0); + + +if 1 + + [ffe_results,equalized_signal] = ffe(eq_ffe, options.M, Scpe_sig, Symbols, Tx_bits, ... + "precode_mode", options.duob_mode, ... + 'showAnalysis', options.debug_plots, ... + "postFFE", [], ... + "eth_style_symbol_mapping", 0); + + ffe_results.metrics.print("description",sprintf('Immediate A2; SIR %d dB',options.dataTable.sir)); + output.conventional_ffe = ffe_results; + +end + +if 1 - eq_core_settings = { ... - "sps", 2, ... - "order", 25, ... - "decide", 0, ... - "adaption_technique", "nlms"}; + showLevelScatter(Scpe_sig_raw, Symbols, ... + "fsym", options.fsym, ... + "fignum", 400, ... + "normalize", true, ... + "debug_plots", false, ... + "showPlot", true, ... + "showStdAnnotations", true, ... + "yLimits", [-2.5 2.5], ... + "xLimits",[0 2.3],... + "scatterAlpha", 0.25, ... + "scatterSize", 1, ... + "avgLineMaxPoints", 500, ... + "avgLineSmoothWindow", 5); - eq_training_settings = { ... - "len_tr", 4096, ... - "epochs_tr", 5, ... - "mu_tr", 0.04}; - eq_dd_settings = { ... - "dd_mode", 1, ... - "epochs_dd", 5, ... - "mu_dd", 0.012}; + exportName = "before_eq"; + exportDir = "C:/Users/Silas/Documents/6971e0b65b380ca6d71c837f/04_Experimental_Evaluation/tikz/mpi"; + pngFile = fullfile(exportDir, exportName + "_img.png"); % absolute filesystem output + tikzFile = fullfile(exportDir, exportName + ".tikz"); % absolute filesystem output + pngTikzPath = "04_Experimental_Evaluation/tikz/mpi/" + exportName + "_img.png"; % path written into TikZ + figure(400); + exportLevelScatterTikz(gca, pngFile, tikzFile, pngTikzPath, ... + "resolutionDpi", 150, ... + "generatedBy", "mpi_recipe_dev.m"); - eq_a1_settings = { ... - "dc_smoothing_a1", 0, ... - "dc_avg_bufferlength_a1", 0}; + showLevelScatter(equalized_signal.*-1, Symbols, ... + "fsym", options.fsym, ... + "fignum", 401, ... + "normalize", true, ... + "debug_plots", false, ... + "showPlot", true, ... + "showStdAnnotations", true, ... + "yLimits", [-2.5 2.5], ... + "xLimits",[0 2.3],... + "scatterAlpha", 0.25, ... + "scatterSize", 1, ... + "avgLineMaxPoints", 500, ... + "avgLineSmoothWindow", 5); - eq_a2_settings = { ... - "dc_smoothing_a2", 0, ... - "dc_level_avg_bufferlength_a2", 0, ... - "dc_level_weights_a2", [0.692 0.979 0.98 0.138]}; + exportName = "after_ffe_eq"; + pngFile = fullfile(exportDir, exportName + "_img.png"); % absolute filesystem output + tikzFile = fullfile(exportDir, exportName + ".tikz"); % absolute filesystem output + pngTikzPath = "04_Experimental_Evaluation/tikz/mpi/" + exportName + "_img.png"; % path written into TikZ + figure(401); + exportLevelScatterTikz(gca, pngFile, tikzFile, pngTikzPath, ... + "resolutionDpi", 150, ... + "generatedBy", "mpi_recipe_dev.m"); - eq_dc_tracking_settings = { ... - "dc_tracking_mu", 0.05, ... - "dc_tracking_adaptive_enabled", 0, ... - "dc_tracking_buffer_len", 1024}; +end - eq_ffe_update_settings = { ... - "ffe_update_buffer_len", 1}; - eq_optimizer_settings = { ... - "optmize_mus", 0, ... - "plot_mu_optimization", options.debug_plots, ... - "optimize_dc_tracking_params", 1, ... - "optimize_a2_level_weights", 0, ... - "a2_level_weight_optimization_len", 2^15, ... - "a2_level_weight_optimization_max_evals", 30, ... - "a2_level_weight_max", 1}; +%% A2 +eq_core_settings = { ... + "sps", 2, ... + "order", 50, ... + "decide", 0, ... + "adaption_technique", "nlms"}; - eq_debug_settings = { ... - "save_debug", true}; +eq_training_settings = { ... + "len_tr", 4096, ... + "epochs_tr", 5, ... + "mu_tr", 0.04}; - eq_settings = [ ... - eq_core_settings, ... - eq_training_settings, ... - eq_dd_settings, ... - eq_a1_settings, ... - eq_a2_settings, ... - eq_dc_tracking_settings, ... - eq_ffe_update_settings, ... - eq_optimizer_settings, ... - eq_debug_settings]; +eq_dd_settings = { ... + "dd_mode", 1, ... + "epochs_dd", 3, ... + "mu_dd", 0.012}; - eq_ffe = FFE(eq_settings{:}); +eq_a1_settings = { ... + "dc_smoothing_a1", 0, ... + "dc_avg_bufferlength_a1", 0, ... + "dc_avg_update_blocklength_a1", 1}; + +eq_a2_settings = { ... + "dc_smoothing_a2", 1, ... + "dc_level_avg_bufferlength_a2", 256, ... + "dc_level_update_blocklength_a2", options.userParameters.block_update, ... + "dc_level_weights_a2", [0.6]}; + +eq_dc_tracking_settings = { ... + "dc_tracking_mu", 0, ... + "dc_tracking_adaptive_enabled", 0, ... + "dc_tracking_persistence_gain", 0, ... + "dc_tracking_buffer_len", 0}; + +eq_ffe_update_settings = { ... + "ffe_update_buffer_len", 1}; + +eq_optimizer_settings = { ... + "optmize_mus", 0, ... + "plot_mu_optimization", options.debug_plots, ... + "optimize_dc_tracking_params", 0, ... + "optimize_a2_level_weights", 1, ... + "a2_level_weight_optimization_len", 2^15, ... + "a2_level_weight_optimization_max_evals", 30, ... + "a2_level_weight_max", 1}; + +eq_debug_settings = { ... + "save_debug", false}; + +eq_settings = [ ... + eq_core_settings, ... + eq_training_settings, ... + eq_dd_settings, ... + eq_a1_settings, ... + eq_a2_settings, ... + eq_dc_tracking_settings, ... + eq_ffe_update_settings, ... + eq_optimizer_settings, ... + eq_debug_settings]; + +eq_ffe = FFE(eq_settings{:}); + +if 0 + + [ffe_results,equalized_signal] = ffe(eq_ffe, options.M, Scpe_sig, Symbols, Tx_bits, ... + "precode_mode", options.duob_mode, ... + 'showAnalysis', options.debug_plots, ... + "postFFE", [], ... + "eth_style_symbol_mapping", 0); + ffe_results.metrics.print("description",sprintf('Immediate A2; SIR %d dB',options.dataTable.sir)); + output.a2_immediate_updates = ffe_results; + +end + +if options.debug_plots showLevelScatter(Scpe_sig_raw, Symbols, ... "fsym", options.fsym, ... "fignum", 400, ... "normalize", true); + [~, avg_for_lvl] = showLevelScatter(equalized_signal, Symbols, ... + "fsym", options.fsym, ... + "fignum", 401, ... + "normalize", true); + +end + +%% A1 + +eq_a1_settings = { ... + "dc_smoothing_a1", 1 ... + "dc_avg_bufferlength_a1", 2048, ... + "dc_avg_update_blocklength_a1", options.userParameters.block_update}; + +eq_a2_settings = { ... + "dc_smoothing_a2", 0, ... + "dc_level_avg_bufferlength_a2", 0, ... + "dc_level_update_blocklength_a2", 0, ... + "dc_level_weights_a2", [0.6 0.6 0.6 0.6]}; + +eq_optimizer_settings = { ... + "optmize_mus", 0, ... + "plot_mu_optimization", options.debug_plots, ... + "optimize_dc_tracking_params", 0, ... + "optimize_a2_level_weights", 0, ... + "a2_level_weight_optimization_len", 2^15, ... + "a2_level_weight_optimization_max_evals", 30, ... + "a2_level_weight_max", 1}; + +eq_debug_settings = { ... + "save_debug", false}; + +eq_settings = [ ... + eq_core_settings, ... + eq_training_settings, ... + eq_dd_settings, ... + eq_a1_settings, ... + eq_a2_settings, ... + eq_dc_tracking_settings, ... + eq_ffe_update_settings, ... + eq_optimizer_settings, ... + eq_debug_settings]; + +eq_ffe = FFE(eq_settings{:}); + +if 0 + + [ffe_results,equalized_signal] = ffe(eq_ffe, options.M, Scpe_sig, Symbols, Tx_bits, ... + "precode_mode", options.duob_mode, ... + 'showAnalysis', options.debug_plots, ... + "postFFE", [], ... + "eth_style_symbol_mapping", 0); + + ffe_results.metrics.print("description",sprintf('Immediate A1; SIR %d dB',options.dataTable.sir)); + output.a1_immediate_updates = ffe_results; + + showLevelScatter(Scpe_sig_raw, Symbols, ... + "fsym", options.fsym, ... + "fignum", 400, ... + "normalize", true); + + [~, avg_for_lvl] = showLevelScatter(equalized_signal, Symbols, ... + "fsym", options.fsym, ... + "fignum", 401, ... + "normalize", true); +end + +%% Tracking fixed + +eq_core_settings = { ... + "sps", 2, ... + "order", 50, ... + "decide", 0, ... + "adaption_technique", "nlms"}; + +eq_training_settings = { ... + "len_tr", 4096, ... + "epochs_tr", 5, ... + "mu_tr", 0.04}; + +eq_dd_settings = { ... + "dd_mode", 1, ... + "epochs_dd", 3, ... + "mu_dd", 0.012}; + +eq_a1_settings = { ... + "dc_smoothing_a1", 0, ... + "dc_avg_bufferlength_a1", 0, ... + "dc_avg_update_blocklength_a1", 1}; + +eq_a2_settings = { ... + "dc_smoothing_a2", 0, ... + "dc_level_avg_bufferlength_a2", 0, ... + "dc_level_update_blocklength_a2", 0, ... + "dc_level_weights_a2", [0.6 0.6 0.6 0.6]}; + +eq_dc_tracking_settings = { ... + "dc_tracking_mu", 0.002, ... + "dc_tracking_adaptive_enabled", 0, ... + "dc_tracking_persistence_gain", 0, ... + "dc_tracking_buffer_len", options.userParameters.block_update}; + +eq_ffe_update_settings = { ... + "ffe_update_buffer_len", 1}; + +eq_optimizer_settings = { ... + "optmize_mus", 0, ... + "plot_mu_optimization", options.debug_plots, ... + "optimize_dc_tracking_params", 1, ... + "optimize_a2_level_weights", 0, ... + "a2_level_weight_optimization_len", 2^15, ... + "a2_level_weight_optimization_max_evals", 30, ... + "a2_level_weight_max", 1}; + +eq_debug_settings = { ... + "save_debug", false}; + +eq_settings = [ ... + eq_core_settings, ... + eq_training_settings, ... + eq_dd_settings, ... + eq_a1_settings, ... + eq_a2_settings, ... + eq_dc_tracking_settings, ... + eq_ffe_update_settings, ... + eq_optimizer_settings, ... + eq_debug_settings]; + +eq_ffe = FFE(eq_settings{:}); + +if 1 + + [ffe_results,equalized_signal] = ffe(eq_ffe, options.M, Scpe_sig, Symbols, Tx_bits, ... + "precode_mode", options.duob_mode, ... + 'showAnalysis', options.debug_plots, ... + "postFFE", [], ... + "eth_style_symbol_mapping", 0); + + ffe_results.metrics.print("description",sprintf('Immediate A2; SIR %d dB',options.dataTable.sir)); + output.tracking_fixmu_immediate_updates = ffe_results; + +end + +%% Tracking Adaptive + +eq_core_settings = { ... + "sps", 2, ... + "order", 50, ... + "decide", 0, ... + "adaption_technique", "nlms"}; + +eq_training_settings = { ... + "len_tr", 4096, ... + "epochs_tr", 5, ... + "mu_tr", 0.04}; + +eq_dd_settings = { ... + "dd_mode", 1, ... + "epochs_dd", 3, ... + "mu_dd", 0.012}; + +eq_a1_settings = { ... + "dc_smoothing_a1", 0, ... + "dc_avg_bufferlength_a1", 0, ... + "dc_avg_update_blocklength_a1", 1}; + +eq_a2_settings = { ... + "dc_smoothing_a2", 0, ... + "dc_level_avg_bufferlength_a2", 0, ... + "dc_level_update_blocklength_a2", 0, ... + "dc_level_weights_a2", [0.6 0.6 0.6 0.6]}; + +eq_dc_tracking_settings = { ... + "dc_tracking_mu", 0.002, ... + "dc_tracking_adaptive_enabled", 1, ... + "dc_tracking_persistence_gain", 0, ... + "dc_tracking_buffer_len", options.userParameters.block_update}; + +eq_ffe_update_settings = { ... + "ffe_update_buffer_len", 1}; + +eq_optimizer_settings = { ... + "optmize_mus", 0, ... + "plot_mu_optimization", options.debug_plots, ... + "optimize_dc_tracking_params", 1, ... + "optimize_a2_level_weights", 0, ... + "a2_level_weight_optimization_len", 2^15, ... + "a2_level_weight_optimization_max_evals", 30, ... + "a2_level_weight_max", 1}; + +eq_debug_settings = { ... + "save_debug", false}; + +eq_settings = [ ... + eq_core_settings, ... + eq_training_settings, ... + eq_dd_settings, ... + eq_a1_settings, ... + eq_a2_settings, ... + eq_dc_tracking_settings, ... + eq_ffe_update_settings, ... + eq_optimizer_settings, ... + eq_debug_settings]; + +eq_ffe = FFE(eq_settings{:}); + +if 1 + + [ffe_results,equalized_signal] = ffe(eq_ffe, options.M, Scpe_sig, Symbols, Tx_bits, ... + "precode_mode", options.duob_mode, ... + 'showAnalysis', options.debug_plots, ... + "postFFE", [], ... + "eth_style_symbol_mapping", 0); + + ffe_results.metrics.print("description",sprintf('Immediate A2; SIR %d dB',options.dataTable.sir)); + output.tracking_adaptive_immediate_updates = ffe_results; + + showLevelScatter(Scpe_sig_raw, Symbols, ... + "fsym", options.fsym, ... + "fignum", 400, ... + "normalize", true); + + [~, avg_for_lvl] = showLevelScatter(equalized_signal, Symbols, ... + "fsym", options.fsym, ... + "fignum", 401, ... + "normalize", true); + + +end - %% NORMAL FFE - if 1 - - [ffe_results,equalized_signal] = ffe(eq_ffe, options.M, Scpe_sig, Symbols, Tx_bits, ... - "precode_mode", options.duob_mode, ... - 'showAnalysis', options.debug_plots, ... - "postFFE", [], ... - "eth_style_symbol_mapping", 0); - - ffe_results.metrics.print("description",sprintf('Normal FFE; SIR %d dB',options.dataTable.sir)); - output.ffe_package = ffe_results; - - [~, avg_for_lvl] = showLevelScatter(equalized_signal, Symbols, ... - "fsym", options.fsym, ... - "fignum", 401, ... - "normalize", true); - - output.ffe_debug = struct(); - dbg = eq_ffe.debug_struct; - - - - smooth_len = 501; - output.ffe_debug.smoothing_window = smooth_len; - error_first_s = movmean(dbg.error_first_epoch(:),smooth_len,"omitnan"); - error_last_s = movmean(dbg.error(:),smooth_len,"omitnan"); - dc_tracking_mu_eff_s = movmean(dbg.dc_tracking_mu_eff(:),smooth_len,"omitnan"); - dc_tracking_est_s = movmean(dbg.dc_tracking_est(:),smooth_len,"omitnan"); - - figure(460); clf; - tiledlayout(2,1,"TileSpacing","compact","Padding","compact"); - - nexttile; - semilogy(error_first_s + eps,"DisplayName","first epoch"); - hold on; - semilogy(error_last_s + eps,"DisplayName","last epoch"); - grid on; - ylabel("error"); - title("FFE debug error"); - legend("Location","best"); - - nexttile;hold on - % plot(eq_sig_mov,"DisplayName","err dc eff"); - plot(dc_tracking_mu_eff_s,"DisplayName","mu DC - either fixed or adaptive"); - e_dc_scale = max(abs(dc_tracking_est_s),[],"omitnan") + eps; - plot(-1.*dc_tracking_est_s./e_dc_scale,"DisplayName","Inverted DC-tracking; Value that is subtracted during EQ"); - avg_lvl_dc=mean(avg_for_lvl,1,"omitnan"); - plot(avg_lvl_dc./0.01,"DisplayName","Smoothed EQ output signal; calc'd by showLevelScatter"); - grid on; - xlabel("Symbol"); - ylabel("mu dc eff"); - title("Effective adaptive DC step"); - legend - - inv_dc_track = -dc_tracking_est_s(:); - avg_lvl_dc = avg_lvl_dc(:); - xcorr_len = min(numel(inv_dc_track),numel(avg_lvl_dc)); - inv_dc_xcorr = inv_dc_track(1:xcorr_len); - avg_lvl_xcorr = avg_lvl_dc(1:xcorr_len); - inv_dc_xcorr = fillmissing(inv_dc_xcorr,"linear","EndValues","nearest"); - avg_lvl_xcorr = fillmissing(avg_lvl_xcorr,"linear","EndValues","nearest"); - inv_dc_xcorr = inv_dc_xcorr - mean(inv_dc_xcorr,"omitnan"); - avg_lvl_xcorr = avg_lvl_xcorr - mean(avg_lvl_xcorr,"omitnan"); - [dc_level_xcorr,dc_level_lags] = xcorr(inv_dc_xcorr,avg_lvl_xcorr,"coeff"); - [dc_level_corr,dc_level_idx] = max(dc_level_xcorr); - dc_level_delay_symbols = dc_level_lags(dc_level_idx); - output.ffe_debug.dc_level_delay_symbols = dc_level_delay_symbols; - output.ffe_debug.dc_level_delay_corr = dc_level_corr; - output.ffe_debug.dc_level_lags = dc_level_lags; - output.ffe_debug.dc_level_xcorr = dc_level_xcorr; - fprintf("FFE debug: xcorr(inv DC tracking, avg level) delay=%d symbols, corr=%6.3f\n", ... - dc_level_delay_symbols,dc_level_corr); - - figure(461); clf; - plot(dc_level_lags,dc_level_xcorr,"DisplayName","xcorr"); - hold on; - plot(dc_level_delay_symbols,dc_level_corr,"ro","DisplayName","max"); - grid on; - xlabel("Lag in symbols"); - ylabel("Correlation coefficient"); - title(sprintf("Delay estimate: %d symbols (corr %.3f)",dc_level_delay_symbols,dc_level_corr)); - legend("Location","best"); - - - end - - - - % - % - % %% MPI Reduction - % if 1 - % - % % dc_buffer_len = 0; - % % ffe_buffer_len = 0; - % % smoothing_buffer_length = options.userParameters.smoothing_length; - % % smoothing_buffer_update = 1; - % - % % eq_ffe_dcr = FFE_DCremoval_adaptive_mu(eq_settings{:}, ... - % % "dc_buffer_len",dc_buffer_len, ... - % % "ffe_buffer_len",ffe_buffer_len,... - % % "smoothing_buffer_length",smoothing_buffer_length,... - % % "smoothing_buffer_update",smoothing_buffer_update); - % - % eq_ = FFE_DCremoval("epochs_tr",5,"epochs_dd",3,"len_tr",4096*2,"mu_dd",... - % 0.0002,"mu_tr",0,"order",25,"sps",2,"decide",0,... - % "mu_dc",0.005,"dc_buffer_len",1); - % - % ffe_results_dcr = ffe(eq_, options.M, Scpe_sig, Symbols, Tx_bits, ... - % "precode_mode", options.duob_mode, ... - % 'showAnalysis', options.debug_plots, ... - % "postFFE", [], ... - % "eth_style_symbol_mapping", 0); - % - % ffe_results_dcr.metrics.print("description",'FFE DCR'); - % output.ffe_dcr_package = ffe_results_dcr; - % end end diff --git a/Functions/EQ_visuals/exportLevelScatterTikz.m b/Functions/EQ_visuals/exportLevelScatterTikz.m new file mode 100644 index 0000000..be0424f --- /dev/null +++ b/Functions/EQ_visuals/exportLevelScatterTikz.m @@ -0,0 +1,320 @@ +function exportLevelScatterTikz(sourceAx, pngFile, tikzFile, pngTikzPath, options) +%EXPORTLEVELSCATTERTIKZ Export showLevelScatter axes as PNG-backed TikZ. +% Rasterizes scatter objects into a PNG layer and writes average lines plus +% text annotations as PGFPlots objects. + +arguments + sourceAx + pngFile (1,1) string + tikzFile (1,1) string + pngTikzPath (1,1) string + options.resolutionDpi (1,1) double {mustBePositive, mustBeInteger} = 150 + options.xLabel (1,1) string = "Time in $\mu$s" + options.yLabel (1,1) string = "Normalized Amplitude" + options.generatedBy (1,1) string = "exportLevelScatterTikz.m" +end + +outputDir = fileparts(pngFile); +if ~exist(outputDir, "dir") + mkdir(outputDir); +end + +exportScatterLayerPng(sourceAx, pngFile, options.resolutionDpi); + +linePlots = collectTikzLinePlots(sourceAx); +textAnnotations = collectTikzTextAnnotations(sourceAx); +writeTikzPngAxis(tikzFile, pngTikzPath, sourceAx, linePlots, textAnnotations, options); +end + +function exportScatterLayerPng(sourceAx, pngFile, resolutionDpi) +sourceFig = ancestor(sourceAx, "figure"); +xLimits = sourceAx.XLim; +yLimits = sourceAx.YLim; +scatterHandles = findGraphicsObjectsByType(sourceAx, "scatter"); + +exportFig = figure( ... + "Visible", "off", ... + "Color", "white", ... + "Units", sourceFig.Units, ... + "Position", sourceFig.Position); +cleanupFigure = onCleanup(@() close(exportFig)); + +exportAx = copyobj(sourceAx, exportFig); +exportAx.Units = "normalized"; +exportAx.Position = [0 0 1 1]; +exportAx.XLim = xLimits; +exportAx.YLim = yLimits; +exportAx.XLimMode = "manual"; +exportAx.YLimMode = "manual"; +stripAxesToScatterOnly(exportAx); +hideAxesInkForRasterExport(exportAx); + +fprintf("Exporting scatter PNG with XLim=[%g %g], YLim=[%g %g], scatter objects=%d: %s\n", ... + xLimits(1), xLimits(2), yLimits(1), yLimits(2), numel(scatterHandles), pngFile); +drawnow; +exportFig.PaperPositionMode = "auto"; +print(exportFig, char(pngFile), "-dpng", sprintf("-r%d", resolutionDpi)); +end + +function stripAxesToScatterOnly(ax) +plotObjects = findall(ax); +for objIdx = 1:numel(plotObjects) + obj = plotObjects(objIdx); + if obj == ax || ~isvalid(obj) + continue + end + + objType = string(get(obj, "Type")); + if lower(objType) ~= "scatter" + delete(obj); + end +end +end + +function hideAxesInkForRasterExport(ax) +ax.Visible = "on"; +ax.Color = "white"; +ax.Box = "off"; +ax.XGrid = "off"; +ax.YGrid = "off"; +ax.XMinorGrid = "off"; +ax.YMinorGrid = "off"; +ax.XTick = []; +ax.YTick = []; +ax.XColor = "white"; +ax.YColor = "white"; +ax.Title.String = ""; +ax.XLabel.String = ""; +ax.YLabel.String = ""; +end + +function handles = findGraphicsObjectsByType(parentHandle, objectType) +allHandles = findall(parentHandle); +matches = false(size(allHandles)); +for handleIdx = 1:numel(allHandles) + currentHandle = allHandles(handleIdx); + if ~isvalid(currentHandle) + continue + end + + matches(handleIdx) = lower(string(get(currentHandle, "Type"))) == lower(string(objectType)); +end + +handles = allHandles(matches); +end + +function linePlots = collectTikzLinePlots(sourceAx) +lineHandles = findGraphicsObjectsByType(sourceAx, "line"); +lineHandles = flipud(lineHandles(:)); + +linePlots = repmat(struct( ... + "xData", [], ... + "yData", [], ... + "color", [], ... + "lineWidth", []), numel(lineHandles), 1); +validLineCount = 0; + +for lineIdx = 1:numel(lineHandles) + lineHandle = lineHandles(lineIdx); + xData = lineHandle.XData(:); + yData = lineHandle.YData(:); + validSamples = isfinite(xData) & isfinite(yData); + + if ~any(validSamples) + continue + end + + validLineCount = validLineCount + 1; + linePlots(validLineCount).xData = xData(validSamples); + linePlots(validLineCount).yData = yData(validSamples); + linePlots(validLineCount).color = lineHandle.Color; + linePlots(validLineCount).lineWidth = lineHandle.LineWidth; +end + +linePlots = linePlots(1:validLineCount); +end + +function textAnnotations = collectTikzTextAnnotations(sourceAx) +textHandles = findGraphicsObjectsByType(sourceAx, "text"); +textHandles = flipud(textHandles(:)); +xLimits = sourceAx.XLim; +yLimits = sourceAx.YLim; + +textAnnotations = repmat(struct( ... + "x", [], ... + "y", [], ... + "text", "", ... + "anchor", "center"), numel(textHandles), 1); +validTextCount = 0; + +for textIdx = 1:numel(textHandles) + textHandle = textHandles(textIdx); + labelText = string(textHandle.String); + if strlength(strtrim(labelText)) == 0 + continue + end + + textPosition = textHandle.Position; + if textPosition(1) < xLimits(1) || textPosition(1) > xLimits(2) || ... + textPosition(2) < yLimits(1) || textPosition(2) > yLimits(2) + continue + end + + validTextCount = validTextCount + 1; + textAnnotations(validTextCount).x = textPosition(1); + textAnnotations(validTextCount).y = textPosition(2); + textAnnotations(validTextCount).text = labelText; + textAnnotations(validTextCount).anchor = matlabTextAlignmentToTikzAnchor( ... + textHandle.HorizontalAlignment, textHandle.VerticalAlignment); +end + +textAnnotations = textAnnotations(1:validTextCount); +end + +function writeTikzPngAxis(tikzFile, pngTikzPath, sourceAx, linePlots, textAnnotations, options) +outputDir = fileparts(tikzFile); +if ~exist(outputDir, "dir") + mkdir(outputDir); +end + +xLimits = sourceAx.XLim; +yLimits = sourceAx.YLim; +fid = fopen(tikzFile, "w"); +if fid < 0 + error("exportLevelScatterTikz:TikzOpenFailed", ... + "Could not open TikZ export file: %s", tikzFile); +end +cleanupFile = onCleanup(@() fclose(fid)); + +axisOptions = { + " every axis/.append style={font=\scriptsize}," + "width=\fwidth," + "height=\fheight," + "at={(0\fwidth,0\fheight)}," + "scale only axis," + "axis on top," + "xmin=" + formatPgfNumber(xLimits(1)) + "," + "xmax=" + formatPgfNumber(xLimits(2)) + "," + "xlabel style={font=\color{white!15!black}\small}," + "xlabel={" + options.xLabel + "}," + "ymin=" + formatPgfNumber(yLimits(1)) + "," + "ymax=" + formatPgfNumber(yLimits(2)) + "," + "ylabel style={font=\color{white!15!black}\small}," + "ylabel={" + options.yLabel + "}," + "axis background/.style={fill=white}," + "xmajorgrids," + "ymajorgrids," + "grid style={line width=0.4pt, dotted, color=black!20}," + "scaled ticks=false," + "tick label style={/pgf/number format/fixed, /pgf/number format/1000 sep={}}," + "legend columns=1" + }; + +fprintf(fid, "%% This file was generated by %s.\n", options.generatedBy); +fprintf(fid, "%%\n"); +for lineIdx = 1:numel(linePlots) + fprintf(fid, "%s\n", formatTikzColorDefinition(lineIdx, linePlots(lineIdx).color)); +end +if ~isempty(linePlots) + fprintf(fid, "%%\n"); +end +fprintf(fid, "%s\n\n", "\begin{tikzpicture}"); +fprintf(fid, "%s\n", "\begin{axis}[%"); +for optionIdx = 1:numel(axisOptions) + fprintf(fid, "%s\n", axisOptions{optionIdx}); +end +fprintf(fid, "%s\n", "]"); + +graphicsLine = sprintf( ... + "\\addplot [forget plot] graphics [xmin=%s, xmax=%s, ymin=%s, ymax=%s] {%s};", ... + formatPgfNumber(xLimits(1)), ... + formatPgfNumber(xLimits(2)), ... + formatPgfNumber(yLimits(1)), ... + formatPgfNumber(yLimits(2)), ... + char(strrep(pngTikzPath, "\", "/"))); +fprintf(fid, "%s\n\n", graphicsLine); +writeTikzLinePlots(fid, linePlots); +writeTikzTextAnnotations(fid, textAnnotations); +fprintf(fid, "%s\n\n", "\end{axis}"); +fprintf(fid, "%s", "\end{tikzpicture}%"); +end + +function valueText = formatPgfNumber(value) +valueText = string(sprintf("%.15g", value)); +end + +function colorDefinition = formatTikzColorDefinition(colorIdx, rgbColor) +colorDefinition = sprintf( ... + "\\definecolor{mycolor%d}{rgb}{%.5f,%.5f,%.5f}%%", ... + colorIdx, rgbColor(1), rgbColor(2), rgbColor(3)); +end + +function writeTikzLinePlots(fid, linePlots) +for lineIdx = 1:numel(linePlots) + fprintf(fid, "\\addplot [color=mycolor%d, line width=%.1fpt, forget plot]\n", ... + lineIdx, linePlots(lineIdx).lineWidth); + fprintf(fid, " table[row sep=crcr]{%%\n"); + + for sampleIdx = 1:numel(linePlots(lineIdx).xData) + fprintf(fid, "%s\t%s\\\\\n", ... + formatPgfNumber(linePlots(lineIdx).xData(sampleIdx)), ... + formatPgfNumber(linePlots(lineIdx).yData(sampleIdx))); + end + + fprintf(fid, "};\n\n"); +end +end + +function writeTikzTextAnnotations(fid, textAnnotations) +for textIdx = 1:numel(textAnnotations) + fprintf(fid, "\\node[font=\\scriptsize, anchor=%s, fill=white, draw=black!40, rounded corners=2pt, inner sep=2pt]\n", ... + textAnnotations(textIdx).anchor); + fprintf(fid, " at (axis cs:%s,%s) {%s};\n\n", ... + formatPgfNumber(textAnnotations(textIdx).x), ... + formatPgfNumber(textAnnotations(textIdx).y), ... + escapeTikzText(textAnnotations(textIdx).text)); +end +end + +function anchor = matlabTextAlignmentToTikzAnchor(horizontalAlignment, verticalAlignment) +horizontalAlignment = string(horizontalAlignment); +verticalAlignment = string(verticalAlignment); + +if horizontalAlignment == "left" + horizontalAnchor = "west"; +elseif horizontalAlignment == "right" + horizontalAnchor = "east"; +else + horizontalAnchor = ""; +end + +if verticalAlignment == "top" + verticalAnchor = "north"; +elseif verticalAlignment == "bottom" + verticalAnchor = "south"; +else + verticalAnchor = ""; +end + +if strlength(horizontalAnchor) > 0 && strlength(verticalAnchor) > 0 + anchor = verticalAnchor + " " + horizontalAnchor; +elseif strlength(horizontalAnchor) > 0 + anchor = horizontalAnchor; +elseif strlength(verticalAnchor) > 0 + anchor = verticalAnchor; +else + anchor = "center"; +end +end + +function textOut = escapeTikzText(textIn) +textOut = char(textIn); +% textOut = strrep(textOut, "\", "\textbackslash{}"); +textOut = strrep(textOut, "{", "\{"); +textOut = strrep(textOut, "}", "\}"); +textOut = strrep(textOut, "%", "\%"); +textOut = strrep(textOut, "&", "\&"); +textOut = strrep(textOut, "#", "\#"); +textOut = strrep(textOut, "_", "\_"); +% textOut = strrep(textOut, "$", "\$"); +end diff --git a/Tests/04_DSP/Equalizer/FFE_test.m b/Tests/04_DSP/Equalizer/FFE_test.m index 3d9896d..7bdfc6f 100644 --- a/Tests/04_DSP/Equalizer/FFE_test.m +++ b/Tests/04_DSP/Equalizer/FFE_test.m @@ -14,10 +14,13 @@ classdef FFE_test < IMDDTestCase testCase.verifyEqual(ffe.epochs_dd, 5); testCase.verifyFalse(logical(ffe.optimize_dc_tracking_params)); testCase.verifyEqual(ffe.dc_tracking_optimization_delay_weight, 1e-3); + testCase.verifyEqual(ffe.dc_tracking_persistence_gain, 0); testCase.verifyEqual(ffe.dc_tracking_power_exponent, 2); testCase.verifyEqual(ffe.dc_avg_bufferlength_a1, 0); + testCase.verifyEqual(ffe.dc_avg_update_blocklength_a1, 0); testCase.verifyEqual(ffe.dc_smoothing_a1, 0); testCase.verifyEqual(ffe.dc_level_avg_bufferlength_a2, 0); + testCase.verifyEqual(ffe.dc_level_update_blocklength_a2, 0); testCase.verifyEqual(ffe.dc_smoothing_a2, 0); testCase.verifyEqual(ffe.dc_level_weights_a2, 0); testCase.verifyEqual(ffe.dc_tracking_buffer_len, 1); @@ -101,6 +104,76 @@ classdef FFE_test < IMDDTestCase testCase.verifyEqual(ffe.e_dc, 0.19, "AbsTol", 1e-12); end + function dcTrackingUnbufferedPathUpdatesEverySymbol(testCase) + x = zeros(4, 1); + d = ones(4, 1); + + ffe = FFE( ... + "sps", 1, ... + "order", 1, ... + "dc_tracking_mu", 0.1, ... + "dc_tracking_buffer_len", 1, ... + "dd_mode", 0, ... + "adaption_technique", adaption_method.lms); + + ffe.constellation = unique(d); + ffe.equalize(x, d, 0, 1, numel(x), true, false); + + testCase.verifyEqual(ffe.e_dc, 1 - 0.9^4, "AbsTol", 1e-12); + end + + function dcTrackingBlockUpdateUsesMeanBlockError(testCase) + ffe = FFE("dc_tracking_mu", 0.1); + + [e_dc_next,stats] = ffe.dcTrackingBlockUpdate(0.5, [1; 2; 3]); + + testCase.verifyEqual(e_dc_next, 0.7, "AbsTol", 1e-12); + testCase.verifyEqual(stats.err_mean, 2, "AbsTol", 1e-12); + testCase.verifyEqual(stats.mu_eff, 0.1, "AbsTol", 1e-12); + testCase.verifyEqual(stats.update, 0.2, "AbsTol", 1e-12); + end + + function dcTrackingBlockUpdateDoesNotSuppressLargeErrorPower(testCase) + ffe = FFE("dc_tracking_mu", 0.1); + + [e_dc_next,stats] = ffe.dcTrackingBlockUpdate(0, [10; 10]); + + testCase.verifyEqual(e_dc_next, 1, "AbsTol", 1e-12); + testCase.verifyEqual(stats.err_mean, 10, "AbsTol", 1e-12); + testCase.verifyEqual(stats.mu_eff, 0.1, "AbsTol", 1e-12); + end + + function dcTrackingBlockUpdateCanBoostPersistentMean(testCase) + ffe = FFE("dc_tracking_mu", 0.1); + + [e_dc_next,stats] = ffe.dcTrackingBlockUpdate(0, [2; 2; 2], ... + "persistence_gain", 1); + + testCase.verifyEqual(e_dc_next, 0.4, "AbsTol", 1e-12); + testCase.verifyEqual(stats.persistence_scale, 1, "AbsTol", 1e-12); + testCase.verifyEqual(stats.mu_eff, 0.2, "AbsTol", 1e-12); + end + + function dcTrackingLoopUsesBlockPersistenceGain(testCase) + x = zeros(4, 1); + d = ones(4, 1); + + ffe = FFE( ... + "sps", 1, ... + "order", 1, ... + "dc_tracking_mu", 0.1, ... + "dc_tracking_adaptive_enabled", true, ... + "dc_tracking_persistence_gain", 1, ... + "dc_tracking_buffer_len", 2, ... + "dd_mode", 0, ... + "adaption_technique", adaption_method.lms); + + ffe.constellation = unique(d); + ffe.equalize(x, d, 0, 1, numel(x), true, false); + + testCase.verifyEqual(ffe.e_dc, 0.36, "AbsTol", 1e-12); + end + function ffeBufferDelaysTapUpdateUntilBlockBoundary(testCase) x = ones(4, 1); d = ones(4, 1); @@ -181,6 +254,28 @@ classdef FFE_test < IMDDTestCase testCase.verifyEqual(ffe.debug_struct.dc_avg_offset, expected_offset, "AbsTol", 1e-12); end + function dcAvgA1UpdateBlocklengthOneRefreshesEverySymbol(testCase) + x = (1:4).'; + d = zeros(4, 1); + expected_offset = [0, 1, 1.5, 2]; + expected_y = [1; 1; 1.5; 2]; + + ffe = FFE( ... + "sps", 1, ... + "order", 1, ... + "dc_avg_bufferlength_a1", 3, ... + "dc_avg_update_blocklength_a1", 1, ... + "save_debug", true, ... + "dd_mode", 0, ... + "adaption_technique", adaption_method.lms); + ffe.e = 1; + + [y,~] = ffe.equalize(x, d, 0, 1, numel(x), true, false); + + testCase.verifyEqual(y, expected_y, "AbsTol", 1e-12); + testCase.verifyEqual(ffe.debug_struct.dc_avg_offset, expected_offset, "AbsTol", 1e-12); + end + function dcAvgA1WarnsWhenCombinedWithDcTracking(testCase) testCase.verifyWarning(@() FFE( ... "dc_avg_bufferlength_a1", 3, ... @@ -211,6 +306,30 @@ classdef FFE_test < IMDDTestCase testCase.verifyEqual(ffe.debug_struct.dc_level_valid_count, [0, 0, 2, 2]); end + function dcLevelAvgA2UpdateBlocklengthOneRefreshesEverySymbol(testCase) + x = (1:4).'; + d = zeros(4, 1); + + ffe = FFE( ... + "sps", 1, ... + "order", 1, ... + "dc_level_avg_bufferlength_a2", 2, ... + "dc_level_update_blocklength_a2", 1, ... + "dc_level_weights_a2", 1, ... + "save_debug", true, ... + "dd_mode", 0, ... + "adaption_technique", adaption_method.lms); + ffe.e = 1; + ffe.constellation = 0; + + [y,~] = ffe.equalize(x, d, 0, 1, numel(x), true, false); + + testCase.verifyEqual(y, [1; 1.5; 1.5; 1.5], "AbsTol", 1e-12); + testCase.verifyEqual(ffe.debug_struct.dc_level_mpi_est, [0, 1, 1.5, 2.5], "AbsTol", 1e-12); + testCase.verifyEqual(ffe.debug_struct.dc_level_weight, [0, 0.5, 1, 1], "AbsTol", 1e-12); + testCase.verifyEqual(ffe.debug_struct.dc_level_valid_count, [0, 1, 2, 2]); + end + function dcLevelAvgA2UsesVectorWeightsInConstellationOrder(testCase) x = [1; 3; 1; 4]; d = [0; 2; 0; 2]; diff --git a/before_eq_-1.png b/before_eq_-1.png new file mode 100644 index 0000000..5b0ec62 Binary files /dev/null and b/before_eq_-1.png differ diff --git a/projects/Diss/MPI_revisit/algorithms/PLOT_ber_vs_sir.m b/projects/Diss/MPI_revisit/algorithms/PLOT_ber_vs_sir.m new file mode 100644 index 0000000..d200677 --- /dev/null +++ b/projects/Diss/MPI_revisit/algorithms/PLOT_ber_vs_sir.m @@ -0,0 +1,237 @@ +%% Simple BER over SIR plot from the recombined warehouses +% This script intentionally keeps the mechanics visible: +% 1) load the run_id warehouse and the grouped config warehouse +% 2) loop over occupied config cells +% 3) read the physical coordinates of the cell +% 4) extract all BER values stored in that cell +% 5) plot individual BER dots and one larger mean marker per SIR + +clear; clc; + +scriptDir = fileparts(mfilename("fullpath")); + +runWhFile = fullfile(scriptDir, "combined_by_run_id_results.mat"); +configWhFile = fullfile(scriptDir, "combined_by_config_results.mat"); + +runWhData = load(runWhFile); +configWhData = load(configWhFile); + +wh_run_id_combined = runWhData.wh_run_id_combined; +wh_config_combined = configWhData.wh_config_combined; +algorithmStorageNames = configWhData.algorithmStorageNames; + +fprintf("Loaded run_id warehouse:\n %s\n", runWhFile); +wh_run_id_combined.showInfo; + +fprintf("\nLoaded config-grouped warehouse:\n %s\n", configWhFile); +wh_config_combined.showInfo; + +%% Plot settings + +pathLengthToPlot = 1; % use 300 to see duplicate run_ids per config; set 1000 for the previous path +selectedPamLevels = 4;%wh_config_combined.parameter.pam_level.values; +selectedAlgorithms = algorithmStorageNames([1,5]); + +useBoundedLines = true; % switch uncertainty bands on/off here +usePolyfit = true; % switch fitted dashed trend lines on/off here +polyfitOrderMax = 4; + + +% SIR vs BER: one subplot per PAM level, all algorithms compared + +if exist("linspecer", "file") + algColors = linspecer(numel(selectedAlgorithms)); +else + algColors = lines(numel(selectedAlgorithms)); +end + +figure(); clf; +tiledlayout(numel(selectedPamLevels), 1, "TileSpacing", "compact"); + +for pamIdx = 1:numel(selectedPamLevels) + pamLevel = selectedPamLevels(pamIdx); + nexttile; hold on; + + for algIdx = 1:numel(selectedAlgorithms) + algorithmName = selectedAlgorithms{algIdx}; + algColor = algColors(algIdx, :); + + berTable = buildBerTable(wh_config_combined, algorithmName, pathLengthToPlot); + berTable = berTable(berTable.pam_level == pamLevel, :); + + % berTable(berTable.sir == 23,:) = []; + + if isempty(berTable) + continue + end + + sirValues = unique(berTable.sir(:).'); + sirValues = sort(sirValues); + + meanBer = nan(size(sirValues)); + minBer = nan(size(sirValues)); + maxBer = nan(size(sirValues)); + + for sirIdx = 1:numel(sirValues) + sirValue = sirValues(sirIdx); + sirRows = berTable(berTable.sir == sirValue, :); + berValues = [sirRows.ber_values{:}]; + berValues = berValues(isfinite(berValues)); + berValues = berValues(berValues<0.1); + berValues = rmoutliers(berValues); + + + if isempty(berValues) + continue + end + + scatter(repmat(sirValue, size(berValues)), berValues, ... + 30, ... + "Marker", ".", ... + "MarkerEdgeColor", algColor, ... + "MarkerFaceColor", algColor, ... + "HandleVisibility", "off"); + + meanBer(sirIdx) = mean(berValues, "omitnan"); + minBer(sirIdx) = min(berValues); + maxBer(sirIdx) = max(berValues); + end + + valid = isfinite(sirValues) & isfinite(meanBer); + if ~any(valid) + continue + end + + if useBoundedLines && exist("boundedline", "file") + yLower = max(meanBer - minBer, 0); + yUpper = max(maxBer - meanBer, 0); + yBounds = [yLower(:), yUpper(:)]; + + [hl, hp] = boundedline(sirValues(valid).', meanBer(valid).', yBounds(valid, :), ... + 'alpha', 'transparency', 0.08, ... + 'cmap', algColor, ... + 'nan', 'fill', ... + 'orientation', 'vert'); + set(hl, "LineStyle", "none", "Marker", "none", "HandleVisibility", "off"); + set(hp, "LineStyle", "none", "HandleVisibility", "off"); + end + + scatter(sirValues(valid), meanBer(valid), ... + 20, ... + "Marker", "o", ... + "MarkerEdgeColor", algColor, ... + "MarkerFaceColor", algColor, ... + "LineWidth", 1, ... + "DisplayName", algorithmName); + + if usePolyfit + fitMask = valid & meanBer > 0; + if nnz(fitMask) >= 2 + fitOrder = min(polyfitOrderMax, nnz(fitMask) - 1); + fitCoeff = polyfit(sirValues(fitMask), log10(meanBer(fitMask)), fitOrder); + xFit = linspace(min(sirValues(fitMask)), max(sirValues(fitMask)), 300); + yFit = 10 .^ polyval(fitCoeff, xFit); + + plot(xFit, yFit, ... + "LineStyle", "--", ... + "LineWidth", 1.1, ... + "Color", algColor, ... + "HandleVisibility", "off"); + end + end + end + + title(sprintf("PAM %.0f, path %.0f m", pamLevel, pathLengthToPlot)); + xlabel("SIR (dB)"); + ylabel("BER"); + set(gca, "YScale", "log"); + xlim([15,45]); + grid on; + box on; + legend("Location", "best", "Interpreter", "none"); +end + +%% Local helpers + +function berTable = buildBerTable(wh, storageName, pathLength) +rows = struct( ... + "pam_level", {}, ... + "symbolrate", {}, ... + "interference_path_length", {}, ... + "sir", {}, ... + "block_update", {}, ... + "run_ids", {}, ... + "n_run_ids", {}, ... + "ber_values", {}, ... + "n_ber_values", {}); + +for linIdx = 1:numel(wh.sto.(storageName)) + storedValue = wh.sto.(storageName){linIdx}; + if isempty(storedValue) + continue + end + + phys = linearIndexToStruct(wh, linIdx); + if phys.interference_path_length ~= pathLength + continue + end + + berValues = extractBerValues(storedValue); + if isempty(berValues) + continue + end + + runIds = wh.sto.run_id{linIdx}; + if isempty(runIds) + runIds = NaN; + end + + row = struct(); + row.pam_level = phys.pam_level; + row.symbolrate = phys.symbolrate; + row.interference_path_length = phys.interference_path_length; + row.sir = phys.sir; + row.block_update = phys.block_update; + row.run_ids = {runIds}; + row.n_run_ids = numel(runIds); + row.ber_values = {berValues}; + row.n_ber_values = numel(berValues); + + rows(end+1) = row; %#ok +end + +if isempty(rows) + berTable = struct2table(rows); +else + berTable = struct2table(rows); + berTable = sortrows(berTable, ["pam_level", "sir", "symbolrate"]); +end +end + +function physStruct = linearIndexToStruct(wh, linIdx) +[physValues, physNames] = wh.getPhysIndicesByLinIndex(linIdx); +physStruct = struct(); + +for paramIdx = 1:numel(physNames) + physStruct.(char(physNames{paramIdx})) = physValues{paramIdx}; +end +end + +function berValues = extractBerValues(value) +berValues = []; + +if isstruct(value) + if isfield(value, "metrics") && hasBer(value.metrics) + berValues(end+1) = value.metrics.BER; + end +elseif iscell(value) + for valueIdx = 1:numel(value) + berValues = [berValues, extractBerValues(value{valueIdx})]; %#ok + end +end +end + +function tf = hasBer(metrics) +tf = (isstruct(metrics) && isfield(metrics, "BER")) || ... + (isobject(metrics) && isprop(metrics, "BER")); +end diff --git a/projects/Diss/MPI_revisit/algorithms/recombine_warehouses.m b/projects/Diss/MPI_revisit/algorithms/recombine_warehouses.m new file mode 100644 index 0000000..85b13b0 --- /dev/null +++ b/projects/Diss/MPI_revisit/algorithms/recombine_warehouses.m @@ -0,0 +1,539 @@ +% === DSP settings === +dsp_options = struct(); +dsp_options.mode = "run_id"; +dsp_options.recipe = @mpi_recipe_dev; +dsp_options.append_to_db = false; +dsp_options.start_occurence = 1; +dsp_options.max_occurences = 15; +dsp_options.debug_plots = false; + +dsp_options.database_type = "mysql"; + +dsp_options.dataBase = "labor"; +dsp_options.storage_path = "W:\labdata\ECOC Silas\ecoc_2025"; + +dsp_options.server = "192.168.178.192"; +dsp_options.port = 3306; +dsp_options.user = "silas"; +dsp_options.password = "silas"; + +db = DBHandler("dataBase", [dsp_options.dataBase], ... + "type", dsp_options.database_type, ... + "server", dsp_options.server, ... + "user", dsp_options.user, "password", dsp_options.password); + +scriptDir = fileparts(mfilename("fullpath")); + +whpam4_loaded = load(fullfile(scriptDir, "pam_4_results.mat")); +whpam6_loaded = load(fullfile(scriptDir, "pam_6_results.mat")); +whpam8_loaded = load(fullfile(scriptDir, "pam_8_results.mat")); + +whpam4 = whpam4_loaded.obj; +whpam6 = whpam6_loaded.obj; +whpam8 = whpam8_loaded.obj; + +whList = {whpam4, whpam6, whpam8}; +algorithmStorageNames = fieldnames(whpam4.sto); + +%% Build a lossless run_id warehouse and cache run metadata + +pamformats = [4, 6, 8]; +baudrates = [112e9, 96e9, 72e9]; +runMetaTable = queryRunMetadata(db, whList, pamformats, baudrates); + +wh_run_id_combined = mergeDataStoragesByUnion(whList); +wh_run_id_combined = addRunMetadataStorages(wh_run_id_combined, runMetaTable); + +runIdSavePath = fullfile(scriptDir, "combined_by_run_id_results.mat"); +save(runIdSavePath, "wh_run_id_combined", "runMetaTable", "algorithmStorageNames", '-v7.3'); + +fprintf("Saved run_id-combined warehouse:\n %s\n", runIdSavePath); +wh_run_id_combined.showInfo; + +%% Build a grouped config warehouse for analysis views + +configFields = ["pam_level", "symbolrate", "interference_path_length", "sir"]; +wh_config_combined = buildConfigWarehouse( ... + wh_run_id_combined, runMetaTable, algorithmStorageNames, configFields); + +configSavePath = fullfile(scriptDir, "combined_by_config_results.mat"); +save(configSavePath, "wh_config_combined", "runMetaTable", "algorithmStorageNames", '-v7.3'); + +fprintf("Saved config-grouped warehouse:\n %s\n", configSavePath); +wh_config_combined.showInfo; + +%% BER over SIR at fixed path length + +plot_options = struct(); +plot_options.path_length = 1000; +plot_options.use_boundedlines = true; +plot_options.boundedline_alpha = 0.10; +plot_options.polyfit_order_max = 3; +plot_options.fec_ber_threshold = 2e-2; +plot_options.figure_base = 7100; +plot_options.xlim = [15 35]; + +plotBerOverSirForPath(wh_config_combined, algorithmStorageNames, plot_options); + +%% Local functions + +function runMetaTable = queryRunMetadata(db, whList, pamformats, baudrates) +metaTables = {}; +metadataFields = db.getTableFieldNames('Runs'); + +for whIdx = 1:numel(whList) + curWh = whList{whIdx}; + curRunIds = curWh.parameter.run_id.values(:); + + fp = QueryFilter(); + fp.where('Runs', 'run_id', 'GREATER_EQUAL', 3153); + fp.where('Runs', 'symbolrate', 'EQUALS', baudrates(whIdx)); + fp.where('Runs', 'fiber_length', 'EQUALS', 0); + fp.where('Runs', 'db_mode', 'EQUALS', '"no_db"'); + fp.where('Runs', 'pam_level', 'EQUALS', pamformats(whIdx)); + + [dataTable, ~] = db.queryDB(fp, metadataFields); + dataTable.run_id = numericColumn(dataTable.run_id); + dataTable = dataTable(ismember(dataTable.run_id, curRunIds), :); + + missingRunIds = setdiff(curRunIds, dataTable.run_id); + if ~isempty(missingRunIds) + warning("recombine_warehouses:MissingMetadata", ... + "%d run_id(s) from PAM %.0f warehouse were not found in the metadata query.", ... + numel(missingRunIds), pamformats(whIdx)); + end + + metaTables{end+1} = dataTable; %#ok +end + +runMetaTable = vertcat(metaTables{:}); +[~, uniqueIdx] = unique(numericColumn(runMetaTable.run_id), "stable"); +runMetaTable = runMetaTable(uniqueIdx, :); + +numericFields = ["run_id", "pam_level", "symbolrate", "interference_path_length", "sir"]; +for fieldIdx = 1:numel(numericFields) + fieldName = numericFields(fieldIdx); + if ismember(fieldName, string(runMetaTable.Properties.VariableNames)) + runMetaTable.(char(fieldName)) = numericColumn(runMetaTable.(char(fieldName))); + end +end + +[~, sortIdx] = sort(runMetaTable.run_id); +runMetaTable = runMetaTable(sortIdx, :); +end + +function x = numericColumn(x) +if isnumeric(x) + x = double(x); +elseif iscell(x) + x = str2double(string(x)); +elseif isstring(x) || ischar(x) || iscategorical(x) + x = str2double(string(x)); +else + x = double(x); +end +x = x(:); +end + +function whMerged = mergeDataStoragesByUnion(whList) +templateWh = whList{1}; +paramNames = cellstr(templateWh.fn); + +for whIdx = 2:numel(whList) + if ~isequal(sort(cellstr(whList{whIdx}.fn)), sort(paramNames)) + error("recombine_warehouses:ParameterMismatch", ... + "All input warehouses must use the same parameter names."); + end +end + +mergedParams = struct(); +for paramIdx = 1:numel(paramNames) + paramName = paramNames{paramIdx}; + values = []; + + for whIdx = 1:numel(whList) + newValues = whList{whIdx}.parameter.(paramName).values(:).'; + values = [values, newValues]; %#ok + end + + values = unique(values, "stable"); + if isnumeric(values) + values = sort(values); + end + mergedParams.(paramName) = values; +end + +whMerged = DataStorage(mergedParams); + +for whIdx = 1:numel(whList) + curWh = whList{whIdx}; + curStorageNames = fieldnames(curWh.sto); + + for storageIdx = 1:numel(curStorageNames) + storageName = curStorageNames{storageIdx}; + if ~isfield(whMerged.sto, storageName) + whMerged.addStorage(storageName); + end + + for linIdx = 1:numel(curWh.sto.(storageName)) + storedValue = curWh.sto.(storageName){linIdx}; + if isempty(storedValue) + continue + end + + targetLinIdx = mapLinearIndex(curWh, whMerged, linIdx); + existingValue = whMerged.sto.(storageName){targetLinIdx}; + whMerged.sto.(storageName){targetLinIdx} = mergeStoredValue(existingValue, storedValue); + end + end +end +end + +function wh = addRunMetadataStorages(wh, runMetaTable) +metadataStorageNames = ["meta_pam_level", "meta_symbolrate", ... + "meta_interference_path_length", "meta_sir"]; +metadataTableFields = ["pam_level", "symbolrate", "interference_path_length", "sir"]; + +for fieldIdx = 1:numel(metadataStorageNames) + if ~isfield(wh.sto, metadataStorageNames(fieldIdx)) + wh.addStorage(char(metadataStorageNames(fieldIdx))); + end +end + +for linIdx = 1:wh.getLastLinIndice() + physStruct = linearIndexToStruct(wh, linIdx); + if ~isfield(physStruct, "run_id") + continue + end + + metaIdx = find(runMetaTable.run_id == physStruct.run_id, 1); + if isempty(metaIdx) + continue + end + + for fieldIdx = 1:numel(metadataStorageNames) + wh.sto.(char(metadataStorageNames(fieldIdx))){linIdx} = ... + runMetaTable.(char(metadataTableFields(fieldIdx)))(metaIdx); + end +end +end + +function whConfig = buildConfigWarehouse(whRun, runMetaTable, algorithmStorageNames, configFields) +runParamNames = cellstr(whRun.fn); +extraParamNames = setdiff(string(runParamNames), "run_id", "stable"); + +configParams = struct(); +for fieldIdx = 1:numel(configFields) + fieldName = configFields(fieldIdx); + values = unique(runMetaTable.(char(fieldName))(:).', "stable"); + if isnumeric(values) + values = sort(values); + end + configParams.(char(fieldName)) = values; +end + +for paramIdx = 1:numel(extraParamNames) + paramName = extraParamNames(paramIdx); + configParams.(char(paramName)) = whRun.parameter.(char(paramName)).values; +end + +whConfig = DataStorage(configParams); +for storageIdx = 1:numel(algorithmStorageNames) + whConfig.addStorage(algorithmStorageNames{storageIdx}); +end +whConfig.addStorage("run_id"); + +for storageIdx = 1:numel(algorithmStorageNames) + storageName = algorithmStorageNames{storageIdx}; + + for linIdx = 1:numel(whRun.sto.(storageName)) + storedValue = whRun.sto.(storageName){linIdx}; + if isempty(storedValue) + continue + end + + sourcePhys = linearIndexToStruct(whRun, linIdx); + metaIdx = find(runMetaTable.run_id == sourcePhys.run_id, 1); + if isempty(metaIdx) + continue + end + + targetArgs = configTargetArgs(whConfig, sourcePhys, runMetaTable, metaIdx, configFields); + targetLinIdx = whConfig.getIndicesByPhys(targetArgs); + + existingValue = whConfig.sto.(storageName){targetLinIdx}; + whConfig.sto.(storageName){targetLinIdx} = mergeStoredValue(existingValue, storedValue); + end +end + +for linIdx = 1:whRun.getLastLinIndice() + sourcePhys = linearIndexToStruct(whRun, linIdx); + metaIdx = find(runMetaTable.run_id == sourcePhys.run_id, 1); + if isempty(metaIdx) + continue + end + + targetArgs = configTargetArgs(whConfig, sourcePhys, runMetaTable, metaIdx, configFields); + targetLinIdx = whConfig.getIndicesByPhys(targetArgs); + whConfig.sto.run_id{targetLinIdx} = mergeUniqueNumeric( ... + whConfig.sto.run_id{targetLinIdx}, sourcePhys.run_id); +end +end + +function targetArgs = configTargetArgs(whConfig, sourcePhys, runMetaTable, metaIdx, configFields) +targetArgs = cell(1, numel(whConfig.fn)); + +for paramIdx = 1:numel(whConfig.fn) + paramName = string(whConfig.fn(paramIdx)); + if any(configFields == paramName) + targetArgs{paramIdx} = runMetaTable.(char(paramName))(metaIdx); + else + targetArgs{paramIdx} = sourcePhys.(char(paramName)); + end +end +end + +function targetLinIdx = mapLinearIndex(sourceWh, targetWh, sourceLinIdx) +sourcePhys = linearIndexToStruct(sourceWh, sourceLinIdx); +targetArgs = cell(1, numel(targetWh.fn)); + +for paramIdx = 1:numel(targetWh.fn) + paramName = char(targetWh.fn(paramIdx)); + targetArgs{paramIdx} = sourcePhys.(paramName); +end + +targetLinIdx = targetWh.getIndicesByPhys(targetArgs); +end + +function physStruct = linearIndexToStruct(wh, linIdx) +[physValues, physNames] = wh.getPhysIndicesByLinIndex(linIdx); +physStruct = struct(); + +for paramIdx = 1:numel(physNames) + physStruct.(char(physNames{paramIdx})) = physValues{paramIdx}; +end +end + +function out = mergeStoredValue(existingValue, newValue) +if isempty(existingValue) + out = newValue; + return +end + +existingPackages = asPackageCell(existingValue); +newPackages = asPackageCell(newValue); +out = [existingPackages(:); newPackages(:)].'; +end + +function packages = asPackageCell(value) +if isempty(value) + packages = {}; +elseif iscell(value) + packages = value(:).'; +else + packages = {value}; +end +end + +function out = mergeUniqueNumeric(existingValue, newValue) +if isempty(existingValue) + out = newValue; +else + out = unique([existingValue(:); newValue(:)].', "stable"); +end +end + +function plotBerOverSirForPath(whConfig, algorithmStorageNames, plotOptions) +paramNames = string(whConfig.fn); + +if ~any(paramNames == "pam_level") || ~any(paramNames == "interference_path_length") || ~any(paramNames == "sir") + error("recombine_warehouses:MissingPlotParameters", ... + "The config warehouse must contain pam_level, interference_path_length, and sir parameters."); +end + +pamLevels = whConfig.parameter.pam_level.values; +hasBlockUpdate = any(paramNames == "block_update"); +if hasBlockUpdate + blockUpdates = whConfig.parameter.block_update.values; +else + blockUpdates = NaN; +end + +if exist("linspecer", "file") + curveColors = linspecer(max(numel(algorithmStorageNames), 1)); +else + curveColors = lines(max(numel(algorithmStorageNames), 1)); +end + +figureIdx = plotOptions.figure_base; +for blockIdx = 1:numel(blockUpdates) + for pamIdx = 1:numel(pamLevels) + pamLevel = pamLevels(pamIdx); + figureIdx = figureIdx + 1; + figure(figureIdx); clf; hold on + + for storageIdx = 1:numel(algorithmStorageNames) + storageName = algorithmStorageNames{storageIdx}; + curveColor = curveColors(storageIdx, :); + + filter = struct(); + filter.pam_level = pamLevel; + filter.interference_path_length = plotOptions.path_length; + if hasBlockUpdate + filter.block_update = blockUpdates(blockIdx); + end + + [sirValues, berGroups] = collectBerBySir(whConfig, storageName, filter); + if isempty(sirValues) + continue + end + + [sirValues, sortIdx] = sort(sirValues); + berGroups = berGroups(sortIdx); + + berAvg = nan(1, numel(sirValues)); + berMin = nan(1, numel(sirValues)); + berMax = nan(1, numel(sirValues)); + + for sirIdx = 1:numel(sirValues) + rawBer = berGroups{sirIdx}; + rawBer = rawBer(isfinite(rawBer)); + if isempty(rawBer) + continue + end + + scatter(repmat(sirValues(sirIdx), size(rawBer)), rawBer, ... + 12, ... + "Marker", ".", ... + "MarkerEdgeColor", curveColor, ... + "MarkerFaceColor", curveColor, ... + "HandleVisibility", "off"); + + berAvg(sirIdx) = mean(rawBer, "omitnan"); + berMin(sirIdx) = min(rawBer); + berMax(sirIdx) = max(rawBer); + end + + validAvg = isfinite(sirValues) & isfinite(berAvg); + if ~any(validAvg) + continue + end + + if plotOptions.use_boundedlines && exist("boundedline", "file") + yLower = max(berAvg - berMin, 0); + yUpper = max(berMax - berAvg, 0); + yBounds = [yLower(:), yUpper(:)]; + + [hl, hp] = boundedline(sirValues(:), berAvg(:), yBounds, ... + "alpha", "transparency", plotOptions.boundedline_alpha, ... + "cmap", curveColor, ... + "nan", "fill", ... + "orientation", "vert"); + set(hl, "LineStyle", "none", "Marker", "none", "HandleVisibility", "off"); + set(hp, "HandleVisibility", "off", "LineStyle", "none"); + end + + scatter(sirValues(validAvg), berAvg(validAvg), ... + 54, ... + "Marker", "o", ... + "MarkerEdgeColor", curveColor, ... + "MarkerFaceColor", "none", ... + "LineWidth", 1.2, ... + "DisplayName", storageName); + + fitMask = validAvg & berAvg > 0; + if nnz(fitMask) >= 2 + fitOrder = min(plotOptions.polyfit_order_max, nnz(fitMask) - 1); + fitCoeff = polyfit(sirValues(fitMask), log10(berAvg(fitMask)), fitOrder); + xFit = linspace(min(sirValues(fitMask)), max(sirValues(fitMask)), 300); + yFit = 10 .^ polyval(fitCoeff, xFit); + + plot(xFit, yFit, ... + "LineWidth", 1.2, ... + "LineStyle", "--", ... + "Marker", "none", ... + "Color", curveColor, ... + "HandleVisibility", "off"); + end + end + + if hasBlockUpdate + title(sprintf("BER over SIR, PAM %.0f, path %.0f m, block update = %g", ... + pamLevel, plotOptions.path_length, blockUpdates(blockIdx))); + else + title(sprintf("BER over SIR, PAM %.0f, path %.0f m", ... + pamLevel, plotOptions.path_length)); + end + xlabel("SIR (dB)"); + ylabel("BER"); + xlim(plotOptions.xlim); + grid on + box on + + if exist("beautifyBERplot", "file") + beautifyBERplot("logscale", true, "setcolors", false, "setmarkers", false); + else + set(gca, "YScale", "log"); + end + legend("Location", "best", "Interpreter", "none"); + end +end +end + +function [sirValues, berGroups] = collectBerBySir(whConfig, storageName, filter) +sirValues = []; +berGroups = {}; + +for linIdx = 1:numel(whConfig.sto.(storageName)) + storedValue = whConfig.sto.(storageName){linIdx}; + if isempty(storedValue) + continue + end + + physStruct = linearIndexToStruct(whConfig, linIdx); + if ~matchesFilter(physStruct, filter) + continue + end + + berValues = extractBerValues(storedValue); + if isempty(berValues) + continue + end + + sirValue = physStruct.sir; + sirIdx = find(sirValues == sirValue, 1); + if isempty(sirIdx) + sirValues(end+1) = sirValue; %#ok + berGroups{end+1} = berValues; %#ok + else + berGroups{sirIdx} = [berGroups{sirIdx}, berValues]; %#ok + end +end +end + +function tf = matchesFilter(physStruct, filter) +tf = true; +filterFields = fieldnames(filter); + +for fieldIdx = 1:numel(filterFields) + fieldName = filterFields{fieldIdx}; + if ~isfield(physStruct, fieldName) || physStruct.(fieldName) ~= filter.(fieldName) + tf = false; + return + end +end +end + +function berValues = extractBerValues(value) +berValues = []; + +if isstruct(value) + if isfield(value, "metrics") && isfield(value.metrics, "BER") + berValues(end+1) = value.metrics.BER; + end +elseif iscell(value) + for valueIdx = 1:numel(value) + berValues = [berValues, extractBerValues(value{valueIdx})]; %#ok + end +end +end diff --git a/projects/Diss/MPI_revisit/investigate_mpi_algorithms.m b/projects/Diss/MPI_revisit/investigate_mpi_algorithms.m index 5eb1a0f..60bc296 100644 --- a/projects/Diss/MPI_revisit/investigate_mpi_algorithms.m +++ b/projects/Diss/MPI_revisit/investigate_mpi_algorithms.m @@ -3,8 +3,8 @@ dsp_options = struct(); dsp_options.mode = "run_id"; dsp_options.recipe = @mpi_recipe_dev; dsp_options.append_to_db = false; -dsp_options.start_occurence = 2; -dsp_options.max_occurences = 1; +dsp_options.start_occurence = 1; +dsp_options.max_occurences = 15; dsp_options.debug_plots = false; dsp_options.database_type = "mysql"; @@ -23,95 +23,58 @@ db = DBHandler("dataBase", [dsp_options.dataBase],... "user", dsp_options.user, "password", dsp_options.password); %% -fp = QueryFilter(); -fp.where('Runs','run_id','GREATER_EQUAL',3153); -fp.where('Runs', 'symbolrate', 'EQUALS', 112e9); -fp.where('Runs', 'fiber_length', 'EQUALS', 0); -fp.where('Runs', 'interference_path_length', 'EQUALS', 1000); -fp.where('Runs', 'sir', 'LESS_EQUAL', 20); -fp.where('Runs', 'db_mode', 'EQUALS', '"no_db"'); -% fp.where('Runs', 'is_mpi', 'EQUALS', 1); -fp.where('Runs', 'pam_level', 'EQUALS', 4); -% fp.where('Runs', 'v_bias', 'EQUALS', 2.65); -[dataTable,~] = db.queryDB(fp, db.getTableFieldNames('Runs')); +pamformats = [4,6,8]; +baudrates = [112e9,96e9,72e9]; -% [~, uniqueSirRows] = unique(dataTable.sir, "stable"); -% dataTable = dataTable(uniqueSirRows, :); +for i = 1:3 + + B = baudrates(i); + M = pamformats(i); + + fp = QueryFilter(); + fp.where('Runs','run_id','GREATER_EQUAL',3153); + fp.where('Runs', 'symbolrate', 'EQUALS', B); % 72 96 112 + fp.where('Runs', 'fiber_length', 'EQUALS', 0); + fp.where('Runs', 'interference_path_length', 'EQUALS', 1000); + fp.where('Runs', 'sir', 'LESS_EQUAL', 20); + fp.where('Runs', 'db_mode', 'EQUALS', '"no_db"'); + % fp.where('Runs', 'is_mpi', 'EQUALS', 1); + fp.where('Runs', 'pam_level', 'EQUALS', M); + % fp.where('Runs', 'v_bias', 'EQUALS', 2.65); + + [dataTable,~] = db.queryDB(fp, db.getTableFieldNames('Runs')); + + [~, sortIdx] = sort(dataTable.sir, 'descend'); + dataTable = dataTable(sortIdx, :); + run_ids = dataTable.run_id; + + %% === Warehouse setup === + dsp_options.userParameters = struct(); + dsp_options.userParameters.block_update = 1;%[1,2,4,8,16,32,64,112,224,448,512,1024,2048,2048*2,2048*4];%%logspace(-3.8,-1,22);%[linspace(2,4096,22)]; + wh = DataStorage(dsp_options.userParameters); + + + %% + n_realizations = (dsp_options.max_occurences - dsp_options.start_occurence + 1); + n_userparams = prod(wh.dim); + n_run_ids = numel(run_ids); + parallel_jobs = n_userparams * n_run_ids; + queried_jobs = n_realizations * n_userparams * n_run_ids; + + fprintf("-> [ %d run_id(s) × %d userParam combination(s) = %d parallel job(s) ] × %d realizations = %d total jobs \n", ... + n_run_ids, n_userparams, parallel_jobs, n_realizations, queried_jobs); + + %% === Run === + % submitJobs returns the raw per-job results and the filled Warehouse. For a + % single run_id and remove_dc = [0, 1], results is a 2-by-1 cell array. + [results, wh] = submitJobs(run_ids, dsp_options, processingMode.serial, ... + "wh", wh, ... + "waitbar", true); + + % savepath = fullfile("C:","Users","Silas","Documents","MATLAB","imdd_simulation","projects","Diss","MPI_revisit","algorithms",sprintf("pam_%d_results.mat",M)); + % save(wh,savepath) -[~, sortIdx] = sort(dataTable.sir, 'descend'); -dataTable = dataTable(sortIdx, :); -dataTable = dataTable(1, :); -run_ids = dataTable.run_id; -if isempty(run_ids) - error("run_minimal_recipe:MissingRunIds", ... - "Set run_ids to one or more known run IDs before running this example."); end -%% === Warehouse setup === -dsp_options.userParameters = struct(); -% dsp_options.userParameters.smoothing_length = [0,linspace(100,1000,10),2500,5000,10000]; -wh = DataStorage(dsp_options.userParameters); - - -%% -n_realizations = (dsp_options.max_occurences - dsp_options.start_occurence + 1); -n_userparams = prod(wh.dim); -n_run_ids = numel(run_ids); -parallel_jobs = n_userparams * n_run_ids; -queried_jobs = n_realizations * n_userparams * n_run_ids; - -fprintf("-> [ %d run_id(s) × %d userParam combination(s) = %d parallel job(s) ] × %d realizations = %d total jobs \n", ... - n_run_ids, n_userparams, parallel_jobs, n_realizations, queried_jobs); - -%% === Run === -% submitJobs returns the raw per-job results and the filled Warehouse. For a -% single run_id and remove_dc = [0, 1], results is a 2-by-1 cell array. -[results, wh] = submitJobs(run_ids, dsp_options, processingMode.serial, ... - "wh", wh, ... - "waitbar", true); - -%% Analyze - -storageNames = fieldnames(wh.sto); -result = wh.sto.(storageNames{1}); -result = result(:).'; - -% Each stored entry is a package cell containing one or more occurrences. -ber_all = cell(1,numel(result)); -for result_idx = 1:numel(result) - packageCell = result{result_idx}; - if isempty(packageCell) - ber_all{result_idx} = NaN; - continue - end - - ber_values = nan(1,numel(packageCell)); - for package_idx = 1:numel(packageCell) - if isstruct(packageCell{package_idx}) && isfield(packageCell{package_idx},"metrics") - ber_values(package_idx) = packageCell{package_idx}.metrics.BER; - end - end - ber_all{result_idx} = ber_values; -end - -maxPackages = max(cellfun(@numel,ber_all)); -ber_mat = nan(maxPackages,numel(ber_all)); -for k = 1:numel(ber_all) - ber_mat(1:numel(ber_all{k}),k) = ber_all{k}; -end -ber = mean(ber_mat,1,"omitnan"); - -x = dataTable.sir(:).'; -if numel(x) ~= numel(ber) - x = 1:numel(ber); -end - -figure(2026); clf; hold on -plot(x,ber); -for k = 1:numel(x) - scatter(repmat(x(k),maxPackages,1),ber_mat(:,k)) -end -beautifyBERplot("logscale",true,"setcolors",false,"setmarkers",true); - diff --git a/projects/Diss/MPI_revisit/parallelization_analysis/RUN_parallel_analysis.m b/projects/Diss/MPI_revisit/parallelization_analysis/RUN_parallel_analysis.m new file mode 100644 index 0000000..f786e4c --- /dev/null +++ b/projects/Diss/MPI_revisit/parallelization_analysis/RUN_parallel_analysis.m @@ -0,0 +1,73 @@ +% === DSP settings === +dsp_options = struct(); +dsp_options.mode = "run_id"; +dsp_options.recipe = @mpi_recipe_dev; +dsp_options.append_to_db = false; +dsp_options.start_occurence = 1; +dsp_options.max_occurences = 15; +dsp_options.debug_plots = false; + +dsp_options.database_type = "mysql"; + +dsp_options.dataBase = "labor"; +dsp_options.storage_path = "W:\labdata\ECOC Silas\ecoc_2025"; + +dsp_options.server = "192.168.178.192"; +dsp_options.port = 3306; +dsp_options.user = "silas"; +dsp_options.password = "silas"; + +db = DBHandler("dataBase", [dsp_options.dataBase],... + "type", dsp_options.database_type,... + "server", dsp_options.server,... + "user", dsp_options.user, "password", dsp_options.password); + +%% +fp = QueryFilter(); +fp.where('Runs','run_id','GREATER_EQUAL',3153); +fp.where('Runs', 'symbolrate', 'EQUALS', 72e9); % 72 96 112 +fp.where('Runs', 'fiber_length', 'EQUALS', 0); +% fp.where('Runs', 'interference_path_length', 'EQUALS', 1000); +% fp.where('Runs', 'sir', 'LESS_EQUAL', 50); +fp.where('Runs', 'db_mode', 'EQUALS', '"no_db"'); +% fp.where('Runs', 'is_mpi', 'EQUALS', 1); +fp.where('Runs', 'pam_level', 'EQUALS', 8); +% fp.where('Runs', 'v_bias', 'EQUALS', 2.65); + +[dataTable,~] = db.queryDB(fp, db.getTableFieldNames('Runs')); + +% [~, uniqueSirRows] = unique(dataTable.sir, "stable"); +% dataTable = dataTable(uniqueSirRows, :); + +[~, sortIdx] = sort(dataTable.sir, 'descend'); +dataTable = dataTable(sortIdx, :); +% dataTable = dataTable(1, :); +run_ids = dataTable.run_id; + +if isempty(run_ids) + error("run_minimal_recipe:MissingRunIds", ... + "Set run_ids to one or more known run IDs before running this example."); +end + +%% === Warehouse setup === +dsp_options.userParameters = struct(); +dsp_options.userParameters.block_update = [1,2,4,8,16,32,64,112,224,448,512,1024,2048,2048*2,2048*4];%%logspace(-3.8,-1,22);%[linspace(2,4096,22)]; +wh = DataStorage(dsp_options.userParameters); + + +%% +n_realizations = (dsp_options.max_occurences - dsp_options.start_occurence + 1); +n_userparams = prod(wh.dim); +n_run_ids = numel(run_ids); +parallel_jobs = n_userparams * n_run_ids; +queried_jobs = n_realizations * n_userparams * n_run_ids; + +fprintf("-> [ %d run_id(s) × %d userParam combination(s) = %d parallel job(s) ] × %d realizations = %d total jobs \n", ... + n_run_ids, n_userparams, parallel_jobs, n_realizations, queried_jobs); + +%% === Run === +% submitJobs returns the raw per-job results and the filled Warehouse. For a +% single run_id and remove_dc = [0, 1], results is a 2-by-1 cell array. +[results, wh] = submitJobs(run_ids, dsp_options, processingMode.parallel, ... + "wh", wh, ... + "waitbar", true); diff --git a/projects/Diss/MPI_revisit/parallelization_analysis/parallelization_pam4.fig b/projects/Diss/MPI_revisit/parallelization_analysis/parallelization_pam4.fig new file mode 100644 index 0000000..b2482da Binary files /dev/null and b/projects/Diss/MPI_revisit/parallelization_analysis/parallelization_pam4.fig differ diff --git a/projects/Diss/MPI_revisit/parallelization_analysis/plot_parallelization_vs_sir_pam4.m b/projects/Diss/MPI_revisit/parallelization_analysis/plot_parallelization_vs_sir_pam4.m new file mode 100644 index 0000000..368a72e --- /dev/null +++ b/projects/Diss/MPI_revisit/parallelization_analysis/plot_parallelization_vs_sir_pam4.m @@ -0,0 +1,360 @@ +%% Analyze +% === DSP settings === +dsp_options = struct(); +dsp_options.mode = "run_id"; +dsp_options.recipe = @mpi_recipe_dev; +dsp_options.append_to_db = false; +dsp_options.start_occurence = 1; +dsp_options.max_occurences = 15; +dsp_options.debug_plots = false; + +dsp_options.database_type = "mysql"; + +dsp_options.dataBase = "labor"; +dsp_options.storage_path = "W:\labdata\ECOC Silas\ecoc_2025"; + +dsp_options.server = "192.168.178.192"; +dsp_options.port = 3306; +dsp_options.user = "silas"; +dsp_options.password = "silas"; + +db = DBHandler("dataBase", [dsp_options.dataBase],... + "type", dsp_options.database_type,... + "server", dsp_options.server,... + "user", dsp_options.user, "password", dsp_options.password); + +%% + +fp = QueryFilter(); +fp.where('Runs','run_id','GREATER_EQUAL',3153); +fp.where('Runs', 'symbolrate', 'EQUALS', 72e9); % 72 96 112 +fp.where('Runs', 'fiber_length', 'EQUALS', 0); +fp.where('Runs', 'interference_path_length', 'EQUALS', 1000); +% fp.where('Runs', 'sir', 'LESS_EQUAL', 50); +fp.where('Runs', 'db_mode', 'EQUALS', '"no_db"'); +% fp.where('Runs', 'is_mpi', 'EQUALS', 1); +fp.where('Runs', 'pam_level', 'EQUALS', 8); +% fp.where('Runs', 'v_bias', 'EQUALS', 2.65); + +[dataTable,~] = db.queryDB(fp, db.getTableFieldNames('Runs')); + +% [~, uniqueSirRows] = unique(dataTable.sir, "stable"); +% dataTable = dataTable(uniqueSirRows, :); + +[~, sortIdx] = sort(dataTable.sir, 'descend'); +dataTable = dataTable(sortIdx, :); +% dataTable = dataTable(1, :); +run_ids = dataTable.run_id; +%% + +wh_ = load("C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\Diss\MPI_revisit\parallelization_analysis\wh_block_update_pam8_combined.mat"); +wh = wh_.wh; + +% wh_ = load("C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\Diss\MPI_revisit\parallelization_analysis\wh_block_update_pam4_short_blocks.mat"); +% wh2 = wh_.wh; +% +% wh = mergeDataStorages(wh1,wh2); + +%% + + +storageNames = fieldnames(wh.sto); +ber_by_storage = struct(); +ber_mat_by_storage = struct(); + +sir_values = dataTable.sir(:).'; +block_updates = wh.parameter.block_update.values; +storage_parameter_names = cellstr(wh.fn); +block_axis = find(strcmp(storage_parameter_names,"block_update"),1); + +if exist('linspecer','file') + curve_colors = linspecer(max(numel(storageNames),1)); +else + curve_colors = lines(max(numel(storageNames),1)); +end + +fec_ber_threshold = 2e-2;%3.8e-3; +fit_coeff_by_storage = struct(); +required_sir_fec_by_storage = struct(); +for storage_idx = 1:numel(storageNames) + storageName = storageNames{storage_idx}; + fit_coeff_by_storage.(storageName) = cell(1,numel(block_updates)); + required_sir_fec_by_storage.(storageName) = nan(1,numel(block_updates)); +end + +for block_idx = 1:numel(block_updates) + block_update = block_updates(block_idx); + + figure(3000 + block_idx); clf; hold on + for storage_idx = 1:numel(storageNames) + storageName = storageNames{storage_idx}; + curve_color = curve_colors(storage_idx,:); + result = wh.sto.(storageName); + + if ~isempty(block_axis) && size(result,block_axis) == numel(block_updates) + result_subs = repmat({':'},1,ndims(result)); + result_subs{block_axis} = block_idx; + result_for_block = result(result_subs{:}); + result_for_block = result_for_block(:).'; + else + result_for_block = result(:).'; + end + + % Each stored entry is a package cell containing one or more occurrences. + ber_all = cell(1,numel(result_for_block)); + for result_idx = 1:numel(result_for_block) + packageCell = result_for_block{result_idx}; + if isempty(packageCell) + ber_all{result_idx} = NaN; + continue + end + + ber_values = nan(1,numel(packageCell)); + for package_idx = 1:numel(packageCell) + if isstruct(packageCell{package_idx}) && isfield(packageCell{package_idx},"metrics") + ber_values(package_idx) = packageCell{package_idx}.metrics.BER; + end + end + ber_all{result_idx} = ber_values; + end + + maxPackages = max(cellfun(@numel,ber_all)); + ber_mat = nan(maxPackages,numel(ber_all)); + for k = 1:numel(ber_all) + ber_mat(1:numel(ber_all{k}),k) = ber_all{k}; + end + + % Replace outliers in ber_mat with NaN (operate column-wise). + for col = 1:size(ber_mat,2) + colData = ber_mat(:,col); + if all(isnan(colData)); continue; end + mask = isoutlier(colData); + ber_mat(mask,col) = NaN; + end + + ber = mean(ber_mat,1,"omitnan"); + ber_min = nan(1,numel(ber)); + ber_max = nan(1,numel(ber)); + for k = 1:numel(ber) + ber_here = ber_mat(:,k); + ber_here = ber_here(isfinite(ber_here)); + if isempty(ber_here) + continue + end + ber_min(k) = min(ber_here); + ber_max(k) = max(ber_here); + end + + x = sir_values; + if numel(x) ~= numel(ber) + x = 1:numel(ber); + end + + ber_by_storage.(storageName){block_idx} = ber; + ber_mat_by_storage.(storageName){block_idx} = ber_mat; + + y_lower = max(ber - ber_min,0); + y_upper = max(ber_max - ber,0); + y_bounds = [y_lower(:), y_upper(:)]; + + [hl, hp] = boundedline(x(:),ber(:),y_bounds, ... + 'alpha', 'transparency', 0.1, ... + 'cmap', curve_color, ... + 'nan', 'fill', ... + 'orientation', 'vert'); + set(hl, ... + 'LineWidth',1.4, ... + 'LineStyle','none', ... + 'Marker','o', ... + 'Markersize',5,... + 'Color',curve_color, ... + 'DisplayName',storageName); + set(hp, ... + 'HandleVisibility','off', ... + 'LineStyle','none'); + + for k = 1:numel(x) + scatter(repmat(x(k),maxPackages,1),ber_mat(:,k), ... + 16, ... + 'Marker','.', ... + 'MarkerEdgeColor',curve_color, ... + 'MarkerFaceColor',curve_color, ... + 'HandleVisibility','off'); + end + + fit_mask = isfinite(x) & isfinite(ber) & ber > 0; + if nnz(fit_mask) >= 2 + fit_order = min(3,nnz(fit_mask)-1); + fit_coeff = polyfit(x(fit_mask),log10(ber(fit_mask)),fit_order); + x_fit = linspace(min(x(fit_mask)),max(x(fit_mask)),300); + y_fit = 10.^polyval(fit_coeff,x_fit); + fit_coeff_by_storage.(storageName){block_idx} = fit_coeff; + + sir_req = NaN; + threshold_mask = isfinite(y_fit) & y_fit <= fec_ber_threshold; + if any(threshold_mask) + first_threshold_idx = find(threshold_mask,1,"first"); + if first_threshold_idx == 1 + sir_req = x_fit(first_threshold_idx); + else + x_pair = x_fit(first_threshold_idx-1:first_threshold_idx); + y_pair = log10(y_fit(first_threshold_idx-1:first_threshold_idx)); + if all(isfinite(y_pair)) && diff(y_pair) ~= 0 + sir_req = interp1(y_pair,x_pair,log10(fec_ber_threshold), ... + "linear","extrap"); + else + sir_req = x_fit(first_threshold_idx); + end + end + end + required_sir_fec_by_storage.(storageName)(block_idx) = sir_req; + + plot(x_fit,y_fit, ... + 'LineWidth',1.2, ... + 'LineStyle','--', ... + 'Marker','none', ... + 'Color',curve_color, ... + 'HandleVisibility','off'); + end + end + + title(sprintf("BER over SIR, block update = %g",block_update)); + xlabel("SIR (dB)"); + ylabel("BER"); + xlim([15 35]); + + beautifyBERplot("logscale",true,"setcolors",false,"setmarkers",false); + legend("Location","best","Interpreter","none"); +end + +%% Required SIR to reach FEC over parallelization + +figure(5000); clf; hold on +for storage_idx = 1:numel(storageNames) + storageName = storageNames{storage_idx}; + curve_color = curve_colors(storage_idx,:); + required_sir = required_sir_fec_by_storage.(storageName); + valid = isfinite(block_updates) & isfinite(required_sir); + + if ~any(valid) + continue + end + + plot(block_updates(valid),required_sir(valid), ... + 'LineWidth',1.4, ... + 'LineStyle','-', ... + 'Marker','o', ... + 'MarkerSize',5, ... + 'Color',curve_color, ... + 'DisplayName',storageName); +end + +set(gca,'XScale','log'); +xticks(block_updates); +xticklabels(string(block_updates)); +grid on +box on +xlabel("Parallelization / block update"); +ylabel(sprintf("Required SIR at BER = %.1e (dB)",fec_ber_threshold)); +title("Required SIR to reach FEC over parallelization"); +legend("Location","best","Interpreter","none"); + + +%% + + +function whMerged = mergeDataStorages(varargin) +% mergeDataStorages merges compatible DataStorage objects by physical indices. +% Existing duplicate entries are concatenated when both values are package cells. + +whList = varargin; +templateWh = whList{1}; +paramNames = cellstr(templateWh.fn); + +mergedParams = struct(); +for param_idx = 1:numel(paramNames) + paramName = paramNames{param_idx}; + values = templateWh.parameter.(paramName).values(:).'; + + for wh_idx = 2:numel(whList) + curWh = whList{wh_idx}; + if ~isequal(sort(cellstr(curWh.fn)),sort(paramNames)) + error("mergeDataStorages:ParameterMismatch", ... + "All warehouses must use the same parameter names."); + end + if ~isfield(curWh.parameter,paramName) + error("mergeDataStorages:ParameterMismatch", ... + "Warehouse %d does not contain parameter '%s'.",wh_idx,paramName); + end + + newValues = curWh.parameter.(paramName).values(:).'; + for value_idx = 1:numel(newValues) + if ~ismember(newValues(value_idx),values) + values(end+1) = newValues(value_idx); %#ok + end + end + end + + if strcmp(paramName,"block_update") && isnumeric(values) + values = sort(values); + end + + mergedParams.(paramName) = values; +end + +whMerged = DataStorage(mergedParams); + +for wh_idx = 1:numel(whList) + curWh = whList{wh_idx}; + curStorageNames = fieldnames(curWh.sto); + + for storage_idx = 1:numel(curStorageNames) + storageName = curStorageNames{storage_idx}; + if ~isfield(whMerged.sto,storageName) + whMerged.addStorage(storageName); + end + + for lin_idx = 1:numel(curWh.sto.(storageName)) + storedValue = curWh.sto.(storageName){lin_idx}; + if isempty(storedValue) + continue + end + + [physValues,physNames] = curWh.getPhysIndicesByLinIndex(lin_idx); + physNames = cellfun(@char,physNames,'UniformOutput',false); + targetSubscripts = cell(1,numel(whMerged.fn)); + + for param_idx = 1:numel(whMerged.fn) + paramName = char(whMerged.fn(param_idx)); + sourceParamIdx = find(strcmp(physNames,paramName),1); + if isempty(sourceParamIdx) + error("mergeDataStorages:ParameterMismatch", ... + "Storage entry is missing parameter '%s'.",paramName); + end + targetSubscripts{param_idx} = whMerged.getIndexByPhys(paramName,physValues{sourceParamIdx}); + end + + if isscalar(targetSubscripts) + targetLinIdx = targetSubscripts{1}; + else + targetLinIdx = sub2ind(whMerged.getStorageSize(),targetSubscripts{:}); + end + + existingValue = whMerged.sto.(storageName){targetLinIdx}; + whMerged.sto.(storageName){targetLinIdx} = mergeStoredValue(existingValue,storedValue); + end + end +end +end + +function out = mergeStoredValue(existingValue,newValue) +if isempty(existingValue) + out = newValue; +elseif iscell(existingValue) && iscell(newValue) + out = [existingValue(:); newValue(:)].'; +else + out = newValue; +end +end + + diff --git a/projects/Diss/MPI_revisit/plot_mpi_timesignal.m b/projects/Diss/MPI_revisit/plot_mpi_timesignal.m index ec3462c..dcd4087 100644 --- a/projects/Diss/MPI_revisit/plot_mpi_timesignal.m +++ b/projects/Diss/MPI_revisit/plot_mpi_timesignal.m @@ -127,7 +127,9 @@ for sirval = [20,30,50] pngTikzPath = "04_Experimental_Evaluation/tikz/mpi/" + exportName + "-1.png"; figure(sirval + pathlen); - exportNakedPlotWithTikz(gca, pngFile, tikzFile, pngTikzPath, 150); + exportLevelScatterTikz(gca, pngFile, tikzFile, pngTikzPath, ... + "resolutionDpi", 150, ... + "generatedBy", "plot_mpi_timesignal.m"); end end @@ -145,309 +147,3 @@ end % var_slope = p(1); % var_offset = p(2);c -function exportNakedPlotWithTikz(sourceAx, pngFile, tikzFile, pngTikzPath, resolutionDpi) -outputDir = fileparts(pngFile); -if ~exist(outputDir, "dir") - mkdir(outputDir); -end - -exportScatterLayerPng(sourceAx, pngFile, resolutionDpi); - -linePlots = collectTikzLinePlots(sourceAx); -textAnnotations = collectTikzTextAnnotations(sourceAx); -writeTikzPngAxis(tikzFile, pngTikzPath, sourceAx, linePlots, textAnnotations); -end - -function exportScatterLayerPng(sourceAx, pngFile, resolutionDpi) -sourceFig = ancestor(sourceAx, "figure"); -xLimits = sourceAx.XLim; -yLimits = sourceAx.YLim; -scatterHandles = findGraphicsObjectsByType(sourceAx, "scatter"); - -exportFig = figure( ... - "Visible", "off", ... - "Color", "white", ... - "Units", sourceFig.Units, ... - "Position", sourceFig.Position); -cleanupFigure = onCleanup(@() close(exportFig)); - -exportAx = copyobj(sourceAx, exportFig); -exportAx.Units = "normalized"; -exportAx.Position = [0 0 1 1]; -exportAx.XLim = xLimits; -exportAx.YLim = yLimits; -exportAx.XLimMode = "manual"; -exportAx.YLimMode = "manual"; -stripAxesToScatterOnly(exportAx); -hideAxesInkForRasterExport(exportAx); - -fprintf("Exporting scatter PNG with XLim=[%g %g], YLim=[%g %g], scatter objects=%d: %s\n", ... - xLimits(1), xLimits(2), yLimits(1), yLimits(2), numel(scatterHandles), pngFile); -drawnow; -exportFig.PaperPositionMode = "auto"; -print(exportFig, char(pngFile), "-dpng", sprintf("-r%d", resolutionDpi)); -end - -function stripAxesToScatterOnly(ax) -plotObjects = findall(ax); -for objIdx = 1:numel(plotObjects) - obj = plotObjects(objIdx); - if obj == ax || ~isvalid(obj) - continue - end - - objType = string(get(obj, "Type")); - if lower(objType) ~= "scatter" - delete(obj); - end -end -end - -function hideAxesInkForRasterExport(ax) -ax.Visible = "on"; -ax.Color = "white"; -ax.Box = "off"; -ax.XGrid = "off"; -ax.YGrid = "off"; -ax.XMinorGrid = "off"; -ax.YMinorGrid = "off"; -ax.XTick = []; -ax.YTick = []; -ax.XColor = "white"; -ax.YColor = "white"; -ax.Title.String = ""; -ax.XLabel.String = ""; -ax.YLabel.String = ""; -end - -function handles = findGraphicsObjectsByType(parentHandle, objectType) -allHandles = findall(parentHandle); -matches = false(size(allHandles)); -for handleIdx = 1:numel(allHandles) - currentHandle = allHandles(handleIdx); - if ~isvalid(currentHandle) - continue - end - - matches(handleIdx) = lower(string(get(currentHandle, "Type"))) == lower(string(objectType)); -end - -handles = allHandles(matches); -end - -function linePlots = collectTikzLinePlots(sourceAx) -lineHandles = findGraphicsObjectsByType(sourceAx, "line"); -lineHandles = flipud(lineHandles(:)); - -linePlots = repmat(struct( ... - "xData", [], ... - "yData", [], ... - "color", [], ... - "lineWidth", []), numel(lineHandles), 1); -validLineCount = 0; - -for lineIdx = 1:numel(lineHandles) - lineHandle = lineHandles(lineIdx); - xData = lineHandle.XData(:); - yData = lineHandle.YData(:); - validSamples = isfinite(xData) & isfinite(yData); - - if ~any(validSamples) - continue - end - - validLineCount = validLineCount + 1; - linePlots(validLineCount).xData = xData(validSamples); - linePlots(validLineCount).yData = yData(validSamples); - linePlots(validLineCount).color = lineHandle.Color; - linePlots(validLineCount).lineWidth = lineHandle.LineWidth; -end - -linePlots = linePlots(1:validLineCount); -end - -function textAnnotations = collectTikzTextAnnotations(sourceAx) -textHandles = findGraphicsObjectsByType(sourceAx, "text"); -textHandles = flipud(textHandles(:)); -xLimits = sourceAx.XLim; -yLimits = sourceAx.YLim; - -textAnnotations = repmat(struct( ... - "x", [], ... - "y", [], ... - "text", "", ... - "anchor", "center"), numel(textHandles), 1); -validTextCount = 0; - -for textIdx = 1:numel(textHandles) - textHandle = textHandles(textIdx); - labelText = string(textHandle.String); - if strlength(strtrim(labelText)) == 0 - continue - end - - textPosition = textHandle.Position; - if textPosition(1) < xLimits(1) || textPosition(1) > xLimits(2) || ... - textPosition(2) < yLimits(1) || textPosition(2) > yLimits(2) - continue - end - - validTextCount = validTextCount + 1; - textAnnotations(validTextCount).x = textPosition(1); - textAnnotations(validTextCount).y = textPosition(2); - textAnnotations(validTextCount).text = labelText; - textAnnotations(validTextCount).anchor = matlabTextAlignmentToTikzAnchor( ... - textHandle.HorizontalAlignment, textHandle.VerticalAlignment); -end - -textAnnotations = textAnnotations(1:validTextCount); -end - -function writeTikzPngAxis(tikzFile, pngTikzPath, sourceAx, linePlots, textAnnotations) -outputDir = fileparts(tikzFile); -if ~exist(outputDir, "dir") - mkdir(outputDir); -end - -xLimits = sourceAx.XLim; -yLimits = sourceAx.YLim; -fid = fopen(tikzFile, "w"); -if fid < 0 - error("plot_mpi_timesignal:TikzOpenFailed", ... - "Could not open TikZ export file: %s", tikzFile); -end -cleanupFile = onCleanup(@() fclose(fid)); - -axisOptions = { - " every axis/.append style={font=\scriptsize}," - "width=\fwidth," - "height=\fheight," - "at={(0\fwidth,0\fheight)}," - "scale only axis," - "axis on top," - "xmin=" + formatPgfNumber(xLimits(1)) + "," - "xmax=" + formatPgfNumber(xLimits(2)) + "," - "xlabel style={font=\color{white!15!black}\small}," - "xlabel={Time in $\mu$s}," - "ymin=" + formatPgfNumber(yLimits(1)) + "," - "ymax=" + formatPgfNumber(yLimits(2)) + "," - "ylabel style={font=\color{white!15!black}\small}," - "ylabel={Normalized Amplitude}," - "axis background/.style={fill=white}," - "xmajorgrids," - "ymajorgrids," - "grid style={line width=0.4pt, dotted, color=black!20}," - "scaled ticks=false," - "tick label style={/pgf/number format/fixed, /pgf/number format/1000 sep={}}," - "legend columns=1" - }; - -fprintf(fid, "%% This file was generated by plot_mpi_timesignal.m.\n"); -fprintf(fid, "%%\n"); -for lineIdx = 1:numel(linePlots) - fprintf(fid, "%s\n", formatTikzColorDefinition(lineIdx, linePlots(lineIdx).color)); -end -if ~isempty(linePlots) - fprintf(fid, "%%\n"); -end -fprintf(fid, "%s\n\n", "\begin{tikzpicture}"); -fprintf(fid, "%s\n", "\begin{axis}[%"); -for optionIdx = 1:numel(axisOptions) - fprintf(fid, "%s\n", axisOptions{optionIdx}); -end -fprintf(fid, "%s\n", "]"); - -graphicsLine = sprintf( ... - "\\addplot [forget plot] graphics [xmin=%s, xmax=%s, ymin=%s, ymax=%s] {%s};", ... - formatPgfNumber(xLimits(1)), ... - formatPgfNumber(xLimits(2)), ... - formatPgfNumber(yLimits(1)), ... - formatPgfNumber(yLimits(2)), ... - char(strrep(pngTikzPath, "\", "/"))); -fprintf(fid, "%s\n\n", graphicsLine); -writeTikzLinePlots(fid, linePlots); -writeTikzTextAnnotations(fid, textAnnotations); -fprintf(fid, "%s\n\n", "\end{axis}"); -fprintf(fid, "%s", "\end{tikzpicture}%"); -end - -function valueText = formatPgfNumber(value) -valueText = string(sprintf("%.15g", value)); -end - -function colorDefinition = formatTikzColorDefinition(colorIdx, rgbColor) -colorDefinition = sprintf( ... - "\\definecolor{mycolor%d}{rgb}{%.5f,%.5f,%.5f}%%", ... - colorIdx, rgbColor(1), rgbColor(2), rgbColor(3)); -end - -function writeTikzLinePlots(fid, linePlots) -for lineIdx = 1:numel(linePlots) - fprintf(fid, "\\addplot [color=mycolor%d, line width=%.1fpt, forget plot]\n", ... - lineIdx, linePlots(lineIdx).lineWidth); - fprintf(fid, " table[row sep=crcr]{%%\n"); - - for sampleIdx = 1:numel(linePlots(lineIdx).xData) - fprintf(fid, "%s\t%s\\\\\n", ... - formatPgfNumber(linePlots(lineIdx).xData(sampleIdx)), ... - formatPgfNumber(linePlots(lineIdx).yData(sampleIdx))); - end - - fprintf(fid, "};\n\n"); -end -end - -function writeTikzTextAnnotations(fid, textAnnotations) -for textIdx = 1:numel(textAnnotations) - fprintf(fid, "\\node[font=\\scriptsize, anchor=%s, fill=white, draw=black!40, rounded corners=2pt, inner sep=2pt]\n", ... - textAnnotations(textIdx).anchor); - fprintf(fid, " at (axis cs:%s,%s) {%s};\n\n", ... - formatPgfNumber(textAnnotations(textIdx).x), ... - formatPgfNumber(textAnnotations(textIdx).y), ... - escapeTikzText(textAnnotations(textIdx).text)); -end -end - -function anchor = matlabTextAlignmentToTikzAnchor(horizontalAlignment, verticalAlignment) -horizontalAlignment = string(horizontalAlignment); -verticalAlignment = string(verticalAlignment); - -if horizontalAlignment == "left" - horizontalAnchor = "west"; -elseif horizontalAlignment == "right" - horizontalAnchor = "east"; -else - horizontalAnchor = ""; -end - -if verticalAlignment == "top" - verticalAnchor = "north"; -elseif verticalAlignment == "bottom" - verticalAnchor = "south"; -else - verticalAnchor = ""; -end - -if strlength(horizontalAnchor) > 0 && strlength(verticalAnchor) > 0 - anchor = verticalAnchor + " " + horizontalAnchor; -elseif strlength(horizontalAnchor) > 0 - anchor = horizontalAnchor; -elseif strlength(verticalAnchor) > 0 - anchor = verticalAnchor; -else - anchor = "center"; -end -end - -function textOut = escapeTikzText(textIn) -textOut = char(textIn); -% textOut = strrep(textOut, "\", "\textbackslash{}"); -textOut = strrep(textOut, "{", "\{"); -textOut = strrep(textOut, "}", "\}"); -textOut = strrep(textOut, "%", "\%"); -textOut = strrep(textOut, "&", "\&"); -textOut = strrep(textOut, "#", "\#"); -textOut = strrep(textOut, "_", "\_"); -% textOut = strrep(textOut, "$", "\$"); -end - diff --git a/projects/ECOC_2025_MPI/plots_from_database/plot_mpi_trial.m b/projects/ECOC_2025_MPI/plots_from_database/plot_mpi_trial.m index 0359b37..5c05de5 100644 --- a/projects/ECOC_2025_MPI/plots_from_database/plot_mpi_trial.m +++ b/projects/ECOC_2025_MPI/plots_from_database/plot_mpi_trial.m @@ -1,49 +1,10 @@ -% dsp_options.database_type = 'mysql'; -% dsp_options.dataBase = 'labor'; -% dsp_options.storage_path = 'Z:\2025\ECOC Silas\ecoc_2025\'; -% database = DBHandler("dataBase", [dsp_options.dataBase], "type", dsp_options.database_type); -% filterParams = database.tables; -% filterParams.Runs.loop_id = 209; -% % filterParams.Configurations = struct( ... -% % 'symbolrate', 112e9, ... %[224,336,360,390,420,448] -% % 'fiber_length', 0, ... -% % 'db_mode', '"no_db"', ... -% % 'interference_attenuation', [], ... -% % 'interference_path_length', 1000, ... -% % 'is_mpi', 1, ... -% % 'pam_level', 4, ... -% % 'wavelength', 1310, ... -% % 'precomp_amp', [], ... -% % 'signal_attenuation', [], ... -% % 'v_awg', [], ... -% % 'v_bias', [] ... -% % ); -% -% % if 1 -% % % filterParams.EqualizerParameters.dc_buffer_len = 1; -% % filterParams.EqualizerParameters.ffe_buffer_len = 1; -% % filterParams.EqualizerParameters.smoothing_buffer_len = 4096; -% % filterParams.EqualizerParameters.smoothing_buffer_update = 224; -% % filterParams.EqualizerParameters.DCmu = 0; -% % end -% a = database.getTableFieldNames('Runs'); -% b = database.getTableFieldNames('Results'); -% c = database.getTableFieldNames('EqualizerParameters'); -% d = [a;b;c]; -% -% [dataTable,~] = database.queryDB(filterParams, d); -% -% selectedFields = {'Configurations.run_id' 'Runs.loop_id' 'Runs.date_of_run' 'Runs.rx_raw_path' 'Runs.bitrate' 'Runs.v_bias' 'Runs.v_awg' 'Runs.precomp_amp' 'Runs.symbolrate' 'Runs.pam_level'... -% 'Runs.db_mode' 'Runs.rop_attenuation' 'Runs.is_mpi' 'Runs.interference_attenuation' 'Runs.interference_path_length' 'Runs.signal_attenuation' ... -% 'EqualizerParameters.equalizer_structure' 'EqualizerParameters.diff_precode' 'EqualizerParameters.eq_id' 'EqualizerParameters.dc_buffer_len' 'EqualizerParameters.ffe_buffer_len' 'EqualizerParameters.smoothing_buffer_len' 'EqualizerParameters.smoothing_buffer_update' 'EqualizerParameters.DCmu' 'Measurements.power_pd_in' ... -% 'Measurements.power_mpi_interference' 'Measurements.power_mpi_signal' 'Results.BER' 'Results.BER_precoded' 'Results.EVM' 'Results.SNR' 'Results.GMI' 'Results.Alpha' 'Results.date_of_processing'}; - db = DBHandler("type","mysql","dataBase",'labor'); fp = QueryFilter(); % fp.where('mpi_superview', 'loop_id','EQUALS', 209); +% fp.where('mpi_superview','run_id','GREATER_EQUAL',3153); fp.where('mpi_superview', 'symbolrate','EQUALS', 112e9); fp.where('mpi_superview', 'pam_level','EQUALS', 4); fn = [db.getTableFieldNames('mpi_superview')]; @@ -67,7 +28,7 @@ figure() tiledlayout(1, 1, 'TileSpacing', 'compact', 'Padding', 'compact'); y_here = 0; figcnt = 0; -for int_len = 1000%[0,50,300,1000] +for int_len = 1%[0,50,300,1000] figcnt = figcnt+1; % figure(int_len+1); nexttile; @@ -139,7 +100,7 @@ for int_len = 1000%[0,50,300,1000] % Modify values in 'interference_path_length' where the condition is met dataTable.interference_path_length(dataTable.interference_path_length < 101 & dataTable.interference_path_length > 1) = 50; - dataTable.interference_path_length(dataTable.interference_path_length == 1) = 0; + % dataTable.interference_path_length(dataTable.interference_path_length == 1) = 0; dataTable = dataTable(dataTable.interference_path_length == int_len, :); diff --git a/workerError.mat b/workerError.mat index e62af78..afabc04 100644 Binary files a/workerError.mat and b/workerError.mat differ