From aa210e5352adca49a324932e107ded71ca7a2f68 Mon Sep 17 00:00:00 2001 From: Silas Oettinghaus Date: Thu, 9 Jul 2026 14:23:16 +0200 Subject: [PATCH] more on MPI tracking --- Classes/04_DSP/Equalizer/FFE.m | 304 +++++++++++++++++++++++--- Functions/EQ_blocks/ffe.m | 3 +- Functions/EQ_recipes/mpi_recipe_dev.m | 8 +- Tests/04_DSP/Equalizer/FFE_test.m | 2 + 4 files changed, 277 insertions(+), 40 deletions(-) diff --git a/Classes/04_DSP/Equalizer/FFE.m b/Classes/04_DSP/Equalizer/FFE.m index b73595d..df8b347 100644 --- a/Classes/04_DSP/Equalizer/FFE.m +++ b/Classes/04_DSP/Equalizer/FFE.m @@ -51,6 +51,23 @@ classdef FFE < handle mu_optimization_len plot_mu_optimization = 0; mu_optimization_fignum = 3010; + optimize_dc_params = 0; + dc_optimization + dc_optimization_iter = 0; + dc_optimization_len + dc_optimization_max_evals + dc_optimization_delay_weight + dc_optimization_smoothing_len + end + + properties (Access = private) + dc_alpha = 0.98 + dc_gamma = 1e-6 + dc_mu_min = 1e-6 + dc_mu_max = 3e-1 + dc_mu_eff_min = 0 + dc_mu_eff_max = inf + dc_power_exponent = 1 end methods @@ -79,6 +96,11 @@ classdef FFE < handle options.mu_optimization_len = 2^15; options.plot_mu_optimization = 0; options.mu_optimization_fignum = 3010; + options.optimize_dc_params = 0; + options.dc_optimization_len = 2^15; + options.dc_optimization_max_evals = 20; + options.dc_optimization_delay_weight = 1e-3; + options.dc_optimization_smoothing_len = 501; end @@ -113,6 +135,13 @@ classdef FFE < handle obj.P = (1/delta) * eye(obj.order); end + if obj.optimize_dc_params + obj.optimizeDcParams(X.signal,D.signal); + obj.e = zeros(obj.order,1); + obj.e_dc = 0; + obj.P = (1/delta) * eye(obj.order); + end + % Training Mode training = 1; showviz = 0; @@ -182,11 +211,14 @@ classdef FFE < handle end P_err = 0; - alpha = 0.98; err_prev = 0; - gamma_dc = 1e-6; - mu_min = 1e-6; - mu_max = 3e-1; + dc_alpha = obj.dc_alpha; + dc_gamma = obj.dc_gamma; + dc_mu_min = obj.dc_mu_min; + dc_mu_max = obj.dc_mu_max; + dc_mu_eff_min = obj.dc_mu_eff_min; + dc_mu_eff_max = obj.dc_mu_eff_max; + dc_power_exponent = obj.dc_power_exponent; debug_enabled = obj.save_debug; if debug_enabled @@ -238,6 +270,16 @@ classdef FFE < handle if training || obj.dd_mode switch obj.adaption_technique + + case adaption_method.nlms + + % mu used as update weight (suggestion: 0.01-0.05; bit higher during tr) + normU = ((U.'*U)) + eps; + weight = mu / normU; + grad = err(symbol) * U; + update = grad * weight; + obj.e = obj.e + update; + case adaption_method.lms @@ -247,14 +289,7 @@ classdef FFE < handle update = grad * weight; obj.e = obj.e + update; - case adaption_method.nlms - - % mu used as update weight (suggestion: 0.01-0.05; bit higher during tr) - normU = ((U.'*U)) + eps; - weight = mu / normU; - grad = err(symbol) * U; - update = grad * weight; - obj.e = obj.e + update; + case adaption_method.rls @@ -275,23 +310,24 @@ classdef FFE < handle if obj.mu_dc ~= 0 if obj.adaptive_dc_enabled - delta_mu = gamma_dc * err(symbol) * err_prev * (U.'*U); - obj.mu_dc = min(max(obj.mu_dc + delta_mu,mu_min),mu_max); + delta_mu = dc_gamma * err(symbol) * err_prev * (U.'*U); + obj.mu_dc = min(max(obj.mu_dc + delta_mu,dc_mu_min),dc_mu_max); err_prev = err(symbol); - P_err = alpha*P_err + (1-alpha)*err(symbol)^2; - mu_dc_eff = obj.mu_dc / (P_err + eps); + P_err = dc_alpha*P_err + (1-dc_alpha)*err(symbol)^2; + mu_dc_eff = obj.mu_dc / ((P_err + eps)^dc_power_exponent); else mu_dc_eff = obj.mu_dc; end + mu_dc_eff = min(max(mu_dc_eff,dc_mu_eff_min),dc_mu_eff_max); obj.e_dc = obj.e_dc + mu_dc_eff * err(symbol); end end - - if debug_enabled && epoch == 1 - obj.debug_struct.error_first_epoch(1,symbol) = err(symbol) * err(symbol)'; - 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)'; @@ -373,13 +409,54 @@ classdef FFE < handle end end - function [x_opt,d_opt,N_opt] = optimizationSignals(obj,x,d) + function optimizeDcParams(obj,x,d) + [x_opt,d_opt,N_opt] = obj.optimizationSignals(x,d,obj.dc_optimization_len); + + vars = optimizableVariable("mu_dc",[1e-5,1e-1],"Transform","log"); + if obj.adaptive_dc_enabled + vars = [vars, ... + optimizableVariable("dc_alpha",[0.85,0.995]), ... + optimizableVariable("dc_gamma",[1e-7,3e-5],"Transform","log"), ... + optimizableVariable("dc_power_exponent",[0,1]), ... + optimizableVariable("dc_mu_eff_max",[1e-3,3e-1],"Transform","log")]; + end + + obj.dc_optimization_iter = 0; + fprintf("FFE DC opt uses fixed mu_tr=%9.3e, mu_dd=%9.3e on %d samples / %d symbols\n", ... + obj.mu_tr,obj.mu_dd,N_opt,numel(d_opt)); + obj.dc_optimization = bayesopt(@(p)obj.dcObjective(p,x_opt,d_opt),vars, ... + "MaxObjectiveEvaluations",obj.dc_optimization_max_evals, ... + "AcquisitionFunctionName","expected-improvement-plus", ... + "IsObjectiveDeterministic",false, ... + "Verbose",0, ... + "PlotFcn",[]); + + best = obj.dc_optimization.XAtMinObjective; + obj.mu_dc = best.mu_dc; + if obj.adaptive_dc_enabled + obj.dc_alpha = best.dc_alpha; + obj.dc_gamma = best.dc_gamma; + obj.dc_power_exponent = best.dc_power_exponent; + obj.dc_mu_eff_max = best.dc_mu_eff_max; + fprintf("\nFFE DC opt done: mu_dc=%9.3e, alpha=%6.3f, gamma=%9.3e, p=%5.2f, mu_eff_max=%9.3e, objective=%9.3e\n", ... + obj.mu_dc,obj.dc_alpha,obj.dc_gamma,obj.dc_power_exponent,obj.dc_mu_eff_max,obj.dc_optimization.MinObjective); + else + fprintf("\nFFE DC opt done: mu_dc=%9.3e, objective=%9.3e\n", ... + obj.mu_dc,obj.dc_optimization.MinObjective); + end + 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(obj.mu_optimization_len) || obj.mu_optimization_len <= 0 || isinf(obj.mu_optimization_len) + if isempty(opt_len) || opt_len <= 0 || isinf(opt_len) N_opt = N_available; else - N_opt = min(N_available,max(obj.len_tr,obj.mu_optimization_len)); + N_opt = min(N_available,max(obj.len_tr,opt_len)); end N_opt = obj.sps * floor(N_opt / obj.sps); @@ -408,17 +485,7 @@ classdef FFE < handle [signal,~] = obj.equalize(x,d,params.mu_dd,obj.epochs_dd,numel(x),0,0); end - 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",10, ... - "skip_end",10, ... - "returnErrorLocation",1); + [ber,errors] = obj.berObjective(signal,d); objective = ber; if ~isfinite(objective) objective = inf; @@ -441,6 +508,175 @@ classdef FFE < handle obj.mu_dc = old_mu_dc; end + function objective = dcObjective(obj,params,x,d) + state = obj.captureObjectiveState(); + cleanup = onCleanup(@()obj.restoreObjectiveState(state)); + + obj.applyDcObjectiveParams(params); + obj.save_debug = 1; + obj.e = zeros(obj.order,1); + obj.e_dc = 0; + obj.P = (1/0.05) * eye(obj.order); + obj.debug_struct = struct(); + + N_tr = min(obj.len_tr,numel(x)); + [signal,~] = obj.equalize(x,d,obj.mu_tr,obj.epochs_tr,N_tr,1,0); + if obj.dd_mode + [signal,~] = obj.equalize(x,d,obj.mu_dd,obj.epochs_dd,numel(x),0,0); + end + + [ber,errors] = obj.berObjective(signal,d); + [delay_symbols,delay_corr] = obj.dcDelayObjective(signal,d); + delay_penalty = obj.dc_optimization_delay_weight * abs(delay_symbols) / max(numel(d),1); + objective = ber + delay_penalty; + if ~isfinite(objective) + objective = inf; + end + + obj.dc_optimization_iter = obj.dc_optimization_iter + 1; + if obj.adaptive_dc_enabled + fprintf("\rFFE DC opt %02d: mu_dc=%9.3e, alpha=%6.3f, gamma=%9.3e, p=%5.2f, mu_eff_max=%9.3e, BER=%9.3e, delay=%7.0f, corr=%6.3f, obj=%9.3e, errors=%d", ... + obj.dc_optimization_iter,params.mu_dc,params.dc_alpha,params.dc_gamma,params.dc_power_exponent,params.dc_mu_eff_max, ... + ber,delay_symbols,delay_corr,objective,errors); + else + fprintf("\rFFE DC opt %02d: mu_dc=%9.3e, BER=%9.3e, delay=%7.0f, corr=%6.3f, obj=%9.3e, errors=%d", ... + obj.dc_optimization_iter,params.mu_dc,ber,delay_symbols,delay_corr,objective,errors); + end + end + + function applyDcObjectiveParams(obj,params) + var_names = string(params.Properties.VariableNames); + if any(var_names == "mu_dc") + obj.mu_dc = params.mu_dc; + end + if any(var_names == "dc_alpha") + obj.dc_alpha = params.dc_alpha; + end + if any(var_names == "dc_gamma") + obj.dc_gamma = params.dc_gamma; + end + if any(var_names == "dc_power_exponent") + obj.dc_power_exponent = params.dc_power_exponent; + end + if any(var_names == "dc_mu_eff_max") + obj.dc_mu_eff_max = params.dc_mu_eff_max; + end + end + + function state = captureObjectiveState(obj) + state.e = obj.e; + state.e_dc = obj.e_dc; + state.P = obj.P; + state.save_debug = obj.save_debug; + state.debug_struct = obj.debug_struct; + state.mu_dc = obj.mu_dc; + state.dc_alpha = obj.dc_alpha; + state.dc_gamma = obj.dc_gamma; + state.dc_power_exponent = obj.dc_power_exponent; + state.dc_mu_eff_max = obj.dc_mu_eff_max; + end + + function restoreObjectiveState(obj,state) + obj.e = state.e; + obj.e_dc = state.e_dc; + obj.P = state.P; + obj.save_debug = state.save_debug; + obj.debug_struct = state.debug_struct; + obj.mu_dc = state.mu_dc; + obj.dc_alpha = state.dc_alpha; + obj.dc_gamma = state.dc_gamma; + obj.dc_power_exponent = state.dc_power_exponent; + obj.dc_mu_eff_max = state.dc_mu_eff_max; + 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",10, ... + "skip_end",10, ... + "returnErrorLocation",1); + end + + function [delay_symbols,delay_corr] = dcDelayObjective(obj,signal,d) + delay_symbols = 0; + delay_corr = 0; + if ~isfield(obj.debug_struct,"e_dc_eff") || isempty(obj.debug_struct.e_dc_eff) + return + end + + smooth_len = obj.dc_optimization_smoothing_len; + e_dc_eff_s = movmean(obj.debug_struct.e_dc_eff(:),smooth_len,"omitnan"); + avg_lvl_dc = obj.averageLevelTrace(signal,d,smooth_len); + + xcorr_len = min(numel(e_dc_eff_s),numel(avg_lvl_dc)); + if xcorr_len < 2 + return + end + + inv_dc_xcorr = -e_dc_eff_s(1:xcorr_len); + avg_lvl_xcorr = avg_lvl_dc(1:xcorr_len); + inv_dc_xcorr = fillmissing(inv_dc_xcorr,"linear","EndValues","nearest"); + avg_lvl_xcorr = fillmissing(avg_lvl_xcorr,"linear","EndValues","nearest"); + inv_dc_xcorr = inv_dc_xcorr - mean(inv_dc_xcorr,"omitnan"); + avg_lvl_xcorr = avg_lvl_xcorr - mean(avg_lvl_xcorr,"omitnan"); + + if rms(inv_dc_xcorr) <= eps || rms(avg_lvl_xcorr) <= eps + return + end + + [dc_level_xcorr,dc_level_lags] = xcorr(inv_dc_xcorr,avg_lvl_xcorr,"coeff"); + [delay_corr,delay_idx] = max(dc_level_xcorr); + delay_symbols = dc_level_lags(delay_idx); + if ~isfinite(delay_corr) + delay_corr = 0; + delay_symbols = 0; + end + end + + function avg_lvl_dc = averageLevelTrace(obj,signal,d,smooth_len) + signal = signal(:); + d = d(:); + n_symbols = min(numel(signal),numel(d)); + signal = signal(1:n_symbols); + d = d(1:n_symbols); + levels = unique(d); + avg_for_lvl = NaN(numel(levels),n_symbols); + + for level_idx = 1:numel(levels) + level_mask = d == levels(level_idx); + level_samples = signal(level_mask); + if isempty(level_samples) + continue + end + + smooth_window = min(smooth_len,numel(level_samples)); + avg_for_lvl(level_idx,level_mask) = movmean(level_samples,smooth_window,"omitnan","Endpoints","shrink"); + avg_for_lvl(level_idx,:) = obj.interpolateMissingAverage(avg_for_lvl(level_idx,:)); + end + + avg_lvl_dc = mean(avg_for_lvl,1,"omitnan").'; + end + + function level_average = interpolateMissingAverage(~,level_average) + valid_samples = isfinite(level_average); + if nnz(valid_samples) == 0 + return + elseif nnz(valid_samples) == 1 + level_average(:) = level_average(valid_samples); + return + end + + t = 1:numel(level_average); + level_average(~valid_samples) = interp1(t(valid_samples),level_average(valid_samples), ... + t(~valid_samples),"linear","extrap"); + end + function plotMuOptimization(obj) if isempty(obj.mu_optimization) return diff --git a/Functions/EQ_blocks/ffe.m b/Functions/EQ_blocks/ffe.m index 6fe7e41..0cbd42b 100644 --- a/Functions/EQ_blocks/ffe.m +++ b/Functions/EQ_blocks/ffe.m @@ -84,6 +84,7 @@ try end try eq_.mu_optimization = []; + eq_.dc_optimization = []; end ffe_results.config = Equalizerstruct(); @@ -196,4 +197,4 @@ function displayAnalysis(eq_noise, eq_signal_sd, rx_signal, eq_, tx_symbols, M, % eq_signal_sd.eye(eq_signal_sd.fs,M,"displayname",'Eye','fignum',105); -end \ No newline at end of file +end diff --git a/Functions/EQ_recipes/mpi_recipe_dev.m b/Functions/EQ_recipes/mpi_recipe_dev.m index 922016c..ea7c5fb 100644 --- a/Functions/EQ_recipes/mpi_recipe_dev.m +++ b/Functions/EQ_recipes/mpi_recipe_dev.m @@ -19,15 +19,15 @@ function output = mpi_recipe_dev(Scpe_sig_raw, Symbols, Tx_bits, options) Scpe_sig = preprocessSignal(Scpe_sig_raw, Symbols, options.fsym, ... "mode", "auto", ... "debug_plots", options.debug_plots); + output = struct(); - eq_settings = { ... "epochs_tr", 5, ... "epochs_dd", 5, ... "len_tr", 4096, ... "mu_tr",0.04, ... "mu_dd",0.012, ... - "mu_dc", 0.012, ... + "mu_dc", 0, ... "adaptive_dc_enabled", 0, ... "order", 25, ... "sps", 2, ... @@ -36,6 +36,7 @@ function output = mpi_recipe_dev(Scpe_sig_raw, Symbols, Tx_bits, options) "dd_mode", 1, ... "adaption_technique", "nlms", ... "plot_mu_optimization", options.debug_plots,... + "optimize_dc_params", false, ... "save_debug",true}; eq_ffe = FFE(eq_settings{:}); @@ -45,9 +46,6 @@ function output = mpi_recipe_dev(Scpe_sig_raw, Symbols, Tx_bits, options) "fignum", 400, ... "normalize", true); - % Scpe_sig_raw.spectrum("normalizeTo0dB",1,"fft_length",4096*4,"fignum",401); - - %options.dataTable.sir; %% NORMAL FFE if 1 diff --git a/Tests/04_DSP/Equalizer/FFE_test.m b/Tests/04_DSP/Equalizer/FFE_test.m index a51e283..8b82541 100644 --- a/Tests/04_DSP/Equalizer/FFE_test.m +++ b/Tests/04_DSP/Equalizer/FFE_test.m @@ -12,6 +12,8 @@ classdef FFE_test < IMDDTestCase testCase.verifyEqual(ffe.dd_mode, 1); testCase.verifyEqual(ffe.mu_dd, 1e-5); testCase.verifyEqual(ffe.epochs_dd, 5); + testCase.verifyFalse(logical(ffe.optimize_dc_params)); + testCase.verifyEqual(ffe.dc_optimization_delay_weight, 1e-3); testCase.verifyFalse(ffe.decide); testCase.verifyEqual(ffe.e, zeros(ffe.order, 1)); testCase.verifyEqual(ffe.error, 0);