Files
imdd_silas/Classes/04_DSP/Equalizer/FFE.m
Silas Oettinghaus f421348e5b before new database
2026-07-13 19:59:57 +02:00

1465 lines
68 KiB
Matlab
Raw Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
classdef FFE < handle
% Implementation of plain and simple FFE.
% 1) Training mode (stable performance when you use NLMS)
% 2) Decision directed mode
%LMS: mu in order of 0.001 for acceptable convergence speed
%NLMS: mu in order of 0.01 for acceptable convergence speed
%RLS: mu is lambda -> 0.99 -> 1 (has a strong dependency on this! use a loop to find out best values)
% 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",adaption_method(adaption),"dd_mode",use_dd_mode);
% eq_ = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",len_tr, ...
% "mu_dd",1e-1,"mu_tr",0.4,"order",50, ...
% "sps",2,"decide",0,"optmize_mus",1,"dd_mode",1, ...
% "adaption_technique","nlms","dc_tracking_mu",1.021e-05);
properties
sps % usually 2
order
e
e_tr
error
len_tr
mu_tr
epochs_tr
adaption_technique % nlms, lms, rls
dd_mode % 1 or 0 to set DD-mode on or off
mu_dd %weight update in dd mode
epochs_dd
dd_len_fraction
% A1 moving-average input suppression
dc_avg_bufferlength_a1
dc_avg_update_blocklength_a1
dc_smoothing_a1
% A2 level-dependent residual suppression
dc_level_avg_bufferlength_a2
dc_level_update_blocklength_a2
dc_smoothing_a2
dc_level_decision_mode_a2
dc_level_weights_a2
% Adaptive DC-tracking loop
dc_tracking_mu
dc_tracking_adaptive_enabled
dc_tracking_persistence_gain
dc_tracking_power_exponent
dc_tracking_buffer_len
% Delayed FFE tap update
ffe_update_buffer_len
e_dc
P % covariance matrix of rls
constellation % symbol constellation
decide %wether to return the (hard) decisions or the result after FFE (soft)
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;
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
optimize_a2_level_weights = 0;
a2_level_weight_optimization
a2_level_weight_optimization_iter = 0;
a2_level_weight_optimization_len
a2_level_weight_optimization_max_evals
a2_level_weight_max
a2_level_weight_initial_stats
end
properties (Access = private)
dc_tracking_alpha = 0.98
dc_tracking_gamma = 1e-6
dc_tracking_mu_min = 1e-6
dc_tracking_mu_max = 3e-1
dc_tracking_mu_eff_min = 0
dc_tracking_mu_eff_max = inf
end
methods
function obj = FFE(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.dc_level_avg_bufferlength_a2 = 0;
options.dc_level_update_blocklength_a2 = 0;
options.dc_smoothing_a2 = 0;
options.dc_level_decision_mode_a2 = "residual";
options.dc_level_weights_a2 = 0;
options.dc_tracking_mu = 0;
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.ffe_update_buffer_len = 1;
options.decide = false;
options.save_debug = 0;
options.optmize_mus = 0;
options.mu_optimization_len = 2^15;
options.plot_mu_optimization = 0;
options.mu_optimization_fignum = 3010;
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;
options.optimize_a2_level_weights = 0;
options.a2_level_weight_optimization_len = 2^15;
options.a2_level_weight_optimization_max_evals = 30;
options.a2_level_weight_max = 1;
end
fn = fieldnames(options);
for n = 1:numel(fn)
obj.(fn{n}) = options.(fn{n});
end
assert(obj.dc_tracking_buffer_len >= 0);
assert(obj.ffe_update_buffer_len >= 0);
assert(obj.dc_avg_bufferlength_a1 >= 0);
assert(obj.dc_avg_update_blocklength_a1 >= 0);
assert(obj.dc_level_avg_bufferlength_a2 >= 0);
assert(obj.dc_level_update_blocklength_a2 >= 0);
obj.e = zeros(obj.order,1);
obj.e_dc = 0;
obj.error = 0;
obj.a2_level_weight_initial_stats = struct();
obj.dc_avg_bufferlength_a1 = floor(obj.dc_avg_bufferlength_a1);
obj.dc_avg_update_blocklength_a1 = floor(obj.dc_avg_update_blocklength_a1);
obj.dc_smoothing_a1 = min(max(obj.dc_smoothing_a1,0),1);
obj.dc_level_avg_bufferlength_a2 = floor(obj.dc_level_avg_bufferlength_a2);
obj.dc_level_update_blocklength_a2 = floor(obj.dc_level_update_blocklength_a2);
obj.dc_smoothing_a2 = min(max(obj.dc_smoothing_a2,0),1);
obj.dc_level_decision_mode_a2 = string(obj.dc_level_decision_mode_a2);
if ~any(obj.dc_level_decision_mode_a2 == ["residual", "tracked_levels"])
builtin("error","FFE:InvalidA2DecisionMode", ...
"dc_level_decision_mode_a2 must be 'residual' or 'tracked_levels'.");
end
obj.dc_tracking_buffer_len = floor(obj.dc_tracking_buffer_len);
obj.ffe_update_buffer_len = floor(obj.ffe_update_buffer_len);
obj.a2_level_weight_max = max(obj.a2_level_weight_max,0);
if obj.dc_avg_bufferlength_a1 > 1 && (obj.dc_tracking_mu ~= 0 || obj.optimize_dc_tracking_params)
warning("FFE:DCAvgWithDcTracking", ...
"A1 moving-average DC suppression and dc_tracking_mu/adaptive DC tracking are alternative DC suppression paths and will be combined.");
end
a2_requested = obj.dc_level_avg_bufferlength_a2 > 1 && ...
(any(obj.dc_level_weights_a2(:) ~= 0) || obj.optimize_a2_level_weights);
if a2_requested && ...
(obj.dc_tracking_mu ~= 0 || obj.optimize_dc_tracking_params || obj.dc_avg_bufferlength_a1 > 1)
warning("FFE:DCLevelAvgWithOtherDcSuppression", ...
"A2 level-dependent MPI suppression is enabled together with another DC/MPI suppression path; both corrections will be combined.");
end
end
function [X,Noi] = process(obj, X, D)
% actual processing of the signal (steps 1. - 3.)
% 1 normalize RMS
X = X.normalize("mode","rms");
obj.constellation = unique(D.signal);
obj.e_dc = 0;
delta = 0.05;
obj.P = (1/delta) * eye(obj.order);
if obj.optmize_mus
obj.optimizeMus(X.signal,D.signal);
obj.e = zeros(obj.order,1);
obj.e_dc = 0;
obj.P = (1/delta) * eye(obj.order);
end
if obj.optimize_dc_tracking_params
obj.optimizeDcTrackingParams(X.signal,D.signal);
obj.e = zeros(obj.order,1);
obj.e_dc = 0;
obj.P = (1/delta) * eye(obj.order);
end
if obj.optimize_a2_level_weights
obj.optimizeA2LevelWeights(X.signal,D.signal);
obj.e = zeros(obj.order,1);
obj.e_dc = 0;
obj.P = (1/delta) * eye(obj.order);
end
% Training Mode
training = 1;
showviz = 0;
obj.equalize(X.signal, D.signal,obj.mu_tr,obj.epochs_tr,obj.len_tr,training,showviz);
obj.e_tr = obj.e;
% Decision Directed Mode
n = X.length;
training = 0;
showviz = 0;
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
% Output Signal
if obj.decide
X.signal = decision;
else
X.signal = signal;
end
X.fs = D.fs; %change sampling frequency of outgoing signal from fdac e.g. 2 sps to symbol spaced = fsym
lbdesc = [num2str(obj.order),' tap FFE'];
X = X.logbookentry(lbdesc); % append to logbook
Noi = X;
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
x = [zeros(floor(obj.order/2),1); x; zeros(obj.order,1)];
lambda = mu;
n_symbols = ceil(N / obj.sps);
y = zeros(n_symbols,1);
d_hat = zeros(n_symbols,1);
err = zeros(n_symbols,1);
true_err = zeros(n_symbols,1);
constellation = obj.constellation;
if obj.adaption_technique == adaption_method.nlms
adaption_code = 1;
elseif obj.adaption_technique == adaption_method.lms
adaption_code = 2;
elseif obj.adaption_technique == adaption_method.rls
adaption_code = 3;
else
builtin("error","FFE:InvalidAdaptionTechnique", ...
"Unsupported FFE adaption technique.");
end
adaption_is_rls = adaption_code == 3;
if training
mask = ones(obj.order,1);
else
mask = zeros(obj.order,1);
mask(900:end) = 1;
mask(ceil(length(obj.e)/2)) = 1;
end
mask = ones(obj.order,1);
maincursor_pos=ceil(length(obj.e)/2);
always_ideal_decision = 0;
grad =0;
weight = 0;
update = 0;
if mu == 0 || (~obj.dd_mode && ~training)
epochs = 1;
end
dc_tracking_mu_min = obj.dc_tracking_mu_min;
dc_tracking_mu_max = obj.dc_tracking_mu_max;
dc_tracking_mu_eff_min = obj.dc_tracking_mu_eff_min;
dc_tracking_mu_eff_max = obj.dc_tracking_mu_eff_max;
dc_tracking_persistence_gain = 0;
if obj.dc_tracking_adaptive_enabled
dc_tracking_persistence_gain = obj.dc_tracking_persistence_gain;
end
dc_tracking_enabled = obj.dc_tracking_mu ~= 0;
if dc_tracking_enabled
obj.dc_tracking_mu = min(max(obj.dc_tracking_mu,dc_tracking_mu_min),dc_tracking_mu_max);
end
dc_tracking_mu = obj.dc_tracking_mu;
dc_tracking_use_persistence = dc_tracking_persistence_gain > 0;
dc_tracking_base_mu_eff = min(max(dc_tracking_mu,dc_tracking_mu_eff_min),dc_tracking_mu_eff_max);
dc_buffer_enabled = dc_tracking_enabled && obj.dc_tracking_buffer_len > 1;
if dc_buffer_enabled
dc_tracking_err_buffer = NaN(obj.dc_tracking_buffer_len,1);
dc_tracking_err_buffer_pos = 0;
dc_tracking_err_sum = 0;
dc_tracking_abs_err_sum = 0;
dc_tracking_valid_count = 0;
end
ffe_buffer_enabled = obj.ffe_update_buffer_len > 1 && ~adaption_is_rls;
if ffe_buffer_enabled
ffe_update_buffer = NaN(obj.order,obj.ffe_update_buffer_len);
end
dc_avg_enabled = obj.dc_avg_bufferlength_a1 > 1;
if dc_avg_enabled
dc_avg_buffer_a1 = NaN(obj.dc_avg_bufferlength_a1,1);
dc_avg_est_a1 = 0;
dc_avg_update_blocklength_a1 = obj.dc_avg_update_blocklength_a1;
if dc_avg_update_blocklength_a1 <= 0
dc_avg_update_blocklength_a1 = obj.dc_avg_bufferlength_a1;
end
dc_avg_update_blocklength_a1 = max(1,floor(dc_avg_update_blocklength_a1));
dc_avg_input_offset_a1 = floor(obj.order/2);
% Hardware-like A1 is causal; dc_smoothing_a1 is reserved for offline variants.
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:MissingConstellation", ...
"A2 level-dependent MPI suppression requires obj.constellation to be set.");
end
dc_level_decision_mode_a2 = obj.dc_level_decision_mode_a2;
dc_level_track_decision_a2 = dc_level_decision_mode_a2 == "tracked_levels";
n_levels = numel(constellation);
if isscalar(obj.dc_level_weights_a2)
dc_level_weight_by_level = repmat(obj.dc_level_weights_a2,n_levels,1);
elseif numel(obj.dc_level_weights_a2) == n_levels
dc_level_weight_by_level = obj.dc_level_weights_a2(:);
else
builtin("error","FFE:InvalidDCLevelWeights", ...
"dc_level_weights_a2 must be scalar or have one entry per constellation level.");
end
dc_level_buffer_len_a2 = obj.dc_level_avg_bufferlength_a2;
dc_level_err_buffer = NaN(n_levels,dc_level_buffer_len_a2);
dc_level_err_buffer_pos_by_level = zeros(n_levels,1);
dc_level_err_sum_by_level = zeros(n_levels,1);
dc_level_buffer_valid_count_by_level = zeros(n_levels,1);
dc_level_mpi_est_by_level = zeros(n_levels,1);
dc_level_valid_count_by_level = zeros(n_levels,1);
dc_level_update_blocklength_a2 = obj.dc_level_update_blocklength_a2;
if dc_level_update_blocklength_a2 <= 0
dc_level_update_blocklength_a2 = dc_level_buffer_len_a2;
end
dc_level_update_blocklength_a2 = max(1,floor(dc_level_update_blocklength_a2));
dc_level_window_future_fraction = obj.dc_smoothing_a2; %#ok<NASGU> % reserved for delayed/offline A2 variants
end
debug_enabled = obj.save_debug;
if debug_enabled
n_symbols_debug = n_symbols;
obj.debug_struct.error = NaN(1,n_symbols_debug);
obj.debug_struct.error_first_epoch = NaN(1,n_symbols_debug);
obj.debug_struct.main_cursor = NaN(1,n_symbols_debug);
obj.debug_struct.mu_nlms = NaN(1,n_symbols_debug);
obj.debug_struct.update_gradient = NaN(1,n_symbols_debug);
obj.debug_struct.dc_tracking_mu_eff = NaN(1,n_symbols_debug);
obj.debug_struct.dc_tracking_est = NaN(1,n_symbols_debug);
obj.debug_struct.dc_avg_offset = NaN(1,n_symbols_debug);
obj.debug_struct.dc_level_mpi_est = NaN(1,n_symbols_debug);
obj.debug_struct.dc_level_weight = NaN(1,n_symbols_debug);
obj.debug_struct.dc_level_valid_count = NaN(1,n_symbols_debug);
obj.debug_struct.dc_level_symbol_idx = NaN(1,n_symbols_debug);
obj.debug_struct.dc_level_decision_level = NaN(1,n_symbols_debug);
obj.debug_struct.dc_level_y_raw = NaN(1,n_symbols_debug);
if training
obj.debug_struct.error_tr = NaN(1,n_symbols_debug);
obj.debug_struct.update_tr = NaN(1,n_symbols_debug);
else
obj.debug_struct.error_dd = NaN(1,n_symbols_debug);
obj.debug_struct.update = NaN(1,n_symbols_debug);
end
end
for epoch = 1 : epochs
symbol = 0;
% obj.e_dc = 0;
for sample = 1 : obj.sps : N
symbol = symbol+1;
dc_tracking_mu_eff = 0;
dc_avg_offset = 0;
dc_level_mpi_est = 0;
dc_level_weight = 0;
dc_level_valid_count = 0;
dc_level_symbol_idx = NaN;
dc_level_decision_level = NaN;
U = x(obj.order+sample-1:-1:sample);
if dc_avg_enabled
dc_avg_offset = dc_avg_est_a1;
U = U - dc_avg_offset;
end
y(symbol,1) = obj.e_dc + (obj.e.*mask).' * U; % Calculating output of LMS __ * |
y_raw = y(symbol);
if dc_avg_enabled
dc_avg_input_idx = dc_avg_input_offset_a1 + sample;
dc_avg_buffer_a1 = circshift(dc_avg_buffer_a1,1);
dc_avg_buffer_a1(1) = x(dc_avg_input_idx);
if mod(symbol,dc_avg_update_blocklength_a1) == 0
dc_avg_est_a1 = mean(dc_avg_buffer_a1,"omitnan");
end
end
if training
d_hat(symbol,1) = d(symbol);
if isempty(constellation)
symbol_idx = NaN;
else
[~,symbol_idx] = min(abs(d_hat(symbol) - constellation));
dc_level_decision_level = constellation(symbol_idx);
end
else
if ~always_ideal_decision
if dc_level_enabled && dc_level_track_decision_a2
% Tracked-level A2 moves the decision constellation; y itself stays unchanged.
dc_level_ramp_by_level = min(dc_level_valid_count_by_level / dc_level_buffer_len_a2,1);
dc_level_decision_constellation = constellation + ...
dc_level_weight_by_level .* dc_level_ramp_by_level .* dc_level_mpi_est_by_level;
[~,symbol_idx] = min(abs(y(symbol) - dc_level_decision_constellation));
dc_level_decision_level = dc_level_decision_constellation(symbol_idx);
else
[~,symbol_idx] = min(abs(y(symbol) - constellation)); % decision for closest constellation point
dc_level_decision_level = constellation(symbol_idx);
end
d_hat(symbol,1) = constellation(symbol_idx);
else
d_hat(symbol,1) = d(symbol);
[~,symbol_idx] = min(abs(d_hat(symbol) - constellation));
dc_level_decision_level = constellation(symbol_idx);
end
end
dc_level_symbol_idx = symbol_idx;
if dc_level_enabled
dc_level_mpi_est = dc_level_mpi_est_by_level(symbol_idx);
dc_level_valid_count = dc_level_valid_count_by_level(symbol_idx);
dc_level_weight = dc_level_weight_by_level(symbol_idx) * ...
min(dc_level_valid_count / dc_level_buffer_len_a2,1);
if dc_level_track_decision_a2
dc_level_buffer_value = y_raw;
y(symbol,1) = y_raw;
if isnan(dc_level_decision_level)
dc_level_decision_level = constellation(symbol_idx) + ...
dc_level_weight * dc_level_mpi_est;
end
else
mpi_err = y_raw - d_hat(symbol);
dc_level_buffer_value = mpi_err;
y(symbol,1) = y_raw - dc_level_weight * dc_level_mpi_est;
if training
d_hat(symbol,1) = d(symbol);
else
if ~always_ideal_decision
[~,symbol_idx] = min(abs(y(symbol) - constellation));
d_hat(symbol,1) = constellation(symbol_idx);
dc_level_decision_level = constellation(symbol_idx);
else
d_hat(symbol,1) = d(symbol);
[~,symbol_idx] = min(abs(d_hat(symbol) - constellation));
dc_level_decision_level = constellation(symbol_idx);
end
end
end
dc_level_err_buffer_pos = dc_level_err_buffer_pos_by_level(dc_level_symbol_idx) + 1;
if dc_level_err_buffer_pos > dc_level_buffer_len_a2
dc_level_err_buffer_pos = 1;
end
dc_level_err_buffer_pos_by_level(dc_level_symbol_idx) = dc_level_err_buffer_pos;
dc_level_old_err = dc_level_err_buffer(dc_level_symbol_idx,dc_level_err_buffer_pos);
if isfinite(dc_level_old_err)
dc_level_err_sum_by_level(dc_level_symbol_idx) = ...
dc_level_err_sum_by_level(dc_level_symbol_idx) - dc_level_old_err;
dc_level_buffer_valid_count_by_level(dc_level_symbol_idx) = ...
dc_level_buffer_valid_count_by_level(dc_level_symbol_idx) - 1;
end
if isfinite(dc_level_buffer_value)
dc_level_err_buffer(dc_level_symbol_idx,dc_level_err_buffer_pos) = dc_level_buffer_value;
dc_level_err_sum_by_level(dc_level_symbol_idx) = ...
dc_level_err_sum_by_level(dc_level_symbol_idx) + dc_level_buffer_value;
dc_level_buffer_valid_count_by_level(dc_level_symbol_idx) = ...
dc_level_buffer_valid_count_by_level(dc_level_symbol_idx) + 1;
else
dc_level_err_buffer(dc_level_symbol_idx,dc_level_err_buffer_pos) = NaN;
end
if mod(symbol,dc_level_update_blocklength_a2) == 0
if dc_level_update_blocklength_a2 == 1
dc_level_valid_count_by_level(dc_level_symbol_idx) = ...
dc_level_buffer_valid_count_by_level(dc_level_symbol_idx);
if dc_level_valid_count_by_level(dc_level_symbol_idx) == 0
dc_level_mpi_est_by_level(dc_level_symbol_idx) = 0;
else
dc_level_mpi_est_by_level(dc_level_symbol_idx) = ...
dc_level_err_sum_by_level(dc_level_symbol_idx) / ...
dc_level_valid_count_by_level(dc_level_symbol_idx);
if dc_level_track_decision_a2
dc_level_mpi_est_by_level(dc_level_symbol_idx) = ...
dc_level_mpi_est_by_level(dc_level_symbol_idx) - ...
constellation(dc_level_symbol_idx);
end
end
else
dc_level_valid_count_by_level = dc_level_buffer_valid_count_by_level;
dc_level_has_valid = dc_level_valid_count_by_level > 0;
dc_level_mpi_est_by_level(:) = 0;
dc_level_mpi_est_by_level(dc_level_has_valid) = ...
dc_level_err_sum_by_level(dc_level_has_valid) ./ ...
dc_level_valid_count_by_level(dc_level_has_valid);
if dc_level_track_decision_a2
dc_level_mpi_est_by_level(dc_level_has_valid) = ...
dc_level_mpi_est_by_level(dc_level_has_valid) - ...
constellation(dc_level_has_valid);
end
end
end
end
% err(symbol) = y(symbol) - d_hat(symbol); % Instantaneous error
err(symbol) = d_hat(symbol) - y(symbol); % Instantaneous error
true_err(symbol) = y(symbol) - d(symbol); % Instantaneous error
if 1 %training || obj.dd_mode
switch adaption_code
case 1
% mu used as update weight (suggestion: 0.01-0.05; bit higher during tr)
normU = ((U.'*U)) + eps;
weight = mu / normU;
grad = err(symbol) * U;
update = grad * weight;
case 2
% mu used as update weight (suggestion: 0.001)
weight = mu;
grad = err(symbol) * U;
update = grad * weight;
case 3
% RLSGain:
denom = lambda + U.' * obj.P * U;
k = (obj.P * U) / denom;
% Gewichtsupdate:
update = k * err(symbol);
% P-MatrixUpdate:
end
if ffe_buffer_enabled
ffe_update_buffer = circshift(ffe_update_buffer,1,2);
ffe_update_buffer(:,1) = update;
if mod(symbol,obj.ffe_update_buffer_len) == 0
obj.e = obj.e + mean(ffe_update_buffer,2,"omitnan");
end
else
obj.e = obj.e + update;
end
if adaption_is_rls
obj.P = (1/lambda) * (obj.P - k * (U.' * obj.P));
end
if dc_tracking_enabled
if dc_buffer_enabled
dc_tracking_err_buffer_pos = dc_tracking_err_buffer_pos + 1;
if dc_tracking_err_buffer_pos > obj.dc_tracking_buffer_len
dc_tracking_err_buffer_pos = 1;
end
dc_tracking_old_err = dc_tracking_err_buffer(dc_tracking_err_buffer_pos);
if isfinite(dc_tracking_old_err)
dc_tracking_err_sum = dc_tracking_err_sum - dc_tracking_old_err;
dc_tracking_valid_count = dc_tracking_valid_count - 1;
if dc_tracking_use_persistence
dc_tracking_abs_err_sum = dc_tracking_abs_err_sum - abs(dc_tracking_old_err);
end
end
dc_tracking_new_err = err(symbol);
if isfinite(dc_tracking_new_err)
dc_tracking_err_buffer(dc_tracking_err_buffer_pos) = dc_tracking_new_err;
dc_tracking_err_sum = dc_tracking_err_sum + dc_tracking_new_err;
dc_tracking_valid_count = dc_tracking_valid_count + 1;
if dc_tracking_use_persistence
dc_tracking_abs_err_sum = dc_tracking_abs_err_sum + abs(dc_tracking_new_err);
end
else
dc_tracking_err_buffer(dc_tracking_err_buffer_pos) = NaN;
end
if mod(symbol,obj.dc_tracking_buffer_len) == 0
if dc_tracking_valid_count == 0
dc_tracking_err_mean = 0;
if dc_tracking_use_persistence
dc_tracking_err_abs_mean = 0;
end
else
dc_tracking_err_mean = dc_tracking_err_sum / dc_tracking_valid_count;
if dc_tracking_use_persistence
dc_tracking_err_abs_mean = dc_tracking_abs_err_sum / dc_tracking_valid_count;
end
end
if dc_tracking_use_persistence
dc_tracking_persistence_scale = abs(dc_tracking_err_mean) / (dc_tracking_err_abs_mean + eps);
dc_tracking_persistence_scale = min(max(dc_tracking_persistence_scale,0),1);
dc_tracking_mu_eff = dc_tracking_mu * ...
(1 + dc_tracking_persistence_gain * dc_tracking_persistence_scale);
dc_tracking_mu_eff = min(max(dc_tracking_mu_eff,dc_tracking_mu_eff_min),dc_tracking_mu_eff_max);
else
dc_tracking_mu_eff = dc_tracking_base_mu_eff;
end
obj.e_dc = obj.e_dc + dc_tracking_mu_eff * dc_tracking_err_mean;
end
else
dc_tracking_err_mean = err(symbol);
if isfinite(dc_tracking_err_mean)
if dc_tracking_use_persistence
dc_tracking_err_abs_mean = abs(dc_tracking_err_mean);
dc_tracking_persistence_scale = abs(dc_tracking_err_mean) / ...
(dc_tracking_err_abs_mean + eps);
else
dc_tracking_mu_eff = dc_tracking_base_mu_eff;
end
else
dc_tracking_err_mean = 0;
if dc_tracking_use_persistence
dc_tracking_persistence_scale = 0;
else
dc_tracking_mu_eff = dc_tracking_base_mu_eff;
end
end
if dc_tracking_use_persistence
dc_tracking_persistence_scale = min(max(dc_tracking_persistence_scale,0),1);
dc_tracking_mu_eff = dc_tracking_mu * ...
(1 + dc_tracking_persistence_gain * dc_tracking_persistence_scale);
dc_tracking_mu_eff = min(max(dc_tracking_mu_eff,dc_tracking_mu_eff_min),dc_tracking_mu_eff_max);
end
obj.e_dc = obj.e_dc + dc_tracking_mu_eff * dc_tracking_err_mean;
end
end
end
%
% 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_tracking_mu_eff(1,symbol) = dc_tracking_mu_eff;
obj.debug_struct.dc_tracking_est(1,symbol) = obj.e_dc;
obj.debug_struct.dc_avg_offset(1,symbol) = dc_avg_offset;
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
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];
end
dc_tracking_mu_range = [1e-5, 1e-1];
mu_tr_var = optimizableVariable("mu_tr",mu_range,"Transform","log");
vars = mu_tr_var;
if obj.dd_mode
vars = [vars, optimizableVariable("mu_dd",mu_range,"Transform","log")];
end
optimize_dc_tracking_mu = obj.dc_tracking_mu ~= 0;
if optimize_dc_tracking_mu
vars = [vars, optimizableVariable("dc_tracking_mu",dc_tracking_mu_range,"Transform","log")];
end
obj.mu_optimization_iter = 0;
fprintf("FFE 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 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 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 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 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");
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")];
end
obj.dc_tracking_optimization_iter = 0;
fprintf("FFE DC opt uses fixed mu_tr=%9.3e, mu_dd=%9.3e on %d samples / %d symbols\n", ...
obj.mu_tr,obj.mu_dd,N_opt,numel(d_opt));
obj.dc_tracking_optimization = bayesopt(@(p)obj.dcTrackingObjective(p,x_opt,d_opt),vars, ...
"MaxObjectiveEvaluations",obj.dc_tracking_optimization_max_evals, ...
"AcquisitionFunctionName","expected-improvement-plus", ...
"IsObjectiveDeterministic",false, ...
"Verbose",0, ...
"PlotFcn",[]);
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 DC opt done: dc_tracking_mu=%9.3e, persistence_gain=%6.3f, mu_eff_max=%9.3e, objective=%9.3e\n", ...
obj.dc_tracking_mu,obj.dc_tracking_persistence_gain,obj.dc_tracking_mu_eff_max,obj.dc_tracking_optimization.MinObjective);
else
fprintf("\nFFE DC opt done: dc_tracking_mu=%9.3e, objective=%9.3e\n", ...
obj.dc_tracking_mu,obj.dc_tracking_optimization.MinObjective);
end
end
function optimizeA2LevelWeights(obj,x,d)
if obj.dc_level_avg_bufferlength_a2 <= 1
warning("FFE:A2OptimizationDisabled", ...
"A2 level-weight optimization requires dc_level_avg_bufferlength_a2 > 1.");
return
end
if obj.a2_level_weight_max <= 0
warning("FFE:A2OptimizationDisabled", ...
"A2 level-weight optimization requires a2_level_weight_max > 0.");
return
end
if isempty(obj.constellation)
obj.constellation = unique(d);
end
[x_opt,d_opt,N_opt] = obj.optimizationSignals(x,d,obj.a2_level_weight_optimization_len);
n_levels = numel(obj.constellation);
var_names = compose("dc_level_weight_a2_%d",1:n_levels);
vars_cell = cell(1,n_levels);
for level_idx = 1:n_levels
vars_cell{level_idx} = optimizableVariable(var_names(level_idx),[0,obj.a2_level_weight_max]);
end
vars = [vars_cell{:}];
[initial_weights,stats] = obj.a2LevelWeightInitialGuess(x_opt,d_opt);
current_weights = obj.expandA2LevelWeights(n_levels);
initial_matrix = [zeros(1,n_levels); initial_weights(:).'];
if any(current_weights ~= 0)
initial_matrix = [initial_matrix; current_weights(:).'];
end
initial_matrix = min(max(initial_matrix,0),obj.a2_level_weight_max);
initial_matrix = unique(initial_matrix,"rows","stable");
initial_x = array2table(initial_matrix,"VariableNames",cellstr(var_names));
[baseline_ber,baseline_errors] = obj.a2LevelWeightBer(initial_x(1,:),x_opt,d_opt);
[x_val,d_val,N_val] = obj.a2ValidationSignals(x,d,N_opt);
stats.baseline_ber = baseline_ber;
stats.baseline_errors = baseline_errors;
obj.a2_level_weight_initial_stats = stats;
obj.a2_level_weight_optimization_iter = 0;
max_evals = max(obj.a2_level_weight_optimization_max_evals,height(initial_x));
fprintf("FFE A2 opt uses fixed mu_tr=%9.3e, mu_dd=%9.3e on %d samples / %d symbols\n", ...
obj.mu_tr,obj.mu_dd,N_opt,numel(d_opt));
fprintf("FFE A2 opt validation uses %d samples / %d symbols\n",N_val,numel(d_val));
fprintf("FFE A2 opt init: baseline BER=%9.3e (%d errors), var_slope=%9.3e, weights=%s\n", ...
baseline_ber,baseline_errors,stats.variance_slope,mat2str(initial_weights(:).',3));
old_rng = rng;
cleanup_rng = onCleanup(@()rng(old_rng));
rng(42,"twister");
obj.a2_level_weight_optimization = bayesopt(@(p)obj.a2LevelWeightObjective(p,x_opt,d_opt,baseline_ber),vars, ...
"MaxObjectiveEvaluations",max_evals, ...
"InitialX",initial_x, ...
"AcquisitionFunctionName","expected-improvement-plus", ...
"IsObjectiveDeterministic",true, ...
"Verbose",0, ...
"PlotFcn",[]);
clear cleanup_rng
[best,best_validation_ber,best_validation_errors] = obj.selectA2LevelWeightsByValidation( ...
obj.a2_level_weight_optimization.XTrace, ...
obj.a2_level_weight_optimization.ObjectiveTrace, ...
x_val,d_val,initial_x);
best_opt = obj.a2_level_weight_optimization.XAtMinObjective;
best_opt_weights = obj.a2LevelWeightsFromParams(best_opt);
obj.dc_level_weights_a2 = obj.a2LevelWeightsFromParams(best);
obj.a2_level_weight_initial_stats.validation_ber = best_validation_ber;
obj.a2_level_weight_initial_stats.validation_errors = best_validation_errors;
obj.a2_level_weight_initial_stats.validation_weights = obj.dc_level_weights_a2;
fprintf("\nFFE A2 opt done: opt_weights=%s, opt_obj=%9.3e, validation_weights=%s, validation_BER=%9.3e (%d errors)\n", ...
mat2str(best_opt_weights(:).',3),obj.a2_level_weight_optimization.MinObjective, ...
mat2str(obj.dc_level_weights_a2(:).',3),best_validation_ber,best_validation_errors);
end
function [x_opt,d_opt,N_opt] = optimizationSignals(obj,x,d,opt_len)
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 [x_val,d_val,N_val] = a2ValidationSignals(obj,x,d,N_opt)
N_available = min(numel(x),numel(d) * obj.sps);
N_val = min(N_opt,N_available);
if N_available <= N_opt
[x_val,d_val,N_val] = obj.optimizationSignals(x,d,N_opt);
return
end
start_symbol = floor((N_available - N_val) / obj.sps) + 1;
start_sample = (start_symbol - 1) * obj.sps + 1;
N_val = obj.sps * floor((N_available - start_sample + 1) / obj.sps);
N_val = max(obj.sps,N_val);
n_symbols = N_val / obj.sps;
x_val = x(start_sample:start_sample+N_val-1);
d_val = d(start_symbol:start_symbol+n_symbols-1);
end
function [best_params,best_ber,best_errors] = selectA2LevelWeightsByValidation(obj,x_trace,objective_trace,x_val,d_val,initial_x)
objective_trace = objective_trace(:);
objective_trace(~isfinite(objective_trace)) = inf;
[~,sort_idx] = sort(objective_trace,"ascend");
n_trace_candidates = min(8,numel(sort_idx));
candidate_x = x_trace(sort_idx(1:n_trace_candidates),:);
candidate_x = [initial_x; candidate_x];
candidate_x = unique(candidate_x,"rows","stable");
n_candidates = height(candidate_x);
validation_ber = inf(n_candidates,1);
validation_errors = nan(n_candidates,1);
for candidate_idx = 1:n_candidates
[validation_ber(candidate_idx),validation_errors(candidate_idx)] = ...
obj.a2LevelWeightBer(candidate_x(candidate_idx,:),x_val,d_val);
end
[best_ber,best_idx] = min(validation_ber);
best_errors = validation_errors(best_idx);
best_params = candidate_x(best_idx,:);
fprintf("FFE A2 validation: checked %d candidates, best weights=%s, BER=%9.3e, errors=%d\n", ...
n_candidates,mat2str(obj.a2LevelWeightsFromParams(best_params).',3),best_ber,best_errors);
end
function objective = muObjective(obj,params,x,d)
old_debug = obj.save_debug;
old_dc_tracking_mu = obj.dc_tracking_mu;
obj.save_debug = 0;
has_dc_tracking_mu = any(strcmp(params.Properties.VariableNames,"dc_tracking_mu"));
if has_dc_tracking_mu
obj.dc_tracking_mu = params.dc_tracking_mu;
end
obj.e = zeros(obj.order,1);
obj.e_dc = 0;
obj.P = (1/0.05) * eye(obj.order);
obj.debug_struct = struct();
N_tr = min(obj.len_tr,numel(x));
obj.equalize(x,d,params.mu_tr,obj.epochs_tr,N_tr,1,0);
if obj.dd_mode
[signal,~] = obj.equalize(x,d,params.mu_dd,obj.epochs_dd,numel(x),0,0);
else
[signal,~] = obj.equalize(x,d,0,1,numel(x),0,0);
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 && has_dc_tracking_mu
fprintf("\rFFE 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 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 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 mu opt %02d: mu_tr=%9.3e, BER=%9.3e, errors=%d", ...
obj.mu_optimization_iter,params.mu_tr,ber,errors);
end
obj.save_debug = old_debug;
obj.dc_tracking_mu = old_dc_tracking_mu;
end
function objective = a2LevelWeightObjective(obj,params,x,d,baseline_ber)
if nargin < 5 || ~isfinite(baseline_ber)
baseline_ber = inf;
end
[ber,errors] = obj.a2LevelWeightBer(params,x,d);
objective = ber;
if isfinite(baseline_ber)
objective = objective + max(0,ber - baseline_ber);
end
if ~isfinite(objective)
objective = inf;
end
obj.a2_level_weight_optimization_iter = obj.a2_level_weight_optimization_iter + 1;
weights = obj.a2LevelWeightsFromParams(params);
fprintf("\rFFE A2 opt %02d: weights=%s, BER=%9.3e, obj=%9.3e, errors=%d", ...
obj.a2_level_weight_optimization_iter,mat2str(weights(:).',3),ber,objective,errors);
end
function [ber,errors] = a2LevelWeightBer(obj,params,x,d)
state = obj.captureObjectiveState();
cleanup = onCleanup(@()obj.restoreObjectiveState(state));
obj.dc_level_weights_a2 = obj.a2LevelWeightsFromParams(params);
obj.save_debug = 0;
obj.e = zeros(obj.order,1);
obj.e_dc = 0;
obj.P = (1/0.05) * eye(obj.order);
obj.debug_struct = struct();
N_tr = min(obj.len_tr,numel(x));
obj.equalize(x,d,obj.mu_tr,obj.epochs_tr,N_tr,1,0);
if obj.dd_mode
[signal,~] = obj.equalize(x,d,obj.mu_dd,obj.epochs_dd,numel(x),0,0);
else
[signal,~] = obj.equalize(x,d,0,1,numel(x),0,0);
end
[ber,errors] = obj.berObjective(signal,d);
end
function objective = dcTrackingObjective(obj,params,x,d)
state = obj.captureObjectiveState();
cleanup = onCleanup(@()obj.restoreObjectiveState(state));
obj.applyDcTrackingObjectiveParams(params);
obj.save_debug = 1;
obj.e = zeros(obj.order,1);
obj.e_dc = 0;
obj.P = (1/0.05) * eye(obj.order);
obj.debug_struct = struct();
N_tr = min(obj.len_tr,numel(x));
obj.equalize(x,d,obj.mu_tr,obj.epochs_tr,N_tr,1,0);
if obj.dd_mode
[signal,~] = obj.equalize(x,d,obj.mu_dd,obj.epochs_dd,numel(x),0,0);
else
[signal,~] = obj.equalize(x,d,0,1,numel(x),0,0);
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 DC opt %02d: dc_tracking_mu=%9.3e, persistence_gain=%6.3f, mu_eff_max=%9.3e, BER=%9.3e, delay=%7.0f, corr=%6.3f, obj=%9.3e, errors=%d", ...
obj.dc_tracking_optimization_iter,params.dc_tracking_mu,params.dc_tracking_persistence_gain,params.dc_tracking_mu_eff_max, ...
ber,delay_symbols,delay_corr,objective,errors);
else
fprintf("\rFFE DC opt %02d: dc_tracking_mu=%9.3e, BER=%9.3e, delay=%7.0f, corr=%6.3f, obj=%9.3e, errors=%d", ...
obj.dc_tracking_optimization_iter,params.dc_tracking_mu,ber,delay_symbols,delay_corr,objective,errors);
end
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 weights = a2LevelWeightsFromParams(~,params)
var_names = string(params.Properties.VariableNames);
weight_names = var_names(startsWith(var_names,"dc_level_weight_a2_"));
weight_idx = extractAfter(weight_names,"dc_level_weight_a2_");
weight_idx = str2double(weight_idx);
[weight_idx,sort_idx] = sort(weight_idx);
weight_names = weight_names(sort_idx);
weights = zeros(numel(weight_names),1);
for n = 1:numel(weight_names)
weights(weight_idx(n),1) = params.(char(weight_names(n)));
end
end
function weights = expandA2LevelWeights(obj,n_levels)
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:InvalidDCLevelWeights", ...
"dc_level_weights_a2 must be scalar or have one entry per constellation level.");
end
weights = min(max(weights,0),obj.a2_level_weight_max);
end
function [weights,stats] = a2LevelWeightInitialGuess(obj,x,d)
levels = obj.constellation(:);
if isempty(levels)
levels = unique(d(:));
end
n_levels = numel(levels);
n_symbols = min(numel(d),floor(numel(x) / obj.sps));
d = d(1:n_symbols);
rx_symbols = x(1:obj.sps:obj.sps*n_symbols);
rx_symbols = rx_symbols(:);
level_mean = NaN(n_levels,1);
level_variance = NaN(n_levels,1);
for level_idx = 1:n_levels
level_samples = rx_symbols(d == levels(level_idx));
level_mean(level_idx) = mean(level_samples,"omitnan");
level_variance(level_idx) = var(level_samples,0,"omitnan");
end
valid = isfinite(level_mean) & isfinite(level_variance);
variance_fit = level_variance;
slope_axis = "mean";
if nnz(valid) >= 2
if numel(unique(level_mean(valid))) >= 2
p = polyfit(level_mean(valid),level_variance(valid),1);
variance_fit = polyval(p,level_mean);
else
slope_axis = "level";
level_axis = (1:n_levels).';
p = polyfit(level_axis(valid),level_variance(valid),1);
variance_fit = polyval(p,(1:n_levels).');
end
else
p = [0, mean(level_variance(valid),"omitnan")];
end
if any(valid)
baseline = min(variance_fit(valid),[],"omitnan");
else
baseline = 0;
end
variance_score = max(variance_fit - baseline,0);
finite_score = isfinite(variance_score);
if ~any(finite_score) || max(variance_score(finite_score),[],"omitnan") <= 0
if any(valid)
baseline = min(level_variance(valid),[],"omitnan");
else
baseline = 0;
end
variance_score = max(level_variance - baseline,0);
finite_score = isfinite(variance_score);
end
weights = zeros(n_levels,1);
usable_score = valid & finite_score;
if any(usable_score)
max_score = max(variance_score(usable_score),[],"omitnan");
else
max_score = 0;
end
if isfinite(max_score) && max_score > 0
initial_weight_ceiling = min(0.7,obj.a2_level_weight_max);
weights(usable_score) = initial_weight_ceiling * variance_score(usable_score) ./ max_score;
end
weights = min(max(weights,0),obj.a2_level_weight_max);
stats = struct( ...
"levels",levels, ...
"mean",level_mean, ...
"variance",level_variance, ...
"variance_fit",variance_fit(:), ...
"variance_slope",p(1), ...
"slope_axis",slope_axis, ...
"initial_weights",weights);
end
function [e_dc_next,stats] = dcTrackingBlockUpdate(obj,e_dc_current,err_block,options)
arguments
obj
e_dc_current (1,1) double
err_block (:,1) double
options.mu_dc (1,1) double = NaN
options.persistence_gain (1,1) double = 0
options.mu_eff_min (1,1) double = NaN
options.mu_eff_max (1,1) double = NaN
end
if isnan(options.mu_dc)
options.mu_dc = obj.dc_tracking_mu;
end
if isnan(options.mu_eff_min)
options.mu_eff_min = obj.dc_tracking_mu_eff_min;
end
if isnan(options.mu_eff_max)
options.mu_eff_max = obj.dc_tracking_mu_eff_max;
end
valid_err = err_block(isfinite(err_block));
if isempty(valid_err)
err_mean = 0;
err_abs_mean = 0;
else
err_mean = mean(valid_err,"omitnan");
err_abs_mean = mean(abs(valid_err),"omitnan");
end
persistence_scale = abs(err_mean) / (err_abs_mean + eps);
persistence_scale = min(max(persistence_scale,0),1);
mu_eff = options.mu_dc * (1 + max(options.persistence_gain,0) * persistence_scale);
mu_eff = min(max(mu_eff,options.mu_eff_min),options.mu_eff_max);
update = mu_eff * err_mean;
e_dc_next = e_dc_current + update;
stats = struct( ...
"err_mean",err_mean, ...
"err_abs_mean",err_abs_mean, ...
"persistence_scale",persistence_scale, ...
"mu_eff",mu_eff, ...
"update",update);
end
function state = captureObjectiveState(obj)
state.e = obj.e;
state.e_dc = obj.e_dc;
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_alpha = obj.dc_tracking_alpha;
state.dc_tracking_gamma = obj.dc_tracking_gamma;
state.dc_tracking_persistence_gain = obj.dc_tracking_persistence_gain;
state.dc_tracking_power_exponent = obj.dc_tracking_power_exponent;
state.dc_tracking_mu_eff_max = obj.dc_tracking_mu_eff_max;
state.dc_level_weights_a2 = obj.dc_level_weights_a2;
state.a2_level_weight_initial_stats = obj.a2_level_weight_initial_stats;
end
function restoreObjectiveState(obj,state)
obj.e = state.e;
obj.e_dc = state.e_dc;
obj.P = state.P;
obj.save_debug = state.save_debug;
obj.debug_struct = state.debug_struct;
obj.dc_tracking_mu = state.dc_tracking_mu;
obj.dc_tracking_alpha = state.dc_tracking_alpha;
obj.dc_tracking_gamma = state.dc_tracking_gamma;
obj.dc_tracking_persistence_gain = state.dc_tracking_persistence_gain;
obj.dc_tracking_power_exponent = state.dc_tracking_power_exponent;
obj.dc_tracking_mu_eff_max = state.dc_tracking_mu_eff_max;
obj.dc_level_weights_a2 = state.dc_level_weights_a2;
obj.a2_level_weight_initial_stats = state.a2_level_weight_initial_stats;
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 [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 plotMuOptimization(obj)
if isempty(obj.mu_optimization)
return
end
X = obj.mu_optimization.XTrace;
objective = obj.mu_optimization.ObjectiveTrace;
objective = objective(:);
valid = isfinite(objective);
if isempty(X) || ~any(valid)
return
end
var_names = X.Properties.VariableNames;
n_vars = numel(var_names);
eval_idx = (1:numel(objective)).';
objective_plot = obj.positiveObjectiveForLogPlot(objective);
best_plot = obj.positiveObjectiveForLogPlot(cummin(objective));
figure(obj.mu_optimization_fignum);
clf;
t = tiledlayout(2,2,"TileSpacing","compact","Padding","compact");
title(t,"FFE Bayesian mu optimization");
nexttile;
h_candidate = semilogy(eval_idx,objective_plot,"o-","DisplayName","candidate");
obj.addOptimizationDataTips(h_candidate,X,objective,objective_plot,eval_idx,var_names);
hold on;
h_best_trace = semilogy(eval_idx,best_plot,"k-","LineWidth",1.2,"DisplayName","best so far");
h_best_trace.DataTipTemplate.DataTipRows(end+1) = dataTipTextRow("Evaluation",eval_idx);
h_best_trace.DataTipTemplate.DataTipRows(end+1) = dataTipTextRow("Best BER",best_plot);
grid on;
xlabel("Evaluation");
ylabel("BER");
legend("Location","best");
if n_vars < 2
return
end
pairs = nchoosek(1:n_vars,2);
n_pair_plots = min(size(pairs,1),3);
[~,best_idx] = min(objective);
for pair_idx = 1:n_pair_plots
nexttile;
x_name = var_names{pairs(pair_idx,1)};
y_name = var_names{pairs(pair_idx,2)};
x_data = X.(x_name);
y_data = X.(y_name);
c_data = log10(objective_plot);
h_scatter = scatter(log10(x_data),log10(y_data),35,c_data,"filled");
obj.addOptimizationDataTips(h_scatter,X,objective,objective_plot,eval_idx,var_names);
hold on;
h_best = plot(log10(x_data(best_idx)),log10(y_data(best_idx)),"kp", ...
"MarkerSize",12, ...
"MarkerFaceColor","y", ...
"DisplayName","best");
obj.addOptimizationDataTips(h_best,X(best_idx,:),objective(best_idx),objective_plot(best_idx),eval_idx(best_idx),var_names);
grid on;
xlabel("log10(" + string(x_name) + ")");
ylabel("log10(" + string(y_name) + ")");
cb = colorbar;
cb.Label.String = "log10(BER)";
title(string(x_name) + " vs " + string(y_name));
end
end
function addOptimizationDataTips(~,plot_handle,X,objective,objective_plot,eval_idx,var_names)
plot_handle.DataTipTemplate.DataTipRows(end+1) = dataTipTextRow("Evaluation",eval_idx);
plot_handle.DataTipTemplate.DataTipRows(end+1) = dataTipTextRow("BER",objective);
plot_handle.DataTipTemplate.DataTipRows(end+1) = dataTipTextRow("BER shown",objective_plot);
for var_idx = 1:numel(var_names)
var_name = var_names{var_idx};
plot_handle.DataTipTemplate.DataTipRows(end+1) = dataTipTextRow(var_name,X.(var_name));
end
end
function objective_plot = positiveObjectiveForLogPlot(~,objective)
objective_plot = objective;
positive_values = objective(isfinite(objective) & objective > 0);
if isempty(positive_values)
floor_value = 1e-12;
else
floor_value = min(positive_values) / 10;
end
objective_plot(~isfinite(objective_plot) | objective_plot <= 0) = floor_value;
end
end
end