From f421348e5b31c9b4e7f9534841199fdfdcc79efb Mon Sep 17 00:00:00 2001 From: Silas Oettinghaus Date: Mon, 13 Jul 2026 19:59:57 +0200 Subject: [PATCH] before new database --- Classes/04_DSP/Equalizer/FFE.m | 73 ++- Classes/04_DSP/Equalizer/FFE_plain.m | 492 ++++++++++++++ .../04_DSP/Equalizer/MPI_reduction/FFE_A1.m | 311 +++++++++ .../Equalizer/MPI_reduction/FFE_A2Residual.m | 347 ++++++++++ .../MPI_reduction/FFE_A2TrackedLevels.m | 361 +++++++++++ .../Equalizer/MPI_reduction/FFE_DCTracking.m | 521 +++++++++++++++ Functions/EQ_recipes/mpi_recipe_dev.m | 601 +++++------------- Functions/EQ_visuals/showLevelScatter.m | 2 +- Tests/04_DSP/Equalizer/FFE_plain_test.m | 105 +++ Tests/04_DSP/Equalizer/FFE_test.m | 171 +++++ .../MPI_revisit/algorithms/PLOT_ber_vs_sir.m | 14 +- .../MPI_revisit/investigate_mpi_algorithms.m | 66 +- .../RUN_parallel_analysis.m | 2 +- .../plot_parallelization_vs_sir_pam4.m | 103 +-- workerError.mat | Bin 1328 -> 1173 bytes 15 files changed, 2682 insertions(+), 487 deletions(-) create mode 100644 Classes/04_DSP/Equalizer/FFE_plain.m create mode 100644 Classes/04_DSP/Equalizer/MPI_reduction/FFE_A1.m create mode 100644 Classes/04_DSP/Equalizer/MPI_reduction/FFE_A2Residual.m create mode 100644 Classes/04_DSP/Equalizer/MPI_reduction/FFE_A2TrackedLevels.m create mode 100644 Classes/04_DSP/Equalizer/MPI_reduction/FFE_DCTracking.m create mode 100644 Tests/04_DSP/Equalizer/FFE_plain_test.m diff --git a/Classes/04_DSP/Equalizer/FFE.m b/Classes/04_DSP/Equalizer/FFE.m index 0ae09ff..7c30919 100644 --- a/Classes/04_DSP/Equalizer/FFE.m +++ b/Classes/04_DSP/Equalizer/FFE.m @@ -42,6 +42,7 @@ classdef FFE < handle dc_level_avg_bufferlength_a2 dc_level_update_blocklength_a2 dc_smoothing_a2 + dc_level_decision_mode_a2 dc_level_weights_a2 % Adaptive DC-tracking loop @@ -119,6 +120,7 @@ classdef FFE < handle options.dc_level_avg_bufferlength_a2 = 0; options.dc_level_update_blocklength_a2 = 0; options.dc_smoothing_a2 = 0; + options.dc_level_decision_mode_a2 = "residual"; options.dc_level_weights_a2 = 0; options.dc_tracking_mu = 0; @@ -170,6 +172,11 @@ classdef FFE < handle 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_level_decision_mode_a2 = string(obj.dc_level_decision_mode_a2); + if ~any(obj.dc_level_decision_mode_a2 == ["residual", "tracked_levels"]) + builtin("error","FFE:InvalidA2DecisionMode", ... + "dc_level_decision_mode_a2 must be 'residual' or 'tracked_levels'."); + end obj.dc_tracking_buffer_len = floor(obj.dc_tracking_buffer_len); obj.ffe_update_buffer_len = floor(obj.ffe_update_buffer_len); obj.a2_level_weight_max = max(obj.a2_level_weight_max,0); @@ -357,6 +364,8 @@ classdef FFE < handle "A2 level-dependent MPI suppression requires obj.constellation to be set."); end + dc_level_decision_mode_a2 = obj.dc_level_decision_mode_a2; + dc_level_track_decision_a2 = dc_level_decision_mode_a2 == "tracked_levels"; 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); @@ -396,6 +405,7 @@ classdef FFE < handle obj.debug_struct.dc_level_weight = NaN(1,n_symbols_debug); obj.debug_struct.dc_level_valid_count = NaN(1,n_symbols_debug); obj.debug_struct.dc_level_symbol_idx = NaN(1,n_symbols_debug); + obj.debug_struct.dc_level_decision_level = NaN(1,n_symbols_debug); obj.debug_struct.dc_level_y_raw = NaN(1,n_symbols_debug); if training @@ -419,6 +429,7 @@ classdef FFE < handle dc_level_weight = 0; dc_level_valid_count = 0; dc_level_symbol_idx = NaN; + dc_level_decision_level = NaN; U = x(obj.order+sample-1:-1:sample); if dc_avg_enabled @@ -444,33 +455,60 @@ classdef FFE < handle symbol_idx = NaN; else [~,symbol_idx] = min(abs(d_hat(symbol) - constellation)); + dc_level_decision_level = constellation(symbol_idx); end else if ~always_ideal_decision - [~,symbol_idx] = min(abs(y(symbol) - constellation)); % decision for closest constellation point + if dc_level_enabled && dc_level_track_decision_a2 + % Tracked-level A2 moves the decision constellation; y itself stays unchanged. + dc_level_ramp_by_level = min(dc_level_valid_count_by_level / dc_level_buffer_len_a2,1); + dc_level_decision_constellation = constellation + ... + dc_level_weight_by_level .* dc_level_ramp_by_level .* dc_level_mpi_est_by_level; + [~,symbol_idx] = min(abs(y(symbol) - dc_level_decision_constellation)); + dc_level_decision_level = dc_level_decision_constellation(symbol_idx); + else + [~,symbol_idx] = min(abs(y(symbol) - constellation)); % decision for closest constellation point + dc_level_decision_level = constellation(symbol_idx); + end d_hat(symbol,1) = constellation(symbol_idx); else d_hat(symbol,1) = d(symbol); + [~,symbol_idx] = min(abs(d_hat(symbol) - constellation)); + dc_level_decision_level = constellation(symbol_idx); end end dc_level_symbol_idx = symbol_idx; if dc_level_enabled - mpi_err = y_raw - d_hat(symbol); 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 / 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); + if dc_level_track_decision_a2 + dc_level_buffer_value = y_raw; + y(symbol,1) = y_raw; + if isnan(dc_level_decision_level) + dc_level_decision_level = constellation(symbol_idx) + ... + dc_level_weight * dc_level_mpi_est; + end else - if ~always_ideal_decision - [~,symbol_idx] = min(abs(y(symbol) - constellation)); - d_hat(symbol,1) = constellation(symbol_idx); - else + mpi_err = y_raw - d_hat(symbol); + dc_level_buffer_value = mpi_err; + 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) - constellation)); + d_hat(symbol,1) = constellation(symbol_idx); + dc_level_decision_level = constellation(symbol_idx); + else + d_hat(symbol,1) = d(symbol); + [~,symbol_idx] = min(abs(d_hat(symbol) - constellation)); + dc_level_decision_level = constellation(symbol_idx); + end end end @@ -488,10 +526,10 @@ classdef FFE < handle 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; + if isfinite(dc_level_buffer_value) + dc_level_err_buffer(dc_level_symbol_idx,dc_level_err_buffer_pos) = dc_level_buffer_value; dc_level_err_sum_by_level(dc_level_symbol_idx) = ... - dc_level_err_sum_by_level(dc_level_symbol_idx) + mpi_err; + dc_level_err_sum_by_level(dc_level_symbol_idx) + dc_level_buffer_value; dc_level_buffer_valid_count_by_level(dc_level_symbol_idx) = ... dc_level_buffer_valid_count_by_level(dc_level_symbol_idx) + 1; else @@ -508,6 +546,11 @@ classdef FFE < handle 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); + if dc_level_track_decision_a2 + dc_level_mpi_est_by_level(dc_level_symbol_idx) = ... + dc_level_mpi_est_by_level(dc_level_symbol_idx) - ... + constellation(dc_level_symbol_idx); + end end else dc_level_valid_count_by_level = dc_level_buffer_valid_count_by_level; @@ -516,6 +559,11 @@ classdef FFE < handle 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); + if dc_level_track_decision_a2 + dc_level_mpi_est_by_level(dc_level_has_valid) = ... + dc_level_mpi_est_by_level(dc_level_has_valid) - ... + constellation(dc_level_has_valid); + end end end end @@ -675,6 +723,7 @@ classdef FFE < handle obj.debug_struct.dc_level_weight(1,symbol) = dc_level_weight; obj.debug_struct.dc_level_valid_count(1,symbol) = dc_level_valid_count; obj.debug_struct.dc_level_symbol_idx(1,symbol) = dc_level_symbol_idx; + obj.debug_struct.dc_level_decision_level(1,symbol) = dc_level_decision_level; obj.debug_struct.dc_level_y_raw(1,symbol) = y_raw; if training diff --git a/Classes/04_DSP/Equalizer/FFE_plain.m b/Classes/04_DSP/Equalizer/FFE_plain.m new file mode 100644 index 0000000..b9d14dd --- /dev/null +++ b/Classes/04_DSP/Equalizer/FFE_plain.m @@ -0,0 +1,492 @@ +classdef FFE_plain < handle + % Plain feed-forward equalizer with training, DD mode, and BER mu tuning. + % + % This class intentionally contains only the core FFE tap adaptation logic. + % MPI-specific DC tracking, A1/A2 suppression, and delayed update buffers + % stay in FFE and the MPI_reduction classes. + + properties + sps + order + e + e_tr + error + + len_tr + mu_tr + epochs_tr + + adaption_technique + dd_mode + mu_dd + epochs_dd + dd_len_fraction + + P + constellation + + decide + + save_debug = 0 + debug_struct + + optmize_mus = 0 + mu_optimization + mu_optimization_iter = 0 + mu_optimization_len + plot_mu_optimization = 0 + mu_optimization_fignum = 3010 + end + + methods + function obj = FFE_plain(options) + arguments(Input) + options.sps = 2 + options.order = 15 + + options.len_tr = 4096 + options.mu_tr = 0 + options.epochs_tr = 5 + + options.adaption_technique adaption_method = adaption_method.lms + options.dd_mode = 1 + options.mu_dd = 1e-5 + options.epochs_dd = 5 + options.dd_len_fraction = 1 + + options.decide = false + + options.save_debug = 0 + options.optmize_mus = 0 + options.mu_optimization_len = 2^15 + options.plot_mu_optimization = 0 + options.mu_optimization_fignum = 3010 + end + + fn = fieldnames(options); + for n = 1:numel(fn) + obj.(fn{n}) = options.(fn{n}); + end + + obj.sps = floor(obj.sps); + obj.order = floor(obj.order); + obj.dd_len_fraction = min(max(obj.dd_len_fraction,0),1); + obj.resetState(); + end + + function [X,Noi] = process(obj, X, D) + X = X.normalize("mode","rms"); + obj.constellation = unique(D.signal); + obj.resetState(); + + if obj.optmize_mus + obj.optimizeMus(X.signal,D.signal); + obj.resetState(); + end + + training = true; + showviz = false; + obj.equalize(X.signal,D.signal,obj.mu_tr,obj.epochs_tr,obj.len_tr,training,showviz); + obj.e_tr = obj.e; + + training = false; + n = X.length; + if obj.dd_mode + [signal,decision] = obj.equalize(X.signal,D.signal,obj.mu_dd,obj.epochs_dd,n,training,showviz); + else + [signal,decision] = obj.equalize(X.signal,D.signal,0,1,n,training,showviz); + end + + if obj.decide + X.signal = decision; + else + X.signal = signal; + end + + X.fs = D.fs; + X = X.logbookentry([num2str(obj.order),' tap FFE']); + + Noi = X - D; + end + + function [y,d_hat] = equalize(obj,x,d,mu,epochs,N,training,showviz) + arguments + obj + x + d + mu + epochs + N + training + showviz + end + + unused_showviz = showviz; %#ok + x = x(:); + d = d(:); + N = obj.validSampleLength(N,x,d); + n_symbols = N / obj.sps; + y = zeros(n_symbols,1); + d_hat = zeros(n_symbols,1); + + if n_symbols == 0 + return + end + + if isempty(obj.constellation) + obj.constellation = unique(d); + end + decision_constellation = obj.constellation(:); + + x = [zeros(floor(obj.order/2),1); x; zeros(obj.order,1)]; + lambda = mu; + mask = ones(obj.order,1); + maincursor_pos = ceil(length(obj.e)/2); + adaption_code = obj.adaptionCode(); + + if mu == 0 || (~obj.dd_mode && ~training) + epochs = 1; + end + + debug_enabled = obj.save_debug; + if debug_enabled + obj.initializeDebug(n_symbols,training); + end + + for epoch = 1:epochs + symbol = 0; + for sample = 1:obj.sps:N + symbol = symbol + 1; + grad = zeros(obj.order,1); + update = zeros(obj.order,1); + weight = 0; + + U = x(obj.order+sample-1:-1:sample); + y(symbol,1) = (obj.e.*mask).' * U; + + if training + d_hat(symbol,1) = d(symbol); + else + [~,symbol_idx] = min(abs(y(symbol) - decision_constellation)); + d_hat(symbol,1) = decision_constellation(symbol_idx); + end + + err = d_hat(symbol) - y(symbol); + + if mu ~= 0 && (training || obj.dd_mode) + switch adaption_code + case 1 + weight = mu / ((U.'*U) + eps); + grad = err * U; + update = grad * weight; + obj.e = obj.e + update; + + case 2 + weight = mu; + grad = err * U; + update = grad * weight; + obj.e = obj.e + update; + + case 3 + denom = lambda + U.' * obj.P * U; + k = (obj.P * U) / denom; + update = k * err; + obj.e = obj.e + update; + obj.P = (1/lambda) * (obj.P - k * (U.' * obj.P)); + end + end + + if debug_enabled && epoch == 1 + obj.debug_struct.error_first_epoch(1,symbol) = err * err'; + end + + if debug_enabled && epoch == epochs + error_power = err * err'; + 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; + obj.debug_struct.update_gradient(1,symbol) = grad.'*grad; + + if training + obj.debug_struct.error_tr(1,symbol) = error_power; + obj.debug_struct.update_tr(1,symbol) = update_power; + else + obj.debug_struct.error_dd(1,symbol) = error_power; + obj.debug_struct.update(1,symbol) = update_power; + end + end + end + end + end + + function optimizeMus(obj,x,d) + [x_opt,d_opt,N_opt] = obj.optimizationSignals(x,d); + + switch obj.adaption_technique + case adaption_method.lms + mu_range = [1e-5, 1e-2]; + case adaption_method.nlms + mu_range = [1e-3, 5e-1]; + case adaption_method.rls + mu_range = [0.98, 0.99999]; + otherwise + builtin("error","FFE_plain:InvalidAdaptionTechnique", ... + "Unsupported FFE adaption technique."); + end + + vars = optimizableVariable("mu_tr",mu_range,"Transform","log"); + if obj.dd_mode + vars = [vars, optimizableVariable("mu_dd",mu_range,"Transform","log")]; + end + + obj.mu_optimization_iter = 0; + fprintf("FFE_plain mu opt uses %d samples / %d symbols\n",N_opt,numel(d_opt)); + obj.mu_optimization = bayesopt(@(p)obj.muObjective(p,x_opt,d_opt),vars, ... + "MaxObjectiveEvaluations",20, ... + "AcquisitionFunctionName","expected-improvement-plus", ... + "IsObjectiveDeterministic",false, ... + "Verbose",0, ... + "PlotFcn",[]); + + obj.mu_tr = obj.mu_optimization.XAtMinObjective.mu_tr; + if obj.dd_mode + obj.mu_dd = obj.mu_optimization.XAtMinObjective.mu_dd; + fprintf("\nFFE_plain mu opt done: mu_tr=%9.3e, mu_dd=%9.3e, BER=%9.3e\n", ... + obj.mu_tr,obj.mu_dd,obj.mu_optimization.MinObjective); + else + fprintf("\nFFE_plain mu opt done: mu_tr=%9.3e, BER=%9.3e\n", ... + obj.mu_tr,obj.mu_optimization.MinObjective); + end + + if obj.plot_mu_optimization + obj.plotMuOptimization(); + end + end + + function objective = muObjective(obj,params,x,d) + old_state = obj.captureObjectiveState(); + cleanup = onCleanup(@()obj.restoreObjectiveState(old_state)); + + obj.save_debug = 0; + obj.resetState(); + N_tr = min(obj.len_tr,numel(x)); + obj.equalize(x,d,params.mu_tr,obj.epochs_tr,N_tr,true,false); + if obj.dd_mode + [signal,~] = obj.equalize(x,d,params.mu_dd,obj.epochs_dd,numel(x),false,false); + else + [signal,~] = obj.equalize(x,d,0,1,numel(x),false,false); + end + + [ber,errors] = obj.berObjective(signal,d); + objective = ber; + if ~isfinite(objective) + objective = inf; + end + + obj.mu_optimization_iter = obj.mu_optimization_iter + 1; + if obj.dd_mode + fprintf("\rFFE_plain mu opt %02d: mu_tr=%9.3e, mu_dd=%9.3e, BER=%9.3e, errors=%d", ... + obj.mu_optimization_iter,params.mu_tr,params.mu_dd,ber,errors); + else + fprintf("\rFFE_plain mu opt %02d: mu_tr=%9.3e, BER=%9.3e, errors=%d", ... + obj.mu_optimization_iter,params.mu_tr,ber,errors); + end + + clear cleanup + end + + function [x_opt,d_opt,N_opt] = optimizationSignals(obj,x,d,opt_len) + if nargin < 4 + opt_len = obj.mu_optimization_len; + end + + N_available = min(numel(x),numel(d) * obj.sps); + + if isempty(opt_len) || opt_len <= 0 || isinf(opt_len) + N_opt = N_available; + else + N_opt = min(N_available,max(obj.len_tr,opt_len)); + end + + N_opt = obj.sps * floor(N_opt / obj.sps); + N_opt = max(obj.sps,N_opt); + + n_symbols = N_opt / obj.sps; + x_opt = x(1:N_opt); + d_opt = d(1:n_symbols); + end + + function [ber,errors] = berObjective(~,signal,d) + M = numel(unique(d)); + mapper = PAMmapper(M,0); + eq_signal_sd = Signal(signal); + eq_signal_hd = mapper.quantize(eq_signal_sd); + tx_symbols = Signal(d); + 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",1000, ... + "skip_end",0, ... + "returnErrorLocation",1); + end + + function plotMuOptimization(obj) + if isempty(obj.mu_optimization) + return + end + + X = obj.mu_optimization.XTrace; + objective = obj.mu_optimization.ObjectiveTrace; + objective = objective(:); + valid = isfinite(objective); + + if isempty(X) || ~any(valid) + return + end + + var_names = X.Properties.VariableNames; + n_vars = numel(var_names); + eval_idx = (1:numel(objective)).'; + objective_plot = obj.positiveObjectiveForLogPlot(objective); + best_plot = obj.positiveObjectiveForLogPlot(cummin(objective)); + + figure(obj.mu_optimization_fignum); + clf; + t = tiledlayout(2,2,"TileSpacing","compact","Padding","compact"); + title(t,"FFE_plain Bayesian mu optimization"); + + nexttile; + h_candidate = semilogy(eval_idx,objective_plot,"o-","DisplayName","candidate"); + obj.addOptimizationDataTips(h_candidate,X,objective,objective_plot,eval_idx,var_names); + hold on; + h_best_trace = semilogy(eval_idx,best_plot,"k-","LineWidth",1.2,"DisplayName","best so far"); + h_best_trace.DataTipTemplate.DataTipRows(end+1) = dataTipTextRow("Evaluation",eval_idx); + h_best_trace.DataTipTemplate.DataTipRows(end+1) = dataTipTextRow("Best BER",best_plot); + grid on; + xlabel("Evaluation"); + ylabel("BER"); + legend("Location","best"); + + if n_vars < 2 + return + end + + pairs = nchoosek(1:n_vars,2); + n_pair_plots = min(size(pairs,1),3); + [~,best_idx] = min(objective); + + for pair_idx = 1:n_pair_plots + nexttile; + x_name = var_names{pairs(pair_idx,1)}; + y_name = var_names{pairs(pair_idx,2)}; + x_data = X.(x_name); + y_data = X.(y_name); + c_data = log10(objective_plot); + + h_scatter = scatter(log10(x_data),log10(y_data),35,c_data,"filled"); + obj.addOptimizationDataTips(h_scatter,X,objective,objective_plot,eval_idx,var_names); + hold on; + h_best = plot(log10(x_data(best_idx)),log10(y_data(best_idx)),"kp", ... + "MarkerSize",12, ... + "MarkerFaceColor","y", ... + "DisplayName","best"); + obj.addOptimizationDataTips(h_best,X(best_idx,:),objective(best_idx),objective_plot(best_idx),eval_idx(best_idx),var_names); + grid on; + xlabel("log10(" + string(x_name) + ")"); + ylabel("log10(" + string(y_name) + ")"); + cb = colorbar; + cb.Label.String = "log10(BER)"; + title(string(x_name) + " vs " + string(y_name)); + end + end + + function addOptimizationDataTips(~,plot_handle,X,objective,objective_plot,eval_idx,var_names) + plot_handle.DataTipTemplate.DataTipRows(end+1) = dataTipTextRow("Evaluation",eval_idx); + plot_handle.DataTipTemplate.DataTipRows(end+1) = dataTipTextRow("BER",objective); + plot_handle.DataTipTemplate.DataTipRows(end+1) = dataTipTextRow("BER shown",objective_plot); + + for var_idx = 1:numel(var_names) + var_name = var_names{var_idx}; + plot_handle.DataTipTemplate.DataTipRows(end+1) = dataTipTextRow(var_name,X.(var_name)); + end + end + + function objective_plot = positiveObjectiveForLogPlot(~,objective) + objective_plot = objective; + positive_values = objective(isfinite(objective) & objective > 0); + + if isempty(positive_values) + floor_value = 1e-12; + else + floor_value = min(positive_values) / 10; + end + + objective_plot(~isfinite(objective_plot) | objective_plot <= 0) = floor_value; + end + + function state = captureObjectiveState(obj) + state.e = obj.e; + state.e_tr = obj.e_tr; + state.error = obj.error; + state.P = obj.P; + state.save_debug = obj.save_debug; + state.debug_struct = obj.debug_struct; + end + + function restoreObjectiveState(obj,state) + obj.e = state.e; + obj.e_tr = state.e_tr; + obj.error = state.error; + obj.P = state.P; + obj.save_debug = state.save_debug; + obj.debug_struct = state.debug_struct; + end + end + + methods (Access = private) + function resetState(obj) + obj.e = zeros(obj.order,1); + obj.e_tr = zeros(obj.order,1); + obj.error = 0; + obj.P = (1/0.05) * eye(obj.order); + obj.debug_struct = struct(); + end + + function N = validSampleLength(obj,N,x,d) + N_available = min(numel(x),numel(d) * obj.sps); + N = min(N,N_available); + N = obj.sps * floor(N / obj.sps); + N = max(0,N); + end + + function adaption_code = adaptionCode(obj) + 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_plain:InvalidAdaptionTechnique", ... + "Unsupported FFE adaption technique."); + end + end + + function initializeDebug(obj,n_symbols,training) + obj.debug_struct.error = NaN(1,n_symbols); + obj.debug_struct.error_first_epoch = NaN(1,n_symbols); + obj.debug_struct.main_cursor = NaN(1,n_symbols); + obj.debug_struct.mu_nlms = NaN(1,n_symbols); + obj.debug_struct.update_gradient = NaN(1,n_symbols); + + if training + obj.debug_struct.error_tr = NaN(1,n_symbols); + obj.debug_struct.update_tr = NaN(1,n_symbols); + else + obj.debug_struct.error_dd = NaN(1,n_symbols); + obj.debug_struct.update = NaN(1,n_symbols); + end + end + end +end diff --git a/Classes/04_DSP/Equalizer/MPI_reduction/FFE_A1.m b/Classes/04_DSP/Equalizer/MPI_reduction/FFE_A1.m new file mode 100644 index 0000000..bec0420 --- /dev/null +++ b/Classes/04_DSP/Equalizer/MPI_reduction/FFE_A1.m @@ -0,0 +1,311 @@ +classdef FFE_A1 < handle + %FFE_A1 Feed-forward equalizer with A1 moving-average input suppression. + % + % This class is intentionally narrow: it contains the normal adaptive FFE + % flow plus the A1 block-rate moving-average input correction. A2, + % DC-tracking, and optimizer logic are kept out of this file. + + properties + sps + order + e + e_tr + error + + len_tr + mu_tr + epochs_tr + + adaption_technique + dd_mode + mu_dd + epochs_dd + dd_len_fraction + + % A1 moving-average input suppression + dc_avg_bufferlength_a1 + dc_avg_update_blocklength_a1 + dc_smoothing_a1 + + P + constellation + + decide + + save_debug = 0 + debug_struct + end + + methods + function obj = FFE_A1(options) + arguments(Input) + options.sps = 2 + options.order = 15 + + options.len_tr = 4096 + options.mu_tr = 0 + options.epochs_tr = 5 + + options.adaption_technique adaption_method = adaption_method.lms + options.dd_mode = 1 + options.mu_dd = 1e-5 + options.epochs_dd = 5 + options.dd_len_fraction = 1 + + options.dc_avg_bufferlength_a1 = 0 + options.dc_avg_update_blocklength_a1 = 0 + options.dc_smoothing_a1 = 0 + + options.decide = false + options.save_debug = 0 + end + + fn = fieldnames(options); + for n = 1:numel(fn) + obj.(fn{n}) = options.(fn{n}); + end + + assert(obj.dc_avg_bufferlength_a1 >= 0); + assert(obj.dc_avg_update_blocklength_a1 >= 0); + + obj.e = zeros(obj.order,1); + obj.error = 0; + 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); + end + + function [X,Noi] = process(obj,X,D) + X = X.normalize("mode","rms"); + + obj.constellation = unique(D.signal); + + delta = 0.05; + obj.P = (1/delta) * eye(obj.order); + + training = true; + showviz = false; + obj.equalize(X.signal,D.signal,obj.mu_tr,obj.epochs_tr,obj.len_tr,training,showviz); + obj.e_tr = obj.e; + + n = X.length; + training = false; + if obj.dd_mode + [signal,decision] = obj.equalize(X.signal,D.signal,obj.mu_dd,obj.epochs_dd,n,training,showviz); + else + [signal,decision] = obj.equalize(X.signal,D.signal,0,1,n,training,showviz); + end + + if obj.decide + X.signal = decision; + else + X.signal = signal; + end + + X.fs = D.fs; + X = X.logbookentry([num2str(obj.order),' tap FFE A1']); + + Noi = X - D; + end + + function [y,d_hat] = equalize(obj,x,d,mu,epochs,N,training,showviz) + arguments + obj + x + d + mu + epochs + N + training + showviz %#ok + end + + 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); + constellation = obj.constellation; + + adaption_code = obj.adaptionCode(); + adaption_is_rls = adaption_code == 3; + + mask = ones(obj.order,1); + maincursor_pos = ceil(length(obj.e)/2); + grad = 0; + weight = 0; + update = 0; + + if mu == 0 || (~obj.dd_mode && ~training) + epochs = 1; + end + + dc_avg_enabled = obj.dc_avg_bufferlength_a1 > 1; + if dc_avg_enabled + dc_avg_buffer = NaN(obj.dc_avg_bufferlength_a1,1); + dc_avg_buffer_pos = 0; + dc_avg_sum = 0; + dc_avg_valid_count = 0; + dc_avg_est = 0; + dc_avg_update_blocklength = obj.dc_avg_update_blocklength_a1; + if dc_avg_update_blocklength <= 0 + dc_avg_update_blocklength = obj.dc_avg_bufferlength_a1; + end + dc_avg_update_blocklength = max(1,floor(dc_avg_update_blocklength)); + dc_avg_input_offset = floor(obj.order/2); + % Hardware-like A1 is causal; dc_smoothing_a1 is reserved for offline variants. + end + + debug_enabled = obj.save_debug; + if debug_enabled + obj.debug_struct.error = NaN(1,n_symbols); + obj.debug_struct.error_first_epoch = NaN(1,n_symbols); + obj.debug_struct.main_cursor = NaN(1,n_symbols); + obj.debug_struct.mu_nlms = NaN(1,n_symbols); + obj.debug_struct.update_gradient = NaN(1,n_symbols); + obj.debug_struct.dc_avg_offset = NaN(1,n_symbols); + + if training + obj.debug_struct.error_tr = NaN(1,n_symbols); + obj.debug_struct.update_tr = NaN(1,n_symbols); + else + obj.debug_struct.error_dd = NaN(1,n_symbols); + obj.debug_struct.update = NaN(1,n_symbols); + end + end + + for epoch = 1:epochs + symbol = 0; + for sample = 1:obj.sps:N + symbol = symbol + 1; + dc_avg_offset = 0; + + U = x(obj.order+sample-1:-1:sample); + if dc_avg_enabled + dc_avg_offset = dc_avg_est; + U = U - dc_avg_offset; + end + + y(symbol,1) = (obj.e.*mask).' * U; + + if dc_avg_enabled + raw_sample_idx = dc_avg_input_offset + sample; + [dc_avg_buffer,dc_avg_buffer_pos,dc_avg_sum,dc_avg_valid_count,dc_avg_est] = ... + obj.updateA1Buffer(dc_avg_buffer,dc_avg_buffer_pos,dc_avg_sum, ... + dc_avg_valid_count,dc_avg_est,x(raw_sample_idx),symbol, ... + dc_avg_update_blocklength); + end + + if training + d_hat(symbol,1) = d(symbol); + else + [~,symbol_idx] = min(abs(y(symbol) - constellation)); + d_hat(symbol,1) = constellation(symbol_idx); + end + + err(symbol) = d_hat(symbol) - y(symbol); + + if training || obj.dd_mode + switch adaption_code + case 1 + normU = (U.'*U) + eps; + weight = mu / normU; + grad = err(symbol) * U; + update = grad * weight; + obj.e = obj.e + update; + + case 2 + weight = mu; + grad = err(symbol) * U; + update = grad * weight; + obj.e = obj.e + update; + + case 3 + denom = lambda + U.' * obj.P * U; + k = (obj.P * U) / denom; + update = k * err(symbol); + obj.e = obj.e + update; + obj.P = (1/lambda) * (obj.P - k * (U.' * obj.P)); + end + end + + if debug_enabled && epoch == 1 + obj.debug_struct.error_first_epoch(1,symbol) = err(symbol) * err(symbol)'; + end + + if debug_enabled && epoch == epochs + error_power = err(symbol) * err(symbol)'; + 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; + obj.debug_struct.update_gradient(1,symbol) = grad.'*grad; + obj.debug_struct.dc_avg_offset(1,symbol) = dc_avg_offset; + + if training + obj.debug_struct.error_tr(1,symbol) = error_power; + obj.debug_struct.update_tr(1,symbol) = update_power; + else + obj.debug_struct.error_dd(1,symbol) = error_power; + obj.debug_struct.update(1,symbol) = update_power; + end + end + end + end + + if ~adaption_is_rls + obj.P = []; + end + end + end + + methods (Access = private) + function adaption_code = adaptionCode(obj) + 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_A1:InvalidAdaptionTechnique", ... + "Unsupported FFE adaption technique."); + end + end + end + + methods (Static, Access = private) + function [buffer,buffer_pos,buffer_sum,valid_count,estimate] = ... + updateA1Buffer(buffer,buffer_pos,buffer_sum,valid_count,estimate, ... + new_value,symbol,update_blocklength) + + buffer_pos = buffer_pos + 1; + if buffer_pos > numel(buffer) + buffer_pos = 1; + end + + old_value = buffer(buffer_pos); + if isfinite(old_value) + buffer_sum = buffer_sum - old_value; + valid_count = valid_count - 1; + end + + if isfinite(new_value) + buffer(buffer_pos) = new_value; + buffer_sum = buffer_sum + new_value; + valid_count = valid_count + 1; + else + buffer(buffer_pos) = NaN; + end + + if mod(symbol,update_blocklength) == 0 + if valid_count == 0 + estimate = 0; + else + estimate = buffer_sum / valid_count; + end + end + end + end +end diff --git a/Classes/04_DSP/Equalizer/MPI_reduction/FFE_A2Residual.m b/Classes/04_DSP/Equalizer/MPI_reduction/FFE_A2Residual.m new file mode 100644 index 0000000..12754cd --- /dev/null +++ b/Classes/04_DSP/Equalizer/MPI_reduction/FFE_A2Residual.m @@ -0,0 +1,347 @@ +classdef FFE_A2Residual < FFE_plain + %FFE_A2Residual Plain FFE with A2 residual MPI suppression. + % + % This class reuses FFE_plain for process flow and BER-based mu + % optimization. Only equalize is specialized to add the paper-style A2 + % residual correction: average y_raw - d_hat per PAM level, subtract the + % level-weighted residual estimate, then redo the DD decision. + + properties + dc_level_avg_bufferlength_a2 + dc_level_update_blocklength_a2 + dc_smoothing_a2 + dc_level_weights_a2 + end + + methods + function obj = FFE_A2Residual(options) + arguments(Input) + options.sps = 2 + options.order = 15 + + options.len_tr = 4096 + options.mu_tr = 0 + options.epochs_tr = 5 + + options.adaption_technique adaption_method = adaption_method.lms + options.dd_mode = 1 + options.mu_dd = 1e-5 + options.epochs_dd = 5 + options.dd_len_fraction = 1 + + 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.decide = false + + options.save_debug = 0 + options.optmize_mus = 0 + options.mu_optimization_len = 2^15 + options.plot_mu_optimization = 0 + options.mu_optimization_fignum = 3010 + end + + obj@FFE_plain( ... + "sps",options.sps, ... + "order",options.order, ... + "len_tr",options.len_tr, ... + "mu_tr",options.mu_tr, ... + "epochs_tr",options.epochs_tr, ... + "adaption_technique",options.adaption_technique, ... + "dd_mode",options.dd_mode, ... + "mu_dd",options.mu_dd, ... + "epochs_dd",options.epochs_dd, ... + "dd_len_fraction",options.dd_len_fraction, ... + "decide",options.decide, ... + "save_debug",options.save_debug, ... + "optmize_mus",options.optmize_mus, ... + "mu_optimization_len",options.mu_optimization_len, ... + "plot_mu_optimization",options.plot_mu_optimization, ... + "mu_optimization_fignum",options.mu_optimization_fignum); + + obj.dc_level_avg_bufferlength_a2 = floor(options.dc_level_avg_bufferlength_a2); + obj.dc_level_update_blocklength_a2 = floor(options.dc_level_update_blocklength_a2); + obj.dc_smoothing_a2 = min(max(options.dc_smoothing_a2,0),1); + obj.dc_level_weights_a2 = options.dc_level_weights_a2; + + assert(obj.dc_level_avg_bufferlength_a2 >= 0); + assert(obj.dc_level_update_blocklength_a2 >= 0); + end + + function [y,d_hat] = equalize(obj,x,d,mu,epochs,N,training,showviz) + arguments + obj + x + d + mu + epochs + N + training + showviz + end + + unused_showviz = showviz; %#ok + x = x(:); + d = d(:); + N = obj.validSampleLength(N,x,d); + n_symbols = N / obj.sps; + y = zeros(n_symbols,1); + d_hat = zeros(n_symbols,1); + + if n_symbols == 0 + return + end + + if isempty(obj.constellation) + obj.constellation = unique(d); + end + decision_constellation = obj.constellation(:); + + x = [zeros(floor(obj.order/2),1); x; zeros(obj.order,1)]; + lambda = mu; + mask = ones(obj.order,1); + maincursor_pos = ceil(length(obj.e)/2); + adaption_code = obj.adaptionCode(); + + if mu == 0 || (~obj.dd_mode && ~training) + epochs = 1; + end + + dc_level_enabled = obj.dc_level_avg_bufferlength_a2 > 1 && any(obj.dc_level_weights_a2(:) ~= 0); + if dc_level_enabled + [dc_level_weight_by_level,n_levels] = obj.expandedLevelWeights(decision_constellation); + dc_level_buffer_len = obj.dc_level_avg_bufferlength_a2; + dc_level_err_buffer = NaN(n_levels,dc_level_buffer_len); + 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 = obj.dc_level_update_blocklength_a2; + if dc_level_update_blocklength <= 0 + dc_level_update_blocklength = dc_level_buffer_len; + end + dc_level_update_blocklength = max(1,floor(dc_level_update_blocklength)); + dc_level_window_future_fraction = obj.dc_smoothing_a2; %#ok + end + + debug_enabled = obj.save_debug; + if debug_enabled + obj.initializeDebug(n_symbols,training); + end + + for epoch = 1:epochs + symbol = 0; + for sample = 1:obj.sps:N + symbol = symbol + 1; + grad = zeros(obj.order,1); + update = zeros(obj.order,1); + weight = 0; + dc_level_mpi_est = 0; + dc_level_weight = 0; + dc_level_valid_count = 0; + + U = x(obj.order+sample-1:-1:sample); + y(symbol,1) = (obj.e.*mask).' * U; + y_raw = y(symbol); + + if training + d_hat(symbol,1) = d(symbol); + [~,symbol_idx] = min(abs(d_hat(symbol) - decision_constellation)); + dc_level_decision_level = decision_constellation(symbol_idx); + else + [~,symbol_idx] = min(abs(y_raw - decision_constellation)); + d_hat(symbol,1) = decision_constellation(symbol_idx); + dc_level_decision_level = decision_constellation(symbol_idx); + end + + dc_level_symbol_idx = symbol_idx; + if dc_level_enabled + dc_level_mpi_est = dc_level_mpi_est_by_level(dc_level_symbol_idx); + dc_level_valid_count = dc_level_valid_count_by_level(dc_level_symbol_idx); + dc_level_weight = dc_level_weight_by_level(dc_level_symbol_idx) * ... + min(dc_level_valid_count / dc_level_buffer_len,1); + + mpi_err = y_raw - d_hat(symbol); + y(symbol,1) = y_raw - dc_level_weight * dc_level_mpi_est; + + if training + d_hat(symbol,1) = d(symbol); + else + [~,symbol_idx] = min(abs(y(symbol) - decision_constellation)); + d_hat(symbol,1) = decision_constellation(symbol_idx); + dc_level_decision_level = decision_constellation(symbol_idx); + end + + [dc_level_err_buffer,dc_level_err_buffer_pos_by_level, ... + dc_level_err_sum_by_level,dc_level_buffer_valid_count_by_level, ... + dc_level_mpi_est_by_level,dc_level_valid_count_by_level] = ... + obj.updateResidualBuffer( ... + dc_level_err_buffer,dc_level_err_buffer_pos_by_level, ... + dc_level_err_sum_by_level,dc_level_buffer_valid_count_by_level, ... + dc_level_mpi_est_by_level,dc_level_valid_count_by_level, ... + mpi_err,dc_level_symbol_idx,symbol,dc_level_update_blocklength); + end + + err = d_hat(symbol) - y(symbol); + + if mu ~= 0 && (training || obj.dd_mode) + switch adaption_code + case 1 + weight = mu / ((U.'*U) + eps); + grad = err * U; + update = grad * weight; + obj.e = obj.e + update; + + case 2 + weight = mu; + grad = err * U; + update = grad * weight; + obj.e = obj.e + update; + + case 3 + denom = lambda + U.' * obj.P * U; + k = (obj.P * U) / denom; + update = k * err; + obj.e = obj.e + update; + obj.P = (1/lambda) * (obj.P - k * (U.' * obj.P)); + end + end + + if debug_enabled && epoch == 1 + obj.debug_struct.error_first_epoch(1,symbol) = err * err'; + end + + if debug_enabled && epoch == epochs + error_power = err * err'; + 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; + obj.debug_struct.update_gradient(1,symbol) = grad.'*grad; + obj.debug_struct.dc_level_mpi_est(1,symbol) = dc_level_mpi_est; + obj.debug_struct.dc_level_weight(1,symbol) = dc_level_weight; + obj.debug_struct.dc_level_valid_count(1,symbol) = dc_level_valid_count; + obj.debug_struct.dc_level_symbol_idx(1,symbol) = dc_level_symbol_idx; + obj.debug_struct.dc_level_decision_level(1,symbol) = dc_level_decision_level; + obj.debug_struct.dc_level_y_raw(1,symbol) = y_raw; + + if training + obj.debug_struct.error_tr(1,symbol) = error_power; + obj.debug_struct.update_tr(1,symbol) = update_power; + else + obj.debug_struct.error_dd(1,symbol) = error_power; + obj.debug_struct.update(1,symbol) = update_power; + end + end + end + end + end + end + + methods (Access = private) + function N = validSampleLength(obj,N,x,d) + N_available = min(numel(x),numel(d) * obj.sps); + N = min(N,N_available); + N = obj.sps * floor(N / obj.sps); + N = max(0,N); + end + + function adaption_code = adaptionCode(obj) + 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_A2Residual:InvalidAdaptionTechnique", ... + "Unsupported FFE adaption technique."); + end + end + + function [weights,n_levels] = expandedLevelWeights(obj,decision_constellation) + n_levels = numel(decision_constellation); + if isscalar(obj.dc_level_weights_a2) + weights = repmat(obj.dc_level_weights_a2,n_levels,1); + elseif numel(obj.dc_level_weights_a2) == n_levels + weights = obj.dc_level_weights_a2(:); + else + builtin("error","FFE_A2Residual:InvalidDCLevelWeights", ... + "dc_level_weights_a2 must be scalar or have one entry per constellation level."); + end + end + + function initializeDebug(obj,n_symbols,training) + obj.debug_struct.error = NaN(1,n_symbols); + obj.debug_struct.error_first_epoch = NaN(1,n_symbols); + obj.debug_struct.main_cursor = NaN(1,n_symbols); + obj.debug_struct.mu_nlms = NaN(1,n_symbols); + obj.debug_struct.update_gradient = NaN(1,n_symbols); + obj.debug_struct.dc_level_mpi_est = NaN(1,n_symbols); + obj.debug_struct.dc_level_weight = NaN(1,n_symbols); + obj.debug_struct.dc_level_valid_count = NaN(1,n_symbols); + obj.debug_struct.dc_level_symbol_idx = NaN(1,n_symbols); + obj.debug_struct.dc_level_decision_level = NaN(1,n_symbols); + obj.debug_struct.dc_level_y_raw = NaN(1,n_symbols); + + if training + obj.debug_struct.error_tr = NaN(1,n_symbols); + obj.debug_struct.update_tr = NaN(1,n_symbols); + else + obj.debug_struct.error_dd = NaN(1,n_symbols); + obj.debug_struct.update = NaN(1,n_symbols); + end + end + end + + methods (Static, Access = private) + function [buffer,buffer_pos_by_level,err_sum_by_level,buffer_valid_count_by_level, ... + mpi_est_by_level,valid_count_by_level] = updateResidualBuffer( ... + buffer,buffer_pos_by_level,err_sum_by_level,buffer_valid_count_by_level, ... + mpi_est_by_level,valid_count_by_level,new_value,symbol_idx,symbol, ... + update_blocklength) + + buffer_len = size(buffer,2); + buffer_pos = buffer_pos_by_level(symbol_idx) + 1; + if buffer_pos > buffer_len + buffer_pos = 1; + end + buffer_pos_by_level(symbol_idx) = buffer_pos; + + old_value = buffer(symbol_idx,buffer_pos); + if isfinite(old_value) + err_sum_by_level(symbol_idx) = err_sum_by_level(symbol_idx) - old_value; + buffer_valid_count_by_level(symbol_idx) = buffer_valid_count_by_level(symbol_idx) - 1; + end + + if isfinite(new_value) + buffer(symbol_idx,buffer_pos) = new_value; + err_sum_by_level(symbol_idx) = err_sum_by_level(symbol_idx) + new_value; + buffer_valid_count_by_level(symbol_idx) = buffer_valid_count_by_level(symbol_idx) + 1; + else + buffer(symbol_idx,buffer_pos) = NaN; + end + + if mod(symbol,update_blocklength) == 0 + if update_blocklength == 1 + valid_count_by_level(symbol_idx) = buffer_valid_count_by_level(symbol_idx); + if valid_count_by_level(symbol_idx) == 0 + mpi_est_by_level(symbol_idx) = 0; + else + mpi_est_by_level(symbol_idx) = ... + err_sum_by_level(symbol_idx) / valid_count_by_level(symbol_idx); + end + else + valid_count_by_level = buffer_valid_count_by_level; + has_valid = valid_count_by_level > 0; + mpi_est_by_level(:) = 0; + mpi_est_by_level(has_valid) = ... + err_sum_by_level(has_valid) ./ valid_count_by_level(has_valid); + end + end + end + end +end diff --git a/Classes/04_DSP/Equalizer/MPI_reduction/FFE_A2TrackedLevels.m b/Classes/04_DSP/Equalizer/MPI_reduction/FFE_A2TrackedLevels.m new file mode 100644 index 0000000..305b2ba --- /dev/null +++ b/Classes/04_DSP/Equalizer/MPI_reduction/FFE_A2TrackedLevels.m @@ -0,0 +1,361 @@ +classdef FFE_A2TrackedLevels < handle + %FFE_A2TrackedLevels FFE with A2 tracked-level decision thresholds. + % + % This class is intentionally narrow: it contains the normal adaptive FFE + % flow plus the A2 tracked-level decision logic. The output sample y stays + % unchanged; only the decision constellation is moved by the tracked + % level-wise offsets. + + properties + sps + order + e + e_tr + error + + len_tr + mu_tr + epochs_tr + + adaption_technique + dd_mode + mu_dd + epochs_dd + dd_len_fraction + + % A2 tracked-level decision + dc_level_avg_bufferlength_a2 + dc_level_update_blocklength_a2 + dc_smoothing_a2 + dc_level_weights_a2 + + P + constellation + + decide + + save_debug = 0 + debug_struct + end + + methods + function obj = FFE_A2TrackedLevels(options) + arguments(Input) + options.sps = 2 + options.order = 15 + + options.len_tr = 4096 + options.mu_tr = 0 + options.epochs_tr = 5 + + options.adaption_technique adaption_method = adaption_method.lms + options.dd_mode = 1 + options.mu_dd = 1e-5 + options.epochs_dd = 5 + options.dd_len_fraction = 1 + + 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.decide = false + options.save_debug = 0 + end + + fn = fieldnames(options); + for n = 1:numel(fn) + obj.(fn{n}) = options.(fn{n}); + end + + assert(obj.dc_level_avg_bufferlength_a2 >= 0); + assert(obj.dc_level_update_blocklength_a2 >= 0); + + obj.e = zeros(obj.order,1); + obj.error = 0; + 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); + end + + function [X,Noi] = process(obj,X,D) + X = X.normalize("mode","rms"); + + obj.constellation = unique(D.signal); + + delta = 0.05; + obj.P = (1/delta) * eye(obj.order); + + training = true; + showviz = false; + obj.equalize(X.signal,D.signal,obj.mu_tr,obj.epochs_tr,obj.len_tr,training,showviz); + obj.e_tr = obj.e; + + n = X.length; + training = false; + if obj.dd_mode + [signal,decision] = obj.equalize(X.signal,D.signal,obj.mu_dd,obj.epochs_dd,n,training,showviz); + else + [signal,decision] = obj.equalize(X.signal,D.signal,0,1,n,training,showviz); + end + + if obj.decide + X.signal = decision; + else + X.signal = signal; + end + + X.fs = D.fs; + X = X.logbookentry([num2str(obj.order),' tap FFE A2 tracked levels']); + + Noi = X - D; + end + + function [y,d_hat] = equalize(obj,x,d,mu,epochs,N,training,showviz) + arguments + obj + x + d + mu + epochs + N + training + showviz %#ok + end + + 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); + constellation = obj.constellation; + + adaption_code = obj.adaptionCode(); + adaption_is_rls = adaption_code == 3; + + mask = ones(obj.order,1); + maincursor_pos = ceil(length(obj.e)/2); + grad = 0; + weight = 0; + update = 0; + + if mu == 0 || (~obj.dd_mode && ~training) + epochs = 1; + end + + dc_level_enabled = obj.dc_level_avg_bufferlength_a2 > 1 && any(obj.dc_level_weights_a2(:) ~= 0); + if dc_level_enabled + if isempty(constellation) + builtin("error","FFE_A2TrackedLevels:MissingConstellation", ... + "A2 tracked-level decision requires obj.constellation to be set."); + end + + 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 + dc_level_weight_by_level = obj.dc_level_weights_a2(:); + else + builtin("error","FFE_A2TrackedLevels:InvalidDCLevelWeights", ... + "dc_level_weights_a2 must be scalar or have one entry per constellation level."); + end + + dc_level_buffer_len = obj.dc_level_avg_bufferlength_a2; + dc_level_buffer = NaN(n_levels,dc_level_buffer_len); + dc_level_buffer_pos_by_level = zeros(n_levels,1); + dc_level_sum_by_level = zeros(n_levels,1); + dc_level_buffer_valid_count_by_level = zeros(n_levels,1); + dc_level_offset_by_level = zeros(n_levels,1); + dc_level_valid_count_by_level = zeros(n_levels,1); + dc_level_update_blocklength = obj.dc_level_update_blocklength_a2; + if dc_level_update_blocklength <= 0 + dc_level_update_blocklength = dc_level_buffer_len; + end + dc_level_update_blocklength = max(1,floor(dc_level_update_blocklength)); + dc_level_window_future_fraction = obj.dc_smoothing_a2; %#ok + end + + debug_enabled = obj.save_debug; + if debug_enabled + obj.debug_struct.error = NaN(1,n_symbols); + obj.debug_struct.error_first_epoch = NaN(1,n_symbols); + obj.debug_struct.main_cursor = NaN(1,n_symbols); + obj.debug_struct.mu_nlms = NaN(1,n_symbols); + obj.debug_struct.update_gradient = NaN(1,n_symbols); + obj.debug_struct.dc_level_mpi_est = NaN(1,n_symbols); + obj.debug_struct.dc_level_weight = NaN(1,n_symbols); + obj.debug_struct.dc_level_valid_count = NaN(1,n_symbols); + obj.debug_struct.dc_level_symbol_idx = NaN(1,n_symbols); + obj.debug_struct.dc_level_decision_level = NaN(1,n_symbols); + obj.debug_struct.dc_level_y_raw = NaN(1,n_symbols); + + if training + obj.debug_struct.error_tr = NaN(1,n_symbols); + obj.debug_struct.update_tr = NaN(1,n_symbols); + else + obj.debug_struct.error_dd = NaN(1,n_symbols); + obj.debug_struct.update = NaN(1,n_symbols); + end + end + + for epoch = 1:epochs + symbol = 0; + for sample = 1:obj.sps:N + symbol = symbol + 1; + dc_level_offset = 0; + dc_level_weight = 0; + dc_level_valid_count = 0; + + U = x(obj.order+sample-1:-1:sample); + y(symbol,1) = (obj.e.*mask).' * U; + y_raw = y(symbol); + + if training + d_hat(symbol,1) = d(symbol); + [~,symbol_idx] = min(abs(d_hat(symbol) - constellation)); + dc_level_decision_level = constellation(symbol_idx); + else + if dc_level_enabled + dc_level_ramp_by_level = min(dc_level_valid_count_by_level / dc_level_buffer_len,1); + dc_level_decision_constellation = constellation + ... + dc_level_weight_by_level .* dc_level_ramp_by_level .* dc_level_offset_by_level; + [~,symbol_idx] = min(abs(y_raw - dc_level_decision_constellation)); + dc_level_decision_level = dc_level_decision_constellation(symbol_idx); + else + [~,symbol_idx] = min(abs(y_raw - constellation)); + dc_level_decision_level = constellation(symbol_idx); + end + d_hat(symbol,1) = constellation(symbol_idx); + end + + dc_level_symbol_idx = symbol_idx; + if dc_level_enabled + dc_level_offset = dc_level_offset_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 / dc_level_buffer_len,1); + + [dc_level_buffer,dc_level_buffer_pos_by_level,dc_level_sum_by_level, ... + dc_level_buffer_valid_count_by_level,dc_level_offset_by_level, ... + dc_level_valid_count_by_level] = obj.updateTrackedLevelBuffer( ... + dc_level_buffer,dc_level_buffer_pos_by_level,dc_level_sum_by_level, ... + dc_level_buffer_valid_count_by_level,dc_level_offset_by_level, ... + dc_level_valid_count_by_level,y_raw,symbol_idx,symbol, ... + dc_level_update_blocklength,constellation); + end + + err(symbol) = d_hat(symbol) - y(symbol); + + if training || obj.dd_mode + switch adaption_code + case 1 + normU = (U.'*U) + eps; + weight = mu / normU; + grad = err(symbol) * U; + update = grad * weight; + obj.e = obj.e + update; + + case 2 + weight = mu; + grad = err(symbol) * U; + update = grad * weight; + obj.e = obj.e + update; + + case 3 + denom = lambda + U.' * obj.P * U; + k = (obj.P * U) / denom; + update = k * err(symbol); + obj.e = obj.e + update; + obj.P = (1/lambda) * (obj.P - k * (U.' * obj.P)); + end + end + + if debug_enabled && epoch == 1 + obj.debug_struct.error_first_epoch(1,symbol) = err(symbol) * err(symbol)'; + end + + if debug_enabled && epoch == epochs + error_power = err(symbol) * err(symbol)'; + 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; + obj.debug_struct.update_gradient(1,symbol) = grad.'*grad; + obj.debug_struct.dc_level_mpi_est(1,symbol) = dc_level_offset; + obj.debug_struct.dc_level_weight(1,symbol) = dc_level_weight; + obj.debug_struct.dc_level_valid_count(1,symbol) = dc_level_valid_count; + obj.debug_struct.dc_level_symbol_idx(1,symbol) = dc_level_symbol_idx; + obj.debug_struct.dc_level_decision_level(1,symbol) = dc_level_decision_level; + obj.debug_struct.dc_level_y_raw(1,symbol) = y_raw; + + if training + obj.debug_struct.error_tr(1,symbol) = error_power; + obj.debug_struct.update_tr(1,symbol) = update_power; + else + obj.debug_struct.error_dd(1,symbol) = error_power; + obj.debug_struct.update(1,symbol) = update_power; + end + end + end + end + + if ~adaption_is_rls + obj.P = []; + end + end + end + + methods (Access = private) + function adaption_code = adaptionCode(obj) + 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_A2TrackedLevels:InvalidAdaptionTechnique", ... + "Unsupported FFE adaption technique."); + end + end + end + + methods (Static, Access = private) + function [buffer,buffer_pos_by_level,level_sum_by_level,buffer_valid_count_by_level, ... + offset_by_level,valid_count_by_level] = updateTrackedLevelBuffer( ... + buffer,buffer_pos_by_level,level_sum_by_level,buffer_valid_count_by_level, ... + offset_by_level,valid_count_by_level,new_value,symbol_idx,symbol, ... + update_blocklength,constellation) + + buffer_pos = buffer_pos_by_level(symbol_idx) + 1; + if buffer_pos > size(buffer,2) + buffer_pos = 1; + end + buffer_pos_by_level(symbol_idx) = buffer_pos; + + old_value = buffer(symbol_idx,buffer_pos); + if isfinite(old_value) + level_sum_by_level(symbol_idx) = level_sum_by_level(symbol_idx) - old_value; + buffer_valid_count_by_level(symbol_idx) = buffer_valid_count_by_level(symbol_idx) - 1; + end + + if isfinite(new_value) + buffer(symbol_idx,buffer_pos) = new_value; + level_sum_by_level(symbol_idx) = level_sum_by_level(symbol_idx) + new_value; + buffer_valid_count_by_level(symbol_idx) = buffer_valid_count_by_level(symbol_idx) + 1; + else + buffer(symbol_idx,buffer_pos) = NaN; + end + + if mod(symbol,update_blocklength) == 0 + valid_count_by_level = buffer_valid_count_by_level; + has_valid = valid_count_by_level > 0; + offset_by_level(:) = 0; + offset_by_level(has_valid) = ... + level_sum_by_level(has_valid) ./ valid_count_by_level(has_valid) - ... + constellation(has_valid); + end + end + end +end diff --git a/Classes/04_DSP/Equalizer/MPI_reduction/FFE_DCTracking.m b/Classes/04_DSP/Equalizer/MPI_reduction/FFE_DCTracking.m new file mode 100644 index 0000000..b730b54 --- /dev/null +++ b/Classes/04_DSP/Equalizer/MPI_reduction/FFE_DCTracking.m @@ -0,0 +1,521 @@ +classdef FFE_DCTracking < FFE_plain + %FFE_DCTracking Plain FFE with the DC-tracking offset loop. + % + % This class reuses the plain FFE process idea and BER mu optimization, + % but keeps only the DC-tracking path from the larger FFE class. A1, A2, + % and delayed FFE tap-update buffers are intentionally excluded. + + properties + dc_tracking_mu + dc_tracking_adaptive_enabled + dc_tracking_persistence_gain + dc_tracking_power_exponent + dc_tracking_buffer_len + e_dc + end + + properties (Access = private) + dc_tracking_mu_min = 1e-6 + dc_tracking_mu_max = 3e-1 + dc_tracking_mu_eff_min = 0 + dc_tracking_mu_eff_max = inf + end + + methods + function obj = FFE_DCTracking(options) + arguments(Input) + options.sps = 2 + options.order = 15 + + options.len_tr = 4096 + options.mu_tr = 0 + options.epochs_tr = 5 + + options.adaption_technique adaption_method = adaption_method.lms + options.dd_mode = 1 + options.mu_dd = 1e-5 + options.epochs_dd = 5 + options.dd_len_fraction = 1 + + 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 + + options.decide = false + + options.save_debug = 0 + options.optmize_mus = 0 + options.mu_optimization_len = 2^15 + options.plot_mu_optimization = 0 + options.mu_optimization_fignum = 3010 + end + + obj@FFE_plain( ... + "sps",options.sps, ... + "order",options.order, ... + "len_tr",options.len_tr, ... + "mu_tr",options.mu_tr, ... + "epochs_tr",options.epochs_tr, ... + "adaption_technique",options.adaption_technique, ... + "dd_mode",options.dd_mode, ... + "mu_dd",options.mu_dd, ... + "epochs_dd",options.epochs_dd, ... + "dd_len_fraction",options.dd_len_fraction, ... + "decide",options.decide, ... + "save_debug",options.save_debug, ... + "optmize_mus",options.optmize_mus, ... + "mu_optimization_len",options.mu_optimization_len, ... + "plot_mu_optimization",options.plot_mu_optimization, ... + "mu_optimization_fignum",options.mu_optimization_fignum); + + obj.dc_tracking_mu = options.dc_tracking_mu; + obj.dc_tracking_adaptive_enabled = options.dc_tracking_adaptive_enabled; + obj.dc_tracking_persistence_gain = options.dc_tracking_persistence_gain; + obj.dc_tracking_power_exponent = options.dc_tracking_power_exponent; + obj.dc_tracking_buffer_len = floor(options.dc_tracking_buffer_len); + obj.e_dc = 0; + + assert(obj.dc_tracking_buffer_len >= 0); + end + + function [X,Noi] = process(obj,X,D) + X = X.normalize("mode","rms"); + obj.constellation = unique(D.signal); + obj.resetTrackingState(); + + if obj.optmize_mus + obj.optimizeMus(X.signal,D.signal); + obj.resetTrackingState(); + end + + training = true; + showviz = false; + obj.equalize(X.signal,D.signal,obj.mu_tr,obj.epochs_tr,obj.len_tr,training,showviz); + obj.e_tr = obj.e; + + n = X.length; + training = false; + if obj.dd_mode + [signal,decision] = obj.equalize(X.signal,D.signal,obj.mu_dd,obj.epochs_dd,n,training,showviz); + else + [signal,decision] = obj.equalize(X.signal,D.signal,0,1,n,training,showviz); + end + + if obj.decide + X.signal = decision; + else + X.signal = signal; + end + + X.fs = D.fs; + X = X.logbookentry([num2str(obj.order),' tap FFE DC tracking']); + + Noi = X - D; + end + + function optimizeMus(obj,x,d) + [x_opt,d_opt,N_opt] = obj.optimizationSignals(x,d); + + switch obj.adaption_technique + case adaption_method.lms + mu_range = [1e-5, 1e-2]; + case adaption_method.nlms + mu_range = [1e-3, 5e-1]; + case adaption_method.rls + mu_range = [0.98, 0.99999]; + otherwise + builtin("error","FFE_DCTracking:InvalidAdaptionTechnique", ... + "Unsupported FFE adaption technique."); + end + + vars = optimizableVariable("mu_tr",mu_range,"Transform","log"); + if obj.dd_mode + vars = [vars, optimizableVariable("mu_dd",mu_range,"Transform","log")]; + end + optimize_dc_tracking_mu = obj.dc_tracking_mu ~= 0; + if optimize_dc_tracking_mu + vars = [vars, optimizableVariable("dc_tracking_mu",[1e-5,1e-1],"Transform","log")]; + end + + obj.mu_optimization_iter = 0; + fprintf("FFE_DCTracking mu opt uses %d samples / %d symbols\n",N_opt,numel(d_opt)); + obj.mu_optimization = bayesopt(@(p)obj.muObjective(p,x_opt,d_opt),vars, ... + "MaxObjectiveEvaluations",20, ... + "AcquisitionFunctionName","expected-improvement-plus", ... + "IsObjectiveDeterministic",false, ... + "Verbose",0, ... + "PlotFcn",[]); + + obj.mu_tr = obj.mu_optimization.XAtMinObjective.mu_tr; + if obj.dd_mode + obj.mu_dd = obj.mu_optimization.XAtMinObjective.mu_dd; + end + if optimize_dc_tracking_mu + obj.dc_tracking_mu = obj.mu_optimization.XAtMinObjective.dc_tracking_mu; + end + + if obj.dd_mode && optimize_dc_tracking_mu + fprintf("\nFFE_DCTracking mu opt done: mu_tr=%9.3e, mu_dd=%9.3e, dc_tracking_mu=%9.3e, BER=%9.3e\n", ... + obj.mu_tr,obj.mu_dd,obj.dc_tracking_mu,obj.mu_optimization.MinObjective); + elseif obj.dd_mode + fprintf("\nFFE_DCTracking mu opt done: mu_tr=%9.3e, mu_dd=%9.3e, BER=%9.3e\n", ... + obj.mu_tr,obj.mu_dd,obj.mu_optimization.MinObjective); + elseif optimize_dc_tracking_mu + fprintf("\nFFE_DCTracking mu opt done: mu_tr=%9.3e, dc_tracking_mu=%9.3e, BER=%9.3e\n", ... + obj.mu_tr,obj.dc_tracking_mu,obj.mu_optimization.MinObjective); + else + fprintf("\nFFE_DCTracking mu opt done: mu_tr=%9.3e, BER=%9.3e\n", ... + obj.mu_tr,obj.mu_optimization.MinObjective); + end + + if obj.plot_mu_optimization + obj.plotMuOptimization(); + end + end + + function objective = muObjective(obj,params,x,d) + old_state = obj.captureTrackingObjectiveState(); + cleanup = onCleanup(@()obj.restoreTrackingObjectiveState(old_state)); + + obj.save_debug = 0; + if any(strcmp(params.Properties.VariableNames,"dc_tracking_mu")) + obj.dc_tracking_mu = params.dc_tracking_mu; + end + obj.resetTrackingState(); + + N_tr = min(obj.len_tr,numel(x)); + obj.equalize(x,d,params.mu_tr,obj.epochs_tr,N_tr,true,false); + if obj.dd_mode + [signal,~] = obj.equalize(x,d,params.mu_dd,obj.epochs_dd,numel(x),false,false); + else + [signal,~] = obj.equalize(x,d,0,1,numel(x),false,false); + end + + [ber,errors] = obj.berObjective(signal,d); + objective = ber; + if ~isfinite(objective) + objective = inf; + end + + obj.mu_optimization_iter = obj.mu_optimization_iter + 1; + has_dc_tracking_mu = any(strcmp(params.Properties.VariableNames,"dc_tracking_mu")); + if obj.dd_mode && has_dc_tracking_mu + fprintf("\rFFE_DCTracking mu opt %02d: mu_tr=%9.3e, mu_dd=%9.3e, dc_tracking_mu=%9.3e, BER=%9.3e, errors=%d", ... + obj.mu_optimization_iter,params.mu_tr,params.mu_dd,params.dc_tracking_mu,ber,errors); + elseif obj.dd_mode + fprintf("\rFFE_DCTracking mu opt %02d: mu_tr=%9.3e, mu_dd=%9.3e, BER=%9.3e, errors=%d", ... + obj.mu_optimization_iter,params.mu_tr,params.mu_dd,ber,errors); + elseif has_dc_tracking_mu + fprintf("\rFFE_DCTracking mu opt %02d: mu_tr=%9.3e, dc_tracking_mu=%9.3e, BER=%9.3e, errors=%d", ... + obj.mu_optimization_iter,params.mu_tr,params.dc_tracking_mu,ber,errors); + else + fprintf("\rFFE_DCTracking mu opt %02d: mu_tr=%9.3e, BER=%9.3e, errors=%d", ... + obj.mu_optimization_iter,params.mu_tr,ber,errors); + end + + clear cleanup + end + + function [y,d_hat] = equalize(obj,x,d,mu,epochs,N,training,showviz) + arguments + obj + x + d + mu + epochs + N + training + showviz + end + + unused_showviz = showviz; %#ok + x = x(:); + d = d(:); + N = obj.validSampleLength(N,x,d); + n_symbols = N / obj.sps; + y = zeros(n_symbols,1); + d_hat = zeros(n_symbols,1); + + if n_symbols == 0 + return + end + + if isempty(obj.constellation) + obj.constellation = unique(d); + end + decision_constellation = obj.constellation(:); + + x = [zeros(floor(obj.order/2),1); x; zeros(obj.order,1)]; + lambda = mu; + mask = ones(obj.order,1); + maincursor_pos = ceil(length(obj.e)/2); + adaption_code = obj.adaptionCode(); + adaption_is_rls = adaption_code == 3; + + if mu == 0 || (~obj.dd_mode && ~training) + epochs = 1; + end + + dc_tracking_mu_eff_min = obj.dc_tracking_mu_eff_min; + dc_tracking_mu_eff_max = obj.dc_tracking_mu_eff_max; + 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, ... + obj.dc_tracking_mu_min),obj.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 + 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; + else + dc_tracking_err_buffer = []; + dc_tracking_err_buffer_pos = 0; + dc_tracking_err_sum = 0; + dc_tracking_abs_err_sum = 0; + dc_tracking_valid_count = 0; + end + + debug_enabled = obj.save_debug; + if debug_enabled + obj.initializeTrackingDebug(n_symbols,training); + end + + for epoch = 1:epochs + symbol = 0; + for sample = 1:obj.sps:N + symbol = symbol + 1; + grad = zeros(obj.order,1); + update = zeros(obj.order,1); + weight = 0; + dc_tracking_mu_eff = 0; + + U = x(obj.order+sample-1:-1:sample); + y(symbol,1) = obj.e_dc + (obj.e.*mask).' * U; + + if training + d_hat(symbol,1) = d(symbol); + else + [~,symbol_idx] = min(abs(y(symbol) - decision_constellation)); + d_hat(symbol,1) = decision_constellation(symbol_idx); + end + + err = d_hat(symbol) - y(symbol); + + switch adaption_code + case 1 + weight = mu / ((U.'*U) + eps); + grad = err * U; + update = grad * weight; + + case 2 + weight = mu; + grad = err * U; + update = grad * weight; + + case 3 + denom = lambda + U.' * obj.P * U; + k = (obj.P * U) / denom; + update = k * err; + end + + obj.e = obj.e + update; + + if adaption_is_rls + obj.P = (1/lambda) * (obj.P - k * (U.' * obj.P)); + end + + if dc_tracking_enabled + [obj.e_dc,dc_tracking_mu_eff,dc_tracking_err_buffer, ... + dc_tracking_err_buffer_pos,dc_tracking_err_sum, ... + dc_tracking_abs_err_sum,dc_tracking_valid_count] = ... + obj.updateDcTracking( ... + err,dc_tracking_mu,dc_tracking_base_mu_eff, ... + dc_tracking_persistence_gain,dc_tracking_use_persistence, ... + dc_buffer_enabled,dc_tracking_mu_eff_min,dc_tracking_mu_eff_max, ... + symbol,dc_tracking_err_buffer,dc_tracking_err_buffer_pos, ... + dc_tracking_err_sum,dc_tracking_abs_err_sum, ... + dc_tracking_valid_count); + end + + if debug_enabled && epoch == epochs + error_power = err * err'; + 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; + obj.debug_struct.update_gradient(1,symbol) = grad.'*grad; + obj.debug_struct.dc_tracking_mu_eff(1,symbol) = dc_tracking_mu_eff; + obj.debug_struct.dc_tracking_est(1,symbol) = obj.e_dc; + + if training + obj.debug_struct.error_tr(1,symbol) = error_power; + obj.debug_struct.update_tr(1,symbol) = update_power; + else + obj.debug_struct.error_dd(1,symbol) = error_power; + obj.debug_struct.update(1,symbol) = update_power; + end + end + end + end + end + end + + methods (Access = private) + function resetTrackingState(obj) + obj.e = zeros(obj.order,1); + obj.e_tr = zeros(obj.order,1); + obj.error = 0; + obj.e_dc = 0; + obj.P = (1/0.05) * eye(obj.order); + obj.debug_struct = struct(); + end + + function state = captureTrackingObjectiveState(obj) + state.e = obj.e; + state.e_tr = obj.e_tr; + state.error = obj.error; + state.e_dc = obj.e_dc; + state.P = obj.P; + state.save_debug = obj.save_debug; + state.debug_struct = obj.debug_struct; + state.dc_tracking_mu = obj.dc_tracking_mu; + end + + function restoreTrackingObjectiveState(obj,state) + obj.e = state.e; + obj.e_tr = state.e_tr; + obj.error = state.error; + obj.e_dc = state.e_dc; + obj.P = state.P; + obj.save_debug = state.save_debug; + obj.debug_struct = state.debug_struct; + obj.dc_tracking_mu = state.dc_tracking_mu; + end + + function N = validSampleLength(obj,N,x,d) + N_available = min(numel(x),numel(d) * obj.sps); + N = min(N,N_available); + N = obj.sps * floor(N / obj.sps); + N = max(0,N); + end + + function adaption_code = adaptionCode(obj) + 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_DCTracking:InvalidAdaptionTechnique", ... + "Unsupported FFE adaption technique."); + end + end + + function initializeTrackingDebug(obj,n_symbols,training) + obj.debug_struct.error = NaN(1,n_symbols); + obj.debug_struct.error_first_epoch = NaN(1,n_symbols); + obj.debug_struct.main_cursor = NaN(1,n_symbols); + obj.debug_struct.mu_nlms = NaN(1,n_symbols); + obj.debug_struct.update_gradient = NaN(1,n_symbols); + obj.debug_struct.dc_tracking_mu_eff = NaN(1,n_symbols); + obj.debug_struct.dc_tracking_est = NaN(1,n_symbols); + + if training + obj.debug_struct.error_tr = NaN(1,n_symbols); + obj.debug_struct.update_tr = NaN(1,n_symbols); + else + obj.debug_struct.error_dd = NaN(1,n_symbols); + obj.debug_struct.update = NaN(1,n_symbols); + end + end + + function [e_dc,mu_eff,err_buffer,err_buffer_pos,err_sum,abs_err_sum,valid_count] = ... + updateDcTracking(obj,err_current,mu_dc,base_mu_eff,persistence_gain, ... + use_persistence,buffer_enabled,mu_eff_min,mu_eff_max,symbol, ... + err_buffer,err_buffer_pos,err_sum,abs_err_sum,valid_count) + + mu_eff = 0; + e_dc = obj.e_dc; + + if buffer_enabled + err_buffer_pos = err_buffer_pos + 1; + if err_buffer_pos > obj.dc_tracking_buffer_len + err_buffer_pos = 1; + end + + old_err = err_buffer(err_buffer_pos); + if isfinite(old_err) + err_sum = err_sum - old_err; + valid_count = valid_count - 1; + if use_persistence + abs_err_sum = abs_err_sum - abs(old_err); + end + end + + if isfinite(err_current) + err_buffer(err_buffer_pos) = err_current; + err_sum = err_sum + err_current; + valid_count = valid_count + 1; + if use_persistence + abs_err_sum = abs_err_sum + abs(err_current); + end + else + err_buffer(err_buffer_pos) = NaN; + end + + if mod(symbol,obj.dc_tracking_buffer_len) == 0 + if valid_count == 0 + err_mean = 0; + err_abs_mean = 0; + else + err_mean = err_sum / valid_count; + err_abs_mean = abs_err_sum / valid_count; + end + + mu_eff = obj.effectiveDcMu(mu_dc,base_mu_eff,persistence_gain, ... + use_persistence,err_mean,err_abs_mean,mu_eff_min,mu_eff_max); + e_dc = e_dc + mu_eff * err_mean; + end + else + if isfinite(err_current) + err_mean = err_current; + err_abs_mean = abs(err_current); + else + err_mean = 0; + err_abs_mean = 0; + end + + mu_eff = obj.effectiveDcMu(mu_dc,base_mu_eff,persistence_gain, ... + use_persistence,err_mean,err_abs_mean,mu_eff_min,mu_eff_max); + e_dc = e_dc + mu_eff * err_mean; + end + end + + function mu_eff = effectiveDcMu(~,mu_dc,base_mu_eff,persistence_gain, ... + use_persistence,err_mean,err_abs_mean,mu_eff_min,mu_eff_max) + + if use_persistence + persistence_scale = abs(err_mean) / (err_abs_mean + eps); + persistence_scale = min(max(persistence_scale,0),1); + mu_eff = mu_dc * (1 + persistence_gain * persistence_scale); + mu_eff = min(max(mu_eff,mu_eff_min),mu_eff_max); + else + mu_eff = base_mu_eff; + end + end + end +end diff --git a/Functions/EQ_recipes/mpi_recipe_dev.m b/Functions/EQ_recipes/mpi_recipe_dev.m index bf50562..8f4022a 100644 --- a/Functions/EQ_recipes/mpi_recipe_dev.m +++ b/Functions/EQ_recipes/mpi_recipe_dev.m @@ -1,8 +1,5 @@ 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. +%mpi_recipe_dev Minimal MPI equalizer comparison recipe. arguments Scpe_sig_raw @@ -16,457 +13,205 @@ arguments options.debug_plots (1,1) logical = false end +%% Execution toggles + +run_conventional_ffe = 1; +run_a2_tracked_levels = 1; +run_a2_residual = 1; +run_a1 = 1; +run_tracking_adaptive = 1; +plot_output_signals = 0; + +%% Shared fixed EQ settings +eq_sps = 2; +eq_order = 50; +eq_adaption = "nlms"; +eq_len_tr = 4096; +eq_epochs_tr = 5; +eq_mu_tr = 0.04; +eq_epochs_dd = 3; +eq_save_debug = false; + +block_update = 1; +if isfield(options.userParameters,"block_update") + block_update = options.userParameters.block_update; +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", 50, ... - "decide", 0, ... - "adaption_technique", "nlms"}; +%% +if plot_output_signals + showLevelScatter(Scpe_sig, Symbols, ... + "fsym", options.fsym, ... + "fignum", 400, ... + "normalize", true); +end +%% Conventional FFE +if run_conventional_ffe + eq_ffe = FFE_plain( ... + "sps", eq_sps, ... + "order", eq_order, ... + "decide", false, ... + "adaption_technique", eq_adaption, ... + "len_tr", eq_len_tr, ... + "epochs_tr", eq_epochs_tr, ... + "mu_tr", eq_mu_tr, ... + "dd_mode", true, ... + "epochs_dd", eq_epochs_dd, ... + "mu_dd", 0.0012, ... + "optmize_mus", false, ... + "plot_mu_optimization", options.debug_plots, ... + "save_debug", eq_save_debug); -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)); + [ffe_results,equalized_signal] = runFfe(eq_ffe, "Conventional FFE", ... + Scpe_sig, Symbols, Tx_bits, options); output.conventional_ffe = ffe_results; + if plot_output_signals + plotEqSignals(equalized_signal,Symbols,options,400,-1); + end end -if 1 - - 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); +%% A2 tracked-level decision +if run_a2_tracked_levels + eq_ffe = FFE_A2TrackedLevels( ... + "sps", eq_sps, ... + "order", eq_order, ... + "decide", true, ... + "adaption_technique", eq_adaption, ... + "len_tr", eq_len_tr, ... + "epochs_tr", eq_epochs_tr, ... + "mu_tr", eq_mu_tr, ... + "dd_mode", false, ... + "epochs_dd", eq_epochs_dd, ... + "mu_dd", 0.012, ... + "dc_smoothing_a2", 1, ... + "dc_level_avg_bufferlength_a2", 112, ... + "dc_level_update_blocklength_a2", block_update, ... + "dc_level_weights_a2", [1], ... + "save_debug", eq_save_debug); + [ffe_results,equalized_signal] = runFfe(eq_ffe, "A2 tracked levels", ... + Scpe_sig, Symbols, Tx_bits, options); + output.a2_adaptive_levels = ffe_results; - 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"); + if plot_output_signals + plotEqSignals(equalized_signal,Symbols,options,410,-1); + end +end - 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); +%% A2 residual correction +if run_a2_residual + eq_ffe = FFE_A2Residual( ... + "sps", eq_sps, ... + "order", eq_order, ... + "decide", false, ... + "adaption_technique", eq_adaption, ... + "len_tr", eq_len_tr, ... + "epochs_tr", eq_epochs_tr, ... + "mu_tr", eq_mu_tr, ... + "dd_mode", true, ... + "epochs_dd", eq_epochs_dd, ... + "mu_dd", 0.012, ... + "dc_smoothing_a2", 1, ... + "dc_level_avg_bufferlength_a2", 112, ... + "dc_level_update_blocklength_a2", block_update, ... + "dc_level_weights_a2", [0.6], ... + "save_debug", eq_save_debug, "optmize_mus",0); - 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"); + [ffe_results,equalized_signal] = runFfe(eq_ffe, "A2 residual", ... + Scpe_sig, Symbols, Tx_bits, options); + output.a2_residual = ffe_results; + + if plot_output_signals + plotEqSignals(equalized_signal,Symbols,options,420,-1); + end +end + +%% A1 moving-average input suppression +if run_a1 + eq_ffe = FFE_A1( ... + "sps", eq_sps, ... + "order", eq_order, ... + "decide", false, ... + "adaption_technique", eq_adaption, ... + "len_tr", eq_len_tr, ... + "epochs_tr", eq_epochs_tr, ... + "mu_tr", eq_mu_tr, ... + "dd_mode", false, ... + "epochs_dd", eq_epochs_dd, ... + "mu_dd", 0.012, ... + "dc_smoothing_a1", 1, ... + "dc_avg_bufferlength_a1", 2048, ... + "dc_avg_update_blocklength_a1", block_update, ... + "save_debug", eq_save_debug); + + [ffe_results,equalized_signal] = runFfe(eq_ffe, "A1", ... + Scpe_sig, Symbols, Tx_bits, options); + output.a1_ff_dc_avg = ffe_results; + + if plot_output_signals + plotEqSignals(equalized_signal,Symbols,options,430,-1); + end +end + +%% DC tracking, adaptive/persistence path +if run_tracking_adaptive + eq_ffe = FFE_DCTracking( ... + "sps", eq_sps, ... + "order", eq_order, ... + "decide", false, ... + "adaption_technique", eq_adaption, ... + "len_tr", eq_len_tr, ... + "epochs_tr", eq_epochs_tr, ... + "mu_tr", eq_mu_tr, ... + "dd_mode", true, ... + "epochs_dd", eq_epochs_dd, ... + "mu_dd", 0.012, ... + "dc_tracking_mu", 0.002, ... + "dc_tracking_adaptive_enabled", true, ... + "dc_tracking_persistence_gain", 0, ... + "dc_tracking_buffer_len", block_update, ... + "optmize_mus", false, ... + "plot_mu_optimization", options.debug_plots, ... + "save_debug", eq_save_debug); + + [ffe_results,equalized_signal] = runFfe(eq_ffe, "DC tracking adaptive", ... + Scpe_sig, Symbols, Tx_bits, options); + output.dc_tracking = ffe_results; + + if plot_output_signals + plotEqSignals(equalized_signal,Symbols,options,450,-1); + end +end end - -%% A2 -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", 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 - +function [ffe_results,equalized_signal] = runFfe(eq_ffe,description,Scpe_sig,Symbols,Tx_bits,options) [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, ... + "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); + ffe_results.metrics.print("description",resultDescription(description,options)); 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; - +function description = resultDescription(prefix,options) + sir = options.dataTable.sir; + if numel(sir) > 1 + sir = sir(1); + end + description = sprintf('%s; SIR %g dB',prefix,sir); end -%% Tracking Adaptive +function plotEqSignals(equalized_signal,Symbols,options,fignum,output_scale) -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, ... + showLevelScatter(equalized_signal .* output_scale, Symbols, ... "fsym", options.fsym, ... - "fignum", 400, ... + "fignum", fignum + 1, ... "normalize", true); - - [~, avg_for_lvl] = showLevelScatter(equalized_signal, Symbols, ... - "fsym", options.fsym, ... - "fignum", 401, ... - "normalize", true); - - -end - end diff --git a/Functions/EQ_visuals/showLevelScatter.m b/Functions/EQ_visuals/showLevelScatter.m index c0bbc3d..22ddfd1 100644 --- a/Functions/EQ_visuals/showLevelScatter.m +++ b/Functions/EQ_visuals/showLevelScatter.m @@ -290,7 +290,7 @@ for levelIdx = 1:numLevels xAxisUs, avg_for_lvl(levelIdx, :), ... options.avgLineMaxPoints, options.avgLineSmoothWindow); plot(ax, xReduced, yReduced, ... - "LineWidth", 2, ... + "LineWidth", 1, ... "Color", cols(levelIdx, :)); end diff --git a/Tests/04_DSP/Equalizer/FFE_plain_test.m b/Tests/04_DSP/Equalizer/FFE_plain_test.m new file mode 100644 index 0000000..f4d513e --- /dev/null +++ b/Tests/04_DSP/Equalizer/FFE_plain_test.m @@ -0,0 +1,105 @@ +classdef FFE_plain_test < IMDDTestCase + methods (Test, TestTags = {'unit', 'fast', 'dsp', 'ffe'}) + function constructorKeepsPlainDefaults(testCase) + ffe = FFE_plain(); + + testCase.verifyEqual(ffe.sps, 2); + testCase.verifyEqual(ffe.order, 15); + testCase.verifyEqual(ffe.len_tr, 4096); + testCase.verifyEqual(ffe.mu_tr, 0); + testCase.verifyEqual(ffe.epochs_tr, 5); + testCase.verifyEqual(ffe.adaption_technique, adaption_method.lms); + testCase.verifyEqual(ffe.dd_mode, 1); + testCase.verifyEqual(ffe.mu_dd, 1e-5); + testCase.verifyEqual(ffe.epochs_dd, 5); + testCase.verifyEqual(ffe.dd_len_fraction, 1); + testCase.verifyFalse(ffe.decide); + testCase.verifyEqual(ffe.e, zeros(ffe.order,1)); + testCase.verifyFalse(isprop(ffe,"dc_tracking_mu")); + testCase.verifyFalse(isprop(ffe,"dc_avg_bufferlength_a1")); + testCase.verifyFalse(isprop(ffe,"dc_level_avg_bufferlength_a2")); + end + + function plainMatchesFullCoreForLms(testCase) + result = runCoreParity(adaption_method.lms,0.02,0.01); + + verifyCoreParity(testCase,result); + end + + function plainMatchesFullCoreForNlms(testCase) + result = runCoreParity(adaption_method.nlms,0.08,0.03); + + verifyCoreParity(testCase,result); + end + + function plainMatchesFullCoreForRls(testCase) + result = runCoreParity(adaption_method.rls,0.995,0.998); + + verifyCoreParity(testCase,result); + end + + function optimizationSignalsUsesShortenedAlignedPrefix(testCase) + ffe = FFE_plain( ... + "sps", 2, ... + "len_tr", 16, ... + "mu_optimization_len", 22); + x = (1:100).'; + d = (1:50).'; + + [x_opt,d_opt,N_opt] = ffe.optimizationSignals(x,d); + + testCase.verifyEqual(N_opt, 22); + testCase.verifyEqual(x_opt, x(1:22)); + testCase.verifyEqual(d_opt, d(1:11)); + end + end +end + +function result = runCoreParity(adaption,mu_tr,mu_dd) + x = [ ... + -2.2; -0.9; 0.8; 2.1; ... + -2.0; -1.1; 1.1; 2.2; ... + -2.1; -0.8; 0.9; 2.0; ... + -2.3; -1.0; 1.0; 2.3]; + d = [ ... + -3; -1; 1; 3; ... + -3; -1; 1; 3; ... + -3; -1; 1; 3; ... + -3; -1; 1; 3]; + settings = { ... + "sps", 1, ... + "order", 3, ... + "len_tr", numel(x), ... + "epochs_tr", 2, ... + "mu_tr", mu_tr, ... + "dd_mode", 1, ... + "epochs_dd", 2, ... + "mu_dd", mu_dd, ... + "adaption_technique", adaption}; + + plain = FFE_plain(settings{:}); + full = FFE(settings{:}); + plain.constellation = unique(d); + full.constellation = unique(d); + plain.P = (1/0.05) * eye(plain.order); + full.P = (1/0.05) * eye(full.order); + + plain.equalize(x,d,plain.mu_tr,plain.epochs_tr,plain.len_tr,true,false); + full.equalize(x,d,full.mu_tr,full.epochs_tr,full.len_tr,true,false); + plain.e_tr = plain.e; + full.e_tr = full.e; + + [result.y_plain,result.dhat_plain] = plain.equalize(x,d,plain.mu_dd,plain.epochs_dd,numel(x),false,false); + [result.y_full,result.dhat_full] = full.equalize(x,d,full.mu_dd,full.epochs_dd,numel(x),false,false); + result.e_plain = plain.e; + result.e_full = full.e; + result.e_tr_plain = plain.e_tr; + result.e_tr_full = full.e_tr; +end + +function verifyCoreParity(testCase,result) + testCase.verifyEqual(result.y_plain,result.y_full,"AbsTol",1e-12); + testCase.verifyEqual(result.dhat_plain,result.dhat_full,"AbsTol",1e-12); + testCase.verifyEqual(result.e_plain,result.e_full,"AbsTol",1e-12); + testCase.verifyEqual(result.e_tr_plain,result.e_tr_full,"AbsTol",1e-12); +end diff --git a/Tests/04_DSP/Equalizer/FFE_test.m b/Tests/04_DSP/Equalizer/FFE_test.m index 7bdfc6f..fa3ba8e 100644 --- a/Tests/04_DSP/Equalizer/FFE_test.m +++ b/Tests/04_DSP/Equalizer/FFE_test.m @@ -22,6 +22,7 @@ classdef FFE_test < IMDDTestCase 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_decision_mode_a2, "residual"); testCase.verifyEqual(ffe.dc_level_weights_a2, 0); testCase.verifyEqual(ffe.dc_tracking_buffer_len, 1); testCase.verifyEqual(ffe.ffe_update_buffer_len, 1); @@ -174,6 +175,64 @@ classdef FFE_test < IMDDTestCase testCase.verifyEqual(ffe.e_dc, 0.36, "AbsTol", 1e-12); end + function ffeDcTrackingClassMatchesFullUnbufferedPath(testCase) + x = zeros(4, 1); + d = ones(4, 1); + settings = { ... + "sps", 1, ... + "order", 1, ... + "dc_tracking_mu", 0.1, ... + "dc_tracking_buffer_len", 1, ... + "save_debug", true, ... + "dd_mode", 0, ... + "adaption_technique", adaption_method.lms}; + + full = FFE(settings{:}); + tracking = FFE_DCTracking(settings{:}); + full.constellation = unique(d); + tracking.constellation = unique(d); + + [y_full,d_hat_full] = full.equalize(x, d, 0, 1, numel(x), true, false); + [y_tracking,d_hat_tracking] = tracking.equalize(x, d, 0, 1, numel(x), true, false); + + testCase.verifyEqual(y_tracking, y_full, "AbsTol", 1e-12); + testCase.verifyEqual(d_hat_tracking, d_hat_full, "AbsTol", 1e-12); + testCase.verifyEqual(tracking.e_dc, full.e_dc, "AbsTol", 1e-12); + testCase.verifyEqual(tracking.debug_struct.dc_tracking_est, ... + full.debug_struct.dc_tracking_est, "AbsTol", 1e-12); + end + + function ffeDcTrackingClassMatchesFullBufferedPersistencePath(testCase) + x = zeros(4, 1); + d = ones(4, 1); + settings = { ... + "sps", 1, ... + "order", 1, ... + "dc_tracking_mu", 0.1, ... + "dc_tracking_adaptive_enabled", true, ... + "dc_tracking_persistence_gain", 1, ... + "dc_tracking_buffer_len", 2, ... + "save_debug", true, ... + "dd_mode", 0, ... + "adaption_technique", adaption_method.lms}; + + full = FFE(settings{:}); + tracking = FFE_DCTracking(settings{:}); + full.constellation = unique(d); + tracking.constellation = unique(d); + + [y_full,d_hat_full] = full.equalize(x, d, 0, 1, numel(x), true, false); + [y_tracking,d_hat_tracking] = tracking.equalize(x, d, 0, 1, numel(x), true, false); + + testCase.verifyEqual(y_tracking, y_full, "AbsTol", 1e-12); + testCase.verifyEqual(d_hat_tracking, d_hat_full, "AbsTol", 1e-12); + testCase.verifyEqual(tracking.e_dc, full.e_dc, "AbsTol", 1e-12); + testCase.verifyEqual(tracking.debug_struct.dc_tracking_mu_eff, ... + full.debug_struct.dc_tracking_mu_eff, "AbsTol", 1e-12); + testCase.verifyEqual(tracking.debug_struct.dc_tracking_est, ... + full.debug_struct.dc_tracking_est, "AbsTol", 1e-12); + end + function ffeBufferDelaysTapUpdateUntilBlockBoundary(testCase) x = ones(4, 1); d = ones(4, 1); @@ -349,6 +408,118 @@ classdef FFE_test < IMDDTestCase testCase.verifyEqual(y, [1; 3; 1; 3.5], "AbsTol", 1e-12); end + function dcLevelAvgA2TrackedLevelsUseAdaptiveDecision(testCase) + x = [0.5; 2.5; 0.5; 2.5; 1.3]; + d = [0; 2; 0; 2; 0]; + + ffe = FFE( ... + "sps", 1, ... + "order", 1, ... + "dc_level_avg_bufferlength_a2", 2, ... + "dc_level_update_blocklength_a2", 1, ... + "dc_level_decision_mode_a2", "tracked_levels", ... + "dc_level_weights_a2", 1, ... + "save_debug", true, ... + "dd_mode", 0, ... + "adaption_technique", adaption_method.lms); + ffe.e = 1; + ffe.constellation = [0; 2]; + + [y,d_hat] = ffe.equalize(x, d, 0, 1, numel(x), false, false); + + testCase.verifyEqual(y, x, "AbsTol", 1e-12); + testCase.verifyEqual(d_hat, d, "AbsTol", 1e-12); + testCase.verifyEqual( ... + ffe.debug_struct.dc_level_decision_level, ... + [0, 2, 0.25, 2.25, 0.5], ... + "AbsTol", 1e-12); + end + + function ffeA2ResidualClassSubtractsScalarWeightedResidual(testCase) + x = (1:4).'; + d = zeros(4, 1); + + ffe = FFE_A2Residual( ... + "sps", 1, ... + "order", 1, ... + "dc_level_avg_bufferlength_a2", 2, ... + "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; 2; 1.5; 2.5], "AbsTol", 1e-12); + testCase.verifyEqual(ffe.debug_struct.dc_level_mpi_est, [0, 0, 1.5, 1.5], "AbsTol", 1e-12); + testCase.verifyEqual(ffe.debug_struct.dc_level_weight, [0, 0, 1, 1], "AbsTol", 1e-12); + testCase.verifyEqual(ffe.debug_struct.dc_level_valid_count, [0, 0, 2, 2]); + end + + function ffeA2ResidualClassMatchesFullResidualPath(testCase) + x = [1; 3; 1; 4; 2; 5]; + d = [0; 2; 0; 2; 0; 2]; + settings = { ... + "sps", 1, ... + "order", 1, ... + "dc_level_avg_bufferlength_a2", 2, ... + "dc_level_update_blocklength_a2", 1, ... + "dc_level_weights_a2", [0.25, 0.75], ... + "save_debug", true, ... + "dd_mode", 0, ... + "adaption_technique", adaption_method.lms}; + + full = FFE(settings{:}); + residual = FFE_A2Residual(settings{:}); + full.e = 1; + residual.e = 1; + full.constellation = [0; 2]; + residual.constellation = [0; 2]; + + [y_full,d_hat_full] = full.equalize(x, d, 0, 1, numel(x), true, false); + [y_residual,d_hat_residual] = residual.equalize(x, d, 0, 1, numel(x), true, false); + + testCase.verifyEqual(y_residual, y_full, "AbsTol", 1e-12); + testCase.verifyEqual(d_hat_residual, d_hat_full, "AbsTol", 1e-12); + testCase.verifyEqual(residual.debug_struct.dc_level_mpi_est, ... + full.debug_struct.dc_level_mpi_est, "AbsTol", 1e-12); + testCase.verifyEqual(residual.debug_struct.dc_level_weight, ... + full.debug_struct.dc_level_weight, "AbsTol", 1e-12); + end + + function dcLevelAvgA2RejectsUnknownDecisionMode(testCase) + testCase.verifyError(@() FFE("dc_level_decision_mode_a2", "unknown"), ... + "FFE:InvalidA2DecisionMode"); + end + + function ffeA2TrackedLevelsClassUsesAdaptiveDecision(testCase) + x = [0.5; 2.5; 0.5; 2.5; 1.3]; + d = [0; 2; 0; 2; 0]; + + ffe = FFE_A2TrackedLevels( ... + "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; 2]; + + [y,d_hat] = ffe.equalize(x, d, 0, 1, numel(x), false, false); + + testCase.verifyEqual(y, x, "AbsTol", 1e-12); + testCase.verifyEqual(d_hat, d, "AbsTol", 1e-12); + testCase.verifyEqual( ... + ffe.debug_struct.dc_level_decision_level, ... + [0, 2, 0.25, 2.25, 0.5], ... + "AbsTol", 1e-12); + end + function dcLevelAvgA2WarnsWhenCombinedWithDcTracking(testCase) testCase.verifyWarning(@() FFE( ... "dc_level_avg_bufferlength_a2", 2, ... diff --git a/projects/Diss/MPI_revisit/algorithms/PLOT_ber_vs_sir.m b/projects/Diss/MPI_revisit/algorithms/PLOT_ber_vs_sir.m index d200677..e8a08d9 100644 --- a/projects/Diss/MPI_revisit/algorithms/PLOT_ber_vs_sir.m +++ b/projects/Diss/MPI_revisit/algorithms/PLOT_ber_vs_sir.m @@ -8,7 +8,8 @@ clear; clc; -scriptDir = fileparts(mfilename("fullpath")); +scriptDir = fileparts("C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\Diss\MPI_revisit\algorithms\"); + runWhFile = fullfile(scriptDir, "combined_by_run_id_results.mat"); configWhFile = fullfile(scriptDir, "combined_by_config_results.mat"); @@ -28,9 +29,9 @@ 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]); +pathLengthToPlot = 1000; % use 300 to see duplicate run_ids per config; set 1000 for the previous path +selectedPamLevels = 8; %wh_config_combined.parameter.pam_level.values; +selectedAlgorithms = algorithmStorageNames([1,2,3,5]); useBoundedLines = true; % switch uncertainty bands on/off here usePolyfit = true; % switch fitted dashed trend lines on/off here @@ -141,6 +142,11 @@ for pamIdx = 1:numel(selectedPamLevels) end end + yline(2.2e-4, 'LineWidth', 1, 'LineStyle', '--', 'HandleVisibility', 'off'); + yline(3.8e-3, 'LineWidth', 1, 'LineStyle', '--', 'HandleVisibility', 'off'); + yline(2e-2, 'LineWidth', 1, 'LineStyle', '--', 'HandleVisibility', 'off'); + ylim([9e-5, 0.1]); + title(sprintf("PAM %.0f, path %.0f m", pamLevel, pathLengthToPlot)); xlabel("SIR (dB)"); ylabel("BER"); diff --git a/projects/Diss/MPI_revisit/investigate_mpi_algorithms.m b/projects/Diss/MPI_revisit/investigate_mpi_algorithms.m index 60bc296..7670ffd 100644 --- a/projects/Diss/MPI_revisit/investigate_mpi_algorithms.m +++ b/projects/Diss/MPI_revisit/investigate_mpi_algorithms.m @@ -27,7 +27,7 @@ db = DBHandler("dataBase", [dsp_options.dataBase],... pamformats = [4,6,8]; baudrates = [112e9,96e9,72e9]; -for i = 1:3 +for i = 1 B = baudrates(i); M = pamformats(i); @@ -37,7 +37,7 @@ for i = 1:3 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', '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); @@ -45,13 +45,18 @@ for i = 1:3 [dataTable,~] = db.queryDB(fp, db.getTableFieldNames('Runs')); + [~, sortIdx] = sort(dataTable.sir, 'descend'); dataTable = dataTable(sortIdx, :); + + % dataTable = dataTable(1,:); + 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)]; + % dsp_options.userParameters.block_update = linspace(1,224,22); wh = DataStorage(dsp_options.userParameters); @@ -68,7 +73,7 @@ for i = 1:3 %% === 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, ... + [results, wh] = submitJobs(run_ids, dsp_options, processingMode.parallel, ... "wh", wh, ... "waitbar", true); @@ -78,3 +83,56 @@ for i = 1:3 end +%% Analyze + +storageNames = fieldnames(wh.sto); + +x_base = dataTable.sir(:).'; +% x_base = dsp_options.userParameters.block_update; + +figure(2026); clf; hold on +for storage_idx = 1:numel(storageNames) + storageName = storageNames{storage_idx}; + result = wh.sto.(storageName); + 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 = x_base; + if numel(x) ~= numel(ber) + x = 1:numel(ber); + end + + h = plot(x,ber,"DisplayName",storageName); + for k = 1:numel(x) + scatter(repmat(x(k),maxPackages,1),ber_mat(:,k), ... + "MarkerEdgeColor",h.Color, ... + "HandleVisibility","off"); + end +end + +beautifyBERplot("logscale",true,"setcolors",false,"setmarkers",true); +legend("Location","best","Interpreter","none"); diff --git a/projects/Diss/MPI_revisit/parallelization_analysis/RUN_parallel_analysis.m b/projects/Diss/MPI_revisit/parallelization_analysis/RUN_parallel_analysis.m index f786e4c..c966986 100644 --- a/projects/Diss/MPI_revisit/parallelization_analysis/RUN_parallel_analysis.m +++ b/projects/Diss/MPI_revisit/parallelization_analysis/RUN_parallel_analysis.m @@ -1,7 +1,7 @@ % === DSP settings === dsp_options = struct(); dsp_options.mode = "run_id"; -dsp_options.recipe = @mpi_recipe_dev; +dsp_options.recipe = @mpi_recipe; dsp_options.append_to_db = false; dsp_options.start_occurence = 1; dsp_options.max_occurences = 15; 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 index 368a72e..d79cf2b 100644 --- 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 @@ -27,13 +27,13 @@ db = DBHandler("dataBase", [dsp_options.dataBase],... fp = QueryFilter(); fp.where('Runs','run_id','GREATER_EQUAL',3153); -fp.where('Runs', 'symbolrate', 'EQUALS', 72e9); % 72 96 112 +fp.where('Runs', 'symbolrate', 'EQUALS', 112e9); % 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', 'pam_level', 'EQUALS', 4); % fp.where('Runs', 'v_bias', 'EQUALS', 2.65); [dataTable,~] = db.queryDB(fp, db.getTableFieldNames('Runs')); @@ -47,7 +47,7 @@ dataTable = dataTable(sortIdx, :); 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_ = load("C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\Diss\MPI_revisit\parallelization_analysis\wh_block_update_pam4_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"); @@ -59,6 +59,7 @@ wh = wh_.wh; storageNames = fieldnames(wh.sto); +storageNames = storageNames([1,2,3,5]); ber_by_storage = struct(); ber_mat_by_storage = struct(); @@ -73,7 +74,10 @@ else curve_colors = lines(max(numel(storageNames),1)); end -fec_ber_threshold = 2e-2;%3.8e-3; +useBoundedLines = true; +usePolyfit = true; +polyfitOrderMax = 4; +fec_ber_threshold = 3.8e-3;%3.8e-3; fit_coeff_by_storage = struct(); required_sir_fec_by_storage = struct(); for storage_idx = 1:numel(storageNames) @@ -124,12 +128,19 @@ for block_idx = 1:numel(block_updates) ber_mat(1:numel(ber_all{k}),k) = ber_all{k}; end - % Replace outliers in ber_mat with NaN (operate column-wise). + % Match the BER-over-SIR plot cleanup: remove high BER points and outliers per SIR. 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; + colData(~isfinite(colData) | colData >= 0.1) = NaN; + validColData = colData(isfinite(colData)); + if isempty(validColData) + ber_mat(:,col) = NaN; + continue + end + outliers = isoutlier(validColData); + validColData(outliers) = NaN; + colData(isfinite(colData)) = validColData; + ber_mat(:,col) = colData; end ber = mean(ber_mat,1,"omitnan"); @@ -153,38 +164,50 @@ for block_idx = 1:numel(block_updates) 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, ... + 30, ... '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); + valid = isfinite(x) & isfinite(ber); + if ~any(valid) + continue + end + + % if useBoundedLines && exist("boundedline","file") + % y_lower = max(ber - ber_min,0); + % y_upper = max(ber_max - ber,0); + % y_bounds = [y_lower(:), y_upper(:)]; + % + % [hl, hp] = boundedline(x(valid).',ber(valid).',y_bounds(valid,:), ... + % 'alpha', 'transparency', 0.08, ... + % 'cmap', curve_color, ... + % 'nan', 'fill', ... + % 'orientation', 'vert'); + % set(hl, ... + % 'LineStyle','none', ... + % 'Marker','none', ... + % 'HandleVisibility','off'); + % set(hp, ... + % 'HandleVisibility','off', ... + % 'LineStyle','none'); + % end + + scatter(x(valid),ber(valid), ... + 20, ... + 'Marker','o', ... + 'MarkerEdgeColor',curve_color, ... + 'MarkerFaceColor',curve_color, ... + 'LineWidth',1, ... + 'DisplayName',storageName); + + fit_mask = valid & ber > 0; + if usePolyfit && nnz(fit_mask) >= 2 + fit_order = min(polyfitOrderMax,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); @@ -210,9 +233,8 @@ for block_idx = 1:numel(block_updates) required_sir_fec_by_storage.(storageName)(block_idx) = sir_req; plot(x_fit,y_fit, ... - 'LineWidth',1.2, ... + 'LineWidth',1.1, ... 'LineStyle','--', ... - 'Marker','none', ... 'Color',curve_color, ... 'HandleVisibility','off'); end @@ -221,10 +243,17 @@ for block_idx = 1:numel(block_updates) 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"); + yline(2.2e-4, 'LineWidth', 1, 'LineStyle', '--', 'HandleVisibility', 'off'); + yline(3.8e-3, 'LineWidth', 1, 'LineStyle', '--', 'HandleVisibility', 'off'); + yline(2e-2, 'LineWidth', 1, 'LineStyle', '--', 'HandleVisibility', 'off'); + ylim([9e-5, 0.1]); + + set(gca, "YScale", "log"); + xlim([15,45]); + grid on; + box on; + legend("Location", "best", "Interpreter", "none"); end %% Required SIR to reach FEC over parallelization diff --git a/workerError.mat b/workerError.mat index afabc042e105ce166ec5de58dc7dd7ce9e6a6345..b33c61c71ed4dd97732ddc29026a029402df862b 100644 GIT binary patch delta 1000 zcmVQFgQ9iG9WTAGB%M>BavVQk#uQ)v;zPD zc%1E6PjAyO6n8^c)@c+tF~Oy`O+uRZBP6s_x3XzNR1|d;6N;3%iJKZ4CyJA5aN@|3 zBgY+g#Xf0qe=&Mi5Mvtk_pQWFP4e?lj4uqY#b_+o@3FRr z`{<(&JrDVYA?$%vg-nqL4SJxEMs#TxW~?a*Nc6x3bc zM`nOs0vfv~Ci0?N6NVIlR4J!-k4Yg#m;08%`S>-V-;tgb#Z?pBCC^LJ^SfY_@5FB@ z2R=^%RFSTyAUn?q7CFzarSYYd@p7DH(J!RN?Vm~ZMa5pm%jJRe-u(%7eEgbX$Mwec z;k)gftu|p(0~2KJg(OOUJ@cbS=Z3Q>i>sxx8BA@Cl+=o{ymTu@m+ok4R5=qWbO!gN2MeM0Ll*;@+^(VAB!t)u!NosVI23KJIGXd38|B z!^@$9<8Urs-nt#5>ANN?J($?KMLok03AQ3ya`94D@Nj+TM)c3c%lExpytEV?twr*} zg4uzA*^6m87f%|VWd#q{6E}15@#j-6J{G%v6h3I!ad5yoY3}v2qTu2B$PGy4^-*$(|4H9}X3xi~OU=hk+xKSAW9?FZ^Keu5+u7}(&y|nceSWrm z{<`YeRol(+f4>m>Y59Mm;O2UE{#S;6dy?O@a}F6-o_bvZ<2bk;|NNilU*#bYnTO=P zXA^KHXE9kg)6NUQn4jXM+OJN*tHvMIBkgMA-Owi6tKaL%cWVcTbI74&9!N4t&Lm?Q zKVBR=r&JCBavVQk#uQ)PXqt} zc%1E7PjA~c6sHw8jnm9+uN%<%81xWi#Re3?usU&KEJeLKaIy`B0F9Pui;ezOq=L4m z9CzGtA0fwngnWb^d)!CZamSJpWz(!|DN_C{O9>Ftqxc?w@AvrU39JJEK0X9skv*4K zt%11Q5^{yrlZWw>7fr~ZvfJT*a-H{+^tBZIq{i*l=vzhf>9k(-5<904CU)+N{+02w z8ujJ!SFEt)|4C)G-)}0@{WZ~6WzS_+|77*OsDHEC zW#4`m)nb*t^%1Zj*+R%|!Ej_zfCtE>c!*H|97G5lAuvEg7}}IjXdVK8gI#=h5yw2~kI%Mk=xX?ks#x4@; zH3R?83VtrPHuM960;2b@4T-+(nei&oW5RUo31E1^8%(ZnU)$eP_?;+n#wy^x+V3D#m`SXW&?v!m$9De#w33$H@ObGR8en{68q` zKRafy&nAcCWASc(4A4=;MV@(4SMYFo+Em&7_?KL~>?=6-OY!o|>l@v8W@C4x4~TCB zp^FFc6v@L&O~J$Eu^G``ZoRZ~@zPUp^v;tP=FAp9n#9V%PUwv z3LhlwJD9RYl6(CuD|om(F#}S0J=^-}M7;c!fn%uP7#84vg`EYjSklURV!}FgciAp* zgbAiPc8?IDcmydj`0fDhkr!ZUIV3rMvht;`6?|NtZS@>@JKc+TTdmK|e`<<8&wrA5 z^1Lg`f1(m4K~ae_Azn<+4jq@jvE-j3&McVVKF;9dzmgYRR-ROR{y%GGo##oB{kU#% zxUBrFcpM&oU1%Ks?|ETy{U1v4{v{tx&qs^vw=OguH#H9~9*^%YG#)oC|1KWK<~7=< zu20@yT>slr@%Zw;QS5)?s^WFodd2-bAN{QOf1}{$^6dC8kL|8Xb}Ncql@FCa)ep6Q zP0~A=Gfr|BB99~Um3ChTrtDNNRDUG{uj+rapJ+h;IUtj(2Px8(_5ch!NGV?-^ zmE=w^F4Nb0Q}=|*MWKWf6rXQ3ZN6Wcc)N=I;jDhmi?`bUtn)jMzMOxg-KYNsWL{%g WzEA&Y&iHBbZPtA{FP{LYU^Ve#j8we<