ML_MLSE (before testing in depth)

minor changes here and there

Classes/04_DSP/Equalizer/ML_MLSE.m

@@ -48,6 +48,11 @@ classdef ML_MLSE < handle
         valid_to_idx
         valid_from_idx
         w
+
+        % --- New: fast state lookup ---
+        state_dict              % containers.Map: key(sequence)->state index
+        key_fmt = '%.8g_';      % key format for sequence strings
+        nSym                      % |constellation|
     end
 
     methods
@@ -79,7 +84,6 @@ classdef ML_MLSE < handle
 
             obj.e = zeros(obj.order,1);
             obj.error = 0;
-
         end
 
         function [X,X_viterbi] = process(obj, X, D)
@@ -88,33 +92,33 @@ classdef ML_MLSE < handle
             % 1 normalize RMS
             X = X.normalize("mode","rms");
 
-            obj.constellation = unique(D.signal);
+            % Use sorted constellation for deterministic mapping
+            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
 
-
             % ==============================================================
             %  INITIALIZATION (only before final epoch and detection mode)
             % ==============================================================
 
-
             % --- Parameters
-            obj.S = numel(unique(D.signal));      % alphabet size
+            obj.S = numel(obj.constellation);      % alphabet size
-            obj.Nf = obj.order*obj.sps;                    % filter length
+            obj.Nf = obj.order*obj.sps;            % filter length
-            % obj.delta = 3;%ceil(obj.Nf/2);         % delay parameter
+            % obj.delta = 3;%ceil(obj.Nf/2);       % delay parameter
             obj.nStates = obj.S^obj.L;
             obj.nFeasible = obj.nStates*obj.S;
 
             % --- Trellis mapping
-            obj.trellis_states = reshape(unique(D.signal),1,[]);
+            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.combs     = fliplr(combvec(pre_comb_cell{:}).');  % rows: states, columns: [x_k, x_{k-1}, ...]
-            obj.first_sym     = obj.combs(:,1);
+            obj.first_sym = obj.combs(:,1);
-            obj.last_sym      = obj.combs(:,end);
+            obj.last_sym  = obj.combs(:,end);
-            obj.nStates       = size(obj.combs,1);
+            obj.nStates   = size(obj.combs,1);
 
             % --- Valid transitions
             obj.valid = false(obj.nStates);
@@ -130,10 +134,17 @@ classdef ML_MLSE < handle
             % --- Allocate vectors and weights
             % !! IF SHAPE FIT, then we already have smth there an we want
             % to start with the existing fitler-set
-            if all(size(obj.w) ~= [obj.Nf+1,obj.nFeasible])
+            if isempty(obj.w) || any(size(obj.w) ~= [obj.Nf+1,obj.nFeasible])
                 obj.w = zeros(obj.Nf+1,obj.nFeasible);     % filter weights per transition + bias tap
+                obj.w = randn(obj.Nf+1,obj.nFeasible);
             end
-            % obj.w = randn(obj.Nf+1,obj.nFeasible);
+
+            % --- Precompute dictionary for fast state lookup (sequence -> state)
+            keys = cell(obj.nStates,1);
+            for i = 1:obj.nStates
+                keys{i} = obj.seq_key(obj.combs(i,:));  % combs row is already [x_k, x_{k-1}, ...]
+            end
+            obj.state_dict = containers.Map(keys, 1:obj.nStates);
 
             % ==============================================================
             %  TRAINING
@@ -165,8 +176,6 @@ classdef ML_MLSE < handle
             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)
@@ -174,24 +183,34 @@ classdef ML_MLSE < handle
             %  FFE + Whitening + ML-Based Branch Metric Estimation + Viterbi
             % ==============================================================
             debug = 0;
+
             % --- Input padding and preallocation
             y = zeros(N,1);
 
+            % number of symbol steps in this block
+            nSymbols = ceil(N/obj.sps);
+
             for epoch = 1:epochs
 
-
+                % state metrics (log-domain costs): keep as column [nStates×1]
-                pm        = zeros(obj.nStates,1);
+                pm        = zeros(obj.nStates,1);  % v_{k-1}(s′)
                 c_hat     = zeros(1,obj.nFeasible);
                 v_tilde   = zeros(1,obj.nFeasible);
-                pred      = zeros(N, obj.nStates, 'uint32');
+                pred      = zeros(nSymbols, obj.nStates, 'uint32');
+                pm_sto    = nan(obj.nStates, nSymbols,'like',pm);
+
+                % --- Initialize "true" trellis state for training (shift-register style)
+                % expect sequences in chronological order [x_{k-L+1}:x_k], but combs rows are [x_k, x_{k-1}, ...]
+                if numel(d) >= obj.L
+                    init_seq = d(1:obj.L);  % [x_1 ... x_L]
+                    true_to_state_idx = obj.state_dict(obj.seq_key(flip(init_seq))); % flip to [x_L, x_{L-1}, ...]
+                else
+                    true_to_state_idx = 1;
+                end
 
-                % ==============================================================
-                %  RUNTIME LOOP
-                % ==============================================================
                 symbol = 0;
                 for sample = 1:obj.sps:N
                     symbol = symbol + 1;
-
                     k = symbol;
 
                     % --- Build Δ-delayed observation window y_k
@@ -201,80 +220,66 @@ classdef ML_MLSE < handle
                     padL = max(0,1 - i1);
                     padR = max(0,i2 - length(x));
                     yk  = [zeros(padL,1); buf(:); zeros(padR,1)];   % Nf×1
-                    yk = [yk;1];
+                    yk  = [yk;1];
 
-                    % --- Predict branch metrics for all feasible transitions
+                    % --- Predict branch metrics for all feasible transitions: c_hat
                     c_hat = (yk.' * obj.w);     % [1×nFeasible]
                     c_hat = c_hat.';            % [nFeasible×1]
 
-                    % --- Extended path metrics
+                    % --- Extended path metrics: v_tilde = pm(from) + c_hat
+                    % normalize pm to avoid growth (invariant to additive const)
+                    pm = pm - min(pm);
                     v_tilde = pm(obj.valid_from_idx) + c_hat;   % [nFeasible×1]
 
                     % ===== Gradient update (Algorithm 1) =====
                     % if training
 
-                    % for current symbol index k -> previous (k-1) and current (k)
                     if k > obj.L
-                        prev_seq = d(k-obj.L:k-1);      % previous state symbols
+                        % shift-register: previous "to" becomes current "from"
-                        true_from_state_idx = find(ismember(obj.combs, flip(prev_seq).', 'rows'));% find state indices in trellis
+                        true_from_state_idx = true_to_state_idx;
 
-                        %if ~(true_from_state_idx == true_to_state_idx), warning('Impossible state transition?!'), end
+                        % current "to" from data window
-
+                        curr_seq = d(k-obj.L+1:k);
-                        curr_seq = d(k-obj.L+1:k);      % next state symbols
+                        key_to   = obj.seq_key(flip(curr_seq));
-                        true_to_state_idx   = find(ismember(obj.combs, flip(curr_seq).', 'rows'));% find state indices in trellis
+                        if isKey(obj.state_dict, key_to)
+                            true_to_state_idx   = obj.state_dict(key_to);
+                        else
+                            % fall back safely (should not happen with proper constellation)
+                            true_to_state_idx   = true_from_state_idx;
+                        end
                     else
                         % not enough history yet
                         true_from_state_idx = 1;
-                        true_to_state_idx   = 1;
+                        true_to_state_idx   = true_to_state_idx; % keep init
                     end
 
                     if 0
                         disp(['FROM: state',char(num2str(true_from_state_idx)),' : symbol transition', char(num2str(obj.combs(true_from_state_idx,:)))]);
                         disp(['TO:   state',char(num2str(true_to_state_idx)),' : symbol transition', char(num2str(obj.combs(true_to_state_idx,:)))]);
                     end
-                    % true_from and true_to are (1) -> (-1)
+
-                    % thus dirac = [1,0,0,0]'
+                    % Dirac delta over correct extended transition (from,to)
                     dirac = zeros(obj.nFeasible,1);
                     dirac(obj.valid_from_idx==true_from_state_idx & obj.valid_to_idx==true_to_state_idx) = 1; % This Dirac delta function δ(s = s∗  k, s′ = s∗  k−1) = 1 if the extended state  (s, s′) corresponds to the true realized states (s∗  k, s∗  k−1), and is  zero otherwise.
-                    
-                    if sum(dirac) == 0
-                        warning('whats happening?!');
-                    end
 
-                    % softmax over -v_tilde, one-hot target t
+                    % softmax over -v_tilde (numerically safe shift)
-                    % 4 feasible transitions: [0,0][0,1][1,0][1,1] #not in
+                    p = exp(-(v_tilde - max(v_tilde)));
-                    % order here!
-                    % first round:
-                    % v_tilde is zero
-                    % -> exp(-0) = 1
-                    % -> 1/(1+1+1+1)
-                    % -> 1/4
-                    % -> each transition is equally likely?
-                    % p = exp(-v_tilde);
-                    p = exp(-(v_tilde-max(v_tilde)));
                     p = p./sum(p); % found in formula (9) and (19)
 
-                    % 1-0.25 = 0.75
+                    % gradient term (t - p)
-                    % 0-0.25 = -0.25
+                    dmp = (dirac - p)';  % 1×nFeasible
-                    % what happens here at the zero? Is this some log
-                    % probability - larger zero is likely; lower zero is
-                    % unlikely? But this is based on known information (training)
-                    dmp = (dirac - p)';
-
 
                     if mod(symbol,128) == 1 && debug
-
                         % --- Normalize and compute probabilities
-                        v_norm = v_tilde - max(v_tilde);          % numerical stability
+                        v_norm = v_tilde - max(v_tilde);
-                        probs_lin = exp(-v_norm);
+                        probs_lin = exp(-v_norm); probs_lin = probs_lin ./ sum(probs_lin);
-                        probs_lin = probs_lin ./ sum(probs_lin);
                         probs_log = -v_norm;
 
                         % --- Map back into nStates×nStates grid
-                        probs_mat = nan(obj.nStates, obj.nStates);
+                        probs_mat    = nan(obj.nStates, obj.nStates);
-                        probs_mat(obj.valid) = probs_lin;         % linear-space probabilities
                         probs_logmat = nan(obj.nStates, obj.nStates);
-                        probs_logmat(obj.valid) = probs_log;      % log-domain scores
+                        probs_mat(obj.valid)    = probs_lin;
+                        probs_logmat(obj.valid) = probs_log;
 
                         % --- Identify the current true transition
                         [to_idx, from_idx] = find(obj.valid);
@@ -282,51 +287,37 @@ classdef ML_MLSE < handle
                         cur_to   = to_idx(cur_idx);
                         cur_from = from_idx(cur_idx);
 
-                        % --- Plot using imagesc (supports hold)
                         figure(11); clf;
-                        imagesc(probs_logmat);
+                        imagesc(probs_logmat); axis xy; colorbar;
-                        axis xy;                                 % origin top-left
+                        xlabel('From state'); ylabel('To state'); set(gca,'FontSize',10);
-                        % colormap(parula);
-                        colorbar;
-                        xlabel('From state');
-                        ylabel('To state');
-                        set(gca,'FontSize',10);
-
-                        % --- Overlay current transition
                         hold on;
-                        plot(cur_from, cur_to, 'rs', ...
+                        plot(cur_from, cur_to, 'rs', 'MarkerSize', 10, 'LineWidth', 2, 'MarkerFaceColor', 'none');
-                            'MarkerSize', 10, 'LineWidth', 2, 'MarkerFaceColor', 'none');
                         hold off;
                     end
 
                     if 1 %training
-                        % the correct state gets a high update, we weight this
+                        % dmp is large for the correct transition -> update emphasizes that branch
-                        % with the input signal, from here the signals are not
-                        % "understandable"
-                        % dmp is large for the correct transition to update
-                        % only this one!
                         dL_Dw = dmp .* (yk); % ∂CE/∂(w) - formula (10)
-                        % from paper: We have observed in simulations that ignoring the derivative of vk−1(s′) during training yields negligible loss after convergence.
+                        % only start with updates when we are inside the signal
-                        % my note: so the update direction is the same for b and w?
-    
-                        %only start with updates when we are inside the signal
                         if k > obj.L
-                            % actual filter training updates:
+                            obj.w = obj.w - mu * dL_Dw; % (Nf+1)×nFeasible
-                            obj.w = obj.w - mu * dL_Dw;% Nf×nFeasible
                         end
-
                     end
 
-                    % compare select
+                    % --- Compare-Select (matrix form, min of costs)
                     v_tilde_mat = inf(obj.nStates, obj.nStates);
                     v_tilde_mat(obj.valid) = v_tilde;
-                    [pm, pred(k,:)] = min(v_tilde_mat, [], 2);
+                    [pm_next, pred(k,:)] = min(v_tilde_mat, [], 2);
 
+                    % re-center to keep metrics bounded (decision-invariant)
+                    pm_next = pm_next - min(pm_next);
+
+                    pm = pm_next;
                     pm_sto(:,symbol) = pm;
-
                 end
 
-                [~, s_end] = max(pm);
+                % --- Traceback (full; you can window with traceback_depth if desired)
+                [~, s_end] = min(pm);
                 viterbi_path = zeros(symbol,1,'uint32');
                 viterbi_path(symbol) = s_end;
                 for n = symbol:-1:2
@@ -334,12 +325,8 @@ classdef ML_MLSE < handle
                 end
 
                 y_vit = obj.first_sym(viterbi_path);
-                y = obj.first_sym(viterbi_path);
+                y     = obj.first_sym(viterbi_path);
 
-
-
-
-                
                 if 1 %debug || training
                     err = sum(y ~= d(1:length(y)));
                     ser = err./length(y);
@@ -357,12 +344,20 @@ classdef ML_MLSE < handle
                     title('Path Metrics (v_tilde)')
 
                     subplot(2,2,4);
-                    % scatter(1:N,pm_sto,1,'.')
+                    scatter(1:symbol,pm_sto,1,'.')
-                    plot(1:symbol,pm_sto)
+                    % plot(1:symbol,pm_sto,'LineStyle','none')
                     title('Path Metric Winners')
                 end
-
             end
         end
     end
+
+    methods (Access=private)
+        function k = seq_key(obj, seq)
+            % Build a stable key string for a sequence row vector in the *same order as combs rows* ([x_k, x_{k-1}, ...])
+            % Use rounding via sprintf to avoid floating-point issues.
+            % seq must be a row vector.
+            k = sprintf(obj.key_fmt, seq);
+        end
+    end
 end

Functions/beautifyBERplot.m

@@ -1,4 +1,7 @@
-function beautifyBERplot()
+function beautifyBERplot(options)
+arguments
+    options.logscale = 1
+end
     % BEAUTIFYBERPLOT Enhances a BER plot for publication-quality figures.
 
     % Set line properties for all current plot lines
@@ -26,8 +29,10 @@ function beautifyBERplot()
     
     % Set logarithmic scale for y-axis, but only if it makes sense.
     % If this is not always desired, you could condition this on the presence of lines or data.
-    set(gca, 'YScale', 'log');
+    if options.logscale
-
+        set(gca, 'YScale', 'log');
+    end
+    
     % Customize grid and box appearance
     set(gca, 'Box', 'on', 'LineWidth', 0.8);  % Thicker border
     grid on;

projects/HighSpeedExperiment_2024/Auswertung_JLT/run_plot_measurements.m

@@ -5,7 +5,7 @@ db = DBHandler("dataBase", [dataBase], "type", database_type);
 M = 4;
 fp = QueryFilter();
 % fp.where('Runs', 'run_id','EQUALS', 987);
-fp.where('Runs', 'pam_level','EQUALS', M);
+% fp.where('Runs', 'pam_level','EQUALS', M);
 % fp.where('Runs', 'symbolrate','EQUALS', 150e9);
 % fp.where('Runs', 'fiber_length','EQUALS', 10);
 fp.where('Runs', 'is_mpi','EQUALS', 0);
@@ -20,10 +20,11 @@ fields =  db.getTableFieldNames('power_state_info');
 fields = [fields; db.getTableFieldNames('dashboard_ungrouped_new')];
 [dataTable,~] = db.queryDB(fp, fields);
 
+%% 
 cfg = struct;
-cfg.x_axis = 'accumulated_dispersion';          % 'symbolrate' | 'baudrate' | 'bitrate' | 'wavelength'
+cfg.x_axis = 'wavelength';          % 'symbolrate' | 'baudrate' | 'bitrate' | 'wavelength'
-cfg.y_axis = 'BER';              % 'BER' | 'GMI' | 'AIR' | ...
+cfg.y_axis = 'power_mzm';              % 'BER' | 'GMI' | 'AIR' | ...
-cfg.group_by = {'equalizer_structure','pre_emph'};
+cfg.group_by = {};
 cfg.filters = struct('is_mpi',0,'pam_level',M,'equalizer_structure',[equalizer_structure.vnle,equalizer_structure.ffe,equalizer_structure.vnle_pf_mlse,equalizer_structure.vnle_db_mlse,equalizer_structure.dfe]);
 
 cfg.y_scale = 'auto';            % auto -> log for BER*, linear otherwise

projects/IMDD_base_system/model_linewidth_evaluation.m

@@ -0,0 +1,181 @@
+%%% 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    = 1300;
+laser_linewidth     = logspace(0,6.2,24);
+
+% Channel
+link_length     = 10;
+
+alpha = 0;
+
+doub_mode = db_mode.no_db;
+cols = linspecer(6);
+rop = [-6];
+bwl = [0.5:0.1:1.5];
+fsym = [200:16:256].*1e9;
+% nonlin_mod = [0.5:0.01:0.75];
+fsym = ones(size(laser_linewidth)).*fsym(1);
+nonlin_mod = ones(size(laser_linewidth)).*0.5;
+
+ffe_results = {};
+mlse_results_lin= {};
+
+parfor r = 1:length(laser_linewidth)
+
+    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",18,"useprbs",1,...
+        "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(r),"randomkey",random_key+1).process(El_sig);
+
+    %%%%%% Fiber %%%%%%
+    mpi = 1;
+    if mpi
+        Combined_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);
+    else
+
+        % 2) ping pong fiber propagation
+        mpi_path = 00;
+        Interference_sig    = Fiber("fsimu",Opt_sig.fs,"fiber_length",mpi_path*2,"alpha",0,"D",0,"lambda0",1310,"gamma",0).process(Opt_sig);
+
+        Interference_sig    = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",-30).process(Interference_sig);
+
+        [Main_sig,dly] = Opt_sig.delay("delay_meter",mpi_path*2);
+
+        % Add
+        Combined_sig = Main_sig + Interference_sig;
+
+        % Cut (due to the delays there is a jump in the signals)
+        if dly == 0;dly = 1;end
+        Combined_sig.signal = Combined_sig.signal(ceil(dly):end);
+
+        % Fiber
+        Combined_sig    = Fiber("fsimu",Combined_sig.fs,"fiber_length",2,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.08).process(Combined_sig);
+    end
+
+
+    %%%%%% ROP %%%%%%
+    Opt_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",rop).process(Combined_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));
+    % Symbols.signal = Symbols.signal(1:Scpe_sig_2sps.length/2);
+
+    % 2sps
+    [~, Scpe_cell, ~, found_sync] = Scpe_sig_2sps.tsynch("reference", Symbols, "fs_ref", fsym(r), "debug_plots", 0);
+    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", 0);
+    Rx_sig_1sps = Scpe_cell_1sps{1};
+    Rx_sig_1sps = Rx_sig_1sps.normalize("mode","rms");
+
+
+    %% RUN DSP
+    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, 3, 3];
+    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",2,"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_2sps, Symbols, Tx_bits, ...
+        "precode_mode", duob_mode,...
+        'showAnalysis', 0, ...
+        "postFFE", [],...
+        "eth_style_symbol_mapping", 0);
+
+
+    ffe_results{r}.metrics.print;
+    mlse_results_lin{r}.metrics.print;
+
+end
+
+
+figure(1);hold on;
+plot(laser_linewidth.*1e-6,cellfun(@(x) x.metrics.BER, ffe_results),'DisplayName','VNLE')
+plot(laser_linewidth.*1e-6,cellfun(@(x) x.metrics.BER, mlse_results_lin),'DisplayName','VNLE+MLSE')
+xlabel('Linewidth [GHz]');
+ylabel('BER')
+set(gca,'YScale','log');
+legend;
+ylim([1e-5 1e-1]);
+beautifyBERplot;
+
+figure();hold on;
+plot(laser_linewidth.*1e-6,cellfun(@(x) x.metrics.AIR.*1e-9, ffe_results),'DisplayName','FFE')
+plot(laser_linewidth.*1e-6,cellfun(@(x) x.metrics.AIR.*1e-9, mlse_results_lin),'DisplayName','FFE+MLSE')
+xlabel('Linewidth [GHz]');
+ylabel('AIR [GBd]')
+% set(gca,'YScale','log');
+legend;
+beautifyBERplot("logscale",0); 
+
+figure(); hold on
+stem(calcWavelengthPlan(16,400e9,1310),ones(16,1),'DisplayName','16x400','Marker','.','LineWidth',1);
+stem(calcWavelengthPlan(8,800e9,1310),ones(8,1),'DisplayName','8x800','Marker','.','LineWidth',1);
+stem(calcWavelengthPlan(16,800e9,1310),ones(16,1),'DisplayName','16x800','Marker','.','LineWidth',1);
+ylim([0,1.2]);
+ylabel('wavelength [nm]');
+xlim([1270, 1350])

projects/IMDD_base_system/simulation_bwl.m

@@ -332,12 +332,14 @@ for r = 1:length(fsym)
         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)
+        for a = 1:numel(alpha_vec)
 
+            alpha_vec(a) = 0.9;
             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);
 
+            showLevelHistogram(eq_signal_prtgt,Symbols_filt,"displayname",'VNLE Out','fignum',201);
 
             if plot_stuff
                 % Plot the response for respective EQ targets

projects/ML_based_MLSE/minimal_model_gpt.m

@@ -0,0 +1,213 @@
+function minimal_model_gpt
+% Minimal binary IM/DD + ISI (M=3), AWGN. RX: ML-based MLSE vs classic MLSE
+% Fixes:
+%  1) Correct (s',s) labeling: s'=[x_{k-1},x_{k-2}], s=[x_k,x_{k-1}]
+%  2) Window standardization (z-score) from pilot
+%  3) Non-negative learned BM: c_hat = (w^T y_std + b).^2
+%  4) Stable training loop bounds; BER length alignment
+
+clear; clc; rng(1);
+
+%% Parameters
+Ktr   = 4000;                 % pilot
+Kte   = 20000;                % test
+M     = 3;                    % true channel memory
+L     = 3;                    % MLSE memory (states keep L-1 symbols)
+h     = [0.55 1.00 0.45];     % ISI FIR
+EsN0dB= 12;
+D     = 5*L;                  % traceback (unused here; full TB)
+Nwin  = 7;                    % BM feature window
+Delta = 0;
+
+%% TX, channel, noise
+map = @(b) 2*b-1;
+b_tr = randi([0 1],Ktr,1);  x_tr = map(b_tr);
+b_te = randi([0 1],Kte,1);  x_te = map(b_te);
+
+convpad = @(x) filter(h,1,[x; zeros(M-1,1)]);
+ytr_clean = convpad(x_tr);
+yte_clean = convpad(x_te);
+
+EsN0   = 10^(EsN0dB/10);
+sigma2 = 1/(2*EsN0);
+awgn   = @(len) sqrt(sigma2)*randn(len,1);
+y_tr   = ytr_clean + awgn(length(ytr_clean));
+y_te   = yte_clean + awgn(length(yte_clean));
+
+% remove tail to keep equal lengths
+trim = M-1;
+y_tr = y_tr(1:end-trim); x_tr = x_tr(1:end-trim); b_tr = b_tr(1:end-trim);
+y_te = y_te(1:end-trim); x_te = x_te(1:end-trim); b_te = b_te(1:end-trim);
+
+%% Trellis
+nStates   = 2^(L-1);
+statesBin = de2bi(0:nStates-1,L-1,'left-msb');
+symOfBit  = @(bit) (2*bit-1);
+
+% transitions (s'->s) with input symbol u in {±1}
+trans = struct('sp',[],'s',[],'u',[],'ubit',[]);
+idx=1;
+for sp=1:nStates
+    prev_bits = statesBin(sp,:); % [x_{k-1}, x_{k-2}] as bits
+    for ubit=[0 1]
+        u  = symOfBit(ubit);
+        nb = [prev_bits(1),ubit];   % new state s = [x_k, x_{k-1}]
+        s  = bi2de(nb,'left-msb')+1;
+        trans(idx).sp=sp; trans(idx).s=s; trans(idx).u=u; trans(idx).ubit=ubit;
+        idx=idx+1;
+    end
+end
+nTrans = numel(trans);
+
+%% Classic BM (Gaussian)
+BM_classic = @(yk, sp, s, u) (( yk - h * [u, symOfBit(statesBin(sp,:))]')^2)/(2*sigma2);
+
+%% ML-based BM: c_hat = (w^T y_std + b)^2  (nonnegative)
+padN   = Nwin-1;
+pad    = @(y) [zeros(Delta+padN,1); y; zeros(max(0,Nwin-1-Delta),1)];
+y_tr_p = pad(y_tr); y_te_p = pad(y_te);
+
+% standardize window features using pilot
+Ytr_mat = im2col_sliding(y_tr_p, Nwin, Delta); % Nwin × T
+mu_win  = mean(Ytr_mat,2);
+sd_win  = std(Ytr_mat,0,2)+1e-8;
+
+W = zeros(Nwin,nTrans);
+b0= zeros(1,nTrans);
+softmax = @(z) exp(z - max(z))./sum(exp(z - max(z)));
+mu = 0.005;           % LR
+nEpoch = 6;           % light training
+Ttr = length(y_tr);
+
+for ep=1:nEpoch
+    v = zeros(nStates,1); v(2:end)=Inf;
+    for k = L : min(Ttr - (L-1), Ttr)   % need x(k-1), x(k-2)
+        % true (s',s):
+        sp_bits = [(x_tr(k-1)<0), (x_tr(k-2)<0)];
+        s_bits  = [(x_tr(k)  <0), (x_tr(k-1)<0)];
+        sp_star = bi2de(sp_bits,'left-msb')+1;
+        s_star  = bi2de(s_bits ,'left-msb')+1;
+
+        % feature window (z-scored)
+        yw = y_tr_p(k-Delta : k-Delta+Nwin-1);
+        yw = (yw - mu_win) ./ sd_win;
+
+        % compute extended PM logits over all (s',s)
+        vtil = -inf(nTrans,1);
+        for t=1:nTrans
+            z = W(:,t).'*yw + b0(t);
+            c_hat = z*z;                         % (.)^2  nonnegative BM
+            vtil(t) = -( v(trans(t).sp) + c_hat);
+        end
+        pext = softmax(vtil).';
+        t_star = find([trans.sp]==sp_star & [trans.s]==s_star,1);
+
+        % CE gradient wrt c_hat, chain rule through square
+        delta = pext; delta(t_star)=delta(t_star)-1; delta = -delta; % sign for vtil=-(...)
+        for t=1:nTrans
+            z = W(:,t).'*yw + b0(t);
+            dc_dz = 2*z;                 % d(z^2)/dz
+            g = delta(t) * dc_dz;        % ∂L/∂z
+            W(:,t) = W(:,t) - mu * g * yw;
+            b0(t)  = b0(t)  - mu * g;
+        end
+
+        % PM update for next step
+        v_next = inf(nStates,1);
+        for t=1:nTrans
+            sp=trans(t).sp; s=trans(t).s;
+            z = W(:,t).'*yw + b0(t);
+            c_hat = z*z;
+            cand = v(sp)+c_hat;
+            if cand < v_next(s)
+                v_next(s)=cand;
+            end
+        end
+        v = v_next;
+    end
+end
+
+%% Decode (full-length traceback)
+[xhat_ml, bhat_ml] = viterbi_decode(y_te, L, trans, @(k) feat_std(y_te_p,k,Nwin,Delta,mu_win,sd_win), ...
+                                    @(yk,sp,s,u) learnedBM(W,b0,yk,sp,s,u), 0);
+[xhat_ex, bhat_ex] = viterbi_decode(y_te, L, trans, @(k) y_te(k), ...
+                                    @(yk,sp,s,u) BM_classic(yk,sp,s,u), 0);
+
+%% BER (align by min length)
+Lmin = min([length(bhat_ml), length(bhat_ex), length(b_te)]);
+BER_ml = mean(bhat_ml(1:Lmin) ~= (b_te(1:Lmin)>0));
+BER_ex = mean(bhat_ex(1:Lmin) ~= (b_te(1:Lmin)>0));
+fprintf('BER (ML-based MLSE):   %.3e\n', BER_ml);
+fprintf('BER (classic MLSE  ):  %.3e\n', BER_ex);
+end
+
+% ----------------- helpers -----------------
+function Y = im2col_sliding(y_pad, Nwin, Delta)
+T = length(y_pad) - (Nwin-1) - Delta;
+Y = zeros(Nwin,T);
+for k=1:T
+    Y(:,k) = y_pad(k-Delta : k-Delta+Nwin-1);
+end
+end
+
+function yw = feat_std(y_pad,k,Nwin,Delta,mu_win,sd_win)
+yw = y_pad(k-Delta : k-Delta+Nwin-1);
+yw = (yw - mu_win) ./ sd_win;
+end
+
+function c = learnedBM(W,b0,yw,sp,s,~)
+% pick parameters of the (sp->s) transition
+persistent map; 
+if isempty(map)
+    % build once: index of W/b0 for each (sp,s)
+    nStates = size(W,1)*0+1; %#ok<NASGU>
+end
+% linear scan is fine at this size:
+% (use first match of (sp,s))
+c = inf;
+for t=1:size(W,2)
+    % suppose we stored (sp,s) order as in training; we cannot access here.
+    % Instead, pass the exact column via a small mapper (build each call):
+end
+% Faster: precompute a table outside; here we reconstruct like in training
+% (rebuild tiny mapper)
+persistent key_sp key_s
+if isempty(key_sp)
+    key_sp = evalin('caller','[trans.sp]');
+    key_s  = evalin('caller','[trans.s]');
+end
+t = find(key_sp==sp & key_s==s,1);
+z = W(:,t).'*yw + b0(t);
+c = z*z;
+end
+
+function [xhat, bhat] = viterbi_decode(y, L, trans, getFeat, BMfun, D_unused)
+nStates = 2^(L-1);
+T = length(y);
+v = inf(nStates,T); v(:,1)=Inf; v(1,1)=0;
+prev = zeros(nStates,T); in_u = zeros(nStates,T);
+
+for k=1:T
+    if k==1, vprev=inf(nStates,1); vprev(1)=0; else, vprev=v(:,k-1); end
+    vcur = inf(nStates,1); pre=zeros(nStates,1); inb=zeros(nStates,1);
+
+    feat = getFeat(k);   % either scalar y(k) or standardized window
+    for t=1:numel(trans)
+        sp=trans(t).sp; s=trans(t).s; u=trans(t).u;
+        ck = BMfun(feat, sp, s, u);
+        cand = vprev(sp)+ck;
+        if cand < vcur(s)
+            vcur(s)=cand; pre(s)=sp; inb(s)=u;
+        end
+    end
+    v(:,k)=vcur; prev(:,k)=pre; in_u(:,k)=inb;
+end
+
+[~,st]=min(v(:,T));
+xhat=zeros(T,1);
+for k=T:-1:1
+    xhat(k)=in_u(st,k);
+    st=prev(st,k); if st==0, st=1; end
+end
+bhat = xhat>0;
+end

projects/ML_based_MLSE/model.m

@@ -22,9 +22,9 @@ alpha = 0;
 
 doub_mode = db_mode.no_db;
 cols = linspecer(6);
-rop = [-6];
+rop = [-5];
 bwl = [0.5:0.1:1.5];
-fsym = [208:16:256].*1e9;
+fsym = [212:16:256].*1e9;
 % nonlin_mod = [0.5:0.01:0.75];
 nonlin_mod = ones(size(fsym)).*0.5;
 
@@ -41,7 +41,7 @@ for r = 1:length(fsym)
     apply_pulsef = 1;
 
     [Digi_sig,Symbols,Tx_bits] = PAMsource(...
-        "fsym",fsym(r),"M",M,"order",17,"useprbs",0,...
+        "fsym",fsym(r),"M",M,"order",18,"useprbs",0,...
         "fs_out",fdac,...
         "applyclipping",0,"clipfactor",1.5,...
         "applypulseform",apply_pulsef,"pulseformer",Pform,...
@@ -94,59 +94,87 @@ for r = 1:length(fsym)
     [~, 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");
+    showLevelHistogram(Rx_sig_1sps,Symbols,"displayname",'ffe','fignum',111);
 
 
-    %% Implement DSP directly here:
+    %%
 
     mu_lms = 0.0005;
-    % tic
+    eq_ = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",2^13,"mu_dd",mu_lms,"mu_tr",mu_lms,"order",50,"sps",2,"dd_mode",1,"adaption_technique","lms");
-    % 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",1);
+    pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1);
-    % [Eq_signal,Vit_signal] = eq.process(Rx_sig_2sps,Symbols);
+    
-    % toc
+    % FFE
+    [y_ffe, ffe_noise] = eq_.process(Rx_sig_2sps, Symbols);
+    
+    Eq_bits = PAMmapper(M, 0, "eth_style", 0).demap(y_ffe);
+    [~, errors, ber_ffe, ~] = calc_ber(Eq_bits.signal, Tx_bits.signal, "skip_front", 0, "skip_end", 0, "returnErrorLocation", 1);
+    fprintf('FFE: %.2e \n',ber_ffe);
 
-    tic
+    % Postfilter
-    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);
+    [y_white,whitened_noise] = pf_.process(y_ffe, ffe_noise);
-    [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;
+    % Sequence Est
-    Vit_signal_.signal = circshift(Vit_signal.signal,0);
+    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,'DIR',pf_.coefficients);
-    Vit_bits = PAMmapper(M, 0, "eth_style", 0).demap(Vit_signal_);
+    [y_mlse] = mlse_.process(y_white,Symbols);
-    [~, errors, ber, errpos] = calc_ber(Vit_bits.signal, Tx_bits.signal, "skip_front", 0, "skip_end", 0, "returnErrorLocation", 1);
+
-    fprintf('Viterbi: %.2e \n',ber);
+    Vit_bits = PAMmapper(M, 0, "eth_style", 0).demap(y_mlse);
+    [~, errors, ber_mlse_normal, errpos] = calc_ber(Vit_bits.signal, Tx_bits.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1);
+    fprintf('MLSE: %.2e \n',ber_mlse_normal);
+
+    showLevelHistogram(y_ffe,Symbols,"displayname",'ffe','fignum',111);
+    % showLevelHistogram(y_white,Symbols,"displayname",'ffe','fignum',111);
+
+    %% optimize length
+
+    mu_lms = 0.15;
+    eq = ML_MLSE("epochs_tr",5,"epochs_dd",5,"len_tr",2^14,...
+        "mu_dd",mu_lms,"mu_tr",mu_lms,"order",5,"sps",2,...
+        "traceback_depth",128,"L",2,"delta",0);
+
+    [y_ml_mlse,Vit_signal] = eq.process(Rx_sig_2sps,Symbols);
+    y_ml_mlse_ = y_ml_mlse;
+    y_ml_mlse_.signal = circshift(y_ml_mlse.signal,0);
+    Vit_bits = PAMmapper(M, 0, "eth_style", 0).demap(y_ml_mlse_);
+    [~, errors, ber, errpos] = calc_ber(Vit_bits.signal, Tx_bits.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1);
+    fprintf('ML MLSE: %.2e \n',ber);
+
+    bursts = count_error_bursts(errpos, 10);
+    e = zeros(size(Vit_bits.signal));
+    e(errpos) = 1;
+    figure(8)
+    stem(e)
+
+    %% optimize delta
+        deltas = [-1:4];
+        ber = zeros(1,length(deltas));
+        parfor m = 1:numel(deltas)
+
+            mu_lms = 0.2;
+            eq = ML_MLSE("epochs_tr",2,"epochs_dd",5,"len_tr",2^13,...
+                        "mu_dd",mu_lms,"mu_tr",mu_lms,"order",4,"sps",1,...
+                        "traceback_depth",128,"L",3,"delta",deltas(m));
+    
+            [y_ml_mlse,Vit_signal] = eq.process(y_ffe,Symbols);
+            y_ml_mlse_ = y_ml_mlse;
+            y_ml_mlse_.signal = circshift(y_ml_mlse.signal,0);
+            Vit_bits = PAMmapper(M, 0, "eth_style", 0).demap(y_ml_mlse_);
+            [~, errors, ber(m), errpos] = calc_ber(Vit_bits.signal, Tx_bits.signal, "skip_front", 0, "skip_end", 0, "returnErrorLocation", 1);
+            fprintf('ML MLSE: %.2e \n',ber(m));
+        end
 
 
+        figure(7); hold on
+        title('ML MLSE')
+        plot(deltas,ber,'DisplayName','BER');
+        yline(ber_ffe,'DisplayName','BER FFE');
+        yline(ber_mlse_normal,'DisplayName','BER MLSE');
+        xlabel('deltas')
+        beautifyBERplot
+        legend
+        ylim([1e-5, 1e-1]);
+        set(gca,'YScale','log');
 
-
+   
-    %% optimize smth.
-    tr_len = 2.^[2:15];
-    tr_len = floor(tr_len);
-    ber = zeros(size(tr_len));
-    for 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
@@ -163,12 +191,12 @@ for r = 1:length(fsym)
     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_ = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",2^13,"mu_dd",mu_lms,"mu_tr",mu_lms,"order",50,"sps",2,"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, ...
+    [ffe_results{r}, mlse_results_lin{r}] = vnle_postfilter_mlse(eq_, pf_, mlse_, M, Rx_sig_2sps, Symbols, Tx_bits, ...
         "precode_mode", duob_mode,...
         'showAnalysis', 0, ...
         "postFFE", [],...