Compare commits

..

6 Commits

Author SHA1 Message Date
Silas Oettinghaus
8c4edf2490 Final MPI Calculations here 2026-07-14 11:50:22 +02:00
Silas Oettinghaus
f421348e5b before new database 2026-07-13 19:59:57 +02:00
Silas Oettinghaus
ef0a74cb7f work on MPI 2026-07-13 15:25:45 +02:00
Silas Oettinghaus
a6cf742121 MPI FFE with A1 and A2 algorithms 2026-07-10 09:14:34 +02:00
Silas Oettinghaus
aa210e5352 more on MPI tracking 2026-07-09 14:23:16 +02:00
Silas Oettinghaus
753bc582b9 MPI analysis, work on DC tracking 2026-07-09 12:29:41 +02:00
33 changed files with 7064 additions and 609 deletions

View File

@@ -24,9 +24,9 @@ classdef Signal
obj.signal = signal; obj.signal = signal;
obj.signal = obj.signal; obj.signal = obj.signal;
obj.fs = options.fs; obj.fs = options.fs;
%
[~,obj.gitSHA] = system('git rev-parse HEAD'); % [~,obj.gitSHA] = system('git rev-parse HEAD');
[~,obj.gitStatus] = system('git status --porcelain'); % [~,obj.gitStatus] = system('git status --porcelain');
% [~,obj.gitPatch] = system('git diff'); % [~,obj.gitPatch] = system('git diff');
%%% Stuff for Logbook %%% %%% Stuff for Logbook %%%

File diff suppressed because it is too large Load Diff

View File

@@ -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<NASGU>
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

View File

@@ -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<INUSA>
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

View File

@@ -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<NASGU>
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<NASGU>
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

View File

@@ -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<INUSA>
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<NASGU>
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

View File

@@ -0,0 +1,752 @@
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
dc_tracking_mu_eff_max
e_dc
optimize_dc_tracking_params = 0
dc_tracking_optimization
dc_tracking_optimization_iter = 0
dc_tracking_optimization_len
dc_tracking_optimization_max_evals
dc_tracking_optimization_delay_weight
dc_tracking_optimization_smoothing_len
end
properties (Access = private)
dc_tracking_mu_min = 1e-6
dc_tracking_mu_max = 3e-1
dc_tracking_mu_eff_min = 0
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.dc_tracking_mu_eff_max = inf
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
options.optimize_dc_tracking_params = 0
options.dc_tracking_optimization_len = 2^15
options.dc_tracking_optimization_max_evals = 20
options.dc_tracking_optimization_delay_weight = 1e-3
options.dc_tracking_optimization_smoothing_len = 501
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.dc_tracking_mu_eff_max = options.dc_tracking_mu_eff_max;
obj.e_dc = 0;
obj.optimize_dc_tracking_params = options.optimize_dc_tracking_params;
obj.dc_tracking_optimization_len = options.dc_tracking_optimization_len;
obj.dc_tracking_optimization_max_evals = max(1, ...
floor(options.dc_tracking_optimization_max_evals));
obj.dc_tracking_optimization_delay_weight = ...
max(0,options.dc_tracking_optimization_delay_weight);
obj.dc_tracking_optimization_smoothing_len = max(1, ...
floor(options.dc_tracking_optimization_smoothing_len));
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
if obj.optimize_dc_tracking_params
obj.optimizeDcTrackingParams(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 && ...
~obj.optimize_dc_tracking_params;
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 optimizeDcTrackingParams(obj,x,d)
[x_opt,d_opt,N_opt] = obj.optimizationSignals( ...
x,d,obj.dc_tracking_optimization_len);
vars = optimizableVariable( ...
"dc_tracking_mu",[1e-5,1e-1],"Transform","log");
initial_values = min(max(obj.dc_tracking_mu,1e-5),1e-1);
initial_names = "dc_tracking_mu";
if obj.dc_tracking_adaptive_enabled
vars = [vars, ...
optimizableVariable("dc_tracking_persistence_gain",[0,2]), ...
optimizableVariable("dc_tracking_mu_eff_max",[1e-3,3e-1], ...
"Transform","log")];
initial_values = [initial_values, ...
min(max(obj.dc_tracking_persistence_gain,0),2), ...
min(max(obj.dc_tracking_mu_eff_max,1e-3),3e-1)];
initial_names = [initial_names, ...
"dc_tracking_persistence_gain", "dc_tracking_mu_eff_max"];
end
initial_x = array2table(initial_values, ...
"VariableNames",cellstr(initial_names));
obj.dc_tracking_optimization_iter = 0;
fprintf("FFE_DCTracking 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));
old_rng = rng;
cleanup_rng = onCleanup(@()rng(old_rng));
rng(42,"twister");
obj.dc_tracking_optimization = bayesopt( ...
@(p)obj.dcTrackingObjective(p,x_opt,d_opt),vars, ...
"MaxObjectiveEvaluations",obj.dc_tracking_optimization_max_evals, ...
"InitialX",initial_x, ...
"AcquisitionFunctionName","expected-improvement-plus", ...
"IsObjectiveDeterministic",true, ...
"Verbose",0, ...
"PlotFcn",[]);
clear cleanup_rng
best = obj.dc_tracking_optimization.XAtMinObjective;
obj.dc_tracking_mu = best.dc_tracking_mu;
if obj.dc_tracking_adaptive_enabled
obj.dc_tracking_persistence_gain = ...
best.dc_tracking_persistence_gain;
obj.dc_tracking_mu_eff_max = best.dc_tracking_mu_eff_max;
fprintf("\nFFE_DCTracking DC opt done: dc_tracking_mu=%9.3e, persistence_gain=%6.3f, mu_eff_max=%9.3e, objective=%9.3e\n", ...
obj.dc_tracking_mu,obj.dc_tracking_persistence_gain, ...
obj.dc_tracking_mu_eff_max, ...
obj.dc_tracking_optimization.MinObjective);
else
fprintf("\nFFE_DCTracking DC opt done: dc_tracking_mu=%9.3e, objective=%9.3e\n", ...
obj.dc_tracking_mu,obj.dc_tracking_optimization.MinObjective);
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 objective = dcTrackingObjective(obj,params,x,d)
old_state = obj.captureTrackingObjectiveState();
cleanup = onCleanup(@()obj.restoreTrackingObjectiveState(old_state));
obj.applyDcTrackingObjectiveParams(params);
obj.save_debug = 1;
obj.resetTrackingState();
N_tr = min(obj.len_tr,numel(x));
obj.equalize(x,d,obj.mu_tr,obj.epochs_tr,N_tr,true,false);
if obj.dd_mode
[signal,~] = obj.equalize( ...
x,d,obj.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);
[delay_symbols,delay_corr] = obj.dcTrackingDelayObjective(signal,d);
delay_penalty = obj.dc_tracking_optimization_delay_weight * ...
abs(delay_symbols) / max(numel(d),1);
objective = ber + delay_penalty;
if ~isfinite(objective)
objective = inf;
end
obj.dc_tracking_optimization_iter = ...
obj.dc_tracking_optimization_iter + 1;
if obj.dc_tracking_adaptive_enabled
fprintf("\rFFE_DCTracking DC opt %02d: dc_tracking_mu=%9.3e, persistence_gain=%6.3f, mu_eff_max=%9.3e, BER=%9.3e, delay=%7.0f, corr=%6.3f, obj=%9.3e, errors=%d", ...
obj.dc_tracking_optimization_iter,params.dc_tracking_mu, ...
params.dc_tracking_persistence_gain, ...
params.dc_tracking_mu_eff_max,ber,delay_symbols,delay_corr, ...
objective,errors);
else
fprintf("\rFFE_DCTracking DC opt %02d: dc_tracking_mu=%9.3e, BER=%9.3e, delay=%7.0f, corr=%6.3f, obj=%9.3e, errors=%d", ...
obj.dc_tracking_optimization_iter,params.dc_tracking_mu, ...
ber,delay_symbols,delay_corr,objective,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<NASGU>
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;
state.dc_tracking_persistence_gain = ...
obj.dc_tracking_persistence_gain;
state.dc_tracking_mu_eff_max = obj.dc_tracking_mu_eff_max;
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;
obj.dc_tracking_persistence_gain = ...
state.dc_tracking_persistence_gain;
obj.dc_tracking_mu_eff_max = state.dc_tracking_mu_eff_max;
end
function applyDcTrackingObjectiveParams(obj,params)
var_names = string(params.Properties.VariableNames);
if any(var_names == "dc_tracking_mu")
obj.dc_tracking_mu = params.dc_tracking_mu;
end
if any(var_names == "dc_tracking_persistence_gain")
obj.dc_tracking_persistence_gain = ...
params.dc_tracking_persistence_gain;
end
if any(var_names == "dc_tracking_mu_eff_max")
obj.dc_tracking_mu_eff_max = params.dc_tracking_mu_eff_max;
end
end
function [delay_symbols,delay_corr] = ...
dcTrackingDelayObjective(obj,signal,d)
delay_symbols = 0;
delay_corr = 0;
if ~isfield(obj.debug_struct,"dc_tracking_est") || ...
isempty(obj.debug_struct.dc_tracking_est)
return
end
smooth_len = obj.dc_tracking_optimization_smoothing_len;
dc_tracking_est_s = movmean( ...
obj.debug_struct.dc_tracking_est(:),smooth_len,"omitnan");
avg_lvl_dc = obj.averageLevelTrace(signal,d,smooth_len);
xcorr_len = min(numel(dc_tracking_est_s),numel(avg_lvl_dc));
if xcorr_len < 2
return
end
inv_dc_xcorr = -dc_tracking_est_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 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

View File

@@ -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 % FFE Processes signals through FFE equalizer
% %
% Inputs: % Inputs:
@@ -74,16 +74,22 @@ end
% Create FFE results structure % Create FFE results structure
ffe_results = struct(); ffe_results = struct();
try config_clear_props = { ...
eq_.e = []; "e", ...
eq_.e2 = []; "e2", ...
eq_.e3 = []; "e3", ...
eq_.b = []; "b", ...
eq_.b2 = []; "b2", ...
eq_.b3 = []; "b3", ...
"mu_optimization", ...
"dc_optimization", ...
"dc_tracking_optimization", ...
"a2_level_weight_optimization"};
for config_clear_idx = 1:numel(config_clear_props)
prop_name = config_clear_props{config_clear_idx};
if isprop(eq_,prop_name)
eq_.(prop_name) = [];
end end
try
eq_.mu_optimization = [];
end end
ffe_results.config = Equalizerstruct(); ffe_results.config = Equalizerstruct();

View File

@@ -21,7 +21,7 @@ function output = dsp_recipe_minimal(Scpe_sig_raw, Symbols, Tx_bits, options)
"debug_plots", options.debug_plots); "debug_plots", options.debug_plots);
eq_ffe = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-1,"mu_tr",0.4,"order",25,... eq_ffe = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-1,"mu_tr",0.4,"order",25,...
"sps",2,"decide",0,"optmize_mus",0,"dd_mode",options.userParameters.dd_mode,"adaption_technique","nlms","mu_dc",1.021e-05); "sps",2,"decide",0,"optmize_mus",0,"dd_mode",options.userParameters.dd_mode,"adaption_technique","nlms","dc_tracking_mu",1.021e-05);
ffe_results = ffe(eq_ffe, options.M, Scpe_sig, Symbols, Tx_bits, ... ffe_results = ffe(eq_ffe, options.M, Scpe_sig, Symbols, Tx_bits, ...
"precode_mode", options.duob_mode, ... "precode_mode", options.duob_mode, ...

View File

@@ -70,8 +70,8 @@ function output = dsp_scope_signal(Scpe_sig_raw, Symbols, Tx_bits, options)
pf_ = Postfilter("ncoeff", pf_ncoeffs, "useBurg", 1); %#ok<NASGU> pf_ = Postfilter("ncoeff", pf_ncoeffs, "useBurg", 1); %#ok<NASGU>
mlse_ = MLSE("duobinary_output", 0, 'M', M, 'trellis_states', PAMmapper(M,0).levels); %#ok<NASGU> mlse_ = MLSE("duobinary_output", 0, 'M', M, 'trellis_states', PAMmapper(M,0).levels); %#ok<NASGU>
eq_post = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",2001,"sps",1,"decide",0,"adaption_technique","lms","mu_dc",mu_dc); %#ok<NASGU> eq_post = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",2001,"sps",1,"decide",0,"adaption_technique","lms","dc_tracking_mu",mu_dc); %#ok<NASGU>
eq_post = FFE("epochs_tr",5,"epochs_dd",2,"len_tr",2^13,"mu_dd",mu_dd,"mu_tr",mu_tr,"order",25,"sps",2,"decide",0, "adaption_technique",adaption_method(adaption),"dd_mode",use_dd_mode,"mu_dc",mu_dc); %#ok<NASGU> eq_post = FFE("epochs_tr",5,"epochs_dd",2,"len_tr",2^13,"mu_dd",mu_dd,"mu_tr",mu_tr,"order",25,"sps",2,"decide",0, "adaption_technique",adaption_method(adaption),"dd_mode",use_dd_mode,"dc_tracking_mu",mu_dc); %#ok<NASGU>
mlse_db_enc = MLSE("DIR", [1,1], "duobinary_output", 0, "M", M, "trellis_states", PAMmapper(M,0).levels); %#ok<NASGU> mlse_db_enc = MLSE("DIR", [1,1], "duobinary_output", 0, "M", M, "trellis_states", PAMmapper(M,0).levels); %#ok<NASGU>
eq_db_enc = EQ("Ne", ffe_order_dbtgt, "Nb", dfe_order_dbtgt, "training_length", len_tr, ... eq_db_enc = EQ("Ne", ffe_order_dbtgt, "Nb", dfe_order_dbtgt, "training_length", len_tr, ...
@@ -94,7 +94,7 @@ function output = dsp_scope_signal(Scpe_sig_raw, Symbols, Tx_bits, options)
if use_ffe if use_ffe
eq_ffe = EQ("Ne",ffe_order_ffe,"Nb",[0,0,0],"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",0); eq_ffe = EQ("Ne",ffe_order_ffe,"Nb",[0,0,0],"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",0);
eq_ffe = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-1,"mu_tr",0.4,"order",ffe_order_ffe(1),... eq_ffe = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-1,"mu_tr",0.4,"order",ffe_order_ffe(1),...
"sps",2,"decide",0,"optmize_mus",0,"dd_mode",1,"adaption_technique","nlms","mu_dc",1.021e-05); "sps",2,"decide",0,"optmize_mus",0,"dd_mode",1,"adaption_technique","nlms","dc_tracking_mu",1.021e-05);
ffe_results = ffe(eq_ffe, M, Scpe_sig, Symbols, Tx_bits, ... ffe_results = ffe(eq_ffe, M, Scpe_sig, Symbols, Tx_bits, ...
"precode_mode", duob_mode, ... "precode_mode", duob_mode, ...

View File

@@ -0,0 +1,404 @@
function output = mpi_recipe_dev(Scpe_sig_raw, Symbols, Tx_bits, options)
%mpi_recipe_dev Minimal example recipe for the DSP job framework.
% This recipe intentionally performs only light preprocessing and records
% summary values. It demonstrates the recipe interface without running a
% full equalizer chain.
arguments
Scpe_sig_raw
Symbols
Tx_bits
options.fsym
options.M
options.duob_mode
options.dataTable table
options.userParameters struct = struct()
options.debug_plots (1,1) logical = false
end
Scpe_sig = preprocessSignal(Scpe_sig_raw, Symbols, options.fsym, ...
"mode", "auto", ...
"debug_plots", options.debug_plots);
output = struct();
%% conventional FFE
eq_core_settings = { ...
"sps", 2, ...
"order", 25, ...
"decide", 0, ...
"adaption_technique", "nlms"};
eq_training_settings = { ...
"len_tr", 4096, ...
"epochs_tr", 5, ...
"mu_tr", 0.04};
eq_dd_settings = { ...
"dd_mode", 1, ...
"epochs_dd", 3, ...
"mu_dd", 0.012};
eq_a1_settings = { ...
"dc_smoothing_a1", 0, ...
"dc_avg_bufferlength_a1", 0, ...
"dc_avg_update_blocklength_a1", 1};
eq_a2_settings = { ...
"dc_smoothing_a2", 0, ...
"dc_level_avg_bufferlength_a2", 0, ...
"dc_level_update_blocklength_a2", 0, ...
"dc_level_weights_a2", [0.6 0.6 0.6 0.6]};
eq_dc_tracking_settings = { ...
"dc_tracking_mu", 0, ...
"dc_tracking_adaptive_enabled", 0, ...
"dc_tracking_persistence_gain", 0, ...
"dc_tracking_buffer_len", 0};
eq_ffe_update_settings = { ...
"ffe_update_buffer_len", 1};
eq_optimizer_settings = { ...
"optmize_mus", 0, ...
"plot_mu_optimization", options.debug_plots, ...
"optimize_dc_tracking_params", 0, ...
"optimize_a2_level_weights", 0, ...
"a2_level_weight_optimization_len", 2^15, ...
"a2_level_weight_optimization_max_evals", 30, ...
"a2_level_weight_max", 1};
eq_debug_settings = { ...
"save_debug", false};
eq_settings = [ ...
eq_core_settings, ...
eq_training_settings, ...
eq_dd_settings, ...
eq_a1_settings, ...
eq_a2_settings, ...
eq_dc_tracking_settings, ...
eq_ffe_update_settings, ...
eq_optimizer_settings, ...
eq_debug_settings];
eq_ffe = FFE(eq_settings{:});
if 1
[ffe_results,equalized_signal] = ffe(eq_ffe, options.M, Scpe_sig, Symbols, Tx_bits, ...
"precode_mode", options.duob_mode, ...
'showAnalysis', options.debug_plots, ...
"postFFE", [], ...
"eth_style_symbol_mapping", 0);
ffe_results.metrics.print("description",sprintf('Immediate A2; SIR %d dB',options.dataTable.sir));
output.conventional_ffe = ffe_results;
end
if options.debug_plots
showLevelScatter(Scpe_sig_raw, Symbols, ...
"fsym", options.fsym, ...
"fignum", 400, ...
"normalize", true);
[~, avg_for_lvl] = showLevelScatter(equalized_signal, Symbols, ...
"fsym", options.fsym, ...
"fignum", 401, ...
"normalize", true);
end
%% A2
eq_core_settings = { ...
"sps", 2, ...
"order", 25, ...
"decide", 0, ...
"adaption_technique", "nlms"};
eq_training_settings = { ...
"len_tr", 4096, ...
"epochs_tr", 5, ...
"mu_tr", 0.04};
eq_dd_settings = { ...
"dd_mode", 1, ...
"epochs_dd", 3, ...
"mu_dd", 0.012};
eq_a1_settings = { ...
"dc_smoothing_a1", 0, ...
"dc_avg_bufferlength_a1", 0, ...
"dc_avg_update_blocklength_a1", 1};
eq_a2_settings = { ...
"dc_smoothing_a2", 1, ...
"dc_level_avg_bufferlength_a2", 256, ...
"dc_level_update_blocklength_a2", options.userParameters.block_update, ...
"dc_level_weights_a2", [0.6 0.6 0.6 0.6]};
eq_dc_tracking_settings = { ...
"dc_tracking_mu", 0, ...
"dc_tracking_adaptive_enabled", 0, ...
"dc_tracking_persistence_gain", 0, ...
"dc_tracking_buffer_len", 0};
eq_ffe_update_settings = { ...
"ffe_update_buffer_len", 1};
eq_optimizer_settings = { ...
"optmize_mus", 0, ...
"plot_mu_optimization", options.debug_plots, ...
"optimize_dc_tracking_params", 0, ...
"optimize_a2_level_weights", 1, ...
"a2_level_weight_optimization_len", 2^15, ...
"a2_level_weight_optimization_max_evals", 30, ...
"a2_level_weight_max", 1};
eq_debug_settings = { ...
"save_debug", false};
eq_settings = [ ...
eq_core_settings, ...
eq_training_settings, ...
eq_dd_settings, ...
eq_a1_settings, ...
eq_a2_settings, ...
eq_dc_tracking_settings, ...
eq_ffe_update_settings, ...
eq_optimizer_settings, ...
eq_debug_settings];
eq_ffe = FFE(eq_settings{:});
if 1
[ffe_results,equalized_signal] = ffe(eq_ffe, options.M, Scpe_sig, Symbols, Tx_bits, ...
"precode_mode", options.duob_mode, ...
'showAnalysis', options.debug_plots, ...
"postFFE", [], ...
"eth_style_symbol_mapping", 0);
ffe_results.metrics.print("description",sprintf('Immediate A2; SIR %d dB',options.dataTable.sir));
output.a2_immediate_updates = ffe_results;
end
if options.debug_plots
showLevelScatter(Scpe_sig_raw, Symbols, ...
"fsym", options.fsym, ...
"fignum", 400, ...
"normalize", true);
[~, avg_for_lvl] = showLevelScatter(equalized_signal, Symbols, ...
"fsym", options.fsym, ...
"fignum", 401, ...
"normalize", true);
end
%% A1
eq_a1_settings = { ...
"dc_smoothing_a1", 1 ...
"dc_avg_bufferlength_a1", 2048, ...
"dc_avg_update_blocklength_a1", options.userParameters.block_update};
eq_a2_settings = { ...
"dc_smoothing_a2", 0, ...
"dc_level_avg_bufferlength_a2", 0, ...
"dc_level_update_blocklength_a2", 0, ...
"dc_level_weights_a2", [0.6 0.6 0.6 0.6]};
eq_optimizer_settings = { ...
"optmize_mus", 0, ...
"plot_mu_optimization", options.debug_plots, ...
"optimize_dc_tracking_params", 0, ...
"optimize_a2_level_weights", 0, ...
"a2_level_weight_optimization_len", 2^15, ...
"a2_level_weight_optimization_max_evals", 30, ...
"a2_level_weight_max", 1};
eq_debug_settings = { ...
"save_debug", false};
eq_settings = [ ...
eq_core_settings, ...
eq_training_settings, ...
eq_dd_settings, ...
eq_a1_settings, ...
eq_a2_settings, ...
eq_dc_tracking_settings, ...
eq_ffe_update_settings, ...
eq_optimizer_settings, ...
eq_debug_settings];
eq_ffe = FFE(eq_settings{:});
if 1
[ffe_results,equalized_signal] = ffe(eq_ffe, options.M, Scpe_sig, Symbols, Tx_bits, ...
"precode_mode", options.duob_mode, ...
'showAnalysis', options.debug_plots, ...
"postFFE", [], ...
"eth_style_symbol_mapping", 0);
ffe_results.metrics.print("description",sprintf('Immediate A1; SIR %d dB',options.dataTable.sir));
output.a1_immediate_updates = ffe_results;
end
%% Tracking
eq_core_settings = { ...
"sps", 2, ...
"order", 25, ...
"decide", 0, ...
"adaption_technique", "nlms"};
eq_training_settings = { ...
"len_tr", 4096, ...
"epochs_tr", 5, ...
"mu_tr", 0.04};
eq_dd_settings = { ...
"dd_mode", 1, ...
"epochs_dd", 3, ...
"mu_dd", 0.012};
eq_a1_settings = { ...
"dc_smoothing_a1", 0, ...
"dc_avg_bufferlength_a1", 0, ...
"dc_avg_update_blocklength_a1", 1};
eq_a2_settings = { ...
"dc_smoothing_a2", 0, ...
"dc_level_avg_bufferlength_a2", 0, ...
"dc_level_update_blocklength_a2", 0, ...
"dc_level_weights_a2", [0.6 0.6 0.6 0.6]};
eq_dc_tracking_settings = { ...
"dc_tracking_mu", 0.002, ...
"dc_tracking_adaptive_enabled", 0, ...
"dc_tracking_persistence_gain", 0, ...
"dc_tracking_buffer_len", options.userParameters.block_update};
eq_ffe_update_settings = { ...
"ffe_update_buffer_len", 1};
eq_optimizer_settings = { ...
"optmize_mus", 0, ...
"plot_mu_optimization", options.debug_plots, ...
"optimize_dc_tracking_params", 1, ...
"optimize_a2_level_weights", 0, ...
"a2_level_weight_optimization_len", 2^15, ...
"a2_level_weight_optimization_max_evals", 30, ...
"a2_level_weight_max", 1};
eq_debug_settings = { ...
"save_debug", false};
eq_settings = [ ...
eq_core_settings, ...
eq_training_settings, ...
eq_dd_settings, ...
eq_a1_settings, ...
eq_a2_settings, ...
eq_dc_tracking_settings, ...
eq_ffe_update_settings, ...
eq_optimizer_settings, ...
eq_debug_settings];
eq_ffe = FFE(eq_settings{:});
if 1
[ffe_results,equalized_signal] = ffe(eq_ffe, options.M, Scpe_sig, Symbols, Tx_bits, ...
"precode_mode", options.duob_mode, ...
'showAnalysis', options.debug_plots, ...
"postFFE", [], ...
"eth_style_symbol_mapping", 0);
ffe_results.metrics.print("description",sprintf('Immediate A2; SIR %d dB',options.dataTable.sir));
output.tracking_fixmu_immediate_updates = ffe_results;
end
%% Tracking Adaptive
eq_core_settings = { ...
"sps", 2, ...
"order", 25, ...
"decide", 0, ...
"adaption_technique", "nlms"};
eq_training_settings = { ...
"len_tr", 4096, ...
"epochs_tr", 5, ...
"mu_tr", 0.04};
eq_dd_settings = { ...
"dd_mode", 1, ...
"epochs_dd", 3, ...
"mu_dd", 0.012};
eq_a1_settings = { ...
"dc_smoothing_a1", 0, ...
"dc_avg_bufferlength_a1", 0, ...
"dc_avg_update_blocklength_a1", 1};
eq_a2_settings = { ...
"dc_smoothing_a2", 0, ...
"dc_level_avg_bufferlength_a2", 0, ...
"dc_level_update_blocklength_a2", 0, ...
"dc_level_weights_a2", [0.6 0.6 0.6 0.6]};
eq_dc_tracking_settings = { ...
"dc_tracking_mu", 0.002, ...
"dc_tracking_adaptive_enabled", 1, ...
"dc_tracking_persistence_gain", 0, ...
"dc_tracking_buffer_len", options.userParameters.block_update};
eq_ffe_update_settings = { ...
"ffe_update_buffer_len", 1};
eq_optimizer_settings = { ...
"optmize_mus", 0, ...
"plot_mu_optimization", options.debug_plots, ...
"optimize_dc_tracking_params", 1, ...
"optimize_a2_level_weights", 0, ...
"a2_level_weight_optimization_len", 2^15, ...
"a2_level_weight_optimization_max_evals", 30, ...
"a2_level_weight_max", 1};
eq_debug_settings = { ...
"save_debug", false};
eq_settings = [ ...
eq_core_settings, ...
eq_training_settings, ...
eq_dd_settings, ...
eq_a1_settings, ...
eq_a2_settings, ...
eq_dc_tracking_settings, ...
eq_ffe_update_settings, ...
eq_optimizer_settings, ...
eq_debug_settings];
eq_ffe = FFE(eq_settings{:});
if 1
[ffe_results,equalized_signal] = ffe(eq_ffe, options.M, Scpe_sig, Symbols, Tx_bits, ...
"precode_mode", options.duob_mode, ...
'showAnalysis', options.debug_plots, ...
"postFFE", [], ...
"eth_style_symbol_mapping", 0);
ffe_results.metrics.print("description",sprintf('Immediate A2; SIR %d dB',options.dataTable.sir));
output.tracking_adaptive_immediate_updates = ffe_results;
end
end

View File

@@ -1,8 +1,5 @@
function output = mpi_recipe_dev(Scpe_sig_raw, Symbols, Tx_bits, options) function output = mpi_recipe_dev(Scpe_sig_raw, Symbols, Tx_bits, options)
%mpi_recipe_dev Minimal example recipe for the DSP job framework. %mpi_recipe_dev Minimal MPI equalizer comparison recipe.
% This recipe intentionally performs only light preprocessing and records
% summary values. It demonstrates the recipe interface without running a
% full equalizer chain.
arguments arguments
Scpe_sig_raw Scpe_sig_raw
@@ -16,77 +13,277 @@ function output = mpi_recipe_dev(Scpe_sig_raw, Symbols, Tx_bits, options)
options.debug_plots (1,1) logical = false options.debug_plots (1,1) logical = false
end 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, ... Scpe_sig = preprocessSignal(Scpe_sig_raw, Symbols, options.fsym, ...
"mode", "auto", ... "mode", "auto", ...
"debug_plots", options.debug_plots); "debug_plots", options.debug_plots);
output = struct();
%%
mu_dc = 0;%1e-5; if plot_output_signals
showLevelScatter(Scpe_sig, Symbols, ...
eq_settings = { ...
"epochs_tr", 5, ...
"epochs_dd", 5, ...
"len_tr", 4096*2, ...
"mu_dd",0.02, ...
"mu_tr",0.2, ...
"order", 25, ...
"sps", 2, ...
"decide", 0, ...
"optmize_mus", 1, ...
"dd_mode", 1, ...
"adaption_technique", "nlms", ...
"mu_dc", mu_dc};
eq_ffe = FFE(eq_settings{:});
showLevelScatter(Scpe_sig_raw, Symbols, ...
"fsym", options.fsym, ... "fsym", options.fsym, ...
"fignum", 400, ... "fignum", 400, ...
"normalize", true); "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);
% Scpe_sig_raw.spectrum("normalizeTo0dB",1,"fft_length",4096*4,"fignum",401); storageName = "conventional_ffe";
[ffe_results,equalized_signal] = runFfe(eq_ffe, "Conventional FFE", ...
Scpe_sig, Symbols, Tx_bits, options);
ffe_results = attachMpiReductionConfig(ffe_results, eq_ffe, storageName, ...
"plain_ffe", "baseline", block_update);
output.(char(storageName)) = ffe_results;
%options.dataTable.sir; if plot_output_signals
plotEqSignals(equalized_signal,Symbols,options,400,-1);
end
end
%% NORMAL FFE %% A2 tracked-level decision
if 1 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", 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", [1], ...
"save_debug", eq_save_debug);
ffe_results = ffe(eq_ffe, options.M, Scpe_sig, Symbols, Tx_bits, ... storageName = "a2_adaptive_levels";
[ffe_results,equalized_signal] = runFfe(eq_ffe, "A2 tracked levels", ...
Scpe_sig, Symbols, Tx_bits, options);
ffe_results = attachMpiReductionConfig(ffe_results, eq_ffe, storageName, ...
"a2_tracked_levels", "tracked_levels", block_update);
output.(char(storageName)) = ffe_results;
if plot_output_signals
plotEqSignals(equalized_signal,Symbols,options,410,-1);
end
end
%% 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);
storageName = "a2_residual";
[ffe_results,equalized_signal] = runFfe(eq_ffe, "A2 residual", ...
Scpe_sig, Symbols, Tx_bits, options);
ffe_results = attachMpiReductionConfig(ffe_results, eq_ffe, storageName, ...
"a2_residual", "residual_correction", block_update);
output.(char(storageName)) = 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);
storageName = "a1_ff_dc_avg";
[ffe_results,equalized_signal] = runFfe(eq_ffe, "A1", ...
Scpe_sig, Symbols, Tx_bits, options);
ffe_results = attachMpiReductionConfig(ffe_results, eq_ffe, storageName, ...
"a1_moving_average", "ff_dc_avg", block_update);
output.(char(storageName)) = 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, ...
"optimize_dc_tracking_params", true, ...
"dc_tracking_optimization_len", 2^15, ...
"dc_tracking_optimization_max_evals", 20, ...
"plot_mu_optimization", options.debug_plots, ...
"save_debug", eq_save_debug);
storageName = "dc_tracking";
[ffe_results,equalized_signal] = runFfe(eq_ffe, "DC tracking adaptive", ...
Scpe_sig, Symbols, Tx_bits, options);
ffe_results = attachMpiReductionConfig(ffe_results, eq_ffe, storageName, ...
"dc_tracking", "adaptive", block_update);
output.(char(storageName)) = ffe_results;
if plot_output_signals
plotEqSignals(equalized_signal,Symbols,options,450,-1);
end
end
end
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, ... "precode_mode", options.duob_mode, ...
'showAnalysis', options.debug_plots, ... "showAnalysis", options.debug_plots, ...
"postFFE", [], ... "postFFE", [], ...
"eth_style_symbol_mapping", 0); "eth_style_symbol_mapping", 0);
ffe_results.metrics.print("description",sprintf('Normal FFE; SIR %d dB',options.dataTable.sir)); ffe_results.metrics.print("description",resultDescription(description,options));
output.ffe_package = ffe_results;
end end
%% MPI Reduction function ffe_results = attachMpiReductionConfig(ffe_results, eq_ffe, storageName, ...
if 1 algorithm, algorithmVariant, block_update)
ffe_results.mpi_reduction_config = struct( ...
"storage_name", char(storageName), ...
"algorithm", char(algorithm), ...
"algorithm_variant", char(algorithmVariant), ...
"eq_class", class(eq_ffe), ...
"params", collectMpiReductionParams(eq_ffe, block_update));
end
% dc_buffer_len = 0; function params = collectMpiReductionParams(eq_ffe, block_update)
% ffe_buffer_len = 0; params = struct();
% smoothing_buffer_length = options.userParameters.smoothing_length; params.block_update = block_update;
% smoothing_buffer_update = 1;
% eq_ffe_dcr = FFE_DCremoval_adaptive_mu(eq_settings{:}, ... whitelistedProps = [ ...
% "dc_buffer_len",dc_buffer_len, ... "sps", ...
% "ffe_buffer_len",ffe_buffer_len,... "order", ...
% "smoothing_buffer_length",smoothing_buffer_length,... "decide", ...
% "smoothing_buffer_update",smoothing_buffer_update); "adaption_technique", ...
"len_tr", ...
"epochs_tr", ...
"mu_tr", ...
"dd_mode", ...
"epochs_dd", ...
"mu_dd", ...
"optmize_mus", ...
"plot_mu_optimization", ...
"save_debug", ...
"dc_smoothing_a1", ...
"dc_avg_bufferlength_a1", ...
"dc_avg_update_blocklength_a1", ...
"dc_smoothing_a2", ...
"dc_level_avg_bufferlength_a2", ...
"dc_level_update_blocklength_a2", ...
"dc_level_weights_a2", ...
"dc_tracking_mu", ...
"dc_tracking_adaptive_enabled", ...
"dc_tracking_persistence_gain", ...
"dc_tracking_buffer_len", ...
"dc_tracking_mu_eff_max", ...
"optimize_dc_tracking_params", ...
"dc_tracking_optimization_len", ...
"dc_tracking_optimization_max_evals", ...
"dc_tracking_optimization_delay_weight", ...
"dc_tracking_optimization_smoothing_len"];
eq_ = FFE_DCremoval("epochs_tr",5,"epochs_dd",3,"len_tr",4096*2,"mu_dd",... for propIdx = 1:numel(whitelistedProps)
0.0002,"mu_tr",0,"order",25,"sps",2,"decide",0,... propName = char(whitelistedProps(propIdx));
"mu_dc",0.005,"dc_buffer_len",1); if isprop(eq_ffe, propName)
params.(propName) = eq_ffe.(propName);
ffe_results_dcr = ffe(eq_, options.M, Scpe_sig, Symbols, Tx_bits, ...
"precode_mode", options.duob_mode, ...
'showAnalysis', options.debug_plots, ...
"postFFE", [], ...
"eth_style_symbol_mapping", 0);
ffe_results_dcr.metrics.print("description",'FFE DCR');
output.ffe_dcr_package = ffe_results_dcr;
end end
end end
end
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
function plotEqSignals(equalized_signal,Symbols,options,fignum,output_scale)
showLevelScatter(equalized_signal .* output_scale, Symbols, ...
"fsym", options.fsym, ...
"fignum", fignum + 1, ...
"normalize", true);
end

View File

@@ -0,0 +1,320 @@
function exportLevelScatterTikz(sourceAx, pngFile, tikzFile, pngTikzPath, options)
%EXPORTLEVELSCATTERTIKZ Export showLevelScatter axes as PNG-backed TikZ.
% Rasterizes scatter objects into a PNG layer and writes average lines plus
% text annotations as PGFPlots objects.
arguments
sourceAx
pngFile (1,1) string
tikzFile (1,1) string
pngTikzPath (1,1) string
options.resolutionDpi (1,1) double {mustBePositive, mustBeInteger} = 150
options.xLabel (1,1) string = "Time in $\mu$s"
options.yLabel (1,1) string = "Normalized Amplitude"
options.generatedBy (1,1) string = "exportLevelScatterTikz.m"
end
outputDir = fileparts(pngFile);
if ~exist(outputDir, "dir")
mkdir(outputDir);
end
exportScatterLayerPng(sourceAx, pngFile, options.resolutionDpi);
linePlots = collectTikzLinePlots(sourceAx);
textAnnotations = collectTikzTextAnnotations(sourceAx);
writeTikzPngAxis(tikzFile, pngTikzPath, sourceAx, linePlots, textAnnotations, options);
end
function exportScatterLayerPng(sourceAx, pngFile, resolutionDpi)
sourceFig = ancestor(sourceAx, "figure");
xLimits = sourceAx.XLim;
yLimits = sourceAx.YLim;
scatterHandles = findGraphicsObjectsByType(sourceAx, "scatter");
exportFig = figure( ...
"Visible", "off", ...
"Color", "white", ...
"Units", sourceFig.Units, ...
"Position", sourceFig.Position);
cleanupFigure = onCleanup(@() close(exportFig));
exportAx = copyobj(sourceAx, exportFig);
exportAx.Units = "normalized";
exportAx.Position = [0 0 1 1];
exportAx.XLim = xLimits;
exportAx.YLim = yLimits;
exportAx.XLimMode = "manual";
exportAx.YLimMode = "manual";
stripAxesToScatterOnly(exportAx);
hideAxesInkForRasterExport(exportAx);
fprintf("Exporting scatter PNG with XLim=[%g %g], YLim=[%g %g], scatter objects=%d: %s\n", ...
xLimits(1), xLimits(2), yLimits(1), yLimits(2), numel(scatterHandles), pngFile);
drawnow;
exportFig.PaperPositionMode = "auto";
print(exportFig, char(pngFile), "-dpng", sprintf("-r%d", resolutionDpi));
end
function stripAxesToScatterOnly(ax)
plotObjects = findall(ax);
for objIdx = 1:numel(plotObjects)
obj = plotObjects(objIdx);
if obj == ax || ~isvalid(obj)
continue
end
objType = string(get(obj, "Type"));
if lower(objType) ~= "scatter"
delete(obj);
end
end
end
function hideAxesInkForRasterExport(ax)
ax.Visible = "on";
ax.Color = "white";
ax.Box = "off";
ax.XGrid = "off";
ax.YGrid = "off";
ax.XMinorGrid = "off";
ax.YMinorGrid = "off";
ax.XTick = [];
ax.YTick = [];
ax.XColor = "white";
ax.YColor = "white";
ax.Title.String = "";
ax.XLabel.String = "";
ax.YLabel.String = "";
end
function handles = findGraphicsObjectsByType(parentHandle, objectType)
allHandles = findall(parentHandle);
matches = false(size(allHandles));
for handleIdx = 1:numel(allHandles)
currentHandle = allHandles(handleIdx);
if ~isvalid(currentHandle)
continue
end
matches(handleIdx) = lower(string(get(currentHandle, "Type"))) == lower(string(objectType));
end
handles = allHandles(matches);
end
function linePlots = collectTikzLinePlots(sourceAx)
lineHandles = findGraphicsObjectsByType(sourceAx, "line");
lineHandles = flipud(lineHandles(:));
linePlots = repmat(struct( ...
"xData", [], ...
"yData", [], ...
"color", [], ...
"lineWidth", []), numel(lineHandles), 1);
validLineCount = 0;
for lineIdx = 1:numel(lineHandles)
lineHandle = lineHandles(lineIdx);
xData = lineHandle.XData(:);
yData = lineHandle.YData(:);
validSamples = isfinite(xData) & isfinite(yData);
if ~any(validSamples)
continue
end
validLineCount = validLineCount + 1;
linePlots(validLineCount).xData = xData(validSamples);
linePlots(validLineCount).yData = yData(validSamples);
linePlots(validLineCount).color = lineHandle.Color;
linePlots(validLineCount).lineWidth = lineHandle.LineWidth;
end
linePlots = linePlots(1:validLineCount);
end
function textAnnotations = collectTikzTextAnnotations(sourceAx)
textHandles = findGraphicsObjectsByType(sourceAx, "text");
textHandles = flipud(textHandles(:));
xLimits = sourceAx.XLim;
yLimits = sourceAx.YLim;
textAnnotations = repmat(struct( ...
"x", [], ...
"y", [], ...
"text", "", ...
"anchor", "center"), numel(textHandles), 1);
validTextCount = 0;
for textIdx = 1:numel(textHandles)
textHandle = textHandles(textIdx);
labelText = string(textHandle.String);
if strlength(strtrim(labelText)) == 0
continue
end
textPosition = textHandle.Position;
if textPosition(1) < xLimits(1) || textPosition(1) > xLimits(2) || ...
textPosition(2) < yLimits(1) || textPosition(2) > yLimits(2)
continue
end
validTextCount = validTextCount + 1;
textAnnotations(validTextCount).x = textPosition(1);
textAnnotations(validTextCount).y = textPosition(2);
textAnnotations(validTextCount).text = labelText;
textAnnotations(validTextCount).anchor = matlabTextAlignmentToTikzAnchor( ...
textHandle.HorizontalAlignment, textHandle.VerticalAlignment);
end
textAnnotations = textAnnotations(1:validTextCount);
end
function writeTikzPngAxis(tikzFile, pngTikzPath, sourceAx, linePlots, textAnnotations, options)
outputDir = fileparts(tikzFile);
if ~exist(outputDir, "dir")
mkdir(outputDir);
end
xLimits = sourceAx.XLim;
yLimits = sourceAx.YLim;
fid = fopen(tikzFile, "w");
if fid < 0
error("exportLevelScatterTikz:TikzOpenFailed", ...
"Could not open TikZ export file: %s", tikzFile);
end
cleanupFile = onCleanup(@() fclose(fid));
axisOptions = {
" every axis/.append style={font=\scriptsize},"
"width=\fwidth,"
"height=\fheight,"
"at={(0\fwidth,0\fheight)},"
"scale only axis,"
"axis on top,"
"xmin=" + formatPgfNumber(xLimits(1)) + ","
"xmax=" + formatPgfNumber(xLimits(2)) + ","
"xlabel style={font=\color{white!15!black}\small},"
"xlabel={" + options.xLabel + "},"
"ymin=" + formatPgfNumber(yLimits(1)) + ","
"ymax=" + formatPgfNumber(yLimits(2)) + ","
"ylabel style={font=\color{white!15!black}\small},"
"ylabel={" + options.yLabel + "},"
"axis background/.style={fill=white},"
"xmajorgrids,"
"ymajorgrids,"
"grid style={line width=0.4pt, dotted, color=black!20},"
"scaled ticks=false,"
"tick label style={/pgf/number format/fixed, /pgf/number format/1000 sep={}},"
"legend columns=1"
};
fprintf(fid, "%% This file was generated by %s.\n", options.generatedBy);
fprintf(fid, "%%\n");
for lineIdx = 1:numel(linePlots)
fprintf(fid, "%s\n", formatTikzColorDefinition(lineIdx, linePlots(lineIdx).color));
end
if ~isempty(linePlots)
fprintf(fid, "%%\n");
end
fprintf(fid, "%s\n\n", "\begin{tikzpicture}");
fprintf(fid, "%s\n", "\begin{axis}[%");
for optionIdx = 1:numel(axisOptions)
fprintf(fid, "%s\n", axisOptions{optionIdx});
end
fprintf(fid, "%s\n", "]");
graphicsLine = sprintf( ...
"\\addplot [forget plot] graphics [xmin=%s, xmax=%s, ymin=%s, ymax=%s] {%s};", ...
formatPgfNumber(xLimits(1)), ...
formatPgfNumber(xLimits(2)), ...
formatPgfNumber(yLimits(1)), ...
formatPgfNumber(yLimits(2)), ...
char(strrep(pngTikzPath, "\", "/")));
fprintf(fid, "%s\n\n", graphicsLine);
writeTikzLinePlots(fid, linePlots);
writeTikzTextAnnotations(fid, textAnnotations);
fprintf(fid, "%s\n\n", "\end{axis}");
fprintf(fid, "%s", "\end{tikzpicture}%");
end
function valueText = formatPgfNumber(value)
valueText = string(sprintf("%.15g", value));
end
function colorDefinition = formatTikzColorDefinition(colorIdx, rgbColor)
colorDefinition = sprintf( ...
"\\definecolor{mycolor%d}{rgb}{%.5f,%.5f,%.5f}%%", ...
colorIdx, rgbColor(1), rgbColor(2), rgbColor(3));
end
function writeTikzLinePlots(fid, linePlots)
for lineIdx = 1:numel(linePlots)
fprintf(fid, "\\addplot [color=mycolor%d, line width=%.1fpt, forget plot]\n", ...
lineIdx, linePlots(lineIdx).lineWidth);
fprintf(fid, " table[row sep=crcr]{%%\n");
for sampleIdx = 1:numel(linePlots(lineIdx).xData)
fprintf(fid, "%s\t%s\\\\\n", ...
formatPgfNumber(linePlots(lineIdx).xData(sampleIdx)), ...
formatPgfNumber(linePlots(lineIdx).yData(sampleIdx)));
end
fprintf(fid, "};\n\n");
end
end
function writeTikzTextAnnotations(fid, textAnnotations)
for textIdx = 1:numel(textAnnotations)
fprintf(fid, "\\node[font=\\scriptsize, anchor=%s, fill=white, draw=black!40, rounded corners=2pt, inner sep=2pt]\n", ...
textAnnotations(textIdx).anchor);
fprintf(fid, " at (axis cs:%s,%s) {%s};\n\n", ...
formatPgfNumber(textAnnotations(textIdx).x), ...
formatPgfNumber(textAnnotations(textIdx).y), ...
escapeTikzText(textAnnotations(textIdx).text));
end
end
function anchor = matlabTextAlignmentToTikzAnchor(horizontalAlignment, verticalAlignment)
horizontalAlignment = string(horizontalAlignment);
verticalAlignment = string(verticalAlignment);
if horizontalAlignment == "left"
horizontalAnchor = "west";
elseif horizontalAlignment == "right"
horizontalAnchor = "east";
else
horizontalAnchor = "";
end
if verticalAlignment == "top"
verticalAnchor = "north";
elseif verticalAlignment == "bottom"
verticalAnchor = "south";
else
verticalAnchor = "";
end
if strlength(horizontalAnchor) > 0 && strlength(verticalAnchor) > 0
anchor = verticalAnchor + " " + horizontalAnchor;
elseif strlength(horizontalAnchor) > 0
anchor = horizontalAnchor;
elseif strlength(verticalAnchor) > 0
anchor = verticalAnchor;
else
anchor = "center";
end
end
function textOut = escapeTikzText(textIn)
textOut = char(textIn);
% textOut = strrep(textOut, "\", "\textbackslash{}");
textOut = strrep(textOut, "{", "\{");
textOut = strrep(textOut, "}", "\}");
textOut = strrep(textOut, "%", "\%");
textOut = strrep(textOut, "&", "\&");
textOut = strrep(textOut, "#", "\#");
textOut = strrep(textOut, "_", "\_");
% textOut = strrep(textOut, "$", "\$");
end

View File

@@ -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. %SHOWLEVELSCATTER Plot received samples separated by reference PAM level.
% Supports plain numeric vectors, Signal objects, synchronized scope cell % Supports plain numeric vectors, Signal objects, synchronized scope cell
% arrays, and raw unsynchronized Signal input. Raw Signal input is % arrays, and raw unsynchronized Signal input. Raw Signal input is
@@ -290,7 +290,7 @@ for levelIdx = 1:numLevels
xAxisUs, avg_for_lvl(levelIdx, :), ... xAxisUs, avg_for_lvl(levelIdx, :), ...
options.avgLineMaxPoints, options.avgLineSmoothWindow); options.avgLineMaxPoints, options.avgLineSmoothWindow);
plot(ax, xReduced, yReduced, ... plot(ax, xReduced, yReduced, ...
"LineWidth", 2, ... "LineWidth", 1, ...
"Color", cols(levelIdx, :)); "Color", cols(levelIdx, :));
end end

View File

@@ -3,6 +3,9 @@ function output = dsp_runid(run_id, options)
arguments arguments
run_id run_id
options.append_to_db = 0; options.append_to_db = 0;
options.append_mpi_reduction_db (1,1) logical = false;
options.mpi_reduction_study_name string = "mpi_reduction_v1";
options.mpi_reduction_writer_path string = "";
options.max_occurences = 4; options.max_occurences = 4;
options.start_occurence = 1; options.start_occurence = 1;
options.userParameters = struct(); options.userParameters = struct();
@@ -24,7 +27,11 @@ try
database = []; database = [];
inputSource = normalizeDspInputSource(options.mode); inputSource = normalizeDspInputSource(options.mode);
if inputSource == "run_id" || options.append_to_db if options.append_mpi_reduction_db && strlength(options.mpi_reduction_writer_path) > 0
addpath(options.mpi_reduction_writer_path);
end
if inputSource == "run_id" || options.append_to_db || options.append_mpi_reduction_db
database = DBHandler("dataBase", [options.dataBase], "type", options.database_type, ... database = DBHandler("dataBase", [options.dataBase], "type", options.database_type, ...
"user", options.user, "password", options.password, ... "user", options.user, "password", options.password, ...
"server", options.server, "port", options.port); "server", options.server, "port", options.port);
@@ -60,6 +67,12 @@ try
if options.append_to_db if options.append_to_db
appendDspOutputToDatabase(database, run_id, dspOutput); appendDspOutputToDatabase(database, run_id, dspOutput);
end end
if options.append_mpi_reduction_db
occurrence_idx = options.start_occurence + r - 1;
appendMpiReductionDspOutput(database, run_id, occurrence_idx, dspOutput, options, ...
"study_name", options.mpi_reduction_study_name, ...
"verbose", false);
end
end end
catch ME catch ME

View File

@@ -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

View File

@@ -12,6 +12,23 @@ classdef FFE_test < IMDDTestCase
testCase.verifyEqual(ffe.dd_mode, 1); testCase.verifyEqual(ffe.dd_mode, 1);
testCase.verifyEqual(ffe.mu_dd, 1e-5); testCase.verifyEqual(ffe.mu_dd, 1e-5);
testCase.verifyEqual(ffe.epochs_dd, 5); testCase.verifyEqual(ffe.epochs_dd, 5);
testCase.verifyFalse(logical(ffe.optimize_dc_tracking_params));
testCase.verifyEqual(ffe.dc_tracking_optimization_delay_weight, 1e-3);
testCase.verifyEqual(ffe.dc_tracking_persistence_gain, 0);
testCase.verifyEqual(ffe.dc_tracking_power_exponent, 2);
testCase.verifyEqual(ffe.dc_avg_bufferlength_a1, 0);
testCase.verifyEqual(ffe.dc_avg_update_blocklength_a1, 0);
testCase.verifyEqual(ffe.dc_smoothing_a1, 0);
testCase.verifyEqual(ffe.dc_level_avg_bufferlength_a2, 0);
testCase.verifyEqual(ffe.dc_level_update_blocklength_a2, 0);
testCase.verifyEqual(ffe.dc_smoothing_a2, 0);
testCase.verifyEqual(ffe.dc_level_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);
testCase.verifyFalse(logical(ffe.optimize_a2_level_weights));
testCase.verifyEqual(ffe.a2_level_weight_optimization_max_evals, 30);
testCase.verifyEqual(ffe.a2_level_weight_max, 1);
testCase.verifyFalse(ffe.decide); testCase.verifyFalse(ffe.decide);
testCase.verifyEqual(ffe.e, zeros(ffe.order, 1)); testCase.verifyEqual(ffe.e, zeros(ffe.order, 1));
testCase.verifyEqual(ffe.error, 0); testCase.verifyEqual(ffe.error, 0);
@@ -69,6 +86,468 @@ classdef FFE_test < IMDDTestCase
testCase.verifyNotEqual(y.signal, zeros(size(x.signal))); testCase.verifyNotEqual(y.signal, zeros(size(x.signal)));
testCase.verifyEqual(length(y.signal), length(x.signal)); testCase.verifyEqual(length(y.signal), length(x.signal));
end end
function dcTrackingBufferDelaysDcUpdateUntilBlockBoundary(testCase)
x = zeros(4, 1);
d = ones(4, 1);
ffe = FFE( ...
"sps", 1, ...
"order", 1, ...
"dc_tracking_mu", 0.1, ...
"dc_tracking_buffer_len", 2, ...
"dd_mode", 0, ...
"adaption_technique", adaption_method.lms);
ffe.constellation = unique(d);
ffe.equalize(x, d, 0, 1, numel(x), true, false);
testCase.verifyEqual(ffe.e_dc, 0.19, "AbsTol", 1e-12);
end
function dcTrackingUnbufferedPathUpdatesEverySymbol(testCase)
x = zeros(4, 1);
d = ones(4, 1);
ffe = FFE( ...
"sps", 1, ...
"order", 1, ...
"dc_tracking_mu", 0.1, ...
"dc_tracking_buffer_len", 1, ...
"dd_mode", 0, ...
"adaption_technique", adaption_method.lms);
ffe.constellation = unique(d);
ffe.equalize(x, d, 0, 1, numel(x), true, false);
testCase.verifyEqual(ffe.e_dc, 1 - 0.9^4, "AbsTol", 1e-12);
end
function dcTrackingBlockUpdateUsesMeanBlockError(testCase)
ffe = FFE("dc_tracking_mu", 0.1);
[e_dc_next,stats] = ffe.dcTrackingBlockUpdate(0.5, [1; 2; 3]);
testCase.verifyEqual(e_dc_next, 0.7, "AbsTol", 1e-12);
testCase.verifyEqual(stats.err_mean, 2, "AbsTol", 1e-12);
testCase.verifyEqual(stats.mu_eff, 0.1, "AbsTol", 1e-12);
testCase.verifyEqual(stats.update, 0.2, "AbsTol", 1e-12);
end
function dcTrackingBlockUpdateDoesNotSuppressLargeErrorPower(testCase)
ffe = FFE("dc_tracking_mu", 0.1);
[e_dc_next,stats] = ffe.dcTrackingBlockUpdate(0, [10; 10]);
testCase.verifyEqual(e_dc_next, 1, "AbsTol", 1e-12);
testCase.verifyEqual(stats.err_mean, 10, "AbsTol", 1e-12);
testCase.verifyEqual(stats.mu_eff, 0.1, "AbsTol", 1e-12);
end
function dcTrackingBlockUpdateCanBoostPersistentMean(testCase)
ffe = FFE("dc_tracking_mu", 0.1);
[e_dc_next,stats] = ffe.dcTrackingBlockUpdate(0, [2; 2; 2], ...
"persistence_gain", 1);
testCase.verifyEqual(e_dc_next, 0.4, "AbsTol", 1e-12);
testCase.verifyEqual(stats.persistence_scale, 1, "AbsTol", 1e-12);
testCase.verifyEqual(stats.mu_eff, 0.2, "AbsTol", 1e-12);
end
function dcTrackingLoopUsesBlockPersistenceGain(testCase)
x = zeros(4, 1);
d = ones(4, 1);
ffe = FFE( ...
"sps", 1, ...
"order", 1, ...
"dc_tracking_mu", 0.1, ...
"dc_tracking_adaptive_enabled", true, ...
"dc_tracking_persistence_gain", 1, ...
"dc_tracking_buffer_len", 2, ...
"dd_mode", 0, ...
"adaption_technique", adaption_method.lms);
ffe.constellation = unique(d);
ffe.equalize(x, d, 0, 1, numel(x), true, false);
testCase.verifyEqual(ffe.e_dc, 0.36, "AbsTol", 1e-12);
end
function 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);
ffe = FFE( ...
"sps", 1, ...
"order", 1, ...
"dc_tracking_mu", 0, ...
"ffe_update_buffer_len", 2, ...
"dd_mode", 0, ...
"adaption_technique", adaption_method.lms);
ffe.constellation = unique(d);
ffe.equalize(x, d, 0.1, 1, numel(x), true, false);
testCase.verifyEqual(ffe.e, 0.19, "AbsTol", 1e-12);
end
function dcAvgA1BufferlengthOneLeavesInputUnchanged(testCase)
x = (1:4).';
d = zeros(4, 1);
ffe = FFE( ...
"sps", 1, ...
"order", 1, ...
"dc_avg_bufferlength_a1", 1, ...
"save_debug", true, ...
"dd_mode", 0, ...
"adaption_technique", adaption_method.lms);
ffe.e = 1;
[y,~] = ffe.equalize(x, d, 0, 1, numel(x), true, false);
testCase.verifyEqual(y, x, "AbsTol", 1e-12);
testCase.verifyEqual(ffe.debug_struct.dc_avg_offset, zeros(1,4), "AbsTol", 1e-12);
end
function dcAvgA1SubtractsCausalMovingAverageByDefault(testCase)
x = (1:4).';
d = zeros(4, 1);
expected_offset = [0, 0, 0, 2];
expected_y = [1; 2; 3; 2];
ffe = FFE( ...
"sps", 1, ...
"order", 1, ...
"dc_avg_bufferlength_a1", 3, ...
"save_debug", true, ...
"dd_mode", 0, ...
"adaption_technique", adaption_method.lms);
ffe.e = 1;
[y,~] = ffe.equalize(x, d, 0, 1, numel(x), true, false);
testCase.verifyEqual(y, expected_y, "AbsTol", 1e-12);
testCase.verifyEqual(ffe.debug_struct.dc_avg_offset, expected_offset, "AbsTol", 1e-12);
end
function dcAvgA1SmoothingIsReservedForOfflineVariants(testCase)
x = (1:4).';
d = zeros(4, 1);
expected_offset = [0, 0, 0, 2];
expected_y = [1; 2; 3; 2];
ffe = FFE( ...
"sps", 1, ...
"order", 1, ...
"dc_avg_bufferlength_a1", 3, ...
"dc_smoothing_a1", 0.5, ...
"save_debug", true, ...
"dd_mode", 0, ...
"adaption_technique", adaption_method.lms);
ffe.e = 1;
[y,~] = ffe.equalize(x, d, 0, 1, numel(x), true, false);
testCase.verifyEqual(y, expected_y, "AbsTol", 1e-12);
testCase.verifyEqual(ffe.debug_struct.dc_avg_offset, expected_offset, "AbsTol", 1e-12);
end
function dcAvgA1UpdateBlocklengthOneRefreshesEverySymbol(testCase)
x = (1:4).';
d = zeros(4, 1);
expected_offset = [0, 1, 1.5, 2];
expected_y = [1; 1; 1.5; 2];
ffe = FFE( ...
"sps", 1, ...
"order", 1, ...
"dc_avg_bufferlength_a1", 3, ...
"dc_avg_update_blocklength_a1", 1, ...
"save_debug", true, ...
"dd_mode", 0, ...
"adaption_technique", adaption_method.lms);
ffe.e = 1;
[y,~] = ffe.equalize(x, d, 0, 1, numel(x), true, false);
testCase.verifyEqual(y, expected_y, "AbsTol", 1e-12);
testCase.verifyEqual(ffe.debug_struct.dc_avg_offset, expected_offset, "AbsTol", 1e-12);
end
function dcAvgA1WarnsWhenCombinedWithDcTracking(testCase)
testCase.verifyWarning(@() FFE( ...
"dc_avg_bufferlength_a1", 3, ...
"dc_tracking_mu", 0.1), "FFE:DCAvgWithDcTracking");
end
function dcLevelAvgA2SubtractsScalarWeightedResidual(testCase)
x = (1:4).';
d = zeros(4, 1);
expected_y = [1; 2; 1.5; 2.5];
ffe = FFE( ...
"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, expected_y, "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 dcLevelAvgA2UpdateBlocklengthOneRefreshesEverySymbol(testCase)
x = (1:4).';
d = zeros(4, 1);
ffe = FFE( ...
"sps", 1, ...
"order", 1, ...
"dc_level_avg_bufferlength_a2", 2, ...
"dc_level_update_blocklength_a2", 1, ...
"dc_level_weights_a2", 1, ...
"save_debug", true, ...
"dd_mode", 0, ...
"adaption_technique", adaption_method.lms);
ffe.e = 1;
ffe.constellation = 0;
[y,~] = ffe.equalize(x, d, 0, 1, numel(x), true, false);
testCase.verifyEqual(y, [1; 1.5; 1.5; 1.5], "AbsTol", 1e-12);
testCase.verifyEqual(ffe.debug_struct.dc_level_mpi_est, [0, 1, 1.5, 2.5], "AbsTol", 1e-12);
testCase.verifyEqual(ffe.debug_struct.dc_level_weight, [0, 0.5, 1, 1], "AbsTol", 1e-12);
testCase.verifyEqual(ffe.debug_struct.dc_level_valid_count, [0, 1, 2, 2]);
end
function dcLevelAvgA2UsesVectorWeightsInConstellationOrder(testCase)
x = [1; 3; 1; 4];
d = [0; 2; 0; 2];
ffe = FFE( ...
"sps", 1, ...
"order", 1, ...
"dc_level_avg_bufferlength_a2", 2, ...
"dc_level_weights_a2", [0, 1], ...
"dd_mode", 0, ...
"adaption_technique", adaption_method.lms);
ffe.e = 1;
ffe.constellation = [0; 2];
[y,~] = ffe.equalize(x, d, 0, 1, numel(x), true, false);
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, ...
"dc_level_weights_a2", 1, ...
"dc_tracking_mu", 0.1), "FFE:DCLevelAvgWithOtherDcSuppression");
end
function a2InitialWeightGuessFollowsLevelVarianceSlope(testCase)
x = [ ...
0.98; 1.02; 1.00; 1.01; ...
1.90; 2.10; 2.00; 2.08; ...
2.50; 3.50; 3.00; 3.40];
d = [ ...
ones(4,1); ...
2*ones(4,1); ...
3*ones(4,1)];
ffe = FFE("sps", 1, "a2_level_weight_max", 1);
ffe.constellation = [1; 2; 3];
[weights,stats] = ffe.a2LevelWeightInitialGuess(x,d);
testCase.verifySize(weights, [3, 1]);
testCase.verifyGreaterThan(weights(3), weights(2));
testCase.verifyGreaterThan(weights(2), weights(1));
testCase.verifyEqual(max(weights), 0.7, "AbsTol", 1e-12);
testCase.verifyGreaterThan(stats.variance_slope, 0);
end
end end
end end

BIN
before_eq_-1.png Normal file

Binary file not shown.

After

Width:  |  Height:  |  Size: 250 KiB

View File

@@ -0,0 +1,243 @@
%% Simple BER over SIR plot from the recombined warehouses
% This script intentionally keeps the mechanics visible:
% 1) load the run_id warehouse and the grouped config warehouse
% 2) loop over occupied config cells
% 3) read the physical coordinates of the cell
% 4) extract all BER values stored in that cell
% 5) plot individual BER dots and one larger mean marker per SIR
clear; clc;
scriptDir = fileparts("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");
runWhData = load(runWhFile);
configWhData = load(configWhFile);
wh_run_id_combined = runWhData.wh_run_id_combined;
wh_config_combined = configWhData.wh_config_combined;
algorithmStorageNames = configWhData.algorithmStorageNames;
fprintf("Loaded run_id warehouse:\n %s\n", runWhFile);
wh_run_id_combined.showInfo;
fprintf("\nLoaded config-grouped warehouse:\n %s\n", configWhFile);
wh_config_combined.showInfo;
%% Plot settings
pathLengthToPlot = 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
polyfitOrderMax = 4;
% SIR vs BER: one subplot per PAM level, all algorithms compared
if exist("linspecer", "file")
algColors = linspecer(numel(selectedAlgorithms));
else
algColors = lines(numel(selectedAlgorithms));
end
figure(); clf;
tiledlayout(numel(selectedPamLevels), 1, "TileSpacing", "compact");
for pamIdx = 1:numel(selectedPamLevels)
pamLevel = selectedPamLevels(pamIdx);
nexttile; hold on;
for algIdx = 1:numel(selectedAlgorithms)
algorithmName = selectedAlgorithms{algIdx};
algColor = algColors(algIdx, :);
berTable = buildBerTable(wh_config_combined, algorithmName, pathLengthToPlot);
berTable = berTable(berTable.pam_level == pamLevel, :);
% berTable(berTable.sir == 23,:) = [];
if isempty(berTable)
continue
end
sirValues = unique(berTable.sir(:).');
sirValues = sort(sirValues);
meanBer = nan(size(sirValues));
minBer = nan(size(sirValues));
maxBer = nan(size(sirValues));
for sirIdx = 1:numel(sirValues)
sirValue = sirValues(sirIdx);
sirRows = berTable(berTable.sir == sirValue, :);
berValues = [sirRows.ber_values{:}];
berValues = berValues(isfinite(berValues));
berValues = berValues(berValues<0.1);
berValues = rmoutliers(berValues);
if isempty(berValues)
continue
end
scatter(repmat(sirValue, size(berValues)), berValues, ...
30, ...
"Marker", ".", ...
"MarkerEdgeColor", algColor, ...
"MarkerFaceColor", algColor, ...
"HandleVisibility", "off");
meanBer(sirIdx) = mean(berValues, "omitnan");
minBer(sirIdx) = min(berValues);
maxBer(sirIdx) = max(berValues);
end
valid = isfinite(sirValues) & isfinite(meanBer);
if ~any(valid)
continue
end
if useBoundedLines && exist("boundedline", "file")
yLower = max(meanBer - minBer, 0);
yUpper = max(maxBer - meanBer, 0);
yBounds = [yLower(:), yUpper(:)];
[hl, hp] = boundedline(sirValues(valid).', meanBer(valid).', yBounds(valid, :), ...
'alpha', 'transparency', 0.08, ...
'cmap', algColor, ...
'nan', 'fill', ...
'orientation', 'vert');
set(hl, "LineStyle", "none", "Marker", "none", "HandleVisibility", "off");
set(hp, "LineStyle", "none", "HandleVisibility", "off");
end
scatter(sirValues(valid), meanBer(valid), ...
20, ...
"Marker", "o", ...
"MarkerEdgeColor", algColor, ...
"MarkerFaceColor", algColor, ...
"LineWidth", 1, ...
"DisplayName", algorithmName);
if usePolyfit
fitMask = valid & meanBer > 0;
if nnz(fitMask) >= 2
fitOrder = min(polyfitOrderMax, nnz(fitMask) - 1);
fitCoeff = polyfit(sirValues(fitMask), log10(meanBer(fitMask)), fitOrder);
xFit = linspace(min(sirValues(fitMask)), max(sirValues(fitMask)), 300);
yFit = 10 .^ polyval(fitCoeff, xFit);
plot(xFit, yFit, ...
"LineStyle", "--", ...
"LineWidth", 1.1, ...
"Color", algColor, ...
"HandleVisibility", "off");
end
end
end
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");
set(gca, "YScale", "log");
xlim([15,45]);
grid on;
box on;
legend("Location", "best", "Interpreter", "none");
end
%% Local helpers
function berTable = buildBerTable(wh, storageName, pathLength)
rows = struct( ...
"pam_level", {}, ...
"symbolrate", {}, ...
"interference_path_length", {}, ...
"sir", {}, ...
"block_update", {}, ...
"run_ids", {}, ...
"n_run_ids", {}, ...
"ber_values", {}, ...
"n_ber_values", {});
for linIdx = 1:numel(wh.sto.(storageName))
storedValue = wh.sto.(storageName){linIdx};
if isempty(storedValue)
continue
end
phys = linearIndexToStruct(wh, linIdx);
if phys.interference_path_length ~= pathLength
continue
end
berValues = extractBerValues(storedValue);
if isempty(berValues)
continue
end
runIds = wh.sto.run_id{linIdx};
if isempty(runIds)
runIds = NaN;
end
row = struct();
row.pam_level = phys.pam_level;
row.symbolrate = phys.symbolrate;
row.interference_path_length = phys.interference_path_length;
row.sir = phys.sir;
row.block_update = phys.block_update;
row.run_ids = {runIds};
row.n_run_ids = numel(runIds);
row.ber_values = {berValues};
row.n_ber_values = numel(berValues);
rows(end+1) = row; %#ok<AGROW>
end
if isempty(rows)
berTable = struct2table(rows);
else
berTable = struct2table(rows);
berTable = sortrows(berTable, ["pam_level", "sir", "symbolrate"]);
end
end
function physStruct = linearIndexToStruct(wh, linIdx)
[physValues, physNames] = wh.getPhysIndicesByLinIndex(linIdx);
physStruct = struct();
for paramIdx = 1:numel(physNames)
physStruct.(char(physNames{paramIdx})) = physValues{paramIdx};
end
end
function berValues = extractBerValues(value)
berValues = [];
if isstruct(value)
if isfield(value, "metrics") && hasBer(value.metrics)
berValues(end+1) = value.metrics.BER;
end
elseif iscell(value)
for valueIdx = 1:numel(value)
berValues = [berValues, extractBerValues(value{valueIdx})]; %#ok<AGROW>
end
end
end
function tf = hasBer(metrics)
tf = (isstruct(metrics) && isfield(metrics, "BER")) || ...
(isobject(metrics) && isprop(metrics, "BER"));
end

View File

@@ -0,0 +1,254 @@
%% BER over SIR from MpiReductionResults
% 1) gather data from the database
% 2) clean data and remove per-curve outliers
% 3) plot BER over SIR grouped by algorithm
clear; clc;
%% 1) Gather data
studyName = "pam4_112_greater_3153";
pathLengthToPlot = 0;
selectedBlockUpdate = 1; % set [] to pool all block_update values
selectedPamLevels = []; % set [] to use all PAM levels in the query result
selectedAlgorithms = []; % set [] to use all algorithms in the query result
useBoundedLines = true;
usePolyfit = true;
polyfitOrderMax = 4;
maxBerForPlot = 0.1;
db = DBHandler( ...
"dataBase", "labor", ...
"type", "mysql", ...
"server", "192.168.178.192", ...
"user", "silas", ...
"password", "silas");
db.refresh();
fp = QueryFilter();
fp.where('Runs', 'interference_path_length', 'EQUALS', pathLengthToPlot);
fp.where('Runs', 'fiber_length', 'EQUALS', 0);
fp.where('Runs', 'db_mode', 'EQUALS', '"no_db"');
fp.where('MpiReductionResults', 'study_name', 'EQUALS', char(studyName));
if ~isempty(selectedBlockUpdate)
fp.where('MpiReductionResults', 'block_update', 'EQUALS', selectedBlockUpdate);
end
selectedFields = { ...
'Runs.run_id', ...
'Runs.sir', ...
'Runs.pam_level', ...
'Runs.symbolrate', ...
'Runs.interference_path_length', ...
'Runs.power_mpi_interference', ...
'MpiReductionResults.occurrence_idx', ...
'MpiReductionResults.storage_name', ...
'MpiReductionResults.algorithm', ...
'MpiReductionResults.algorithm_variant', ...
'MpiReductionResults.eq_class', ...
'MpiReductionResults.block_update', ...
'MpiReductionResults.BER', ...
'MpiReductionResults.BER_precoded', ...
'MpiReductionResults.SNR', ...
'MpiReductionResults.GMI', ...
'MpiReductionResults.AIR'};
selectedFields = selectedFields(:);
[rawData, query] = db.queryDB(fp, selectedFields);
disp(query);
fprintf("Fetched %d MPI reduction rows.\n", height(rawData));
%% 2) Clean data
data = rawData;
numericFields = ["run_id", "sir", "pam_level", "symbolrate", ...
"interference_path_length", "occurrence_idx", "block_update", ...
"power_mpi_interference", "BER", "BER_precoded", "SNR", "GMI", "AIR"];
for fieldIdx = 1:numel(numericFields)
fieldName = numericFields(fieldIdx);
if ismember(fieldName, string(data.Properties.VariableNames))
data.(char(fieldName)) = numericColumn(data.(char(fieldName)));
end
end
stringFields = ["storage_name", "algorithm", "algorithm_variant", "eq_class"];
for fieldIdx = 1:numel(stringFields)
fieldName = stringFields(fieldIdx);
if ismember(fieldName, string(data.Properties.VariableNames))
data.(char(fieldName)) = stringColumn(data.(char(fieldName)));
end
end
data = data(isfinite(data.BER) & data.BER > 0 & data.BER < maxBerForPlot, :);
data.sir_exact = -7 - data.power_mpi_interference;
if isempty(data)
warning("PLOT_mpi_reduction_db_ber_vs_sir:NoRows", ...
"No rows remain after initial BER/path/study/block filtering.");
return
end
if isempty(selectedPamLevels)
selectedPamLevels = unique(data.pam_level(isfinite(data.pam_level))).';
else
data = data(ismember(data.pam_level, selectedPamLevels), :);
end
if isempty(selectedAlgorithms)
selectedAlgorithms = unique(data.algorithm, "stable").';
else
selectedAlgorithms = string(selectedAlgorithms);
data = data(ismember(data.algorithm, selectedAlgorithms), :);
end
if isempty(data)
warning("PLOT_mpi_reduction_db_ber_vs_sir:NoSelectedRows", ...
"No rows remain after selectedPamLevels/selectedAlgorithms filtering.");
return
end
data.clean_keep = true(height(data), 1);
[groupId, groupPam, groupAlgorithm, groupSir] = findgroups(data.pam_level, data.algorithm, data.sir);
for groupIdx = 1:max(groupId)
rowMask = groupId == groupIdx;
berValues = data.BER(rowMask);
if nnz(rowMask) > 3
data.clean_keep(rowMask) = ~isoutlier(berValues);
end
end
cleanData = data(data.clean_keep, :);
summaryTable = groupsummary(cleanData, ["pam_level", "algorithm", "sir"], ...
{"mean", "min", "max"}, "BER");
sirExactTable = groupsummary(cleanData, ["pam_level", "algorithm", "sir"], ...
"median", "sir_exact");
summaryTable.sir_exact = sirExactTable.median_sir_exact;
summaryTable = sortrows(summaryTable, ["pam_level", "algorithm", "sir"]);
fprintf("Cleaned to %d rows across %d PAM/algorithm/SIR groups.\n", ...
height(cleanData), height(summaryTable));
disp(groupcounts(cleanData, ["pam_level", "algorithm"]));
%% 3) Plot data
if exist("linspecer", "file")
algColors = linspecer(numel(selectedAlgorithms));
else
algColors = lines(numel(selectedAlgorithms));
end
figure(); clf;
tiledlayout(numel(selectedPamLevels), 1, "TileSpacing", "compact");
for pamIdx = 1:numel(selectedPamLevels)
pamLevel = selectedPamLevels(pamIdx);
nexttile; hold on;
for algIdx = 1:numel(selectedAlgorithms)
algorithmName = selectedAlgorithms(algIdx);
algColor = algColors(algIdx, :);
rawMask = cleanData.pam_level == pamLevel & cleanData.algorithm == algorithmName;
curveMask = summaryTable.pam_level == pamLevel & summaryTable.algorithm == algorithmName;
if ~any(curveMask)
continue
end
scatter(cleanData.sir_exact(rawMask), cleanData.BER(rawMask), ...
30, ...
"Marker", ".", ...
"MarkerEdgeColor", algColor, ...
"HandleVisibility", "off");
sirValues = summaryTable.sir_exact(curveMask).';
meanBer = summaryTable.mean_BER(curveMask).';
minBer = summaryTable.min_BER(curveMask).';
maxBer = summaryTable.max_BER(curveMask).';
valid = isfinite(sirValues) & isfinite(meanBer) & meanBer > 0;
if useBoundedLines && exist("boundedline", "file") && any(valid)
yLower = max(meanBer - minBer, 0);
yUpper = max(maxBer - meanBer, 0);
yBounds = [yLower(:), yUpper(:)];
[hl, hp] = boundedline(sirValues(valid).', meanBer(valid).', yBounds(valid, :), ...
'alpha', 'transparency', 0.08, ...
'cmap', algColor, ...
'nan', 'fill', ...
'orientation', 'vert');
set(hl, "LineStyle", "none", "Marker", "none", "HandleVisibility", "off");
set(hp, "LineStyle", "none", "HandleVisibility", "off");
end
plot(sirValues(valid), meanBer(valid), ...
"LineStyle", "none", ...
"Marker", "o", ...
"MarkerSize", 5, ...
"LineWidth", 1.2, ...
"Color", algColor, ...
"MarkerFaceColor", algColor, ...
"DisplayName", char(algorithmName));
if usePolyfit
fitMask = valid & meanBer > 0;
if nnz(fitMask) >= 2
fitOrder = min(polyfitOrderMax, nnz(fitMask) - 1);
fitCoeff = polyfit(sirValues(fitMask), log10(meanBer(fitMask)), fitOrder);
xFit = linspace(min(sirValues(fitMask)), max(sirValues(fitMask)), 300);
yFit = 10 .^ polyval(fitCoeff, xFit);
plot(xFit, yFit, ...
"LineStyle", "--", ...
"LineWidth", 1.1, ...
"Color", algColor, ...
"HandleVisibility", "off");
end
end
end
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");
title(sprintf("PAM %.0f, path %.0f m, block update %s", ...
pamLevel, pathLengthToPlot, blockUpdateLabel(selectedBlockUpdate)));
xlabel("SIR (dB)");
ylabel("BER");
set(gca, "YScale", "log");
ylim([9e-5, maxBerForPlot]);
grid on;
box on;
legend("Location", "best", "Interpreter", "none");
end
if exist("beautifyBERplot", "file")
beautifyBERplot("logscale", true, "setcolors", false, "setmarkers", false);
end
%% Local helpers
function values = numericColumn(values)
if iscell(values)
values = string(values);
end
if isstring(values) || ischar(values)
values = str2double(values);
end
values = double(values);
end
function values = stringColumn(values)
if iscell(values)
values = string(values);
elseif ischar(values)
values = string(values);
end
values = strip(string(values));
end
function label = blockUpdateLabel(selectedBlockUpdate)
if isempty(selectedBlockUpdate)
label = "pooled";
else
label = string(selectedBlockUpdate);
end
end

View File

@@ -0,0 +1,539 @@
% === DSP settings ===
dsp_options = struct();
dsp_options.mode = "run_id";
dsp_options.recipe = @mpi_recipe_dev;
dsp_options.append_to_db = false;
dsp_options.start_occurence = 1;
dsp_options.max_occurences = 15;
dsp_options.debug_plots = false;
dsp_options.database_type = "mysql";
dsp_options.dataBase = "labor";
dsp_options.storage_path = "W:\labdata\ECOC Silas\ecoc_2025";
dsp_options.server = "192.168.178.192";
dsp_options.port = 3306;
dsp_options.user = "silas";
dsp_options.password = "silas";
db = DBHandler("dataBase", [dsp_options.dataBase], ...
"type", dsp_options.database_type, ...
"server", dsp_options.server, ...
"user", dsp_options.user, "password", dsp_options.password);
scriptDir = fileparts(mfilename("fullpath"));
whpam4_loaded = load(fullfile(scriptDir, "pam_4_results.mat"));
whpam6_loaded = load(fullfile(scriptDir, "pam_6_results.mat"));
whpam8_loaded = load(fullfile(scriptDir, "pam_8_results.mat"));
whpam4 = whpam4_loaded.obj;
whpam6 = whpam6_loaded.obj;
whpam8 = whpam8_loaded.obj;
whList = {whpam4, whpam6, whpam8};
algorithmStorageNames = fieldnames(whpam4.sto);
%% Build a lossless run_id warehouse and cache run metadata
pamformats = [4, 6, 8];
baudrates = [112e9, 96e9, 72e9];
runMetaTable = queryRunMetadata(db, whList, pamformats, baudrates);
wh_run_id_combined = mergeDataStoragesByUnion(whList);
wh_run_id_combined = addRunMetadataStorages(wh_run_id_combined, runMetaTable);
runIdSavePath = fullfile(scriptDir, "combined_by_run_id_results.mat");
save(runIdSavePath, "wh_run_id_combined", "runMetaTable", "algorithmStorageNames", '-v7.3');
fprintf("Saved run_id-combined warehouse:\n %s\n", runIdSavePath);
wh_run_id_combined.showInfo;
%% Build a grouped config warehouse for analysis views
configFields = ["pam_level", "symbolrate", "interference_path_length", "sir"];
wh_config_combined = buildConfigWarehouse( ...
wh_run_id_combined, runMetaTable, algorithmStorageNames, configFields);
configSavePath = fullfile(scriptDir, "combined_by_config_results.mat");
save(configSavePath, "wh_config_combined", "runMetaTable", "algorithmStorageNames", '-v7.3');
fprintf("Saved config-grouped warehouse:\n %s\n", configSavePath);
wh_config_combined.showInfo;
%% BER over SIR at fixed path length
plot_options = struct();
plot_options.path_length = 1000;
plot_options.use_boundedlines = true;
plot_options.boundedline_alpha = 0.10;
plot_options.polyfit_order_max = 3;
plot_options.fec_ber_threshold = 2e-2;
plot_options.figure_base = 7100;
plot_options.xlim = [15 35];
plotBerOverSirForPath(wh_config_combined, algorithmStorageNames, plot_options);
%% Local functions
function runMetaTable = queryRunMetadata(db, whList, pamformats, baudrates)
metaTables = {};
metadataFields = db.getTableFieldNames('Runs');
for whIdx = 1:numel(whList)
curWh = whList{whIdx};
curRunIds = curWh.parameter.run_id.values(:);
fp = QueryFilter();
fp.where('Runs', 'run_id', 'GREATER_EQUAL', 3153);
fp.where('Runs', 'symbolrate', 'EQUALS', baudrates(whIdx));
fp.where('Runs', 'fiber_length', 'EQUALS', 0);
fp.where('Runs', 'db_mode', 'EQUALS', '"no_db"');
fp.where('Runs', 'pam_level', 'EQUALS', pamformats(whIdx));
[dataTable, ~] = db.queryDB(fp, metadataFields);
dataTable.run_id = numericColumn(dataTable.run_id);
dataTable = dataTable(ismember(dataTable.run_id, curRunIds), :);
missingRunIds = setdiff(curRunIds, dataTable.run_id);
if ~isempty(missingRunIds)
warning("recombine_warehouses:MissingMetadata", ...
"%d run_id(s) from PAM %.0f warehouse were not found in the metadata query.", ...
numel(missingRunIds), pamformats(whIdx));
end
metaTables{end+1} = dataTable; %#ok<AGROW>
end
runMetaTable = vertcat(metaTables{:});
[~, uniqueIdx] = unique(numericColumn(runMetaTable.run_id), "stable");
runMetaTable = runMetaTable(uniqueIdx, :);
numericFields = ["run_id", "pam_level", "symbolrate", "interference_path_length", "sir"];
for fieldIdx = 1:numel(numericFields)
fieldName = numericFields(fieldIdx);
if ismember(fieldName, string(runMetaTable.Properties.VariableNames))
runMetaTable.(char(fieldName)) = numericColumn(runMetaTable.(char(fieldName)));
end
end
[~, sortIdx] = sort(runMetaTable.run_id);
runMetaTable = runMetaTable(sortIdx, :);
end
function x = numericColumn(x)
if isnumeric(x)
x = double(x);
elseif iscell(x)
x = str2double(string(x));
elseif isstring(x) || ischar(x) || iscategorical(x)
x = str2double(string(x));
else
x = double(x);
end
x = x(:);
end
function whMerged = mergeDataStoragesByUnion(whList)
templateWh = whList{1};
paramNames = cellstr(templateWh.fn);
for whIdx = 2:numel(whList)
if ~isequal(sort(cellstr(whList{whIdx}.fn)), sort(paramNames))
error("recombine_warehouses:ParameterMismatch", ...
"All input warehouses must use the same parameter names.");
end
end
mergedParams = struct();
for paramIdx = 1:numel(paramNames)
paramName = paramNames{paramIdx};
values = [];
for whIdx = 1:numel(whList)
newValues = whList{whIdx}.parameter.(paramName).values(:).';
values = [values, newValues]; %#ok<AGROW>
end
values = unique(values, "stable");
if isnumeric(values)
values = sort(values);
end
mergedParams.(paramName) = values;
end
whMerged = DataStorage(mergedParams);
for whIdx = 1:numel(whList)
curWh = whList{whIdx};
curStorageNames = fieldnames(curWh.sto);
for storageIdx = 1:numel(curStorageNames)
storageName = curStorageNames{storageIdx};
if ~isfield(whMerged.sto, storageName)
whMerged.addStorage(storageName);
end
for linIdx = 1:numel(curWh.sto.(storageName))
storedValue = curWh.sto.(storageName){linIdx};
if isempty(storedValue)
continue
end
targetLinIdx = mapLinearIndex(curWh, whMerged, linIdx);
existingValue = whMerged.sto.(storageName){targetLinIdx};
whMerged.sto.(storageName){targetLinIdx} = mergeStoredValue(existingValue, storedValue);
end
end
end
end
function wh = addRunMetadataStorages(wh, runMetaTable)
metadataStorageNames = ["meta_pam_level", "meta_symbolrate", ...
"meta_interference_path_length", "meta_sir"];
metadataTableFields = ["pam_level", "symbolrate", "interference_path_length", "sir"];
for fieldIdx = 1:numel(metadataStorageNames)
if ~isfield(wh.sto, metadataStorageNames(fieldIdx))
wh.addStorage(char(metadataStorageNames(fieldIdx)));
end
end
for linIdx = 1:wh.getLastLinIndice()
physStruct = linearIndexToStruct(wh, linIdx);
if ~isfield(physStruct, "run_id")
continue
end
metaIdx = find(runMetaTable.run_id == physStruct.run_id, 1);
if isempty(metaIdx)
continue
end
for fieldIdx = 1:numel(metadataStorageNames)
wh.sto.(char(metadataStorageNames(fieldIdx))){linIdx} = ...
runMetaTable.(char(metadataTableFields(fieldIdx)))(metaIdx);
end
end
end
function whConfig = buildConfigWarehouse(whRun, runMetaTable, algorithmStorageNames, configFields)
runParamNames = cellstr(whRun.fn);
extraParamNames = setdiff(string(runParamNames), "run_id", "stable");
configParams = struct();
for fieldIdx = 1:numel(configFields)
fieldName = configFields(fieldIdx);
values = unique(runMetaTable.(char(fieldName))(:).', "stable");
if isnumeric(values)
values = sort(values);
end
configParams.(char(fieldName)) = values;
end
for paramIdx = 1:numel(extraParamNames)
paramName = extraParamNames(paramIdx);
configParams.(char(paramName)) = whRun.parameter.(char(paramName)).values;
end
whConfig = DataStorage(configParams);
for storageIdx = 1:numel(algorithmStorageNames)
whConfig.addStorage(algorithmStorageNames{storageIdx});
end
whConfig.addStorage("run_id");
for storageIdx = 1:numel(algorithmStorageNames)
storageName = algorithmStorageNames{storageIdx};
for linIdx = 1:numel(whRun.sto.(storageName))
storedValue = whRun.sto.(storageName){linIdx};
if isempty(storedValue)
continue
end
sourcePhys = linearIndexToStruct(whRun, linIdx);
metaIdx = find(runMetaTable.run_id == sourcePhys.run_id, 1);
if isempty(metaIdx)
continue
end
targetArgs = configTargetArgs(whConfig, sourcePhys, runMetaTable, metaIdx, configFields);
targetLinIdx = whConfig.getIndicesByPhys(targetArgs);
existingValue = whConfig.sto.(storageName){targetLinIdx};
whConfig.sto.(storageName){targetLinIdx} = mergeStoredValue(existingValue, storedValue);
end
end
for linIdx = 1:whRun.getLastLinIndice()
sourcePhys = linearIndexToStruct(whRun, linIdx);
metaIdx = find(runMetaTable.run_id == sourcePhys.run_id, 1);
if isempty(metaIdx)
continue
end
targetArgs = configTargetArgs(whConfig, sourcePhys, runMetaTable, metaIdx, configFields);
targetLinIdx = whConfig.getIndicesByPhys(targetArgs);
whConfig.sto.run_id{targetLinIdx} = mergeUniqueNumeric( ...
whConfig.sto.run_id{targetLinIdx}, sourcePhys.run_id);
end
end
function targetArgs = configTargetArgs(whConfig, sourcePhys, runMetaTable, metaIdx, configFields)
targetArgs = cell(1, numel(whConfig.fn));
for paramIdx = 1:numel(whConfig.fn)
paramName = string(whConfig.fn(paramIdx));
if any(configFields == paramName)
targetArgs{paramIdx} = runMetaTable.(char(paramName))(metaIdx);
else
targetArgs{paramIdx} = sourcePhys.(char(paramName));
end
end
end
function targetLinIdx = mapLinearIndex(sourceWh, targetWh, sourceLinIdx)
sourcePhys = linearIndexToStruct(sourceWh, sourceLinIdx);
targetArgs = cell(1, numel(targetWh.fn));
for paramIdx = 1:numel(targetWh.fn)
paramName = char(targetWh.fn(paramIdx));
targetArgs{paramIdx} = sourcePhys.(paramName);
end
targetLinIdx = targetWh.getIndicesByPhys(targetArgs);
end
function physStruct = linearIndexToStruct(wh, linIdx)
[physValues, physNames] = wh.getPhysIndicesByLinIndex(linIdx);
physStruct = struct();
for paramIdx = 1:numel(physNames)
physStruct.(char(physNames{paramIdx})) = physValues{paramIdx};
end
end
function out = mergeStoredValue(existingValue, newValue)
if isempty(existingValue)
out = newValue;
return
end
existingPackages = asPackageCell(existingValue);
newPackages = asPackageCell(newValue);
out = [existingPackages(:); newPackages(:)].';
end
function packages = asPackageCell(value)
if isempty(value)
packages = {};
elseif iscell(value)
packages = value(:).';
else
packages = {value};
end
end
function out = mergeUniqueNumeric(existingValue, newValue)
if isempty(existingValue)
out = newValue;
else
out = unique([existingValue(:); newValue(:)].', "stable");
end
end
function plotBerOverSirForPath(whConfig, algorithmStorageNames, plotOptions)
paramNames = string(whConfig.fn);
if ~any(paramNames == "pam_level") || ~any(paramNames == "interference_path_length") || ~any(paramNames == "sir")
error("recombine_warehouses:MissingPlotParameters", ...
"The config warehouse must contain pam_level, interference_path_length, and sir parameters.");
end
pamLevels = whConfig.parameter.pam_level.values;
hasBlockUpdate = any(paramNames == "block_update");
if hasBlockUpdate
blockUpdates = whConfig.parameter.block_update.values;
else
blockUpdates = NaN;
end
if exist("linspecer", "file")
curveColors = linspecer(max(numel(algorithmStorageNames), 1));
else
curveColors = lines(max(numel(algorithmStorageNames), 1));
end
figureIdx = plotOptions.figure_base;
for blockIdx = 1:numel(blockUpdates)
for pamIdx = 1:numel(pamLevels)
pamLevel = pamLevels(pamIdx);
figureIdx = figureIdx + 1;
figure(figureIdx); clf; hold on
for storageIdx = 1:numel(algorithmStorageNames)
storageName = algorithmStorageNames{storageIdx};
curveColor = curveColors(storageIdx, :);
filter = struct();
filter.pam_level = pamLevel;
filter.interference_path_length = plotOptions.path_length;
if hasBlockUpdate
filter.block_update = blockUpdates(blockIdx);
end
[sirValues, berGroups] = collectBerBySir(whConfig, storageName, filter);
if isempty(sirValues)
continue
end
[sirValues, sortIdx] = sort(sirValues);
berGroups = berGroups(sortIdx);
berAvg = nan(1, numel(sirValues));
berMin = nan(1, numel(sirValues));
berMax = nan(1, numel(sirValues));
for sirIdx = 1:numel(sirValues)
rawBer = berGroups{sirIdx};
rawBer = rawBer(isfinite(rawBer));
if isempty(rawBer)
continue
end
scatter(repmat(sirValues(sirIdx), size(rawBer)), rawBer, ...
12, ...
"Marker", ".", ...
"MarkerEdgeColor", curveColor, ...
"MarkerFaceColor", curveColor, ...
"HandleVisibility", "off");
berAvg(sirIdx) = mean(rawBer, "omitnan");
berMin(sirIdx) = min(rawBer);
berMax(sirIdx) = max(rawBer);
end
validAvg = isfinite(sirValues) & isfinite(berAvg);
if ~any(validAvg)
continue
end
if plotOptions.use_boundedlines && exist("boundedline", "file")
yLower = max(berAvg - berMin, 0);
yUpper = max(berMax - berAvg, 0);
yBounds = [yLower(:), yUpper(:)];
[hl, hp] = boundedline(sirValues(:), berAvg(:), yBounds, ...
"alpha", "transparency", plotOptions.boundedline_alpha, ...
"cmap", curveColor, ...
"nan", "fill", ...
"orientation", "vert");
set(hl, "LineStyle", "none", "Marker", "none", "HandleVisibility", "off");
set(hp, "HandleVisibility", "off", "LineStyle", "none");
end
scatter(sirValues(validAvg), berAvg(validAvg), ...
54, ...
"Marker", "o", ...
"MarkerEdgeColor", curveColor, ...
"MarkerFaceColor", "none", ...
"LineWidth", 1.2, ...
"DisplayName", storageName);
fitMask = validAvg & berAvg > 0;
if nnz(fitMask) >= 2
fitOrder = min(plotOptions.polyfit_order_max, nnz(fitMask) - 1);
fitCoeff = polyfit(sirValues(fitMask), log10(berAvg(fitMask)), fitOrder);
xFit = linspace(min(sirValues(fitMask)), max(sirValues(fitMask)), 300);
yFit = 10 .^ polyval(fitCoeff, xFit);
plot(xFit, yFit, ...
"LineWidth", 1.2, ...
"LineStyle", "--", ...
"Marker", "none", ...
"Color", curveColor, ...
"HandleVisibility", "off");
end
end
if hasBlockUpdate
title(sprintf("BER over SIR, PAM %.0f, path %.0f m, block update = %g", ...
pamLevel, plotOptions.path_length, blockUpdates(blockIdx)));
else
title(sprintf("BER over SIR, PAM %.0f, path %.0f m", ...
pamLevel, plotOptions.path_length));
end
xlabel("SIR (dB)");
ylabel("BER");
xlim(plotOptions.xlim);
grid on
box on
if exist("beautifyBERplot", "file")
beautifyBERplot("logscale", true, "setcolors", false, "setmarkers", false);
else
set(gca, "YScale", "log");
end
legend("Location", "best", "Interpreter", "none");
end
end
end
function [sirValues, berGroups] = collectBerBySir(whConfig, storageName, filter)
sirValues = [];
berGroups = {};
for linIdx = 1:numel(whConfig.sto.(storageName))
storedValue = whConfig.sto.(storageName){linIdx};
if isempty(storedValue)
continue
end
physStruct = linearIndexToStruct(whConfig, linIdx);
if ~matchesFilter(physStruct, filter)
continue
end
berValues = extractBerValues(storedValue);
if isempty(berValues)
continue
end
sirValue = physStruct.sir;
sirIdx = find(sirValues == sirValue, 1);
if isempty(sirIdx)
sirValues(end+1) = sirValue; %#ok<AGROW>
berGroups{end+1} = berValues; %#ok<AGROW>
else
berGroups{sirIdx} = [berGroups{sirIdx}, berValues]; %#ok<AGROW>
end
end
end
function tf = matchesFilter(physStruct, filter)
tf = true;
filterFields = fieldnames(filter);
for fieldIdx = 1:numel(filterFields)
fieldName = filterFields{fieldIdx};
if ~isfield(physStruct, fieldName) || physStruct.(fieldName) ~= filter.(fieldName)
tf = false;
return
end
end
end
function berValues = extractBerValues(value)
berValues = [];
if isstruct(value)
if isfield(value, "metrics") && isfield(value.metrics, "BER")
berValues(end+1) = value.metrics.BER;
end
elseif iscell(value)
for valueIdx = 1:numel(value)
berValues = [berValues, extractBerValues(value{valueIdx})]; %#ok<AGROW>
end
end
end

View File

@@ -0,0 +1,99 @@
function [summary, rows] = appendMpiReductionDspOutput(db, run_id, occurrence_idx, dspOutput, dsp_options, options)
%appendMpiReductionDspOutput Store one DSP occurrence in MpiReductionResults.
arguments
db
run_id
occurrence_idx
dspOutput struct
dsp_options struct = struct()
options.study_name string = "mpi_reduction_v1"
options.dry_run (1,1) logical = false
options.verbose (1,1) logical = false
end
summary = emptySummary();
rows = {};
packageNames = fieldnames(dspOutput);
if isempty(packageNames)
return
end
block_update = resolveBlockUpdate(dspOutput, dsp_options);
wh = DataStorage(struct( ...
"run_id", double(run_id), ...
"block_update", double(block_update)));
hasPackages = false;
for packageIdx = 1:numel(packageNames)
packageName = packageNames{packageIdx};
package = dspOutput.(packageName);
if isempty(package)
continue
end
if iscell(package)
packageCell = package;
else
packageCell = {package};
end
wh.addStorage(packageName);
wh.addValueToStorageByLinIdx(packageCell, packageName, 1);
hasPackages = true;
end
if ~hasPackages
return
end
dsp_options_local = dsp_options;
dsp_options_local.start_occurence = occurrence_idx;
dsp_options_local.userParameters.block_update = block_update;
[summary, rows] = appendMpiReductionWarehouse(db, wh, run_id, dsp_options_local, ...
"study_name", options.study_name, ...
"dry_run", options.dry_run, ...
"verbose", options.verbose);
end
function summary = emptySummary()
summary = struct( ...
"rows_considered", 0, ...
"rows_inserted", 0, ...
"rows_skipped_duplicate", 0, ...
"rows_skipped_existing_key", 0, ...
"rows_skipped_empty", 0);
end
function block_update = resolveBlockUpdate(dspOutput, dsp_options)
if isfield(dsp_options, "userParameters") && ...
isfield(dsp_options.userParameters, "block_update") && ...
isscalar(dsp_options.userParameters.block_update)
block_update = dsp_options.userParameters.block_update;
return
end
packageNames = fieldnames(dspOutput);
for packageIdx = 1:numel(packageNames)
package = dspOutput.(packageNames{packageIdx});
if isempty(package)
continue
end
if iscell(package)
package = package{1};
end
if isstruct(package) && ...
isfield(package, "mpi_reduction_config") && ...
isfield(package.mpi_reduction_config, "params") && ...
isfield(package.mpi_reduction_config.params, "block_update")
block_update = package.mpi_reduction_config.params.block_update;
return
end
end
error("appendMpiReductionDspOutput:MissingBlockUpdate", ...
"Could not resolve scalar block_update from dsp_options or package metadata.");
end

View File

@@ -0,0 +1,307 @@
function [summary, rows] = appendMpiReductionWarehouse(db, wh, run_ids, dsp_options, options)
%appendMpiReductionWarehouse Store MPI reduction warehouse entries in MySQL.
arguments
db
wh
run_ids
dsp_options struct
options.study_name string = "mpi_reduction_v1"
options.dry_run (1,1) logical = false
options.verbose (1,1) logical = true
end
ensureTableVisible(db);
storageNames = fieldnames(wh.sto);
summary = struct( ...
"rows_considered", 0, ...
"rows_inserted", 0, ...
"rows_skipped_duplicate", 0, ...
"rows_skipped_existing_key", 0, ...
"rows_skipped_empty", 0);
rows = {};
for storageIdx = 1:numel(storageNames)
storageName = storageNames{storageIdx};
storageArray = wh.sto.(storageName);
for linIdx = 1:numel(storageArray)
packageCell = storageArray{linIdx};
if isempty(packageCell)
summary.rows_skipped_empty = summary.rows_skipped_empty + 1;
continue
end
phys = physicalCoordinates(wh, linIdx);
run_id = resolveRunId(phys, run_ids);
block_update = resolveBlockUpdate(phys, dsp_options);
packageCell = normalizePackageCell(packageCell);
for packageIdx = 1:numel(packageCell)
package = packageCell{packageIdx};
if isempty(package) || ~isstruct(package) || ~isfield(package, "metrics")
summary.rows_skipped_empty = summary.rows_skipped_empty + 1;
continue
end
occurrence_idx = startOccurrence(dsp_options) + packageIdx - 1;
row = buildRow(db, package, storageName, run_id, occurrence_idx, ...
block_update, options.study_name);
summary.rows_considered = summary.rows_considered + 1;
rows{end+1,1} = row; %#ok<AGROW>
if options.dry_run
continue
end
if resultHashExists(db, row.result_hash)
summary.rows_skipped_duplicate = summary.rows_skipped_duplicate + 1;
if options.verbose
fprintf("MPI reduction DB: skipped duplicate %s, run_id=%d, occurrence=%d, block_update=%d\n", ...
row.storage_name, row.run_id, row.occurrence_idx, row.block_update);
end
continue
end
if uniqueResultKeyExists(db, row)
summary.rows_skipped_existing_key = summary.rows_skipped_existing_key + 1;
if options.verbose
fprintf("MPI reduction DB: skipped existing key %s, run_id=%d, occurrence=%d, block_update=%d\n", ...
row.storage_name, row.run_id, row.occurrence_idx, row.block_update);
end
continue
end
db.appendToTable("MpiReductionResults", row);
summary.rows_inserted = summary.rows_inserted + 1;
end
end
end
if options.verbose
fprintf("MPI reduction DB: considered %d rows, inserted %d, duplicates %d, existing keys %d, empty %d\n", ...
summary.rows_considered, summary.rows_inserted, ...
summary.rows_skipped_duplicate, summary.rows_skipped_existing_key, ...
summary.rows_skipped_empty);
end
end
function ensureTableVisible(db)
if ~isfield(db.tables, "MpiReductionResults")
db.refresh();
end
if ~isfield(db.tables, "MpiReductionResults")
error("appendMpiReductionWarehouse:MissingTable", ...
"Table MpiReductionResults is not visible to DBHandler. Create the table and call db.refresh().");
end
end
function phys = physicalCoordinates(wh, linIdx)
phys = struct();
[physValues, physNames] = wh.getPhysIndicesByLinIndex(linIdx);
for physIdx = 1:numel(physNames)
name = char(physNames{physIdx});
phys.(name) = physValues{physIdx};
end
end
function run_id = resolveRunId(phys, run_ids)
if isfield(phys, "run_id")
run_id = phys.run_id;
return
end
if isscalar(run_ids)
run_id = run_ids;
return
end
error("appendMpiReductionWarehouse:MissingRunId", ...
"Warehouse has no run_id axis, but run_ids is not scalar.");
end
function block_update = resolveBlockUpdate(phys, dsp_options)
if isfield(phys, "block_update")
block_update = phys.block_update;
return
end
if isfield(dsp_options, "userParameters") && ...
isfield(dsp_options.userParameters, "block_update") && ...
isscalar(dsp_options.userParameters.block_update)
block_update = dsp_options.userParameters.block_update;
return
end
error("appendMpiReductionWarehouse:MissingBlockUpdate", ...
"Could not resolve block_update from warehouse coordinates or scalar userParameters.");
end
function packageCell = normalizePackageCell(packageCell)
if ~iscell(packageCell)
packageCell = {packageCell};
end
end
function occurrence = startOccurrence(dsp_options)
occurrence = 1;
if isfield(dsp_options, "start_occurence")
occurrence = dsp_options.start_occurence;
end
end
function row = buildRow(db, package, storageName, run_id, occurrence_idx, block_update, study_name)
config = fieldOr(package, "mpi_reduction_config", fallbackConfig(storageName, block_update));
params = fieldOr(config, "params", struct());
if isempty(params)
params = struct();
end
metrics = package.metrics;
if isobject(metrics) && ismethod(metrics, "toStruct")
metrics = metrics.toStruct();
end
param_json = jsonencode(params);
if strlength(string(param_json)) == 0
param_json = "{}";
end
row = struct();
row.run_id = double(run_id);
row.occurrence_idx = double(occurrence_idx);
row.study_name = char(study_name);
row.storage_name = char(fieldOr(config, "storage_name", storageName));
row.algorithm = char(fieldOr(config, "algorithm", "unknown"));
row.algorithm_variant = char(fieldOr(config, "algorithm_variant", ""));
row.eq_class = char(fieldOr(config, "eq_class", "unknown"));
row.block_update = double(block_update);
row.sps = scalarParam(params, "sps");
row.eq_order = scalarParam(params, "order");
row.len_tr = scalarParam(params, "len_tr");
row.epochs_tr = scalarParam(params, "epochs_tr");
row.mu_tr = scalarParam(params, "mu_tr");
row.dd_mode = logicalParam(params, "dd_mode");
row.epochs_dd = scalarParam(params, "epochs_dd");
row.mu_dd = scalarParam(params, "mu_dd");
row.numBits = metricScalar(metrics, "numBits", 0);
row.numBitErr = metricScalar(metrics, "numBitErr", 0);
row.BER = metricScalar(metrics, "BER", NaN);
row.numBitErr_precoded = metricScalar(metrics, "numBitErr_precoded", NaN);
row.BER_precoded = metricScalar(metrics, "BER_precoded", NaN);
row.SNR = metricScalar(metrics, "SNR", NaN);
row.SNR_level = metricJson(metrics, "SNR_level");
row.GMI = metricScalar(metrics, "GMI", NaN);
row.AIR = metricScalar(metrics, "AIR", NaN);
row.EVM = metricScalar(metrics, "EVM", NaN);
row.EVM_level = metricJson(metrics, "EVM_level");
row.STD = metricScalar(metrics, "STD", NaN);
row.STD_level = metricJson(metrics, "STD_level");
row.STDrx = metricScalar(metrics, "STDrx", NaN);
row.STDrx_level = metricJson(metrics, "STDrx_level");
row.Alpha = firstFiniteScalar(fieldOr(metrics, "Alpha", NaN));
row.param_json = char(param_json);
row.param_hash = db.calcHash(params);
row.result_hash = db.calcHash(row);
end
function config = fallbackConfig(storageName, block_update)
params = struct("block_update", block_update);
config = struct( ...
"storage_name", storageName, ...
"algorithm", "unknown", ...
"algorithm_variant", "", ...
"eq_class", "unknown", ...
"params", params);
end
function value = fieldOr(s, fieldName, defaultValue)
fieldName = char(fieldName);
if isstruct(s) && isfield(s, fieldName) && ~isempty(s.(fieldName))
value = s.(fieldName);
else
value = defaultValue;
end
end
function value = scalarParam(params, fieldName)
value = firstFiniteScalar(fieldOr(params, fieldName, NaN));
end
function value = logicalParam(params, fieldName)
value = firstFiniteScalar(fieldOr(params, fieldName, false));
if isnan(value)
value = false;
else
value = logical(value);
end
end
function value = metricScalar(metrics, fieldName, defaultValue)
value = firstFiniteScalar(fieldOr(metrics, fieldName, defaultValue));
end
function value = metricJson(metrics, fieldName)
metricValue = fieldOr(metrics, fieldName, []);
value = char(jsonencode(metricValue));
end
function value = firstFiniteScalar(valueIn)
if islogical(valueIn)
valueIn = double(valueIn);
end
if isempty(valueIn)
value = NaN;
return
end
if isnumeric(valueIn)
finiteValues = valueIn(isfinite(valueIn));
if isempty(finiteValues)
value = NaN;
else
value = double(finiteValues(1));
end
return
end
value = str2double(string(valueIn));
if isnan(value)
value = NaN;
end
end
function exists = resultHashExists(db, resultHash)
query = sprintf("SELECT mpi_result_id FROM MpiReductionResults WHERE result_hash = '%s' LIMIT 1", ...
char(resultHash));
existing = db.fetch(query);
exists = ~isempty(existing);
end
function exists = uniqueResultKeyExists(db, row)
query = sprintf( ...
"SELECT mpi_result_id FROM MpiReductionResults " + ...
"WHERE run_id = %d " + ...
"AND occurrence_idx = %d " + ...
"AND study_name = '%s' " + ...
"AND storage_name = '%s' " + ...
"AND block_update = %d " + ...
"AND param_hash = '%s' LIMIT 1", ...
row.run_id, row.occurrence_idx, ...
sqlString(row.study_name), sqlString(row.storage_name), ...
row.block_update, sqlString(row.param_hash));
existing = db.fetch(query);
exists = ~isempty(existing);
end
function value = sqlString(value)
value = strrep(char(value), '''', '''''');
end

View File

@@ -4,8 +4,14 @@ dsp_options.mode = "run_id";
dsp_options.recipe = @mpi_recipe_dev; dsp_options.recipe = @mpi_recipe_dev;
dsp_options.append_to_db = false; dsp_options.append_to_db = false;
dsp_options.start_occurence = 1; dsp_options.start_occurence = 1;
dsp_options.max_occurences = 10; dsp_options.max_occurences = 15;
dsp_options.debug_plots = true; dsp_options.debug_plots = false;
write_mpi_reduction_db = 1;
stream_mpi_reduction_db = 1;
mpi_reduction_study_name = "block_update_sweep";
mpi_reduction_writer_path = fullfile(fileparts(mfilename('fullpath')),"db");
addpath(mpi_reduction_writer_path);
dsp_options.database_type = "mysql"; dsp_options.database_type = "mysql";
@@ -16,6 +22,9 @@ dsp_options.server = "192.168.178.192";
dsp_options.port = 3306; dsp_options.port = 3306;
dsp_options.user = "silas"; dsp_options.user = "silas";
dsp_options.password = "silas"; dsp_options.password = "silas";
dsp_options.append_mpi_reduction_db = write_mpi_reduction_db && stream_mpi_reduction_db;
dsp_options.mpi_reduction_study_name = mpi_reduction_study_name;
dsp_options.mpi_reduction_writer_path = mpi_reduction_writer_path;
db = DBHandler("dataBase", [dsp_options.dataBase],... db = DBHandler("dataBase", [dsp_options.dataBase],...
"type", dsp_options.database_type,... "type", dsp_options.database_type,...
@@ -23,37 +32,41 @@ db = DBHandler("dataBase", [dsp_options.dataBase],...
"user", dsp_options.user, "password", dsp_options.password); "user", dsp_options.user, "password", dsp_options.password);
%% %%
pamformats = [4,6,8];
baudrates = [112e9,96e9,72e9];
for i = 1
B = baudrates(i);
M = pamformats(i);
fp = QueryFilter(); fp = QueryFilter();
fp.where('Runs','run_id','GREATER_EQUAL',3153); fp.where('Runs','run_id','GREATER_EQUAL',3153);
fp.where('Runs', 'symbolrate', 'EQUALS', 112e9); fp.where('Runs', 'symbolrate', 'EQUALS', B); % 72 96 112
fp.where('Runs', 'fiber_length', 'EQUALS', 0); fp.where('Runs', 'fiber_length', 'EQUALS', 0);
fp.where('Runs', 'interference_path_length', 'EQUALS', 1000); fp.where('Runs', 'interference_path_length', 'EQUALS', 1000);
fp.where('Runs', 'sir', 'LESS_EQUAL', 30); % fp.where('Runs', 'sir', 'LESS_EQUAL', 20);
fp.where('Runs', 'db_mode', 'EQUALS', '"no_db"'); fp.where('Runs', 'db_mode', 'EQUALS', '"no_db"');
% fp.where('Runs', 'is_mpi', 'EQUALS', 1); % fp.where('Runs', 'is_mpi', 'EQUALS', 1);
fp.where('Runs', 'pam_level', 'EQUALS', 4); fp.where('Runs', 'pam_level', 'EQUALS', M);
% fp.where('Runs', 'v_bias', 'EQUALS', 2.65); % fp.where('Runs', 'v_bias', 'EQUALS', 2.65);
[dataTable,~] = db.queryDB(fp, db.getTableFieldNames('Runs')); [dataTable,~] = db.queryDB(fp, db.getTableFieldNames('Runs'));
% [~, uniqueSirRows] = unique(dataTable.sir, "stable");
% dataTable = dataTable(uniqueSirRows, :);
[~, sortIdx] = sort(dataTable.sir, 'descend'); [~, sortIdx] = sort(dataTable.sir, 'descend');
dataTable = dataTable(sortIdx, :); dataTable = dataTable(sortIdx, :);
run_ids = dataTable.run_id;
if isempty(run_ids) % dataTable = dataTable(1,:);
error("run_minimal_recipe:MissingRunIds", ...
"Set run_ids to one or more known run IDs before running this example."); run_ids = dataTable.run_id;
end
%% === Warehouse setup === %% === Warehouse setup ===
dsp_options.userParameters = struct(); dsp_options.userParameters = struct();
% dsp_options.userParameters.smoothing_length = [0,linspace(100,1000,10),2500,5000,10000]; dsp_options.userParameters.block_update = [1,2,4,8,16,32,64,112,224,448,448*2,1024,2048,4096,8192,16384];%%logspace(-3.8,-1,22);%[linspace(2,4096,22)];
% dsp_options.userParameters.block_update = linspace(1,224,22);
wh = DataStorage(dsp_options.userParameters); wh = DataStorage(dsp_options.userParameters);
%% %%
n_realizations = (dsp_options.max_occurences - dsp_options.start_occurence + 1); n_realizations = (dsp_options.max_occurences - dsp_options.start_occurence + 1);
n_userparams = prod(wh.dim); n_userparams = prod(wh.dim);
@@ -67,38 +80,74 @@ fprintf("-> [ %d run_id(s) × %d userParam combination(s) = %d parallel job(s) ]
%% === Run === %% === Run ===
% submitJobs returns the raw per-job results and the filled Warehouse. For a % 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. % 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, ... "wh", wh, ...
"waitbar", true); "waitbar", true);
%%
if write_mpi_reduction_db && ~stream_mpi_reduction_db
db.refresh();
appendMpiReductionWarehouse(db, wh, run_ids, dsp_options, ...
"study_name", mpi_reduction_study_name);
end
% savepath = fullfile("C:","Users","Silas","Documents","MATLAB","imdd_simulation","projects","Diss","MPI_revisit","algorithms",sprintf("pam_%d_results.mat",M));
% save(wh,savepath)
end
%% Analyze %% Analyze
storageNames = fieldnames(wh.sto); storageNames = fieldnames(wh.sto);
x_vars = run_ids;%dsp_options.userParameters.smoothing_length;
x_vars = reshape(x_vars,1,[]);
figure(2026);hold on x_base = dataTable.sir(:).';
for i = 1:numel(run_ids) % x_base = dsp_options.userParameters.block_update;
disp(run_ids(i))
result = wh.getStoValue(storageNames{1}, x_vars); figure(2026); clf; hold on
% Each cell in result is a cell array (packageCell) containing multiple packages. for storage_idx = 1:numel(storageNames)
% Extract BERs from all packages and, for example, average them per x_var. storageName = storageNames{storage_idx};
ber_all = cellfun(@(packageCell) cellfun(@(pkg) pkg.metrics.BER, packageCell), result, 'UniformOutput', false); result = wh.sto.(storageName);
% Convert to numeric matrix (rows = packages, cols = x_vars) by padding if needed 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)); maxPackages = max(cellfun(@numel,ber_all));
ber_mat = nan(maxPackages,numel(ber_all)); ber_mat = nan(maxPackages,numel(ber_all));
for k = 1:numel(ber_all) for k = 1:numel(ber_all)
ber_mat(1:numel(ber_all{k}),k) = ber_all{k}; ber_mat(1:numel(ber_all{k}),k) = ber_all{k};
end end
% Choose aggregation: mean across packages (ignore NaNs) ber = mean(ber_mat,1,"omitnan");
ber = mean(ber_mat, 1, 'omitnan');
x = x_base;
x = dataTable.sir; if numel(x) ~= numel(ber)
y = ber; x = 1:numel(ber);
plot(x,ber);
scatter(x,ber_mat)
beautifyBERplot("logscale",true,"setcolors",false,"setmarkers",true);
end 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");

View File

@@ -0,0 +1,73 @@
% === DSP settings ===
dsp_options = struct();
dsp_options.mode = "run_id";
dsp_options.recipe = @mpi_recipe;
dsp_options.append_to_db = false;
dsp_options.start_occurence = 1;
dsp_options.max_occurences = 15;
dsp_options.debug_plots = false;
dsp_options.database_type = "mysql";
dsp_options.dataBase = "labor";
dsp_options.storage_path = "W:\labdata\ECOC Silas\ecoc_2025";
dsp_options.server = "192.168.178.192";
dsp_options.port = 3306;
dsp_options.user = "silas";
dsp_options.password = "silas";
db = DBHandler("dataBase", [dsp_options.dataBase],...
"type", dsp_options.database_type,...
"server", dsp_options.server,...
"user", dsp_options.user, "password", dsp_options.password);
%%
fp = QueryFilter();
fp.where('Runs','run_id','GREATER_EQUAL',3153);
fp.where('Runs', 'symbolrate', 'EQUALS', 72e9); % 72 96 112
fp.where('Runs', 'fiber_length', 'EQUALS', 0);
% fp.where('Runs', 'interference_path_length', 'EQUALS', 1000);
% fp.where('Runs', 'sir', 'LESS_EQUAL', 50);
fp.where('Runs', 'db_mode', 'EQUALS', '"no_db"');
% fp.where('Runs', 'is_mpi', 'EQUALS', 1);
fp.where('Runs', 'pam_level', 'EQUALS', 8);
% fp.where('Runs', 'v_bias', 'EQUALS', 2.65);
[dataTable,~] = db.queryDB(fp, db.getTableFieldNames('Runs'));
% [~, uniqueSirRows] = unique(dataTable.sir, "stable");
% dataTable = dataTable(uniqueSirRows, :);
[~, sortIdx] = sort(dataTable.sir, 'descend');
dataTable = dataTable(sortIdx, :);
% dataTable = dataTable(1, :);
run_ids = dataTable.run_id;
if isempty(run_ids)
error("run_minimal_recipe:MissingRunIds", ...
"Set run_ids to one or more known run IDs before running this example.");
end
%% === Warehouse setup ===
dsp_options.userParameters = struct();
dsp_options.userParameters.block_update = [1,2,4,8,16,32,64,112,224,448,512,1024,2048,2048*2,2048*4];%%logspace(-3.8,-1,22);%[linspace(2,4096,22)];
wh = DataStorage(dsp_options.userParameters);
%%
n_realizations = (dsp_options.max_occurences - dsp_options.start_occurence + 1);
n_userparams = prod(wh.dim);
n_run_ids = numel(run_ids);
parallel_jobs = n_userparams * n_run_ids;
queried_jobs = n_realizations * n_userparams * n_run_ids;
fprintf("-> [ %d run_id(s) × %d userParam combination(s) = %d parallel job(s) ] × %d realizations = %d total jobs \n", ...
n_run_ids, n_userparams, parallel_jobs, n_realizations, queried_jobs);
%% === Run ===
% submitJobs returns the raw per-job results and the filled Warehouse. For a
% single run_id and remove_dc = [0, 1], results is a 2-by-1 cell array.
[results, wh] = submitJobs(run_ids, dsp_options, processingMode.parallel, ...
"wh", wh, ...
"waitbar", true);

View File

@@ -0,0 +1,389 @@
%% Analyze
% === DSP settings ===
dsp_options = struct();
dsp_options.mode = "run_id";
dsp_options.recipe = @mpi_recipe_dev;
dsp_options.append_to_db = false;
dsp_options.start_occurence = 1;
dsp_options.max_occurences = 15;
dsp_options.debug_plots = false;
dsp_options.database_type = "mysql";
dsp_options.dataBase = "labor";
dsp_options.storage_path = "W:\labdata\ECOC Silas\ecoc_2025";
dsp_options.server = "192.168.178.192";
dsp_options.port = 3306;
dsp_options.user = "silas";
dsp_options.password = "silas";
db = DBHandler("dataBase", [dsp_options.dataBase],...
"type", dsp_options.database_type,...
"server", dsp_options.server,...
"user", dsp_options.user, "password", dsp_options.password);
%%
fp = QueryFilter();
fp.where('Runs','run_id','GREATER_EQUAL',3153);
fp.where('Runs', 'symbolrate', 'EQUALS', 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', 4);
% fp.where('Runs', 'v_bias', 'EQUALS', 2.65);
[dataTable,~] = db.queryDB(fp, db.getTableFieldNames('Runs'));
% [~, uniqueSirRows] = unique(dataTable.sir, "stable");
% dataTable = dataTable(uniqueSirRows, :);
[~, sortIdx] = sort(dataTable.sir, 'descend');
dataTable = dataTable(sortIdx, :);
% dataTable = dataTable(1, :);
run_ids = dataTable.run_id;
%%
wh_ = load("C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\Diss\MPI_revisit\parallelization_analysis\wh_block_update_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");
% wh2 = wh_.wh;
%
% wh = mergeDataStorages(wh1,wh2);
%%
storageNames = fieldnames(wh.sto);
storageNames = storageNames([1,2,3,5]);
ber_by_storage = struct();
ber_mat_by_storage = struct();
sir_values = dataTable.sir(:).';
block_updates = wh.parameter.block_update.values;
storage_parameter_names = cellstr(wh.fn);
block_axis = find(strcmp(storage_parameter_names,"block_update"),1);
if exist('linspecer','file')
curve_colors = linspecer(max(numel(storageNames),1));
else
curve_colors = lines(max(numel(storageNames),1));
end
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)
storageName = storageNames{storage_idx};
fit_coeff_by_storage.(storageName) = cell(1,numel(block_updates));
required_sir_fec_by_storage.(storageName) = nan(1,numel(block_updates));
end
for block_idx = 1:numel(block_updates)
block_update = block_updates(block_idx);
figure(3000 + block_idx); clf; hold on
for storage_idx = 1:numel(storageNames)
storageName = storageNames{storage_idx};
curve_color = curve_colors(storage_idx,:);
result = wh.sto.(storageName);
if ~isempty(block_axis) && size(result,block_axis) == numel(block_updates)
result_subs = repmat({':'},1,ndims(result));
result_subs{block_axis} = block_idx;
result_for_block = result(result_subs{:});
result_for_block = result_for_block(:).';
else
result_for_block = result(:).';
end
% Each stored entry is a package cell containing one or more occurrences.
ber_all = cell(1,numel(result_for_block));
for result_idx = 1:numel(result_for_block)
packageCell = result_for_block{result_idx};
if isempty(packageCell)
ber_all{result_idx} = NaN;
continue
end
ber_values = nan(1,numel(packageCell));
for package_idx = 1:numel(packageCell)
if isstruct(packageCell{package_idx}) && isfield(packageCell{package_idx},"metrics")
ber_values(package_idx) = packageCell{package_idx}.metrics.BER;
end
end
ber_all{result_idx} = ber_values;
end
maxPackages = max(cellfun(@numel,ber_all));
ber_mat = nan(maxPackages,numel(ber_all));
for k = 1:numel(ber_all)
ber_mat(1:numel(ber_all{k}),k) = ber_all{k};
end
% 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);
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");
ber_min = nan(1,numel(ber));
ber_max = nan(1,numel(ber));
for k = 1:numel(ber)
ber_here = ber_mat(:,k);
ber_here = ber_here(isfinite(ber_here));
if isempty(ber_here)
continue
end
ber_min(k) = min(ber_here);
ber_max(k) = max(ber_here);
end
x = sir_values;
if numel(x) ~= numel(ber)
x = 1:numel(ber);
end
ber_by_storage.(storageName){block_idx} = ber;
ber_mat_by_storage.(storageName){block_idx} = ber_mat;
for k = 1:numel(x)
scatter(repmat(x(k),maxPackages,1),ber_mat(:,k), ...
30, ...
'Marker','.', ...
'MarkerEdgeColor',curve_color, ...
'MarkerFaceColor',curve_color, ...
'HandleVisibility','off');
end
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);
fit_coeff_by_storage.(storageName){block_idx} = fit_coeff;
sir_req = NaN;
threshold_mask = isfinite(y_fit) & y_fit <= fec_ber_threshold;
if any(threshold_mask)
first_threshold_idx = find(threshold_mask,1,"first");
if first_threshold_idx == 1
sir_req = x_fit(first_threshold_idx);
else
x_pair = x_fit(first_threshold_idx-1:first_threshold_idx);
y_pair = log10(y_fit(first_threshold_idx-1:first_threshold_idx));
if all(isfinite(y_pair)) && diff(y_pair) ~= 0
sir_req = interp1(y_pair,x_pair,log10(fec_ber_threshold), ...
"linear","extrap");
else
sir_req = x_fit(first_threshold_idx);
end
end
end
required_sir_fec_by_storage.(storageName)(block_idx) = sir_req;
plot(x_fit,y_fit, ...
'LineWidth',1.1, ...
'LineStyle','--', ...
'Color',curve_color, ...
'HandleVisibility','off');
end
end
title(sprintf("BER over SIR, block update = %g",block_update));
xlabel("SIR (dB)");
ylabel("BER");
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
figure(5000); clf; hold on
for storage_idx = 1:numel(storageNames)
storageName = storageNames{storage_idx};
curve_color = curve_colors(storage_idx,:);
required_sir = required_sir_fec_by_storage.(storageName);
valid = isfinite(block_updates) & isfinite(required_sir);
if ~any(valid)
continue
end
plot(block_updates(valid),required_sir(valid), ...
'LineWidth',1.4, ...
'LineStyle','-', ...
'Marker','o', ...
'MarkerSize',5, ...
'Color',curve_color, ...
'DisplayName',storageName);
end
set(gca,'XScale','log');
xticks(block_updates);
xticklabels(string(block_updates));
grid on
box on
xlabel("Parallelization / block update");
ylabel(sprintf("Required SIR at BER = %.1e (dB)",fec_ber_threshold));
title("Required SIR to reach FEC over parallelization");
legend("Location","best","Interpreter","none");
%%
function whMerged = mergeDataStorages(varargin)
% mergeDataStorages merges compatible DataStorage objects by physical indices.
% Existing duplicate entries are concatenated when both values are package cells.
whList = varargin;
templateWh = whList{1};
paramNames = cellstr(templateWh.fn);
mergedParams = struct();
for param_idx = 1:numel(paramNames)
paramName = paramNames{param_idx};
values = templateWh.parameter.(paramName).values(:).';
for wh_idx = 2:numel(whList)
curWh = whList{wh_idx};
if ~isequal(sort(cellstr(curWh.fn)),sort(paramNames))
error("mergeDataStorages:ParameterMismatch", ...
"All warehouses must use the same parameter names.");
end
if ~isfield(curWh.parameter,paramName)
error("mergeDataStorages:ParameterMismatch", ...
"Warehouse %d does not contain parameter '%s'.",wh_idx,paramName);
end
newValues = curWh.parameter.(paramName).values(:).';
for value_idx = 1:numel(newValues)
if ~ismember(newValues(value_idx),values)
values(end+1) = newValues(value_idx); %#ok<AGROW>
end
end
end
if strcmp(paramName,"block_update") && isnumeric(values)
values = sort(values);
end
mergedParams.(paramName) = values;
end
whMerged = DataStorage(mergedParams);
for wh_idx = 1:numel(whList)
curWh = whList{wh_idx};
curStorageNames = fieldnames(curWh.sto);
for storage_idx = 1:numel(curStorageNames)
storageName = curStorageNames{storage_idx};
if ~isfield(whMerged.sto,storageName)
whMerged.addStorage(storageName);
end
for lin_idx = 1:numel(curWh.sto.(storageName))
storedValue = curWh.sto.(storageName){lin_idx};
if isempty(storedValue)
continue
end
[physValues,physNames] = curWh.getPhysIndicesByLinIndex(lin_idx);
physNames = cellfun(@char,physNames,'UniformOutput',false);
targetSubscripts = cell(1,numel(whMerged.fn));
for param_idx = 1:numel(whMerged.fn)
paramName = char(whMerged.fn(param_idx));
sourceParamIdx = find(strcmp(physNames,paramName),1);
if isempty(sourceParamIdx)
error("mergeDataStorages:ParameterMismatch", ...
"Storage entry is missing parameter '%s'.",paramName);
end
targetSubscripts{param_idx} = whMerged.getIndexByPhys(paramName,physValues{sourceParamIdx});
end
if isscalar(targetSubscripts)
targetLinIdx = targetSubscripts{1};
else
targetLinIdx = sub2ind(whMerged.getStorageSize(),targetSubscripts{:});
end
existingValue = whMerged.sto.(storageName){targetLinIdx};
whMerged.sto.(storageName){targetLinIdx} = mergeStoredValue(existingValue,storedValue);
end
end
end
end
function out = mergeStoredValue(existingValue,newValue)
if isempty(existingValue)
out = newValue;
elseif iscell(existingValue) && iscell(newValue)
out = [existingValue(:); newValue(:)].';
else
out = newValue;
end
end

View File

@@ -127,7 +127,9 @@ for sirval = [20,30,50]
pngTikzPath = "04_Experimental_Evaluation/tikz/mpi/" + exportName + "-1.png"; pngTikzPath = "04_Experimental_Evaluation/tikz/mpi/" + exportName + "-1.png";
figure(sirval + pathlen); figure(sirval + pathlen);
exportNakedPlotWithTikz(gca, pngFile, tikzFile, pngTikzPath, 150); exportLevelScatterTikz(gca, pngFile, tikzFile, pngTikzPath, ...
"resolutionDpi", 150, ...
"generatedBy", "plot_mpi_timesignal.m");
end end
end end
@@ -145,309 +147,3 @@ end
% var_slope = p(1); % var_slope = p(1);
% var_offset = p(2);c % var_offset = p(2);c
function exportNakedPlotWithTikz(sourceAx, pngFile, tikzFile, pngTikzPath, resolutionDpi)
outputDir = fileparts(pngFile);
if ~exist(outputDir, "dir")
mkdir(outputDir);
end
exportScatterLayerPng(sourceAx, pngFile, resolutionDpi);
linePlots = collectTikzLinePlots(sourceAx);
textAnnotations = collectTikzTextAnnotations(sourceAx);
writeTikzPngAxis(tikzFile, pngTikzPath, sourceAx, linePlots, textAnnotations);
end
function exportScatterLayerPng(sourceAx, pngFile, resolutionDpi)
sourceFig = ancestor(sourceAx, "figure");
xLimits = sourceAx.XLim;
yLimits = sourceAx.YLim;
scatterHandles = findGraphicsObjectsByType(sourceAx, "scatter");
exportFig = figure( ...
"Visible", "off", ...
"Color", "white", ...
"Units", sourceFig.Units, ...
"Position", sourceFig.Position);
cleanupFigure = onCleanup(@() close(exportFig));
exportAx = copyobj(sourceAx, exportFig);
exportAx.Units = "normalized";
exportAx.Position = [0 0 1 1];
exportAx.XLim = xLimits;
exportAx.YLim = yLimits;
exportAx.XLimMode = "manual";
exportAx.YLimMode = "manual";
stripAxesToScatterOnly(exportAx);
hideAxesInkForRasterExport(exportAx);
fprintf("Exporting scatter PNG with XLim=[%g %g], YLim=[%g %g], scatter objects=%d: %s\n", ...
xLimits(1), xLimits(2), yLimits(1), yLimits(2), numel(scatterHandles), pngFile);
drawnow;
exportFig.PaperPositionMode = "auto";
print(exportFig, char(pngFile), "-dpng", sprintf("-r%d", resolutionDpi));
end
function stripAxesToScatterOnly(ax)
plotObjects = findall(ax);
for objIdx = 1:numel(plotObjects)
obj = plotObjects(objIdx);
if obj == ax || ~isvalid(obj)
continue
end
objType = string(get(obj, "Type"));
if lower(objType) ~= "scatter"
delete(obj);
end
end
end
function hideAxesInkForRasterExport(ax)
ax.Visible = "on";
ax.Color = "white";
ax.Box = "off";
ax.XGrid = "off";
ax.YGrid = "off";
ax.XMinorGrid = "off";
ax.YMinorGrid = "off";
ax.XTick = [];
ax.YTick = [];
ax.XColor = "white";
ax.YColor = "white";
ax.Title.String = "";
ax.XLabel.String = "";
ax.YLabel.String = "";
end
function handles = findGraphicsObjectsByType(parentHandle, objectType)
allHandles = findall(parentHandle);
matches = false(size(allHandles));
for handleIdx = 1:numel(allHandles)
currentHandle = allHandles(handleIdx);
if ~isvalid(currentHandle)
continue
end
matches(handleIdx) = lower(string(get(currentHandle, "Type"))) == lower(string(objectType));
end
handles = allHandles(matches);
end
function linePlots = collectTikzLinePlots(sourceAx)
lineHandles = findGraphicsObjectsByType(sourceAx, "line");
lineHandles = flipud(lineHandles(:));
linePlots = repmat(struct( ...
"xData", [], ...
"yData", [], ...
"color", [], ...
"lineWidth", []), numel(lineHandles), 1);
validLineCount = 0;
for lineIdx = 1:numel(lineHandles)
lineHandle = lineHandles(lineIdx);
xData = lineHandle.XData(:);
yData = lineHandle.YData(:);
validSamples = isfinite(xData) & isfinite(yData);
if ~any(validSamples)
continue
end
validLineCount = validLineCount + 1;
linePlots(validLineCount).xData = xData(validSamples);
linePlots(validLineCount).yData = yData(validSamples);
linePlots(validLineCount).color = lineHandle.Color;
linePlots(validLineCount).lineWidth = lineHandle.LineWidth;
end
linePlots = linePlots(1:validLineCount);
end
function textAnnotations = collectTikzTextAnnotations(sourceAx)
textHandles = findGraphicsObjectsByType(sourceAx, "text");
textHandles = flipud(textHandles(:));
xLimits = sourceAx.XLim;
yLimits = sourceAx.YLim;
textAnnotations = repmat(struct( ...
"x", [], ...
"y", [], ...
"text", "", ...
"anchor", "center"), numel(textHandles), 1);
validTextCount = 0;
for textIdx = 1:numel(textHandles)
textHandle = textHandles(textIdx);
labelText = string(textHandle.String);
if strlength(strtrim(labelText)) == 0
continue
end
textPosition = textHandle.Position;
if textPosition(1) < xLimits(1) || textPosition(1) > xLimits(2) || ...
textPosition(2) < yLimits(1) || textPosition(2) > yLimits(2)
continue
end
validTextCount = validTextCount + 1;
textAnnotations(validTextCount).x = textPosition(1);
textAnnotations(validTextCount).y = textPosition(2);
textAnnotations(validTextCount).text = labelText;
textAnnotations(validTextCount).anchor = matlabTextAlignmentToTikzAnchor( ...
textHandle.HorizontalAlignment, textHandle.VerticalAlignment);
end
textAnnotations = textAnnotations(1:validTextCount);
end
function writeTikzPngAxis(tikzFile, pngTikzPath, sourceAx, linePlots, textAnnotations)
outputDir = fileparts(tikzFile);
if ~exist(outputDir, "dir")
mkdir(outputDir);
end
xLimits = sourceAx.XLim;
yLimits = sourceAx.YLim;
fid = fopen(tikzFile, "w");
if fid < 0
error("plot_mpi_timesignal:TikzOpenFailed", ...
"Could not open TikZ export file: %s", tikzFile);
end
cleanupFile = onCleanup(@() fclose(fid));
axisOptions = {
" every axis/.append style={font=\scriptsize},"
"width=\fwidth,"
"height=\fheight,"
"at={(0\fwidth,0\fheight)},"
"scale only axis,"
"axis on top,"
"xmin=" + formatPgfNumber(xLimits(1)) + ","
"xmax=" + formatPgfNumber(xLimits(2)) + ","
"xlabel style={font=\color{white!15!black}\small},"
"xlabel={Time in $\mu$s},"
"ymin=" + formatPgfNumber(yLimits(1)) + ","
"ymax=" + formatPgfNumber(yLimits(2)) + ","
"ylabel style={font=\color{white!15!black}\small},"
"ylabel={Normalized Amplitude},"
"axis background/.style={fill=white},"
"xmajorgrids,"
"ymajorgrids,"
"grid style={line width=0.4pt, dotted, color=black!20},"
"scaled ticks=false,"
"tick label style={/pgf/number format/fixed, /pgf/number format/1000 sep={}},"
"legend columns=1"
};
fprintf(fid, "%% This file was generated by plot_mpi_timesignal.m.\n");
fprintf(fid, "%%\n");
for lineIdx = 1:numel(linePlots)
fprintf(fid, "%s\n", formatTikzColorDefinition(lineIdx, linePlots(lineIdx).color));
end
if ~isempty(linePlots)
fprintf(fid, "%%\n");
end
fprintf(fid, "%s\n\n", "\begin{tikzpicture}");
fprintf(fid, "%s\n", "\begin{axis}[%");
for optionIdx = 1:numel(axisOptions)
fprintf(fid, "%s\n", axisOptions{optionIdx});
end
fprintf(fid, "%s\n", "]");
graphicsLine = sprintf( ...
"\\addplot [forget plot] graphics [xmin=%s, xmax=%s, ymin=%s, ymax=%s] {%s};", ...
formatPgfNumber(xLimits(1)), ...
formatPgfNumber(xLimits(2)), ...
formatPgfNumber(yLimits(1)), ...
formatPgfNumber(yLimits(2)), ...
char(strrep(pngTikzPath, "\", "/")));
fprintf(fid, "%s\n\n", graphicsLine);
writeTikzLinePlots(fid, linePlots);
writeTikzTextAnnotations(fid, textAnnotations);
fprintf(fid, "%s\n\n", "\end{axis}");
fprintf(fid, "%s", "\end{tikzpicture}%");
end
function valueText = formatPgfNumber(value)
valueText = string(sprintf("%.15g", value));
end
function colorDefinition = formatTikzColorDefinition(colorIdx, rgbColor)
colorDefinition = sprintf( ...
"\\definecolor{mycolor%d}{rgb}{%.5f,%.5f,%.5f}%%", ...
colorIdx, rgbColor(1), rgbColor(2), rgbColor(3));
end
function writeTikzLinePlots(fid, linePlots)
for lineIdx = 1:numel(linePlots)
fprintf(fid, "\\addplot [color=mycolor%d, line width=%.1fpt, forget plot]\n", ...
lineIdx, linePlots(lineIdx).lineWidth);
fprintf(fid, " table[row sep=crcr]{%%\n");
for sampleIdx = 1:numel(linePlots(lineIdx).xData)
fprintf(fid, "%s\t%s\\\\\n", ...
formatPgfNumber(linePlots(lineIdx).xData(sampleIdx)), ...
formatPgfNumber(linePlots(lineIdx).yData(sampleIdx)));
end
fprintf(fid, "};\n\n");
end
end
function writeTikzTextAnnotations(fid, textAnnotations)
for textIdx = 1:numel(textAnnotations)
fprintf(fid, "\\node[font=\\scriptsize, anchor=%s, fill=white, draw=black!40, rounded corners=2pt, inner sep=2pt]\n", ...
textAnnotations(textIdx).anchor);
fprintf(fid, " at (axis cs:%s,%s) {%s};\n\n", ...
formatPgfNumber(textAnnotations(textIdx).x), ...
formatPgfNumber(textAnnotations(textIdx).y), ...
escapeTikzText(textAnnotations(textIdx).text));
end
end
function anchor = matlabTextAlignmentToTikzAnchor(horizontalAlignment, verticalAlignment)
horizontalAlignment = string(horizontalAlignment);
verticalAlignment = string(verticalAlignment);
if horizontalAlignment == "left"
horizontalAnchor = "west";
elseif horizontalAlignment == "right"
horizontalAnchor = "east";
else
horizontalAnchor = "";
end
if verticalAlignment == "top"
verticalAnchor = "north";
elseif verticalAlignment == "bottom"
verticalAnchor = "south";
else
verticalAnchor = "";
end
if strlength(horizontalAnchor) > 0 && strlength(verticalAnchor) > 0
anchor = verticalAnchor + " " + horizontalAnchor;
elseif strlength(horizontalAnchor) > 0
anchor = horizontalAnchor;
elseif strlength(verticalAnchor) > 0
anchor = verticalAnchor;
else
anchor = "center";
end
end
function textOut = escapeTikzText(textIn)
textOut = char(textIn);
% textOut = strrep(textOut, "\", "\textbackslash{}");
textOut = strrep(textOut, "{", "\{");
textOut = strrep(textOut, "}", "\}");
textOut = strrep(textOut, "%", "\%");
textOut = strrep(textOut, "&", "\&");
textOut = strrep(textOut, "#", "\#");
textOut = strrep(textOut, "_", "\_");
% textOut = strrep(textOut, "$", "\$");
end

Binary file not shown.

View File

@@ -133,7 +133,7 @@ for i = 1:numel(fsym_values)
eq_ = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",len_tr, ... eq_ = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",len_tr, ...
"mu_dd",1e-3,"mu_tr",0.4,"order",50, ... "mu_dd",1e-3,"mu_tr",0.4,"order",50, ...
"sps",2,"decide",0,"optmize_mus",0,"dd_mode",1, ... "sps",2,"decide",0,"optmize_mus",0,"dd_mode",1, ...
"adaption_technique","nlms","mu_dc",1e-3); "adaption_technique","nlms","dc_tracking_mu",1e-3);
pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1);
mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels);

View File

@@ -134,7 +134,7 @@ for i = 1:numel(fsym_values)
eq_ = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",len_tr, ... eq_ = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",len_tr, ...
"mu_dd",1e-1,"mu_tr",0.4,"order",50, ... "mu_dd",1e-1,"mu_tr",0.4,"order",50, ...
"sps",2,"decide",0,"optmize_mus",1,"dd_mode",1, ... "sps",2,"decide",0,"optmize_mus",1,"dd_mode",1, ...
"adaption_technique","nlms","mu_dc",1.021e-05); "adaption_technique","nlms","dc_tracking_mu",1.021e-05);
case "VNLE" case "VNLE"
eq_ = VNLE("epochs_tr",5,"epochs_dd",5,"len_tr",len_tr, ... eq_ = VNLE("epochs_tr",5,"epochs_dd",5,"len_tr",len_tr, ...
"mu_dd",[0.0004 0.0005 0.0006],"mu_tr",0, ... "mu_dd",[0.0004 0.0005 0.0006],"mu_tr",0, ...

View File

@@ -1,49 +1,10 @@
% dsp_options.database_type = 'mysql';
% dsp_options.dataBase = 'labor';
% dsp_options.storage_path = 'Z:\2025\ECOC Silas\ecoc_2025\';
% database = DBHandler("dataBase", [dsp_options.dataBase], "type", dsp_options.database_type);
% filterParams = database.tables;
% filterParams.Runs.loop_id = 209;
% % filterParams.Configurations = struct( ...
% % 'symbolrate', 112e9, ... %[224,336,360,390,420,448]
% % 'fiber_length', 0, ...
% % 'db_mode', '"no_db"', ...
% % 'interference_attenuation', [], ...
% % 'interference_path_length', 1000, ...
% % 'is_mpi', 1, ...
% % 'pam_level', 4, ...
% % 'wavelength', 1310, ...
% % 'precomp_amp', [], ...
% % 'signal_attenuation', [], ...
% % 'v_awg', [], ...
% % 'v_bias', [] ...
% % );
%
% % if 1
% % % filterParams.EqualizerParameters.dc_buffer_len = 1;
% % filterParams.EqualizerParameters.ffe_buffer_len = 1;
% % filterParams.EqualizerParameters.smoothing_buffer_len = 4096;
% % filterParams.EqualizerParameters.smoothing_buffer_update = 224;
% % filterParams.EqualizerParameters.DCmu = 0;
% % end
% a = database.getTableFieldNames('Runs');
% b = database.getTableFieldNames('Results');
% c = database.getTableFieldNames('EqualizerParameters');
% d = [a;b;c];
%
% [dataTable,~] = database.queryDB(filterParams, d);
%
% selectedFields = {'Configurations.run_id' 'Runs.loop_id' 'Runs.date_of_run' 'Runs.rx_raw_path' 'Runs.bitrate' 'Runs.v_bias' 'Runs.v_awg' 'Runs.precomp_amp' 'Runs.symbolrate' 'Runs.pam_level'...
% 'Runs.db_mode' 'Runs.rop_attenuation' 'Runs.is_mpi' 'Runs.interference_attenuation' 'Runs.interference_path_length' 'Runs.signal_attenuation' ...
% 'EqualizerParameters.equalizer_structure' 'EqualizerParameters.diff_precode' 'EqualizerParameters.eq_id' 'EqualizerParameters.dc_buffer_len' 'EqualizerParameters.ffe_buffer_len' 'EqualizerParameters.smoothing_buffer_len' 'EqualizerParameters.smoothing_buffer_update' 'EqualizerParameters.DCmu' 'Measurements.power_pd_in' ...
% 'Measurements.power_mpi_interference' 'Measurements.power_mpi_signal' 'Results.BER' 'Results.BER_precoded' 'Results.EVM' 'Results.SNR' 'Results.GMI' 'Results.Alpha' 'Results.date_of_processing'};
db = DBHandler("type","mysql","dataBase",'labor'); db = DBHandler("type","mysql","dataBase",'labor');
fp = QueryFilter(); fp = QueryFilter();
% fp.where('mpi_superview', 'loop_id','EQUALS', 209); % fp.where('mpi_superview', 'loop_id','EQUALS', 209);
% fp.where('mpi_superview','run_id','GREATER_EQUAL',3153);
fp.where('mpi_superview', 'symbolrate','EQUALS', 112e9); fp.where('mpi_superview', 'symbolrate','EQUALS', 112e9);
fp.where('mpi_superview', 'pam_level','EQUALS', 4); fp.where('mpi_superview', 'pam_level','EQUALS', 4);
fn = [db.getTableFieldNames('mpi_superview')]; fn = [db.getTableFieldNames('mpi_superview')];
@@ -67,7 +28,7 @@ figure()
tiledlayout(1, 1, 'TileSpacing', 'compact', 'Padding', 'compact'); tiledlayout(1, 1, 'TileSpacing', 'compact', 'Padding', 'compact');
y_here = 0; y_here = 0;
figcnt = 0; figcnt = 0;
for int_len = 1000%[0,50,300,1000] for int_len = 1%[0,50,300,1000]
figcnt = figcnt+1; figcnt = figcnt+1;
% figure(int_len+1); % figure(int_len+1);
nexttile; nexttile;
@@ -139,7 +100,7 @@ for int_len = 1000%[0,50,300,1000]
% Modify values in 'interference_path_length' where the condition is met % Modify values in 'interference_path_length' where the condition is met
dataTable.interference_path_length(dataTable.interference_path_length < 101 & dataTable.interference_path_length > 1) = 50; dataTable.interference_path_length(dataTable.interference_path_length < 101 & dataTable.interference_path_length > 1) = 50;
dataTable.interference_path_length(dataTable.interference_path_length == 1) = 0; % dataTable.interference_path_length(dataTable.interference_path_length == 1) = 0;
dataTable = dataTable(dataTable.interference_path_length == int_len, :); dataTable = dataTable(dataTable.interference_path_length == int_len, :);

Binary file not shown.