diff --git a/Classes/04_DSP/Equalizer/FFE.m b/Classes/04_DSP/Equalizer/FFE.m index f655098..b73595d 100644 --- a/Classes/04_DSP/Equalizer/FFE.m +++ b/Classes/04_DSP/Equalizer/FFE.m @@ -33,6 +33,7 @@ classdef FFE < handle epochs_dd dd_len_fraction mu_dc + adaptive_dc_enabled e_dc P % covariance matrix of rls @@ -48,6 +49,8 @@ classdef FFE < handle mu_optimization mu_optimization_iter = 0; mu_optimization_len + plot_mu_optimization = 0; + mu_optimization_fignum = 3010; end methods @@ -67,12 +70,15 @@ classdef FFE < handle options.epochs_dd = 5; options.dd_len_fraction = 1; options.mu_dc = 0; + options.adaptive_dc_enabled = false; 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 @@ -175,12 +181,40 @@ classdef FFE < handle epochs = 1; end + P_err = 0; + alpha = 0.98; + err_prev = 0; + gamma_dc = 1e-6; + mu_min = 1e-6; + mu_max = 3e-1; + + debug_enabled = obj.save_debug; + if debug_enabled + n_symbols_debug = ceil(N / obj.sps); + obj.debug_struct.error = NaN(1,n_symbols_debug); + obj.debug_struct.error_first_epoch = NaN(1,n_symbols_debug); + obj.debug_struct.main_cursor = NaN(1,n_symbols_debug); + obj.debug_struct.mu_nlms = NaN(1,n_symbols_debug); + obj.debug_struct.update_gradient = NaN(1,n_symbols_debug); + obj.debug_struct.mu_dc_eff = NaN(1,n_symbols_debug); + obj.debug_struct.e_dc_eff = NaN(1,n_symbols_debug); + + if training + obj.debug_struct.error_tr = NaN(1,n_symbols_debug); + obj.debug_struct.update_tr = NaN(1,n_symbols_debug); + else + obj.debug_struct.error_dd = NaN(1,n_symbols_debug); + obj.debug_struct.update = NaN(1,n_symbols_debug); + end + end + for epoch = 1 : epochs symbol = 0; % obj.e_dc = 0; for sample = 1 : obj.sps : N symbol = symbol+1; + mu_dc_eff = 0; U = x(obj.order+sample-1:-1:sample); @@ -240,31 +274,42 @@ classdef FFE < handle end if obj.mu_dc ~= 0 - obj.e_dc = obj.e_dc + obj.mu_dc * err(symbol); + 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); + err_prev = err(symbol); + P_err = alpha*P_err + (1-alpha)*err(symbol)^2; + mu_dc_eff = obj.mu_dc / (P_err + eps); + else + mu_dc_eff = obj.mu_dc; + end + + obj.e_dc = obj.e_dc + mu_dc_eff * err(symbol); end end - if obj.save_debug - obj.debug_struct.error(epoch,symbol) = err(symbol) * err(symbol)'; % Instantaneous square error - - if training - obj.debug_struct.error_tr(epoch,symbol) = err(symbol) * err(symbol)'; % Instantaneous square error - obj.debug_struct.update_tr(epoch,symbol) = update.'*update ./ rms(obj.e); - else - % if symbol == length(d) - % M = numel(unique(d)); - % bits_ref = PAMmapper(M,0).demap(d); - % bits_eq = PAMmapper(M,0).demap(d_hat); - % [~,~,obj.debug_struct.ber(epoch)] = calc_ber(bits_eq,bits_ref); - % end - end + if debug_enabled && epoch == 1 + obj.debug_struct.error_first_epoch(1,symbol) = err(symbol) * err(symbol)'; + end - % obj.debug_struct.main_cursor(epoch,symbol) = abs(obj.e(maincursor_pos)); - % obj.debug_struct.mu_nlms(epoch,symbol) = weight; - % obj.debug_struct.update_gradient(epoch,symbol) = grad.'*grad; - % obj.debug_struct.update(epoch,symbol) = update.'*update ./ rms(obj.e); - + if debug_enabled && epoch == epochs + error_power = err(symbol) * err(symbol)'; + update_power = update.'*update ./ (rms(obj.e) + 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.mu_dc_eff(1,symbol) = mu_dc_eff; + obj.debug_struct.e_dc_eff(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 @@ -272,17 +317,6 @@ classdef FFE < handle end - function N_dd = ddLength(obj,N) - if isempty(obj.dd_len_fraction) || obj.dd_len_fraction <= 0 || obj.dd_len_fraction >= 1 - N_dd = N; - return - end - - N_dd = floor(N * obj.dd_len_fraction); - N_dd = max(obj.sps,N_dd); - N_dd = min(N,N_dd); - end - function optimizeMus(obj,x,d) [x_opt,d_opt,N_opt] = obj.optimizationSignals(x,d); @@ -333,6 +367,10 @@ classdef FFE < handle fprintf("\nFFE 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 [x_opt,d_opt,N_opt] = optimizationSignals(obj,x,d) @@ -402,5 +440,99 @@ classdef FFE < handle obj.save_debug = old_debug; obj.mu_dc = old_mu_dc; 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 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 end end diff --git a/Functions/EQ_blocks/ffe.m b/Functions/EQ_blocks/ffe.m index 964905e..6fe7e41 100644 --- a/Functions/EQ_blocks/ffe.m +++ b/Functions/EQ_blocks/ffe.m @@ -1,4 +1,4 @@ -function [ffe_results] = ffe(eq_, M, rx_signal, tx_symbols, tx_bits, options) +function [ffe_results,eq_signal_sd] = ffe(eq_, M, rx_signal, tx_symbols, tx_bits, options) % FFE Processes signals through FFE equalizer % % Inputs: diff --git a/Functions/EQ_recipes/mpi_recipe_dev.m b/Functions/EQ_recipes/mpi_recipe_dev.m index 7babb83..922016c 100644 --- a/Functions/EQ_recipes/mpi_recipe_dev.m +++ b/Functions/EQ_recipes/mpi_recipe_dev.m @@ -20,22 +20,23 @@ function output = mpi_recipe_dev(Scpe_sig_raw, Symbols, Tx_bits, options) "mode", "auto", ... "debug_plots", options.debug_plots); - - mu_dc = 0;%1e-5; - + eq_settings = { ... "epochs_tr", 5, ... "epochs_dd", 5, ... - "len_tr", 4096*2, ... - "mu_dd",0.02, ... - "mu_tr",0.2, ... + "len_tr", 4096, ... + "mu_tr",0.04, ... + "mu_dd",0.012, ... + "mu_dc", 0.012, ... + "adaptive_dc_enabled", 0, ... "order", 25, ... "sps", 2, ... "decide", 0, ... - "optmize_mus", 1, ... + "optmize_mus", 0, ... "dd_mode", 1, ... "adaption_technique", "nlms", ... - "mu_dc", mu_dc}; + "plot_mu_optimization", options.debug_plots,... + "save_debug",true}; eq_ffe = FFE(eq_settings{:}); @@ -51,7 +52,7 @@ function output = mpi_recipe_dev(Scpe_sig_raw, Symbols, Tx_bits, options) %% NORMAL FFE if 1 - ffe_results = ffe(eq_ffe, options.M, Scpe_sig, Symbols, Tx_bits, ... + [ffe_results,equalized_signal] = ffe(eq_ffe, options.M, Scpe_sig, Symbols, Tx_bits, ... "precode_mode", options.duob_mode, ... 'showAnalysis', options.debug_plots, ... "postFFE", [], ... @@ -60,8 +61,84 @@ function output = mpi_recipe_dev(Scpe_sig_raw, Symbols, Tx_bits, options) ffe_results.metrics.print("description",sprintf('Normal FFE; SIR %d dB',options.dataTable.sir)); output.ffe_package = ffe_results; + [~, avg_for_lvl] = showLevelScatter(equalized_signal, Symbols, ... + "fsym", options.fsym, ... + "fignum", 401, ... + "normalize", true); + + output.ffe_debug = struct(); + dbg = eq_ffe.debug_struct; + + + + smooth_len = 501; + output.ffe_debug.smoothing_window = smooth_len; + error_first_s = movmean(dbg.error_first_epoch(:),smooth_len,"omitnan"); + error_last_s = movmean(dbg.error(:),smooth_len,"omitnan"); + mu_dc_eff_s = movmean(dbg.mu_dc_eff(:),smooth_len,"omitnan"); + e_dc_eff_s = movmean(dbg.e_dc_eff(:),smooth_len,"omitnan"); + + figure(460); clf; + tiledlayout(2,1,"TileSpacing","compact","Padding","compact"); + + nexttile; + semilogy(error_first_s + eps,"DisplayName","first epoch"); + hold on; + semilogy(error_last_s + eps,"DisplayName","last epoch"); + grid on; + ylabel("error"); + title("FFE debug error"); + legend("Location","best"); + + nexttile;hold on + % plot(eq_sig_mov,"DisplayName","err dc eff"); + plot(mu_dc_eff_s,"DisplayName","mu DC - either fixed or adaptive"); + e_dc_scale = max(abs(e_dc_eff_s),[],"omitnan") + eps; + plot(-1.*e_dc_eff_s./e_dc_scale,"DisplayName","Inverted DC-tracking; Value that is subtracted during EQ"); + avg_lvl_dc=mean(avg_for_lvl,1,"omitnan"); + plot(avg_lvl_dc./0.01,"DisplayName","Smoothed EQ output signal; calc'd by showLevelScatter"); + grid on; + xlabel("Symbol"); + ylabel("mu dc eff"); + title("Effective adaptive DC step"); + legend + + inv_dc_track = -e_dc_eff_s(:); + avg_lvl_dc = avg_lvl_dc(:); + xcorr_len = min(numel(inv_dc_track),numel(avg_lvl_dc)); + inv_dc_xcorr = inv_dc_track(1:xcorr_len); + avg_lvl_xcorr = avg_lvl_dc(1:xcorr_len); + inv_dc_xcorr = fillmissing(inv_dc_xcorr,"linear","EndValues","nearest"); + avg_lvl_xcorr = fillmissing(avg_lvl_xcorr,"linear","EndValues","nearest"); + inv_dc_xcorr = inv_dc_xcorr - mean(inv_dc_xcorr,"omitnan"); + avg_lvl_xcorr = avg_lvl_xcorr - mean(avg_lvl_xcorr,"omitnan"); + [dc_level_xcorr,dc_level_lags] = xcorr(inv_dc_xcorr,avg_lvl_xcorr,"coeff"); + [dc_level_corr,dc_level_idx] = max(dc_level_xcorr); + dc_level_delay_symbols = dc_level_lags(dc_level_idx); + output.ffe_debug.dc_level_delay_symbols = dc_level_delay_symbols; + output.ffe_debug.dc_level_delay_corr = dc_level_corr; + output.ffe_debug.dc_level_lags = dc_level_lags; + output.ffe_debug.dc_level_xcorr = dc_level_xcorr; + fprintf("FFE debug: xcorr(inv DC tracking, avg level) delay=%d symbols, corr=%6.3f\n", ... + dc_level_delay_symbols,dc_level_corr); + + figure(461); clf; + plot(dc_level_lags,dc_level_xcorr,"DisplayName","xcorr"); + hold on; + plot(dc_level_delay_symbols,dc_level_corr,"ro","DisplayName","max"); + grid on; + xlabel("Lag in symbols"); + ylabel("Correlation coefficient"); + title(sprintf("Delay estimate: %d symbols (corr %.3f)",dc_level_delay_symbols,dc_level_corr)); + legend("Location","best"); + + end + + + + %% MPI Reduction if 1 diff --git a/Functions/EQ_visuals/showLevelScatter.m b/Functions/EQ_visuals/showLevelScatter.m index 7b3b62d..c0bbc3d 100644 --- a/Functions/EQ_visuals/showLevelScatter.m +++ b/Functions/EQ_visuals/showLevelScatter.m @@ -1,4 +1,4 @@ -function [symbols_for_lvl, avg_for_lvl, info] = showLevelScatter(rxInput, refSymbols, options) +function [symbols_for_lvl, avg_for_lvl, xAxisUs, info] = showLevelScatter(rxInput, refSymbols, options) %SHOWLEVELSCATTER Plot received samples separated by reference PAM level. % Supports plain numeric vectors, Signal objects, synchronized scope cell % arrays, and raw unsynchronized Signal input. Raw Signal input is diff --git a/projects/Diss/MPI_revisit/investigate_mpi_algorithms.m b/projects/Diss/MPI_revisit/investigate_mpi_algorithms.m index 31f45e4..9d5633d 100644 --- a/projects/Diss/MPI_revisit/investigate_mpi_algorithms.m +++ b/projects/Diss/MPI_revisit/investigate_mpi_algorithms.m @@ -3,9 +3,9 @@ dsp_options = struct(); dsp_options.mode = "run_id"; dsp_options.recipe = @mpi_recipe_dev; dsp_options.append_to_db = false; -dsp_options.start_occurence = 1; +dsp_options.start_occurence = 2; dsp_options.max_occurences = 10; -dsp_options.debug_plots = true; +dsp_options.debug_plots = false; dsp_options.database_type = "mysql"; diff --git a/workerError.mat b/workerError.mat index d0ec048..a8ca784 100644 Binary files a/workerError.mat and b/workerError.mat differ