From 7085ba093133c39566c335f7c78f773a1b1efd5c Mon Sep 17 00:00:00 2001 From: Silas Oettinghaus Date: Fri, 24 Oct 2025 16:58:06 +0200 Subject: [PATCH] Gif stuff nonlinear MLSE investigation trying hard to implment ML-based pre Equalization to find the branch metrics --- Classes/00_signals/Signal.m | 52 ++-- Classes/04_DSP/Equalizer/FFE.m | 2 +- Classes/04_DSP/Equalizer/FFE_MLSE.m | 233 ++++++++++++++ Classes/04_DSP/Equalizer/ML_MLSE.m | 246 +++++++++++++++ Classes/04_DSP/Sequence Detection/MLSE.m | 114 ++++++- Classes/GifWriter.m | 146 +++++++++ Functions/EQ_structures/dsp_runid.m | 6 +- Functions/EQ_visuals/showEQNoisePSD.m | 2 +- Functions/Theory/CCDM/ccdm_gpt_example.m | 27 ++ Functions/beautifyBERplot.m | 13 +- .../Auswertung_JLT/run_dsp_from_db.m | 16 +- .../Auswertung_JLT/run_plot_measurements.m | 2 +- projects/IMDD_base_system/simulation_bwl.m | 286 ++++++++++++------ projects/IMDD_base_system/simulation_bwl_2.m | 8 +- projects/ML_based_MLSE/model.m | 194 ++++++++++++ .../Nonlinear_MLSE/simulation_nonlin_dsp.m | 229 ++++++++++++++ 16 files changed, 1414 insertions(+), 162 deletions(-) create mode 100644 Classes/04_DSP/Equalizer/FFE_MLSE.m create mode 100644 Classes/04_DSP/Equalizer/ML_MLSE.m create mode 100644 Classes/GifWriter.m create mode 100644 Functions/Theory/CCDM/ccdm_gpt_example.m create mode 100644 projects/ML_based_MLSE/model.m create mode 100644 projects/Nonlinear_MLSE/simulation_nonlin_dsp.m diff --git a/Classes/00_signals/Signal.m b/Classes/00_signals/Signal.m index c84ebf6..159b4e3 100644 --- a/Classes/00_signals/Signal.m +++ b/Classes/00_signals/Signal.m @@ -344,7 +344,7 @@ classdef Signal arguments obj - options.fignum + options.fignum = 2025 options.displayname = ""; options.color = []; options.normalizeToNyquist = 0; @@ -449,21 +449,25 @@ classdef Signal ylabel(ylab); + % --- Y-Axis scaling (auto with margin) --- + y_min = min(p_dbm(:)); + y_max = max(p_dbm(:)); + + % Add 5% dynamic range margin on both sides + y_range = y_max - y_min; + if y_range == 0 + y_range = 10; % fallback if flat + end + + y_margin = 0.05 * y_range; + ylim([y_min - y_margin, y_max + y_margin]); + + % Set ticks automatically, avoid overpopulation try - ylim([max(min(floor(min(p_dbm))-3, ax.YLim(1)),-40), min(max(ceil(max(p_dbm))+3, ax.YLim(2)),10)]); - catch - ylim([floor(min(p_dbm,[],'all'))-3, ceil(max(p_dbm,[],'all'))+3]); + yticks(round(linspace(y_min, y_max, min(10, max(4, ceil(y_range/10)))))); end - - if options.normalizeTo0dB - ylim([floor(min(p_dbm,[],'all'))-3, ceil(max(p_dbm,[],'all'))+3]); - else - ylim([floor(min(p_dbm,[],'all'))-3, ceil(max(p_dbm,[],'all'))+3]); - end - - yticks(-200:10:10); grid on; - legend + end @@ -961,8 +965,8 @@ classdef Signal maxA = max(sig(100:end-100))*1.3; minA = min(sig(100:end-100))*1.3; - % maxA = 0.0015; - % minA = 0; + maxA = 0.0025; + minA = 0; difference= maxA-minA; @@ -1012,7 +1016,7 @@ classdef Signal % add information - if 1 + if 0 pwr_dbm = round(obj.power,3); pwr_lin = obj.power("unit",power_notation.W); @@ -1121,20 +1125,20 @@ classdef Signal end - yticks(linspace(0,histpoints,6)); - y_tickstring = sprintfc('%.2f', y_tickstring); - yticklabels(y_tickstring); - - xticks(linspace(0,histpoints_horizontal,6)) - x_tickstring = sprintfc('%.2f', linspace(0, 2/fsym, 8) .* 1e12); - xticklabels(x_tickstring); - grid off + end + yticks(linspace(0,histpoints,6)); + y_tickstring = sprintfc('%.2f', y_tickstring); + yticklabels(y_tickstring); + xticks(linspace(0,histpoints_horizontal,6)) + x_tickstring = sprintfc('%.2f', linspace(0, 2/fsym, 8) .* 1e12); + xticklabels(x_tickstring); + % end % disp('h'); diff --git a/Classes/04_DSP/Equalizer/FFE.m b/Classes/04_DSP/Equalizer/FFE.m index 5920e89..18442d6 100644 --- a/Classes/04_DSP/Equalizer/FFE.m +++ b/Classes/04_DSP/Equalizer/FFE.m @@ -131,7 +131,7 @@ classdef FFE < handle mask = ones(obj.order,1); maincursor_pos=ceil(length(obj.e)/2); always_ideal_decision = 0; - save_debug = 1; + save_debug = 0; grad =0; weight = 0; update = 0; diff --git a/Classes/04_DSP/Equalizer/FFE_MLSE.m b/Classes/04_DSP/Equalizer/FFE_MLSE.m new file mode 100644 index 0000000..e7a4bbd --- /dev/null +++ b/Classes/04_DSP/Equalizer/FFE_MLSE.m @@ -0,0 +1,233 @@ +classdef ML_MLSE < handle + % Implementation of plain and simple FFE. + % 1) Training mode (stable performance when you use NLMS) + % 2) Decision directed mode + + %LMS: mu in order of 0.0001 for acceptable convergence speed + %NLMS: mu in order of 0.01 for acceptable convergence speed + %RLS: mu is lambda -> 0.99 -> 1 (has a strong dependency on this! use a loop to find out best values) + + % FFE("epochs_tr",5,"epochs_dd",2,"len_tr",2^13,"mu_dd",mu_dd,"mu_tr",mu_tr,"order",25,"sps",2,"decide",0, "adaption",adaption_method(adaption),"dd_mode",use_dd_mode); + + properties + sps % usually 2 + order + e + e_tr + error + + len_tr + mu_tr + epochs_tr + + dd_mode % 1 or 0 to set DD-mode on or off + mu_dd %weight update in dd mode + epochs_dd + + constellation + + L %viterbi memory length + + alpha + DIR + DIR_flip + trellis_states + + traceback_depth + end + + methods + function obj = ML_MLSE(options) + arguments(Input) + + options.sps = 2; + options.order = 15; + + options.len_tr = 4096; + options.mu_tr = 0; + options.epochs_tr = 5; + + options.dd_mode = 1; + options.mu_dd = 1e-5; + options.epochs_dd = 5; + + 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 + + function [X,X_viterbi] = process(obj, X, D) + + % actual processing of the signal (steps 1. - 3.) + % 1 normalize RMS + X = X.normalize("mode","rms"); + + obj.constellation = unique(D.signal); + + if length(X)/length(D) ~= obj.sps + warning('Signal length does not fit to reference!'); + end + + % Training Mode + n = obj.len_tr; + training = 1; + obj.equalize(X.signal, D.signal,obj.mu_tr,obj.epochs_dd,n,training); + obj.e_tr = obj.e; + + % Decision Directed Mode + n = X.length; + training = 0; + [y,y_vit]=obj.equalize(X.signal, D.signal,obj.mu_dd,obj.epochs_dd,n,training); + + X_viterbi = X; + + X.signal = y; + X.fs = D.fs; %change sampling frequency of outgoing signal from fdac e.g. 2 sps to symbol spaced = fsym + lbdesc = [num2str(obj.order),' tap FFE']; + X = X.logbookentry(lbdesc); % append to logbook + + X_viterbi.signal = y_vit; + X_viterbi.fs = D.fs; %change sampling frequency of outgoing signal from fdac e.g. 2 sps to symbol spaced = fsym + lbdesc = [num2str(obj.order),'order FFE + PF + Viterbi']; + X_viterbi = X_viterbi.logbookentry(lbdesc); % append to logbook + + + end + + function [y,y_vit] = equalize(obj,x,d,mu,epochs,N,training) + % ============================================================== + % FFE + Whitening + Viterbi Equalizer (reference implementation) + % ============================================================== + + % --- Input padding and preallocation + x = [zeros(floor(obj.order/2),1); x; zeros(obj.order,1)]; + N_ = N / obj.sps; + y = zeros(N_,1); + y_white = zeros(N_,1); + + for epoch = 1:epochs + + % ============================================================== + % INITIALIZATION (only before final epoch and detection mode) + % ============================================================== + if epoch == epochs && ~training + + % --- Parameters + S = numel(unique(d)); + L = obj.L; + nStates = S^L; + nFeasible = S^(L-1)*S; + + % --- Trellis setup + obj.DIR = arburg(y-d, L); + obj.DIR_flip = flip(obj.DIR); + obj.trellis_states = reshape(unique(d),1,[]); + + pre_comb_mat = repmat(obj.trellis_states, L, 1); + pre_comb_cell = mat2cell(pre_comb_mat, ones(1,L), size(pre_comb_mat,2)); + combs = fliplr(combvec(pre_comb_cell{:}).'); + first_sym = combs(:,1); + last_sym = combs(:,end); + nStates = size(combs,1); + + noise_free_received = inf(nStates,nStates); + valid = false(nStates); + + for from = 1:nStates + for to = 1:nStates + if all(combs(to,2:end) == combs(from,1:end-1)) + noise_free_received(to,from) = ... + dot(combs(to,:), obj.DIR_flip(end:-1:2)) + last_sym(from)*obj.DIR_flip(1); + valid(to,from) = true; + end + end + end + + nf_vec = noise_free_received(valid); + [valid_to, valid_from] = find(valid); + from_per_to = arrayfun(@(to)find(valid(to,:)), 1:nStates, 'UniformOutput', false); + + % --- Noise stats + y_ideal = conv(d(:), obj.DIR(:), "same"); + sigma2 = mean(abs(y - y_ideal).^2); + inv2s2 = 1/(2*sigma2); + + % --- Vector initialization + bm_vec = zeros(1,nFeasible); + pm = zeros(nStates,1); + pm_next = zeros(nStates,1); + bm_fw = zeros(nStates,nStates,length(y)); + zi = zeros(max(numel(obj.DIR)-1,0),1); + end + + % ============================================================== + % RUNTIME LOOP (FFE update + Viterbi detection in last epoch) + % ============================================================== + symbol = 0; + for sample = 1:obj.sps:N + symbol = symbol + 1; + + % --- FFE output + U = x(obj.order+sample-1:-1:sample); + y(symbol,1) = obj.e.' * U; + + % --- Decision + if training + d_hat(symbol,1) = d(symbol); + else + [~,symbol_idx] = min(abs(y(symbol) - obj.constellation)); + d_hat(symbol,1) = obj.constellation(symbol_idx); + end + + % --- LMS weight update + err(symbol) = d_hat(symbol) - y(symbol); + obj.e = obj.e + mu * (err(symbol) * U); + + % --- Whitening + Viterbi (final epoch only) + if epoch == epochs && ~training + + [y_white(symbol), zi] = filter(obj.DIR,1,y(symbol), zi); + + if symbol == 1 + pm = -inf(nStates,nStates); + pm(:,1:nStates) = 0; + else + bm_vec = -(y_white(symbol) - nf_vec).^2 * inv2s2; + bm_mat = -inf(nStates,nStates); + bm_mat(valid) = bm_vec; + + pm_new = pm + bm_mat; + [pm_survive(:,symbol), pm_survivor_fw_idx(:,symbol)] = max(pm_new,[],2); + pm = repmat(pm_survive(:,symbol).', nStates,1); + end + + % --- Traceback + if mod(symbol,obj.traceback_depth) == 0 + [~,viterbi_path(symbol)] = max(pm_survive(:,symbol)); + for n = symbol:-1:symbol-obj.traceback_depth+2 + viterbi_path(n-1) = pm_survivor_fw_idx(viterbi_path(n),n); + end + end + end + end + + % --- Final reconstruction + if epoch == epochs && ~training + y_vit = first_sym(viterbi_path); + end + end + end + + end +end \ No newline at end of file diff --git a/Classes/04_DSP/Equalizer/ML_MLSE.m b/Classes/04_DSP/Equalizer/ML_MLSE.m new file mode 100644 index 0000000..db6c686 --- /dev/null +++ b/Classes/04_DSP/Equalizer/ML_MLSE.m @@ -0,0 +1,246 @@ +classdef ML_MLSE < handle + % Implementation of plain and simple FFE. + % 1) Training mode (stable performance when you use NLMS) + % 2) Decision directed mode + + %LMS: mu in order of 0.0001 for acceptable convergence speed + %NLMS: mu in order of 0.01 for acceptable convergence speed + %RLS: mu is lambda -> 0.99 -> 1 (has a strong dependency on this! use a loop to find out best values) + + % FFE("epochs_tr",5,"epochs_dd",2,"len_tr",2^13,"mu_dd",mu_dd,"mu_tr",mu_tr,"order",25,"sps",2,"decide",0, "adaption",adaption_method(adaption),"dd_mode",use_dd_mode); + + properties + sps % usually 2 + order + e + e_tr + error + + len_tr + mu_tr + epochs_tr + + dd_mode % 1 or 0 to set DD-mode on or off + mu_dd %weight update in dd mode + epochs_dd + + constellation + + L %viterbi memory length + + alpha + DIR + DIR_flip + trellis_states + + traceback_depth + end + + methods + function obj = ML_MLSE(options) + arguments(Input) + + options.sps = 2; + options.order = 15; + + options.len_tr = 4096; + options.mu_tr = 0; + options.epochs_tr = 5; + + options.dd_mode = 1; + options.mu_dd = 1e-5; + options.epochs_dd = 5; + + 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 + + function [X,X_viterbi] = process(obj, X, D) + + % actual processing of the signal (steps 1. - 3.) + % 1 normalize RMS + X = X.normalize("mode","rms"); + + obj.constellation = unique(D.signal); + + if length(X)/length(D) ~= obj.sps + warning('Signal length does not fit to reference!'); + end + + % Training Mode + n = obj.len_tr; + training = 1; + obj.equalize(X.signal, D.signal,obj.mu_tr,obj.epochs_dd,n,training); + obj.e_tr = obj.e; + + % Decision Directed Mode + n = X.length; + training = 0; + [y,y_vit]=obj.equalize(X.signal, D.signal,obj.mu_dd,obj.epochs_dd,n,training); + + X_viterbi = X; + + X.signal = y; + X.fs = D.fs; %change sampling frequency of outgoing signal from fdac e.g. 2 sps to symbol spaced = fsym + lbdesc = [num2str(obj.order),' tap FFE']; + X = X.logbookentry(lbdesc); % append to logbook + + X_viterbi.signal = y_vit; + X_viterbi.fs = D.fs; %change sampling frequency of outgoing signal from fdac e.g. 2 sps to symbol spaced = fsym + lbdesc = [num2str(obj.order),'order FFE + PF + Viterbi']; + X_viterbi = X_viterbi.logbookentry(lbdesc); % append to logbook + + + end + + function [y,y_vit] = equalize(obj,x,d,mu,epochs,N,training) + % ============================================================== + % FFE + Whitening + ML-Based Branch Metric Estimation + Viterbi + % ============================================================== + + % --- Input padding and preallocation + x = [zeros(floor(obj.order/2),1); x; zeros(obj.order,1)]; + N_ = N / obj.sps; + y = zeros(N_,1); + y_white = zeros(N_,1); + + for epoch = 1:epochs + + % ============================================================== + % INITIALIZATION (only before final epoch and detection mode) + % ============================================================== + if epoch == epochs && ~training + + % --- Parameters + S = numel(unique(d)); % alphabet size + L = obj.L; % MLSE memory + Nf = L; % filter length + Delta = ceil(L/2); % delay parameter + nStates = S^L; + nFeasible = S^(L-1)*S; + + % --- Trellis mapping + obj.DIR = arburg(y-d, L); + obj.DIR_flip = flip(obj.DIR); + obj.trellis_states = reshape(unique(d),1,[]); + pre_comb_mat = repmat(obj.trellis_states, L, 1); + pre_comb_cell = mat2cell(pre_comb_mat, ones(1,L), size(pre_comb_mat,2)); + combs = fliplr(combvec(pre_comb_cell{:}).'); + first_sym = combs(:,1); + last_sym = combs(:,end); + nStates = size(combs,1); + + % --- Valid transitions + valid = false(nStates); + for from = 1:nStates + for to = 1:nStates + if all(combs(to,2:end) == combs(from,1:end-1)) + valid(to,from) = true; + end + end + end + [valid_to, valid_from] = find(valid); + + % --- Noise estimation + y_ideal = conv(d(:), obj.DIR(:), "same"); + sigma2 = mean(abs(y - y_ideal).^2); + inv2s2 = 1/(2*sigma2); + + % --- Allocate vectors and weights + pm = zeros(nStates,1); + w = zeros(Nf,nFeasible); % filter weights per transition + b = zeros(1,nFeasible); % bias terms + v_hat = zeros(1,nFeasible); + v_tilde = zeros(1,nFeasible); + bm_vec = zeros(1,nFeasible); + zi = zeros(max(numel(obj.DIR)-1,0),1); + end + + % ============================================================== + % RUNTIME LOOP + % ============================================================== + symbol = 0; + for sample = 1:obj.sps:N + symbol = symbol + 1; + + % --- FFE output + U = x(obj.order+sample-1:-1:sample); + y(symbol,1) = obj.e.' * U; + + % --- Decision / FFE adaptation + if training + d_hat = d(symbol); + else + [~,idx] = min(abs(y(symbol) - obj.constellation)); + d_hat = obj.constellation(idx); + end + err = d_hat - y(symbol); + obj.e = obj.e + mu * (err * U); + + % --- Whitening + MLSE in last epoch + if epoch == epochs && ~training + + [y_white(symbol), zi] = filter(obj.DIR,1,y(symbol), zi); + k = symbol; + + % --- Build Δ-delayed observation window y_k + i1 = k - Nf + 1 + Delta; + i2 = k + Delta; + buf = y_white(max(1,i1):min(length(y_white),i2)); + padL = max(0,1 - i1); + padR = max(0,i2 - length(y_white)); + yk = [zeros(padL,1); buf(:); zeros(padR,1)]; % Nf×1 + + % --- Predict branch metrics for all feasible transitions + v_hat = (yk.' * w) + b; % [1×nFeasible] + v_hat = v_hat.'; % [nFeasible×1] + + % --- Extended path metrics + v_tilde = pm(valid_from) + v_hat; % [nFeasible×1] + + % --- Compute branch metrics (distance) + bm_vec = -(y_white(k) - v_hat).^2 * inv2s2; % 1×nFeasible + + % --- Survivor selection (vector aggregation) + pm_new_vec = pm(valid_from) + bm_vec.'; % nFeasible×1 + pm_next = -inf(nStates,1); + surv_idx = zeros(nStates,1); + for t = 1:nStates + mask = (valid_to==t); + [pm_next(t), arg] = max(pm_new_vec(mask)); + surv_idx(t) = valid_from(find(mask,1,'first')-1+arg); + end + pm = pm_next; + + % --- Traceback + if mod(symbol,obj.traceback_depth) == 0 + [~,viterbi_path(symbol)] = max(pm); + for n = symbol:-1:symbol-obj.traceback_depth+2 + viterbi_path(n-1) = surv_idx(viterbi_path(n)); + end + end + end + end + + % --- Output reconstructed path + if epoch == epochs && ~training + y_vit = first_sym(viterbi_path); + end + end + end + + + end +end \ No newline at end of file diff --git a/Classes/04_DSP/Sequence Detection/MLSE.m b/Classes/04_DSP/Sequence Detection/MLSE.m index 3a8e288..fabee6f 100644 --- a/Classes/04_DSP/Sequence Detection/MLSE.m +++ b/Classes/04_DSP/Sequence Detection/MLSE.m @@ -277,14 +277,107 @@ classdef MLSE < handle if debug - alpha_ = alpha - min(alpha) + eps; - figure();hold on; - n = 10; - scatter(1:n,obj.trellis_states(repmat([1:numel(obj.trellis_states)]',1,n)),abs(alpha_(:,end-n+1:end)),'Marker','o','LineWidth',1); - scatter(1:n,obj.trellis_states(viterbi_path(end-n+1:end)),500,'Marker','x','LineWidth',1,'MarkerEdgeColor','green'); - % scatter(1:n,data_ref(end-n+1:end),500,'Marker','x','LineWidth',1,'MarkerEdgeColor','red'); - yticks(obj.trellis_states); - ylim([min(obj.trellis_states)-1 max(obj.trellis_states)+1]); + % alpha_ = alpha - min(alpha) + eps; + % figure();hold on; + % n = 10; + % scatter(1:n,obj.trellis_states(repmat([1:numel(obj.trellis_states)]',1,n)),abs(alpha_(:,end-n+1:end)),'Marker','o','LineWidth',1); + % scatter(1:n,obj.trellis_states(viterbi_path(end-n+1:end)),500,'Marker','x','LineWidth',1,'MarkerEdgeColor','green'); + % % scatter(1:n,data_ref(end-n+1:end),500,'Marker','x','LineWidth',1,'MarkerEdgeColor','red'); + % yticks(obj.trellis_states); + % ylim([min(obj.trellis_states)-1 max(obj.trellis_states)+1]); + + %% ----- Build true path from known sequence (data_ref) ----- + % Memory length used by VA (L = length(obj.DIR)-1 symbols stored in state) + L = size(combs,2); % each row of combs is the L-tap state vector + N = length(data_in); + + % Map each trellis level to an index (1..M) + [~, level_to_idx] = ismember(levels, levels); %#ok % identity map + [ok_ref, ref_idx] = ismember(data_ref(:).', levels); + if ~all(ok_ref) + warning('Some data_ref symbols are not in "levels". True-path build may fail.'); + end + + % Precompute next_state(from_state, u_idx) LUT such that: + % combs(next_state,:) == [ combs(from_state,2:end) , levels(u_idx) ] + next_state = zeros(size(combs,1), numel(levels), 'uint32'); + for from = 1:size(combs,1) + prefix = combs(from,2:end); % what must match in 'to' for a valid transition + for ui = 1:numel(levels) + target = [prefix, levels(ui)]; + % find the unique 'to' whose state vector equals target + to = find(all(bsxfun(@eq, combs, target), 2), 1, 'first'); + if isempty(to), to = 0; end + next_state(from, ui) = to; + end + end + + % Initialize true path at n=1: pick any state with last_sym == data_ref(1) + % Prefer one whose suffix matches the first available history if L>1. + cand = find(last_sym == data_ref(1)); + if isempty(cand) + % fallback: choose closest in amplitude (should not happen if levels match) + [~,ix] = min(abs(last_sym - data_ref(1))); + cand = ix; + end + true_state_path = zeros(1,N,'uint32'); + true_state_path(1) = cand(1); + + % Propagate forward using the known inputs data_ref(n) + for n = 2:N + ui = ref_idx(n); % index of the actual transmitted level at time n + from = true_state_path(n-1); + if from==0 || ui==0 + true_state_path(n) = 0; + else + true_state_path(n) = next_state(from, ui); + if true_state_path(n)==0 + % Safety fallback: if no valid transition found (should not happen) + % choose any to-state whose vector matches shift+current symbol + target = [combs(from,2:end), levels(ui)]; + to = find(all(bsxfun(@eq, combs, target), 2), 1, 'first'); + if isempty(to), to = from; end + true_state_path(n) = to; + end + end + end + + %% ----- Collect branch metrics along decoded vs. true path ----- + bm_decoded = nan(1,N); + bm_true = nan(1,N); + + % n=1 in your code stores pm into bm_fw(:,:,1); real BMs start at n>=2 + for n = 2:N + % Decoded path: to = viterbi_path(n), from = survivor that fed it + to_d = viterbi_path(n); + from_d = pm_survivor_fw_idx(to_d, n); + bm_decoded(n) = bm_fw(to_d, from_d, n); + + % True path: transition true_state_path(n-1) -> true_state_path(n) + to_t = true_state_path(n); + from_t = true_state_path(n-1); + if to_t>0 && from_t>0 + bm_true(n) = bm_fw(to_t, from_t, n); + end + end + + % Convert to "cost" for intuitive plotting (your BM is a log-likelihood) + cost_dec = -bm_decoded; + cost_true= -bm_true; + + %% ----- Plot a short window for clarity ----- + win = max(2, N-20000):N; % last 200 samples (adjust as needed) + figure('Color','w'); hold on; grid on; box on; + plot(win, cost_true(win), 'LineWidth',1.2, 'DisplayName','True path cost (−BM)'); + plot(win, cost_dec(win), 'LineWidth',1.2, 'DisplayName','Decoded path cost (−BM)'); + xlabel('Time index n'); ylabel('Branch cost'); title('Branch metrics along true vs decoded path'); + legend('Location','best'); + + %% ----- Optional: overlay symbol levels for the same window ----- + yyaxis right + plot(win, data_in(win), ':', 'LineWidth',0.8, 'DisplayName','y(n)'); + ylabel('Amplitude'); + end VITERBI_ESTIMATION_SYMBOLS(1:length(data_in)) = first_sym(viterbi_path); @@ -564,11 +657,6 @@ classdef MLSE < handle end - function s = logsumexp(a,dim) - % returns log(sum(exp(a),dim)) safely - amax = max(a,[],dim); - s = amax + log(sum(exp(a - amax), dim)); - end end end diff --git a/Classes/GifWriter.m b/Classes/GifWriter.m new file mode 100644 index 0000000..c6483db --- /dev/null +++ b/Classes/GifWriter.m @@ -0,0 +1,146 @@ +classdef GifWriter < handle + %GIFWRITER Simple class to create GIFs from figures (parallel-safe) + % + % Example: + % g = GifWriter('Name','mySim','Parallel',true); + % parfor i = 1:10 + % plot(rand(10,1)); + % g.addFrame(1,i); + % end + % g.compile(1); + + properties + Name (1,:) char = 'default' % GIF base name + DelayTime (1,1) double = 0.1 % Frame delay in seconds + Parallel (1,1) logical = false % Enable parallel-safe mode + BaseDir (1,:) char % Base directory for temp frames + OutputDir (1,:) char % Final GIF output directory + end + + methods + %% Constructor + function obj = GifWriter(varargin) + % Parse name/value pairs + p = inputParser; + addParameter(p, 'Name', 'default', @ischar); + addParameter(p, 'DelayTime', 0.1, @isnumeric); + addParameter(p, 'Parallel', false, @islogical); + addParameter(p, 'OutputDir', fullfile(pwd, 'gif_output'), @ischar); + parse(p, varargin{:}); + + obj.Name = p.Results.Name; + obj.DelayTime = p.Results.DelayTime; + obj.Parallel = p.Results.Parallel; + obj.OutputDir = p.Results.OutputDir; + obj.BaseDir = fullfile(obj.OutputDir, 'tmp', obj.Name); + + if ~exist(obj.BaseDir, 'dir'), mkdir(obj.BaseDir); end + if ~exist(obj.OutputDir, 'dir'), mkdir(obj.OutputDir); end + end + + %% + function addFrame(obj, figInput, pos) + %ADDFRAME Add a figure frame to the GIF (supports parallel mode) + % + % Usage: + % obj.addFrame(figHandle) + % obj.addFrame(figNum) + % obj.addFrame(figHandle, pos) % parallel mode + % obj.addFrame(figNum, pos) + % + % In parallel mode, 'pos' must be a unique integer (loop index). + + if nargin < 3, pos = []; end + + % --- Resolve figure handle --- + if isnumeric(figInput) + % User passed a figure number + if ~ishandle(figInput) + warning('GifWriter:addFrame', 'Figure %d not found.', figInput); + return; + end + figHandle = figure(figInput); + elseif isa(figInput, 'matlab.ui.Figure') + figHandle = figInput; + else + error('GifWriter:addFrame:InvalidInput', ... + 'Input must be a figure handle or figure number.'); + end + + % --- Parallel-safe frame writing --- + if obj.Parallel + if isempty(pos) + error('GifWriter:ParallelMode', ... + 'In parallel mode, provide a unique ''pos'' identifier.'); + end + + % Directory for this figure number + frameDir = fullfile(obj.BaseDir, sprintf('fig_%d', figHandle.Number)); + if ~exist(frameDir, 'dir') + mkdir(frameDir); + end + + % File path for this frame + frameFile = fullfile(frameDir, sprintf('frame_%05d.png', pos)); + + % Export to PNG (headless-safe) + exportgraphics(figHandle, frameFile, 'Resolution', 150); + + else + % --- Serial mode: append directly to GIF --- + gifFile = fullfile(obj.OutputDir, ... + sprintf('%s_fig_%d.gif', obj.Name, figHandle.Number)); + + % Export frame temporarily + tmpFile = [tempname, '.png']; + exportgraphics(figHandle, tmpFile, 'Resolution', 150); + img = imread(tmpFile); + delete(tmpFile); + + % Append to GIF + [A, map] = rgb2ind(img, 256); + if ~isfile(gifFile) + imwrite(A, map, gifFile, 'gif', ... + 'LoopCount', Inf, 'DelayTime', obj.DelayTime); + else + imwrite(A, map, gifFile, 'gif', ... + 'WriteMode', 'append', 'DelayTime', obj.DelayTime); + end + end + end + + + + %% Compile all PNGs into a GIF (and clean up) + function compile(obj, fignum) + figDir = fullfile(obj.BaseDir, sprintf('fig_%d', fignum)); + gifFile = fullfile(obj.OutputDir, sprintf('%s_fig_%d.gif', obj.Name, fignum)); + + frames = dir(fullfile(figDir, 'frame_*.png')); + if isempty(frames) + warning('GifWriter:NoFrames', 'No frames found for figure %d.', fignum); + return; + end + + % Sort by frame name + [~, idx] = sort({frames.name}); + frames = frames(idx); + + % Combine into a GIF + for i = 1:numel(frames) + img = imread(fullfile(frames(i).folder, frames(i).name)); + [A, map] = rgb2ind(img, 256); + if i == 1 + imwrite(A, map, gifFile, 'gif', ... + 'LoopCount', Inf, 'DelayTime', obj.DelayTime); + else + imwrite(A, map, gifFile, 'gif', ... + 'WriteMode', 'append', 'DelayTime', obj.DelayTime); + end + end + + % Clean up temporary frames + rmdir(figDir, 's'); + end + end +end diff --git a/Functions/EQ_structures/dsp_runid.m b/Functions/EQ_structures/dsp_runid.m index 132cc75..1a20686 100644 --- a/Functions/EQ_structures/dsp_runid.m +++ b/Functions/EQ_structures/dsp_runid.m @@ -96,10 +96,10 @@ try adaption= 1; use_dd_mode = 1; - use_ffe = 1; - use_dfe = 1; + use_ffe = 0; + use_dfe = 0; use_vnle_mlse = 1; - use_dbtgt = 1; + use_dbtgt = 0; use_dbenc = 0; addProcessingResultToDatabase = 0; diff --git a/Functions/EQ_visuals/showEQNoisePSD.m b/Functions/EQ_visuals/showEQNoisePSD.m index 64f002d..8eeed21 100644 --- a/Functions/EQ_visuals/showEQNoisePSD.m +++ b/Functions/EQ_visuals/showEQNoisePSD.m @@ -43,6 +43,6 @@ end end xlim([-eq_noise.fs/2* 1e-9 eq_noise.fs/2* 1e-9]); - ylim([-15, 0]); + % ylim([-15, 0]); end diff --git a/Functions/Theory/CCDM/ccdm_gpt_example.m b/Functions/Theory/CCDM/ccdm_gpt_example.m new file mode 100644 index 0000000..bea3bf3 --- /dev/null +++ b/Functions/Theory/CCDM/ccdm_gpt_example.m @@ -0,0 +1,27 @@ +%% Target source entropy for PS-PAM8 +clear; clc; + +M = 8; +a = -(M-1):2:(M-1); % PAM-8 amplitude levels: [-7 -5 -3 -1 1 3 5 7] +H_target = 2.79; % desired entropy [bits/symbol] + +% Objective: find nu such that H(PA) = H_target +f = @(nu) entropy_MB(a,nu) - H_target; +nu_opt = fzero(f, [0, 2]); % search ν in reasonable range + +% Compute final distribution +P = exp(-nu_opt*a.^2); +P = P/sum(P); +H = -sum(P .* log2(P)); + +fprintf('Shaping parameter ν = %.4f\n', nu_opt); +fprintf('Entropy H(A) = %.3f bits/symbol\n', H); +disp('Probability vector (P_A):'); +disp(P.'); + +%% Helper: entropy function +function H = entropy_MB(a,nu) + P = exp(-nu*a.^2); + P = P/sum(P); + H = -sum(P .* log2(P)); +end diff --git a/Functions/beautifyBERplot.m b/Functions/beautifyBERplot.m index 0d0a665..9e799e2 100644 --- a/Functions/beautifyBERplot.m +++ b/Functions/beautifyBERplot.m @@ -8,12 +8,13 @@ function beautifyBERplot() for i = 1:length(lines) lines(i).LineWidth = 1.3; % Thicker line width - %lines(i).LineStyle = '-'; % Solid lines for simplicity - % if string(lines(i).Marker) == "none" - % lines(i).Marker = markers{mod(i-1, num_markers) + 1}; % Assign markers cyclically - % end - lines(i).MarkerSize = 4; % Marker size - lines(i).MarkerFaceColor = 'auto'; % Use line color for marker face + lines(i).LineStyle = '-'; % Solid lines for simplicity + if string(lines(i).Marker) == "none" + lines(i).Marker = markers{mod(i-1, num_markers) + 1}; % Assign markers cyclically + end + lines(i).MarkerSize = 7; % Marker size + lines(i).MarkerFaceColor = lines(i).Color; % Use line color for marker face + lines(i).MarkerEdgeColor = 'white'; end % Change all text interpreters to LaTeX diff --git a/projects/HighSpeedExperiment_2024/Auswertung_JLT/run_dsp_from_db.m b/projects/HighSpeedExperiment_2024/Auswertung_JLT/run_dsp_from_db.m index 8da64d8..ff9eac0 100644 --- a/projects/HighSpeedExperiment_2024/Auswertung_JLT/run_dsp_from_db.m +++ b/projects/HighSpeedExperiment_2024/Auswertung_JLT/run_dsp_from_db.m @@ -40,18 +40,18 @@ end fp = QueryFilter(); % fp.where('Runs', 'run_id','EQUALS', 987); -M = 6; -% fp.where('Runs', 'pam_level','EQUALS', M); -% fp.where('Runs', 'bitrate','EQUALS', 480e9); +M = 4; +fp.where('Runs', 'pam_level','EQUALS', M); +fp.where('Runs', 'bitrate','EQUALS', 420e9);%360,390 % fp.where('Runs', 'symbolrate','EQUALS', 162e9); -% fp.where('Runs', 'fiber_length','EQUALS', 1); +fp.where('Runs', 'fiber_length','EQUALS', 2); fp.where('Runs', 'is_mpi','EQUALS', 0); % fp.where('Runs', 'interference_path_length','EQUALS', 1000); % fp.where('Runs', 'loop_id','GREATER_THAN', 11); % fp.where('Runs', 'sir','EQUALS',18); -fp.where('Runs', 'wavelength','LESS_THAN', 1311); -% fp.where('Runs', 'db_mode','EQUALS', 0); -% fp.where('Runs', 'rop_attenuation','NOT_EQUAL', 0); +fp.where('Runs', 'wavelength','EQUAL', 1310); +fp.where('Runs', 'db_mode','EQUALS', 0); +fp.where('Runs', 'rop_attenuation','EQUAL', 0); % fp.where('Runs', 'power_pd_in','GREATER_THAN', 7); [dataTable,~] = db.queryDB(fp, db.getTableFieldNames('Runs')); @@ -69,7 +69,7 @@ wh.addStorage("dbenc_package"); % === RUN IT === -[results,wh] = submitJobs(dataTable.run_id(:), dsp_options, "parallel", 'wh', wh, 'waitbar', true); +[results,wh] = submitJobs(dataTable.run_id(:), dsp_options, "serial", 'wh', wh, 'waitbar', true); results_db = results(dataTable.db_mode==1); results_nodb = results(dataTable.db_mode==0); diff --git a/projects/HighSpeedExperiment_2024/Auswertung_JLT/run_plot_measurements.m b/projects/HighSpeedExperiment_2024/Auswertung_JLT/run_plot_measurements.m index e7dd62f..ebcfef1 100644 --- a/projects/HighSpeedExperiment_2024/Auswertung_JLT/run_plot_measurements.m +++ b/projects/HighSpeedExperiment_2024/Auswertung_JLT/run_plot_measurements.m @@ -6,7 +6,7 @@ M = 4; fp = QueryFilter(); % fp.where('Runs', 'run_id','EQUALS', 987); fp.where('Runs', 'pam_level','EQUALS', M); -fp.where('Runs', 'symbolrate','EQUALS', 150e9); +% fp.where('Runs', 'symbolrate','EQUALS', 150e9); % fp.where('Runs', 'fiber_length','EQUALS', 10); fp.where('Runs', 'is_mpi','EQUALS', 0); % fp.where('Runs', 'interference_path_length','EQUALS', 1000); diff --git a/projects/IMDD_base_system/simulation_bwl.m b/projects/IMDD_base_system/simulation_bwl.m index 79fc4c7..1ac779f 100644 --- a/projects/IMDD_base_system/simulation_bwl.m +++ b/projects/IMDD_base_system/simulation_bwl.m @@ -47,16 +47,17 @@ doub_mode = db_mode.no_db; cols = linspecer(6); rop = [-6]; bwl = [0.5:0.1:1.5]; -fsym = [120:8:256].*1e9; -fsym =210e9; +fsym = [192:16:256].*1e9; +fsym = 208e9; -ber_vnle = []; +ber_vnle = []; ber_mlse = []; +ber_mlse_burg = []; ber_viterbi = []; ber_db = []; ber_db_diff_precoded = []; gmi_vnle_bitwise = []; -gmi_mlse = []; +gmi_mlse = []; gmi_mlse_db = []; for r = 1:length(fsym) @@ -69,16 +70,16 @@ for r = 1:length(fsym) apply_pulsef = 1; [Digi_sig,Symbols,Tx_bits] = PAMsource(... - "fsym",fsym(r),"M",M,"order",18,"useprbs",0,... + "fsym",fsym(r),"M",M,"order",19,"useprbs",0,... "fs_out",fdac,... "applyclipping",0,"clipfactor",1.5,... "applypulseform",apply_pulsef,"pulseformer",Pform,... "randkey",random_key,... "db_precode",db_precode,"db_encode",db_encode,... "mrds_code",0,"mrds_blocklength",512,"duobinary_mode",duob_mode).process(); - + % El_sig = AWG("fdac",fdac,"f_cutoff",fsym(r),"lpf_active",0,"kover",kover,"bit_resolution",12,"upsampling_method","samplehold","precomp_sinc_rolloff",0).process(Digi_sig); - El_sig = M8199B("kover",kover).process(Digi_sig); + El_sig = M8199B("kover",kover).process(Digi_sig); % AWG("fdac",fdac,"f_cutoff",fsym(r),"lpf_active",0,"kover",kover,"bit_resolution",12,"upsampling_method","samplehold","precomp_sinc_rolloff",0).process(Digi_sig); %%%%% Low-pass el. components %%%%%% @@ -95,15 +96,17 @@ for r = 1:length(fsym) vbias = -u_pi*0.5; [Opt_sig] = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs,"lambda",laser_wavelength,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",random_key+1).process(El_sig); - figure(15); - hold on - scatter(El_sig.signal(1:100000)+vbias,(abs(Opt_sig.signal(1:100000)).^2)*1e3,0.1,'.','DisplayName','Modulator TF') - xlabel('Input in V') - ylabel('abs(Eopt)2 in mW','Interpreter','latex') - ylim([0 2]); - xlim([-3.2 0]); - - Opt_sig.eye(fsym(r),M,"fignum",103837); + if 0 + figure(15); + hold on + scatter(El_sig.signal(1:100000)+vbias,(abs(Opt_sig.signal(1:100000)).^2)*1e3,0.1,'.','DisplayName','Modulator TF') + xlabel('Input in V') + ylabel('abs(Eopt)2 in mW','Interpreter','latex') + ylim([0 2]); + xlim([-3.2 0]); + + Opt_sig.eye(fsym(r),M,"fignum",103837); + end %%%%%% Fiber %%%%%% Opt_sig = Fiber("fsimu",Opt_sig.fs,"fiber_length",link_length/1000,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.07).process(Opt_sig); @@ -113,50 +116,50 @@ for r = 1:length(fsym) % Opt_sig.eye(fsym(r),M,"fignum",103838); - % % Opt_sig.signal = Opt_sig.signal + 5*abs(mean(Opt_sig.signal)); + % % Opt_sig.signal = Opt_sig.signal + 5*abs(mean(Opt_sig.signal)); % Opt_sig.move_it_spectrum("displayname",'Opt Sig after Amp','fignum',1223323); % Pc = abs(mean(Opt_sig.signal)).^2; % carrier power % Ptot = mean(abs(Opt_sig.signal).^2); % total power - % Ps = max(Ptot - Pc, eps); + % Ps = max(Ptot - Pc, eps); % Pcdb = 10*log10(Pc); % Psdb = 10*log10(Ps); - % + % % cspr_dB = 10*log10(Pc / Ps); - % + % % % Minimal in-place CSPR set (real, nonnegative field constraint) % E = Opt_sig.signal; % real field samples % target_cspr_dB = 20; % <-- set your target CSPR (dB) - % + % % % Decompose into DC + zero-mean waveform % m = mean(E); % x0 = E - m; % zero-mean modulation % Ps0 = mean(x0.^2); % sideband power (fixed if shape kept) - % + % % % Current CSPR (for reference) % Pc_cur = m^2; % Ptot_cur = mean(E.^2); % Ps_cur = max(Ptot_cur - Pc_cur, eps); % cspr_in = 10*log10(Pc_cur / Ps_cur); - % + % % % Bias needed for target CSPR, and minimal bias to keep E>=0 % R_tgt = 10^(target_cspr_dB/10); % Pc/Ps % a_req = sqrt(R_tgt * Ps0); % required DC bias % a_min = -min(x0); % to avoid negatives everywhere % a = max(a_req, a_min); % if infeasible, lands at CSPR_min - % + % % % Apply bias (preserves waveform shape) % E_new = a + x0; - % + % % % Achieved CSPR % Pc_new = mean(E_new)^2; % Ptot_new = mean(E_new.^2); % Ps_new = max(Ptot_new - Pc_new, eps); % cspr_out = 10*log10(Pc_new / Ps_new); - % + % % % (Optional) show feasibility info % cspr_min = 10*log10((a_min^2)/max(Ps0,eps)); % disp(table(cspr_in, target_cspr_dB, cspr_min, cspr_out)); - % + % % % Use E_new as your adjusted field % Opt_sig.signal = E_new; @@ -166,57 +169,56 @@ for r = 1:length(fsym) %%%%%% Low-pass RX (PD, El. Connectors and Scope %%%%%% rx_bwl = 70e9; PD_sig = Filter('filtdegree',4,"f_cutoff",rx_bwl,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(PD_sig); - + % %%%%%% Low-pass Scope %%%%%% Lp_scpe = Filter('filtdegree',4,"f_cutoff",110e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true); - + %%%%%% Scope %%%%%% Scpe_sig = Scope("fsimu",fdac*kover,"fadc",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',Lp_scpe).process(PD_sig); - + Scpe_sig_2sps = Scpe_sig.resample("fs_out",2*fsym(r)); % Scpe_sig_resampled.signal = Scpe_sig_resampled.signal(1:2*length(Symbols)); - + [~, Scpe_cell, ~, found_sync] = Scpe_sig_2sps.tsynch("reference", Symbols, "fs_ref", fsym(r), "debug_plots", 0); Rx_sig = Scpe_cell{1}; Rx_sig = Rx_sig.normalize("mode","rms"); - if 0 + if 1 %Duobinary Targeting eq_ = EQ("Ne",[vnle_order1,vnle_order2,vnle_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); - db_ref_sequence = Duobinary().encode(Symbols); db_ref_constellation = unique(db_ref_sequence.signal); [eq_signal, eq_noise] = eq_.process(Rx_sig,db_ref_sequence); - + viterbi = 0; - if viterbi + if viterbi mlse_ = MLSE_viterbi("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); mlse_.DIR = [1,1]; - [mlse_sig_sd] = mlse_.process(eq_signal); + [eq_signal_whitened] = mlse_.process(eq_signal); else mlse_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",M,"trellis_states",PAMmapper(M,0).get_levels ./ PAMmapper(M,0).get_scaling); mlse_.DIR = [1,1]; - [mlse_sig_sd,LLR,gmi_mlse_db(r)] = mlse_.process(eq_signal,Symbols); + [eq_signal_whitened,LLR,gmi_mlse_db(r)] = mlse_.process(eq_signal,Symbols); end - - mlse_sig_hd = PAMmapper(M,0,"eth_style",0).quantize(mlse_sig_sd); + + mlse_sig_hd = PAMmapper(M,0,"eth_style",0).quantize(eq_signal_whitened); mlse_sig_hd_precoded = Duobinary().encode(mlse_sig_hd,"M",M); mlse_sig_hd_precoded = Duobinary().decode(mlse_sig_hd_precoded,"M",M); - + tx_symbols_precoded = Duobinary().encode(Symbols); tx_symbols_precoded = Duobinary().decode(tx_symbols_precoded); - + tx_bits_precoded = PAMmapper(M,0,"eth_style",0).demap(tx_symbols_precoded); - + rx_bits_mlse = PAMmapper(M,0,"eth_style",0).demap(mlse_sig_hd_precoded); [~,errors_db_diff_precoded,ber_db_diff_precoded(r),a] = calc_ber(rx_bits_mlse.signal,tx_bits_precoded.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); burst_db_pre(r,:) = count_error_bursts(a, 15)./numel(Tx_bits.signal); - + %B) Just determine BER rx_bits_mlse = PAMmapper(M,0,"eth_style",0).demap(mlse_sig_hd); [bits_mlse,errors_db,ber_db(r),a] = calc_ber(rx_bits_mlse.signal,Tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); @@ -225,39 +227,37 @@ for r = 1:length(fsym) fprintf('BER ber_db_diff_precoded: %.2e \n',ber_db_diff_precoded(r)); fprintf('BER Vber_dbNLE: %.2e \n',ber_db(r)); % figure();hold on;stem(1:15,burst_db(r,:),'LineWidth',1,'Color',cols(1,:));stem(1:15,burst_db_pre(r,:),'LineWidth',1,'Color',cols(2,:));set(gca, 'yscale', 'log'); - + end % FFE or VNLE eq_ = EQ("Ne",[vnle_order1,vnle_order2,vnle_order3],"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.00,"FFEmu",0,"plotfinal",0,"ideal_dfe",0); % eq = VNLE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",[0.0004 0.0005 0.0006],"mu_tr",0,"order",[50,2,2],"sps",2,"decide",0); - - [eq_signal_sd, eq_noise] = eq_.process(Rx_sig, Symbols); + + [eq_signal_fullresp, eq_noise] = eq_.process(Rx_sig, Symbols); showEQNoisePSD(eq_noise, "fignum",1273876,"displayname",'noise after EQ'); - [mi_gomez(r)] = calc_air(eq_signal_sd, Symbols, "skip_front", 100, "skip_end", 100); - [gmi_vnle_bitwise(r)] = calc_ngmi(eq_signal_sd,Symbols); - [gmi_bitwise_2(r)] = calc_gmi_bitwise(eq_signal_sd,Symbols); - snr_vnle(r) = calc_snr(Symbols, eq_signal_sd-Symbols); - - eq_signal_sd.plot("displayname",'bla','fignum',199); - eq_signal_sd.eye(fsym(r),M,"fignum",103837); - + [mi_gomez(r)] = calc_air(eq_signal_fullresp, Symbols, "skip_front", 100, "skip_end", 100); + [gmi_vnle_bitwise(r)] = calc_ngmi(eq_signal_fullresp,Symbols); + [gmi_bitwise_2(r)] = calc_gmi_bitwise(eq_signal_fullresp,Symbols); + snr_vnle(r) = calc_snr(Symbols, eq_signal_fullresp-Symbols); + + % eq_signal_fullresp.plot("displayname",'bla','fignum',199); + % eq_signal_fullresp.eye(fsym(r),M,"fignum",103837); % Hard decision on VNLE output - eq_signal_hd = PAMmapper(M, 0).quantize(eq_signal_sd); + eq_signal_hd = PAMmapper(M, 0).quantize(eq_signal_fullresp); rx_bits = PAMmapper(M,0,"eth_style",0).demap(eq_signal_hd); [~,tot_err,ber_vnle(r),a] = calc_ber(rx_bits.signal,Tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); burst_vnle(r,:) = count_error_bursts(a, 10)./tot_err; - showLevelConfusionMatrix(eq_signal_hd,Symbols,"M",M,"fignum",200,"displayname",'bla'); - showLevelScatter(eq_signal_sd,Symbols,"displayname",'VNLE Out','f_sym',fsym(r),'fignum',201); - show2Dconstellation(eq_signal_sd,Symbols,"displayname",'VNLE Out','fignum',2241); - + % showLevelConfusionMatrix(eq_signal_hd,Symbols,"M",M,"fignum",200,"displayname",'bla'); + % showLevelScatter(eq_signal_fullresp,Symbols,"displayname",'VNLE Out','f_sym',fsym(r),'fignum',201); + % show2Dconstellation(eq_signal_fullresp,Symbols,"displayname",'VNLE Out','fignum',2241); + fprintf('BER VNLE: %.2e \n',ber_vnle(r)); fprintf('NGMI VNLE: %.2f \n',gmi_vnle_bitwise(r)./m); - if 1 % Process through postfilter and MLSE @@ -268,45 +268,133 @@ for r = 1:length(fsym) pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1,"coefficients",[1,0.85]); end - mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).get_levels ./ PAMmapper(M,0).get_scaling); - [mlse_sig_sd,whitened_noise] = pf_.process(eq_signal_sd, eq_noise); - - mlse_.DIR = pf_.coefficients; + % showEQNoisePSD(eq_noise,"postfilter_taps",pf_.coefficients,"displayname",'Postfilter Burg based'); alpha(r) = pf_.coefficients(2); - [signalclass_hd,LLR,gmi_mlse(r)] = mlse_.process(mlse_sig_sd,Symbols); - mlse_sig_hd = PAMmapper(M, 0, "eth_style", 0).quantize(signalclass_hd); - rx_bits = PAMmapper(M,0,"eth_style",0).demap(mlse_sig_hd); + alpha_vec = max(0,round(alpha(r),2)-0.2):0.025:round(alpha(r),2)+0.4; + alpha_vec = unique(sort([alpha_vec, 1, alpha(r)])); - [~,tot_err,ber_mlse(r),a] = calc_ber(rx_bits.signal,Tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - burst_mlse(r,:) = count_error_bursts(a, 10); + gmi_mlse_ = zeros(size(alpha_vec)); + ber_mlse_ = zeros(size(alpha_vec)); + parfor a=1:numel(alpha_vec) - showLevelConfusionMatrix(mlse_sig_hd,Symbols,"M",M,"fignum",300,"displayname",'bla'); + mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).get_levels ./ PAMmapper(M,0).get_scaling,'DIR',[1,alpha_vec(a)]); - fprintf('BER MLSE: %.2e \n',ber_mlse(r)); - fprintf('NGMI MLSE: %.5f \n',gmi_mlse(r)./m); + pf_ = Postfilter("ncoeff",1,"useBurg",0,"coefficients",[1,alpha_vec(a)]); + [eq_signal_whitened,whitened_noise] = pf_.process(eq_signal_fullresp, eq_noise); + + [signalclass_hd,LLR,gmi_mlse_(a)] = mlse_.process(eq_signal_whitened,Symbols); - levels = sort(unique(Symbols.signal(:)).'); % 1×6 - pairs = reshape(mlse_sig_hd.signal,2,[]).'; - isedge = ismember(pairs, [levels(1) levels(end)]); - isforbidden = sum(isedge,2)==2; - fprintf('Found %d forbidden transitions (even→odd edges).\n', nnz(isforbidden)); + mlse_sig_hd = PAMmapper(M, 0, "eth_style", 0).quantize(signalclass_hd); + + rx_bits = PAMmapper(M,0,"eth_style",0).demap(mlse_sig_hd); + + [~,tot_err,ber_mlse_(a),errpos] = calc_ber(rx_bits.signal,Tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + % burst_mlse(r,:) = count_error_bursts(errpos, 10); + + % if 0 + % fprintf('BER MLSE: %.2e \n',ber_mlse(r)); + % fprintf('NGMI MLSE: %.5f \n',gmi_mlse(r)./m); + % + % showLevelConfusionMatrix(mlse_sig_hd,Symbols,"M",M,"fignum",300,"displayname",'bla'); + % + % levels = sort(unique(Symbols.signal(:)).'); % 1×6 + % pairs = reshape(mlse_sig_hd.signal,2,[]).'; + % isedge = ismember(pairs, [levels(1) levels(end)]); + % isforbidden = sum(isedge,2)==2; + % fprintf('Found %d forbidden transitions (even→odd edges).\n', nnz(isforbidden)); + % + % + % % Process through postfilter and MLSE + % pf_ncoeffs = 1; + % pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); + % mlse_ = MLSE_viterbi("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); + % [eq_signal_whitened,whitened_noise] = pf_.process(eq_signal_fullresp, eq_noise); + % mlse_.DIR = pf_.coefficients; + % mlse_output = mlse_.process(eq_signal_whitened); + % mlse_sig_hd = PAMmapper(M, 0, "eth_style", 0).quantize(mlse_output); + % rx_bits = PAMmapper(M,0,"eth_style",0).demap(mlse_sig_hd); + % [~,~,ber_viterbi(r),~] = calc_ber(rx_bits.signal,Tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + % fprintf('Viterbi BER: %.2e \n',ber_viterbi(r)); + % end + end + + [ber_mlse(r),idx] = min(ber_mlse_); + gmi_mlse(r) = gmi_mlse_(idx); + ber_mlse_burg(r) = ber_mlse_(alpha_vec==alpha(r)); + best_alpha(r) = alpha_vec(idx); + + end + + % IR target in EQ + if 1 + + % alpha_vec = max(0,round(alpha(r),2)-0.1):0.01:min(1,round(alpha(r),2)+0.1); + plot_stuff = 0; + gmi_mlse_pr_tgt_ = zeros(size(alpha_vec)); + ber_mlse_pr_tgt_ = zeros(size(alpha_vec)); + parfor a = 1:numel(alpha_vec) + + eq_ = EQ("Ne",[vnle_order1,vnle_order2,vnle_order3],"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.00,"FFEmu",0,"plotfinal",0,"ideal_dfe",0); + Symbols_filt = Symbols.filter([1,alpha_vec(a)],1); + [eq_signal_prtgt, eq_noise] = eq_.process(Rx_sig, Symbols_filt); - % Process through postfilter and MLSE - pf_ncoeffs = 1; - pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); - mlse_ = MLSE_viterbi("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); - [mlse_sig_sd,whitened_noise] = pf_.process(eq_signal_sd, eq_noise); - mlse_.DIR = pf_.coefficients; - mlse_sig_sd = mlse_.process(mlse_sig_sd); - mlse_sig_hd = PAMmapper(M, 0, "eth_style", 0).quantize(mlse_sig_sd); - rx_bits = PAMmapper(M,0,"eth_style",0).demap(mlse_sig_hd); - [~,~,ber_viterbi(r),~] = calc_ber(rx_bits.signal,Tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - fprintf('Viterbi BER: %.2e \n',ber_viterbi(r)); + if plot_stuff + % Plot the response for respective EQ targets + Symbols_filt.spectrum("displayname",'IDEAL Filtered Reference','fignum',240587); + eq_signal_whitened.spectrum("displayname",'Full tgt. EQ + PF','fignum',240587); + eq_signal_prtgt.spectrum("displayname",'Partial Resp. Target EQ','fignum',240587); + + noise_pf_out = Symbols_filt-eq_signal_whitened; + noise_pr_tgt = Symbols_filt-eq_signal_prtgt; + + noise_pf_out.spectrum("displayname",'Ideal PR - Whitening Out','fignum',240588); + noise_pr_tgt.spectrum("displayname",'Ideal PR - PR Target Out','fignum',240588); + end + + mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).get_levels ./ PAMmapper(M,0).get_scaling,'DIR',[1,alpha_vec(a)],'debug',0); + + [signalclass_hd,LLR,gmi_mlse_pr_tgt_(a)] = mlse_.process(eq_signal_prtgt,Symbols); + + mlse_sig_hd = PAMmapper(M, 0, "eth_style", 0).quantize(signalclass_hd); + + rx_bits = PAMmapper(M,0,"eth_style",0).demap(mlse_sig_hd); + + [~,tot_err,ber_mlse_pr_tgt_(a),errpos] = calc_ber(rx_bits.signal,Tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + + % burst_mlse_(a,:) = count_error_bursts(errpos, 10); + + % fprintf('BER MLSE: %.2e \n',ber_mlse_pr_tgt_(a)); + % fprintf('NGMI MLSE: %.5f \n',gmi_mlse_pr_tgt_(a)./m); + + end + + [ber_mlse_pr_tgt(r),idx] = min(ber_mlse_pr_tgt_); + gmi_mlse_pr_tgt(r) = gmi_mlse_pr_tgt_(idx); + best_alpha_pr_tgt(r) = alpha_vec(idx); end + cols = cbrewer2('paired',8); + figure(); hold on + title(sprintf('%d GBd',fsym(r).*1e-9)); + scatter(alpha_vec,ber_mlse_,15,'Marker','o','LineWidth',1,'DisplayName','MLSE','MarkerEdgeColor',cols(1,:)); + scatter(best_alpha(r),ber_mlse(r),15,'Marker','o','LineWidth',2,'DisplayName','MLSE','MarkerEdgeColor',cols(2,:)); + scatter(alpha(r),ber_mlse_burg(r),25,'Marker','+','LineWidth',2,'DisplayName','MLSE','MarkerEdgeColor',cols(2,:)); + + scatter(1,ber_db_diff_precoded(r),15,'Marker','diamond','LineWidth',2,'DisplayName','Duobinary','MarkerEdgeColor',cols(4,:)); + scatter(1,ber_db(r),15,'Marker','diamond','LineWidth',2,'DisplayName','Duobinary','MarkerEdgeColor',cols(4,:)); + + scatter(alpha_vec,ber_mlse_pr_tgt_,15,'Marker','x','LineWidth',1,'DisplayName','MLSE Partial Resp tgt','MarkerEdgeColor',cols(5,:)); + scatter(best_alpha_pr_tgt(r),ber_mlse_pr_tgt(r),25,'Marker','x','LineWidth',2,'DisplayName','MLSE','MarkerEdgeColor',cols(6,:)); + + set(gca,"YScale","log"); + % ylim([1e-6 0.5]); + % xlim([0.1 1]); + drawnow; + + end @@ -413,9 +501,9 @@ m = floor(log2(M)*10)/10; figure(113+M); clf; hold on; netrates_vnle = tp.calculateNetRate(fsym.* m, ... - 'NGMI', gmi_vnle_bitwise./m, ... - 'BER', ber_vnle); -% + 'NGMI', gmi_vnle_bitwise./m, ... + 'BER', ber_vnle); +% plot(xGHz, gmi_vnle_bitwise.*xGHz, ... 'DisplayName','GMI*R VNLE', ... mk.VNLE{:}, 'LineStyle','-','LineWidth',LINE_WIDTH,'Color',cm.VNLE); @@ -428,15 +516,15 @@ plot(xGHz, netrates_vnle.HD.NetRate.*1e-9, ... mk.VNLE{:}, 'LineStyle','-.','LineWidth',LINE_WIDTH,'Color',cm.VNLE); -% +% % MLSE symbol-wise (if present) plot(xGHz, gmi_mlse.*xGHz, ... 'DisplayName','GMI*R MLSE', ... mk.MLSE{:}, 'LineStyle','-','LineWidth',LINE_WIDTH,'Color',cm.MLSE); netrates_mlse = tp.calculateNetRate(fsym.* m, ... - 'NGMI', gmi_mlse./m, ... - 'BER', ber_mlse); + 'NGMI', gmi_mlse./m, ... + 'BER', ber_mlse); plot(xGHz, netrates_mlse.SDHD.NetRate.*1e-9, ... 'DisplayName','SD+HD MLSE', ... mk.MLSE{:}, 'LineStyle',':','LineWidth',LINE_WIDTH,'Color',cm.MLSE); @@ -451,8 +539,8 @@ plot(xGHz, gmi_mlse_db.*xGHz, ... mk.DB{:}, 'LineStyle','-','LineWidth',LINE_WIDTH,'Color',cm.DB); netrates_db = tp.calculateNetRate(fsym.* m, ... - 'NGMI', gmi_mlse_db./m, ... - 'BER', ber_db); + 'NGMI', gmi_mlse_db./m, ... + 'BER', ber_db); plot(xGHz, netrates_db.SDHD.NetRate.*1e-9, ... 'DisplayName','SD+HD DB', ... @@ -473,7 +561,7 @@ xlim([184, 256]) % Auxiliary nested helper for numerically stable log-sum-exp function s = logsumexp(a) - % LOGSUMEXP Compute log(sum(exp(a))) in a numerically stable way - m = max(a); - s = m + log(sum(exp(a - m))); +% LOGSUMEXP Compute log(sum(exp(a))) in a numerically stable way +m = max(a); +s = m + log(sum(exp(a - m))); end diff --git a/projects/IMDD_base_system/simulation_bwl_2.m b/projects/IMDD_base_system/simulation_bwl_2.m index 8dff4a2..27a71cb 100644 --- a/projects/IMDD_base_system/simulation_bwl_2.m +++ b/projects/IMDD_base_system/simulation_bwl_2.m @@ -136,9 +136,5 @@ eq_signal_hd = PAMmapper(M, 0).quantize(eq_signal_sd); % Process through postfilter and MLSE [mlse_sig_sd,whitened_noise] = pf_.process(eq_signal_sd, eq_noise); mlse_.DIR = pf_.coefficients; - -constellation = [-3, -1, 1, 3]; -chatgpt_answer(mlse_sig_sd.signal, Symbols.signal,mlse_.DIR,constellation); - -% mlse_sig_sd = mlse_.process(mlse_sig_sd,Symbols); -% mlse_sig_hd = PAMmapper(M, 0, "eth_style", options.eth_style_symbol_mapping).quantize(mlse_sig_sd); +mlse_sig_sd = mlse_.process(mlse_sig_sd,Symbols); +mlse_sig_hd = PAMmapper(M, 0, "eth_style", 0).quantize(mlse_sig_sd); diff --git a/projects/ML_based_MLSE/model.m b/projects/ML_based_MLSE/model.m new file mode 100644 index 0000000..f74e5fe --- /dev/null +++ b/projects/ML_based_MLSE/model.m @@ -0,0 +1,194 @@ +%%% Run parameters +% TX +M = 4; + +apply_pulsef = 1; +fdac = 256e9; +fadc = 256e9; +random_key = 2; + +rcalpha = 0.05; +kover = 8; +vbias_rel = 0.5; +u_pi = 3.2; +vbias = -vbias_rel*u_pi; +laser_wavelength = 1310; +laser_linewidth = 1e6; + +% Channel +link_length = 0; + +alpha = 0; + +doub_mode = db_mode.no_db; +cols = linspecer(6); +rop = [-6]; +bwl = [0.5:0.1:1.5]; +fsym = [208:16:256].*1e9; +% nonlin_mod = [0.5:0.01:0.75]; +nonlin_mod = ones(size(fsym)).*0.5; + +ffe_results = {}; +mlse_results_lin= {}; + +for r = 1:length(fsym) + + Pform = Pulseformer("fsym",fsym(r),"fdac",4*fsym(r),"pulse","rc","pulselength",16,"alpha",rcalpha); + + db_precode = 0; + db_encode = 0; + duob_mode = db_mode.no_db; + apply_pulsef = 1; + + [Digi_sig,Symbols,Tx_bits] = PAMsource(... + "fsym",fsym(r),"M",M,"order",17,"useprbs",0,... + "fs_out",fdac,... + "applyclipping",0,"clipfactor",1.5,... + "applypulseform",apply_pulsef,"pulseformer",Pform,... + "randkey",random_key,... + "db_precode",db_precode,"db_encode",db_encode,... + "mrds_code",0,"mrds_blocklength",512,"duobinary_mode",duob_mode).process(); + + El_sig = M8199B("kover",kover).process(Digi_sig); + + %%%%% Electrical Driver Amplifier %%%%%% + El_sig = El_sig.normalize("mode","oneone"); + + %%%%% MODULATE E/O CONVERSION %%%%% + u_pi = 3.2; + vbias = -u_pi*nonlin_mod(r); + [Opt_sig] = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs,"lambda",laser_wavelength,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",random_key+1).process(El_sig); + + %%%%%% Fiber %%%%%% + Opt_sig = Fiber("fsimu",Opt_sig.fs,"fiber_length",link_length/1000,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.07).process(Opt_sig); + + %%%%%% ROP %%%%%% + Opt_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",rop).process(Opt_sig); + + %%%%%% PD Square Law %%%%%% + PD_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20,"nep",1.8e-11,"randomkey",random_key).process(Opt_sig); + + %%%%%% Low-pass RX (PD, El. Connectors and Scope %%%%%% + rx_bwl = 70e9; + PD_sig = Filter('filtdegree',4,"f_cutoff",rx_bwl,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(PD_sig); + + % %%%%%% Low-pass Scope %%%%%% + Lp_scpe = Filter('filtdegree',4,"f_cutoff",110e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true); + + %%%%%% Scope %%%%%% + Scpe_sig = Scope("fsimu",fdac*kover,"fadc",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',Lp_scpe).process(PD_sig); + + Scpe_sig_2sps = Scpe_sig.resample("fs_out",2*fsym(r)); + + % 2sps + [~, Scpe_cell, ~, found_sync] = Scpe_sig_2sps.tsynch("reference", Symbols, "fs_ref", fsym(r), "debug_plots", 1); + Rx_sig_2sps = Scpe_cell{1}; + Rx_sig_2sps = Rx_sig_2sps.normalize("mode","rms"); + + % 1sps + Scpe_sig_1sps = Scpe_sig.resample("fs_out",1*fsym(r)); + + [~, Scpe_cell_1sps, ~, found_sync] = Scpe_sig_1sps.tsynch("reference", Symbols, "fs_ref", fsym(r), "debug_plots", 1); + Rx_sig_1sps = Scpe_cell_1sps{1}; + Rx_sig_1sps = Rx_sig_1sps.normalize("mode","rms"); + + + %% Implement DSP directly here: + + mu_lms = 0.0005; + tic + eq = ML_MLSE("epochs_tr",5,"epochs_dd",5,"len_tr",2^13,"mu_dd",mu_lms,"mu_tr",mu_lms,"order",50,"sps",2,"traceback_depth",32,"L",1); + [Eq_signal,Vit_signal] = eq.process(Rx_sig_2sps,Symbols); + toc + + % tic + % eq = FFE_MLSE("epochs_tr",5,"epochs_dd",5,"len_tr",2^13,"mu_dd",mu_lms,"mu_tr",mu_lms,"order",50,"sps",2,"traceback_depth",32,"L",3); + % [Eq_signal,Vit_signal] = eq.process(Rx_sig_2sps,Symbols); + % toc + + Eq_bits = PAMmapper(M, 0, "eth_style", 0).demap(Eq_signal); + [~, errors, ber, ~] = calc_ber(Eq_bits.signal, Tx_bits.signal, "skip_front", 0, "skip_end", 0, "returnErrorLocation", 1); + fprintf('FFE: %.2e \n',ber); + + Vit_signal_ = Vit_signal; + Vit_signal_.signal = circshift(Vit_signal.signal,0); + Vit_bits = PAMmapper(M, 0, "eth_style", 0).demap(Vit_signal_); + [~, errors, ber, errpos] = calc_ber(Vit_bits.signal, Tx_bits.signal, "skip_front", 0, "skip_end", 0, "returnErrorLocation", 1); + fprintf('Viterbi: %.2e \n',ber); + + + + + + %% optimize smth. + tr_len = 2.^[2:15]; + tr_len = floor(tr_len); + ber = zeros(size(tr_len)); + parfor m = 1:numel(tr_len) + + mu_lms = 0.0005; + eq = ML_MLSE("epochs_tr",5,"epochs_dd",5,"len_tr",2^13,"mu_dd",mu_lms,"mu_tr",mu_lms,"order",50,"sps",2,"traceback_depth",tr_len(m)); + [Eq_signal,Vit_signal] = eq.process(Rx_sig_2sps,Symbols); + + % Eq_bits = PAMmapper(M, 0, "eth_style", 0).demap(Eq_signal); + % [~, errors, ber(m), ~] = calc_ber(Eq_bits.signal, Tx_bits.signal, "skip_front", 0, "skip_end", 0, "returnErrorLocation", 1); + + Vit_bits = PAMmapper(M, 0, "eth_style", 0).demap(Vit_signal); + [~, errors, ber(m), errpos] = calc_ber(Vit_bits.signal, Tx_bits.signal, "skip_front", 0, "skip_end", 0, "returnErrorLocation", 1); + end + + figure(5); hold on + title('1-SPS') + plot(tr_len,ber,'DisplayName','BER'); + xlabel('Traceback Length') + beautifyBERplot + legend + ylim([1e-4, 1e-1]); + set(gca,'YScale','log'); + + + + + + + + + + + + + + + + %% RUN Comparison + len_tr = 4096*2; + + mu_ffe1 = 0.0001; + mu_ffe2 = 0.0008; + mu_ffe3 = 0.001; + mu_dc = 0.005; + % mu_dc = 0; + + mu_ffe = [mu_ffe1 mu_ffe2 mu_ffe3]; + mu_dfe = 0.0004; + pf_ncoeffs = 1; + ffe_order = [50, 0, 0]; + mu_lms = 0.0005; + eq_ = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",2^13,"mu_dd",mu_lms,"mu_tr",mu_lms,"order",50,"sps",1,"dd_mode",1,"adaption_technique","lms"); + % eq_ = EQ("Ne",ffe_order,"Nb",[0,0,0],"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",1,"DCmu",0.00,"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,'scale_mode',0,'trellis_exclusion',0,'trellis_state_mode',2,'debug',0); + + [ffe_results{r}, mlse_results_lin{r}] = vnle_postfilter_mlse(eq_, pf_, mlse_, M, Rx_sig_1sps, Symbols, Tx_bits, ... + "precode_mode", duob_mode,... + 'showAnalysis', 0, ... + "postFFE", [],... + "eth_style_symbol_mapping", 0); + + mlse_results_lin{r}.metrics.print; + ffe_results{r}.metrics.print; + +end \ No newline at end of file diff --git a/projects/Nonlinear_MLSE/simulation_nonlin_dsp.m b/projects/Nonlinear_MLSE/simulation_nonlin_dsp.m new file mode 100644 index 0000000..1896dd4 --- /dev/null +++ b/projects/Nonlinear_MLSE/simulation_nonlin_dsp.m @@ -0,0 +1,229 @@ +%%% Run parameters +% TX +M = 4; +m = floor(log2(M)*10)/10; +fsym = 224e9; + +apply_pulsef = 1; +fdac = 256e9; +fadc = 256e9; +random_key = 2; + +rcalpha = 0.05; +kover = 8; +vbias_rel = 0.5; +u_pi = 3.2; +vbias = -vbias_rel*u_pi; +laser_wavelength = 1310; +laser_linewidth = 1e6; + + +% Channel +link_length = 0; + +vnle_order1 = 50; +vnle_order2 = 0; +vnle_order3 = 0; + +vnle_order=[vnle_order1,vnle_order2,vnle_order3]; +dfe_order = [0 0 0]; + +alpha = 0; + +len_tr = 4096*2; + +mu_ffe1 = 0.0001; +mu_ffe2 = 0.0008; +mu_ffe3 = 0.001; +mu_dc = 0.005; +% mu_dc = 0; + +mu_ffe = [mu_ffe1 mu_ffe3 mu_ffe3]; +mu_dfe = 0.0004; + +dfe_ = sum(dfe_order)>0; + +doub_mode = db_mode.no_db; +cols = linspecer(6); +rop = [-8]; +bwl = [0.5:0.1:1.5]; +fsym = [208:16:256].*1e9; +nonlin_mod = [0.5:0.01:0.75]; +fsym = ones(size(nonlin_mod)).*fsym(1); + +ffe_results = {}; +mlse_results_lin= {}; +vnle_results= {}; +mlse_results_nonlin= {}; +mlse_results_nonlin_states= {}; + +g_eye = GifWriter('Name','eye','Parallel',true); +g_mod = GifWriter('Name','modulator','Parallel',true); + +for r = 1:length(nonlin_mod) + + Pform = Pulseformer("fsym",fsym(r),"fdac",4*fsym(r),"pulse","rc","pulselength",16,"alpha",rcalpha); + + db_precode = 0; + db_encode = 0; + duob_mode = db_mode.no_db; + apply_pulsef = 1; + + [Digi_sig,Symbols,Tx_bits] = PAMsource(... + "fsym",fsym(r),"M",M,"order",19,"useprbs",0,... + "fs_out",fdac,... + "applyclipping",0,"clipfactor",1.5,... + "applypulseform",apply_pulsef,"pulseformer",Pform,... + "randkey",random_key,... + "db_precode",db_precode,"db_encode",db_encode,... + "mrds_code",0,"mrds_blocklength",512,"duobinary_mode",duob_mode).process(); + + % El_sig = AWG("fdac",fdac,"f_cutoff",fsym(r),"lpf_active",0,"kover",kover,"bit_resolution",12,"upsampling_method","samplehold","precomp_sinc_rolloff",0).process(Digi_sig); + El_sig = M8199B("kover",kover).process(Digi_sig); + % AWG("fdac",fdac,"f_cutoff",fsym(r),"lpf_active",0,"kover",kover,"bit_resolution",12,"upsampling_method","samplehold","precomp_sinc_rolloff",0).process(Digi_sig); + + %%%%% Low-pass el. components %%%%%% + % tx_bwl = 100e9; + % El_sig = Filter('filtdegree',3,"f_cutoff",tx_bwl,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(El_sig); + + %%%%% Electrical Driver Amplifier %%%%%% + El_sig = El_sig.normalize("mode","oneone"); + % El_sig = El_sig.setPower(1,"dBm"); + % figure;histogram(El_sig.signal); + + %%%%% MODULATE E/O CONVERSION %%%%% + u_pi = 3.2; + vbias = -u_pi*nonlin_mod(r); + [Opt_sig] = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs,"lambda",laser_wavelength,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",random_key+1).process(El_sig); + + if 1 + figure(15); + hold on + scatter(El_sig.signal(1:100000)+vbias,(abs(Opt_sig.signal(1:100000)).^2)*1e3,0.1,'.','DisplayName','Modulator TF') + xlabel('Input in V') + ylabel('abs(Eopt)2 in mW','Interpreter','latex') + ylim([0 2]); + xlim([-3.2 0]); + g_mod.addFrame(15, r); + + Opt_sig.eye(fsym(r), M, "fignum", 103837); + g_eye.addFrame(103837, r); + end + + %%%%%% Fiber %%%%%% + Opt_sig = Fiber("fsimu",Opt_sig.fs,"fiber_length",link_length/1000,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.07).process(Opt_sig); + + %%%%%% ROP %%%%%% + Opt_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",rop).process(Opt_sig); + + %%%%%% PD Square Law %%%%%% + PD_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20,"nep",1.8e-11,"randomkey",random_key).process(Opt_sig); + + %%%%%% Low-pass RX (PD, El. Connectors and Scope %%%%%% + rx_bwl = 70e9; + PD_sig = Filter('filtdegree',4,"f_cutoff",rx_bwl,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(PD_sig); + + % %%%%%% Low-pass Scope %%%%%% + Lp_scpe = Filter('filtdegree',4,"f_cutoff",110e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true); + + %%%%%% Scope %%%%%% + Scpe_sig = Scope("fsimu",fdac*kover,"fadc",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',Lp_scpe).process(PD_sig); + + Scpe_sig_2sps = Scpe_sig.resample("fs_out",2*fsym(r)); + % Scpe_sig_resampled.signal = Scpe_sig_resampled.signal(1:2*length(Symbols)); + + [~, Scpe_cell, ~, found_sync] = Scpe_sig_2sps.tsynch("reference", Symbols, "fs_ref", fsym(r), "debug_plots", 0); + Rx_sig = Scpe_cell{1}; + Rx_sig = Rx_sig.normalize("mode","rms"); + + if 1 + + %% FFE + % ffe_order = [50, 0, 0]; + % eq_ffe = 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); + % + % ffe_results = ffe(eq_ffe,M,Rx_sig,Symbols,Tx_bits,... + % "precode_mode",duob_mode,... + % 'showAnalysis',0,... + % "postFFE",[],... + % "eth_style_symbol_mapping",0); + % + % ffe_results.metrics.print; + % ffe_results.config.equalizer_structure = "ffe"; + % + + %% MLSE linear + + pf_ncoeffs = 1; + ffe_order = [50, 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",1); + pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); + mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels,'scale_mode',2,'trellis_exclusion',0,'trellis_state_mode',2,'debug',1); + + [ffe_results{r}, mlse_results_lin{r}] = vnle_postfilter_mlse(eq_, pf_, mlse_, M, Rx_sig, Symbols, Tx_bits, ... + "precode_mode", duob_mode,... + 'showAnalysis', 0, ... + "postFFE", [],... + "eth_style_symbol_mapping", 0); + + mlse_results_lin{r}.metrics.print; + ffe_results{r}.metrics.print; + + %% MLSE nonlinear pre + + pf_ncoeffs = 1; + ffe_order = [50, 1, 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",1); + pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); + mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels,'scale_mode',2,'trellis_exclusion',0,'trellis_state_mode',2); + + [vnle_results{r}, mlse_results_nonlin{r}] = vnle_postfilter_mlse(eq_, pf_, mlse_, M, Rx_sig, Symbols, Tx_bits, ... + "precode_mode", duob_mode,... + 'showAnalysis', 0, ... + "postFFE", [],... + "eth_style_symbol_mapping", 0); + + mlse_results_nonlin{r}.metrics.print; + vnle_results{r}.metrics.print; + + %% nonlinear states MLSE linear pre + + pf_ncoeffs = 1; + ffe_order = [50, 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",1); + pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); + mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels,'scale_mode',2,'trellis_exclusion',0,'trellis_state_mode',3); + + [~, mlse_results_nonlin_states{r}] = vnle_postfilter_mlse(eq_, pf_, mlse_, M, Rx_sig, Symbols, Tx_bits, ... + "precode_mode", duob_mode,... + 'showAnalysis', 0, ... + "postFFE", [],... + "eth_style_symbol_mapping", 0); + + mlse_results_nonlin_states{r}.metrics.print; + + end +end + +g_mod.compile(15); +g_eye.compile(103837); + +%% + +figure();hold on; +plot(nonlin_mod,cellfun(@(x) x.metrics.BER, ffe_results),'DisplayName','FFE') +plot(nonlin_mod,cellfun(@(x) x.metrics.BER, vnle_results),'DisplayName','VNLE') +plot(nonlin_mod,cellfun(@(x) x.metrics.BER, mlse_results_lin),'DisplayName','FFE+MLSE') +plot(nonlin_mod,cellfun(@(x) x.metrics.BER, mlse_results_nonlin_states),'DisplayName','FFE+nonlin. states MLSE') +plot(nonlin_mod,cellfun(@(x) x.metrics.BER, mlse_results_nonlin),'DisplayName','VNLE+MLSE') +xlabel('Nonlinear Driving'); +ylabel('BER') +set(gca,'YScale','log'); +legend; +ylim([1e-4 1e-1]); +beautifyBERplot; + +