diff --git a/Classes/04_DSP/Equalizer/ML_MLSE_DUOBINARY.m b/Classes/04_DSP/Equalizer/ML_MLSE_DUOBINARY.m new file mode 100644 index 0000000..2a9c206 --- /dev/null +++ b/Classes/04_DSP/Equalizer/ML_MLSE_DUOBINARY.m @@ -0,0 +1,409 @@ +classdef ML_MLSE_DUOBINARY < handle + % --------------------------------------------------------------------- + % W. Lanneer and Y. Lefevre, + % “Machine Learning-Based Pre-Equalizers for Maximum Likelihood + % Sequence Estimation in High-Speed PONs,” EUSIPCO 2023 + % --------------------------------------------------------------------- + % This implementation reproduces the closed-loop ML-based + % pre-equalizer training for MLSE, supporting both training and + % detection (decision-directed) modes. + % --------------------------------------------------------------------- + + properties + sps + order + e + e_tr + error + + len_tr + mu_tr + epochs_tr + + dd_mode + mu_dd + epochs_dd + adaptive_mu + + constellation + L + alpha + DIR + DIR_flip + trellis_states + traceback_depth + delta + + % Internal variables + S + Nf + nStates + nFeasible + combs + first_sym + last_sym + valid + valid_to_idx + valid_from_idx + w + + % Fast lookup + nSym + key_table + trans_index + true_to_state_idx + + % Debug metrics + ber = [] + ber_dd = [] + ce = ones(1,1) + end + + methods + function obj = ML_MLSE_DUOBINARY(options) + arguments(Input) + options.sps = 2; + options.order = 15; + options.len_tr = 4096; + options.mu_tr = 0.001; + options.epochs_tr = 5; + options.dd_mode = 1; + options.mu_dd = 1e-5; + options.epochs_dd = 5; + options.adaptive_mu = 1; + options.delta = 0; + options.traceback_depth = 1024; + options.L = 1; + end + + fn = fieldnames(options); + for n = 1:numel(fn) + obj.(fn{n}) = options.(fn{n}); + end + + obj.e = zeros(obj.order,1); + obj.error = 0; + end + + % ============================================================== + % PROCESS + % ============================================================== + function [X,X_viterbi] = process(obj, X, D) + % Normalize input RMS + X = X.normalize("mode","rms"); + obj.constellation = sort(unique(D.signal),'ascend'); + obj.nSym = numel(obj.constellation); + + if length(X)/length(D) ~= obj.sps + warning('Signal length does not fit to reference!'); + end + + % --- Parameters + obj.S = obj.nSym; + obj.Nf = obj.order * obj.sps; + obj.nStates = obj.S^obj.L; + obj.nFeasible = obj.nStates * obj.S; + + % --- Trellis mapping + obj.trellis_states = reshape(obj.constellation,1,[]); + pre_comb_mat = repmat(obj.trellis_states, obj.L, 1); + pre_comb_cell = mat2cell(pre_comb_mat, ones(1,obj.L), size(pre_comb_mat,2)); + obj.combs = fliplr(combvec(pre_comb_cell{:}).'); + obj.first_sym = obj.combs(:,1); + obj.last_sym = obj.combs(:,end); + obj.nStates = size(obj.combs,1); + + % --- Valid transitions + obj.valid = false(obj.nStates); + for from = 1:obj.nStates + for to = 1:obj.nStates + if all(obj.combs(to,2:end) == obj.combs(from,1:end-1)) + obj.valid(to,from) = true; + end + end + end + [obj.valid_to_idx,obj.valid_from_idx] = find(obj.valid); + + % --- Initialize weights + if isempty(obj.w) || any(size(obj.w) ~= [obj.Nf+1,obj.nFeasible]) + % obj.w = randn(obj.Nf+1,obj.nFeasible); + obj.w = zeros(obj.Nf+1,obj.nFeasible); + end + + % --- Fast lookup tables + [~, sym_idx_mat] = ismember(obj.combs, obj.constellation); + key_vals = 1 + sum((sym_idx_mat - 1) .* (obj.nSym .^ (0:obj.L-1)), 2); + max_key = obj.nSym^obj.L; + obj.key_table = zeros(max_key,1,'uint32'); + obj.key_table(key_vals) = 1:obj.nStates; + + obj.trans_index = sparse(obj.nStates,obj.nStates); + for i = 1:length(obj.valid_from_idx) + f = obj.valid_from_idx(i); + t = obj.valid_to_idx(i); + obj.trans_index(t,f) = i; + end + + % ============================================================== + % TRAINING + % ============================================================== + fprintf('\n--- Training mode ---\n'); + obj.equalize(X.signal, D.signal, obj.mu_tr, obj.epochs_tr, obj.len_tr, true); + obj.e_tr = obj.e; + + % ============================================================== + % DECISION-DIRECTED / TESTING + % ============================================================== + fprintf('--- Decision-directed / detection mode ---\n'); + [y, y_vit] = obj.equalize(X.signal, D.signal, obj.mu_dd, obj.epochs_dd, X.length, false); + + X_viterbi = X; + X.signal = y; + X_viterbi.signal = y_vit; + end + + % ============================================================== + % EQUALIZE + % ============================================================== + function [y,y_ref] = equalize(obj,x,d,mu,epochs,N,training) + debug = 0; + showPlots = 0; + y = zeros(N,1); + nSymbols = ceil(N/obj.sps); + + for epoch = 1:epochs + pm = zeros(obj.nStates,1); + pred = zeros(nSymbols,obj.nStates,'uint32'); + pm_sto = nan(obj.nStates,nSymbols,'like',pm); + CE_accum = 0; + + start_sample = 1; + end_sample = N; + start_symbol = 1 + floor((start_sample - 1)/obj.sps); + + % --- initialize true state + if numel(d) >= obj.L && start_symbol >= obj.L + init_seq = d(start_symbol-obj.L+1:start_symbol); + key_init = obj.seq2key(init_seq); + true_to_state_idx = obj.key_table(key_init); + if true_to_state_idx==0, true_to_state_idx=1; end + else + true_to_state_idx = uint32(1); + end + + for sample = start_sample:obj.sps:end_sample + symbol = (sample - start_sample)/obj.sps + 1; + sym_idx = start_symbol + (symbol - 1); + + % --- Observation window (with delta) + i1 = sample - obj.Nf + 1 + obj.delta; + i2 = sample + obj.delta; + buf = x(max(1,i1):min(length(x),i2)); + padL = max(0,1 - i1); + padR = max(0,i2 - length(x)); + yk = [zeros(padL,1); buf(:); zeros(padR,1)]; + yk = [yk;1]; + + % --- Branch metrics + c_hat = (yk.' * obj.w).'; + pm = pm - min(pm); + v_tilde = pm(obj.valid_from_idx) + c_hat; + + % --- allocate once + if epoch==1 && symbol==1 + obj.true_to_state_idx = ones(ceil(N/obj.sps),1,'uint32'); + end + + % --- previous "to" becomes "from" + if symbol>1 + true_from_state_idx = obj.true_to_state_idx(symbol-1); + else + true_from_state_idx = 1; + end + + % --- compute or reuse "to" state + if epoch==1 + if sym_idx>=obj.L + key_to = obj.seq2key(d(sym_idx-obj.L+1:sym_idx)); + state_idx = obj.key_table(key_to); + if state_idx==0 + state_idx = true_from_state_idx; + end + obj.true_to_state_idx(symbol) = state_idx; + else + obj.true_to_state_idx(symbol) = true_from_state_idx; + end + end + true_to_state_idx = obj.true_to_state_idx(symbol); + + % --- fast Dirac creation + dirac = zeros(obj.nFeasible,1); + trans_idx = obj.trans_index(true_to_state_idx,true_from_state_idx); + if trans_idx~=0 + dirac(trans_idx)=1; + end + + % --- ensure valid (from,to) + if ~any(dirac) + mask = obj.valid_from_idx==true_from_state_idx & ... + obj.valid_to_idx ==true_to_state_idx; + if any(mask) + dirac(mask) = 1; + else + idx = find(obj.valid_from_idx==true_from_state_idx,1,'first'); + dirac(idx) = 1; + obj.true_to_state_idx(symbol) = obj.valid_to_idx(idx); + end + end + + % =================================================================== + % TRAINING MODE (weight update) + % =================================================================== + if training + % --- Softmax and CE + v_shift = -(v_tilde - min(v_tilde)); + v_shift = min(v_shift,100); + expv = exp(v_shift); + p = expv./(sum(expv)+eps); + CE_symbol(symbol) = -log(p(dirac==1)+eps); + + % --- CE smoothing and adaptive μ + if sym_idx>obj.L + CE_smooth(symbol)=0.01*CE_symbol(symbol)+0.99*CE_symbol(symbol-1); + else + CE_smooth(symbol)=CE_symbol(symbol); + end + CE_accum=CE_accum+CE_symbol(symbol); + + % --- Gradient update + dmp=(dirac-p)'; + dL_Dw=(yk).*dmp; + if sym_idx>=obj.L + if obj.adaptive_mu + mu_eff=CE_smooth(symbol); + mu_eff=max(min(mu_eff,0.2),1e-4); + else + mu_eff=mu; + end + obj.w=obj.w - mu_eff.*dL_Dw; + end + end + + % =================================================================== + % DECODING MODE (Viterbi only) + % =================================================================== + % Compare-Select (always executed) + vmat=inf(obj.nStates,obj.nStates); + vmat(obj.valid)=v_tilde; + [pm_next,pred(symbol,:)]=min(vmat,[],2); + pm_next=pm_next-min(pm_next); + pm=pm_next; + pm_sto(:,symbol)=pm; + end + + % --- Traceback + [~,s_end]=min(pm); + vpath=zeros(symbol,1,'uint32'); + vpath(symbol)=s_end; + for n=symbol:-1:2 + vpath(n-1)=pred(n,vpath(n)); + end + + y_ref=d(start_symbol:end); + y=obj.first_sym(vpath); + + % --- BER/CE reporting and plots + if training + err=sum(y~=y_ref(1:length(y))); + ser=err/length(y); + [ber, ~] = obj.calculateDuobinaryBer(y, y_ref); + if isfinite(ber) + fprintf('Epoch %d - BER: %.2e\n',epoch,ber); + obj.ber(epoch)=ber; + else + fprintf('Epoch %d - SER: %.2e\n',epoch,ser); + obj.ber(epoch)=ser; + end + obj.ce(epoch)=CE_accum/symbol; + + if debug && mod(epoch,10)==1 && showPlots + figure(10);clf + subplot(3,2,1:2); + imagesc(obj.w);axis xy;colorbar;title('Filter W'); + subplot(3,2,3); + vtilde_mat=NaN(obj.nStates,obj.nStates); + vtilde_mat(obj.valid)=v_tilde; + imagesc(vtilde_mat);axis xy;colorbar;title('Path Metrics (v\_tilde)'); + subplot(3,2,4); + plot(1:symbol,pm_sto);title('Path Metric Evolution'); + subplot(3,2,5);hold on; + scatter(1:symbol,CE_symbol,1,'.'); + scatter(1:symbol,CE_smooth,1,'.'); + title('Cross Entropy'); + subplot(3,2,6);hold on; + yyaxis left + scatter(1:length(obj.ce),obj.ce,10,'s','filled'); + ylabel('Cross Entropy'); + yyaxis right + scatter(1:length(obj.ber),obj.ber,10,'d','filled'); + set(gca,'YScale','log'); + ylabel('BER (log)'); + xlabel('Epoch');grid on; + title('Convergence'); + drawnow; + end + else + [ber, ser] = obj.calculateDuobinaryBer(y, y_ref); + if isfinite(ber) + fprintf('DD epoch %d - BER: %.2e\n',epoch,ber); + obj.ber_dd(epoch)=ber; + else + fprintf('DD epoch %d - SER: %.2e\n',epoch,ser); + obj.ber_dd(epoch)=ser; + end + end + end + end + + end + + methods (Access=private) + % ============================================================== + % Helper: convert a detected precoded sequence back to PAM data + % ============================================================== + function data = invertDuobinaryPrecoder(obj, precoded) + encoded = Duobinary().encode(precoded,"M",obj.S); + encoded_signal = Signal(encoded); + decoded_signal = Duobinary().decode(encoded_signal,"M",obj.S); + data = decoded_signal.signal; + end + + % ============================================================== + % Helper: calculate BER after inverting the duobinary precoder + % ============================================================== + function [ber, ser] = calculateDuobinaryBer(obj, detected, reference) + n = min(numel(detected), numel(reference)); + detected = detected(1:n); + reference = reference(1:n); + + ser = sum(detected ~= reference) / n; + + detected_data = obj.invertDuobinaryPrecoder(detected); + reference_data = obj.invertDuobinaryPrecoder(reference); + mapper = PAMmapper(obj.S, 0); + detected_bits = mapper.demap(detected_data); + reference_bits = mapper.demap(reference_data); + [~,~,ber,~] = calc_ber(reference_bits, detected_bits, ... + "skip_front",10,"skip_end",10,"returnErrorLocation",1); + end + + % ============================================================== + % Helper: Sequence → key (always scalar) + % ============================================================== + function key = seq2key(obj, seq) + [~, idx] = ismember(flip(seq), obj.constellation); + pow = (obj.nSym .^ (0:obj.L-1)).'; + key = 1 + sum((idx(:) - 1) .* pow); + end + end +end diff --git a/Classes/04_DSP/Equalizer/VNLE.m b/Classes/04_DSP/Equalizer/VNLE.m index 51d96a5..1281dd5 100644 --- a/Classes/04_DSP/Equalizer/VNLE.m +++ b/Classes/04_DSP/Equalizer/VNLE.m @@ -30,6 +30,9 @@ classdef VNLE < handle optmize_mus = 0; mu_optimization mu_optimization_iter = 0; + mu_optimization_len + plot_mu_optimization = 0 + mu_optimization_fignum = 3020; x_norm ce @@ -55,6 +58,9 @@ classdef VNLE < handle 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 = 3020; end @@ -110,7 +116,6 @@ classdef VNLE < handle lbdesc = [num2str(obj.order),' tap FFE']; X = X.logbookentry(lbdesc); % append to logbook - N = X; N = X - D; @@ -138,13 +143,18 @@ classdef VNLE < handle ones(1,obj.ce(3))*mu(3) ]); end + x = x(:); + d = d(:); x = [zeros(floor(obj.order(1)/2),1); x; zeros(obj.order(1),1)]; + n_symbols = floor(N / obj.sps); + y = zeros(n_symbols,1); + d_hat = zeros(n_symbols,1); if showviz - f = figure(111); + figure(111); subplot(2,2,1:2); hold on - a = scatter(1:numel(x),x,1,'.'); + scatter(1:numel(x),x,1,'.'); a2 = scatter(1,1,1,'.'); a3 = scatter(1,1,2,'.'); a4 = xline(1); @@ -215,40 +225,50 @@ classdef VNLE < handle mu_range = [1e-5, 1e-2]; mu_dc_range = [1e-5, 1e-1]; - vars = [optimizableVariable("mu_tr",mu_range,"Transform","log"), ... - optimizableVariable("mu_dd",mu_range,"Transform","log")]; + [x_opt,d_opt,N_opt] = obj.optimizationSignals(x,d); + + vars = obj.muOptimizableVariables("mu_tr",mu_range); + vars = [vars, obj.muOptimizableVariables("mu_dd",mu_range)]; optimize_mu_dc = obj.mu_dc ~= 0; if optimize_mu_dc vars = [vars, optimizableVariable("mu_dc",mu_dc_range,"Transform","log")]; end obj.mu_optimization_iter = 0; - obj.mu_optimization = bayesopt(@(p)obj.muObjective(p,x,d),vars, ... - "MaxObjectiveEvaluations",10, ... + fprintf("VNLE 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; - obj.mu_dd = obj.mu_optimization.XAtMinObjective.mu_dd; + + obj.mu_tr = obj.muVectorFromParams(obj.mu_optimization.XAtMinObjective,"mu_tr"); + obj.mu_dd = obj.muVectorFromParams(obj.mu_optimization.XAtMinObjective,"mu_dd"); if optimize_mu_dc obj.mu_dc = obj.mu_optimization.XAtMinObjective.mu_dc; end - objective_db = 10*log10(obj.mu_optimization.MinObjective); if optimize_mu_dc - fprintf("\nVNLE mu opt done: mu_tr=%9.3e, mu_dd=%9.3e, mu_dc=%9.3e, MSE=%9.3e, MSE_dB=%7.2f dB\n", ... - obj.mu_tr,obj.mu_dd,obj.mu_dc,obj.mu_optimization.MinObjective,objective_db); + fprintf("\nVNLE mu opt done: mu_tr=[%s], mu_dd=[%s], mu_dc=%9.3e, BER=%9.3e\n", ... + obj.formatMuVector(obj.mu_tr),obj.formatMuVector(obj.mu_dd), ... + obj.mu_dc,obj.mu_optimization.MinObjective); else - fprintf("\nVNLE mu opt done: mu_tr=%9.3e, mu_dd=%9.3e, MSE=%9.3e, MSE_dB=%7.2f dB\n", ... - obj.mu_tr,obj.mu_dd,obj.mu_optimization.MinObjective,objective_db); + fprintf("\nVNLE mu opt done: mu_tr=[%s], mu_dd=[%s], BER=%9.3e\n", ... + obj.formatMuVector(obj.mu_tr),obj.formatMuVector(obj.mu_dd), ... + obj.mu_optimization.MinObjective); + end + + if obj.plot_mu_optimization + obj.plotMuOptimization(); end end function objective = muObjective(obj,params,x,d) - old_debug = obj.save_debug; - old_mu_dc = obj.mu_dc; - obj.save_debug = 1; + old_state = obj.captureObjectiveState(); + cleanup = onCleanup(@()obj.restoreObjectiveState(old_state)); + + obj.save_debug = 0; optimize_mu_dc = ismember("mu_dc",string(params.Properties.VariableNames)); if optimize_mu_dc obj.mu_dc = params.mu_dc; @@ -257,34 +277,217 @@ classdef VNLE < handle obj.e = zeros(sum(obj.ce),1); obj.e_dc = 0; obj.debug_struct = struct(); - obj.equalize(x,d,params.mu_tr,obj.epochs_tr,obj.len_tr,1,0); - obj.equalize(x,d,params.mu_dd,obj.epochs_dd,numel(x),0,0); + muTrCandidate = obj.muVectorFromParams(params,"mu_tr"); + muDdCandidate = obj.muVectorFromParams(params,"mu_dd"); + N_tr = min(obj.len_tr,numel(x)); + obj.equalize(x,d,muTrCandidate,obj.epochs_tr,N_tr,1,0); + [signal,~] = obj.equalize(x,d,muDdCandidate,obj.epochs_dd,numel(x),0,0); - objective = mean(obj.debug_struct.error(end,:),"omitnan"); + [ber,errors] = obj.berObjective(signal,d); + objective = ber; if ~isfinite(objective) objective = inf; end - objective_db = 10*log10(objective); obj.mu_optimization_iter = obj.mu_optimization_iter + 1; if optimize_mu_dc - fprintf("\rVNLE mu opt %02d: mu_tr=%9.3e, mu_dd=%9.3e, mu_dc=%9.3e, MSE=%9.3e, MSE_dB=%7.2f dB", ... - obj.mu_optimization_iter,params.mu_tr,params.mu_dd,params.mu_dc,objective,objective_db); + fprintf("\rVNLE mu opt %02d: mu_tr=[%s], mu_dd=[%s], mu_dc=%9.3e, BER=%9.3e, errors=%d", ... + obj.mu_optimization_iter,obj.formatMuVector(muTrCandidate), ... + obj.formatMuVector(muDdCandidate),params.mu_dc,ber,errors); else - fprintf("\rVNLE mu opt %02d: mu_tr=%9.3e, mu_dd=%9.3e, MSE=%9.3e, MSE_dB=%7.2f dB", ... - obj.mu_optimization_iter,params.mu_tr,params.mu_dd,objective,objective_db); + fprintf("\rVNLE mu opt %02d: mu_tr=[%s], mu_dd=[%s], BER=%9.3e, errors=%d", ... + obj.mu_optimization_iter,obj.formatMuVector(muTrCandidate), ... + obj.formatMuVector(muDdCandidate),ber,errors); end - obj.save_debug = old_debug; - obj.mu_dc = old_mu_dc; + + 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); + skip_front = min(1000,max(0,floor(numel(rx_bits.signal) / 4))); + [~,errors,ber,~] = calc_ber(rx_bits.signal,tx_bits.signal, ... + "skip_front",skip_front, ... + "skip_end",0, ... + "returnErrorLocation",1); + end + + function vars = muOptimizableVariables(obj,prefix,mu_range) + vars = optimizableVariable.empty; + activeOrders = find(obj.order > 0); + + for idx = 1:numel(activeOrders) + orderIdx = activeOrders(idx); + varName = sprintf("%s_%d",prefix,orderIdx); + vars = [vars, optimizableVariable(varName,mu_range,"Transform","log")]; %#ok + end + end + + function mu = muVectorFromParams(obj,params,prefix) + mu = zeros(1,3); + for orderIdx = 1:numel(mu) + varName = sprintf("%s_%d",prefix,orderIdx); + if ismember(varName,string(params.Properties.VariableNames)) + mu(orderIdx) = params.(varName); + elseif numel(obj.(char(prefix))) >= orderIdx + mu(orderIdx) = obj.(char(prefix))(orderIdx); + else + mu(orderIdx) = obj.(char(prefix))(1); + end + end + 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,"VNLE 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_dc = obj.e_dc; + state.error = obj.error; + state.mu_dc = obj.mu_dc; + state.save_debug = obj.save_debug; + state.debug_struct = obj.debug_struct; + end + + function restoreObjectiveState(obj,state) + obj.e = state.e; + obj.e_dc = state.e_dc; + obj.error = state.error; + obj.mu_dc = state.mu_dc; + obj.save_debug = state.save_debug; + obj.debug_struct = state.debug_struct; + end + + function s = formatMuVector(~,mu) + s = strtrim(sprintf("%9.3e ",mu)); end %% Functions needed During Adaption function x_in_vnle_format = calcVNLENonlinVecs(~,x_in_block,I_2,I_3,N_,norm_) % These are the second and third order input signal products of the VNLE EQ % ∑ h1 x_in(k-n1) + ∑∑ h2 x_in(k-n1)*x_in(k-n2) + ∑∑∑ h3 x_in(k-n1)*x_in(k-n2)*x_in(k-n3) + x_in_block = x_in_block(:); l1=length(x_in_block); - l2=length(I_2); - l3=length(I_3); + l2=size(I_2,1); + l3=size(I_3,1); final_length = l1+l2+l3; x_in_vnle_format = zeros(final_length,1); @@ -382,6 +585,7 @@ classdef VNLE < handle end function powerNorm = calcPowerNormalization(~,v) + v = v(:); powerNorm(1) = sqrt(mean(abs(v ).^2)); powerNorm(2) = sqrt(mean(abs(v.^2).^2)); diff --git a/Classes/DataBaseHandler/Metricstruct.m b/Classes/DataBaseHandler/Metricstruct.m index fe9a818..b755718 100644 --- a/Classes/DataBaseHandler/Metricstruct.m +++ b/Classes/DataBaseHandler/Metricstruct.m @@ -13,6 +13,9 @@ classdef Metricstruct BER_precoded (1,1) double {mustBeNumeric, mustBeNonnegative} = 0 numBitErr_precoded (1,1) double {mustBeInteger, mustBeNonnegative} = 0 + % BER_DB_memoryless (1,1) double {mustBeNumeric} = NaN + % BER_DB_sequencedetection (1,1) double {mustBeNumeric} = NaN + SNR (1,1) double {mustBeNumeric} = NaN SNR_level (:,1) double {mustBeNumeric} = [] STD (1,1) double {mustBeNumeric} = NaN diff --git a/Functions/EQ_blocks/duobinary_signaling.m b/Functions/EQ_blocks/duobinary_signaling.m index 7f13c35..5c342c7 100644 --- a/Functions/EQ_blocks/duobinary_signaling.m +++ b/Functions/EQ_blocks/duobinary_signaling.m @@ -21,10 +21,13 @@ arguments tx_symbols tx_bits options.precode_mode db_mode + options.decoding_mode db_decoder = db_decoder.sequencedetection; + options.showAnalysis = 0 options.eth_style_symbol_mapping = 0 options.postFFE = [] options.database = [] + end %% Process signals through equalizer @@ -35,43 +38,72 @@ if ~isempty(options.postFFE) [eq_signal, eq_noise] = options.postFFE.process(eq_signal, tx_symbols); end -if isa(mlse_,'MLSE_viterbi') - [mlse_signal] = mlse_.process(eq_signal); -else - % Aufpassen mit welcher Sequenz man hier vergleicht für LLR stuff... - % gespeichtere "Symbols" sind schon DB codiert, das wollen wir hier - % nicht! Sondern die precoded aber nicht db-encoded müssen als ref in - % die LLR berechnung gehen! - ref_sym = PAMmapper(M,0).map(tx_bits); %ist klar - ref_sym_dpc = Duobinary().precode(ref_sym); % precoded - % ref_sym_dbenc = Duobinary().encode(ref_sym_dpc); %encoded - das wurde gesendet! - % ref_sym_dec = Duobinary().decode(ref_sym_dbenc); %ref_sym wieder zurück! - mlse_.trellis_states = PAMmapper(M,0).levels; - mlse_.trellis_state_mode = 1; - [mlse_signal,LLR,GMI_MLSE] = mlse_.process(eq_signal,ref_sym_dpc); +run_both_detection_schemes = 1; +if options.decoding_mode == db_decoder.sequencedetection || run_both_detection_schemes %MLSE + + + if isa(mlse_,'MLSE_viterbi') + [mlse_signal] = mlse_.process(eq_signal); + else + + % Aufpassen mit welcher Sequenz man hier vergleicht für LLR stuff... + % gespeichtere "Symbols" sind schon DB codiert, das wollen wir hier + % nicht! Sondern die precoded aber nicht db-encoded müssen als ref in + % die LLR berechnung gehen! + ref_sym = PAMmapper(M,0).map(tx_bits); %ist klar + ref_sym_dpc = Duobinary().precode(ref_sym); % precoded + % ref_sym_dbenc = Duobinary().encode(ref_sym_dpc); %encoded - das wurde gesendet! + % ref_sym_dec = Duobinary().decode(ref_sym_dbenc); %ref_sym wieder zurück! + + mlse_.trellis_states = PAMmapper(M,0).levels; + % mlse_.scale_mode = 2; + if M == 6 + mlse_.trellis_exclusion = 0; + end + mlse_.debug = 0; + [mlse_signal,LLR,GMI_MLSE] = mlse_.process(eq_signal.normalize("mode","rms"),ref_sym_dpc); + end + + % Apply duobinary encoding and decoding + mlse_signal = Duobinary().encode(mlse_signal,"M",M); + + mlse_signal = Duobinary().decode(mlse_signal,"M",M); + + % Demap symbols to bits + rx_bits = PAMmapper(M, 0, "eth_style", options.eth_style_symbol_mapping).demap(mlse_signal); + + [bits_db, ~, ber_sequencedetection, ~] = calc_ber(rx_bits.signal, tx_bits.signal, "skip_front", 100, "skip_end", 150, "returnErrorLocation", 1); + +end + +if options.decoding_mode == db_decoder.memoryless || run_both_detection_schemes + + + db_ref_constellation = unique(tx_symbols.signal); + eq_signal_hd = PAMmapper(M, 0).quantize(eq_signal,'custom_const',db_ref_constellation.'); + eq_signal_hd = Duobinary().decode(eq_signal_hd,"M",M); + + % Demap + rx_bits = PAMmapper(M, 0, "eth_style", options.eth_style_symbol_mapping).demap(eq_signal_hd); + + [bits_db, ~, ber_memoryless, ~] = calc_ber(rx_bits.signal, tx_bits.signal, "skip_front", 100, "skip_end", 150, "returnErrorLocation", 1); + end -% tx_symbols_ = Duobinary().decode(tx_symbols); -% [mlse_signal,~,GMI_MLSE] = mlse_.process(eq_signal,tx_symbols); - -% Apply duobinary encoding and decoding -mlse_signal = Duobinary().encode(mlse_signal); -mlse_signal = Duobinary().decode(mlse_signal); - -% Demap symbols to bits -rx_bits = PAMmapper(M, 0, "eth_style", options.eth_style_symbol_mapping).demap(mlse_signal); - -%% Calculate BER and metrics -[bits_db, errors_db, ber_db, error_pos] = calc_ber(rx_bits.signal, tx_bits.signal, "skip_front", 100, "skip_end", 150, "returnErrorLocation", 1); +% Calculate BER and metrics % Calculate performance metrics after duobinary FFE! [snr, snr_lvl] = calc_snr(tx_symbols.signal, eq_noise.signal); %SNR of duobinary sequence - not directly comparable to -[gmi] = calc_air(eq_signal, tx_symbols, "skip_front", 10000, "skip_end", 10000); -air = tx_symbols.fs .* floor(log2(double(M))*10)/10 .* gmi ./ log2(double(M)); +% [gmi] = calc_air(eq_signal, tx_symbols, "skip_front", 10000, "skip_end", 10000); % Not working for Duobinary == Channel with memory +% air = tx_symbols.fs .* floor(log2(double(M))*10)/10 .* gmi ./ log2(double(M)); +if options.decoding_mode == db_decoder.sequencedetection || run_both_detection_schemes + [gmi] = GMI_MLSE; + air = tx_symbols.fs .* floor(log2(double(M))*10)/10 .* gmi ./ log2(double(M)); +end [evm_total, evm_lvl] = calc_evm(eq_signal, tx_symbols); [std_total, std_lvl] = calc_std(eq_signal, tx_symbols); [std_rxraw_total, std_rxraw_lvl] = calc_std(rx_signal.resample("fs_out", tx_symbols.fs), tx_symbols); @@ -91,9 +123,12 @@ db_results.metrics.result_id = NaN; db_results.metrics.run_id = NaN; db_results.metrics.eqParam_id = NaN; db_results.metrics.date_of_processing = datetime('now'); -db_results.metrics.BER = ber_db; +db_results.metrics.BER = ber_sequencedetection; % THIS IS A CONVENTION +db_results.metrics.BER_precoded = ber_memoryless; % THIS IS A CONVENTION +% db_results.metrics.BER_DB_memoryless = ber_memoryless; +% db_results.metrics.BER_DB_sequencedetection = ber_sequencedetection; db_results.metrics.numBits = bits_db; -db_results.metrics.numBitErr = errors_db; +% db_results.metrics.numBitErr = NaN; db_results.metrics.SNR = snr; db_results.metrics.SNR_level = snr_lvl; db_results.metrics.STD = std_total; diff --git a/Functions/EQ_recipes/dsp_400g_recipe.m b/Functions/EQ_recipes/dsp_400g_recipe.m new file mode 100644 index 0000000..751534c --- /dev/null +++ b/Functions/EQ_recipes/dsp_400g_recipe.m @@ -0,0 +1,449 @@ +function output = dsp_400g_recipe(Scpe_sig_raw, Symbols, Tx_bits, options) +%dsp_400g_recipe Run 400G equalizer schemes from a synchronized scope signal. + +arguments + Scpe_sig_raw + Symbols + Tx_bits + options.fsym + options.M + options.duob_mode + options.dataTable table + options.userParameters struct = struct() + options.preprocess_mode string = "auto" + options.tx_pulseformer = [] + options.debug_plots (1,1) logical = false +end + +p = defaultRecipeParameters(); +p.preprocess_mode = options.preprocess_mode; +p = applyUserParameters(p, options.userParameters); + +output = struct(); + +Scpe_sig = preprocessSignal(Scpe_sig_raw, Symbols, options.fsym, ... + "mode", p.preprocess_mode, ... + "tx_pulseformer", options.tx_pulseformer, ... + "debug_plots", options.debug_plots); +[Scpe_sig, Symbols, Tx_bits] = alignDspInputs(Scpe_sig, Symbols, Tx_bits, p.eq_sps); +if isempty(p.ml_mlse_len_tr) + p.ml_mlse_len_tr = floor(length(Scpe_sig) / 4); +end + +if p.plot_input_signal + showLevelScatter(Scpe_sig.normalize("mode","rms"), Symbols, ... + "fsym", options.fsym, ... + "fignum", p.input_plot_fignum, ... + "normalize", true); +end + +if options.duob_mode ~= db_mode.db_encoded + if p.run_ffe + eq_ffe = FFE("epochs_tr", p.epochs_tr, ... + "epochs_dd", p.epochs_dd, ... + "len_tr", p.len_tr, ... + "mu_dd", p.ffe_mu_dd, ... + "mu_tr", p.ffe_mu_tr, ... + "order", p.ffe_order(1), ... + "sps", p.eq_sps, ... + "decide", false, ... + "optmize_mus", p.optimize_mus, ... + "dd_mode", p.dd_mode, ... + "adaption_technique", p.ffe_adaption, ... + "dc_tracking_mu", p.mu_dc); + + [ffe_results, equalized_signal] = runFfe(eq_ffe, "FFE", ... + Scpe_sig, Symbols, Tx_bits, options); + ffe_results.config.equalizer_structure = equalizer_structure.ffe; + ffe_results.recipe_config = collectRecipeConfig("ffe", eq_ffe, p, options); + output.ffe_package = ffe_results; + + if p.plot_output_signals + plotEqSignals(equalized_signal, Symbols, options, p.output_plot_fignum, -1); + end + end + + + if p.run_vnle + % eq_vnle = VNLE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2, ... + % "mu_dd",[0.0004 0.0005 0.0006],"mu_tr",[0.0001 0.0008 0.001], ... + % "order",[150,5,5],"sps",2,"decide",1, ... + % "optmize_mus",1,"mu_optimization_len",2^15); + + eq_vnle = EQ("Ne", p.vnle_ffe_order, ... + "Nb", p.vnle_dfe_order, ... + "training_length", p.len_tr, ... + "training_loops", p.epochs_tr, ... + "dd_loops", p.epochs_dd, ... + "K", p.eq_K, ... + "DCmu", p.mu_dc, ... + "DDmu", [p.eq_mu_ffe p.mu_dfe], ... + "DFEmu", p.dfe_mu_feedback, ... + "FFEmu", 0, ... + "plotfinal", 0, ... + "ideal_dfe", false); + + [vnle_results, equalized_signal] = runFfe(eq_vnle, "VNLE", ... + Scpe_sig, Symbols, Tx_bits, options); + vnle_results.config.equalizer_structure = equalizer_structure.ffe; + ffe_results.recipe_config = collectRecipeConfig("vnle", eq_ffe, p, options); + output.vnle_package = vnle_results; + + if p.plot_output_signals + plotEqSignals(equalized_signal, Symbols, options, p.output_plot_fignum, -1); + end + end + + + if p.run_dfe + eq_dfe = EQ("Ne", p.dfe_ffe_order, ... + "Nb", p.dfe_feedback_order, ... + "training_length", p.len_tr, ... + "training_loops", p.epochs_tr, ... + "dd_loops", p.epochs_dd, ... + "K", p.eq_K, ... + "DCmu", p.mu_dc, ... + "DDmu", [p.eq_mu_ffe p.mu_dfe], ... + "DFEmu", p.dfe_mu_feedback, ... + "FFEmu", 0, ... + "plotfinal", 0, ... + "ideal_dfe", false); + + dfe_results = ffe(eq_dfe, options.M, Scpe_sig, Symbols, Tx_bits, ... + "precode_mode", options.duob_mode, ... + "showAnalysis", options.debug_plots, ... + "postFFE", [], ... + "eth_style_symbol_mapping", 0); + dfe_results.config.equalizer_structure = equalizer_structure.dfe; + dfe_results.metrics.print("description", resultDescription("DFE", options)); + dfe_results.recipe_config = collectRecipeConfig("dfe", eq_dfe, p, options); + output.dfe_package = dfe_results; + end + + if p.run_vnle_mlse + eq_vnle = EQ("Ne", p.vnle_ffe_order, ... + "Nb", p.vnle_dfe_order, ... + "training_length", p.len_tr, ... + "training_loops", p.epochs_tr, ... + "dd_loops", p.epochs_dd, ... + "K", p.eq_K, ... + "DCmu", p.mu_dc, ... + "DDmu", [p.eq_mu_ffe p.mu_dfe], ... + "DFEmu", p.dfe_mu_feedback, ... + "FFEmu", 0, ... + "plotfinal", 0, ... + "ideal_dfe", false); + pf = Postfilter("ncoeff", p.pf_ncoeffs, "useBurg", true); + mlse = buildMlse(options.M, options.duob_mode, p, "pf_mlse"); + + [vnle_results, mlse_results] = vnle_postfilter_mlse(eq_vnle, pf, mlse, ... + options.M, Scpe_sig, Symbols, Tx_bits, ... + "precode_mode", options.duob_mode, ... + "showAnalysis", options.debug_plots, ... + "postFFE", [], ... + "eth_style_symbol_mapping", 0); + + vnle_results.config.equalizer_structure = equalizer_structure.vnle; + vnle_results.recipe_config = collectRecipeConfig("vnle", eq_vnle, p, options); + mlse_results.config.equalizer_structure = equalizer_structure.vnle_pf_mlse; + mlse_results.recipe_config = collectRecipeConfig("vnle_pf_mlse", eq_vnle, p, options); + vnle_results.metrics.print("description", resultDescription("VNLE", options)); + mlse_results.metrics.print("description", resultDescription("VNLE + PF + MLSE", options)); + + output.vnle_package = vnle_results; + output.mlse_package = mlse_results; + end + + if p.run_dbtgt + eq_dbtgt = EQ("Ne", p.dbtgt_ffe_order, ... + "Nb", p.dbtgt_dfe_order, ... + "training_length", p.len_tr, ... + "training_loops", p.epochs_tr, ... + "dd_loops", p.epochs_dd, ... + "K", p.eq_K, ... + "DCmu", p.mu_dc, ... + "DDmu", [p.eq_mu_ffe p.mu_dfe], ... + "DFEmu", p.dfe_mu_feedback, ... + "FFEmu", 0, ... + "plotfinal", 0, ... + "ideal_dfe", true); + mlse_db = buildMlse(options.M, options.duob_mode, p, "db_target"); + + dbtgt_results = runDuobinaryTarget(eq_dbtgt, mlse_db, ... + Scpe_sig, Symbols, Tx_bits, options, p); + dbtgt_results.config.equalizer_structure = equalizer_structure.vnle_db_mlse; + dbtgt_results.recipe_config = collectRecipeConfig("vnle_db_mlse", eq_dbtgt, p, options); + dbtgt_results.metrics.print("description", resultDescription("VNLE DB target + MLSE", options)); + output.dbtgt_package = dbtgt_results; + end + + if p.run_ml_mlse + ml_mlse_equalizer = ML_MLSE("epochs_tr", p.ml_mlse_epochs_tr, ... + "epochs_dd", p.ml_mlse_epochs_dd, ... + "len_tr", p.ml_mlse_len_tr, ... + "mu_dd", p.ml_mlse_mu_dd, ... + "mu_tr", p.ml_mlse_mu_tr, ... + "order", p.ml_mlse_order, ... + "sps", p.eq_sps, ... + "traceback_depth", p.ml_mlse_traceback_depth, ... + "L", p.ml_mlse_L, ... + "delta", p.ml_mlse_delta, ... + "adaptive_mu", p.ml_mlse_adaptive_mu); + + ml_mlse_results = ml_mlse(ml_mlse_equalizer, options.M, ... + Scpe_sig, Symbols, Tx_bits, ... + "precode_mode", options.duob_mode); + ml_mlse_results.config.equalizer_structure = equalizer_structure.ml_mlse; + ml_mlse_results.recipe_config = collectRecipeConfig("ml_mlse", ml_mlse_equalizer, p, options); + output.mlmlse_package = ml_mlse_results; + end +else + if p.run_ml_mlse_db + + ml_mlse_db_equalizer = ML_MLSE_DUOBINARY("epochs_tr", p.ml_mlse_epochs_tr, ... + "epochs_dd", p.ml_mlse_epochs_dd, ... + "len_tr", p.ml_mlse_len_tr, ... + "mu_dd", p.ml_mlse_mu_dd, ... + "mu_tr", p.ml_mlse_mu_tr, ... + "order", p.ml_mlse_order, ... + "sps", p.eq_sps, ... + "traceback_depth", p.ml_mlse_traceback_depth, ... + "L", p.ml_mlse_L, ... + "delta", p.ml_mlse_delta, ... + "adaptive_mu", p.ml_mlse_adaptive_mu); + + % Use a precoded reference for the detector; the received waveform remains encoded. + ref_sym = PAMmapper(options.M,0).map(Tx_bits); + Symbols_precoded = Duobinary().precode(ref_sym); % precoded + + + ml_mlse_db_results = ml_mlse(ml_mlse_db_equalizer, options.M, ... + Scpe_sig, Symbols_precoded, Tx_bits, ... + "precode_mode", db_mode.db_precoded); + + ml_mlse_db_results.config.equalizer_structure = equalizer_structure.ml_mlse; + ml_mlse_db_results.config.comment = 'function: ML-based MLSE; duobinary encoded'; + ml_mlse_db_results.recipe_config = collectRecipeConfig("ml_mlse_db", ml_mlse_db_equalizer, p, options); + output.mlmlse_db_package = ml_mlse_db_results; + end + + if p.run_mlse_db + eq_db_enc = EQ("Ne", p.dbtgt_ffe_order, ... + "Nb", p.dbtgt_dfe_order, ... + "training_length", p.len_tr, ... + "training_loops", p.epochs_tr, ... + "dd_loops", p.epochs_dd, ... + "K", p.eq_K, ... + "DCmu", p.mu_dc, ... + "DDmu", [p.eq_mu_ffe p.mu_dfe], ... + "DFEmu", p.dfe_mu_feedback, ... + "FFEmu", 0, ... + "plotfinal", 0, ... + "ideal_dfe", true); + mlse_db_enc = MLSE("DIR", [1,1], ... + "duobinary_output", 0, ... + "M", options.M, ... + "trellis_states", PAMmapper(options.M,0).levels); + + if isempty(p.decoding_mode) + mlse_db_results = duobinary_signaling(eq_db_enc, mlse_db_enc, ... + options.M, Scpe_sig, Symbols, Tx_bits, ... + "precode_mode", options.duob_mode, ... + "showAnalysis", options.debug_plots, ... + "postFFE", []); + mlse_db_results.metrics.print("description", resultDescription(["DB Encoded; MLSE"], options)); + else + mlse_db_results = duobinary_signaling(eq_db_enc, mlse_db_enc, ... + options.M, Scpe_sig, Symbols, Tx_bits, ... + "precode_mode", options.duob_mode, ... + "showAnalysis", options.debug_plots, ... + "postFFE", [],"decoding_mode",p.decoding_mode); + mlse_db_results.metrics.print("description", resultDescription(["DB Encoded "+string(p.decoding_mode)], options)); + end + mlse_db_results.config.equalizer_structure = equalizer_structure.db_encoded; + mlse_db_results.recipe_config = collectRecipeConfig("mlse_db", eq_db_enc, p, options); + output.mlse_db_package = mlse_db_results; + + end +end + +end + +function p = defaultRecipeParameters() +p = struct(); + +p.run_ffe = false; +p.run_vnle = false; +p.run_dfe = false; +p.run_vnle_mlse = false; +p.run_dbtgt = false; +p.run_ml_mlse = true; % non-encoded and precoded branches +p.run_ml_mlse_db = true; % db_encoded: ML-based MLSE +p.run_mlse_db = true; % db_encoded: conventional MLSE + + +p.preprocess_mode = "auto"; +p.plot_input_signal = false; +p.plot_output_signals = false; +p.input_plot_fignum = 400; +p.output_plot_fignum = 410; + +p.eq_sps = 2; +p.eq_K = 2; +p.len_tr = 4096*2; +p.epochs_tr = 5; +p.epochs_dd = 5; +p.dd_mode = true; +p.optimize_mus = true; + +p.ffe_order = [50, 0, 0]; +p.dfe_ffe_order = [50, 5, 5]; +p.vnle_ffe_order = [50, 5, 5]; +p.dbtgt_ffe_order = [50, 5, 5]; +p.vnle_dfe_order = [0, 0, 0]; +p.dbtgt_dfe_order = [0, 0, 0]; +p.dfe_feedback_order = [2, 0, 0]; + +p.eq_mu_ffe = [0.0001, 0.0008, 0.001]; +p.ffe_mu_tr = 0.4; +p.ffe_mu_dd = 0.1; +p.ffe_adaption = "nlms"; +p.mu_dfe = 0.0004; +p.mu_dc = 1.021e-05; +p.dfe_mu_feedback = 0.005; + +p.pf_ncoeffs = 1; +p.use_viterbi = false; +p.mlse_scale_mode = 2; +p.mlse_trellis_state_mode = 2; +p.dbtgt_trellis_state_mode = 3; +p.decoding_mode = []; + +p.ml_mlse_mu_tr = 0.03; +p.ml_mlse_mu_dd = 0.03; +p.ml_mlse_epochs_tr = 100; +p.ml_mlse_epochs_dd = 1; +p.ml_mlse_len_tr = []; +p.ml_mlse_order = 11; +p.ml_mlse_traceback_depth = 256; +p.ml_mlse_L = 1; +p.ml_mlse_delta = 4; +p.ml_mlse_adaptive_mu = false; +end + +function p = applyUserParameters(p, userParameters) +if isempty(userParameters) + return +end + +paramNames = fieldnames(userParameters); +for paramIdx = 1:numel(paramNames) + paramName = paramNames{paramIdx}; + if ~isfield(p, paramName) + warning("dsp_400g_recipe:UnknownUserParameter", ... + "Ignoring unknown user parameter '%s'.", paramName); + continue + end + p.(paramName) = userParameters.(paramName); +end +end + +function [Scpe_sig, Symbols, Tx_bits] = alignDspInputs(Scpe_sig, Symbols, Tx_bits, sps) +nSymbols = min(length(Symbols), floor(length(Scpe_sig) / sps)); +if nSymbols <= 0 + error("dsp_400g_recipe:EmptyAlignedSignal", ... + "No overlapping samples remain after preprocessing and synchronization."); +end + +Scpe_sig.signal = real(Scpe_sig.signal(1:sps*nSymbols)); +Symbols.signal = Symbols.signal(1:nSymbols,:); + +% if isprop(Tx_bits, "signal") +% Tx_bits.signal = Tx_bits.signal(1:nSymbols,:); <- THIS IS WRONG!! +% 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, ... + "showAnalysis", options.debug_plots, ... + "postFFE", [], ... + "eth_style_symbol_mapping", 0); + +ffe_results.metrics.print("description", resultDescription(description, options)); +end + +function dbtgt_results = runDuobinaryTarget(eq_dbtgt, mlse_db, ... + Scpe_sig, Symbols, Tx_bits, options, p) +if isempty(p.decoding_mode) + dbtgt_results = duobinary_target(eq_dbtgt, mlse_db, options.M, ... + Scpe_sig, Symbols, Tx_bits, ... + "precode_mode", options.duob_mode, ... + "showAnalysis", options.debug_plots, ... + "postFFE", []); +else + dbtgt_results = duobinary_target(eq_dbtgt, mlse_db, options.M, ... + Scpe_sig, Symbols, Tx_bits, ... + "precode_mode", options.duob_mode, ... + "showAnalysis", options.debug_plots, ... + "postFFE", [], ... + "decoding_mode", p.decoding_mode); +end +end + +function mlse = buildMlse(M, duobMode, p, mode) +if p.use_viterbi + mlse = MLSE_viterbi("duobinary_output", 0, ... + "M", M, ... + "trellis_states", PAMmapper(M,0).levels); + return +end + +if duobMode == db_mode.no_db && M == 6 + trellisExclusion = true; +else + trellisExclusion = false; +end + +switch string(mode) + case "db_target" + mlse = MLSE("DIR", [1,1], ... + "duobinary_output", 0, ... + "M", M, ... + "trellis_states", PAMmapper(M,0).levels, ... + "scale_mode", p.mlse_scale_mode, ... + "trellis_exclusion", trellisExclusion, ... + "trellis_state_mode", p.dbtgt_trellis_state_mode); + otherwise + mlse = MLSE("duobinary_output", 0, ... + "M", M, ... + "trellis_states", PAMmapper(M,0).levels, ... + "scale_mode", p.mlse_scale_mode, ... + "trellis_exclusion", trellisExclusion, ... + "trellis_state_mode", p.mlse_trellis_state_mode); +end +end + +function config = collectRecipeConfig(storageName, eqObject, p, options) +config = struct(); +config.storage_name = char(storageName); +config.recipe = "dsp_400g_recipe"; +config.eq_class = class(eqObject); +config.run_id = options.dataTable.run_id; +config.user_parameters = options.userParameters; +config.parameters = p; +end + +function description = resultDescription(prefix, options) +dt = options.dataTable; +description = sprintf('%s; run %d; PAM-%d; %.0f GBd; %.0f km', ... + prefix, dt.run_id, dt.pam_level, dt.symbolrate * 1e-9, dt.fiber_length); +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 diff --git a/Functions/EQ_recipes/mpi_recipe_dev.m b/Functions/EQ_recipes/mpi_recipe_dev.m index e056ff9..358a1b1 100644 --- a/Functions/EQ_recipes/mpi_recipe_dev.m +++ b/Functions/EQ_recipes/mpi_recipe_dev.m @@ -22,7 +22,7 @@ run_a2_tracked_levels = 0; run_a2_residual = 0; run_a1 = 0; run_tracking_adaptive = 1; -plot_output_signals = 1; +plot_output_signals = options.debug_plots; %% Shared fixed EQ settings eq_sps = 2; @@ -42,11 +42,12 @@ end Scpe_sig = preprocessSignal(Scpe_sig_raw, Symbols, options.fsym, ... "mode", "auto", ... "debug_plots", options.debug_plots); +[Scpe_sig, Symbols, Tx_bits] = alignDspInputs(Scpe_sig, Symbols, Tx_bits, eq_sps); output = struct(); %% if plot_output_signals - showLevelScatter(Scpe_sig, Symbols, ... + showLevelScatter(Scpe_sig.normalize("mode","rms"), Symbols, ... "fsym", options.fsym, ... "fignum", 400, ... "normalize", true); @@ -283,9 +284,9 @@ if run_tracking_adaptive "dc_tracking_persistence_gain", 0, ... "dc_tracking_buffer_len", block_update, ... "optmize_mus", false, ... - "optimize_dc_tracking_params", false, ... + "optimize_dc_tracking_params", true, ... "dc_tracking_optimization_len", 2^15, ... - "dc_tracking_optimization_max_evals", 30, ... + "dc_tracking_optimization_max_evals", 20, ... "plot_mu_optimization", options.debug_plots, ... "save_debug", eq_save_debug); @@ -296,16 +297,18 @@ if run_tracking_adaptive "dc_tracking", "adaptive", block_update); output.(char(storageName)) = ffe_results; - fignum = 106; + if options.debug_plots + fignum = 106; - eq_noise = equalized_signal - Symbols; - dn = sprintf("FFE DCT; SIR: %d dB",options.dataTable.sir); - showEQNoisePSD(eq_noise, "fignum", fignum, "displayname", dn,"colormode","diverging"); - ylim([-70 -30]); - % - % dn = sprintf("FFE only; SIR: %d dB",options.dataTable.sir); - % showEQNoisePSD(eq_noise_ffe, "fignum", fignum+1, "displayname", dn,"colormode","diverging"); - % ylim([-70 -30]); + eq_noise = equalized_signal - Symbols; + dn = sprintf("FFE DCT; SIR: %d dB",options.dataTable.sir); + showEQNoisePSD(eq_noise, "fignum", fignum, "displayname", dn,"colormode","diverging"); + ylim([-70 -30]); + % + dn = sprintf("FFE only; SIR: %d dB",options.dataTable.sir); + showEQNoisePSD(eq_noise_ffe, "fignum", fignum+1, "displayname", dn,"colormode","diverging"); + ylim([-70 -30]); + end if plot_output_signals @@ -387,6 +390,22 @@ function description = resultDescription(prefix,options) description = sprintf('%s; SIR %g dB',prefix,sir); end +function [Scpe_sig, Symbols, Tx_bits] = alignDspInputs(Scpe_sig, Symbols, Tx_bits, sps) + nSymbols = min(length(Symbols), floor(length(Scpe_sig) / sps)); + if nSymbols <= 0 + error("mpi_recipe_dev:EmptyAlignedSignal", ... + "No overlapping samples remain after preprocessing and synchronization."); + end + + Scpe_sig.signal = Scpe_sig.signal(1:sps*nSymbols); + Scpe_sig.signal = real(Scpe_sig.signal); + Symbols.signal = Symbols.signal(1:nSymbols,:); + + if isprop(Tx_bits, "signal") + Tx_bits.signal = Tx_bits.signal(1:nSymbols,:); + end +end + function plotEqSignals(equalized_signal,Symbols,options,fignum,output_scale) showLevelScatter(equalized_signal .* output_scale, Symbols, ... diff --git a/Functions/EQ_visuals/showMarkovDiagram.m b/Functions/EQ_visuals/showMarkovDiagram.m index 49068e3..634d9b4 100644 --- a/Functions/EQ_visuals/showMarkovDiagram.m +++ b/Functions/EQ_visuals/showMarkovDiagram.m @@ -8,7 +8,8 @@ sequence = sequence(2+mod(1,length(sequence)):end); %filtered sequences often ha x = sequence; -levels = sort(unique(x)).'; % or provide known 1x6 level values +levels = sort(unique(x)).'; + [~,ix] = min(abs(x - levels),[],2); x = levels(ix); diff --git a/Functions/Metrics/calc_air.m b/Functions/Metrics/calc_air.m index 2b10827..dd6498e 100644 --- a/Functions/Metrics/calc_air.m +++ b/Functions/Metrics/calc_air.m @@ -29,8 +29,18 @@ end % CALC AIR %%% new implementation of AIR constellation = unique(reference_signal); +% map reference symbols to constellation indices reference_idx = arrayfun(@(x) find(constellation == x, 1), reference_signal); -ach_inf_rate = air_garcia_implementation(constellation',test_signal',reference_idx'); + +% compute probability mass function of symbol indices +[unique_idx,~,ic] = unique(reference_idx); +counts = accumarray(ic,1); +pmf = zeros(size(constellation)); +pmf(unique_idx) = counts / sum(counts); + +% ensure pmf corresponds to full set of constellation indices 1:N +% (pmf already aligned because unique_idx are indices into constellation) +ach_inf_rate = air_garcia_implementation(constellation',test_signal',reference_idx',pmf); function [data_,reference_]=trimseq(data,reference,skipstart,skip_end) diff --git a/Tests/04_DSP/Equalizer/ML_MLSE_DUOBINARY_test.m b/Tests/04_DSP/Equalizer/ML_MLSE_DUOBINARY_test.m new file mode 100644 index 0000000..0c1786c --- /dev/null +++ b/Tests/04_DSP/Equalizer/ML_MLSE_DUOBINARY_test.m @@ -0,0 +1,84 @@ +classdef ML_MLSE_DUOBINARY_test < IMDDTestCase + methods (Test, TestTags = {'unit', 'fast', 'dsp', 'ml_mlse', 'duobinary'}) + function constructorStoresDuobinaryConfiguration(testCase) + eq = ML_MLSE_DUOBINARY(); + + testCase.verifyEqual(class(eq), 'ML_MLSE_DUOBINARY'); + testCase.verifyEqual(eq.sps, 2); + testCase.verifyEqual(eq.order, 15); + testCase.verifyEqual(eq.ber, []); + testCase.verifyEqual(eq.ber_dd, []); + end + + function processReturnsPrecodeDomainAndFiniteDiagnostics(testCase) + [~, ~, txPrecoded, txEncoded] = makePam4Fixture(512); + + eq = ML_MLSE_DUOBINARY( ... + "sps", 1, ... + "order", 1, ... + "len_tr", length(txEncoded), ... + "epochs_tr", 1, ... + "epochs_dd", 1, ... + "mu_tr", 0.01, ... + "mu_dd", 0.01, ... + "adaptive_mu", false, ... + "L", 1); + + [detected, detectedViterbi] = eq.process(txEncoded, txPrecoded); + + testCase.verifyEqual(length(detected), length(txEncoded)); + testCase.verifyEqual(length(detectedViterbi), length(txEncoded)); + testCase.verifyTrue(all(ismembertol(detected.signal, unique(txPrecoded.signal), 1e-12))); + testCase.verifySize(eq.ber, [1 1]); + testCase.verifySize(eq.ber_dd, [1 1]); + testCase.verifyTrue(isfinite(eq.ber(1))); + testCase.verifyTrue(isfinite(eq.ber_dd(1))); + + end + + function encodeDecodeInvertsPrecodingAwayFromInitialState(testCase) + [~, txSymbols, txPrecoded, ~] = makePam4Fixture(512); + + encoded = Duobinary().encode(txPrecoded, "M", 4); + decoded = Duobinary().decode(encoded, "M", 4); + + testCase.verifyEqual(decoded.signal(11:end-10), ... + txSymbols.signal(11:end-10), "AbsTol", 1e-12); + end + + function resultWrapperReportsBothPrecodedAndOriginalBer(testCase) + [txBits, ~, txPrecoded, txEncoded] = makePam4Fixture(32000); + + eq = ML_MLSE_DUOBINARY( ... + "sps", 1, ... + "order", 1, ... + "len_tr", length(txEncoded), ... + "epochs_tr", 1, ... + "epochs_dd", 1, ... + "mu_tr", 0.01, ... + "mu_dd", 0.01, ... + "adaptive_mu", false, ... + "L", 1); + + results = ml_mlse(eq, 4, txEncoded, txPrecoded, txBits, ... + "precode_mode", db_mode.db_precoded); + + testCase.verifyTrue(isfinite(results.metrics.BER)); + testCase.verifyTrue(isfinite(results.metrics.BER_precoded)); + testCase.verifyGreaterThanOrEqual(results.metrics.numBits, 0); + testCase.verifyGreaterThanOrEqual(results.metrics.numBitErr, 0); + testCase.verifyGreaterThanOrEqual(results.metrics.numBitErr_precoded, 0); + end + end +end + +function [txBits, txSymbols, txPrecoded, txEncoded] = makePam4Fixture(nSymbols) + pattern = [0 0; 0 1; 1 1; 1 0]; + bits = repmat(pattern, ceil(nSymbols / size(pattern, 1)), 1); + bits = bits(1:nSymbols, :); + + txBits = Informationsignal(bits); + txSymbols = PAMmapper(4, 0).map(txBits); + txPrecoded = Duobinary().precode(txSymbols); + txEncoded = Duobinary().encode(txPrecoded, "M", 4); +end diff --git a/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/power_vs_wavelength.m b/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/power_vs_wavelength.m index d018b47..b7e38ad 100644 --- a/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/power_vs_wavelength.m +++ b/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/power_vs_wavelength.m @@ -5,12 +5,12 @@ db = DBHandler("dataBase", [dataBase], "type", database_type); fp = QueryFilter(); fp.where('power_state_info', 'pam_level','EQUALS', 4); -fp.where('power_state_info', 'db_mode','EQUALS', 1); +fp.where('power_state_info', 'db_mode','EQUALS', 0); % fp.where('power_state_info', 'fiber_length','EQUALS', 1); fp.where('power_state_info', 'is_mpi','EQUALS', 0); fields = db.getTableFieldNames('power_state_info'); -% [dataTable,~] = db.queryDB(fp, fields); +[dataTable,~] = db.queryDB(fp, fields); fiber_len = unique(dataTable.fiber_length); cnt = 0; diff --git a/projects/Diss/400G_revisit/PLOT_BER_BEST_ALGOS_WITH_DUOBINARY_SIGNALING.m b/projects/Diss/400G_revisit/PLOT_BER_BEST_ALGOS_WITH_DUOBINARY_SIGNALING.m new file mode 100644 index 0000000..7522ccf --- /dev/null +++ b/projects/Diss/400G_revisit/PLOT_BER_BEST_ALGOS_WITH_DUOBINARY_SIGNALING.m @@ -0,0 +1,343 @@ +%% 400G BER over bitrate: best normal algorithms plus duobinary signaling +% Normal algorithms are reduced to the best BER per gross rate across +% db_mode 0/1, pre-emphasis on/off, and BER/BER_precoded result variants. +% Duobinary signaling uses db_mode = 2 and only the sequence-detection BER +% stored in the BER field. + +clear; clc; + +%% 1) Query data + +selectedPamLevel = 8; +selectedFiberLengthKm = 10; +selectedWavelengthNm = 1310; +selectedRopAttenuation = 0; % set [] to use all ROP attenuation values +selectedIsMpi = 0; % set [] to use all entries + +normalDbModes = [double(db_mode.no_db), double(db_mode.db_precoded)]; +duobinaryDbMode = double(db_mode.db_encoded); + +maxBerForPlot = 0.5; +showRawEntries = false; +showBestLine = true; + +algoStyles = defaultAlgorithmStyles(); + +db = DBHandler( ... + "dataBase", "labor_highspeed", ... + "type", "mysql", ... + "server", "192.168.178.192", ... + "user", "silas", ... + "password", "silas"); +db.refresh(); + +fp = QueryFilter(); +fp.where('Runs', 'fiber_length', 'EQUALS', selectedFiberLengthKm); +fp.where('Runs', 'pam_level', 'EQUALS', selectedPamLevel); +fp.where('Runs', 'wavelength', 'EQUALS', selectedWavelengthNm); +if ~isempty(selectedRopAttenuation) + fp.where('Runs', 'rop_attenuation', 'EQUALS', selectedRopAttenuation); +end +if ~isempty(selectedIsMpi) + fp.where('Runs', 'is_mpi', 'EQUALS', selectedIsMpi); +end + +selectedFields = db.getTableFieldNames('dashboard_ungrouped_alltime'); +selectedFields = appendMissingFields(selectedFields, ... + {'Runs.precomp_amp'; 'Runs.is_mpi'}); +selectedFields = selectedFields(:); + +[rawData, query] = db.queryDB(fp, selectedFields); +disp(query); +fprintf("Fetched %d 400G result rows.\n", height(rawData)); + +%% 2) Clean data and build the five plotted curves + +data = rawData; +numericFields = ["result_id", "run_id", "eq_id", "bitrate", "grossrate", ... + "symbolrate", "pam_level", "wavelength", "fiber_length", "db_mode", ... + "rop_attenuation", "precomp_amp", "is_mpi", "numBits", "numBitErr", ... + "BER", "numBitErr_precoded", "BER_precoded", "STD", "STDrx", ... + "GMI", "AIR", "NGMI", "EVM", "Alpha"]; +for fieldIdx = 1:numel(numericFields) + fieldName = numericFields(fieldIdx); + if ismember(fieldName, string(data.Properties.VariableNames)) + data.(char(fieldName)) = numericColumn(data.(char(fieldName))); + end +end + +if ~ismember("precomp_amp", string(data.Properties.VariableNames)) + warning("plot_best_algos:NoPrecompAmp", ... + "Runs.precomp_amp was not returned. Falling back to pre_emphasis = (db_mode == 0)."); + data.pre_emphasis = data.db_mode == double(db_mode.no_db); +else + data.pre_emphasis = derivePreEmphasis(data.precomp_amp, data.db_mode); +end + +normalRows = data(ismember(data.db_mode, normalDbModes), :); +normalPlotData = buildNormalMetricRows(normalRows); +normalPlotData = normalPlotData(isfinite(normalPlotData.BER_plot) & ... + normalPlotData.BER_plot > 0 & normalPlotData.BER_plot < maxBerForPlot, :); + +duobinaryRows = data(data.db_mode == duobinaryDbMode, :); +if ismember("equalizer_structure", string(duobinaryRows.Properties.VariableNames)) + duobinaryRows = duobinaryRows( ... + equalizerMask(duobinaryRows.equalizer_structure, ... + equalizer_structure.db_encoded), :); +end +duobinaryPlotData = buildDuobinarySignalingRows(duobinaryRows); +duobinaryPlotData = duobinaryPlotData(isfinite(duobinaryPlotData.BER_plot) & ... + duobinaryPlotData.BER_plot > 0 & ... + duobinaryPlotData.BER_plot < maxBerForPlot, :); + +plotData = [normalPlotData; duobinaryPlotData]; +if isempty(plotData) + warning("plot_best_algos:NoRows", ... + "No rows remain after length/PAM/wavelength/BER filtering."); + return +end + +plotData.bitrate_Gbps = plotData.bitrate .* 1e-9; +plotData.grossrate_Gbps = plotData.grossrate .* 1e-9; + +fprintf("Remaining candidate BER rows: %d\n", height(plotData)); +disp(groupcounts(plotData, ["algorithm_key", "db_mode", "pre_emphasis", "precode"])); + +bestPlotData = bestBerByAlgorithmAndGrossRate(plotData); +fprintf("Keeping %d best-BER rows across algorithm/gross-rate groups.\n", ... + height(bestPlotData)); +disp(groupcounts(bestPlotData, "algorithm_key")); + +%% 3) Plot one figure with five lines + +availableStyles = algoStyles(hasAlgorithmRows(bestPlotData, algoStyles), :); +if isempty(availableStyles) + warning("plot_best_algos:NoSelectedAlgorithms", ... + "None of the configured algorithm styles match the queried rows."); + return +end + +fig = figure(); clf; +ax = axes(fig); hold(ax, "on"); + +for styleIdx = 1:height(availableStyles) + style = availableStyles(styleIdx, :); + rowMask = bestPlotData.algorithm_key == style.algorithm_key; + if ~any(rowMask) + continue + end + + algoData = sortrows(bestPlotData(rowMask, :), "grossrate_Gbps"); + + if showRawEntries + scatter(ax, algoData.grossrate_Gbps, algoData.BER_plot, ... + 9, ... + "Marker", ".", ... + "MarkerEdgeColor", style.color, ... + "MarkerFaceColor", style.color, ... + "MarkerEdgeAlpha", 0.25, ... + "MarkerFaceAlpha", 0.25, ... + "HandleVisibility", "off"); + end + + if showBestLine + plot(ax, algoData.grossrate_Gbps, algoData.BER_plot, ... + "LineStyle", style.lineStyle, ... + "Marker", style.marker, ... + "MarkerSize", 5, ... + "LineWidth", 1.5, ... + "Color", style.color, ... + "MarkerFaceColor", style.markerFaceColor, ... + "MarkerEdgeColor", style.color, ... + "DisplayName", style.name); + end +end + +yline(ax, [2.2e-4, 4.85e-3, 2e-2], ... + "LineWidth", 1, ... + "LineStyle", "--", ... + "Color", [0.25 0.25 0.25], ... + "HandleVisibility", "off"); + +title(ax, sprintf("PAM-%d, %.0f km, %.0f nm", ... + selectedPamLevel, selectedFiberLengthKm, selectedWavelengthNm)); +xlabel(ax, "Gross rate [Gb/s]"); +ylabel(ax, "BER"); +set(ax, "YScale", "log"); +ylim(ax, [1e-5, maxBerForPlot]); +grid(ax, "on"); +box(ax, "on"); + +xTicks = unique(bestPlotData.grossrate_Gbps(isfinite(bestPlotData.grossrate_Gbps))); +if ~isempty(xTicks) + xticks(ax, xTicks); + xlim(ax, [min(xTicks), max(xTicks)]); +end + +legend(ax, "Location", "best", "Interpreter", "none"); +if exist("beautifyBERplot", "file") + beautifyBERplot("logscale", true, "setcolors", false, ... + "setmarkers", false, "changemarkers", false); +end + +set(fig, "Position", 1e3 .* [0.1000 0.5500 0.7200 0.4200]); + +%% Local helpers + +function fields = appendMissingFields(fields, extraFields) +fields = cellstr(fields); +extraFields = cellstr(extraFields); +for idx = 1:numel(extraFields) + if ~any(strcmp(fields, extraFields{idx})) + fields{end+1, 1} = extraFields{idx}; %#ok + end +end +end + +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 styles = defaultAlgorithmStyles() +styles = table( ... + ["vnle"; ... + "vnle_pf_mlse"; ... + "vnle_db_mlse"; ... + "ml_mlse"; ... + "db_encoded"], ... + [equalizer_structure.vnle; ... + equalizer_structure.vnle_pf_mlse; ... + equalizer_structure.vnle_db_mlse; ... + equalizer_structure.ml_mlse; ... + equalizer_structure.db_encoded], ... + ["VNLE"; ... + "VNLE + PF + MLSE"; ... + "VNLE DBt. + MLSE"; ... + "ML pre-EQ + Viterbi"; ... + "Duobinary signaling"], ... + ["o"; "square"; "diamond"; "^"; "v"], ... + ["-"; "-"; "-"; "-"; "-"], ... + ["w"; "w"; "w"; "w"; "w"], ... + [clr.Paired.red; ... + clr.Paired.green; ... + clr.Paired.blue; ... + clr.Paired.purple; ... + clr.Paired.orange], ... + 'VariableNames', ["algorithm_key", "eq", "name", "marker", ... + "lineStyle", "markerFaceColor", "color"]); +end + +function preEmphasis = derivePreEmphasis(precompAmp, dbMode) +preEmphasis = false(size(dbMode)); + +validPrecomp = isfinite(precompAmp); +preEmphasis(validPrecomp) = precompAmp(validPrecomp) > -45; + +missingPrecomp = ~validPrecomp; +preEmphasis(missingPrecomp) = dbMode(missingPrecomp) == double(db_mode.no_db); +end + +function plotData = buildNormalMetricRows(data) +baseRows = data(isfinite(data.BER), :); +baseRows.precode = false(height(baseRows), 1); +baseRows.BER_plot = baseRows.BER; +baseRows.algorithm_key = algorithmKeyFromEqualizer(baseRows.equalizer_structure); +baseRows = baseRows(baseRows.algorithm_key ~= "", :); + +if ismember("BER_precoded", string(data.Properties.VariableNames)) + precodedRows = data(isfinite(data.BER_precoded), :); + precodedRows.precode = true(height(precodedRows), 1); + precodedRows.BER_plot = precodedRows.BER_precoded; + precodedRows.algorithm_key = algorithmKeyFromEqualizer( ... + precodedRows.equalizer_structure); + precodedRows = precodedRows(precodedRows.algorithm_key ~= "", :); + plotData = [baseRows; precodedRows]; +else + warning("plot_best_algos:NoPrecodedBer", ... + "BER_precoded was not returned. Plotting only BER rows for normal algorithms."); + plotData = baseRows; +end +end + +function plotData = buildDuobinarySignalingRows(data) +plotData = data(isfinite(data.BER), :); +plotData.precode = false(height(plotData), 1); +plotData.BER_plot = plotData.BER; +plotData.algorithm_key = repmat("db_encoded", height(plotData), 1); +end + +function algorithmKey = algorithmKeyFromEqualizer(equalizerColumn) +eqNumeric = equalizerNumeric(equalizerColumn); +algorithmKey = strings(size(eqNumeric)); + +algorithmKey(eqNumeric == enumValue(equalizer_structure.vnle)) = "vnle"; +algorithmKey(eqNumeric == enumValue(equalizer_structure.vnle_pf_mlse)) = ... + "vnle_pf_mlse"; +algorithmKey(eqNumeric == enumValue(equalizer_structure.vnle_db_mlse)) = ... + "vnle_db_mlse"; +algorithmKey(eqNumeric == enumValue(equalizer_structure.ml_mlse)) = "ml_mlse"; +end + +function mask = equalizerMask(equalizerColumn, eqValue) +eqNumeric = equalizerNumeric(equalizerColumn); +mask = eqNumeric == enumValue(eqValue); +end + +function eqNumeric = equalizerNumeric(equalizerColumn) +if isa(equalizerColumn, "equalizer_structure") + eqNumeric = double(equalizerColumn); +elseif isnumeric(equalizerColumn) + eqNumeric = double(equalizerColumn); +else + equalizerString = string(equalizerColumn); + eqNumeric = str2double(equalizerString); + + enumNames = ["vnle", "ffe", "dfe", "vnle_pf_mlse", ... + "vnle_db_mlse", "db_encoded", "ml_mlse"]; + enumValues = [ ... + enumValue(equalizer_structure.vnle), ... + enumValue(equalizer_structure.ffe), ... + enumValue(equalizer_structure.dfe), ... + enumValue(equalizer_structure.vnle_pf_mlse), ... + enumValue(equalizer_structure.vnle_db_mlse), ... + enumValue(equalizer_structure.db_encoded), ... + enumValue(equalizer_structure.ml_mlse)]; + + for idx = 1:numel(enumNames) + missingNumeric = isnan(eqNumeric); + eqNumeric(missingNumeric & equalizerString == enumNames(idx)) = ... + enumValues(idx); + end +end +end + +function value = enumValue(enumEntry) +value = double(enumEntry); +end + +function bestData = bestBerByAlgorithmAndGrossRate(data) +groupVars = ["algorithm_key", "grossrate_Gbps"]; +groupId = findgroups(data(:, groupVars)); +keepIdx = NaN(max(groupId), 1); + +for curGroup = 1:max(groupId) + rowIdx = find(groupId == curGroup); + [~, localBestIdx] = min(data.BER_plot(rowIdx)); + keepIdx(curGroup) = rowIdx(localBestIdx); +end + +bestData = sortrows(data(keepIdx, :), groupVars); +end + +function keep = hasAlgorithmRows(data, algoStyles) +keep = false(height(algoStyles), 1); +for idx = 1:height(algoStyles) + keep(idx) = any(data.algorithm_key == algoStyles.algorithm_key(idx)); +end +end diff --git a/projects/Diss/400G_revisit/PLOT_BER_VS_ALGO.m b/projects/Diss/400G_revisit/PLOT_BER_VS_ALGO.m new file mode 100644 index 0000000..eed6bc1 --- /dev/null +++ b/projects/Diss/400G_revisit/PLOT_BER_VS_ALGO.m @@ -0,0 +1,294 @@ +%% 400G BER over bitrate from labor_highspeed.dashboard_ungrouped_alltime +% 1) gather all BER entries for one PAM format and fiber length +% 2) derive pre-emphasis and precoding groups +% 3) plot one bitrate-vs-BER tile per equalizer structure + +clear; clc; + +%% 1) Gather data + +selectedPamLevel = 4; +selectedFiberLengthKm = 10; +selectedWavelength =1310; % set [] to use all wavelengths +selectedRopAttenuation = []; % set [] to use all ROP attenuation values +selectedIsMpi = []; % set [] to use all entries + +maxBerForPlot = 0.5; +showRawEntries = true; +showMedianLine = true; + +eqStyles = defaultEqualizerStyles(); +comboStyles = defaultCombinationStyles(); + +db = DBHandler( ... + "dataBase", "labor_highspeed", ... + "type", "mysql"); +db.refresh(); + +fp = QueryFilter(); +fp.where('Runs', 'fiber_length', 'EQUALS', selectedFiberLengthKm); +fp.where('Runs', 'pam_level', 'EQUALS', selectedPamLevel); +if ~isempty(selectedWavelength) + fp.where('Runs', 'wavelength', 'EQUALS', selectedWavelength); +end +if ~isempty(selectedRopAttenuation) + fp.where('Runs', 'rop_attenuation', 'EQUALS', selectedRopAttenuation); +end +if ~isempty(selectedIsMpi) + fp.where('Runs', 'is_mpi', 'EQUALS', selectedIsMpi); +end + +selectedFields = db.getTableFieldNames('dashboard_ungrouped_alltime'); +selectedFields = [selectedFields; {'Runs.precomp_amp'; 'Runs.is_mpi'}]; +selectedFields = selectedFields(:); + +[rawData, query] = db.queryDB(fp, selectedFields); +disp(query); +fprintf("Fetched %d 400G result rows.\n", height(rawData)); + +%% 2) Clean data and derive analysis groups + +data = rawData; +numericFields = ["result_id", "run_id", "eq_id", "bitrate", "grossrate", ... + "symbolrate", "pam_level", "wavelength", "fiber_length", "db_mode", ... + "rop_attenuation", "precomp_amp", "is_mpi", "numBits", "numBitErr", ... + "BER", "numBitErr_precoded", "BER_precoded", "STD", "STDrx", ... + "GMI", "AIR", "NGMI", "EVM", "Alpha"]; +for fieldIdx = 1:numel(numericFields) + fieldName = numericFields(fieldIdx); + if ismember(fieldName, string(data.Properties.VariableNames)) + data.(char(fieldName)) = numericColumn(data.(char(fieldName))); + end +end + +if ~ismember("precomp_amp", string(data.Properties.VariableNames)) + warning("analyze_db:NoPrecompAmp", ... + "Runs.precomp_amp was not returned. Falling back to pre_emphasis = (db_mode == 0)."); + data.pre_emphasis = data.db_mode == 0; +else + data.pre_emphasis = derivePreEmphasis(data.precomp_amp, data.db_mode); +end + +plotData = buildBerMetricRows(data); +plotData = plotData(isfinite(plotData.BER_plot) & ... + plotData.BER_plot > 0 & plotData.BER_plot < maxBerForPlot, :); + +if isempty(plotData) + warning("analyze_db:NoRows", ... + "No rows remain after fiber/PAM/BER filtering."); + return +end + +plotData.bitrate_Gbps = plotData.bitrate .* 1e-9; +plotData.grossrate_Gbps = plotData.grossrate .* 1e-9; + +fprintf("Remaining plotted BER rows: %d\n", height(plotData)); +disp(groupcounts(plotData, ["equalizer_structure", "pre_emphasis", "precode"])); + +bestPlotData = bestBerByBitrateAndGroup(plotData); +fprintf("Keeping %d best-BER rows across bitrate/EQ/pre-emphasis/precode groups.\n", ... + height(bestPlotData)); + +%% 3) Plot bitrate versus BER + +availableEqStyles = eqStyles(hasEqualizerRows(bestPlotData, eqStyles), :); +if isempty(availableEqStyles) + warning("analyze_db:NoSelectedEqualizers", ... + "None of the configured equalizer styles match the queried rows."); + return +end + +fig = figure(401); clf; +tiledlayout(1, height(availableEqStyles), ... + "TileSpacing", "compact", ... + "Padding", "compact"); + +for eqIdx = 1:height(availableEqStyles) + eqStyle = availableEqStyles(eqIdx, :); + ax = nexttile; hold(ax, "on"); + eqMask = equalizerMask(bestPlotData.equalizer_structure, eqStyle.eq); + + for comboIdx = 1:height(comboStyles) + comboStyle = comboStyles(comboIdx, :); + rowMask = eqMask & ... + bestPlotData.pre_emphasis == comboStyle.pre_emphasis & ... + bestPlotData.precode == comboStyle.precode; + + if ~any(rowMask) + continue + end + + comboData = sortrows(bestPlotData(rowMask, :), "bitrate_Gbps"); + comboColor = emphasisColor(eqStyle.color, comboStyle.pre_emphasis); + label = sprintf("%s, %s", eqStyle.name, comboStyle.name); + + if showRawEntries + scatter(ax, comboData.bitrate_Gbps, comboData.BER_plot, ... + 5, ... + "Marker", '.', ... + "MarkerEdgeColor", comboColor, ... + "MarkerFaceColor", comboColor, ... + "MarkerEdgeAlpha", 0.25, ... + "MarkerFaceAlpha", 0.25, ... + "HandleVisibility", "off"); + end + + if showMedianLine + summaryTable = summarizeBerByBitrate(comboData); + plot(ax, summaryTable.bitrate_Gbps, summaryTable.median_BER_plot, ... + "LineStyle", comboStyle.lineStyle, ... + "Marker", eqStyle.marker, ... + "MarkerSize", 4, ... + "LineWidth", 1.4, ... + "Color", comboColor, ... + "MarkerFaceColor", comboStyle.markerFaceColor, ... + "MarkerEdgeColor", comboColor, ... + "DisplayName", label); + end + end + + yline(ax, [2.2e-4, 4.85e-3, 2e-2], ... + "LineWidth", 1, ... + "LineStyle", "--", ... + "Color", [0.25 0.25 0.25], ... + "HandleVisibility", "off"); + + title(ax, sprintf("PAM-%d, %s", selectedPamLevel, eqStyle.name)); + xlabel(ax, "Gross rate [Gb/s]"); + ylabel(ax, "BER"); + set(ax, "YScale", "log"); + ylim(ax, [1e-5, maxBerForPlot]); + grid(ax, "on"); + box(ax, "on"); + + xTicks = unique(bestPlotData.bitrate_Gbps(isfinite(bestPlotData.bitrate_Gbps))); + if ~isempty(xTicks) + xticks(ax, xTicks); + xlim(ax, [min(xTicks), max(xTicks)]); + end + + legend(ax, "Location", "best", "Interpreter", "none"); + if exist("beautifyBERplot", "file") + beautifyBERplot("logscale", true, "setcolors", false, ... + "setmarkers", false, "changemarkers", false); + end +end + +set(fig, "Position", 1e3 .* [0.1000 0.5500 1.4113 0.3200]); + +%% 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 styles = defaultEqualizerStyles() + styles = table( ... + [equalizer_structure.vnle; ... + equalizer_structure.vnle_pf_mlse; ... + equalizer_structure.vnle_db_mlse; ... + equalizer_structure.ml_mlse], ... + ["VNLE"; ... + "VNLE + PF + MLSE"; ... + "VNLE DBt. + MLSE"; ... + "ML pre-EQ + Viterbi"], ... + ["o"; "square"; "diamond"; "^"], ... + [clr.Paired.red; ... + clr.Paired.green; ... + clr.Paired.blue; ... + clr.Paired.purple], ... + 'VariableNames', ["eq", "name", "marker", "color"]); +end + +function styles = defaultCombinationStyles() + styles = table( ... + [true; true; false; false], ... + [true; false; true; false], ... + ["w/ pre-emph., w/ precode"; ... + "w/ pre-emph., w/o precode"; ... + "w/o pre-emph., w/ precode"; ... + "w/o pre-emph., w/o precode"], ... + ["--"; "-"; "--"; "-"], ... + ["w"; "w"; "none"; "none"], ... + 'VariableNames', ["pre_emphasis", "precode", "name", ... + "lineStyle", "markerFaceColor"]); +end + +function preEmphasis = derivePreEmphasis(precompAmp, dbMode) + preEmphasis = false(size(dbMode)); + + validPrecomp = isfinite(precompAmp); + % In the 400G measurement scripts, -50 dB is the low/no-pre-emphasis + % setting, while -38/-37/-34 dB are the active pre-emphasis settings. + preEmphasis(validPrecomp) = precompAmp(validPrecomp) > -45; + + missingPrecomp = ~validPrecomp; + preEmphasis(missingPrecomp) = dbMode(missingPrecomp) == 0; +end + +function plotData = buildBerMetricRows(data) + baseRows = data(isfinite(data.BER), :); + baseRows.precode = false(height(baseRows), 1); + baseRows.BER_plot = baseRows.BER; + + if ismember("BER_precoded", string(data.Properties.VariableNames)) + precodedRows = data(isfinite(data.BER_precoded), :); + precodedRows.precode = true(height(precodedRows), 1); + precodedRows.BER_plot = precodedRows.BER_precoded; + plotData = [baseRows; precodedRows]; + else + warning("analyze_db:NoPrecodedBer", ... + "BER_precoded was not returned. Plotting only precode = 0 rows."); + plotData = baseRows; + end +end + +function mask = equalizerMask(equalizerColumn, eqValue) + if isa(equalizerColumn, "equalizer_structure") + mask = equalizerColumn == eqValue; + elseif isnumeric(equalizerColumn) + mask = double(equalizerColumn) == double(int32(eqValue)); + else + mask = string(equalizerColumn) == string(eqValue); + end +end + +function keep = hasEqualizerRows(data, eqStyles) + keep = false(height(eqStyles), 1); + for idx = 1:height(eqStyles) + keep(idx) = any(equalizerMask(data.equalizer_structure, eqStyles.eq(idx))); + end +end + +function summaryTable = summarizeBerByBitrate(data) + summaryTable = groupsummary(data, "bitrate_Gbps", "median", "BER_plot"); + summaryTable = sortrows(summaryTable, "bitrate_Gbps"); +end + +function bestData = bestBerByBitrateAndGroup(data) + groupVars = ["equalizer_structure", "pre_emphasis", "precode", "bitrate_Gbps"]; + groupId = findgroups(data(:, groupVars)); + keepIdx = NaN(max(groupId), 1); + + for curGroup = 1:max(groupId) + rowIdx = find(groupId == curGroup); + [~, localBestIdx] = min(data.BER_plot(rowIdx)); + keepIdx(curGroup) = rowIdx(localBestIdx); + end + + bestData = sortrows(data(keepIdx, :), groupVars); +end + +function color = emphasisColor(baseColor, preEmphasis) + if preEmphasis + color = 0.65 .* baseColor + 0.35; + else + color = baseColor; + end +end diff --git a/projects/Diss/400G_revisit/PLOT_DUOBINARY_DETECTION_BER_VS_BITRATE.m b/projects/Diss/400G_revisit/PLOT_DUOBINARY_DETECTION_BER_VS_BITRATE.m new file mode 100644 index 0000000..9c52776 --- /dev/null +++ b/projects/Diss/400G_revisit/PLOT_DUOBINARY_DETECTION_BER_VS_BITRATE.m @@ -0,0 +1,242 @@ +%% Duobinary transmission: BER over bitrate for detection algorithms +% DB transmission means the sequence was precoded and encoded at the Tx +% (Runs.db_mode = db_mode.db_encoded). In this stored-result convention: +% BER -> VNLE + MLSE +% BER_precoded -> VNLE + memoryless detection + +clear; clc; + +%% 1) Query data + +selectedPamLevels = [4, 6, 8]; +selectedFiberLengthKm = 10; +selectedWavelengthNm = 1310; +selectedRopAttenuation = 0; +selectedIsMpi = 0; +selectedDbMode = db_mode.db_encoded; +selectedEqualizerStructure = equalizer_structure.db_encoded; + +maxBerForPlot = 0.5; +showRawEntries = false; +showBestLine = true; + +db = DBHandler( ... + "dataBase", "labor_highspeed", ... + "type", "mysql", ... + "server", "192.168.178.192", ... + "user", "silas", ... + "password", "silas"); +db.refresh(); + +fp = QueryFilter(); +fp.where('Runs', 'fiber_length', 'EQUALS', selectedFiberLengthKm); +fp.where('Runs', 'wavelength', 'EQUALS', selectedWavelengthNm); +fp.where('Runs', 'rop_attenuation', 'EQUALS', selectedRopAttenuation); +fp.where('Runs', 'is_mpi', 'EQUALS', selectedIsMpi); +fp.where('Runs', 'db_mode', 'EQUALS', double(selectedDbMode)); + +selectedFields = db.getTableFieldNames('dashboard_ungrouped_alltime'); +selectedFields = appendMissingFields(selectedFields, {'Runs.is_mpi'}); +selectedFields = selectedFields(:); + +[rawData, query] = db.queryDB(fp, selectedFields); +disp(query); +fprintf("Fetched %d duobinary result rows.\n", height(rawData)); + +%% 2) Clean and reshape BER metrics + +data = rawData; +numericFields = ["result_id", "run_id", "eq_id", "bitrate", "grossrate", ... + "symbolrate", "pam_level", "wavelength", "fiber_length", "db_mode", ... + "rop_attenuation", "is_mpi", "numBits", "numBitErr", "BER", ... + "numBitErr_precoded", "BER_precoded", "STD", "STDrx", ... + "GMI", "AIR", "NGMI", "EVM", "Alpha"]; +for fieldIdx = 1:numel(numericFields) + fieldName = numericFields(fieldIdx); + if ismember(fieldName, string(data.Properties.VariableNames)) + data.(char(fieldName)) = numericColumn(data.(char(fieldName))); + end +end + +data = data(ismember(data.pam_level, selectedPamLevels), :); +if ismember("equalizer_structure", string(data.Properties.VariableNames)) && ... + ~isempty(selectedEqualizerStructure) + data = data(equalizerMask(data.equalizer_structure, selectedEqualizerStructure), :); +end + +plotData = buildDetectionMetricRows(data); +plotData = plotData(isfinite(plotData.BER_plot) & ... + plotData.BER_plot > 0 & plotData.BER_plot < maxBerForPlot, :); + +if isempty(plotData) + warning("plot_duobinary_detection:NoRows", ... + "No rows remain after duobinary/PAM/wavelength/BER filtering."); + return +end + +plotData.bitrate_Gbps = plotData.bitrate .* 1e-9; +plotData = sortrows(plotData, ... + ["pam_level", "wavelength", "detection_type", "bitrate_Gbps", "run_id"]); + +fprintf("Remaining plotted BER rows: %d\n", height(plotData)); +disp(groupcounts(plotData, ["pam_level", "wavelength", "detection_type"])); + +%% 3) Plot BER versus bitrate + +detectionStyles = defaultDetectionStyles(); +availablePamLevels = selectedPamLevels(ismember(selectedPamLevels, unique(plotData.pam_level).')); + +fig = figure(430); clf; +tiledlayout(1, numel(availablePamLevels), ... + "TileSpacing", "compact", ... + "Padding", "compact"); + +for pamIdx = 1:numel(availablePamLevels) + pamLevel = availablePamLevels(pamIdx); + ax = nexttile; hold(ax, "on"); + pamMask = plotData.pam_level == pamLevel; + + for styleIdx = 1:height(detectionStyles) + style = detectionStyles(styleIdx, :); + rowMask = pamMask & plotData.detection_type == style.detection_type; + if ~any(rowMask) + continue + end + + detectionData = sortrows(plotData(rowMask, :), "bitrate_Gbps"); + summaryTable = summarizeBerByBitrate(detectionData); + + if showRawEntries + scatter(ax, detectionData.bitrate_Gbps, detectionData.BER_plot, ... + 9, ... + "Marker", ".", ... + "MarkerEdgeColor", style.color, ... + "MarkerFaceColor", style.color, ... + "MarkerEdgeAlpha", 0.25, ... + "MarkerFaceAlpha", 0.25, ... + "HandleVisibility", "off"); + end + + if showBestLine + plot(ax, summaryTable.bitrate_Gbps, summaryTable.BER_plot, ... + "LineStyle", style.lineStyle, ... + "Marker", style.marker, ... + "MarkerSize", 5, ... + "LineWidth", 1.4, ... + "Color", style.color, ... + "MarkerFaceColor", style.markerFaceColor, ... + "MarkerEdgeColor", style.color, ... + "DisplayName", style.name); + end + end + + yline(ax, [2.2e-4, 4.85e-3, 2e-2], ... + "LineWidth", 1, ... + "LineStyle", "--", ... + "Color", [0.25 0.25 0.25], ... + "HandleVisibility", "off"); + + title(ax, sprintf("PAM-%d, %.0f nm", pamLevel, selectedWavelengthNm)); + xlabel(ax, "Bitrate [Gb/s]"); + ylabel(ax, "BER"); + set(ax, "YScale", "log"); + ylim(ax, [1e-5, maxBerForPlot]); + grid(ax, "on"); + box(ax, "on"); + + xTicks = unique(plotData.bitrate_Gbps(pamMask & isfinite(plotData.bitrate_Gbps))); + if ~isempty(xTicks) + xticks(ax, xTicks); + xlim(ax, [min(xTicks), max(xTicks)]); + end + + legend(ax, "Location", "best", "Interpreter", "none"); + if exist("beautifyBERplot", "file") + beautifyBERplot("logscale", true, "setcolors", false, ... + "setmarkers", false, "changemarkers", false); + end +end + +set(fig, "Position", 1e3 .* [0.1000 0.5500 1.4113 0.3200]); + +%% Local helpers + +function fields = appendMissingFields(fields, extraFields) +fields = cellstr(fields); +extraFields = cellstr(extraFields); +for idx = 1:numel(extraFields) + if ~any(strcmp(fields, extraFields{idx})) + fields{end+1, 1} = extraFields{idx}; %#ok + end +end +end + +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 plotData = buildDetectionMetricRows(data) +baseRows = data(isfinite(data.BER), :); +baseRows.detection_type = repmat("VNLE + MLSE", height(baseRows), 1); +baseRows.BER_plot = baseRows.BER; + +if ismember("BER_precoded", string(data.Properties.VariableNames)) + memorylessRows = data(isfinite(data.BER_precoded), :); + memorylessRows.detection_type = repmat("VNLE + memoryless", ... + height(memorylessRows), 1); + memorylessRows.BER_plot = memorylessRows.BER_precoded; + plotData = [baseRows; memorylessRows]; +else + warning("plot_duobinary_detection:NoPrecodedBer", ... + "BER_precoded was not returned. Plotting only VNLE + MLSE rows."); + plotData = baseRows; +end +end + +function styles = defaultDetectionStyles() +styles = table( ... + ["VNLE + MLSE"; "VNLE + memoryless"], ... + ["VNLE + MLSE"; "VNLE + memoryless"], ... + ["o"; "square"], ... + ["-"; "--"], ... + ["w"; "none"], ... + [clr.Paired.blue; clr.Paired.orange], ... + 'VariableNames', ["detection_type", "name", "marker", ... + "lineStyle", "markerFaceColor", "color"]); +end + +function mask = equalizerMask(equalizerColumn, eqValue) +if isa(equalizerColumn, "equalizer_structure") + mask = equalizerColumn == eqValue; +elseif isnumeric(equalizerColumn) + mask = double(equalizerColumn) == enumValue(eqValue); +else + equalizerString = string(equalizerColumn); + numericEqualizer = str2double(equalizerString); + mask = equalizerString == string(eqValue) | numericEqualizer == enumValue(eqValue); +end +end + +function value = enumValue(enumEntry) +value = double(enumEntry); +end + +function summaryTable = summarizeBerByBitrate(data) +groupId = findgroups(data.bitrate_Gbps); +keepIdx = NaN(max(groupId), 1); + +for curGroup = 1:max(groupId) + rowIdx = find(groupId == curGroup); + [~, localBestIdx] = min(data.BER_plot(rowIdx)); + keepIdx(curGroup) = rowIdx(localBestIdx); +end + +summaryTable = data(keepIdx, :); +summaryTable = sortrows(summaryTable, "bitrate_Gbps"); +end diff --git a/projects/Diss/400G_revisit/duobinary_partly_failed.mat b/projects/Diss/400G_revisit/duobinary_partly_failed.mat new file mode 100644 index 0000000..3979c55 Binary files /dev/null and b/projects/Diss/400G_revisit/duobinary_partly_failed.mat differ diff --git a/projects/Diss/400G_revisit/investigate_400g_algorithms.m b/projects/Diss/400G_revisit/investigate_400g_algorithms.m new file mode 100644 index 0000000..56d5856 --- /dev/null +++ b/projects/Diss/400G_revisit/investigate_400g_algorithms.m @@ -0,0 +1,361 @@ +% === 400G DSP settings === +dsp_options = struct(); +dsp_options.mode = "run_id"; +dsp_options.recipe = @dsp_400g_recipe; +dsp_options.append_to_db = false; +% dsp_options.append_mpi_reduction_db = false; +dsp_options.start_occurence = 1; +dsp_options.max_occurences = 1; +dsp_options.debug_plots = false; + +dsp_options.database_type = "mysql"; +dsp_options.dataBase = "labor_highspeed"; + +if ismac + dsp_options.storage_path = "/Volumes/media/labdata/sioe_labor"; +else + dsp_options.storage_path = "W:\labdata\sioe_labor"; +end +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); + +%% Select runs + +maxRunIds = 1; % keep small until the recipe settings are settled + +fp = QueryFilter(); +fp.where('Runs','fiber_length','EQUALS', 10); +fp.where('Runs','wavelength','EQUALS', 1310); +fp.where('Runs','bitrate','EQUALS', 330e9); +fp.where('Runs','pam_level','EQUALS', 4); +fp.where('Runs','rop_attenuation','EQUALS', 0); +fp.where('Runs','is_mpi','EQUALS', 0); +fp.where('Runs', 'db_mode','EQUALS', 2); + +fields = db.getTableFieldNames('Runs'); +[dataTable, query] = db.queryDB(fp, fields); +disp(query); + +dataTable = sortrows(dataTable, {'bitrate', 'run_id'}); + +% dataTable = dataTable(1:maxRunIds, :); + +run_ids = dataTable.run_id(:).'; + +if isempty(run_ids) + error("investigate_400g_algorithms:MissingRunIds", ... + "No 400G runs match the current filters."); +end + +fprintf("Selected %d run_id(s): %s\n", numel(run_ids), mat2str(run_ids)); + +%% Parameter sweep + +dsp_options.userParameters = struct(); +dsp_options.userParameters.run_ml_mlse_db = true; +dsp_options.userParameters.run_mlse_db = true; + +% Enable/disable equalizer branches. +% dsp_options.userParameters.run_ffe = false; +% dsp_options.userParameters.run_vnle = false; +% dsp_options.userParameters.run_dfe = false; +% dsp_options.userParameters.run_vnle_mlse = true; +% dsp_options.userParameters.run_dbtgt = true; + +% Examples for parameter loops. DataStorage expands every vector-valued field. +% dsp_options.userParameters.len_tr = 4096*2; +% dsp_options.userParameters.pf_ncoeffs = 1; + +% dsp_options.userParameters.decoding_mode = [db_decoder.memoryless,db_decoder.sequencedetection]; +% dsp_options.userParameters.pf_ncoeffs = [1, 2, 3]; +% dsp_options.userParameters.mu_dc = [0, 1e-5, 1e-4]; +% dsp_options.userParameters.run_ml_mlse_db = [false, true]; +% dsp_options.userParameters.run_mlse_db = [false, true]; + +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) x %d userParam combination(s) = %d job(s) ] x %d realizations = %d total jobs \n", ... + n_run_ids, n_userparams, parallel_jobs, n_realizations, queried_jobs); + +%% Run + +[results, wh] = submitJobs(run_ids, dsp_options, processingMode.serial, ... + "wh", wh, ... + "waitbar", true); + + +%% Quick result overview + +printBerSummary(wh); +plotBerVsBitrateQuick(wh, dataTable); + +function printBerSummary(wh) +storageNames = fieldnames(wh.sto); +if isempty(storageNames) + fprintf("No non-empty recipe outputs were stored.\n"); + return +end + +fprintf("\nBER summary by stored package:\n"); +for storageIdx = 1:numel(storageNames) + storageName = storageNames{storageIdx}; + values = wh.sto.(storageName)(:).'; + berValues = extractBerValues(values); + + if isempty(berValues) + fprintf(" %-18s no BER values\n", storageName); + else + fprintf(" %-18s min %.3e | median %.3e | n %d\n", ... + storageName, min(berValues), median(berValues), numel(berValues)); + end +end +end + +function berValues = extractBerValues(values) +berValues = []; +for valueIdx = 1:numel(values) + packageCell = values{valueIdx}; + if isempty(packageCell) + continue + end + if ~iscell(packageCell) + packageCell = {packageCell}; + end + + for packageIdx = 1:numel(packageCell) + package = packageCell{packageIdx}; + if isstruct(package) && isfield(package, "metrics") + metrics = package.metrics; + if isprop(metrics, "BER") + berValues(end+1) = metrics.BER; %#ok + elseif isstruct(metrics) && isfield(metrics, "BER") + berValues(end+1) = metrics.BER; %#ok + end + end + end +end +berValues = berValues(isfinite(berValues)); +end + +function plotBerVsBitrateQuick(wh, dataTable) +plotData = buildQuickBerTable(wh, dataTable); +if isempty(plotData) + fprintf("No BER values available for quick BER-vs-bitrate plot.\n"); + return +end + +storageNames = unique(plotData.storage_name, "stable"); +decodingModes = [db_decoder.memoryless, db_decoder.sequencedetection]; +berMetrics = ["BER", "BER_precoded"]; +lineStyles = ["-", ":"]; +rawMarkers = [".", "x"]; +markers = ["o", "square", "diamond", "^", "v", ">"]; + +fig = figure(402); clf; +ax = axes(fig); hold(ax, "on"); + +for storageIdx = 1:numel(storageNames) + storageName = storageNames(storageIdx); + for modeIdx = 1:numel(decodingModes) + decodingMode = decodingModes(modeIdx); + for metricIdx = 1:numel(berMetrics) + metricName = berMetrics(metricIdx); + rowMask = plotData.storage_name == storageName & ... + plotData.decoding_mode == decodingMode & ... + plotData.metric_name == metricName; + if ~any(rowMask) + continue + end + + modeData = sortrows(plotData(rowMask, :), "bitrate_Gbps"); + summaryData = groupsummary(modeData, "bitrate_Gbps", "median", "BER"); + summaryData = sortrows(summaryData, "bitrate_Gbps"); + color = quickPlotColor(modeIdx); + marker = markers(1 + mod(storageIdx - 1, numel(markers))); + label = sprintf("%s, %s, %s", storageName, ... + decodingModeLabel(decodingMode), metricName); + + scatter(ax, modeData.bitrate_Gbps, modeData.BER, ... + 12, ... + "Marker", rawMarkers(metricIdx), ... + "MarkerEdgeColor", color, ... + "MarkerEdgeAlpha", 0.25, ... + "HandleVisibility", "off"); + + plot(ax, summaryData.bitrate_Gbps, summaryData.median_BER, ... + "LineStyle", lineStyles(metricIdx), ... + "Marker", marker, ... + "MarkerSize", 5, ... + "LineWidth", 1.4, ... + "Color", color, ... + "DisplayName", label); + end + end +end + +yline(ax, [2.2e-4, 4.85e-3, 2e-2], ... + "LineWidth", 1, ... + "LineStyle", "--", ... + "Color", [0.25 0.25 0.25], ... + "HandleVisibility", "off"); + +xlabel(ax, "Bitrate [Gb/s]"); +ylabel(ax, "BER"); +title(ax, "Quick BER vs bitrate"); +set(ax, "YScale", "log"); +grid(ax, "on"); +box(ax, "on"); +legend(ax, "Location", "best", "Interpreter", "none"); + +if exist("beautifyBERplot", "file") + beautifyBERplot("logscale", true, "setcolors", false, ... + "setmarkers", false, "changemarkers", false); +end +end + +function plotData = buildQuickBerTable(wh, dataTable) +storageNames = fieldnames(wh.sto); +if isempty(storageNames) + plotData = table(); + return +end + +runIds = dataTable.run_id(:); +bitrates = dataTable.bitrate(:); + +storageCol = strings(0, 1); +runIdCol = zeros(0, 1); +bitrateCol = zeros(0, 1); +decodingCol = db_decoder.empty(0, 1); +metricCol = strings(0, 1); +berCol = zeros(0, 1); + +for storageIdx = 1:numel(storageNames) + storageName = storageNames{storageIdx}; + storageValues = wh.sto.(storageName); + for linIdx = 1:numel(storageValues) + [phys, storedValue] = wh.getPhysAndValueByLinIndex(storageName, linIdx); + metricRows = extractBerMetricRows({storedValue}); + if isempty(metricRows) || ~isfield(phys, "decoding_mode") + continue + end + + runId = resolveRunId(phys, runIds); + bitrate = resolveBitrate(runId, runIds, bitrates); + if ~isfinite(bitrate) + continue + end + + nRows = height(metricRows); + storageCol(end+1:end+nRows, 1) = string(storageName); + runIdCol(end+1:end+nRows, 1) = double(runId); + bitrateCol(end+1:end+nRows, 1) = double(bitrate); + decodingCol(end+1:end+nRows, 1) = phys.decoding_mode; + metricCol(end+1:end+nRows, 1) = metricRows.metric_name; + berCol(end+1:end+nRows, 1) = metricRows.BER; + end +end + +plotData = table(storageCol, runIdCol, bitrateCol, decodingCol, metricCol, berCol, ... + 'VariableNames', ["storage_name", "run_id", "bitrate", ... + "decoding_mode", "metric_name", "BER"]); +if ~isempty(plotData) + plotData = plotData(isfinite(plotData.BER) & plotData.BER > 0, :); + plotData.bitrate_Gbps = plotData.bitrate .* 1e-9; +end +end + +function metricRows = extractBerMetricRows(values) +metricNames = strings(0, 1); +berValues = zeros(0, 1); +requestedMetrics = ["BER", "BER_precoded"]; + +for valueIdx = 1:numel(values) + packageCell = values{valueIdx}; + if isempty(packageCell) + continue + end + if ~iscell(packageCell) + packageCell = {packageCell}; + end + + for packageIdx = 1:numel(packageCell) + package = packageCell{packageIdx}; + if ~isstruct(package) || ~isfield(package, "metrics") + continue + end + + metrics = package.metrics; + for metricIdx = 1:numel(requestedMetrics) + metricName = requestedMetrics(metricIdx); + value = readMetricValue(metrics, metricName); + if isfinite(value) + metricNames(end+1, 1) = metricName; %#ok + berValues(end+1, 1) = value; %#ok + end + end + end +end + +metricRows = table(metricNames, berValues, ... + 'VariableNames', ["metric_name", "BER"]); +end + +function value = readMetricValue(metrics, metricName) +value = NaN; +fieldName = char(metricName); +if isstruct(metrics) && isfield(metrics, fieldName) + value = metrics.(fieldName); +elseif isobject(metrics) && isprop(metrics, fieldName) + value = metrics.(fieldName); +end +end + +function runId = resolveRunId(phys, runIds) +if isfield(phys, "run_id") + runId = phys.run_id; +else + runId = runIds(1); +end +end + +function bitrate = resolveBitrate(runId, runIds, bitrates) +rowIdx = find(double(runIds) == double(runId), 1, "first"); +if isempty(rowIdx) + bitrate = NaN; +else + bitrate = bitrates(rowIdx); +end +end + +function color = quickPlotColor(modeIdx) +colors = [ ... + 0.1059 0.6196 0.4667; ... + 0.8510 0.3725 0.0078]; +color = colors(1 + mod(modeIdx - 1, size(colors, 1)), :); +end + +function label = decodingModeLabel(decodingMode) +switch decodingMode + case db_decoder.memoryless + label = "memoryless"; + case db_decoder.sequencedetection + label = "sequence detection"; + otherwise + label = string(decodingMode); +end +end diff --git a/projects/Diss/MPI_revisit/algorithms/PLOT_mpi_reduction_db_delay_regimes.m b/projects/Diss/MPI_revisit/algorithms/PLOT_mpi_reduction_db_delay_regimes.m index dfb99a3..007b4c4 100644 --- a/projects/Diss/MPI_revisit/algorithms/PLOT_mpi_reduction_db_delay_regimes.m +++ b/projects/Diss/MPI_revisit/algorithms/PLOT_mpi_reduction_db_delay_regimes.m @@ -9,7 +9,7 @@ clear; clc; studyName = "block_update_sweep"; selectedBlockUpdate = 1; % set [] to pool all block_update values -selectedPamLevels = 6; % set [] to use all PAM levels in the query result +selectedPamLevels = 4; % set [] to use all PAM levels in the query result algorithmSelection = table( ... ["plain_ffe"; ... @@ -17,7 +17,7 @@ algorithmSelection = table( ... "a2_residual"; ... "a1_moving_average"; ... "dc_tracking"], ... - [true; true; true; true; true], ... + [true; false; false; false; true], ... 'VariableNames', ["algorithm", "enabled"]); selectedAlgorithms = algorithmSelection.algorithm(algorithmSelection.enabled); @@ -286,6 +286,81 @@ for regimeIdx = 1:numel(regimeNames) % mat2tikz_improved("C:/Users/Silas/Documents/6971e0b65b380ca6d71c837f/04_Experimental_Evaluation/tikz/mpi/ber_vs_sir_" + regimeName + ".tikz","cleanfigure",1); end +%% 4) Plot BER spread over SIR by delay/coherence regime and algorithm + +spreadTable = berSpreadSummary(cleanData, groupVars); +fprintf("Calculated BER spread for %d regime/PAM/algorithm/SIR groups.\n", ... + height(spreadTable)); + +figure(); clf; hold on; + +for regimeIdx = 1:numel(regimeNames) + regimeName = regimeNames(regimeIdx); + lineStyle = pathRegimeLineStyle(regimeName); + + for algIdx = 1:numel(selectedAlgorithms) + algorithmName = selectedAlgorithms(algIdx); + algColor = algorithmColor(algorithmName); + marker = algorithmMarker(algorithmName, algorithmMarkers); + displayName = algorithmDisplayName(algorithmName); + + curveMask = spreadTable.path_regime == regimeName & ... + spreadTable.algorithm == algorithmName; + + if ~any(curveMask) + continue + end + + curveTable = sortrows(spreadTable(curveMask, :), ... + ["pam_level", "sir_exact"]); + + for pamIdx = 1:numel(selectedPamLevels) + pamLevel = selectedPamLevels(pamIdx); + pamMask = curveTable.pam_level == pamLevel; + if ~any(pamMask) + continue + end + + x = curveTable.sir_exact(pamMask).'; + y = curveTable.std_log10_BER(pamMask).'; + valid = isfinite(x) & isfinite(y); + + if ~any(valid) + continue + end + + if isscalar(selectedPamLevels) + legendText = sprintf("%s, %s", displayName, regimeName); + else + legendText = sprintf("%s, %s, PAM %.0f", ... + displayName, regimeName, pamLevel); + end + + plot(x(valid), y(valid), ... + "LineStyle", lineStyle, ... + "Marker", marker, ... + "MarkerSize", 3.5, ... + "LineWidth", 1, ... + "Color", algColor, ... + "MarkerFaceColor", "w", ... + "MarkerEdgeColor", algColor, ... + "DisplayName", legendText); + end + end +end + +xlabel("SIR (dB)"); +ylabel("Std. dev. of log_{10}(BER)"); +title("BER spread over SIR by delay regime and algorithm"); +xlim([15, 45]); +grid on; +box on; +legend("Location", "northeast", "Interpreter", "none"); + +if exist("beautifyBERplot", "file") + beautifyBERplot("logscale", false, "setcolors", false, "setmarkers", false); +end + %% Local helpers function values = numericColumn(values) @@ -487,6 +562,38 @@ function yFit = fitLogBer(x, y, xFit, fitOrder) yFit = 10 .^ polyval(coeff, xFit); end +function spreadTable = berSpreadSummary(cleanData, groupVars) + summaryGroups = groupsummary(cleanData, groupVars); + sirExactTable = groupsummary(cleanData, groupVars, "median", "sir_exact"); + summaryGroups = sortrows(summaryGroups, groupVars); + sirExactTable = sortrows(sirExactTable, groupVars); + summaryGroups.sir_exact = sirExactTable.median_sir_exact; + + stdLogBer = NaN(height(summaryGroups), 1); + stdBer = NaN(height(summaryGroups), 1); + for groupIdx = 1:height(summaryGroups) + rowMask = true(height(cleanData), 1); + for varIdx = 1:numel(groupVars) + varName = groupVars(varIdx); + rowMask = rowMask & cleanData.(char(varName)) == ... + summaryGroups.(char(varName))(groupIdx); + end + + berValues = cleanData.BER(rowMask); + berValues = berValues(isfinite(berValues) & berValues > 0); + if isempty(berValues) + continue + end + + stdLogBer(groupIdx) = std(log10(berValues), 0, "omitnan"); + stdBer(groupIdx) = std(berValues, 0, "omitnan"); + end + + spreadTable = summaryGroups; + spreadTable.std_log10_BER = stdLogBer; + spreadTable.std_BER = stdBer; +end + function label = algorithmDisplayName(algorithmName) algorithmName = string(algorithmName); switch algorithmName @@ -532,3 +639,18 @@ function marker = algorithmMarker(algorithmName, algorithmMarkers) end marker = algorithmMarkers{mod(markerIdx - 1, numel(algorithmMarkers)) + 1}; end + +function lineStyle = pathRegimeLineStyle(regimeName) + switch string(regimeName) + case "0-1 m" + lineStyle = "-"; + case "10-100 m" + lineStyle = "--"; + case "300 m" + lineStyle = ":"; + case "1000 m" + lineStyle = "-."; + otherwise + lineStyle = "-"; + end +end diff --git a/projects/Diss/MPI_revisit/parallelization_analysis/PLOT_mpi_reduction_db_parallelization_vs_sir.m b/projects/Diss/MPI_revisit/parallelization_analysis/PLOT_mpi_reduction_db_parallelization_vs_sir.m index a7ca35f..dff72c7 100644 --- a/projects/Diss/MPI_revisit/parallelization_analysis/PLOT_mpi_reduction_db_parallelization_vs_sir.m +++ b/projects/Diss/MPI_revisit/parallelization_analysis/PLOT_mpi_reduction_db_parallelization_vs_sir.m @@ -132,7 +132,7 @@ summaryTable = addBerStdBounds(summaryTable, cleanData, groupVars); blockUpdates = unique(cleanData.block_update(isfinite(cleanData.block_update))).'; blockUpdates = sort(blockUpdates); -helperBlockUpdates = unique(blockUpdates([1, end]), "stable"); +helperBlockUpdates = unique(blockUpdates([1, end-3]), "stable"); requiredSirRows = table(); for pamIdx = 1:numel(selectedPamLevels) @@ -214,14 +214,14 @@ for pamIdx = 1:numel(selectedPamLevels) sirValues, boundCenterBer, boundLowerBer, boundUpperBer, ... boundMode, boundaryPolyfitOrderMax); - [hl, hp] = boundedline(xBand, centerBand, yBounds, ... - 'alpha', 'transparency', 0.1, ... - 'cmap', algColor, ... - 'nan', 'fill', ... - 'orientation', 'vert'); - set(hl, "LineStyle", "none", 'LineWidth', 1, "Marker", "none", ... - "HandleVisibility", "off", "DisplayName", char(displayName)); - set(hp, "LineStyle", "-", "HandleVisibility", "off", "Marker", "none"); + % [hl, hp] = boundedline(xBand, centerBand, yBounds, ... + % 'alpha', 'transparency', 0.1, ... + % 'cmap', algColor, ... + % 'nan', 'fill', ... + % 'orientation', 'vert'); + % set(hl, "LineStyle", "none", 'LineWidth', 1, "Marker", "none", ... + % "HandleVisibility", "off", "DisplayName", char(displayName)); + % set(hp, "LineStyle", "-", "HandleVisibility", "off", "Marker", "none"); end plot(sirValues(valid), meanBer(valid), ... @@ -242,11 +242,11 @@ for pamIdx = 1:numel(selectedPamLevels) 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"); + plot(xFit, yFit, ... + "LineStyle", "--", ... + "LineWidth", 1.1, ... + "Color", algColor, ... + "HandleVisibility", "off"); end end end diff --git a/projects/Diss/MPI_revisit/parallelization_analysis/parallalization_algos_at_1.fig b/projects/Diss/MPI_revisit/parallelization_analysis/parallalization_algos_at_1.fig new file mode 100644 index 0000000..45561f9 Binary files /dev/null and b/projects/Diss/MPI_revisit/parallelization_analysis/parallalization_algos_at_1.fig differ diff --git a/projects/Diss/MPI_revisit/parallelization_analysis/parallalization_algos_at_2048.fig b/projects/Diss/MPI_revisit/parallelization_analysis/parallalization_algos_at_2048.fig new file mode 100644 index 0000000..9b561c8 Binary files /dev/null and b/projects/Diss/MPI_revisit/parallelization_analysis/parallalization_algos_at_2048.fig differ diff --git a/projects/Diss/MPI_revisit/parallelization_analysis/parallalization_analysis.fig b/projects/Diss/MPI_revisit/parallelization_analysis/parallalization_analysis.fig new file mode 100644 index 0000000..7f9fa22 Binary files /dev/null and b/projects/Diss/MPI_revisit/parallelization_analysis/parallalization_analysis.fig differ diff --git a/projects/Diss/MPI_revisit/simulation/PLOT_mpi_simulation_linewidth_vs_sir.m b/projects/Diss/MPI_revisit/simulation/PLOT_mpi_simulation_linewidth_vs_sir.m new file mode 100644 index 0000000..ff25551 --- /dev/null +++ b/projects/Diss/MPI_revisit/simulation/PLOT_mpi_simulation_linewidth_vs_sir.m @@ -0,0 +1,608 @@ +%% BER over SIR from saved MPI simulation warehouses grouped by linewidth +% Loads a saved simulation warehouse and plots one BER-vs-SIR curve per +% linewidth for each algorithm. + +% clear; +% clc; + +%% Load data +resultFile = ""; +if resultFile == "" + resultDir = fullfile(fileparts(mfilename("fullpath")), "results"); + files = dir(fullfile(resultDir, "mpi_simulation_*.mat")); + if isempty(files) + error("PLOT_mpi_simulation_linewidth_vs_sir:NoResultFiles", ... + "No mpi_simulation_*.mat files found in %s.", resultDir); + end + [~, newestIdx] = max([files.datenum]); + resultFile = fullfile(files(newestIdx).folder, files(newestIdx).name); +end + +loaded = load(resultFile, "wh", "simulation_config"); +wh = loaded.wh; +fprintf("Loaded MPI simulation warehouse:\n%s\n", resultFile); +wh.showInfo; + +%% Plot settings +useBoundedLines = true; +usePolyfit = true; +polyfitOrderMax = 4; +boundaryPolyfitOrderMax = 4; +fecBerThreshold = 3.8e-3; +maxBerForPlot = 0.1; +crossingSirWindow = [15 40]; +boundMode = "fitStd"; +algorithmMarkers = {'o','square','diamond','^','v','>','<','pentagram'}; + +%% Collect and clean +cleanData = collectMpiSimulationBerRows(wh); +cleanData = normalizeBerColumnName(cleanData); +cleanData = cleanData(isfinite(cleanData.BER) & cleanData.BER > 0 & ... + cleanData.BER < maxBerForPlot, :); + +if isempty(cleanData) + warning("PLOT_mpi_simulation_linewidth_vs_sir:NoRows", ... + "No valid BER rows remain for plotting."); + return +end + +if all(~isfinite(cleanData.laser_linewidth)) + if isfield(loaded, "simulation_config") && isfield(loaded.simulation_config, "laser_linewidth") + cleanData.laser_linewidth(:) = loaded.simulation_config.laser_linewidth; + else + cleanData.laser_linewidth(:) = 0; + end +end + +cleanData.clean_keep = true(height(cleanData), 1); +groupId = findgroups(cleanData.storage_name, cleanData.block_update, ... + cleanData.laser_linewidth, cleanData.sir); +for curGroup = unique(groupId(isfinite(groupId))).' + rowMask = groupId == curGroup; + berValues = cleanData.BER(rowMask); + if nnz(rowMask) > 3 + cleanData.clean_keep(rowMask) = ~isoutlier(berValues); + end +end +cleanData = cleanData(cleanData.clean_keep, :); + +groupVars = ["storage_name", "algorithm", "block_update", "laser_linewidth", "sir"]; +summaryTable = groupsummary(cleanData, groupVars, {"mean", "min", "max"}, "BER"); +summaryTable = sortrows(summaryTable, groupVars); +summaryTable.sir_exact = summaryTable.sir; +summaryTable = addBerStdBounds(summaryTable, cleanData, groupVars); + +selectedBlockUpdates = unique(cleanData.block_update(isfinite(cleanData.block_update))).'; +selectedAlgorithms = unique(cleanData.storage_name, "stable").'; +selectedLinewidths = unique(cleanData.laser_linewidth(isfinite(cleanData.laser_linewidth))).'; +selectedLinewidths = sort(selectedLinewidths); + +requiredSirRows = table(); +for blockIdx = 1:numel(selectedBlockUpdates) + blockUpdate = selectedBlockUpdates(blockIdx); + for algIdx = 1:numel(selectedAlgorithms) + storageName = selectedAlgorithms(algIdx); + algorithmName = string(summaryTable.algorithm(find(summaryTable.storage_name == storageName, 1, "first"))); + for linewidthIdx = 1:numel(selectedLinewidths) + laserLinewidth = selectedLinewidths(linewidthIdx); + curveMask = summaryTable.storage_name == storageName & ... + summaryTable.block_update == blockUpdate & ... + summaryTable.laser_linewidth == laserLinewidth; + + sirValues = summaryTable.sir_exact(curveMask).'; + meanBer = summaryTable.mean_BER(curveMask).'; + [~, ~, requiredSir, fitOrder] = fitBerAtFec( ... + sirValues, meanBer, polyfitOrderMax, fecBerThreshold, crossingSirWindow); + + newRow = table(storageName, algorithmName, blockUpdate, ... + laserLinewidth, requiredSir, fitOrder, nnz(isfinite(sirValues) & isfinite(meanBer)), ... + 'VariableNames', {'storage_name', 'algorithm', 'block_update', ... + 'laser_linewidth', 'required_sir', 'fit_order', 'n_points'}); + requiredSirRows = [requiredSirRows; newRow]; %#ok + end + end +end + +fprintf("Cleaned to %d simulation BER rows across %d storage/block/linewidth/SIR groups.\n", ... + height(cleanData), height(summaryTable)); +disp(groupcounts(cleanData, ["storage_name", "block_update", "laser_linewidth"])); +disp(requiredSirRows); + +%% Plot BER vs SIR, linewidth as curve family +for blockIdx = 1:numel(selectedBlockUpdates) + blockUpdate = selectedBlockUpdates(blockIdx); + + figure(); + clf; + tiledlayout(numel(selectedAlgorithms), 1, "TileSpacing", "compact"); + + for algIdx = 1:numel(selectedAlgorithms) + storageName = selectedAlgorithms(algIdx); + nexttile; + hold on; + + for linewidthIdx = 1:numel(selectedLinewidths) + laserLinewidth = selectedLinewidths(linewidthIdx); + lineColor = linewidthColor(linewidthIdx, numel(selectedLinewidths)); + marker = algorithmMarker(storageName, algorithmMarkers); + displayName = sprintf("%s, %s", ... + algorithmDisplayName(storageName), linewidthLabel(laserLinewidth)); + + rawMask = cleanData.storage_name == storageName & ... + cleanData.block_update == blockUpdate & ... + cleanData.laser_linewidth == laserLinewidth; + curveMask = summaryTable.storage_name == storageName & ... + summaryTable.block_update == blockUpdate & ... + summaryTable.laser_linewidth == laserLinewidth; + + if ~any(curveMask) + continue + end + + scatterRows = cleanData(rawMask, :); + scatter(scatterRows.sir, scatterRows.BER, ... + 22, ... + "Marker", ".", ... + "MarkerEdgeColor", lineColor, ... + "MarkerFaceColor", lineColor, ... + "HandleVisibility", "off"); + + sirValues = summaryTable.sir_exact(curveMask).'; + meanBer = summaryTable.mean_BER(curveMask).'; + boundCenterBer = summaryTable.std_center_BER(curveMask).'; + boundLowerBer = summaryTable.std_lower_BER(curveMask).'; + boundUpperBer = summaryTable.std_upper_BER(curveMask).'; + valid = isfinite(sirValues) & isfinite(meanBer) & meanBer > 0; + + if useBoundedLines && exist("boundedline", "file") && any(valid) + [xBand, centerBand, yBounds] = berStdBounds( ... + sirValues, boundCenterBer, boundLowerBer, boundUpperBer, ... + boundMode, boundaryPolyfitOrderMax); + + [hl, hp] = boundedline(xBand, centerBand, yBounds, ... + 'alpha', 'transparency', 0.08, ... + 'cmap', lineColor, ... + 'nan', 'fill', ... + 'orientation', 'vert'); + set(hl, "LineStyle", "none", "LineWidth", 1, "Marker", "none", ... + "HandleVisibility", "off", "DisplayName", displayName); + set(hp, "LineStyle", "-", "HandleVisibility", "off", "Marker", "none"); + end + + plot(sirValues(valid), meanBer(valid), ... + "LineStyle", "-", ... + "Marker", marker, ... + "MarkerSize", 3, ... + "LineWidth", 1, ... + "Color", lineColor, ... + "MarkerFaceColor", "w", ... + "MarkerEdgeColor", lineColor, ... + "DisplayName", displayName, ... + "HandleVisibility", "on"); + + 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", lineColor, ... + "HandleVisibility", "off"); + end + end + end + + yline(2.2e-4, "LineWidth", 1, "LineStyle", "--", "HandleVisibility", "off"); + yline(fecBerThreshold, "LineWidth", 1, "LineStyle", "--", "HandleVisibility", "off"); + yline(2e-2, "LineWidth", 1, "LineStyle", "--", "HandleVisibility", "off"); + + title(sprintf("%s, block update %.0f", algorithmDisplayName(storageName), blockUpdate)); + xlabel("SIR (dB)"); + ylabel("BER"); + set(gca, "YScale", "log"); + ylim([9e-5, maxBerForPlot]); + xlim(crossingSirWindow); + grid on; + box on; + legend("Location", "northeast", "Interpreter", "none"); + + if exist("beautifyBERplot", "file") + beautifyBERplot("logscale", true, "setcolors", false, ... + "setmarkers", false, "changemarkers", false); + end + end +end + +%% Plot required SIR at FEC over linewidth +for blockIdx = 1:numel(selectedBlockUpdates) + blockUpdate = selectedBlockUpdates(blockIdx); + + figure(); + clf; + hold on; + + for algIdx = 1:numel(selectedAlgorithms) + storageName = selectedAlgorithms(algIdx); + algColor = algorithmColor(storageName); + marker = algorithmMarker(storageName, algorithmMarkers); + displayName = algorithmDisplayName(storageName); + rowMask = requiredSirRows.storage_name == storageName & ... + requiredSirRows.block_update == blockUpdate; + + x = requiredSirRows.laser_linewidth(rowMask).'; + y = requiredSirRows.required_sir(rowMask).'; + [x, sortIdx] = sort(x); + y = y(sortIdx); + valid = isfinite(x) & isfinite(y); + + if ~any(valid) + continue + end + + plot(x(valid), y(valid), ... + "LineWidth", 1.4, ... + "LineStyle", "-", ... + "Marker", marker, ... + "MarkerSize", 5, ... + "Color", algColor, ... + "MarkerFaceColor", "w", ... + "MarkerEdgeColor", algColor, ... + "DisplayName", char(displayName)); + end + + set(gca, "XScale", "log"); + xticks(selectedLinewidths); + xticklabels(arrayfun(@linewidthLabel, selectedLinewidths, "UniformOutput", false)); + ylim(crossingSirWindow); + grid on; + box on; + xlabel("Laser linewidth"); + ylabel(sprintf("Required SIR at BER = %.1e (dB)", fecBerThreshold)); + title(sprintf("Required SIR over linewidth, block update %.0f", blockUpdate)); + legend("Location", "best", "Interpreter", "none"); + + if exist("beautifyBERplot", "file") + beautifyBERplot("logscale", false, "setcolors", false, ... + "setmarkers", false, "changemarkers", false); + end +end + +%% Local helpers +function data = collectMpiSimulationBerRows(wh) +storageNames = string(fieldnames(wh.sto)); +data = table(); + +for storageIdx = 1:numel(storageNames) + storageName = storageNames(storageIdx); + storage = wh.sto.(char(storageName)); + + for linIdx = 1:numel(storage) + package = storage{linIdx}; + ber = extractPackageBer(package); + if ~isfinite(ber) + continue + end + + [physValues, physNames] = wh.getPhysIndicesByLinIndex(linIdx); + phys = struct(); + for physIdx = 1:numel(physNames) + phys.(char(physNames{physIdx})) = physValues{physIdx}; + end + + algorithm = extractPackageAlgorithm(package, storageName); + newRow = table( ... + storageName, ... + algorithm, ... + readPhysValue(phys, "sir", NaN), ... + readPhysValue(phys, "block_update", NaN), ... + readPhysValue(phys, "random_key", NaN), ... + readPhysValue(phys, "laser_linewidth", NaN), ... + ber, ... + 'VariableNames', {'storage_name', 'algorithm', 'sir', ... + 'block_update', 'random_key', 'laser_linewidth', 'BER'}); + data = [data; newRow]; %#ok + end +end +end + +function data = normalizeBerColumnName(data) +variableNames = string(data.Properties.VariableNames); +if ismember("BER", variableNames) + return +end + +if ismember("ber", variableNames) + data.Properties.VariableNames(variableNames == "ber") = {'BER'}; + return +end + +error("PLOT_mpi_simulation_linewidth_vs_sir:MissingBerColumn", ... + "Could not find a BER or ber column in the simulation BER table."); +end + +function value = readPhysValue(phys, name, defaultValue) +if isfield(phys, name) + value = phys.(name); +else + value = defaultValue; +end +end + +function ber = extractPackageBer(package) +ber = NaN; +if isempty(package) + return +end + +if iscell(package) + package = package{1}; +end + +if isstruct(package) && isfield(package, "metrics") + metrics = package.metrics; + if isstruct(metrics) && isfield(metrics, "BER") + ber = metrics.BER; + elseif isobject(metrics) && isprop(metrics, "BER") + ber = metrics.BER; + end +end +end + +function algorithm = extractPackageAlgorithm(package, fallbackName) +algorithm = fallbackName; +if isempty(package) + return +end + +if iscell(package) + package = package{1}; +end + +if isstruct(package) && isfield(package, "mpi_reduction_config") && ... + isfield(package.mpi_reduction_config, "algorithm") + algorithm = string(package.mpi_reduction_config.algorithm); +end +end + +function summaryTable = addBerStdBounds(summaryTable, cleanData, groupVars) +nGroups = height(summaryTable); +stdCenter = NaN(nGroups, 1); +stdLower = NaN(nGroups, 1); +stdUpper = NaN(nGroups, 1); + +for groupIdx = 1:nGroups + rowMask = true(height(cleanData), 1); + for varIdx = 1:numel(groupVars) + varName = groupVars(varIdx); + rowMask = rowMask & cleanData.(char(varName)) == summaryTable.(char(varName))(groupIdx); + end + + [stdCenter(groupIdx), stdLower(groupIdx), stdUpper(groupIdx)] = ... + logBerMeanStdInterval(cleanData.BER(rowMask)); +end + +summaryTable.std_center_BER = stdCenter; +summaryTable.std_lower_BER = stdLower; +summaryTable.std_upper_BER = stdUpper; +end + +function [centerBer, lowerBer, upperBer] = logBerMeanStdInterval(berValues) +berValues = berValues(isfinite(berValues) & berValues > 0); +if isempty(berValues) + centerBer = NaN; + lowerBer = NaN; + upperBer = NaN; + return +end + +logBer = log10(berValues(:)); +centerLog = mean(logBer, "omitnan"); +stdLog = std(logBer, 0, "omitnan"); + +centerBer = 10 .^ centerLog; +lowerBer = 10 .^ (centerLog - stdLog); +upperBer = 10 .^ (centerLog + stdLog); +end + +function [xBand, centerBand, yBounds] = berStdBounds(sirValues, centerBer, lowerBer, upperBer, boundMode, maxOrder) +if boundMode == "directStd" + [xBand, centerBand, yBounds] = directBerBounds(sirValues, centerBer, lowerBer, upperBer); + return +end + +[xBand, centerBand, yBounds] = fittedBerBounds(sirValues, centerBer, lowerBer, upperBer, maxOrder); +end + +function [xBand, centerBand, yBounds] = directBerBounds(sirValues, centerBer, lowerBer, upperBer) +valid = isfinite(sirValues) & isfinite(centerBer) & isfinite(lowerBer) & ... + isfinite(upperBer) & centerBer > 0 & lowerBer > 0 & upperBer > 0; + +xBand = sirValues(valid).'; +centerBand = centerBer(valid).'; +lowerBand = lowerBer(valid).'; +upperBand = upperBer(valid).'; + +[xBand, orderIdx] = sort(xBand(:)); +centerBand = centerBand(orderIdx); +lowerBand = lowerBand(orderIdx); +upperBand = upperBand(orderIdx); + +lowerTmp = min(lowerBand, upperBand); +upperBand = max(lowerBand, upperBand); +lowerBand = lowerTmp; + +centerBand = min(max(centerBand, lowerBand), upperBand); +yBounds = [max(centerBand - lowerBand, 0), max(upperBand - centerBand, 0)]; +end + +function [xBand, centerBand, yBounds] = fittedBerBounds(sirValues, meanBer, minBer, maxBer, maxOrder) +valid = isfinite(sirValues) & isfinite(meanBer) & isfinite(minBer) & ... + isfinite(maxBer) & meanBer > 0 & minBer > 0 & maxBer > 0; + +x = sirValues(valid); +yMean = meanBer(valid); +yMin = minBer(valid); +yMax = maxBer(valid); + +if numel(x) < 2 + xBand = x(:); + centerBand = yMean(:); + yBounds = [max(yMean(:) - yMin(:), 0), max(yMax(:) - yMean(:), 0)]; + return +end + +[x, orderIdx] = sort(x(:)); +yMean = yMean(orderIdx); +yMin = yMin(orderIdx); +yMax = yMax(orderIdx); + +xBand = linspace(min(x), max(x), 300).'; +fitOrder = min(maxOrder, numel(unique(x)) - 1); + +if fitOrder < 1 + centerBand = interp1(x, yMean, xBand, "linear", "extrap"); + lowerBand = interp1(x, yMin, xBand, "linear", "extrap"); + upperBand = interp1(x, yMax, xBand, "linear", "extrap"); +else + centerBand = fitLogBer(x, yMean, xBand, fitOrder); + lowerBand = fitLogBer(x, yMin, xBand, fitOrder); + upperBand = fitLogBer(x, yMax, xBand, fitOrder); +end + +lowerTmp = min(lowerBand, upperBand); +upperBand = max(lowerBand, upperBand); +lowerBand = lowerTmp; + +centerBand = min(max(centerBand, lowerBand), upperBand); +yBounds = [max(centerBand - lowerBand, 0), max(upperBand - centerBand, 0)]; +end + +function yFit = fitLogBer(x, y, xFit, fitOrder) +coeff = polyfit(x, log10(y), fitOrder); +yFit = 10 .^ polyval(coeff, xFit); +end + +function [xFit, yFit, requiredSir, fitOrder] = fitBerAtFec( ... + sirValues, meanBer, polyfitOrderMax, fecBerThreshold, crossingSirWindow) + +xFit = NaN; +yFit = NaN; +requiredSir = NaN; +fitOrder = NaN; + +valid = isfinite(sirValues) & isfinite(meanBer) & meanBer > 0 & ... + sirValues >= crossingSirWindow(1) & sirValues <= crossingSirWindow(2); +if nnz(valid) < 2 + return +end + +sirValues = sirValues(valid); +meanBer = meanBer(valid); +[sirValues, sortIdx] = sort(sirValues); +meanBer = meanBer(sortIdx); + +fitOrder = min(polyfitOrderMax, nnz(valid) - 1); +fitCoeff = polyfit(sirValues, log10(meanBer), fitOrder); +xFit = linspace(max(min(sirValues), crossingSirWindow(1)), ... + min(max(sirValues), crossingSirWindow(2)), 300); +yFit = 10 .^ polyval(fitCoeff, xFit); + +thresholdMask = isfinite(yFit) & yFit <= fecBerThreshold; +if ~any(thresholdMask) + return +end + +firstThresholdIdx = find(thresholdMask, 1, "first"); +if firstThresholdIdx == 1 + requiredSir = xFit(firstThresholdIdx); + return +end + +xPair = xFit(firstThresholdIdx - 1:firstThresholdIdx); +yPair = log10(yFit(firstThresholdIdx - 1:firstThresholdIdx)); +if all(isfinite(yPair)) && diff(yPair) ~= 0 + requiredSir = interp1(yPair, xPair, log10(fecBerThreshold), ... + "linear", "extrap"); +else + requiredSir = xFit(firstThresholdIdx); +end + +if requiredSir < crossingSirWindow(1) || requiredSir > crossingSirWindow(2) + requiredSir = NaN; +end +end + +function label = algorithmDisplayName(algorithmName) +algorithmName = string(algorithmName); +switch algorithmName + case {"plain_ffe", "conventional_ffe"} + label = "FFE only"; + case "a2_tracked_levels" + label = "ACT"; + case "a2_residual" + label = "L-DCA"; + case "a1_moving_average" + label = "DCA"; + case "dc_tracking" + label = "DCT"; + otherwise + label = algorithmName; +end +end + +function color = algorithmColor(algorithmName) +algorithmName = string(algorithmName); +switch algorithmName + case {"plain_ffe", "conventional_ffe"} + color = [0.3467 0.5360 0.6907]; + case "a2_tracked_levels" + color = [0.9153 0.2816 0.2878]; + case "a2_residual" + color = [0.4416 0.7490 0.4322]; + case "a1_moving_average" + color = [1.0000 0.5984 0.2000]; + case "dc_tracking" + color = [0.6769 0.4447 0.7114]; + otherwise + color = [0 0 0]; +end +end + +function color = linewidthColor(linewidthIdx, nLinewidths) +if nLinewidths <= 1 + color = [0.3467 0.5360 0.6907]; + return +end + +if exist("cbrewer2", "file") + colors = cbrewer2("Set1", max(nLinewidths, 3)); +else + colors = lines(nLinewidths); +end + +color = colors(linewidthIdx, :); +end + +function label = linewidthLabel(laserLinewidth) +if abs(laserLinewidth) >= 1e6 + label = sprintf("%.3g MHz", laserLinewidth * 1e-6); +elseif abs(laserLinewidth) >= 1e3 + label = sprintf("%.3g kHz", laserLinewidth * 1e-3); +else + label = sprintf("%.3g Hz", laserLinewidth); +end +end + +function marker = algorithmMarker(algorithmName, algorithmMarkers) +algorithmOrder = ["conventional_ffe", "dc_tracking", "a2_tracked_levels", ... + "a2_residual", "a1_moving_average"]; +markerIdx = find(algorithmOrder == string(algorithmName), 1); +if isempty(markerIdx) + markerIdx = 1; +end +marker = algorithmMarkers{mod(markerIdx - 1, numel(algorithmMarkers)) + 1}; +end diff --git a/projects/Diss/MPI_revisit/simulation/mpi_simulation_worker.m b/projects/Diss/MPI_revisit/simulation/mpi_simulation_worker.m new file mode 100644 index 0000000..11a5051 --- /dev/null +++ b/projects/Diss/MPI_revisit/simulation/mpi_simulation_worker.m @@ -0,0 +1,214 @@ +function output = mpi_simulation_worker(userParameters, simulation_config) +%MPI_SIMULATION_WORKER Build one decorrelated MPI point and run the DSP recipe. + +arguments + userParameters struct + simulation_config struct +end + +config = applyDefaults(simulation_config); + +[Scpe_sig_raw, Symbols, Tx_bits, dataTable] = buildMpiScopeSignal(userParameters, config); + +output = config.recipe(Scpe_sig_raw, Symbols, Tx_bits, ... + "fsym", config.fsym, ... + "M", config.M, ... + "duob_mode", config.duob_mode, ... + "dataTable", dataTable, ... + "userParameters", userParameters, ... + "debug_plots", config.debug_plots); +end + +function config = applyDefaults(config) +defaults = struct( ... + "M", 4, ... + "fsym", 112e9, ... + "mpi_path_meter", 1000, ... + "laser_linewidth", 150e3, ... + "recipe", @mpi_recipe_dev, ... + "debug_plots", false, ... + "fdac", 256e9, ... + "fadc", 256e9, ... + "kover", 16, ... + "random_key", 1, ... + "rcalpha", 0.05, ... + "duob_mode", db_mode.no_db, ... + "vbias_rel", 0.5, ... + "u_pi", 3, ... + "laser_wavelength", 1293, ... + "link_length_km", 1, ... + "rop", -9, ... + "rx_bwl", 80e9, ... + "scope_bwl", 110e9, ... + "alpha", 0); + +names = fieldnames(defaults); +for nameIdx = 1:numel(names) + name = names{nameIdx}; + if ~isfield(config, name) || isempty(config.(name)) + config.(name) = defaults.(name); + end +end +end + +function [Scpe_sig, Symbols, Tx_bits, dataTable] = buildMpiScopeSignal(userParameters, config) +sir = readUserParameter(userParameters, "sir", 30); +randomKey = readUserParameter(userParameters, "random_key", config.random_key); +laserLinewidth = readUserParameter(userParameters, "laser_linewidth", config.laser_linewidth); +config.laser_linewidth = laserLinewidth; + +Pform = Pulseformer( ... + "fsym", config.fsym, ... + "fdac", 4 * config.fsym, ... + "pulse", "rrc", ... + "pulselength", 16, ... + "alpha", config.rcalpha); + +vbias = -config.vbias_rel * config.u_pi; +mainSource = buildOpticalSource(Pform, config, randomKey, randomKey + 1, vbias); +interferenceSource = buildOpticalSource(Pform, config, randomKey + 100000, ... + randomKey + 100001, vbias); + +Symbols = mainSource.Symbols; +Tx_bits = mainSource.Tx_bits; +Opt_sig = applyDecorrelatedMpi(mainSource.Opt_sig, interferenceSource.Opt_sig, ... + sir, config.link_length_km); + +Rx_sig = Amplifier( ... + "amp_mode", "ideal_no_noise", ... + "gain_mode", "output_power", ... + "amplification_db", config.rop).process(Opt_sig); + +Rx_sig = Photodiode( ... + "fsimu", config.fdac * config.kover, ... + "dark_current", 2e-08, ... + "responsivity", 1, ... + "temperature", 20, ... + "nep", 1.8e-11).process(Rx_sig); +Rx_sig.signal = real(Rx_sig.signal); + +Rx_sig = Filter( ... + "filtdegree", 4, ... + "f_cutoff", config.rx_bwl, ... + "fs", config.fdac * config.kover, ... + "filterType", filtertypes.butterworth, ... + "active", true).process(Rx_sig); +Rx_sig.signal = real(Rx_sig.signal); + +scopeFilter = Filter( ... + "filtdegree", 4, ... + "f_cutoff", config.scope_bwl, ... + "fs", config.fadc, ... + "filterType", filtertypes.butterworth, ... + "active", true); + +Scpe_sig = Scope( ... + "fsimu", config.fdac * config.kover, ... + "fadc", config.fadc, ... + "delay", 0, ... + "fixed_delay", 0, ... + "filtertype", filtertypes.butterworth, ... + "samplingdelay", 0, ... + "rand_samplingdelay", 0, ... + "freq_offset", 0, ... + "samp_jitter", 0, ... + "adcresolution", 8, ... + "quantbuffer", 0.1, ... + "block_dc", 1, ... + "lpf_active", 1, ... + "H_lpf", scopeFilter).process(Rx_sig); +Scpe_sig.signal = real(Scpe_sig.signal); + +HighpassFilter = Filter( ... + "filtdegree", 6, ... + "f_cutoff", 1e6, ... + "fs", config.fadc, ... + "filterType", filtertypes.butterworth, ... + "active", true, "lowpass",0); + +Scpe_sig = HighpassFilter.process(Scpe_sig); + +dataTable = table( ... + sir, ... + config.fsym, ... + config.M, ... + config.mpi_path_meter, ... + laserLinewidth, ... + randomKey, ... + 'VariableNames', {'sir', 'symbolrate', 'pam_level', ... + 'interference_path_length', 'laser_linewidth', 'random_key'}); +end + +function source = buildOpticalSource(Pform, config, sourceRandomKey, laserRandomKey, vbias) +[Digi_sig, Symbols, Tx_bits] = PAMsource( ... + "fsym", config.fsym, ... + "M", config.M, ... + "order", 18, ... + "useprbs", 0, ... + "fs_out", config.fdac, ... + "applyclipping", 0, ... + "clipfactor", 1.5, ... + "applypulseform", 1, ... + "pulseformer", Pform, ... + "randkey", sourceRandomKey, ... + "duobinary_mode", config.duob_mode, ... + "mrds_code", 0, ... + "mrds_blocklength", 512).process(); + +El_sig = M8199A("kover", config.kover).process(Digi_sig); + + +El_sig = Filter("f_cutoff",65e9,"filterType","butterworth","filtdegree",4,"fs",El_sig.fs).process(El_sig); + +El_sig = El_sig.normalize("mode", "oneone"); + +scaling = 0.6 * (config.u_pi / 2 - abs(vbias - config.u_pi / 2)); +El_sig = El_sig .* scaling; + +Opt_sig = EML( ... + "mode", eml_mode.im_cosinus, ... + "power", 3, ... + "fsimu", El_sig.fs, ... + "lambda", config.laser_wavelength, ... + "bias", vbias, ... + "u_pi", config.u_pi, ... + "linewidth", config.laser_linewidth, ... + "randomkey", laserRandomKey, ... + "alpha", config.alpha).process(El_sig); + +source = struct( ... + "Opt_sig", Opt_sig, ... + "Symbols", Symbols, ... + "Tx_bits", Tx_bits); +end + +function Opt_sig = applyDecorrelatedMpi(main_sig, interference_sig, sir, linkLengthKm) +interference_sig = Amplifier( ... + "amp_mode", "ideal_no_noise", ... + "gain_mode", "output_power", ... + "amplification_db", main_sig.power - sir).process(interference_sig); + +if numel(main_sig.signal) ~= numel(interference_sig.signal) + minLength = min(numel(main_sig.signal), numel(interference_sig.signal)); + main_sig.signal = main_sig.signal(1:minLength); + interference_sig.signal = interference_sig.signal(1:minLength); +end +combined_sig = main_sig + interference_sig; + +Opt_sig = Fiber( ... + "fsimu", combined_sig.fs, ... + "fiber_length", linkLengthKm, ... + "alpha", 0.3, ... + "D", 0, ... + "lambda0", 1310, ... + "gamma", 0, ... + "Dslope", 0.07).process(combined_sig); +end + +function value = readUserParameter(userParameters, name, defaultValue) +if isfield(userParameters, name) + value = userParameters.(name); +else + value = defaultValue; +end +end diff --git a/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_091910.mat b/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_091910.mat new file mode 100644 index 0000000..06ceb50 Binary files /dev/null and b/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_091910.mat differ diff --git a/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_092424.mat b/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_092424.mat new file mode 100644 index 0000000..74f5143 Binary files /dev/null and b/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_092424.mat differ diff --git a/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_095244.mat b/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_095244.mat new file mode 100644 index 0000000..b4deaa5 Binary files /dev/null and b/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_095244.mat differ diff --git a/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_100202.mat b/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_100202.mat new file mode 100644 index 0000000..43847eb Binary files /dev/null and b/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_100202.mat differ diff --git a/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_101144.mat b/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_101144.mat new file mode 100644 index 0000000..b26653a Binary files /dev/null and b/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_101144.mat differ diff --git a/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_104628.mat b/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_104628.mat new file mode 100644 index 0000000..9477d83 Binary files /dev/null and b/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_104628.mat differ diff --git a/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_104748.mat b/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_104748.mat new file mode 100644 index 0000000..771e348 Binary files /dev/null and b/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_104748.mat differ diff --git a/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_104851.mat b/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_104851.mat new file mode 100644 index 0000000..8b8e83f Binary files /dev/null and b/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_104851.mat differ diff --git a/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_105718.mat b/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_105718.mat new file mode 100644 index 0000000..439c00f Binary files /dev/null and b/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_105718.mat differ diff --git a/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_105936.mat b/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_105936.mat new file mode 100644 index 0000000..bfb3aca Binary files /dev/null and b/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_105936.mat differ diff --git a/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_110122.mat b/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_110122.mat new file mode 100644 index 0000000..07ce9c0 Binary files /dev/null and b/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_110122.mat differ diff --git a/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_110410.mat b/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_110410.mat new file mode 100644 index 0000000..20ab49f Binary files /dev/null and b/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_110410.mat differ diff --git a/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_110802.mat b/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_110802.mat new file mode 100644 index 0000000..300193d Binary files /dev/null and b/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_110802.mat differ diff --git a/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_111645.mat b/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_111645.mat new file mode 100644 index 0000000..8600f11 Binary files /dev/null and b/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_111645.mat differ diff --git a/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_122431.mat b/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_122431.mat new file mode 100644 index 0000000..d2afda5 Binary files /dev/null and b/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_122431.mat differ diff --git a/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_130744.mat b/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_130744.mat new file mode 100644 index 0000000..370ea2d Binary files /dev/null and b/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_130744.mat differ diff --git a/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_134103.mat b/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_134103.mat new file mode 100644 index 0000000..79e0500 Binary files /dev/null and b/projects/Diss/MPI_revisit/simulation/results/mpi_simulation_20260717_134103.mat differ diff --git a/projects/Diss/MPI_revisit/simulation/run_mpi_simulation_recipe.m b/projects/Diss/MPI_revisit/simulation/run_mpi_simulation_recipe.m new file mode 100644 index 0000000..7a92282 --- /dev/null +++ b/projects/Diss/MPI_revisit/simulation/run_mpi_simulation_recipe.m @@ -0,0 +1,476 @@ +%RUN_MPI_SIMULATION_RECIPE Queue MPI simulation points and store DSP packages. + +clear; +% clc; + +%% Fixed experiment-like simulation setup +simulation_config = struct(); +simulation_config.M = 4; +simulation_config.fsym = 112e9; +simulation_config.mpi_path_meter = 1000; +simulation_config.laser_linewidth = 150e3;%[150e3 150e3 250e3 500e3 750e3 1e6 5e6 10e6 20e6 50e6]; +simulation_config.laser_linewidths = simulation_config.laser_linewidth; +simulation_config.recipe = @mpi_recipe_dev; +simulation_config.debug_plots = true; +simulation_config.waitbar = true; +simulation_config.processing_mode = processingMode.parallel; +simulation_config.num_workers = 0; +simulation_config.random_keys = 1:10; + +simulation_config.fdac = 256e9; +simulation_config.fadc = 256e9; +simulation_config.kover = 16; +simulation_config.random_key = 1; +simulation_config.rcalpha = 0.05; +simulation_config.duob_mode = db_mode.no_db; +simulation_config.vbias_rel = 0.5; +simulation_config.u_pi = 3; +simulation_config.laser_wavelength = 1310; +simulation_config.link_length_km = 0.5; +simulation_config.rop = -7.5; +simulation_config.rx_bwl = 70e9; +simulation_config.scope_bwl = 70e9; +simulation_config.alpha = 0; + +%% Warehouse sweep setup +sweep_params = struct(); +sweep_params.sir = 15:3:45; +sweep_params.block_update = 1; +sweep_params.random_key = simulation_config.random_keys; +sweep_params.laser_linewidth = simulation_config.laser_linewidths; + +wh = DataStorage(sweep_params); + +fprintf("Requested %d MPI simulation job(s).\n", wh.getLastLinIndice()); + +%% Run queued simulations +[results, wh] = submitMpiSimulationJobs(wh, simulation_config, ... + "mode", simulation_config.processing_mode, ... + "waitbar", simulation_config.waitbar, ... + "numWorkers", simulation_config.num_workers); + +%% Save result artifact +result_dir = fullfile(fileparts(mfilename("fullpath")), "results"); +if ~exist(result_dir, "dir") + mkdir(result_dir); +end + +timestamp = string(datetime("now", "Format", "yyyyMMdd_HHmmss")); +result_file = fullfile(result_dir, "mpi_simulation_" + timestamp + ".mat"); +save(result_file, "wh", "simulation_config", "results"); +fprintf("Saved MPI simulation warehouse to:\n%s\n", result_file); + +%% BER inspection +plotMpiSimulationBer(wh); + +function plotMpiSimulationBer(wh) +storageNames = string(fieldnames(wh.sto)); +if isempty(storageNames) + warning("run_mpi_simulation_recipe:NoStorage", ... + "The warehouse does not contain any stored DSP packages."); + return +end + +useBoundedLines = true; +usePolyfit = true; +polyfitOrderMax = 4; +boundaryPolyfitOrderMax = 4; +fecBerThreshold = 3.8e-3; +maxBerForPlot = 0.1; +boundMode = "fitStd"; +algorithmMarkers = {'o','square','diamond','^','v','>','<','pentagram'}; + +cleanData = collectMpiSimulationBerRows(wh); +cleanData = normalizeBerColumnName(cleanData); +cleanData = cleanData(isfinite(cleanData.BER) & cleanData.BER > 0 & ... + cleanData.BER < maxBerForPlot, :); + +if isempty(cleanData) + warning("run_mpi_simulation_recipe:NoBerRows", ... + "No valid BER rows remain for plotting."); + return +end + +cleanData.clean_keep = true(height(cleanData), 1); +groupId = findgroups(cleanData.storage_name, cleanData.block_update, cleanData.sir); +for curGroup = unique(groupId(isfinite(groupId))).' + rowMask = groupId == curGroup; + berValues = cleanData.BER(rowMask); + if nnz(rowMask) > 3 + cleanData.clean_keep(rowMask) = ~isoutlier(berValues); + end +end +cleanData = cleanData(cleanData.clean_keep, :); + +groupVars = ["storage_name", "algorithm", "block_update", "sir"]; +summaryTable = groupsummary(cleanData, groupVars, {"mean", "min", "max"}, "BER"); +summaryTable = sortrows(summaryTable, groupVars); +summaryTable.sir_exact = summaryTable.sir; +summaryTable = addBerStdBounds(summaryTable, cleanData, groupVars); + +selectedBlockUpdates = unique(cleanData.block_update(isfinite(cleanData.block_update))).'; +selectedAlgorithms = unique(cleanData.storage_name, "stable").'; + +fprintf("Cleaned to %d simulation BER rows across %d storage/block/SIR groups.\n", ... + height(cleanData), height(summaryTable)); +disp(groupcounts(cleanData, ["storage_name", "block_update"])); + +figure(); +clf; +tiledlayout(numel(selectedBlockUpdates), 1, "TileSpacing", "compact"); + +for blockIdx = 1:numel(selectedBlockUpdates) + blockUpdate = selectedBlockUpdates(blockIdx); + nexttile; hold on; + + for algIdx = 1:numel(selectedAlgorithms) + storageName = selectedAlgorithms(algIdx); + algColor = algorithmColor(storageName); + marker = algorithmMarker(storageName, algorithmMarkers); + displayName = algorithmDisplayName(storageName); + + rawMask = cleanData.storage_name == storageName & ... + cleanData.block_update == blockUpdate; + curveMask = summaryTable.storage_name == storageName & ... + summaryTable.block_update == blockUpdate; + + if ~any(curveMask) + continue + end + + scatterRows = cleanData(rawMask, :); + scatter(scatterRows.sir, scatterRows.BER, ... + 26, ... + "Marker", ".", ... + "MarkerEdgeColor", algColor, ... + "MarkerFaceColor", algColor, ... + "HandleVisibility", "off"); + + sirValues = summaryTable.sir_exact(curveMask).'; + meanBer = summaryTable.mean_BER(curveMask).'; + boundCenterBer = summaryTable.std_center_BER(curveMask).'; + boundLowerBer = summaryTable.std_lower_BER(curveMask).'; + boundUpperBer = summaryTable.std_upper_BER(curveMask).'; + valid = isfinite(sirValues) & isfinite(meanBer) & meanBer > 0; + + if useBoundedLines && exist("boundedline", "file") && any(valid) + [xBand, centerBand, yBounds] = berStdBounds( ... + sirValues, boundCenterBer, boundLowerBer, boundUpperBer, ... + boundMode, boundaryPolyfitOrderMax); + + [hl, hp] = boundedline(xBand, centerBand, yBounds, ... + 'alpha', 'transparency', 0.1, ... + 'cmap', algColor, ... + 'nan', 'fill', ... + 'orientation', 'vert'); + set(hl, "LineStyle", "none", "LineWidth", 1, "Marker", "none", ... + "HandleVisibility", "off", "DisplayName", char(displayName)); + set(hp, "LineStyle", "-", "HandleVisibility", "off", "Marker", "none"); + end + + plot(sirValues(valid), meanBer(valid), ... + "LineStyle", "-", ... + "Marker", marker, ... + "MarkerSize", 3, ... + "LineWidth", 1, ... + "Color", algColor, ... + "MarkerFaceColor", "w", ... + "MarkerEdgeColor", algColor, ... + "DisplayName", char(displayName), ... + "HandleVisibility", "on"); + + 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(fecBerThreshold, "LineWidth", 1, "LineStyle", "--", "HandleVisibility", "off"); + yline(2e-2, "LineWidth", 1, "LineStyle", "--", "HandleVisibility", "off"); + + title(sprintf("Simulated MPI, block update %.0f", blockUpdate)); + xlabel("SIR (dB)"); + ylabel("BER"); + set(gca, "YScale", "log"); + ylim([9e-5, maxBerForPlot]); + xlim([min(cleanData.sir) - 1, max(cleanData.sir) + 1]); + grid on; + box on; + legend("Location", "northeast", "Interpreter", "none"); + + if exist("beautifyBERplot", "file") + beautifyBERplot("logscale", true, "setcolors", false, "setmarkers", false); + end +end +end + +function data = collectMpiSimulationBerRows(wh) +storageNames = string(fieldnames(wh.sto)); +data = table(); + +for storageIdx = 1:numel(storageNames) + storageName = storageNames(storageIdx); + storage = wh.sto.(char(storageName)); + + for linIdx = 1:numel(storage) + package = storage{linIdx}; + ber = extractPackageBer(package); + if ~isfinite(ber) + continue + end + + [physValues, physNames] = wh.getPhysIndicesByLinIndex(linIdx); + phys = struct(); + for physIdx = 1:numel(physNames) + phys.(char(physNames{physIdx})) = physValues{physIdx}; + end + + algorithm = extractPackageAlgorithm(package, storageName); + newRow = table( ... + storageName, ... + algorithm, ... + readPhysValue(phys, "sir", NaN), ... + readPhysValue(phys, "block_update", NaN), ... + readPhysValue(phys, "random_key", NaN), ... + readPhysValue(phys, "laser_linewidth", NaN), ... + ber, ... + 'VariableNames', {'storage_name', 'algorithm', 'sir', ... + 'block_update', 'random_key', 'laser_linewidth', 'BER'}); + data = [data; newRow]; %#ok + end +end +end + +function data = normalizeBerColumnName(data) +variableNames = string(data.Properties.VariableNames); +if ismember("BER", variableNames) + return +end + +if ismember("ber", variableNames) + data.Properties.VariableNames(variableNames == "ber") = {'BER'}; + return +end + +error("run_mpi_simulation_recipe:MissingBerColumn", ... + "Could not find a BER or ber column in the simulation BER table."); +end + +function value = readPhysValue(phys, name, defaultValue) +if isfield(phys, name) + value = phys.(name); +else + value = defaultValue; +end +end + +function ber = extractPackageBer(package) +ber = NaN; +if isempty(package) + return +end + +if iscell(package) + package = package{1}; +end + +if isstruct(package) && isfield(package, "metrics") + metrics = package.metrics; + if isstruct(metrics) && isfield(metrics, "BER") + ber = metrics.BER; + elseif isobject(metrics) && isprop(metrics, "BER") + ber = metrics.BER; + end +end +end + +function algorithm = extractPackageAlgorithm(package, fallbackName) +algorithm = fallbackName; +if isempty(package) + return +end + +if iscell(package) + package = package{1}; +end + +if isstruct(package) && isfield(package, "mpi_reduction_config") && ... + isfield(package.mpi_reduction_config, "algorithm") + algorithm = string(package.mpi_reduction_config.algorithm); +end +end + +function summaryTable = addBerStdBounds(summaryTable, cleanData, groupVars) +nGroups = height(summaryTable); +stdCenter = NaN(nGroups, 1); +stdLower = NaN(nGroups, 1); +stdUpper = NaN(nGroups, 1); + +for groupIdx = 1:nGroups + rowMask = true(height(cleanData), 1); + for varIdx = 1:numel(groupVars) + varName = groupVars(varIdx); + rowMask = rowMask & cleanData.(char(varName)) == summaryTable.(char(varName))(groupIdx); + end + + [stdCenter(groupIdx), stdLower(groupIdx), stdUpper(groupIdx)] = ... + logBerMeanStdInterval(cleanData.BER(rowMask)); +end + +summaryTable.std_center_BER = stdCenter; +summaryTable.std_lower_BER = stdLower; +summaryTable.std_upper_BER = stdUpper; +end + +function [centerBer, lowerBer, upperBer] = logBerMeanStdInterval(berValues) +berValues = berValues(isfinite(berValues) & berValues > 0); +if isempty(berValues) + centerBer = NaN; + lowerBer = NaN; + upperBer = NaN; + return +end + +logBer = log10(berValues(:)); +centerLog = mean(logBer, "omitnan"); +stdLog = std(logBer, 0, "omitnan"); + +centerBer = 10 .^ centerLog; +lowerBer = 10 .^ (centerLog - stdLog); +upperBer = 10 .^ (centerLog + stdLog); +end + +function [xBand, centerBand, yBounds] = berStdBounds(sirValues, centerBer, lowerBer, upperBer, boundMode, maxOrder) +if boundMode == "directStd" + [xBand, centerBand, yBounds] = directBerBounds(sirValues, centerBer, lowerBer, upperBer); + return +end + +[xBand, centerBand, yBounds] = fittedBerBounds(sirValues, centerBer, lowerBer, upperBer, maxOrder); +end + +function [xBand, centerBand, yBounds] = directBerBounds(sirValues, centerBer, lowerBer, upperBer) +valid = isfinite(sirValues) & isfinite(centerBer) & isfinite(lowerBer) & ... + isfinite(upperBer) & centerBer > 0 & lowerBer > 0 & upperBer > 0; + +xBand = sirValues(valid).'; +centerBand = centerBer(valid).'; +lowerBand = lowerBer(valid).'; +upperBand = upperBer(valid).'; + +[xBand, orderIdx] = sort(xBand(:)); +centerBand = centerBand(orderIdx); +lowerBand = lowerBand(orderIdx); +upperBand = upperBand(orderIdx); + +lowerTmp = min(lowerBand, upperBand); +upperBand = max(lowerBand, upperBand); +lowerBand = lowerTmp; + +centerBand = min(max(centerBand, lowerBand), upperBand); +yBounds = [max(centerBand - lowerBand, 0), max(upperBand - centerBand, 0)]; +end + +function [xBand, centerBand, yBounds] = fittedBerBounds(sirValues, meanBer, minBer, maxBer, maxOrder) +valid = isfinite(sirValues) & isfinite(meanBer) & isfinite(minBer) & ... + isfinite(maxBer) & meanBer > 0 & minBer > 0 & maxBer > 0; + +x = sirValues(valid); +yMean = meanBer(valid); +yMin = minBer(valid); +yMax = maxBer(valid); + +if numel(x) < 2 + xBand = x(:); + centerBand = yMean(:); + yBounds = [max(yMean(:) - yMin(:), 0), max(yMax(:) - yMean(:), 0)]; + return +end + +[x, orderIdx] = sort(x(:)); +yMean = yMean(orderIdx); +yMin = yMin(orderIdx); +yMax = yMax(orderIdx); + +xBand = linspace(min(x), max(x), 300).'; +fitOrder = min(maxOrder, numel(unique(x)) - 1); + +if fitOrder < 1 + centerBand = interp1(x, yMean, xBand, "linear", "extrap"); + lowerBand = interp1(x, yMin, xBand, "linear", "extrap"); + upperBand = interp1(x, yMax, xBand, "linear", "extrap"); +else + centerBand = fitLogBer(x, yMean, xBand, fitOrder); + lowerBand = fitLogBer(x, yMin, xBand, fitOrder); + upperBand = fitLogBer(x, yMax, xBand, fitOrder); +end + +lowerTmp = min(lowerBand, upperBand); +upperBand = max(lowerBand, upperBand); +lowerBand = lowerTmp; + +centerBand = min(max(centerBand, lowerBand), upperBand); +yBounds = [max(centerBand - lowerBand, 0), max(upperBand - centerBand, 0)]; +end + +function yFit = fitLogBer(x, y, xFit, fitOrder) +coeff = polyfit(x, log10(y), fitOrder); +yFit = 10 .^ polyval(coeff, xFit); +end + +function label = algorithmDisplayName(algorithmName) +algorithmName = string(algorithmName); +switch algorithmName + case {"plain_ffe", "conventional_ffe"} + label = "FFE only"; + case "a2_tracked_levels" + label = "ACT"; + case "a2_residual" + label = "L-DCA"; + case "a1_moving_average" + label = "DCA"; + case "dc_tracking" + label = "DCT"; + otherwise + label = algorithmName; +end +end + +function color = algorithmColor(algorithmName) +algorithmName = string(algorithmName); +switch algorithmName + case {"plain_ffe", "conventional_ffe"} + color = [0.3467 0.5360 0.6907]; + case "a2_tracked_levels" + color = [0.9153 0.2816 0.2878]; + case "a2_residual" + color = [0.4416 0.7490 0.4322]; + case "a1_moving_average" + color = [1.0000 0.5984 0.2000]; + case "dc_tracking" + color = [0.6769 0.4447 0.7114]; + otherwise + color = [0 0 0]; +end +end + +function marker = algorithmMarker(algorithmName, algorithmMarkers) +algorithmOrder = ["conventional_ffe", "dc_tracking", "a2_tracked_levels", ... + "a2_residual", "a1_moving_average"]; +markerIdx = find(algorithmOrder == string(algorithmName), 1); +if isempty(markerIdx) + markerIdx = 1; +end +marker = algorithmMarkers{mod(markerIdx - 1, numel(algorithmMarkers)) + 1}; +end diff --git a/projects/Diss/MPI_revisit/simulation/submitMpiSimulationJobs.m b/projects/Diss/MPI_revisit/simulation/submitMpiSimulationJobs.m new file mode 100644 index 0000000..f6fafc0 --- /dev/null +++ b/projects/Diss/MPI_revisit/simulation/submitMpiSimulationJobs.m @@ -0,0 +1,87 @@ +function [results, wh] = submitMpiSimulationJobs(wh, simulation_config, options) +%SUBMITMPISIMULATIONJOBS Execute MPI simulation warehouse points via runBatch. + +arguments + wh DataStorage + simulation_config struct + options.mode = processingMode.serial + options.waitbar (1,1) logical = true + options.numWorkers (1,1) double {mustBeNonnegative, mustBeInteger} = 0 + options.idleTimeout (1,1) double {mustBePositive} = 300 + options.cancelExistingQueue (1,1) logical = true +end + +nJobs = wh.getLastLinIndice(); +jobs = repmat(struct("args", {{}}, "label", "", "meta", struct()), 1, nJobs); + +for linIdx = 1:nJobs + userParameters = buildUserParameters(wh, linIdx); + jobs(linIdx).args = {userParameters, simulation_config}; + jobs(linIdx).label = buildJobLabel(userParameters, linIdx); + jobs(linIdx).meta.lin_idx = linIdx; + jobs(linIdx).meta.userParameters = userParameters; +end + +results = runBatch(@mpi_simulation_worker, jobs, ... + "mode", options.mode, ... + "waitbar", options.waitbar, ... + "waitbarMessage", "Processing MPI simulations...", ... + "numWorkers", options.numWorkers, ... + "idleTimeout", options.idleTimeout, ... + "cancelExistingQueue", options.cancelExistingQueue, ... + "resultHandler", @storeResult, ... + "errorHandler", @handleError); + + function userParameters = buildUserParameters(storageWh, linIdx) + userParameters = struct(); + if isempty(storageWh.getDimension()) + return + end + + [values, names] = storageWh.getPhysIndicesByLinIndex(linIdx); + for paramIdx = 1:numel(names) + userParameters.(char(names{paramIdx})) = values{paramIdx}; + end + end + + function label = buildJobLabel(userParameters, linIdx) + label = sprintf("MPI sim job %d", linIdx); + if isfield(userParameters, "sir") + label = sprintf("%s, SIR %g dB", label, userParameters.sir); + end + if isfield(userParameters, "block_update") + label = sprintf("%s, block %g", label, userParameters.block_update); + end + end + + function storeResult(val, job, ~) + if isempty(val) || ~isstruct(val) + return + end + + storageNames = fieldnames(val); + for storageIdx = 1:numel(storageNames) + storageName = storageNames{storageIdx}; + if isempty(val.(storageName)) + continue + end + + ensureStorage(storageName); + wh.addValueToStorageByLinIdx(val.(storageName), storageName, job.meta.lin_idx); + end + end + + function ensureStorage(storageName) + if ~isfield(wh.sto, storageName) + wh.addStorage(storageName); + end + end + + function handleError(ME, job, ~) + fprintf("[%s] ERROR [%s]: %s\n", job.label, ME.identifier, ME.message); + for st = ME.stack' + fprintf(" %s:%d (%s)\n", st.file, st.line, st.name); + end + fprintf("Full report:\n%s\n", getReport(ME, "extended")); + end +end diff --git a/projects/IMDD_base_system/minimal_example.m b/projects/IMDD_base_system/minimal_example.m index 2b97947..b0891d5 100644 --- a/projects/IMDD_base_system/minimal_example.m +++ b/projects/IMDD_base_system/minimal_example.m @@ -133,14 +133,19 @@ Scpe_sig.signal = Scpe_sig.signal(1:2*length(Symbols)); if 1 % -------------------- FFE -------------------- ffe_order = [150, 0, 0]; - eq_ = EQ("Ne",ffe_order,"Nb",[2,0,0], ... + eq_ = EQ("Ne",ffe_order,"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); % mu_tr_rls = 9.805e-01; mu_dd_rls = 0.999989348903919; eq_ = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",1.933e-04,"order",ffe_order(1),"sps",2,"decide",0,"optmize_mus",1,"dd_mode",1,"adaption_technique","nlms"); - + eq_ = FFE_plain("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",1.933e-04,"order",ffe_order(1),"sps",2,"decide",0,"optmize_mus",1,"dd_mode",1,"adaption_technique","nlms"); + + eq_ = VNLE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2, ... + "mu_dd",[0.0004 0.0005 0.0006],"mu_tr",[0.0001 0.0008 0.001], ... + "order",[150,5,5],"sps",2,"decide",1, ... + "optmize_mus",1,"mu_optimization_len",2^15); % eq_ = FFE("epochs_tr",4,"epochs_dd",5,"len_tr",4096,"mu_dd",0.01,"mu_tr",0.01,"order",50,"sps",2,"decide",0, "adaption",adaption_method.nlms,"dd_mode",1); % eq_ = FFE_DFE("epochs_tr",5,"epochs_dd",5,"len_tr",512,"ffe_mu_dd",1e-4,"dfe_mu_dd",5e-4,"ffe_mu_tr",0,"dfe_mu_tr",0,"ffe_order",99,"dfe_order",99,"sps",2,"decide",0); @@ -149,37 +154,37 @@ if 1 "precode_mode",duob_mode,'showAnalysis',1,"postFFE",[], ... "eth_style_symbol_mapping",0); - output.ffe_results.metrics.print("description",'DFE'); + output.ffe_results.metrics.print("description",'VNLE'); end %% - -% -------------------- VNLE + MLSE -------------------- -pf_ncoeffs = 1; -ffe_order3 = [50, 5, 5]; -eq_v = EQ("Ne",ffe_order3,"Nb",dfe_order, ... - "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",1); -pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); - -mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); - -[output.vnle_results, output.mlse_results] = vnle_postfilter_mlse(eq_v, pf_, mlse_, M, Scpe_sig, Symbols, Tx_bits, ... - "precode_mode", duob_mode, 'showAnalysis', 1, "postFFE", [], "eth_style_symbol_mapping", 0); - -output.mlse_results.metrics.print("description",'MLSE'); - +if 0 + % -------------------- VNLE + MLSE -------------------- + pf_ncoeffs = 1; + ffe_order3 = [50, 5, 5]; + eq_v = EQ("Ne",ffe_order3,"Nb",dfe_order, ... + "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",1); + pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); + + mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); + + [output.vnle_results, output.mlse_results] = vnle_postfilter_mlse(eq_v, pf_, mlse_, M, Scpe_sig, Symbols, Tx_bits, ... + "precode_mode", duob_mode, 'showAnalysis', 1, "postFFE", [], "eth_style_symbol_mapping", 0); + + output.mlse_results.metrics.print("description",'MLSE'); +end %% - -% -------------------- DB target -------------------- -mlse_db_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",M,'trellis_states',PAMmapper(M,0).levels); -ffe_order = [50, 5, 5]; -eq_ = EQ("Ne",ffe_order,"Nb",dfe_order,"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",1); -output.dbt_results = duobinary_target(eq_,mlse_db_, M, Scpe_sig, Symbols, Tx_bits, ... - "precode_mode", duob_mode, 'showAnalysis', 0, "postFFE", []); - -output.dbt_results.metrics.print("description",'Duobinary'); - +if 0 + % -------------------- DB target -------------------- + mlse_db_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",M,'trellis_states',PAMmapper(M,0).levels); + ffe_order = [50, 5, 5]; + eq_ = EQ("Ne",ffe_order,"Nb",dfe_order,"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",1); + output.dbt_results = duobinary_target(eq_,mlse_db_, M, Scpe_sig, Symbols, Tx_bits, ... + "precode_mode", duob_mode, 'showAnalysis', 0, "postFFE", []); + + output.dbt_results.metrics.print("description",'Duobinary'); +end diff --git a/workerError.mat b/workerError.mat index 34c4a9f..f03b730 100644 Binary files a/workerError.mat and b/workerError.mat differ