From 16846d4eb75f9543ae7964caee7d69c69abb023d Mon Sep 17 00:00:00 2001
From: Silas Oettinghaus <silas.oettinghaus@email.uni-kiel.de>
Date: Mon, 12 Jan 2026 11:36:15 +0100
Subject: [PATCH 1/7] stuff

---
 .../a_minimal_example.m                       | 125 +++++++++---------
 projects/WDM/WDM_auswertung.m                 |  14 +-
 projects/WDM/WDM_model_16x200.m               |  30 +++--
 3 files changed, 91 insertions(+), 78 deletions(-)

diff --git a/projects/HighSpeedExperiment_2024/a_minimal_example.m b/projects/HighSpeedExperiment_2024/a_minimal_example.m
index 966d5b9..2a94e13 100644
--- a/projects/HighSpeedExperiment_2024/a_minimal_example.m
+++ b/projects/HighSpeedExperiment_2024/a_minimal_example.m
@@ -4,70 +4,77 @@ dsp_options.storage_path    = 'Z:\2024\sioe_labor\';
 dsp_options.max_occurences = 1;
 db    = DBHandler("dataBase", 'labor_highspeed', "type",  'mysql' );
 
-fp = QueryFilter();
-
-fp.where('Runs','fiber_length','EQUALS', 2);
-fp.where('Runs','wavelength','EQUALS', 1310);
-fp.where('Runs','bitrate','EQUALS', 300e9);
-fp.where('Runs','pam_level','EQUALS', 4);
-fp.where('Runs','rop_attenuation','EQUALS', 0);
-fp.where('Runs','is_mpi','EQUALS', 0);
-fp.where('Runs', 'db_mode','EQUALS', 0); 
-% fields = db.getTableFieldNames('Runs');
-% [dataTable,~] = db.queryDB(fp, fields);
-fields = db.getTableFieldNames('power_state_info');
-fields = [fields; db.getTableFieldNames('Runs')];
-fields = [fields; db.getTableFieldNames('dashboard_ungrouped_alltime')]; %dashboard_ungrouped_after_nov_2025  dashboard_ungrouped_aug_nov_2025 
-fields = unique(fields);
-[dataTable,~] = db.queryDB(fp, fields);
-
-
-fsym        = dataTable(1,:).symbolrate;
-M           = double(dataTable(1,:).pam_level);
-duob_mode   = db_mode(strrep(dataTable(1,:).db_mode,'"',''));
-
-
-% Load and Sync signal data from DB
-[Tx_bits, Symbols, Scpe_cell, ~] = loadAndSyncSignalDataFromDb(dataTable(1,:), dsp_options);
-
-% Preprocess signal
-Scpe_sig = preprocessSignal(Scpe_cell{1}, Symbols, fsym);
-
-% Show spectrum
-Scpe_sig.spectrum("fignum",1,"displayname",'Rx')
-
-% meta = struct();
-% meta.varnames = dataTable.Properties.VariableNames;
+% fp = QueryFilter();
 % 
-% for k = 1:numel(meta.varnames)
-%     v = meta.varnames{k};
-%     col = dataTable.(v);
+% fp.where('Runs','fiber_length','EQUALS', 2);
+% fp.where('Runs','wavelength','EQUALS', 1310);
+% fp.where('Runs','bitrate','EQUALS', 300e9);
+% fp.where('Runs','pam_level','EQUALS', 4);
+% fp.where('Runs','rop_attenuation','EQUALS', 0);
+% fp.where('Runs','is_mpi','EQUALS', 0);
+% fp.where('Runs', 'db_mode','EQUALS', 0); 
+% % fields = db.getTableFieldNames('Runs');
+% % [dataTable,~] = db.queryDB(fp, fields);
 % 
-%     if isnumeric(col) || islogical(col)
-%         meta.(v) = col;
-%     elseif isstring(col)
-%         meta.(v) = cellstr(col);
-%     elseif iscellstr(col)
-%         meta.(v) = col;
-%     else
-%         error("Unsupported table column type: %s", class(col))
-%     end
-% end
-% exp_data.metadata = meta;
 
-exp_data = struct();
-exp_data.metadata = dataTable;
-exp_data.tx_bits = Tx_bits.signal;
-exp_data.tx_signal = Symbols.signal;
-exp_data.rx_signal_2sps = Scpe_sig.signal;
+run_ids = [ 993 1205 1413 1623 1833 2043 2253 2628 2836 2958 3000 3042 3323 5098 5225 5267 5309];
 
-fname = dataTable(1,:).rx_raw_path;
-[~, filename, ext] = fileparts(fname);
-filename = strrep(filename,"_raw_signal","");
-filename = filename + ext;
-savepath = fullfile('F:\2024\sioe_labor\export_skuehl\',filename);
-save(savepath,'exp_data','-v7.3');
+for id = run_ids
+    
+    fp = QueryFilter();
+    fp.where('Runs', 'run_id','EQUALS', id); 
+    fields = db.getTableFieldNames('power_state_info');
+    fields = [fields; db.getTableFieldNames('Runs')];
+    fields = [fields; db.getTableFieldNames('dashboard_ungrouped_alltime')]; %dashboard_ungrouped_after_nov_2025  dashboard_ungrouped_aug_nov_2025 
+    fields = unique(fields);
+    [dataTable,~] = db.queryDB(fp, fields);
+    
+    fsym        = dataTable(1,:).symbolrate;
+    M           = double(dataTable(1,:).pam_level);
+    duob_mode   = db_mode(strrep(dataTable(1,:).db_mode,'"',''));
+    
+    % Load and Sync signal data from DB
+    [Tx_bits, Symbols, Scpe_cell, ~] = loadAndSyncSignalDataFromDb(dataTable(1,:), dsp_options);
+    
+    % Preprocess signal
+    Scpe_sig = preprocessSignal(Scpe_cell{1}, Symbols, fsym);
+    
+    % Show spectrum
+    Scpe_sig.spectrum("fignum",1,"displayname",'Rx')
+    
+    meta = struct();
+    meta.varnames = dataTable.Properties.VariableNames;
+    
+    for k = 1:numel(meta.varnames)
+        v = meta.varnames{k};
+        col = dataTable.(v);
+    
+        if isnumeric(col) || islogical(col)
+            meta.(v) = col;
+        elseif isstring(col)
+            meta.(v) = cellstr(col);
+        elseif iscellstr(col)
+            meta.(v) = col;
+        else
+            error("Unsupported table column type: %s", class(col))
+        end
+    end
+    exp_data.metadata = meta;
+    
+    exp_data = struct();
+    exp_data.metadata = dataTable;
+    exp_data.tx_bits = Tx_bits.signal;
+    exp_data.tx_signal = Symbols.signal;
+    exp_data.rx_signal_2sps = Scpe_sig.signal;
+    
+    fname = dataTable(1,:).rx_raw_path;
+    [~, filename, ext] = fileparts(fname);
+    filename = strrep(filename,"_raw_signal","");
+    filename = filename + ext;
+    savepath = fullfile('F:\2024\sioe_labor\export_skuehl\',filename);
+    save(savepath,'exp_data','-v7.3');
 
+end
 
 
 
diff --git a/projects/WDM/WDM_auswertung.m b/projects/WDM/WDM_auswertung.m
index d13c97f..1f3b949 100644
--- a/projects/WDM/WDM_auswertung.m
+++ b/projects/WDM/WDM_auswertung.m
@@ -21,7 +21,7 @@ Smlse = cell(1,N);
 Sdbt = cell(1,N);
 
 % Choose your quantile band. For your old style, use 0.04/0.99:
-qLow  = 0.0;                  % lower quantile (e.g., 0.04 for old script)
+qLow  = 0.0;                  % lower quantile (e.g.,s 0.04 for old script)
 qHigh = 1;                  % upper quantile (e.g., 0.99 for old script)
 cols = linspecer(N);           % one color per wavelength (Ch)
 cols = cbrewer2('set1',N); 
@@ -51,10 +51,10 @@ for l = 1:N
     showLegend = 1; % one legend entry per technique
 
     % Plot shaded band + mean line with boundedline
-    % plotBandMeanBL(rop, ffe_mat,  cols(l,:), sprintf('FFE @ %.1d nm',round(wavelengthplan(l))),  qLow, qHigh, '--s', showLegend);
+    plotBandMeanBL(rop, ffe_mat,  cols(l,:), sprintf('FFE @ %.1d nm',round(wavelengthplan(l))),  qLow, qHigh, '--s', showLegend);
     % scatter(Sffe,fec.*ones(size(Sffe)),20,'v','MarkerFaceColor','black');
 
-    plotBandMeanBL(rop, vnle_mat, cols(l,:), sprintf('VNLE @ %.1d nm',round(wavelengthplan(l))), qLow, qHigh, '--x', showLegend);
+    % plotBandMeanBL(rop, vnle_mat, cols(l,:), sprintf('VNLE @ %.1d nm',round(wavelengthplan(l))), qLow, qHigh, '--x', showLegend);
 
     % plotBandMeanBL(rop, mlse_mat, cols(l,:), sprintf('VNLE+PF+MLSE @ %.1d nm',round(wavelengthplan(l))), qLow, qHigh, '-o',  showLegend);
 
@@ -73,7 +73,7 @@ ylim([1e-5 0.3]);
 grid on;
 legend show;
 
-
+%%
 S_cell = Sdbt;
 S_cell =Smlse;
 S_cell = {Svnle,Smlse,Sdbt};
@@ -106,7 +106,7 @@ for i = 1:length(S_cell)
         'Bandwidth', 0.05);
     
     ylim([floor(min(S_mat,[],'all')), ceil(max(S_mat,[],'all'))])
-     ylim([-8 0]);
+    % ylim([-8 0]);
     ylabel('ROP at FEC crossing');
     title(sprintf('RROP to cross BER %.2e', fec));
     grid on; box on;
@@ -118,10 +118,6 @@ end
 
 
 
-
-
-
-
 %% ================= helper =================
 function plotBandMeanBL(x, Y, color, techLabel, qLow, qHigh, lineSpec, showLegend)
     % Y: (nPoints x nRealizations)
diff --git a/projects/WDM/WDM_model_16x200.m b/projects/WDM/WDM_model_16x200.m
index 7000966..3d6dac9 100644
--- a/projects/WDM/WDM_model_16x200.m
+++ b/projects/WDM/WDM_model_16x200.m
@@ -6,9 +6,9 @@ run(fullfile(fileparts(mfilename('fullpath')),'WDM_settings.m'));
 s.num_realiz      = 2;
 s.wavelengthplan  = calcWavelengthPlan(16,400e9,1310);
 % wavelengthplan  = [1295,1305,1315,1325];
-link_length     = 0;
-pmd             = 0.1;
-gamma           = 0.0023;
+link_length     = 1;
+s.pmd             = 0.1;
+s.gamma           = 0.0023;
 
 s.M               = 4;
 m = floor(log2(s.M)*10)/10;
@@ -52,7 +52,7 @@ f_plan = physconst('lightspeed')./(s.wavelengthplan.*1e-9);
 margin = 25e12; % some THz left and right
 f_span = (max(f_plan)+margin)-(min(f_plan)-margin);
 f_nyq = f_span/2;
-kover = 8;
+kover = 4;
 upsample_required = f_nyq./(fdac*kover/2);
 upsample_pow = 2^nextpow2(upsample_required);
 upsample_ceil = ceil(upsample_required);
@@ -64,13 +64,23 @@ signal_cell = {};
 Symbols = {};
 Tx_bits = {};
 
-s.rop = -6:0.75:-0.75;
+s.rop = -12:0.75:-0.75;
 
 output_ffe = cell(length(s.wavelengthplan),length(s.rop),s.num_realiz);
 output_vnle = cell(length(s.wavelengthplan),length(s.rop),s.num_realiz);
 output_mlse = cell(length(s.wavelengthplan),length(s.rop),s.num_realiz);
 output_dbt = cell(length(s.wavelengthplan),length(s.rop),s.num_realiz);
 
+s.p = "pair";
+switch s.p
+    case "co"
+        pol_rot = 100.*ones(1,length(s.wavelengthplan));
+    case "pair"
+        pol_rot = repmat([100,100,0,0],1,length(s.wavelengthplan)/4);
+    case "alt"
+        pol_rot = repmat([100,0,100,0],1,length(s.wavelengthplan)/4);
+end
+
 for realiz = 1:s.num_realiz
 
   
@@ -100,7 +110,7 @@ for realiz = 1:s.num_realiz
         %%%%% s.MODULATE E/O CONVERSION %%%%%
         Eml_out = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs,"lambda",s.wavelengthplan(l),"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",s.random_key+l+realiz).process(El_sig);
         
-        signal_cell{l} = Polarization_Controller("mode","rot_power","desired_power",100).process(Eml_out);
+        signal_cell{l} = Polarization_Controller("mode","rot_power","desired_power",pol_rot(l)).process(Eml_out);
     end
     
     Opt_sig_wdm = Optical_Multiplex("fs_in",fdac*kover,"fs_out",upsample_pow*fdac*kover,...
@@ -108,7 +118,7 @@ for realiz = 1:s.num_realiz
     
     Opt_sig_wdm = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",3+10*log10(N)).process(Opt_sig_wdm);
     
-    % Opt_sig_wdm.spectrum("fignum",101,"displayname",'bla','normalizeTo0dB',0,'lambda0_nm',1310,'useWavelengthAxis',0);
+    Opt_sig_wdm.spectrum("fignum",101,"displayname",'bla','normalizeTo0dB',0,'lambda0_nm',1310,'useWavelengthAxis',0);
     
     % Opt_sig_wdm.spectrum("fignum",101,"displayname",'bla','normalizeTo0dB',1,'max_num_lines',2);
     
@@ -123,7 +133,7 @@ for realiz = 1:s.num_realiz
     Dvec = getDispersionVector(nSegments, D_local, zdw, randomize_D, s.random_key+realiz);
     for seg = 1:nSegments
         
-        Opt_sig_wdm_fib = DP_Fiber("L",segment_length,"D",Dvec(s),"Dpmd",pmd,"Ds",0.07,...
+        Opt_sig_wdm_fib = DP_Fiber("L",segment_length,"D",Dvec(seg),"Dpmd",s.pmd,"Ds",0.07,...
                         "beat_len",10,"corr_len",100,"dz",1,"manakov",0,...
                         "gamma",s.gamma,"lambda",zdw,"n_waveplates",10,"SS_dphimax",0.01,...
                         "SS_dzmax",50,"SS_dzmin",10,"X_alpha",0.3,"X_beta",0,"rng",1).process(Opt_sig_wdm_fib);
@@ -135,11 +145,11 @@ for realiz = 1:s.num_realiz
     % Opt_sig_wdm_fib.move_it_spectrum("fignum",100212,"displayname",'bla');
     
     % Opt_sig    = Fiber("fsimu",Opt_sig.fs,"fiber_length",s.link_length/1000,"alpha",0.3,"D",0,"lambda0",1310,"s.gamma",0,"Dslope",0.07).process(Opt_sig)
-    Opt_sig_wdm_demux = Optical_Demultiplex("attenuation",0,"B",200e9,"filtype",1,"fs_out",Opt_sig_wdm_rx.fs/upsample_pow,"fs_in",Opt_sig_wdm_rx.fs,"lambda_center",1310).process(Opt_sig_wdm_fib);
+    Opt_sig_wdm_demux = Optical_Demultiplex("attenuation",0,"B",200e9,"filtype",1,"fs_out",fdac*kover,"fs_in",fdac*kover*upsample_pow,"lambda_center",1310).process(Opt_sig_wdm_fib);
         
     for ri = 1:length(s.rop)
     
-        parfor l = 1:N
+        for l = 1:N
 
             %%%%%% ROP %%%%%%
             Opt_sig_wdm_rx = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",s.rop(ri)).process(Opt_sig_wdm_demux{l}); % rop+10*log10(N)

From 719e5508e776c18e3dcf72a954c71bdedda179e0 Mon Sep 17 00:00:00 2001
From: Silas Oettinghaus <silas.oettinghaus@email.uni-kiel.de>
Date: Tue, 20 Jan 2026 10:49:59 +0100
Subject: [PATCH 2/7] FSO analysis: - minimal changes in several files - new
 Timing Rec from Magnus

---
 Classes/00_signals/Signal.m                   |   6 +-
 Classes/04_DSP/Timing_Recovery.m              |  10 +-
 Functions/EQ_structures/dsp_runid.m           |   7 +-
 Functions/EQ_structures/ffe.m                 |  25 +-
 Functions/EQ_structures/ml_mlse.m             |   6 +-
 .../EQ_structures/vnle_postfilter_mlse.m      |   2 +-
 .../{violinplot.m => violinplot_community.m}  |   0
 projects/FSO_transmission/first_analysis.m    | 215 ++++++++++++------
 .../Auswertung_JLT/final/FIGURE_EYES.m        |  36 ++-
 .../Auswertung_JLT/run_dsp_from_db.m          |   8 +-
 .../baudrate_evaluation.m                     |   2 +-
 projects/IMDD_base_system/minimal_example.m   |  50 ++--
 .../lambda_vs_spectrum.m                      |  88 +++++++
 projects/WDM/WDM_auswertung.m                 | 186 ++++++++++++++-
 14 files changed, 511 insertions(+), 130 deletions(-)
 rename Libs/Violinplot-Matlab-master/{violinplot.m => violinplot_community.m} (100%)
 create mode 100644 projects/Lab_analysis/data_analysis_and_plot_scripts/lambda_vs_spectrum.m

diff --git a/Classes/00_signals/Signal.m b/Classes/00_signals/Signal.m
index 5905556..1377ddc 100644
--- a/Classes/00_signals/Signal.m
+++ b/Classes/00_signals/Signal.m
@@ -971,8 +971,8 @@ classdef Signal
                 maxA = max(sig(100:end-100))*1.3;
                 minA = min(sig(100:end-100))*1.3;
 
-                % maxA = 0.12;
-                % minA = -0.08;
+                maxA = 0.12;
+                minA = -0.08;
 
                 difference= maxA-minA;
 
@@ -993,7 +993,9 @@ classdef Signal
                 % beautify
                 cm=flip(cbrewer2("RdYlBu",4096));
                 % cm=flip(cbrewer2("RdBu",4096));
+                % cm=flip(cbrewer2("Blues",4096));
                 cm(1,:) = [1,1,1]; % set zeros to white => clean background
+
                 colormap(cm);
                 % colormap('turbo');
                 % ax.CLim =  [0   50];
diff --git a/Classes/04_DSP/Timing_Recovery.m b/Classes/04_DSP/Timing_Recovery.m
index ba16f10..23cc331 100644
--- a/Classes/04_DSP/Timing_Recovery.m
+++ b/Classes/04_DSP/Timing_Recovery.m
@@ -9,7 +9,7 @@ classdef Timing_Recovery < handle
     end
 
     methods(Access=public)
-        function obj = FFE(options)
+        function obj = Timing_Recovery(options)
             arguments(Input)
 
                 options.timing_error_detector = 'Gardner';
@@ -25,12 +25,11 @@ classdef Timing_Recovery < handle
                 obj.(fn{n}) = options.(fn{n});
             end
 
-            obj.e = zeros(obj.order,1);
-            obj.error = 0;
+ 
 
         end
 
-        function data_out = process(obj, data_in)
+        function [data_out,timing_error] = process(obj, data_in)
 
             timing_synchronization = comm.SymbolSynchronizer( ...
                 "TimingErrorDetector", obj.timing_error_detector, ...
@@ -39,7 +38,8 @@ classdef Timing_Recovery < handle
                 "NormalizedLoopBandwidth", obj.normalized_loop_bandwidth, ...
                 "DetectorGain",        obj.detector_gain);
 
-            data_out.signal = timing_synchronization(data_in.signal);
+            data_out = data_in;
+            [data_out.signal,timing_error] = timing_synchronization(data_in.signal);
 
         end
     end
diff --git a/Functions/EQ_structures/dsp_runid.m b/Functions/EQ_structures/dsp_runid.m
index 4d207bd..c10cb5c 100644
--- a/Functions/EQ_structures/dsp_runid.m
+++ b/Functions/EQ_structures/dsp_runid.m
@@ -95,7 +95,7 @@ try
     adaption= 1;
     use_dd_mode = 1;
 
-    use_ffe         = 0;
+    use_ffe         = 1;
     use_dfe         = 0;
     use_vnle_mlse   = 0;
     use_dbtgt       = 0;
@@ -149,6 +149,9 @@ try
         % Scpe_sig.spectrum("fignum",22233,"normalizeTo0dB",0,"displayname",'Rx');
         % Scpe_sig.eye(fsym,M,"fignum",1024);
 
+        Scpe_sig.signal = Scpe_sig.signal(1:2*Symbols.length);
+        Scpe_sig.signal = real(Scpe_sig.signal);
+
         if duob_mode ~= db_mode.db_encoded
 
             if use_ffe
@@ -246,7 +249,7 @@ try
                 %ML-based MLSE (L=2)
                 mu_ml = 0.01; training_epochs = 100;
                 ml_mlse_equalizer = ML_MLSE("epochs_tr",training_epochs,"epochs_dd",1, ...
-                    "len_tr",length(Scpe_sig),"mu_dd",mu_ml,"mu_tr",mu_ml,"order",11,"sps",2, ...
+                    "len_tr",length(Scpe_sig)/4,"mu_dd",mu_ml,"mu_tr",mu_ml,"order",11,"sps",2, ...
                     "traceback_depth",128,"L",1,"delta",4,"adaptive_mu",0);
 
                 [ml_mlse_results] = ml_mlse(ml_mlse_equalizer, M, Scpe_sig, Symbols, Tx_bits,"precode_mode",duob_mode);
diff --git a/Functions/EQ_structures/ffe.m b/Functions/EQ_structures/ffe.m
index df4138e..c7e7c52 100644
--- a/Functions/EQ_structures/ffe.m
+++ b/Functions/EQ_structures/ffe.m
@@ -157,7 +157,7 @@ function displayAnalysis(eq_noise, eq_signal_sd, rx_signal, eq_, tx_symbols, M,
 
     % Initialize figure handles
     % Corrected line - added tx_symbols as second positional argument
-    showLevelScatter(rx_signal.resample("fs_out", tx_symbols.fs), tx_symbols, "fignum", 100);
+    % showLevelScatter(rx_signal.resample("fs_out", tx_symbols.fs), tx_symbols, "fignum", 100);
 
     warning off
     showLevelScatter(eq_signal_sd, tx_symbols, "fignum", 101);
@@ -174,19 +174,20 @@ function displayAnalysis(eq_noise, eq_signal_sd, rx_signal, eq_, tx_symbols, M,
     end
 
     try
-        figure(339);
-        showEQfilter(eq_.e_tr, eq_signal_sd.fs.*2,"displayname",'training','fignum',339);
-        showEQfilter(eq_.e, eq_signal_sd.fs.*2,"displayname",'dec. directed','fignum',339);
+        figure(339);hold on
+        showEQfilter(eq_.e, eq_signal_sd.fs.*2,"displayname",'training','fignum',339);
+        % showEQfilter(eq_.e, eq_signal_sd.fs.*2,"displayname",'dec. directed','fignum',339);
         legend on
     end
 
-    try
-        figure(240); hold on; plot(pow2db(movmean(eq_.debug_struct.error_tr',100)));ylim([-30,3]);title('error training');
-    
-        figure(241); hold on; plot(pow2db(movmean(eq_.debug_struct.update_tr',100)));title('update step training');
-    
-        figure(242); hold on; plot(pow2db(movmean(eq_.debug_struct.update',1000)));title('update step dd');
-    end
-    % eq_signal_sd.eye(eq_signal_sd.fs,M,"displayname",'Eye','fignum',105);
+    % try
+    %     figure(240); hold on; plot(pow2db(movmean(eq_.debug_struct.error_tr',100)));ylim([-30,3]);title('error training');
+    % 
+    %     figure(241); hold on; plot(pow2db(movmean(eq_.debug_struct.update_tr',100)));title('update step training');
+    % 
+    %     figure(242); hold on; plot(pow2db(movmean(eq_.debug_struct.update',1000)));title('update step dd');
+    % end
+
+    eq_signal_sd.eye(eq_signal_sd.fs,M,"displayname",'Eye','fignum',105);
 
 end
\ No newline at end of file
diff --git a/Functions/EQ_structures/ml_mlse.m b/Functions/EQ_structures/ml_mlse.m
index f2d468d..ce2112c 100644
--- a/Functions/EQ_structures/ml_mlse.m
+++ b/Functions/EQ_structures/ml_mlse.m
@@ -45,8 +45,7 @@ end
 
 ml_mlse_results.config = Equalizerstruct();
 
-eq_small = strip_eq(eq_, 10);
-json_str = jsonencode(eq_small);
+eq_ = strip_eq(eq_, 10);
 
 ml_mlse_results.config.eq = jsonencode(eq_);
 ml_mlse_results.config.equalizer_structure = int32(equalizer_structure.ml_mlse);
@@ -135,6 +134,9 @@ function eq_out = strip_eq(eq_, max_elems)
                 eq_.(props{i}) = [];
             end
         end
+        if issparse(val)
+            eq_.(props{i}) = find(eq_.(props{i}));
+        end
     end
     eq_out = eq_;
 end
diff --git a/Functions/EQ_structures/vnle_postfilter_mlse.m b/Functions/EQ_structures/vnle_postfilter_mlse.m
index c7f10e1..0e6aea5 100644
--- a/Functions/EQ_structures/vnle_postfilter_mlse.m
+++ b/Functions/EQ_structures/vnle_postfilter_mlse.m
@@ -283,7 +283,7 @@ if ~isempty(postFFE)
     showEQcoefficients('n1', postFFE.e, "displayname", 'Coefficients', 'fignum', 338);
 end
 
-showEQcoefficients('n1', eq_.e,'n2', eq_.e2,'n3', eq_.e3, "displayname", 'Coefficients', 'fignum', 339);
+% showEQcoefficients('n1', eq_.e1,'n2', eq_.e2,'n3', eq_.e3, "displayname", 'Coefficients', 'fignum', 339);
 showEQfilter(eq_.e, eq_signal_sd.fs.*2);
 
 figure(340); clf;
diff --git a/Libs/Violinplot-Matlab-master/violinplot.m b/Libs/Violinplot-Matlab-master/violinplot_community.m
similarity index 100%
rename from Libs/Violinplot-Matlab-master/violinplot.m
rename to Libs/Violinplot-Matlab-master/violinplot_community.m
diff --git a/projects/FSO_transmission/first_analysis.m b/projects/FSO_transmission/first_analysis.m
index c47a2b1..d8361bc 100644
--- a/projects/FSO_transmission/first_analysis.m
+++ b/projects/FSO_transmission/first_analysis.m
@@ -1,15 +1,15 @@
 
-base = "C:\Users\Silas\Nextcloud\Dokumente\02_Ablage_Office\FSO_FP_QCL_60umUTC";
+base = "C:\Users\Silas\Nextcloud4\Dokumente\02_Ablage_Office\FSO_FP_QCL_60umUTC";
 mode = 0; %0 oder 1
 M = 2;
 
 all_files = dir(fullfile(base, "**/*.mat"));
 
 if M == 2
-    tx_data = load("C:\Users\Silas\Nextcloud\Dokumente\02_Ablage_Office\FSO_FP_QCL_60umUTC\14G_PAM2\tx_info\tx_info_PAM2_14Gbd0.75RRC.mat");
+    tx_data = load("C:\Users\Silas\Nextcloud4\Dokumente\02_Ablage_Office\FSO_FP_QCL_60umUTC\14G_PAM2\tx_info\tx_info_PAM2_14Gbd0.75RRC.mat");
     filename = fullfile(base, "14G_PAM2\M=2_Rs=1.4e10_Fs=8e10_I=265mA_RoP=46.3mW_L=31m_PS=RRC_rolloff=0.75_Mode=Rise.mat");
 elseif M == 4
-    tx_data = load("C:\Users\Silas\Nextcloud\Dokumente\02_Ablage_Office\FSO_FP_QCL_60umUTC\6G_PAM4\tx_info\tx_info_PAM4_6Gbd0.6RRC.mat");
+    tx_data = load("C:\Users\Silas\Nextcloud4\Dokumente\02_Ablage_Office\FSO_FP_QCL_60umUTC\6G_PAM4\tx_info\tx_info_PAM4_6Gbd0.6RRC.mat");
     filename = fullfile(base, "6G_PAM4\M=4_Rs=6e9_Fs=8e10_I=255mA_RoP=42.3mW_L=31m_PS=RRC_rolloff=0.6_Mode=Rise.mat");
 end
 
@@ -53,10 +53,14 @@ Bits_ = PM.demap(Symbols);
 assert(ber == 0);
 
 %% For comparison, apply pulsef on Tx Symbols
-Pform = Pulseformer("fsym",fsym,"fdac",fs,"pulse","rc","pulselength",16,"alpha",rolloff);
-Digi_sig_compare = Pform.process(Symbols);
-MF = Pulseformer("fsym",fsym,"fdac",2*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff);
-Rx_sig_compare = MF.process(Digi_sig_compare);
+Pform = Pulseformer("fsym",fsym,"fdac",fs,"pulse","rrc","pulselength",16,"alpha",rolloff);
+Digi_sig_tx_compare = Pform.process(Symbols);
+
+Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff);
+Digi_sync = Pform.process(Symbols);
+
+MF = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff);
+Rx_sig_compare = MF.process(Digi_sig_tx_compare);
 
 %%
 
@@ -68,53 +72,47 @@ scoperead_volts = (traceData.RawData - traceData.YReference) * traceData.YIncrem
 demystified = isequal(traceData.YData,scoperead_volts);
 assert(demystified);
 
+timesig_compare = [0:1:datas.tr.lastData(1).trace.ch3.Points-1] ./ fs;
+timesig = datas.tr.lastData(1).trace.ch3.XData;
+assert(isequal(traceData.YData,scoperead_volts));
+
 Scope_sig = Electricalsignal(traceData.YData,"fs",fs);
 
-Scope_sig.plot("displayname",'raw','fignum',100);
-Scope_sig.spectrum("displayname",'raw','fignum',101)
-
+%%
 % 1) matched filter
-% pulse is symmetric, hence we can use pulsef firectly as matched filter.
+% pulse is symmetric, hence we can use pulsef directly as matched filter.
 % It feels off (bit I think correct) that the fsym is now the output freq.!!
 % -> output 2 sps to omit timing recovery!?
-Pform = Pulseformer("fsym",fsym,"fdac",2*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff,"matched",1);
-Rx_matched = Pform.process(Scope_sig);
-Rx_matched.spectrum("displayname",'Signal after matched filter','fignum',1);
+apply_matched_filter = 0;
+k = 4;
+if apply_matched_filter
+    
+    Pform = Pulseformer("fsym",fsym,"fdac",k*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff,"matched",1);
+    Rx_matched = Pform.process(Scope_sig);
+else
 
+    Rx_matched = Filter('filtdegree',4,"f_cutoff",fsym*0.5,"fs",Scope_sig.fs,"filterType",filtertypes.gaussian,"active",true).process(Scope_sig);
+    Rx_matched = Rx_matched.resample("fs_out",k*fsym);
+end
+Rx_matched.spectrum();
 
-%%
-
-sys = comm.SymbolSynchronizer('TimingErrorDetector', 'Gardner (non-data-aided)', ...
-                              'SamplesPerSymbol', 2, ...
-                              'DampingFactor', 0.7, ...
-                              'NormalizedLoopBandwidth', 0.01);
-Rx_symbolsync = Rx_matched;
-[Rx_symbolsync.signal, timing_error] = sys(Rx_matched.signal);
-
-plot(timing_error); % If this is a ramp, you have drift!
-
-%% timing sync -> at this point we still have no symbol timing recovery, we
-% % try to do this with 2sps EQ!
-
-[~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_symbolsync.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1);
-
-
-%% not working.. 
-Rx_synced = Rx_synced_cell{1};
-len_tr  = 4096*2;
-mu_ffe1 = 0.0001;
-mu_ffe2 = 0.0008;
-mu_ffe3 = 0.001;
-mu_dc = 0.005;
-mu_ffe = [mu_ffe1 mu_ffe3 mu_ffe3];
-mu_dfe = 0.0004;
-duob_mode = db_mode.no_db;
-
-Rx_synced.plot("displayname",'RX: Matched+Sync+2sps','fignum',2);
-
-Digi_sig_compare.normalize("mode","rms").spectrum("displayname",'Tx: RC-shaped','fignum',1,'normalizeTo0dB',0);
-Rx_sig_compare.normalize("mode","rms").spectrum("displayname",'Tx: RC-shaped + matched filtered ','fignum',1,'normalizeTo0dB',0);
-Rx_synced.normalize("mode","rms").spectrum("displayname",'RX: matched filtered + synced','fignum',1,'normalizeTo0dB',0);
+%% Timing Rec
+apply_timing_rec = 1;
+if apply_timing_rec
+    [Rx_symbolsync, timing_error] = Timing_Recovery("timing_error_detector",'Gardner (non-data-aided)','sps',k,'damping_factor',0.1,'normalized_loop_bandwidth',0.1,'detector_gain',2.7).process(Rx_matched);
+    figure();plot(timing_error);
+    Rx_symbolsync.fs = fsym;
+    % Tsynch
+    [~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_symbolsync.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 0);
+    Rx_synced = Rx_synced_cell{1};
+    sps = 1;
+else
+    % Tsynch
+    [~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_matched.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 0);
+    Rx_synced = Rx_synced_cell{1};
+    Rx_synced = Rx_synced.resample("fs_out",2*fsym);
+    sps = 2;
+end
 
 if M == 2
     ber_in_paper = 10^(-2.6); %fig 3a) 4 Gb/s MWIR FSO Transmission using Directly Modulated QCL and an Uncooled UTC-PD at Room-Temperature
@@ -122,50 +120,121 @@ elseif M == 4
     ber_in_paper = 10^(-2.5);
 end
 
-%% -------------------- FFE --------------------
-% requires some more digging what is going on :-) 
-eq_ffe      = EQ("Ne",[50, 5, 5],"Nb",[2,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, ...
+
+%%
+
+Rx_synced = Rx_synced_cell{1};
+
+
+% -------------------- FFE --------------------
+% requires some more digging what is going on :-)
+eq_ffe      = EQ("Ne",[50, 1, 1],"Nb",[2,0,0], ...
+    "training_length",512,"training_loops",5,"dd_loops",5, ...
+    "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ...
     "FFEmu",0,"plotfinal",0,"ideal_dfe",1);
 
-ffe_results = ffe(eq_ffe,M,Rx_synced,Symbols,Bits, ...
-    "precode_mode",duob_mode,'showAnalysis',0,"postFFE",[], ...
-    "eth_style_symbol_mapping",mapping_style);
+vars = logspace(-4,-3,36);
 
-% ffe_results.metrics.print
-fprintf('My EQ: %.1e \n',ffe_results.metrics.BER);
+parfor i = 1:numel(vars)
+
+    len_tr  = 4096;
+    mu_ffe1 = 0.01;% mus(i);%0.0001;
+    mu_ffe2 = 0.0008;
+    mu_ffe3 = 0.001;
+    mu_dc = 0.005;
+    mu_ffe = [mu_ffe1 mu_ffe2 mu_ffe3];
+    mu_dfe = vars(i);
+    duob_mode = db_mode.no_db;
+
+    % requires some more digging what is going on :-)
+    eq_ffe_1      = EQ("Ne",[150, 1, 0],"Nb",[50,0,0], ...
+        "training_length",len_tr,"training_loops",5,"dd_loops",5, ...
+        "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ...
+        "FFEmu",0,"plotfinal",0,"ideal_dfe",0);
+
+    eq_ffe_2 = FFE("epochs_tr",1,"epochs_dd",vars(i),"len_tr",4096,"mu_dd",vars(i),"mu_tr",vars(i),"order",999,"sps",1,"decide",0, "adaption",adaption_method.nlms,"dd_mode",0);
+    % eq_ffe_2 = FFE_DFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"ffe_mu_dd",1e-5,"dfe_mu_dd",mus(i),"ffe_mu_tr",0,"dfe_mu_tr",0,"ffe_order",50,"dfe_order",10,"sps",1,"decide",1);
+
+
+    ffe_results = ffe(eq_ffe_1,M,Rx_synced,Symbols,Bits, ...
+        "precode_mode",duob_mode,'showAnalysis',1,"postFFE",[], ...
+        "eth_style_symbol_mapping",mapping_style);
+
+    ffe_results.metrics.BER
+    bers(i) = ffe_results.metrics.BER;
+end
+
+figure();
+plot(vars,bers);
+yline(ber_in_paper)
+beautifyBERplot();
+% 
 fprintf('Paper: %.1e \n \n',ber_in_paper);
+ffe_results.metrics.print("description",'FFE');
+fprintf('FFE: %.1e \n',ffe_results.metrics.BER);
+
 
 %% -------------------- VNLE + MLSE --------------------
+len_tr  = 4096;
+mu_ffe1 = 0.0001;% mus(i);%0.0001;
+mu_ffe2 = 0.0008;
+mu_ffe3 = 0.001;
+mu_dc = 0.005;
+mu_ffe = [mu_ffe1 mu_ffe2 mu_ffe3];
+mu_dfe = 0.0004;
+duob_mode = db_mode.no_db;
 
-pf_ncoeffs  = 1;
-eq_v        = EQ("Ne",[100, 5, 5],"Nb",[0, 0, 0], ...
+pf_ncoeffs  = 4;
+
+vars = 1:7;
+bers = zeros(size(vars));
+parfor i = 1:numel(vars)
+    
+    eqv  = EQ("Ne",[200, 1, 0],"Nb",[2, 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, ...
+    "K",sps,"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,"eth_style",mapping_style).levels);
-
-[vnle_results, mlse_results] = vnle_postfilter_mlse(eq_v, pf_, mlse_, M, Rx_synced, Symbols, Bits, ...
-    "precode_mode", duob_mode, 'showAnalysis', 1, "postFFE", [], "eth_style_symbol_mapping", mapping_style);
-
-mlse_results.metrics.print("description",'MLSE')
-fprintf('My EQ: %.1e \n',mlse_results.metrics.BER);
-fprintf('Paper: %.1e \n \n',ber_in_paper);
+    pf_ncoeffs  = vars(i);
+    pf      = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1);
+    mlse_   = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0,"eth_style",mapping_style).levels);
+    [vnle_results, mlse_results] = vnle_postfilter_mlse(eqv, pf, mlse_, M, Rx_synced, Symbols, Bits, ...
+        "precode_mode", duob_mode, 'showAnalysis', 1, "postFFE", [], "eth_style_symbol_mapping", mapping_style);
+    fprintf('Paper: %.1e \n \n',ber_in_paper);
+    vnle_results.metrics.print("description",'VNLE');
+    mlse_results.metrics.print("description",'MLSE');
+    bers(i) = mlse_results.metrics.BER;
+end
+%%
+figure();hold on
+plot(vars,bers_ffe,'DisplayName','FFE [200,0,0] + PF + MLSE');
+plot(vars,bers_vnle,'DisplayName','VNLE [200,1,0] + PF + MLSE');
+plot(vars,bers_vnledfe,'DisplayName','VNLE [200,1,0] + DFE [2] + PF + MLSE');
+plot(vars,bers_vnledfe_ideal,'DisplayName','VNLE [200,1,0] + ideal DFE [2] + PF + MLSE');
+yline(ber_in_paper);
+yline([2e-2, 4.85e-3, 3.8e-3, 2,2e-4],'LineWidth',2,'Color',[0.8,0.8,0.8],'LineStyle',':','HandleVisibility','off');
+ylim([1e-5,0.1]);
+beautifyBERplot();
 
 %% -------------------- DB target --------------------
 mlse_db_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",M,'trellis_states',PAMmapper(M,0).levels);
 
 eq_       = EQ("Ne",[50, 5, 5],"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);
+    "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1);
 
 dbt_results = duobinary_target(eq_,mlse_db_, M, Rx_synced, Symbols, Bits, ...
     "precode_mode", duob_mode, 'showAnalysis', 0, "postFFE", [],"eth_style_symbol_mapping",mapping_style);
 
-dbt_results.metrics.print("description",'Duobinary');
-mlse_results.metrics.print
-fprintf('My EQ: %.1e \n',dbt_results.metrics.BER);
 fprintf('Paper: %.1e \n \n',ber_in_paper);
+dbt_results.metrics.print("description",'Duobinary');
 
 
+%%
+
+%ML-based MLSE (L=2)
+mu_ml = 0.01; training_epochs = 100;
+ml_mlse_equalizer = ML_MLSE("epochs_tr",training_epochs,"epochs_dd",1, ...
+    "len_tr",len_tr,"mu_dd",mu_ml,"mu_tr",mu_ml,"order",11,"sps",1, ...
+    "traceback_depth",128,"L",3,"delta",4,"adaptive_mu",0);
+
+[ml_mlse_results] = ml_mlse(ml_mlse_equalizer, M, Rx_synced, Symbols, Bits,"precode_mode",duob_mode);
+ml_mlse_results.metrics.print("description",'ML ');
\ No newline at end of file
diff --git a/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/FIGURE_EYES.m b/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/FIGURE_EYES.m
index c22e90f..0d953af 100644
--- a/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/FIGURE_EYES.m
+++ b/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/FIGURE_EYES.m
@@ -2,7 +2,7 @@ dsp_options.storage_path    = 'Z:\2024\sioe_labor\';
 dsp_options.max_occurences = 1;
 database    = DBHandler("dataBase", 'labor_highspeed', "type",  'mysql' );
 
-rate = [300e9];
+
 cols = cbrewer2('BuPu',25);
 cols = [cols(end-10:2:end,:)];
 cols = cbrewer2('Set1',6);
@@ -12,10 +12,10 @@ fig=figure(fignum);clf;
 
 dbmode = 0;
 
-
 % 1 - PAM 4 with preemphasis
 fp = QueryFilter();
-M = 6;
+M = 8;
+rate = [360e9];
 fp.where('Runs', 'pam_level','EQUALS', M);
 fp.where('Runs', 'bitrate','EQUALS', rate);%360,390
 fp.where('Runs', 'fiber_length','EQUALS', 2);
@@ -23,7 +23,7 @@ fp.where('Runs', 'wavelength','EQUALS', 1310);
 fp.where('Runs', 'db_mode','EQUALS', dbmode);
 fp.where('Runs', 'rop_attenuation','EQUAL', 0);
 
-[dataTable,~] = db.queryDB(fp, database.getTableFieldNames('Runs'));
+[dataTable,~] = database.queryDB(fp, database.getTableFieldNames('Runs'));
 
 dataTable   = queryRunid(dataTable.run_id, database);
 fsym        = dataTable.symbolrate;
@@ -33,14 +33,36 @@ duob_mode   = db_mode(strrep(dataTable.db_mode,'"',''));
 % Load and Sync signal data from DB
 [Tx_bits, Symbols, Scpe_cell, ~] = loadAndSyncSignalDataFromDb(dataTable, dsp_options);
 
+Scpe_sig_syncd = Scpe_cell{1};
+Scpe_sig_syncd.eye(fsym,M,"fignum",rate.*1e-9*M+1,"displayname",' Eye of Signal');
+%%%%%% SNR CHEAT - Avges the measured signal occurences found after correlation in "tsynch" %%%%%%
+average_signals = 1;
+if average_signals
+    Scpe_sig_avg = Scpe_sig_syncd;
+    scope_mean = zeros(size(Scpe_cell{1}.signal));
+    for n=1:numel(Scpe_cell)
+        scope_mean = scope_mean + Scpe_cell{n}.signal;
+    end
+    scope_mean = scope_mean ./ n;
+    Scpe_sig_avg.signal = scope_mean;
+
+    Scpe_sig_avg.spectrum("displayname","Scope PSD","fignum",20,"normalizeTo0dB",1);
+    Scpe_sig_avg.plot("displayname","Scope raw signal","fignum",27,"clear",1);
+    Scpe_sig_avg = Scpe_sig_avg.*1.25;
+    Scpe_sig_avg.eye(fsym,M,"fignum",rate.*1e-9*M,"displayname",' Eye of AVG Signal');
+end
+
 % Preprocess signal
-Scpe_sig = preprocessSignal(Scpe_cell{1}, Symbols, fsym);
+Scpe_sig = preprocessSignal(Scpe_sig_avg, Symbols, fsym);
 
 Scpe_sig.eye(fsym,M,"fignum",M*10);
 
+
+
 %% === EXPORT TO TIKZ ===
-% outfile = ['C:\Users\Silas\Documents\latex\JLT_400G copy\media\matlab2tikz\eye_pam_',num2str(M),'.tikz'];
-% outfile = ['C:\Users\Silas\Documents\latex\JLT_400G copy\media\matlab2tikz\vnle_optimization.tikz'];
+
+% outfile = ['C:\Users\Silas\Documents\latex\JLT_400G_submission\media\matlab2tikz\eye_pam_',num2str(M),'-2.tikz'];
+% % outfile = ['C:\Users\Silas\Documents\latex\JLT_400G copy\media\matlab2tikz\vnle_optimization.tikz'];
 % matlab2tikz(outfile, ...
 %     'width','\fwidth', ...
 %     'height','\fheight', ...
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 8f96ba2..ab31451 100644
--- a/projects/HighSpeedExperiment_2024/Auswertung_JLT/run_dsp_from_db.m
+++ b/projects/HighSpeedExperiment_2024/Auswertung_JLT/run_dsp_from_db.m
@@ -42,14 +42,14 @@ fp = QueryFilter();
 % fp.where('Runs', 'run_id','EQUALS', 2776);
 M = 6;
 fp.where('Runs', 'pam_level','EQUALS', M);     
-% fp.where('Runs', 'bitrate','EQUALS', 390e9);%360,390
+fp.where('Runs', 'bitrate','EQUALS', 360e9);%360,390
 % fp.where('Runs', 'symbolrate','EQUALS', 195e9);
-fp.where('Runs', 'fiber_length','EQUALS', 10);   
+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','EQUALS', 1310);         
+fp.where('Runs', 'wavelength','EQUALS', 1310);         
 fp.where('Runs', 'db_mode','EQUALS', 0);              
 fp.where('Runs', 'rop_attenuation','EQUAL', 0);      
 % fp.where('Runs', 'power_pd_in','LESS_THAN', 7);    
@@ -70,7 +70,7 @@ wh.addStorage("mlmlse_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);
 
 
 
diff --git a/projects/HighSpeedExperiment_2024/baudrate_evaluation.m b/projects/HighSpeedExperiment_2024/baudrate_evaluation.m
index 0a48e57..100fa1d 100644
--- a/projects/HighSpeedExperiment_2024/baudrate_evaluation.m
+++ b/projects/HighSpeedExperiment_2024/baudrate_evaluation.m
@@ -1,5 +1,5 @@
 
-filename = "C:\Users\sioe\Documents\High_Speed_Measurement_2024\baudrate_sweep_b2b\PAMX_b2b_baudrate20241024_210648_wh.mat";
+filename = "F:\2024\sioe\High Speed Messungen Oktober\baudrate_sweep_b2b\PAMX_b2b_baudrate20241024_210648_wh_final.mat";
 
 a = load(filename);
 wh = a.obj;
diff --git a/projects/IMDD_base_system/minimal_example.m b/projects/IMDD_base_system/minimal_example.m
index 89f24ff..32e1ef8 100644
--- a/projects/IMDD_base_system/minimal_example.m
+++ b/projects/IMDD_base_system/minimal_example.m
@@ -8,7 +8,7 @@ fdac            = 256e9;
 fadc            = 256e9;
 random_key      = 1;
 
-rcalpha         = 0.05;
+rcalpha         = 0.6;
 kover           = 16;
 
 duob_mode = db_mode.no_db;
@@ -25,7 +25,7 @@ tx_bw_nyquist   = 0.8;
 link_length     = 1;
 
 % RX
-rop = -8;
+rop = -9;
 rx_bw_nyquist = 0.8;
 
 vnle_order1 = 50;
@@ -51,7 +51,7 @@ mu_ffe = [mu_ffe1 mu_ffe3 mu_ffe3];
 mu_dfe = 0.0004;
 
 
-Pform       = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",16,"alpha",rcalpha);
+Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"alpha",rcalpha);
 
 [Digi_sig,Symbols,Tx_bits] = PAMsource(...
     "fsym",fsym,"M",M,"order",18,"useprbs",0,...
@@ -65,7 +65,7 @@ Pform       = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",1
 %%%%% AWG
 El_sig = M8199A("kover",kover).process(Digi_sig);
 % El_sig = AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",0,"kover",kover,"bit_resolution",12,"upsampling_method","samplehold","precomp_sinc_rolloff",1).process(Digi_sig);
-El_sig.spectrum("displayname",'Digi Spectrum','fignum',100,'normalizeTo0dB',0);
+El_sig.spectrum("displayname",'Digi Spectrum','fignum',1,'normalizeTo0dB',1);
 % El_sig = El_sig.setPower(0,"dBm");
 
 %%%%% Electrical Driver Amplifier %%%%%%
@@ -101,32 +101,53 @@ 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(Rx_sig);
+%%
 
-%%%%%% Sample to 2x fsym %%%%%%
-Scpe_sig =  Scpe_sig.resample("fs_out",2*fsym);
-Scpe_sig.signal = Scpe_sig.signal(1:2*length(Symbols));
+% 1) matched filter
+% pulse is symmetric, hence we can use pulsef firectly as matched filter.
+% It feels off (bit I think correct) that the fsym is now the output freq.!!
+% -> output 2 sps to omit timing recovery!?
+Pform = Pulseformer("fsym",fsym,"fdac",2*fsym,"pulse","rrc","pulselength",16,"alpha",rcalpha,"matched",1);
+Scpe_sig = Pform.process(Scpe_sig);
+Scpe_sig.spectrum("displayname",'Signal after matched filter','fignum',1,'normalizeTo0dB',1);
+% 
 
+% %%
+% %%%%%% Sample to 2x fsym %%%%%%
+% Scpe_sig =  Scpe_sig.resample("fs_out",2*fsym);
+% Scpe_sig.signal = Scpe_sig.signal(1:2*length(Symbols));
+
+%%
 %%%%%% Sync Rx signal with reference %%%%%%
-[Scpe_sig,~] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym,"debug_plots",0);
+[Scpe_sig,~] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym,"debug_plots",1);
 Scpe_sig.spectrum("displayname",'Opt Spectrum','fignum',11,'normalizeTo0dB',1);
 
-Scpe_sig = Filter('filtdegree',4,"f_cutoff",Symbols.fs.*0.5,"fs",Scpe_sig.fs,"filterType",filtertypes.gaussian,"active",true).process(Scpe_sig);
+% Scpe_sig = Filter('filtdegree',4,"f_cutoff",Symbols.fs.*0.5,"fs",Scpe_sig.fs,"filterType",filtertypes.gaussian,"active",true).process(Scpe_sig);
 
 Scpe_sig = Scpe_sig - mean(Scpe_sig.signal);
+Scpe_sig.signal = Scpe_sig.signal(1:2*length(Symbols));
+
+%%
 
 % -------------------- FFE --------------------
 ffe_order   = [50, 0, 0];
-eq_ffe      = EQ("Ne",ffe_order,"Nb",[0,0,0], ...
+eq_      = EQ("Ne",ffe_order,"Nb",[2,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);
 
-output.ffe_results = ffe(eq_ffe,M,Scpe_sig,Symbols,Tx_bits, ...
-    "precode_mode",duob_mode,'showAnalysis',0,"postFFE",[], ...
+% eq_ = FFE("epochs_tr",4,"epochs_dd",5,"len_tr",4096,"mu_dd",0.01,"mu_tr",0.01,"order",50,"sps",2,"decide",0, "adaption",adaption_method.nlms,"dd_mode",1);
+eq_ = FFE_DFE("epochs_tr",5,"epochs_dd",5,"len_tr",512,"ffe_mu_dd",1e-4,"dfe_mu_dd",5e-4,"ffe_mu_tr",0,"dfe_mu_tr",0,"ffe_order",99,"dfe_order",99,"sps",2,"decide",0);
+
+
+output.ffe_results = ffe(eq_,M,Scpe_sig,Symbols,Tx_bits, ...
+    "precode_mode",duob_mode,'showAnalysis',1,"postFFE",[], ...
     "eth_style_symbol_mapping",0);
 
 output.ffe_results.metrics.print
 
+%%
+
 % -------------------- VNLE + MLSE --------------------
 pf_ncoeffs  = 1;
 ffe_order3  = [50, 5, 5];
@@ -139,8 +160,11 @@ pf_         = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1);
 mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels);
 
 [output.vnle_results, output.mlse_results] = vnle_postfilter_mlse(eq_v, pf_, mlse_, M, Scpe_sig, Symbols, Tx_bits, ...
-    "precode_mode", duob_mode, 'showAnalysis', 0, "postFFE", [], "eth_style_symbol_mapping", 0);
+    "precode_mode", duob_mode, 'showAnalysis', 1, "postFFE", [], "eth_style_symbol_mapping", 0);
 
+output.mlse_results.metrics.print
+
+%%
 
 % -------------------- DB target --------------------
 mlse_db_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",M,'trellis_states',PAMmapper(M,0).levels);
diff --git a/projects/Lab_analysis/data_analysis_and_plot_scripts/lambda_vs_spectrum.m b/projects/Lab_analysis/data_analysis_and_plot_scripts/lambda_vs_spectrum.m
new file mode 100644
index 0000000..c14d43d
--- /dev/null
+++ b/projects/Lab_analysis/data_analysis_and_plot_scripts/lambda_vs_spectrum.m
@@ -0,0 +1,88 @@
+wh_aeon = load("C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\Lab_analysis\aeon_soa_measurement_lambda_plaser_pump.mat");
+wh_aeon = wh_aeon.wh;
+wh_thor = load("C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\Lab_analysis\thorlabs_pdfa_measurement_lambda_plaser_pump.mat");
+wh_thor = wh_thor.wh;
+
+% Silas' custom "warehouse" datatype
+wh_aeon.showInfo;
+wh_thor.showInfo;
+
+%% PLOT OSA SPECTRA
+
+figure();hold on
+cols = linspecer(4);
+cols = cbrewer2('Paired',8);
+subplot(1,2,1);hold on;
+wh = wh_thor;
+ccnt = 1;
+
+pumps = wh.parameter.pump.values;
+plasers = wh.parameter.laserpower.values;
+scnt = 1;
+for p = 1:numel(plasers)
+    
+    for pmp = 1:numel(pumps)
+        subplot(numel(plasers),numel(pumps),scnt);hold on;
+        title(sprintf('Laser: %d dBm; Pump: %d %',plasers(p),pumps(pmp)),"Interpreter","tex")
+        for lambda = wh.parameter.lambda.values
+    
+            spectrum_osa = wh_thor.getStoValue('spectrum_osa',plasers(p),lambda,pumps(pmp));
+            wavelength_osa = wh_thor.getStoValue('wavelength_osa',plasers(p),lambda,pumps(pmp));
+            plot(wavelength_osa,spectrum_osa,'DisplayName',sprintf('P_{in}: %d dB ',plasers(p)),'Color',cols(ccnt,:));
+            
+
+            spectrum_osa = wh_aeon.getStoValue('spectrum_osa',plasers(p),lambda,pumps(pmp));
+            wavelength_osa = wh_aeon.getStoValue('wavelength_osa',plasers(p),lambda,pumps(pmp));
+            plot(wavelength_osa,spectrum_osa,'DisplayName',sprintf('P_{in}: %d dB ',plasers(p)),'Color',cols(ccnt+1,:));
+
+
+            ylim([-60, 20]);
+            beautifyBERplot("logscale",false,"setmarkers",0,"setcolors",0);
+        end
+        scnt = scnt +1;
+    end
+    ccnt = ccnt+2;
+end
+
+
+xlabel('Wavelength [nm]');
+ylabel('OSNR [dB]')
+
+
+%%
+
+figure();hold on
+cols = linspecer(4);
+subplot(1,2,1);hold on;
+wh = wh_thor;
+ccnt = 1;
+
+pumps = wh.parameter.pump.values;
+plasers = wh.parameter.laserpower.values;
+scnt = 1;
+for p = 1:numel(plasers)
+    
+    for pmp = 1:numel(pumps)
+        subplot(numel(plasers),numel(pumps),scnt);hold on;
+        
+
+        ase_noise = wh_thor.getStoValue('pase_osa',plasers(p),wh.parameter.lambda.values,pumps(pmp));
+        signal = wh_thor.getStoValue('psig_osa',plasers(p),wh.parameter.lambda.values,pumps(pmp));
+        plot(wh.parameter.lambda.values,ase_noise,'DisplayName',sprintf('P_{in}: %d dB ',plasers(p)),'Color',cols(ccnt,:));
+
+
+        ase_noise = wh_aeon.getStoValue('pase_osa',plasers(p),wh.parameter.lambda.values,pumps(pmp));
+        signal = wh_aeon.getStoValue('psig_osa',plasers(p),wh.parameter.lambda.values,pumps(pmp));
+        plot(wh.parameter.lambda.values,ase_noise,'DisplayName',sprintf('P_{in}: %d dB ',plasers(p)),'Color',cols(ccnt,:),'LineStyle','-');
+
+        ylim([-60, 20]);
+        beautifyBERplot("logscale",false,"setmarkers",0,"setcolors",0);
+
+        scnt = scnt +1;
+    end
+    ccnt = ccnt+1;
+end
+
+
+xlabel('Wavelength [nm]');
+ylabel('OSNR [dB]')
diff --git a/projects/WDM/WDM_auswertung.m b/projects/WDM/WDM_auswertung.m
index 98b42b9..1435d1f 100644
--- a/projects/WDM/WDM_auswertung.m
+++ b/projects/WDM/WDM_auswertung.m
@@ -1,5 +1,5 @@
 
-base = "C:\Users\Silas\Nextcloud\Cluster";
+base = "C:\Users\Silas\Nextcloud4\Cluster";
 all_files = dir(fullfile(base, "**/*.mat"));
 
 schemes = ["co","pair","alt","seg"];
@@ -13,8 +13,8 @@ T = table('Size',[0 10], ...
 rx = "^WDM_(?<date>\d{8})_(?<time>\d{6})_n(?<node>\d+)_(?<jobid>\d+)_" + ...
      "(?<L>\d+)km_(?<Nch>\d+)ch_(?<df>\d+)ghz_(?<scheme>[a-z]+)_alpha(?<alpha>\d+(?:_\d+)?)\.mat$";
 
-for k = 1:numel(all_files)
-    f = all_files(k);
+for i = 1:numel(all_files)
+    f = all_files(i);
     folder = string(f.folder);
     file   = string(f.name);
 
@@ -44,6 +44,8 @@ for k = 1:numel(all_files)
 end
 
 %%
+
+
 idx = strcmp(T.scheme,"co") & ...
       T.alpha == 0.4 & ...
       T.date >= datetime(2026,1,9) & ...
@@ -51,13 +53,181 @@ idx = strcmp(T.scheme,"co") & ...
 
 T_sel = T(idx,:);
 
-i = 2;
-res = load(fullfile(T_sel.folder(i),T_sel.file(i)),'res');
-res = res.res;
+% ---- Load all res ----
+R = cell(numel(T_sel.file),1);
+for i = 1:numel(T_sel.file)
+    S = load(fullfile(T_sel.folder(i),T_sel.file(i)),'res');
+    S.res = strip_config_from_res(S.res);
+    R{i} = S.res;
+end
+
+res_all = combine_res_list(R);
+res_all = drop_empty_realizations(res_all);
+
+
+%%
+
+
+
+%%
+
 
 %%
 % Routine A: plot BER curves and compute crossings
-S = plot_BER_vs_ROP(res, 'fec', 3.8e-3);
+S = plot_BER_vs_ROP(res_all, 'fec', 3.8e-3);
 
 %% Routine B: violin plot (independent)
-plot_FEC_violin(res, 'tech','VNLE', 'fec',3.8e-3, 'ylim',[-10 0]);
+plot_FEC_violin(res_all, 'tech','VNLE', 'fec',3.8e-3, 'ylim',[-10 0],'eval_ptr',5);
+
+
+%%
+
+function res = strip_config_from_res(res)
+% Remove the "config" payload from every non-empty result cell
+% Works whether entries are structs-with-field or objects-with-property.
+
+techs = {'ffe','dfe','vnle','mlse','dbt'};
+
+for t = 1:numel(techs)
+    fn = techs{t};
+    if ~isfield(res, fn) || isempty(res.(fn)), continue; end
+
+    C = res.(fn);
+    if ~iscell(C), continue; end
+
+    idx = find(~cellfun('isempty', C));  % fast builtin
+    for k = 1:numel(idx)
+        x = C{idx(k)};
+
+        % Case 1: struct entry with field "config"
+        if isstruct(x) && isfield(x,'config')
+            x = rmfield(x,'config');
+
+        % Case 2: object entry with property "config"
+        elseif isobject(x) && isprop(x,'config')
+            try
+                x.config = [];  % lighter than keeping Equalizerstruct JSON
+            catch
+                % ignore if class forbids assignment
+            end
+        end
+
+        C{idx(k)} = x;
+    end
+
+    res.(fn) = C;
+end
+end
+
+
+% ===================== Local helper functions =====================
+
+function res_out = combine_res_list(R)
+% Concatenate along realization dimension (dim=3) for all techniques.
+% Requires consistent sizes in dims 1,2,4 and consistent eval_dist_km.
+
+fieldsTech = ["ffe","dfe","vnle","mlse","dbt"];
+
+% Start with first
+res_out = R{1};
+
+% --- Build common distance axis (union, stable) ---
+dist_all = res_out.eval_dist_km(:).';
+for k = 2:numel(R)
+    dist_all = unique([dist_all, R{k}.eval_dist_km(:).'], 'stable');
+end
+
+% --- If res_out not already on dist_all, expand it ---
+if ~isequal(res_out.eval_dist_km(:).', dist_all)
+    res_out = expand_res_dist(res_out, dist_all);
+end
+
+% --- For each file: expand to dist_all, then merge realizations (your logic) ---
+for k = 2:numel(R)
+    R{k} = expand_res_dist(R{k}, dist_all);
+end
+
+
+% Concatenate technique cell arrays along dim=3 (realizations)
+for f = fieldsTech
+    A = res_out.(f);
+    for k = 2:numel(R)
+        B = R{k}.(f);
+        A = cat(3, A, B);
+    end
+    res_out.(f) = A;
+end
+
+% Update num_realiz in settings to match concatenated dim
+dims = size(res_out.ffe); if numel(dims)<3, dims(3)=1; end
+res_out.settings.num_realiz = dims(3);
+end
+
+function res = expand_res_dist(res, dist_all)
+techs = {'ffe','dfe','vnle','mlse','dbt'};
+
+old = res.eval_dist_km(:).';
+new = dist_all(:).';
+[~,loc] = ismember(old, new);              % mapping old -> new positions
+if any(loc==0)
+    error('expand_res_dist: internal: old distances not found in union.');
+end
+
+% sizes from an existing field (prefer ffe)
+Cref = res.ffe;
+sz = size(Cref); sz(end+1:4) = 1;          % ensure 4 dims
+Nch=sz(1); Nrop=sz(2); Nreal=sz(3); Nnew=numel(new);
+
+for t = 1:numel(techs)
+    fn = techs{t};
+    if ~isfield(res,fn) || isempty(res.(fn)), continue; end
+    Cold = res.(fn);
+    sz2 = size(Cold); sz2(end+1:4) = 1;
+
+    Cnew = cell(sz2(1), sz2(2), sz2(3), Nnew);
+    Cnew(:,:,:,loc) = Cold;                % place old distances into new axis
+    res.(fn) = Cnew;
+end
+
+res.eval_dist_km = new;
+end
+
+
+function res_out = drop_empty_realizations(res_in)
+% Removes realizations where ALL entries are empty across ALL techniques
+% (across channels, rop, distances).
+
+    res_out = res_in;
+    fieldsTech = ["ffe","dfe","vnle","mlse","dbt"];
+    
+    % Determine sizes from one field
+    dims = size(res_in.ffe);
+    if numel(dims) < 4, dims(end+1:4) = 1; end
+    Nreal = dims(3);
+    
+    % For each realization, check if there is at least one non-empty result anywhere
+    keep = false(1, Nreal);
+    for r = 1:Nreal
+        hasAny = false;
+        for f = fieldsTech
+            C = res_in.(f); % cell array
+            % Slice all ch,rop,dist for this realization
+            slice = C(:,:,r,:);                % still a cell array
+            hasAny = any(~cellfun(@isempty, slice(:)));
+            if hasAny, break; end
+        end
+        keep(r) = hasAny;
+    end
+    
+    fprintf('drop_empty_realizations: keeping %d/%d realizations (dropping %d)\n', ...
+        sum(keep), Nreal, sum(~keep));
+    
+    % Apply keep mask
+    for f = fieldsTech
+        res_out.(f) = res_out.(f)(:,:,keep,:);
+    end
+    
+    % Update settings
+    res_out.settings.num_realiz = sum(keep);
+end
+

From 616945bf2700ee43675fe80160a3e3d6cf13515a Mon Sep 17 00:00:00 2001
From: magf <magf@tf.uni-kiel.de>
Date: Tue, 20 Jan 2026 13:36:21 +0100
Subject: [PATCH 3/7] - Fixed 'ML_MLSE' - Added weighted DFE to 'EQ' (work in
 progress). - Added the folders 'Evaluation Scripts' and 'Data' in
 'projects\FSO_transmission', which contain different evaluations and measured
 data regarding the FSO transmission data.

---
 Classes/00_signals/Signal.m                   |   6 +-
 Classes/04_DSP/Equalizer/EQ.m                 |  48 ++
 Classes/04_DSP/Equalizer/ML_MLSE.m            | 513 +++++++++++++++++-
 Classes/04_DSP/Timing_Recovery.m              |  20 +-
 Functions/EQ_structures/ml_mlse.m             |  21 +-
 projects/FSO_transmission/first_analysis.m    |  16 +-
 .../theoretic_channel_evaluation.m            |  60 +-
 test/duobinary_minimal_example.m              |   3 +
 8 files changed, 635 insertions(+), 52 deletions(-)

diff --git a/Classes/00_signals/Signal.m b/Classes/00_signals/Signal.m
index 5905556..8fef781 100644
--- a/Classes/00_signals/Signal.m
+++ b/Classes/00_signals/Signal.m
@@ -172,9 +172,11 @@ classdef Signal
 
             hold on;
             if isempty(options.color)
-                plot(t* 1e6, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Marker', '.', 'LineStyle','none', 'MarkerSize', 0.1);
+                % plot(t* 1e6, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Marker', '.', 'LineStyle','none', 'MarkerSize', 0.1);
+                plot(t* 1e6, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1);
             else
-                plot(t* 1e6, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Marker', '.', 'LineStyle','none', 'MarkerSize', 0.1,'Color',options.color);
+                % plot(t* 1e6, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Marker', '.', 'LineStyle','none', 'MarkerSize', 0.1,'Color',options.color);
+                plot(t* 1e6, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Color',options.color);
             end
             % 2 c)
             % - xlabel if not already here: time in readable format (1 ms and not 1e-3 s)
diff --git a/Classes/04_DSP/Equalizer/EQ.m b/Classes/04_DSP/Equalizer/EQ.m
index a03d66b..f9d58d2 100644
--- a/Classes/04_DSP/Equalizer/EQ.m
+++ b/Classes/04_DSP/Equalizer/EQ.m
@@ -10,6 +10,10 @@ classdef EQ < handle
         training_length %Number of training symbols
         training_loops %Number of loops through sequence for training mode
         ideal_dfe %Error free DFE decisions
+        weighted_DFE %Weighted DFE on/off
+        weighted_DFE_mode %Weighted DFE mode
+        weighted_DFE_d_min %d_min threshold parameter for weighted DFE mode 1
+        weighted_DFE_I_mode %[a_s, b_s, I_max]-parameters for the weighted DFE
 
         DB_aim %Aim at duobinary output sequence
 
@@ -68,6 +72,10 @@ classdef EQ < handle
                 options.training_length = 1024 %Number of training symbols
                 options.training_loops = 1 %Number of loops through sequence for training mode
                 options.ideal_dfe = 0 %Error free DFE decisions
+                options.weighted_DFE = 0;
+                options.weighted_DFE_mode = 'R1';
+                options.weighted_DFE_d_min = 0.5;
+                options.weighted_DFE_I_mode = [5,0.5,0.6];
 
                 options.DB_aim %Aim at duobinary output sequence
 
@@ -438,6 +446,43 @@ classdef EQ < handle
                             dd_out(k) = constellation_in_(dd_idx);
                         end
 
+                        % Implementation of a weighted DFE in
+                        % order to prevent error propagation.
+                        % For further details, study [1], chapter 3.2.2 -
+                        % Modifications of DFE
+
+                        if obj.weighted_DFE
+                            % define new constellations
+                            const = unique(ref_in);
+
+                            % determine reliability factor gamma_k
+                            if output_vec(m) > min(const) && output_vec(m) < max(const)
+                                gamma_k = 1 - abs(output_vec(m) - dd_out(k));
+                            else
+                                gamma_k = 1;
+                            end
+
+                            % select mode
+                            if strcmp(obj.weighted_DFE_mode,'R1')
+                                if gamma_k >= obj.weighted_DFE_d_min
+                                    f_gamma_k = 1;
+                                else
+                                    f_gamma_k = 0;
+                                end
+                            elseif strcmp(obj.weighted_DFE_mode,'R2')
+                                f_gamma_k = gamma_k;
+                            elseif strcmp(obj.weighted_DFE_mode,'I1')
+                                nom = 1-exp(-obj.weighted_DFE_I_mode(1)*((gamma_k/obj.weighted_DFE_I_mode(2))-1));
+                                denom = 1+exp(-obj.weighted_DFE_I_mode(1)*((gamma_k/obj.weighted_DFE_I_mode(2))-1));
+                                f_gamma_k = (1/2)*((nom/denom) - 1);
+                            elseif strcmp(obj.weighted_DFE_mode,'I2')
+                                nom = 1-exp(-obj.weighted_DFE_I_mode(1)*((gamma_k/obj.weighted_DFE_I_mode(2))-1));
+                                denom = 1+exp(-obj.weighted_DFE_I_mode(1)*((gamma_k/obj.weighted_DFE_I_mode(2))-1));
+                                f_gamma_k = (obj.weighted_DFE_I_mode(3)/2)*((nom/denom) - 1);
+                            end
+                            % calculate weighted output
+                            output_vec(m) = f_gamma_k.*dd_out(k)+(1-f_gamma_k).*output_vec(m);
+                        end
 
                         if obj.Nb(1) > 0
                             dd_DFE(2:end) = dd_DFE(1:end-1);
@@ -739,3 +784,6 @@ classdef EQ < handle
     end
 end
 
+% References
+% [1] T. J. Wettlin, “Experimental Evaluation of Advanced Digital Signal Processing for Intra-Datacenter Systems using Direct-Detection,” 2023. [Online]. Available: https://nbn-resolving.org/urn:nbn:de:gbv:8:3-2023-00703-8
+
diff --git a/Classes/04_DSP/Equalizer/ML_MLSE.m b/Classes/04_DSP/Equalizer/ML_MLSE.m
index 3f584c3..9965525 100644
--- a/Classes/04_DSP/Equalizer/ML_MLSE.m
+++ b/Classes/04_DSP/Equalizer/ML_MLSE.m
@@ -1,3 +1,499 @@
+% classdef ML_MLSE < handle
+%     % ALGORITHM DESCRIBED IN:
+%     % W. Lanneer and Y. Lefevre, “Machine Learning-Based Pre-Equalizers for
+%     % Maximum Likelihood Sequence Estimation in High-Speed PONs,”
+%     % in 2023 31st European Signal Processing Conference
+% 
+%     % Further ML Refs:
+%     % https://machinelearningmastery.com/cross-entropy-for-machine-learning/
+%     % https://docs.pytorch.org/docs/stable/generated/torch.nn.CrossEntropyLoss.html
+% 
+%     % The central idea is to overcome the (white-) noise assumption within the previously described
+%     % Viterbi algorithm, more precisely a closed-loop optimization is proposed that finds a suitable
+%     % filter-set to directly compute the branch metrics c_k (s,s^' ). These can directly be used to
+%     % carry out the conventional Viterbi algorithm. The system consists of S^L S=F linear FIR filters,
+%     % combined with one bias coefficient respectively. These filters take the received input samples to
+%     % compute the branch metrics estimates (c_k ) ̂(s,s^' ) according toThe central idea is to overcome
+%     % the (white-) noise assumption within the previously described Viterbi algorithm, more precisely
+%     % a closed-loop optimization is proposed that finds a suitable filter-set to directly compute the
+%     % branch metrics c_k (s,s^' ). These can directly be used to carry out the conventional Viterbi
+%     % algorithm. The system consists of S^L S=F linear FIR filters, combined with one bias coefficient
+%     % respectively. These filters take the received input samples to compute the branch metrics
+%     % estimates. Finally, the usual Viterbi is carried out...
+% 
+%     % Recommended Settings and some findings:
+% 
+%     % Requires many training epochs. According to ML people, 100,200 or
+%     % even up to 1000 epochs are normal for ML-convergence
+% 
+%     % The mu parameter _can_ be adaptive - using the cross entropy and when
+%     % analyzing the isolated training it looks very promisig. However, is
+%     % later use I found this is not as stable as a fixed learning rate.
+%     % mu = 0.1 worked good for me
+% 
+%     % Longer orders/ filter length are not always better. For me order=11
+%     % was good.
+% 
+%     % Delay factor (delta) is good when the order is also increased. With
+%     % order = 11, a delta of =4 shows good results
+% 
+%     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
+% 
+%         adaptive_mu
+% 
+%         constellation
+% 
+%         L %viterbi memory length
+% 
+%         alpha
+%         DIR
+%         DIR_flip
+%         trellis_states
+% 
+%         traceback_depth
+% 
+%         % --- Added internal class variables used later ---
+%         S
+%         Nf
+%         delta
+%         nStates
+%         nFeasible
+%         combs
+%         first_sym
+%         last_sym
+%         valid
+%         valid_to_idx
+%         valid_from_idx
+%         w
+% 
+%         % --- New: fast state lookup ---
+%         true_to_state_idx
+%         state_dict              % containers.Map: key(sequence)->state index
+%         key_fmt = '%.8g_';      % key format for sequence strings
+%         nSym                      % |constellation|
+% 
+%         ber = []
+%         ce  = ones(1,1);
+%     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.adaptive_mu = 1;
+% 
+%                 options.delta = 0;
+%                 options.traceback_depth = 1024;
+% 
+%                 options.L = 1
+% 
+%             end
+% 
+%             fn = fieldnames(options);
+%             for n = 1:numel(fn)
+%                 obj.(fn{n}) = options.(fn{n});
+%             end
+% 
+%             obj.e = zeros(obj.order,1);
+%             obj.error = 0;
+%         end
+% 
+%         function [X,X_viterbi] = process(obj, X, D)
+% 
+%             % actual processing of the signal (steps 1. - 3.)
+%             % 1 normalize RMS
+%             X = X.normalize("mode","rms");
+% 
+%             % 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(obj.constellation);      % alphabet size
+%             obj.Nf = obj.order*obj.sps;            % filter length
+%             % 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(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{:}).');  % rows: states, columns: [x_k, x_{k-1}, ...]
+%             obj.first_sym = obj.combs(:,1);
+%             obj.last_sym  = obj.combs(:,end);
+%             obj.nStates   = size(obj.combs,1);
+% 
+%             % --- Valid transitions
+%             obj.valid = false(obj.nStates);
+%             for from = 1:obj.nStates
+%                 for to = 1:obj.nStates
+%                     if all(obj.combs(to,2:end) == obj.combs(from,1:end-1))
+%                         obj.valid(to,from) = true;
+%                     end
+%                 end
+%             end
+%             [obj.valid_to_idx, obj.valid_from_idx] = find(obj.valid);
+% 
+%             % --- Allocate vectors and weights
+%             % !! IF SHAPE FIT, then we already have smth there an we want
+%             % to start with the existing fitler-set
+%             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
+% 
+%             % --- 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
+%             % ==============================================================
+% 
+%             % Training Mode
+%             n = obj.len_tr;
+%             training = 1;
+%             obj.equalize(X.signal, D.signal,obj.mu_tr,obj.epochs_tr,n,training);
+%             obj.e_tr = obj.e;
+% 
+%             % ==============================================================
+%             %  DD-Mode / Fixed Mode
+%             % ==============================================================
+% 
+%             % 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_ref] = equalize(obj,x,d,mu,epochs,N,training)
+%             % ==============================================================
+%             %  FFE + Whitening + ML-Based Branch Metric Estimation + Viterbi
+%             % ==============================================================
+%             debug = 1;
+%             showPlots = 1;
+% 
+%             % --- 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);  % v_{k-1}(s′)
+%                 c_hat     = zeros(1,obj.nFeasible);
+%                 v_tilde   = zeros(1,obj.nFeasible);
+%                 pred      = zeros(nSymbols, obj.nStates, 'uint32');
+%                 pm_sto    = nan(obj.nStates, nSymbols,'like',pm);
+%                 CE_accum = 0;
+% 
+% 
+%                 %%% START IDX 
+%                 if training
+%                     max_start = length(x) - ( (ceil(N/obj.sps)-1)*obj.sps + 1 );
+%                     max_start = max(1, max_start);  % safety
+%                     start_sample = randi([1, max_start], 1); %rnd training; not really good
+%                     start_sample = 1;
+%                     end_sample   = start_sample + (ceil(N/obj.sps)-1)*obj.sps;
+%                 else
+%                     start_sample = 1;%obj.len_tr;
+%                     end_sample   = N;
+%                 end
+% 
+%                 start_symbol = 1 + floor((start_sample - 1)/obj.sps);   % ABSOLUTE symbol index
+% 
+%                 if numel(d) >= obj.L && start_symbol >= obj.L
+%                     init_seq = d(start_symbol-obj.L+1 : start_symbol);      % [d_k-L+1 ... d_k]
+%                     true_to_state_idx = obj.state_dict(obj.seq_key(flip(init_seq))); % [d_k ... d_k-L+1]
+%                 else
+%                     % Not enough history – fall back to state 1
+%                     true_to_state_idx = uint32(1);
+%                 end
+% 
+%                 symbol = 0;
+%                 for sample = start_sample:obj.sps:end_sample
+%                     symbol = symbol + 1;
+%                     k = symbol;
+%                     sym_idx = start_symbol + (symbol - 1);
+% 
+%                     % --- Build Δ-delayed observation window y_k
+%                     i1  = sample - obj.Nf + 1 + obj.delta;
+%                     i2  = sample + obj.delta;
+%                     buf = x(max(1,i1):min(length(x),i2));
+%                     padL = max(0,1 - i1);
+%                     padR = max(0,i2 - length(x));
+%                     yk  = [zeros(padL,1); buf(:); zeros(padR,1)];   % Nf×1
+%                     yk  = [yk;1];
+% 
+%                     % --- 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: 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 1 %training
+%                         % --- allocate storage once
+%                         if epoch == 1 && symbol == 1
+%                             obj.true_to_state_idx = ones(ceil(N/obj.sps),1,'uint32');
+%                         end
+% 
+%                         % --- previous "to" becomes current "from"
+%                         if symbol > 1
+%                             true_from_state_idx = obj.true_to_state_idx(symbol-1);
+%                         else
+%                             true_from_state_idx = 1;
+%                         end
+% 
+%                         % --- compute or reuse "to" state
+%                         if epoch == 1
+%                             % only compute in first epoch
+%                             if sym_idx >= obj.L
+%                                 key_to = obj.seq_key(flip(d(sym_idx-obj.L+1 : sym_idx)));
+%                                 if isKey(obj.state_dict, key_to)
+%                                     obj.true_to_state_idx(symbol) = obj.state_dict(key_to);
+%                                 else
+%                                     obj.true_to_state_idx(symbol) = true_from_state_idx;
+%                                 end
+%                             else
+%                                 obj.true_to_state_idx(symbol) = true_from_state_idx;
+%                             end
+%                         end
+% 
+%                         % --- reuse cached state from second epoch onward
+%                         true_to_state_idx = obj.true_to_state_idx(symbol);
+% 
+%                         % --- ensure valid (from,to)
+%                         dirac = zeros(obj.nFeasible,1);
+%                         mask  = obj.valid_from_idx==true_from_state_idx & ...
+%                             obj.valid_to_idx  ==true_to_state_idx;
+%                         if any(mask)
+%                             dirac(mask) = 1;
+%                         else
+%                             idx = find(obj.valid_from_idx==true_from_state_idx,1,'first');
+%                             dirac(idx) = 1;
+%                             obj.true_to_state_idx(symbol) = obj.valid_to_idx(idx);
+%                         end
+% 
+% 
+% 
+% 
+%                         % softmax over -v_tilde (numerically safe shift)
+%                         v_shift = -(v_tilde - min(v_tilde));         % shift to small positive numbers
+%                         v_shift = min(v_shift, 100);                  % clamp exponent argument (≈ exp(50)=3e21)
+%                         expv = exp(v_shift);
+%                         p = expv ./ (sum(expv) + eps);
+% 
+%                         % for logging only:
+%                         CE_symbol(symbol) = -log(p(dirac==1) + eps);
+% 
+%                         if sym_idx > obj.L
+%                             CE_smooth(symbol) = 0.01*CE_symbol(symbol) + 0.99*CE_smooth(symbol-1);
+%                         else
+%                             if epoch > 1
+%                                 CE_smooth(symbol) = obj.ce(end); %use ce from last epoch or =1 for very first round?!
+%                             else
+%                                 CE_smooth(symbol) = CE_symbol(symbol);
+%                             end
+%                         end
+% 
+%                         CE_accum = CE_symbol(symbol) + CE_accum;
+% 
+% 
+%                         % gradient term (t - p)
+%                         dmp = (dirac - p)';  % 1×nFeasible
+% 
+%                         % Per-feature gradient; implicit expansion gives (Nf+1)×nFeasible
+%                         dL_Dw = (yk) .* dmp;
+% 
+%                         % Start updates only when the ABSOLUTE symbol index has ≥ L history
+%                         if sym_idx >= obj.L
+%                             if obj.adaptive_mu
+%                                 mu_eff = CE_smooth(sym_idx);
+%                                 mu_eff = max(min(mu_eff, 0.2), 1e-4);
+%                             else
+%                                 mu_eff = mu;
+%                             end
+% 
+%                             obj.w = obj.w - mu_eff .* dL_Dw;   % (Nf+1)×nFeasible
+%                         end
+% 
+%                         % if debug && epoch > 2
+%                         %     figure(100);
+%                         %     subplot(4,1,1);
+%                         %     heatmap(p');
+%                         %     title('Probs')
+%                         %     subplot(4,1,2);
+%                         %     heatmap(dmp);
+%                         %     title('Update')
+%                         %     subplot(4,1,3);
+%                         %     heatmap(dL_Dw);
+%                         %     title('Update')
+%                         %     subplot(4,1,4);
+%                         %     heatmap(bj.w);
+%                         %     title('Update')
+%                         % 
+%                         % end
+% 
+%                     end
+% 
+% 
+% 
+%                     % --- Compare-Select (matrix form, min of costs)
+%                     v_tilde_mat = inf(obj.nStates, obj.nStates);
+%                     v_tilde_mat(obj.valid) = v_tilde;
+%                     [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
+% 
+%                 % --- 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
+%                     viterbi_path(n-1) = pred(n, viterbi_path(n));
+%                 end
+% 
+%                 y_ref = d(start_symbol:end);
+%                 y     = obj.first_sym(viterbi_path);
+% 
+%                 if debug && training
+%                     sym_start = start_symbol;
+%                     sym_end   = start_symbol + symbol - 1;
+%                     ref_slice = d(sym_start : sym_end);
+%                     err = sum(y ~= ref_slice(1:numel(y)));
+% 
+%                     try
+%                         ref_bits = PAMmapper(obj.S,0).demap(ref_slice);
+%                         eq_bits = PAMmapper(obj.S,0).demap(y);
+%                         [~, ~, ber, ~] = calc_ber(ref_bits, eq_bits, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1);
+%                         fprintf('Epoch: %d - BER: %.1e  \n',epoch, ber);
+%                         obj.ber(epoch) = ber;
+%                     catch
+%                         ser = err./length(y);
+%                         fprintf('Epoch: %d - SER: %.1e  \n',epoch, ser);
+%                     end
+% 
+%                     obj.ce(epoch) = CE_accum./symbol;
+% 
+%                     if showPlots
+%                         figure(10);clf
+%                         subplot(3,2,1:2);
+%                         heatmap(obj.w);
+%                         title('Filter')
+% 
+%                         subplot(3,2,3);
+%                         v_tildemat = NaN(obj.nStates, obj.nStates);
+%                         v_tildemat(obj.valid) = v_tilde;      % log-domain scores
+%                         heatmap(v_tildemat);
+%                         title('Path Metrics (v_tilde)')
+% 
+%                         subplot(3,2,4);
+%                         scatter(1:symbol,pm_sto,1,'.')
+%                         title('Path Metric Winners')
+% 
+%                         subplot(3,2,5);hold on
+%                         scatter(1:symbol,CE_symbol,1,'.');
+%                          scatter(1:symbol,CE_smooth,1,'.')
+%                         title('Cross Entropy')
+% 
+%                         subplot(3,2,6); hold on
+% 
+%                         % Left y-axis: Cross Entropy (linear)
+%                         yyaxis left
+%                         scatter(1:length(obj.ce), obj.ce, 10, 's', 'filled')
+%                         ylabel('Cross Entropy')
+% 
+%                         % Right y-axis: BER (logarithmic)
+%                         yyaxis right
+%                         scatter(1:length(obj.ber), obj.ber, 10, 'd', 'filled')
+%                         set(gca, 'YScale', 'log')
+%                         ylabel('BER (log scale)')
+% 
+%                         xlim([1, epochs])
+%                         xlabel('Epoch')
+%                         title('Cross Entropy // BER')
+%                         grid on
+% 
+%                         drawnow
+%                     end
+%                 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
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
 classdef ML_MLSE < handle
     % ---------------------------------------------------------------------
     % W. Lanneer and Y. Lefevre,
@@ -164,8 +660,8 @@ classdef ML_MLSE < handle
         % EQUALIZE
         % ==============================================================
         function [y,y_ref] = equalize(obj,x,d,mu,epochs,N,training)
-            debug = 0;
-            showPlots = 0;
+            debug = 1;
+            showPlots = 1;
             y = zeros(N,1);
             nSymbols = ceil(N/obj.sps);
 
@@ -241,6 +737,19 @@ classdef ML_MLSE < handle
                         dirac(trans_idx)=1;
                     end
 
+                    % --- ensure valid (from,to)
+                    if ~any(dirac)
+                        mask  = obj.valid_from_idx==true_from_state_idx & ...
+                            obj.valid_to_idx  ==true_to_state_idx;
+                        if any(mask)
+                            dirac(mask) = 1;
+                        else
+                            idx = find(obj.valid_from_idx==true_from_state_idx,1,'first');
+                            dirac(idx) = 1;
+                            obj.true_to_state_idx(symbol) = obj.valid_to_idx(idx);
+                        end
+                    end
+
                     % ===================================================================
                     % TRAINING MODE (weight update)
                     % ===================================================================
diff --git a/Classes/04_DSP/Timing_Recovery.m b/Classes/04_DSP/Timing_Recovery.m
index ba16f10..3999ad7 100644
--- a/Classes/04_DSP/Timing_Recovery.m
+++ b/Classes/04_DSP/Timing_Recovery.m
@@ -1,6 +1,7 @@
 classdef Timing_Recovery < handle
 
     properties(Access=public)
+        modulation
         timing_error_detector
         sps
         damping_factor
@@ -9,14 +10,15 @@ classdef Timing_Recovery < handle
     end
 
     methods(Access=public)
-        function obj = FFE(options)
+        function obj = Timing_Recovery(options)
             arguments(Input)
 
-                options.timing_error_detector = 'Gardner';
+                options.modulation = 'PAM/PSK/QAM';
+                options.timing_error_detector = 'Gardner (non-data-aided)';
                 options.sps = 2;
                 options.damping_factor = 1.0;
-                options.normalized_loop_bandwidth = 0.005;
-                options.detector_gain = 1;
+                options.normalized_loop_bandwidth = 0.01;
+                options.detector_gain = 2.7;
 
             end
 
@@ -25,21 +27,21 @@ classdef Timing_Recovery < handle
                 obj.(fn{n}) = options.(fn{n});
             end
 
-            obj.e = zeros(obj.order,1);
-            obj.error = 0;
+            %obj-Initialization here%
 
         end
 
-        function data_out = process(obj, data_in)
+        function [data_out, timing_error] = process(obj, data_in)
 
             timing_synchronization = comm.SymbolSynchronizer( ...
+                "Modulation",          obj.modulation,...
                 "TimingErrorDetector", obj.timing_error_detector, ...
                 "SamplesPerSymbol",    obj.sps, ...
                 "DampingFactor",       obj.damping_factor, ...
                 "NormalizedLoopBandwidth", obj.normalized_loop_bandwidth, ...
                 "DetectorGain",        obj.detector_gain);
-
-            data_out.signal = timing_synchronization(data_in.signal);
+            data_out = data_in;
+            [data_out.signal, timing_error] = timing_synchronization(data_in.signal);
 
         end
     end
diff --git a/Functions/EQ_structures/ml_mlse.m b/Functions/EQ_structures/ml_mlse.m
index f2d468d..80e28cf 100644
--- a/Functions/EQ_structures/ml_mlse.m
+++ b/Functions/EQ_structures/ml_mlse.m
@@ -46,8 +46,20 @@ end
 ml_mlse_results.config = Equalizerstruct();
 
 eq_small = strip_eq(eq_, 10);
+fn = fieldnames(eq_small);
+for k = 1:numel(fn)
+    if issparse(eq_small.(fn{k}))
+        eq_small.(fn{k}) = full(eq_small.(fn{k}));
+    end
+end
 json_str = jsonencode(eq_small);
 
+fn = fieldnames(eq_);
+for k = 1:numel(fn)
+    if issparse(eq_.(fn{k}))
+        eq_.(fn{k}) = full(eq_.(fn{k}));
+    end
+end
 ml_mlse_results.config.eq = jsonencode(eq_);
 ml_mlse_results.config.equalizer_structure = int32(equalizer_structure.ml_mlse);
 ml_mlse_results.config.comment = 'function: ML-based MLSE';
@@ -81,6 +93,7 @@ end
 function [bits, errors, ber, error_pos, errors_precoded, ber_precoded] = calculateBER(eq_signal_hd, tx_symbols, tx_bits, precode_mode, M, eth_style)
 % Calculate BER based on precoding mode
 mapper = PAMmapper(M, 0, "eth_style", eth_style);
+skip_front = 150;
 
 switch precode_mode
     case db_mode.no_db
@@ -95,11 +108,11 @@ switch precode_mode
         tx_bits_precoded = mapper.demap(tx_symbols_precoded);
         
         rx_bits = mapper.demap(eq_signal_hd_precoded);
-        [~, errors_precoded, ber_precoded, ~] = calc_ber(rx_bits.signal, tx_bits_precoded.signal, "skip_front", 30000, "skip_end", 150, "returnErrorLocation", 1);
+        [~, errors_precoded, ber_precoded, ~] = calc_ber(rx_bits.signal, tx_bits_precoded.signal, "skip_front", skip_front, "skip_end", 150, "returnErrorLocation", 1);
         
         % B) Just determine BER
         rx_bits = mapper.demap(eq_signal_hd);
-        [bits, errors, ber, error_pos] = calc_ber(rx_bits.signal, tx_bits.signal, "skip_front", 30000, "skip_end", 150, "returnErrorLocation", 1);
+        [bits, errors, ber, error_pos] = calc_ber(rx_bits.signal, tx_bits.signal, "skip_front", skip_front, "skip_end", 150, "returnErrorLocation", 1);
         
     case db_mode.db_precoded
         % Data is precoded on TX side
@@ -107,12 +120,12 @@ switch precode_mode
         eq_signal_hd_decoded = Duobinary().encode(eq_signal_hd, "M", M);
         eq_signal_hd_decoded = Duobinary().decode(eq_signal_hd_decoded, "M", M);
         rx_bits_decoded = mapper.demap(eq_signal_hd_decoded);
-        [~, errors_precoded, ber_precoded, ~] = calc_ber(rx_bits_decoded.signal, tx_bits.signal, "skip_front", 30000, "skip_end", 150, "returnErrorLocation", 1);
+        [~, errors_precoded, ber_precoded, ~] = calc_ber(rx_bits_decoded.signal, tx_bits.signal, "skip_front", skip_front, "skip_end", 150, "returnErrorLocation", 1);
         
         % B) Omit the Coding by comparing with demapped TX symbol sequence
         tx_bits_demapped = mapper.demap(tx_symbols);
         rx_bits = mapper.demap(eq_signal_hd);
-        [bits, errors, ber, error_pos] = calc_ber(rx_bits.signal, tx_bits_demapped.signal, "skip_front", 30000, "skip_end", 150, "returnErrorLocation", 1);
+        [bits, errors, ber, error_pos] = calc_ber(rx_bits.signal, tx_bits_demapped.signal, "skip_front", skip_front, "skip_end", 150, "returnErrorLocation", 1);
 end
 end
 
diff --git a/projects/FSO_transmission/first_analysis.m b/projects/FSO_transmission/first_analysis.m
index c47a2b1..42ea7b7 100644
--- a/projects/FSO_transmission/first_analysis.m
+++ b/projects/FSO_transmission/first_analysis.m
@@ -1,15 +1,15 @@
 
-base = "C:\Users\Silas\Nextcloud\Dokumente\02_Ablage_Office\FSO_FP_QCL_60umUTC";
+base = "C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\FSO Equalizer\FSO_FP_QCL_60umUTC\";
 mode = 0; %0 oder 1
 M = 2;
 
 all_files = dir(fullfile(base, "**/*.mat"));
 
 if M == 2
-    tx_data = load("C:\Users\Silas\Nextcloud\Dokumente\02_Ablage_Office\FSO_FP_QCL_60umUTC\14G_PAM2\tx_info\tx_info_PAM2_14Gbd0.75RRC.mat");
+    tx_data_path = fullfile(base, "14G_PAM2\tx_info\tx_info_PAM2_14Gbd0.75RRC.mat");
     filename = fullfile(base, "14G_PAM2\M=2_Rs=1.4e10_Fs=8e10_I=265mA_RoP=46.3mW_L=31m_PS=RRC_rolloff=0.75_Mode=Rise.mat");
 elseif M == 4
-    tx_data = load("C:\Users\Silas\Nextcloud\Dokumente\02_Ablage_Office\FSO_FP_QCL_60umUTC\6G_PAM4\tx_info\tx_info_PAM4_6Gbd0.6RRC.mat");
+    tx_data_path = fullfile(base, "6G_PAM4\tx_info\tx_info_PAM4_6Gbd0.6RRC.mat");
     filename = fullfile(base, "6G_PAM4\M=4_Rs=6e9_Fs=8e10_I=255mA_RoP=42.3mW_L=31m_PS=RRC_rolloff=0.6_Mode=Rise.mat");
 end
 
@@ -20,7 +20,7 @@ if mode == 1
     end
 end
 
-
+tx_data = load(tx_data_path);
 datas = load(filename);
 
 %%
@@ -97,9 +97,15 @@ plot(timing_error); % If this is a ramp, you have drift!
 % % try to do this with 2sps EQ!
 
 [~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_symbolsync.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1);
-
+% Rx_matched_1 = Rx_synced_cell{1};
+% 
+% % Rx_Time_Rec = Rx_matched;
+% [Rx_Time_Rec, Timing_Error] = Timing_Recovery("modulation", 'PAM/PSK/QAM', "timing_error_detector",'Gardner (non-data-aided)','sps',2,'damping_factor',1,'normalized_loop_bandwidth',0.01,'detector_gain',2.7).process(Rx_matched_1);
+% figure;plot(Timing_Error);
+% % Rx_Time_Rec.signal = resample(Rx_Time_Rec.signal, 28e9, 14e9);
 
 %% not working.. 
+% Rx_synced = Rx_Time_Rec;
 Rx_synced = Rx_synced_cell{1};
 len_tr  = 4096*2;
 mu_ffe1 = 0.0001;
diff --git a/projects/ML_based_MLSE/theoretic_channel_evaluation.m b/projects/ML_based_MLSE/theoretic_channel_evaluation.m
index 3eb6f09..3e1efcc 100644
--- a/projects/ML_based_MLSE/theoretic_channel_evaluation.m
+++ b/projects/ML_based_MLSE/theoretic_channel_evaluation.m
@@ -48,36 +48,36 @@ for i = 1:numel(SNR_dB)
     symbols_noi = symbols_filt;
     symbols_noi.signal = awgn(symbols_filt.signal, SNR_dB(i), 'measured');  % AWGN with given SNR
     
-    % Sequence Est L=5
-    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',h);
-    mlse_.DIR = h;
-    [y_mlse] = mlse_.process(symbols_noi,Symbols);
-    mlse_bits = PAMmapper(M, 0, "eth_style", 0).demap(y_mlse);
-    [~, ~, ber_mlse_l5(i), ~] = calc_ber(mlse_bits.signal, Bits.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1);
-    fprintf('MLSE L5: %.2e \n',ber_mlse_l5(i));
-
-    % 2nd Approach 
-    mu_lms = 0.0005;
-    pf_ncoeffs = 1;
-    eq_ = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",2^13,"mu_dd",mu_lms,"mu_tr",mu_lms,"order",16,"sps",1,"dd_mode",1,"adaption_technique","lms");
-    pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1);
-
-    % FFE
-    [y_ffe, ffe_noise] = eq_.process(symbols_noi, Symbols);
-
-    Eq_bits = PAMmapper(M, 0, "eth_style", 0).demap(y_ffe);
-    [~, ~, ber_ffe(i), ~] = calc_ber(Eq_bits.signal, Bits.signal, "skip_front", 0, "skip_end", 0, "returnErrorLocation", 1);
-    fprintf('FFE: %.2e \n',ber_ffe(i));
-
-    % Postfilter
-    [y_white,~] = pf_.process(y_ffe, ffe_noise);
-
-    % Sequence Est
-    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);
-    [y_mlse] = mlse_.process(y_white,Symbols);
-    mlse_bits = PAMmapper(M, 0, "eth_style", 0).demap(y_mlse);
-    [~, errors, ber_nwf_mlse_l2(i), errpos] = calc_ber(mlse_bits.signal, Bits.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1);
-    fprintf('MLSE: %.2e \n',ber_nwf_mlse_l2(i));
+    % % Sequence Est L=5
+    % 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',h);
+    % mlse_.DIR = h;
+    % [y_mlse] = mlse_.process(symbols_noi,Symbols);
+    % mlse_bits = PAMmapper(M, 0, "eth_style", 0).demap(y_mlse);
+    % [~, ~, ber_mlse_l5(i), ~] = calc_ber(mlse_bits.signal, Bits.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1);
+    % fprintf('MLSE L5: %.2e \n',ber_mlse_l5(i));
+    % 
+    % % 2nd Approach 
+    % mu_lms = 0.0005;
+    % pf_ncoeffs = 1;
+    % eq_ = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",2^13,"mu_dd",mu_lms,"mu_tr",mu_lms,"order",16,"sps",1,"dd_mode",1,"adaption_technique","lms");
+    % pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1);
+    % 
+    % % FFE
+    % [y_ffe, ffe_noise] = eq_.process(symbols_noi, Symbols);
+    % 
+    % Eq_bits = PAMmapper(M, 0, "eth_style", 0).demap(y_ffe);
+    % [~, ~, ber_ffe(i), ~] = calc_ber(Eq_bits.signal, Bits.signal, "skip_front", 0, "skip_end", 0, "returnErrorLocation", 1);
+    % fprintf('FFE: %.2e \n',ber_ffe(i));
+    % 
+    % % Postfilter
+    % [y_white,~] = pf_.process(y_ffe, ffe_noise);
+    % 
+    % % Sequence Est
+    % 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);
+    % [y_mlse] = mlse_.process(y_white,Symbols);
+    % mlse_bits = PAMmapper(M, 0, "eth_style", 0).demap(y_mlse);
+    % [~, errors, ber_nwf_mlse_l2(i), errpos] = calc_ber(mlse_bits.signal, Bits.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1);
+    % fprintf('MLSE: %.2e \n',ber_nwf_mlse_l2(i));
 
     % ML-base MLSE L=2
     adaptive_mu = 0;
diff --git a/test/duobinary_minimal_example.m b/test/duobinary_minimal_example.m
index c39183e..d61c044 100644
--- a/test/duobinary_minimal_example.m
+++ b/test/duobinary_minimal_example.m
@@ -29,6 +29,9 @@ para.skip =0;
 para.bruijn = 0;
 para.reset_prms = 0;
 para.method = 1;
+para.pcs = 0;
+para.shape_para = 0;
+para.rng_num = 0;
 
 data_in = [];
 global loop;

From 1bf17ec6bee092ca4f64960b44beabefc1d77cf2 Mon Sep 17 00:00:00 2001
From: Silas Oettinghaus <silas.oettinghaus@email.uni-kiel.de>
Date: Tue, 20 Jan 2026 13:48:15 +0100
Subject: [PATCH 4/7] FSO: - try to analyze spectrum of FSO and use the coeffs
 in the ML analysis. simple channel coeff + awgn example

---
 projects/FSO_transmission/first_analysis.m     | 16 ++++++++++++++--
 projects/ML_based_MLSE/analyze_filter_length.m | 18 +++++++++++-------
 2 files changed, 25 insertions(+), 9 deletions(-)

diff --git a/projects/FSO_transmission/first_analysis.m b/projects/FSO_transmission/first_analysis.m
index d8361bc..9a327ab 100644
--- a/projects/FSO_transmission/first_analysis.m
+++ b/projects/FSO_transmission/first_analysis.m
@@ -83,8 +83,8 @@ Scope_sig = Electricalsignal(traceData.YData,"fs",fs);
 % pulse is symmetric, hence we can use pulsef directly as matched filter.
 % It feels off (bit I think correct) that the fsym is now the output freq.!!
 % -> output 2 sps to omit timing recovery!?
-apply_matched_filter = 0;
-k = 4;
+apply_matched_filter = 1;
+k = 1;
 if apply_matched_filter
     
     Pform = Pulseformer("fsym",fsym,"fdac",k*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff,"matched",1);
@@ -96,6 +96,18 @@ else
 end
 Rx_matched.spectrum();
 
+
+%%
+
+coefficients = arburg(Rx_matched.signal,25);
+
+figure()
+[h,w] = freqz(1,coefficients,Rx_matched.length,"whole",Rx_matched.fs);
+h = h/max(abs(h));
+hold on
+w_ = (w - Rx_matched.fs/2);
+plot(w_.*1e-9,20*log10(fftshift(abs(h))),'DisplayName',['Burg Coeffs: ', num2str(round(coefficients,2)), ' '],'LineWidth',2);
+
 %% Timing Rec
 apply_timing_rec = 1;
 if apply_timing_rec
diff --git a/projects/ML_based_MLSE/analyze_filter_length.m b/projects/ML_based_MLSE/analyze_filter_length.m
index 616b715..bcf2a91 100644
--- a/projects/ML_based_MLSE/analyze_filter_length.m
+++ b/projects/ML_based_MLSE/analyze_filter_length.m
@@ -5,9 +5,9 @@ M = 4;
 randkey = 1;
 
 % --- Parameter sweep
-order_range = 2:3:11;        % FFE order
-delta_range = 0:2:4;        % delta
-SNR_dB = 20;
+order_range = 5:5:50;        % FFE order
+delta_range = 0:5:20;        % delta
+SNR_dB = 30;
 
 % --- Prepare bit sequence
 order_bits = 19;
@@ -21,10 +21,14 @@ Symbols = PAMmapper(M,0).map(Bits);
 Symbols.fs = 200e9;
 
 % --- Channel (minimal ISI + AWGN)
-h = [0.3 0.9 0.3]; h = h/norm(h);
-symbols_filt = Symbols.filter(h,1);
+h = abs([0.3 0.9 0.3]); h = h/norm(h);
+
+% h = [1	-1.67085330039878	1.17918163282514	-0.805210559745616	0.571564213123367	-0.296337147529674	0.00649773445209780	0.0854177610195952	-0.0576009020965258	0.0520994427061551	-0.0624586034913656	0.0553280962699552	-0.00705582559925755	-0.0336399056707792	0.0706903719452810	-0.0334124287931977	0.0131699455037966	0.0587431373842994	-0.0515902976066452	0.00647904355473619	0.0137506750904990	-0.0547974515885928	0.00994735499340592	-0.0135513582534086	-0.00463322575007739	0.0277311946101940];
+% h = h/norm(h);
+symbols_filt = Symbols.filter(1,h);
 symbols_noi = symbols_filt;
 symbols_noi.signal = awgn(symbols_filt.signal,SNR_dB,'measured');
+symbols_noi.spectrum;
 
 % --- Generate all parameter pairs
 [O,D] = ndgrid(order_range, delta_range);
@@ -51,7 +55,7 @@ parfor k = 1:size(pairs,1)
     try
         ml = ML_MLSE("epochs_tr",training_len,"epochs_dd",1,"len_tr",2^15, ...
             "mu_dd",0.1,"mu_tr",0.1,"order",order_k,"sps",1, ...
-            "traceback_depth",128,"L",3,"delta",delta_k,"adaptive_mu",0);
+            "traceback_depth",128,"L",2,"delta",delta_k,"adaptive_mu",0);
 
         [y_ml,y_ref] = ml.process(symbols_noi,Symbols);
         ref_bits = PAMmapper(M,0).demap(y_ref);
@@ -89,7 +93,7 @@ end
 beautifyBERplot
 ylabel('BER'); xlabel('Filter Order [N]');
 title('BER vs. Filter order');
-ylim([1e-4, 0.1]);
+% ylim([1e-4, 0.1]);
 yline(3.8e-3,'HandleVisibility','off');
 yline(2.2e-4,'HandleVisibility','off');
 

From 72c1cc39f8d30b5ca961c0368ba3e72a9b694a0f Mon Sep 17 00:00:00 2001
From: magf <magf@tf.uni-kiel.de>
Date: Tue, 20 Jan 2026 14:16:02 +0100
Subject: [PATCH 5/7] - Resolved previous conflict

---
 projects/FSO_transmission/first_analysis.m | 293 ---------------------
 1 file changed, 293 deletions(-)
 delete mode 100644 projects/FSO_transmission/first_analysis.m

diff --git a/projects/FSO_transmission/first_analysis.m b/projects/FSO_transmission/first_analysis.m
deleted file mode 100644
index b2a6027..0000000
--- a/projects/FSO_transmission/first_analysis.m
+++ /dev/null
@@ -1,293 +0,0 @@
-
-<<<<<<< HEAD
-base = "C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\FSO Equalizer\FSO_FP_QCL_60umUTC\";
-=======
-base = "C:\Users\Silas\Nextcloud4\Dokumente\02_Ablage_Office\FSO_FP_QCL_60umUTC";
->>>>>>> 719e5508e776c18e3dcf72a954c71bdedda179e0
-mode = 0; %0 oder 1
-M = 2;
-
-all_files = dir(fullfile(base, "**/*.mat"));
-
-if M == 2
-<<<<<<< HEAD
-    tx_data_path = fullfile(base, "14G_PAM2\tx_info\tx_info_PAM2_14Gbd0.75RRC.mat");
-    filename = fullfile(base, "14G_PAM2\M=2_Rs=1.4e10_Fs=8e10_I=265mA_RoP=46.3mW_L=31m_PS=RRC_rolloff=0.75_Mode=Rise.mat");
-elseif M == 4
-    tx_data_path = fullfile(base, "6G_PAM4\tx_info\tx_info_PAM4_6Gbd0.6RRC.mat");
-=======
-    tx_data = load("C:\Users\Silas\Nextcloud4\Dokumente\02_Ablage_Office\FSO_FP_QCL_60umUTC\14G_PAM2\tx_info\tx_info_PAM2_14Gbd0.75RRC.mat");
-    filename = fullfile(base, "14G_PAM2\M=2_Rs=1.4e10_Fs=8e10_I=265mA_RoP=46.3mW_L=31m_PS=RRC_rolloff=0.75_Mode=Rise.mat");
-elseif M == 4
-    tx_data = load("C:\Users\Silas\Nextcloud4\Dokumente\02_Ablage_Office\FSO_FP_QCL_60umUTC\6G_PAM4\tx_info\tx_info_PAM4_6Gbd0.6RRC.mat");
->>>>>>> 719e5508e776c18e3dcf72a954c71bdedda179e0
-    filename = fullfile(base, "6G_PAM4\M=4_Rs=6e9_Fs=8e10_I=255mA_RoP=42.3mW_L=31m_PS=RRC_rolloff=0.6_Mode=Rise.mat");
-end
-
-if mode == 1
-    [f, p] = uigetfile(fullfile(base, "**/*.mat")); 
-    if f~=0
-        filename = fullfile(p,f);
-    end
-end
-
-tx_data = load(tx_data_path);
-datas = load(filename);
-
-%%
-str = filename;
-M_       = str2double(regexp(str, 'M=([^_]+)',       'tokens', 'once'));
-assert(M==M_);
-fsym      = str2double(regexp(str, 'Rs=([^_]+)',      'tokens', 'once'));
-fs      = str2double(regexp(str, 'Fs=([^_]+)',      'tokens', 'once'));
-I       = sscanf(char(regexp(str, 'I=([^_]+)',       'tokens', 'once')), '%f');
-rop     = sscanf(char(regexp(str, 'RoP=([^_]+)',     'tokens', 'once')), '%f');
-L       = sscanf(char(regexp(str, 'L=([^_]+)',       'tokens', 'once')), '%f');
-pulseshape = string(    regexp(str, 'PS=([^_]+)',      'tokens', 'once'));
-rolloff = str2double(regexp(str, 'rolloff=([^_]+)', 'tokens', 'once'));
-mode    = string(    regexp(str, 'Mode=([^\.]+)',   'tokens', 'once'));
-
-%%
-% Tx data
-
-Bits = Informationsignal(tx_data.tx_data,"fs",fsym);
-Symbols = Informationsignal(real(tx_data.tx_PAM_sym),"fs",fsym);
-
-mapping_style = M==4; % Pam2 is like move-it; PAM-4 is different, same mapping like ETH peopled used in Zurich... hence the "eth_style" argument here and there
-PM = PAMmapper(M,0,"eth_style",mapping_style); % one should rename "eth style" as this is simply a different mapping scheme
-
-Symbols_ = PM.map(Bits) .* PM.scaling;
-assert(isequal(Symbols.signal,Symbols_.signal));
-
-Bits_ = PM.demap(Symbols);
-[bits,errors,ber,errorIndice] = calc_ber(Bits_.signal,Bits.signal);
-assert(ber == 0);
-
-%% For comparison, apply pulsef on Tx Symbols
-Pform = Pulseformer("fsym",fsym,"fdac",fs,"pulse","rrc","pulselength",16,"alpha",rolloff);
-Digi_sig_tx_compare = Pform.process(Symbols);
-
-Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff);
-Digi_sync = Pform.process(Symbols);
-
-MF = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff);
-Rx_sig_compare = MF.process(Digi_sig_tx_compare);
-
-%%
-
-% Rx Data
-traceData = datas.tr.lastData(2).trace.ch3;
-
-%FYI: Voltage=(RawData−YReference)×YIncrement+YOrigin
-scoperead_volts = (traceData.RawData - traceData.YReference) * traceData.YIncrement + traceData.YOrigin;
-demystified = isequal(traceData.YData,scoperead_volts);
-assert(demystified);
-
-timesig_compare = [0:1:datas.tr.lastData(1).trace.ch3.Points-1] ./ fs;
-timesig = datas.tr.lastData(1).trace.ch3.XData;
-assert(isequal(traceData.YData,scoperead_volts));
-
-Scope_sig = Electricalsignal(traceData.YData,"fs",fs);
-
-%%
-% 1) matched filter
-% pulse is symmetric, hence we can use pulsef directly as matched filter.
-% It feels off (bit I think correct) that the fsym is now the output freq.!!
-% -> output 2 sps to omit timing recovery!?
-apply_matched_filter = 0;
-k = 4;
-if apply_matched_filter
-    
-    Pform = Pulseformer("fsym",fsym,"fdac",k*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff,"matched",1);
-    Rx_matched = Pform.process(Scope_sig);
-else
-
-    Rx_matched = Filter('filtdegree',4,"f_cutoff",fsym*0.5,"fs",Scope_sig.fs,"filterType",filtertypes.gaussian,"active",true).process(Scope_sig);
-    Rx_matched = Rx_matched.resample("fs_out",k*fsym);
-end
-Rx_matched.spectrum();
-
-<<<<<<< HEAD
-%%
-
-sys = comm.SymbolSynchronizer('TimingErrorDetector', 'Gardner (non-data-aided)', ...
-                              'SamplesPerSymbol', 2, ...
-                              'DampingFactor', 0.7, ...
-                              'NormalizedLoopBandwidth', 0.01);
-Rx_symbolsync = Rx_matched;
-[Rx_symbolsync.signal, timing_error] = sys(Rx_matched.signal);
-
-plot(timing_error); % If this is a ramp, you have drift!
-
-%% timing sync -> at this point we still have no symbol timing recovery, we
-% % try to do this with 2sps EQ!
-
-[~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_symbolsync.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1);
-% Rx_matched_1 = Rx_synced_cell{1};
-% 
-% % Rx_Time_Rec = Rx_matched;
-% [Rx_Time_Rec, Timing_Error] = Timing_Recovery("modulation", 'PAM/PSK/QAM', "timing_error_detector",'Gardner (non-data-aided)','sps',2,'damping_factor',1,'normalized_loop_bandwidth',0.01,'detector_gain',2.7).process(Rx_matched_1);
-% figure;plot(Timing_Error);
-% % Rx_Time_Rec.signal = resample(Rx_Time_Rec.signal, 28e9, 14e9);
-
-%% not working.. 
-% Rx_synced = Rx_Time_Rec;
-Rx_synced = Rx_synced_cell{1};
-len_tr  = 4096*2;
-mu_ffe1 = 0.0001;
-mu_ffe2 = 0.0008;
-mu_ffe3 = 0.001;
-mu_dc = 0.005;
-mu_ffe = [mu_ffe1 mu_ffe3 mu_ffe3];
-mu_dfe = 0.0004;
-duob_mode = db_mode.no_db;
-
-Rx_synced.plot("displayname",'RX: Matched+Sync+2sps','fignum',2);
-
-Digi_sig_compare.normalize("mode","rms").spectrum("displayname",'Tx: RC-shaped','fignum',1,'normalizeTo0dB',0);
-Rx_sig_compare.normalize("mode","rms").spectrum("displayname",'Tx: RC-shaped + matched filtered ','fignum',1,'normalizeTo0dB',0);
-Rx_synced.normalize("mode","rms").spectrum("displayname",'RX: matched filtered + synced','fignum',1,'normalizeTo0dB',0);
-=======
-%% Timing Rec
-apply_timing_rec = 1;
-if apply_timing_rec
-    [Rx_symbolsync, timing_error] = Timing_Recovery("timing_error_detector",'Gardner (non-data-aided)','sps',k,'damping_factor',0.1,'normalized_loop_bandwidth',0.1,'detector_gain',2.7).process(Rx_matched);
-    figure();plot(timing_error);
-    Rx_symbolsync.fs = fsym;
-    % Tsynch
-    [~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_symbolsync.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 0);
-    Rx_synced = Rx_synced_cell{1};
-    sps = 1;
-else
-    % Tsynch
-    [~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_matched.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 0);
-    Rx_synced = Rx_synced_cell{1};
-    Rx_synced = Rx_synced.resample("fs_out",2*fsym);
-    sps = 2;
-end
->>>>>>> 719e5508e776c18e3dcf72a954c71bdedda179e0
-
-if M == 2
-    ber_in_paper = 10^(-2.6); %fig 3a) 4 Gb/s MWIR FSO Transmission using Directly Modulated QCL and an Uncooled UTC-PD at Room-Temperature
-elseif M == 4
-    ber_in_paper = 10^(-2.5);
-end
-
-
-%%
-
-Rx_synced = Rx_synced_cell{1};
-
-
-% -------------------- FFE --------------------
-% requires some more digging what is going on :-)
-eq_ffe      = EQ("Ne",[50, 1, 1],"Nb",[2,0,0], ...
-    "training_length",512,"training_loops",5,"dd_loops",5, ...
-    "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ...
-    "FFEmu",0,"plotfinal",0,"ideal_dfe",1);
-
-vars = logspace(-4,-3,36);
-
-parfor i = 1:numel(vars)
-
-    len_tr  = 4096;
-    mu_ffe1 = 0.01;% mus(i);%0.0001;
-    mu_ffe2 = 0.0008;
-    mu_ffe3 = 0.001;
-    mu_dc = 0.005;
-    mu_ffe = [mu_ffe1 mu_ffe2 mu_ffe3];
-    mu_dfe = vars(i);
-    duob_mode = db_mode.no_db;
-
-    % requires some more digging what is going on :-)
-    eq_ffe_1      = EQ("Ne",[150, 1, 0],"Nb",[50,0,0], ...
-        "training_length",len_tr,"training_loops",5,"dd_loops",5, ...
-        "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ...
-        "FFEmu",0,"plotfinal",0,"ideal_dfe",0);
-
-    eq_ffe_2 = FFE("epochs_tr",1,"epochs_dd",vars(i),"len_tr",4096,"mu_dd",vars(i),"mu_tr",vars(i),"order",999,"sps",1,"decide",0, "adaption",adaption_method.nlms,"dd_mode",0);
-    % eq_ffe_2 = FFE_DFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"ffe_mu_dd",1e-5,"dfe_mu_dd",mus(i),"ffe_mu_tr",0,"dfe_mu_tr",0,"ffe_order",50,"dfe_order",10,"sps",1,"decide",1);
-
-
-    ffe_results = ffe(eq_ffe_1,M,Rx_synced,Symbols,Bits, ...
-        "precode_mode",duob_mode,'showAnalysis',1,"postFFE",[], ...
-        "eth_style_symbol_mapping",mapping_style);
-
-    ffe_results.metrics.BER
-    bers(i) = ffe_results.metrics.BER;
-end
-
-figure();
-plot(vars,bers);
-yline(ber_in_paper)
-beautifyBERplot();
-% 
-fprintf('Paper: %.1e \n \n',ber_in_paper);
-ffe_results.metrics.print("description",'FFE');
-fprintf('FFE: %.1e \n',ffe_results.metrics.BER);
-
-
-%% -------------------- VNLE + MLSE --------------------
-len_tr  = 4096;
-mu_ffe1 = 0.0001;% mus(i);%0.0001;
-mu_ffe2 = 0.0008;
-mu_ffe3 = 0.001;
-mu_dc = 0.005;
-mu_ffe = [mu_ffe1 mu_ffe2 mu_ffe3];
-mu_dfe = 0.0004;
-duob_mode = db_mode.no_db;
-
-pf_ncoeffs  = 4;
-
-vars = 1:7;
-bers = zeros(size(vars));
-parfor i = 1:numel(vars)
-    
-    eqv  = EQ("Ne",[200, 1, 0],"Nb",[2, 0, 0], ...
-    "training_length",len_tr,"training_loops",5,"dd_loops",5, ...
-    "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ...
-    "FFEmu",0,"plotfinal",0,"ideal_dfe",1);
-    pf_ncoeffs  = vars(i);
-    pf      = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1);
-    mlse_   = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0,"eth_style",mapping_style).levels);
-    [vnle_results, mlse_results] = vnle_postfilter_mlse(eqv, pf, mlse_, M, Rx_synced, Symbols, Bits, ...
-        "precode_mode", duob_mode, 'showAnalysis', 1, "postFFE", [], "eth_style_symbol_mapping", mapping_style);
-    fprintf('Paper: %.1e \n \n',ber_in_paper);
-    vnle_results.metrics.print("description",'VNLE');
-    mlse_results.metrics.print("description",'MLSE');
-    bers(i) = mlse_results.metrics.BER;
-end
-%%
-figure();hold on
-plot(vars,bers_ffe,'DisplayName','FFE [200,0,0] + PF + MLSE');
-plot(vars,bers_vnle,'DisplayName','VNLE [200,1,0] + PF + MLSE');
-plot(vars,bers_vnledfe,'DisplayName','VNLE [200,1,0] + DFE [2] + PF + MLSE');
-plot(vars,bers_vnledfe_ideal,'DisplayName','VNLE [200,1,0] + ideal DFE [2] + PF + MLSE');
-yline(ber_in_paper);
-yline([2e-2, 4.85e-3, 3.8e-3, 2,2e-4],'LineWidth',2,'Color',[0.8,0.8,0.8],'LineStyle',':','HandleVisibility','off');
-ylim([1e-5,0.1]);
-beautifyBERplot();
-
-%% -------------------- DB target --------------------
-mlse_db_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",M,'trellis_states',PAMmapper(M,0).levels);
-
-eq_       = EQ("Ne",[50, 5, 5],"Nb",[0,0,0],"training_length",len_tr,"training_loops",5,"dd_loops",5, ...
-    "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1);
-
-dbt_results = duobinary_target(eq_,mlse_db_, M, Rx_synced, Symbols, Bits, ...
-    "precode_mode", duob_mode, 'showAnalysis', 0, "postFFE", [],"eth_style_symbol_mapping",mapping_style);
-
-fprintf('Paper: %.1e \n \n',ber_in_paper);
-dbt_results.metrics.print("description",'Duobinary');
-
-
-%%
-
-%ML-based MLSE (L=2)
-mu_ml = 0.01; training_epochs = 100;
-ml_mlse_equalizer = ML_MLSE("epochs_tr",training_epochs,"epochs_dd",1, ...
-    "len_tr",len_tr,"mu_dd",mu_ml,"mu_tr",mu_ml,"order",11,"sps",1, ...
-    "traceback_depth",128,"L",3,"delta",4,"adaptive_mu",0);
-
-[ml_mlse_results] = ml_mlse(ml_mlse_equalizer, M, Rx_synced, Symbols, Bits,"precode_mode",duob_mode);
-ml_mlse_results.metrics.print("description",'ML ');
\ No newline at end of file

From 3886524b0e27fdc128e1df3baaa919aca6b87295 Mon Sep 17 00:00:00 2001
From: magf <magf@tf.uni-kiel.de>
Date: Tue, 20 Jan 2026 14:49:36 +0100
Subject: [PATCH 6/7] - Resolved conflict (hopefully)

---
 Classes/04_DSP/Timing_Recovery.m  | 14 --------------
 Functions/EQ_structures/ml_mlse.m |  4 ----
 2 files changed, 18 deletions(-)

diff --git a/Classes/04_DSP/Timing_Recovery.m b/Classes/04_DSP/Timing_Recovery.m
index 5e33b6b..3999ad7 100644
--- a/Classes/04_DSP/Timing_Recovery.m
+++ b/Classes/04_DSP/Timing_Recovery.m
@@ -27,19 +27,11 @@ classdef Timing_Recovery < handle
                 obj.(fn{n}) = options.(fn{n});
             end
 
-<<<<<<< HEAD
             %obj-Initialization here%
 
         end
 
         function [data_out, timing_error] = process(obj, data_in)
-=======
- 
-
-        end
-
-        function [data_out,timing_error] = process(obj, data_in)
->>>>>>> 719e5508e776c18e3dcf72a954c71bdedda179e0
 
             timing_synchronization = comm.SymbolSynchronizer( ...
                 "Modulation",          obj.modulation,...
@@ -48,14 +40,8 @@ classdef Timing_Recovery < handle
                 "DampingFactor",       obj.damping_factor, ...
                 "NormalizedLoopBandwidth", obj.normalized_loop_bandwidth, ...
                 "DetectorGain",        obj.detector_gain);
-<<<<<<< HEAD
             data_out = data_in;
             [data_out.signal, timing_error] = timing_synchronization(data_in.signal);
-=======
-
-            data_out = data_in;
-            [data_out.signal,timing_error] = timing_synchronization(data_in.signal);
->>>>>>> 719e5508e776c18e3dcf72a954c71bdedda179e0
 
         end
     end
diff --git a/Functions/EQ_structures/ml_mlse.m b/Functions/EQ_structures/ml_mlse.m
index 17e25ca..84ff5cb 100644
--- a/Functions/EQ_structures/ml_mlse.m
+++ b/Functions/EQ_structures/ml_mlse.m
@@ -45,7 +45,6 @@ end
 
 ml_mlse_results.config = Equalizerstruct();
 
-<<<<<<< HEAD
 eq_small = strip_eq(eq_, 10);
 fn = fieldnames(eq_small);
 for k = 1:numel(fn)
@@ -54,9 +53,6 @@ for k = 1:numel(fn)
     end
 end
 json_str = jsonencode(eq_small);
-=======
-eq_ = strip_eq(eq_, 10);
->>>>>>> 719e5508e776c18e3dcf72a954c71bdedda179e0
 
 fn = fieldnames(eq_);
 for k = 1:numel(fn)

From 7eb33648146b49212f6a8e0b67972ae7c96f2a37 Mon Sep 17 00:00:00 2001
From: Silas Oettinghaus <silas.oettinghaus@email.uni-kiel.de>
Date: Thu, 29 Jan 2026 16:49:50 +0100
Subject: [PATCH 7/7] Minimal changes Add Theory plots for Silas Diss

---
 Classes/04_DSP/Equalizer/ML_MLSE.m            | 515 +++++++++++++++++-
 .../Dissertation/mach_zehnder_modulator.m     | 140 +++++
 Functions/mat2tikz_improved.m                 |  23 +
 Libs/wesanderson_colors/WesPalette.m          | 127 +++--
 .../minimal_example_wespalette.m              |   4 +-
 .../final/Copy_of_FIGURE_WAVELENGTH.m         |  26 +-
 .../final/FIGURE_introduction.m               | 132 +++--
 .../ML_based_MLSE/analyze_filter_length.m     |   8 +-
 8 files changed, 860 insertions(+), 115 deletions(-)
 create mode 100644 Functions/Theory/Dissertation/mach_zehnder_modulator.m
 create mode 100644 Functions/mat2tikz_improved.m

diff --git a/Classes/04_DSP/Equalizer/ML_MLSE.m b/Classes/04_DSP/Equalizer/ML_MLSE.m
index 3f584c3..168b837 100644
--- a/Classes/04_DSP/Equalizer/ML_MLSE.m
+++ b/Classes/04_DSP/Equalizer/ML_MLSE.m
@@ -1,3 +1,499 @@
+% classdef ML_MLSE < handle
+%     % ALGORITHM DESCRIBED IN:
+%     % W. Lanneer and Y. Lefevre, “Machine Learning-Based Pre-Equalizers for
+%     % Maximum Likelihood Sequence Estimation in High-Speed PONs,”
+%     % in 2023 31st European Signal Processing Conference
+% 
+%     % Further ML Refs:
+%     % https://machinelearningmastery.com/cross-entropy-for-machine-learning/
+%     % https://docs.pytorch.org/docs/stable/generated/torch.nn.CrossEntropyLoss.html
+% 
+%     % The central idea is to overcome the (white-) noise assumption within the previously described
+%     % Viterbi algorithm, more precisely a closed-loop optimization is proposed that finds a suitable
+%     % filter-set to directly compute the branch metrics c_k (s,s^' ). These can directly be used to
+%     % carry out the conventional Viterbi algorithm. The system consists of S^L S=F linear FIR filters,
+%     % combined with one bias coefficient respectively. These filters take the received input samples to
+%     % compute the branch metrics estimates (c_k ) ̂(s,s^' ) according toThe central idea is to overcome
+%     % the (white-) noise assumption within the previously described Viterbi algorithm, more precisely
+%     % a closed-loop optimization is proposed that finds a suitable filter-set to directly compute the
+%     % branch metrics c_k (s,s^' ). These can directly be used to carry out the conventional Viterbi
+%     % algorithm. The system consists of S^L S=F linear FIR filters, combined with one bias coefficient
+%     % respectively. These filters take the received input samples to compute the branch metrics
+%     % estimates. Finally, the usual Viterbi is carried out...
+% 
+%     % Recommended Settings and some findings:
+% 
+%     % Requires many training epochs. According to ML people, 100,200 or
+%     % even up to 1000 epochs are normal for ML-convergence
+% 
+%     % The mu parameter _can_ be adaptive - using the cross entropy and when
+%     % analyzing the isolated training it looks very promisig. However, is
+%     % later use I found this is not as stable as a fixed learning rate.
+%     % mu = 0.1 worked good for me
+% 
+%     % Longer orders/ filter length are not always better. For me order=11
+%     % was good.
+% 
+%     % Delay factor (delta) is good when the order is also increased. With
+%     % order = 11, a delta of =4 shows good results
+% 
+%     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
+% 
+%         adaptive_mu
+% 
+%         constellation
+% 
+%         L %viterbi memory length
+% 
+%         alpha
+%         DIR
+%         DIR_flip
+%         trellis_states
+% 
+%         traceback_depth
+% 
+%         % --- Added internal class variables used later ---
+%         S
+%         Nf
+%         delta
+%         nStates
+%         nFeasible
+%         combs
+%         first_sym
+%         last_sym
+%         valid
+%         valid_to_idx
+%         valid_from_idx
+%         w
+% 
+%         % --- New: fast state lookup ---
+%         true_to_state_idx
+%         state_dict              % containers.Map: key(sequence)->state index
+%         key_fmt = '%.8g_';      % key format for sequence strings
+%         nSym                      % |constellation|
+% 
+%         ber = []
+%         ce  = ones(1,1);
+%     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.adaptive_mu = 1;
+% 
+%                 options.delta = 0;
+%                 options.traceback_depth = 1024;
+% 
+%                 options.L = 1
+% 
+%             end
+% 
+%             fn = fieldnames(options);
+%             for n = 1:numel(fn)
+%                 obj.(fn{n}) = options.(fn{n});
+%             end
+% 
+%             obj.e = zeros(obj.order,1);
+%             obj.error = 0;
+%         end
+% 
+%         function [X,X_viterbi] = process(obj, X, D)
+% 
+%             % actual processing of the signal (steps 1. - 3.)
+%             % 1 normalize RMS
+%             X = X.normalize("mode","rms");
+% 
+%             % 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(obj.constellation);      % alphabet size
+%             obj.Nf = obj.order*obj.sps;            % filter length
+%             % 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(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{:}).');  % rows: states, columns: [x_k, x_{k-1}, ...]
+%             obj.first_sym = obj.combs(:,1);
+%             obj.last_sym  = obj.combs(:,end);
+%             obj.nStates   = size(obj.combs,1);
+% 
+%             % --- Valid transitions
+%             obj.valid = false(obj.nStates);
+%             for from = 1:obj.nStates
+%                 for to = 1:obj.nStates
+%                     if all(obj.combs(to,2:end) == obj.combs(from,1:end-1))
+%                         obj.valid(to,from) = true;
+%                     end
+%                 end
+%             end
+%             [obj.valid_to_idx, obj.valid_from_idx] = find(obj.valid);
+% 
+%             % --- Allocate vectors and weights
+%             % !! IF SHAPE FIT, then we already have smth there an we want
+%             % to start with the existing fitler-set
+%             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
+% 
+%             % --- 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
+%             % ==============================================================
+% 
+%             % Training Mode
+%             n = obj.len_tr;
+%             training = 1;
+%             obj.equalize(X.signal, D.signal,obj.mu_tr,obj.epochs_tr,n,training);
+%             obj.e_tr = obj.e;
+% 
+%             % ==============================================================
+%             %  DD-Mode / Fixed Mode
+%             % ==============================================================
+% 
+%             % 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_ref] = equalize(obj,x,d,mu,epochs,N,training)
+%             % ==============================================================
+%             %  FFE + Whitening + ML-Based Branch Metric Estimation + Viterbi
+%             % ==============================================================
+%             debug = 1;
+%             showPlots = 1;
+% 
+%             % --- 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);  % v_{k-1}(s′)
+%                 c_hat     = zeros(1,obj.nFeasible);
+%                 v_tilde   = zeros(1,obj.nFeasible);
+%                 pred      = zeros(nSymbols, obj.nStates, 'uint32');
+%                 pm_sto    = nan(obj.nStates, nSymbols,'like',pm);
+%                 CE_accum = 0;
+% 
+% 
+%                 %%% START IDX 
+%                 if training
+%                     max_start = length(x) - ( (ceil(N/obj.sps)-1)*obj.sps + 1 );
+%                     max_start = max(1, max_start);  % safety
+%                     start_sample = randi([1, max_start], 1); %rnd training; not really good
+%                     start_sample = 1;
+%                     end_sample   = start_sample + (ceil(N/obj.sps)-1)*obj.sps;
+%                 else
+%                     start_sample = 1;%obj.len_tr;
+%                     end_sample   = N;
+%                 end
+% 
+%                 start_symbol = 1 + floor((start_sample - 1)/obj.sps);   % ABSOLUTE symbol index
+% 
+%                 if numel(d) >= obj.L && start_symbol >= obj.L
+%                     init_seq = d(start_symbol-obj.L+1 : start_symbol);      % [d_k-L+1 ... d_k]
+%                     true_to_state_idx = obj.state_dict(obj.seq_key(flip(init_seq))); % [d_k ... d_k-L+1]
+%                 else
+%                     % Not enough history – fall back to state 1
+%                     true_to_state_idx = uint32(1);
+%                 end
+% 
+%                 symbol = 0;
+%                 for sample = start_sample:obj.sps:end_sample
+%                     symbol = symbol + 1;
+%                     k = symbol;
+%                     sym_idx = start_symbol + (symbol - 1);
+% 
+%                     % --- Build Δ-delayed observation window y_k
+%                     i1  = sample - obj.Nf + 1 + obj.delta;
+%                     i2  = sample + obj.delta;
+%                     buf = x(max(1,i1):min(length(x),i2));
+%                     padL = max(0,1 - i1);
+%                     padR = max(0,i2 - length(x));
+%                     yk  = [zeros(padL,1); buf(:); zeros(padR,1)];   % Nf×1
+%                     yk  = [yk;1];
+% 
+%                     % --- 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: 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 1 %training
+%                         % --- allocate storage once
+%                         if epoch == 1 && symbol == 1
+%                             obj.true_to_state_idx = ones(ceil(N/obj.sps),1,'uint32');
+%                         end
+% 
+%                         % --- previous "to" becomes current "from"
+%                         if symbol > 1
+%                             true_from_state_idx = obj.true_to_state_idx(symbol-1);
+%                         else
+%                             true_from_state_idx = 1;
+%                         end
+% 
+%                         % --- compute or reuse "to" state
+%                         if epoch == 1
+%                             % only compute in first epoch
+%                             if sym_idx >= obj.L
+%                                 key_to = obj.seq_key(flip(d(sym_idx-obj.L+1 : sym_idx)));
+%                                 if isKey(obj.state_dict, key_to)
+%                                     obj.true_to_state_idx(symbol) = obj.state_dict(key_to);
+%                                 else
+%                                     obj.true_to_state_idx(symbol) = true_from_state_idx;
+%                                 end
+%                             else
+%                                 obj.true_to_state_idx(symbol) = true_from_state_idx;
+%                             end
+%                         end
+% 
+%                         % --- reuse cached state from second epoch onward
+%                         true_to_state_idx = obj.true_to_state_idx(symbol);
+% 
+%                         % --- ensure valid (from,to)
+%                         dirac = zeros(obj.nFeasible,1);
+%                         mask  = obj.valid_from_idx==true_from_state_idx & ...
+%                             obj.valid_to_idx  ==true_to_state_idx;
+%                         if any(mask)
+%                             dirac(mask) = 1;
+%                         else
+%                             idx = find(obj.valid_from_idx==true_from_state_idx,1,'first');
+%                             dirac(idx) = 1;
+%                             obj.true_to_state_idx(symbol) = obj.valid_to_idx(idx);
+%                         end
+% 
+% 
+% 
+% 
+%                         % softmax over -v_tilde (numerically safe shift)
+%                         v_shift = -(v_tilde - min(v_tilde));         % shift to small positive numbers
+%                         v_shift = min(v_shift, 100);                  % clamp exponent argument (≈ exp(50)=3e21)
+%                         expv = exp(v_shift);
+%                         p = expv ./ (sum(expv) + eps);
+% 
+%                         % for logging only:
+%                         CE_symbol(symbol) = -log(p(dirac==1) + eps);
+% 
+%                         if sym_idx > obj.L
+%                             CE_smooth(symbol) = 0.01*CE_symbol(symbol) + 0.99*CE_smooth(symbol-1);
+%                         else
+%                             if epoch > 1
+%                                 CE_smooth(symbol) = obj.ce(end); %use ce from last epoch or =1 for very first round?!
+%                             else
+%                                 CE_smooth(symbol) = CE_symbol(symbol);
+%                             end
+%                         end
+% 
+%                         CE_accum = CE_symbol(symbol) + CE_accum;
+% 
+% 
+%                         % gradient term (t - p)
+%                         dmp = (dirac - p)';  % 1×nFeasible
+% 
+%                         % Per-feature gradient; implicit expansion gives (Nf+1)×nFeasible
+%                         dL_Dw = (yk) .* dmp;
+% 
+%                         % Start updates only when the ABSOLUTE symbol index has ≥ L history
+%                         if sym_idx >= obj.L
+%                             if obj.adaptive_mu
+%                                 mu_eff = CE_smooth(sym_idx);
+%                                 mu_eff = max(min(mu_eff, 0.2), 1e-4);
+%                             else
+%                                 mu_eff = mu;
+%                             end
+% 
+%                             obj.w = obj.w - mu_eff .* dL_Dw;   % (Nf+1)×nFeasible
+%                         end
+% 
+%                         % if debug && epoch > 2
+%                         %     figure(100);
+%                         %     subplot(4,1,1);
+%                         %     heatmap(p');
+%                         %     title('Probs')
+%                         %     subplot(4,1,2);
+%                         %     heatmap(dmp);
+%                         %     title('Update')
+%                         %     subplot(4,1,3);
+%                         %     heatmap(dL_Dw);
+%                         %     title('Update')
+%                         %     subplot(4,1,4);
+%                         %     heatmap(bj.w);
+%                         %     title('Update')
+%                         % 
+%                         % end
+% 
+%                     end
+% 
+% 
+% 
+%                     % --- Compare-Select (matrix form, min of costs)
+%                     v_tilde_mat = inf(obj.nStates, obj.nStates);
+%                     v_tilde_mat(obj.valid) = v_tilde;
+%                     [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
+% 
+%                 % --- 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
+%                     viterbi_path(n-1) = pred(n, viterbi_path(n));
+%                 end
+% 
+%                 y_ref = d(start_symbol:end);
+%                 y     = obj.first_sym(viterbi_path);
+% 
+%                 if debug && training
+%                     sym_start = start_symbol;
+%                     sym_end   = start_symbol + symbol - 1;
+%                     ref_slice = d(sym_start : sym_end);
+%                     err = sum(y ~= ref_slice(1:numel(y)));
+% 
+%                     try
+%                         ref_bits = PAMmapper(obj.S,0).demap(ref_slice);
+%                         eq_bits = PAMmapper(obj.S,0).demap(y);
+%                         [~, ~, ber, ~] = calc_ber(ref_bits, eq_bits, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1);
+%                         fprintf('Epoch: %d - BER: %.1e  \n',epoch, ber);
+%                         obj.ber(epoch) = ber;
+%                     catch
+%                         ser = err./length(y);
+%                         fprintf('Epoch: %d - SER: %.1e  \n',epoch, ser);
+%                     end
+% 
+%                     obj.ce(epoch) = CE_accum./symbol;
+% 
+%                     if showPlots
+%                         figure(10);clf
+%                         subplot(3,2,1:2);
+%                         heatmap(obj.w);
+%                         title('Filter')
+% 
+%                         subplot(3,2,3);
+%                         v_tildemat = NaN(obj.nStates, obj.nStates);
+%                         v_tildemat(obj.valid) = v_tilde;      % log-domain scores
+%                         heatmap(v_tildemat);
+%                         title('Path Metrics (v_tilde)')
+% 
+%                         subplot(3,2,4);
+%                         scatter(1:symbol,pm_sto,1,'.')
+%                         title('Path Metric Winners')
+% 
+%                         subplot(3,2,5);hold on
+%                         scatter(1:symbol,CE_symbol,1,'.');
+%                          scatter(1:symbol,CE_smooth,1,'.')
+%                         title('Cross Entropy')
+% 
+%                         subplot(3,2,6); hold on
+% 
+%                         % Left y-axis: Cross Entropy (linear)
+%                         yyaxis left
+%                         scatter(1:length(obj.ce), obj.ce, 10, 's', 'filled')
+%                         ylabel('Cross Entropy')
+% 
+%                         % Right y-axis: BER (logarithmic)
+%                         yyaxis right
+%                         scatter(1:length(obj.ber), obj.ber, 10, 'd', 'filled')
+%                         set(gca, 'YScale', 'log')
+%                         ylabel('BER (log scale)')
+% 
+%                         xlim([1, epochs])
+%                         xlabel('Epoch')
+%                         title('Cross Entropy // BER')
+%                         grid on
+% 
+%                         drawnow
+%                     end
+%                 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
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
 classdef ML_MLSE < handle
     % ---------------------------------------------------------------------
     % W. Lanneer and Y. Lefevre,
@@ -101,7 +597,7 @@ classdef ML_MLSE < handle
             obj.S = obj.nSym;
             obj.Nf = obj.order * obj.sps;
             obj.nStates = obj.S^obj.L;
-            obj.nFeasible = obj.nStates * obj.S;
+            obj.nFeasible = obj.nStates * obj.S; %feasible state transitions
 
             % --- Trellis mapping
             obj.trellis_states = reshape(obj.constellation,1,[]);
@@ -164,8 +660,8 @@ classdef ML_MLSE < handle
         % EQUALIZE
         % ==============================================================
         function [y,y_ref] = equalize(obj,x,d,mu,epochs,N,training)
-            debug = 0;
-            showPlots = 0;
+            debug = 1;
+            showPlots = 1;
             y = zeros(N,1);
             nSymbols = ceil(N/obj.sps);
 
@@ -241,6 +737,19 @@ classdef ML_MLSE < handle
                         dirac(trans_idx)=1;
                     end
 
+                    % --- ensure valid (from,to)
+                    if ~any(dirac)
+                        mask  = obj.valid_from_idx==true_from_state_idx & ...
+                            obj.valid_to_idx  ==true_to_state_idx;
+                        if any(mask)
+                            dirac(mask) = 1;
+                        else
+                            idx = find(obj.valid_from_idx==true_from_state_idx,1,'first');
+                            dirac(idx) = 1;
+                            obj.true_to_state_idx(symbol) = obj.valid_to_idx(idx);
+                        end
+                    end
+
                     % ===================================================================
                     % TRAINING MODE (weight update)
                     % ===================================================================
diff --git a/Functions/Theory/Dissertation/mach_zehnder_modulator.m b/Functions/Theory/Dissertation/mach_zehnder_modulator.m
new file mode 100644
index 0000000..7b379c2
--- /dev/null
+++ b/Functions/Theory/Dissertation/mach_zehnder_modulator.m
@@ -0,0 +1,140 @@
+% Minimal MZM transfer-function demo (sinusoidal drive) — aligned with your notation
+%
+% Implements exactly:
+%   E_out(t) = E0 * exp(j*w0*t) * exp(-j*w0*L*n_eff/c0) * 1/2 * [ exp(-j*phi1(t)) + rho*exp(-j*phi2(t)) ]
+% with  phi_{1,2}(t) = pi * v_{1,2}(t)/Vpi
+%
+% Push-pull:
+%   v1(t) = +v_drive(t)/2 ,  v2(t) = -v_drive(t)/2   =>  phi1 = +pi/2 * v_drive/Vpi, phi2 = -pi/2 * v_drive/Vpi
+%
+% And the ideal TF (rho=1):
+%   E_out/E_in = exp(-j*w0*L*n_eff/c0) * cos( (pi/2) * v_drive/Vpi )
+%
+% Note: E_in(t) = E0 * exp(j*w0*t) in this script.
+
+% Parameters
+c0    = physconst('lightspeed');      % [m/s]
+lambda0 = 1310e-9;      % [m]
+omega0  = 2*pi*c0/lambda0;
+
+L     = 5e-3;           % [m] effective phase section length (set as needed)
+n_eff = 2.2;            % [-] effective index (set as needed)
+
+E0    = 1;              % field amplitude (arbitrary)
+Vpi   = 3.2;            % [V] half-wave voltage (your V_pi)
+
+% Drive
+f0    = 1e9;            % [Hz]
+fs    = 100e9;          % [Hz]
+Nper  = 1;              % number of periods
+Vpp   = 0.5*Vpi;        % [V] peak-to-peak of v_drive(t)
+
+biasV = 2;            % [V] differential bias added to v_drive
+
+% Analytic
+v_ = linspace(-1,2, 2001);
+% Field transfer function (amplitude)
+Field_mzm_analytic = cos((pi/2)*v_);
+% Power transfer function (intensity)
+P_mzm_analytic = Field_mzm_analytic.^2;
+
+% Imbalance factor in YOUR notation:
+rho   = 1;            % rho=1 -> ideal balanced MZM (collapses to ideal TF)
+
+% Time axis + differential drive voltage v_drive(t)
+T  = Nper/f0;
+t  = (0:1/fs:T-1/fs).';
+
+v_drive = biasV + (Vpp/2)*sin(2*pi*f0*t);  % v_drive(t) (peak = Vpp/2)
+
+% Push-pull branch voltages (consistent with v_drive = v1 - v2)
+v1 = +0.5*v_drive;   % arm 1
+v2 = -0.5*v_drive;   % arm 2
+
+% Phases phi1, phi2
+phi1 = pi * v1 / Vpi;
+phi2 = pi * v2 / Vpi;
+
+% Fields: E_in and E_out (exactly your Eq. (mzm_e_field))
+E_in  = E0 .* exp(1i*omega0*t);
+
+common_phase = exp(-1i * (omega0*L*n_eff/c0));  % exp(-j*omega0*L*n_eff/c0)
+
+E_out = E0 .* exp(1i*omega0*t) .* common_phase .* 0.5 .* ...
+        ( exp(-1i*phi1) + rho .* exp(-1i*phi2) );
+
+% Transfer function (numerical): E_out/E_in
+H_num = E_out ./ E_in;
+
+% Power (normalized)
+Pnorm_num   = abs(H_num).^2;     % since |E_out/E_in|^2
+
+
+
+% Ideal TF (analytic) for comparison (rho=1, push-pull)
+H_ideal = common_phase .* cos( (pi/2) * (v_drive./Vpi) );
+
+Pnorm_ideal = abs(H_ideal).^2;
+Pnorm_math = cos( (pi/2) * (v_drive./Vpi) ).^2;
+
+
+
+
+
+set(groot, 'defaultLegendInterpreter', 'tex');
+set(groot, 'defaultAxesTickLabelInterpreter', 'tex');
+set(groot, 'defaultTextInterpreter', 'tex');
+
+% Normalized voltage axis (multiples of Vpi)
+v_norm = v_drive./Vpi;
+
+colfield = [0,0,0]; %is black
+colpow = linspecer(2);
+colpow = colpow(1,:);
+colvdrive = linspecer(2);
+colvdrive = colvdrive(2,:);
+
+%% SIGNAL IN
+figure(1); clf
+plot(v_norm,t*1e9, 'LineWidth', 1.0,'Color',colvdrive); grid on;
+ylabel('t [ns]'); xlabel('v_{drive}(t)/V_\pi'); 
+title('Drive voltage (normalized)');
+xlim([min(v_) max(v_)]);
+
+%% IN/OUT (static transfer) — normalized x-axis + analytic curve
+figure(2); clf
+plot(v_, Field_mzm_analytic, 'LineWidth', 1.2,'LineStyle','--','Color',colfield); hold on;% analytic power TF
+plot(v_, P_mzm_analytic, 'LineWidth', 1.2, 'Color',colpow); hold on;% analytic power TF
+% show input time signal
+plot(v_norm,-1+t*1e9, 'LineWidth', 1.0,'Color',colvdrive); grid on;
+% show output time signal
+plot(2+t*1e9, Pnorm_num,   'LineWidth', 1.0,'DisplayName','Intensity', 'Color',colvdrive); hold on;
+plot(2+t*1e9, real(H_ideal), '--', 'LineWidth', 1.0,'DisplayName','Field','Color',colfield); hold on;
+scatter(v_norm, Pnorm_num, 12, '.', 'LineWidth', 1,'MarkerEdgeColor',colvdrive); 
+scatter(biasV./Vpi,(cos((pi/2)*biasV./Vpi)^2),10,'Marker','o');
+line([min(v_drive), min(v_drive)]./Vpi,[(cos((pi/2)*min(v_drive)./Vpi)^2), -2],'linewidth',0.5,'color','black','linestyle','--');
+line([max(v_drive) max(v_drive)]./Vpi,[(cos((pi/2)*max(v_drive)./Vpi)^2), -2],'linewidth',0.5,'color','black','linestyle','--');
+xline([min(v_norm) max(v_norm)])
+
+grid on;
+xlabel('v_{drive}(t)/V_\pi'); ylabel('|E_{out}/E_{in}|^2');
+% legend
+xlim([min(v_) max(v_)+1]);
+ylim([-1 1]);
+
+% mat2tikz_improved('C:\Users\Silas\Documents\6971e0b65b380ca6d71c837f\02_IMDD_System\tikz\mzm.tex');
+
+
+%%
+% % FIELD TF (only field here; do not mix power into this figure)
+figure(3); clf
+% plot(t*1e9, real(H_num), 'LineWidth', 1.0); hold on;
+% plot(t*1e9, real(H_ideal), '--', 'LineWidth', 1.0,'DisplayName','Field','Color',colfield); hold on;
+plot(t*1e9, Pnorm_num,   'LineWidth', 1.0,'DisplayName','Intensity', 'Color',colpow); hold on;
+grid on;
+xlabel('t [ns]'); ylabel('Re\{E_{out}/E_{in}\}');
+legend
+mat2tikz_improved('C:\Users\Silas\Documents\6971e0b65b380ca6d71c837f\02_IMDD_System\tikz\mzm_out.tex');
+
+
+
diff --git a/Functions/mat2tikz_improved.m b/Functions/mat2tikz_improved.m
new file mode 100644
index 0000000..a36b1d1
--- /dev/null
+++ b/Functions/mat2tikz_improved.m
@@ -0,0 +1,23 @@
+function mat2tikz_improved(filename)
+arguments
+    % Default to the path in your example if no argument is provided
+    filename (1,1) string = 'C:\Users\Silas\Documents\Dissertation\00_Examples\tikz\textfig.tikz';
+end
+cleanfigure;
+matlab2tikz(char(filename), ...
+    'width','\fwidth', ...
+    'height','\fheight', ...
+    'showInfo',false, ...
+    'extraAxisOptions',{ ...
+        'legend style={font=\footnotesize}', ...
+        'xlabel style={font=\color{white!15!black},font=\small},',...
+        'ylabel style={font=\color{white!15!black},font=\small},',...
+        'legend columns=1', ...
+        'every axis/.append style={font=\scriptsize}',...
+    'legend columns=1',...
+    'legend style={at={(0.02,0.98)},font=\footnotesize,draw=black!60,rounded corners=2pt,inner sep=1pt,fill=white,column sep=6pt,anchor= north west}',...
+    'legend style={at={(0.02,0.98)},draw=white!0!white,font=\scriptsize,inner sep=0.1pt,fill=white,column sep=1pt,anchor= north west}',...
+     'every axis/.append style={font=\scriptsize}',...
+    });
+
+end
\ No newline at end of file
diff --git a/Libs/wesanderson_colors/WesPalette.m b/Libs/wesanderson_colors/WesPalette.m
index 7cd1dcb..13e976f 100644
--- a/Libs/wesanderson_colors/WesPalette.m
+++ b/Libs/wesanderson_colors/WesPalette.m
@@ -1,10 +1,13 @@
 classdef WesPalette
     % WESPALETTE  Wes Anderson color palettes with auto-completion
     % Usage:
-    %   cmap = WesPalette.Zissou1.rgb()
-    %   cmap = WesPalette.Zissou1.rgb(3)
-
-    % https://github.com/karthik/wesanderson?tab=readme-ov-file
+    %   cmap = WesPalette.Zissou1.rgb()                     % full palette
+    %   cmap = WesPalette.Zissou1.rgb(3)                    % 3 colors (discrete default)
+    %   cmap = WesPalette.Zissou1.rgb(12,"discrete")        % any n, no interpolation
+    %   cmap = WesPalette.Zissou1.rgb(256,"continuous")     % smooth colormap (Lab interpolation)
+    %
+    % Requires:
+    %   - colorspace.m (Pascal Getreuer) on MATLAB path for "continuous" mode
 
     enumeration
         BottleRocket1
@@ -34,21 +37,33 @@ classdef WesPalette
     end
 
     methods
-        function cmap = rgb(obj, n)
+        function cmap = rgb(obj, n, mode)
             % Return palette as Nx3 RGB colormap [0–1]
+            %
+            % n    : number of requested colors (optional)
+            % mode : "discrete" (default) or "continuous"
 
-            hex = obj.hex();
+            base_hex = obj.hex();
+            base_rgb = WesPalette.hex2rgb(base_hex);
 
-            rgb = hex2rgb(hex);
+            if nargin < 2 || isempty(n)
+                cmap = base_rgb;
+                return;
+            end
+            if nargin < 3 || isempty(mode)
+                mode = "discrete";
+            end
+            mode = lower(string(mode));
 
-            if nargin == 2
-                if n > size(rgb,1)
-                    error('Requested %d colors, but only %d available.', ...
-                          n, size(rgb,1))
-                end
-                cmap = rgb(1:n,:);
+            validateattributes(n, {'numeric'}, {'scalar','integer','positive'}, mfilename, 'n');
+            if mode ~= "discrete" && mode ~= "continuous"
+                error('mode must be "discrete" or "continuous".');
+            end
+
+            if mode == "discrete"
+                cmap = WesPalette.sample_discrete(base_rgb, n);
             else
-                cmap = rgb;
+                cmap = WesPalette.interpolate_continuous_lab(base_rgb, n);
             end
         end
     end
@@ -56,78 +71,116 @@ classdef WesPalette
     methods (Access = private)
         function hex = hex(obj)
             % Internal HEX storage
-
             switch obj
                 case WesPalette.BottleRocket1
                     hex = {'#A42820','#5F5647','#9B110E','#3F5151','#4E2A1E','#550307','#0C1707'};
-
                 case WesPalette.BottleRocket2
                     hex = {'#FAD510','#CB2314','#273046','#354823','#1E1E1E'};
-
                 case {WesPalette.Rushmore1, WesPalette.Rushmore}
                     hex = {'#E1BD6D','#EABE94','#0B775E','#35274A','#F2300F'};
-
                 case WesPalette.Royal1
                     hex = {'#899DA4','#C93312','#FAEFD1','#DC863B'};
-
                 case WesPalette.Royal2
                     hex = {'#9A8822','#F5CDB4','#F8AFA8','#FDDDA0','#74A089'};
-
                 case WesPalette.Zissou1
                     hex = {'#3B9AB2','#78B7C5','#EBCC2A','#E1AF00','#F21A00'};
-
                 case WesPalette.Zissou1Continuous
                     hex = {'#3A9AB2','#6FB2C1','#91BAB6','#A5C2A3','#BDC881', ...
                            '#DCCB4E','#E3B710','#E79805','#EC7A05','#EF5703','#F11B00'};
-
                 case WesPalette.Darjeeling1
                     hex = {'#FF0000','#00A08A','#F2AD00','#F98400','#5BBCD6'};
-
                 case WesPalette.Darjeeling2
                     hex = {'#ECCBAE','#046C9A','#D69C4E','#ABDDDE','#000000'};
-
                 case WesPalette.Chevalier1
                     hex = {'#446455','#FDD262','#D3DDDC','#C7B19C'};
-
                 case WesPalette.FantasticFox1
                     hex = {'#DD8D29','#E2D200','#46ACC8','#E58601','#B40F20'};
-
                 case WesPalette.Moonrise1
                     hex = {'#F3DF6C','#CEAB07','#D5D5D3','#24281A'};
-
                 case WesPalette.Moonrise2
                     hex = {'#798E87','#C27D38','#CCC591','#29211F'};
-
                 case WesPalette.Moonrise3
                     hex = {'#85D4E3','#F4B5BD','#9C964A','#CDC08C','#FAD77B'};
-
                 case WesPalette.Cavalcanti1
                     hex = {'#D8B70A','#02401B','#A2A475','#81A88D','#972D15'};
-
                 case WesPalette.GrandBudapest1
                     hex = {'#F1BB7B','#FD6467','#5B1A18','#D67236'};
-
                 case WesPalette.GrandBudapest2
                     hex = {'#E6A0C4','#C6CDF7','#D8A499','#7294D4'};
-
                 case WesPalette.IsleofDogs1
                     hex = {'#9986A5','#79402E','#CCBA72','#0F0D0E','#D9D0D3','#8D8680'};
-
                 case WesPalette.IsleofDogs2
                     hex = {'#EAD3BF','#AA9486','#B6854D','#39312F','#1C1718'};
-
                 case WesPalette.FrenchDispatch
                     hex = {'#90D4CC','#BD3027','#B0AFA2','#7FC0C6','#9D9C85'};
-
                 case WesPalette.AsteroidCity1
                     hex = {'#0A9F9D','#CEB175','#E54E21','#6C8645','#C18748'};
-
                 case WesPalette.AsteroidCity2
                     hex = {'#C52E19','#AC9765','#54D8B1','#B67C3B','#175149','#AF4E24'};
-
                 case WesPalette.AsteroidCity3
                     hex = {'#FBA72A','#D3D4D8','#CB7A5C','#5785C1'};
             end
         end
     end
+
+    methods (Static, Access = private)
+        function rgb = hex2rgb(hex)
+            % hex: cellstr like {'#RRGGBB', ...}
+            if isstring(hex), hex = cellstr(hex); end
+            n = numel(hex);
+            rgb = zeros(n,3);
+            for i = 1:n
+                h = char(hex{i});
+                if startsWith(h,'#'), h = h(2:end); end
+                if numel(h) ~= 6
+                    error('Invalid HEX color: %s', hex{i});
+                end
+                rgb(i,1) = hex2dec(h(1:2))/255;
+                rgb(i,2) = hex2dec(h(3:4))/255;
+                rgb(i,3) = hex2dec(h(5:6))/255;
+            end
+        end
+
+        function cmap = sample_discrete(base_rgb, n)
+            % No interpolation; allow any n by sampling/repeating.
+            k = size(base_rgb,1);
+
+            if n <= k
+                idx = round(linspace(1, k, n));      % spread across palette
+                idx = max(1, min(k, idx));
+                cmap = base_rgb(idx,:);
+            else
+                reps = floor(n / k);
+                rmd  = mod(n, k);
+                cmap = [repmat(base_rgb, reps, 1); base_rgb(1:rmd,:)];
+            end
+        end
+
+        function cmap = interpolate_continuous_lab(base_rgb, n)
+            % Smooth interpolation in Lab using colorspace().
+            % Requires colorspace.m by Pascal Getreuer on MATLAB path.
+
+            k = size(base_rgb,1);
+            if k == 1
+                cmap = repmat(base_rgb, n, 1);
+                return;
+            end
+
+            % Convert to Lab, interpolate each channel, convert back
+            lab = colorspace('Lab<-RGB', base_rgb);
+
+            t_base = linspace(0, 1, k);
+            t_new  = linspace(0, 1, n);
+
+            lab_new = zeros(n,3);
+            for c = 1:3
+                lab_new(:,c) = interp1(t_base, lab(:,c), t_new, 'linear');
+            end
+
+            rgb_new = colorspace('RGB<-Lab', lab_new);
+
+            % Clamp to displayable gamut
+            cmap = min(max(rgb_new, 0), 1);
+        end
+    end
 end
diff --git a/Libs/wesanderson_colors/minimal_example_wespalette.m b/Libs/wesanderson_colors/minimal_example_wespalette.m
index 9fcd99b..f308f27 100644
--- a/Libs/wesanderson_colors/minimal_example_wespalette.m
+++ b/Libs/wesanderson_colors/minimal_example_wespalette.m
@@ -5,8 +5,8 @@ y2 = 1e0 ./ (1 + exp(-0.4*(x-12)));  % NLPN
 y3 = 1e-6 * 10.^(0.45*x);             % RP on gamma
 y4 = 1e-2 * 10.^(0.18*(x-8));         % RP on beta2
 
-cmap = WesPalette.AsteroidCity1.rgb(4);
-cmap = linspecer(4);
+cmap = WesPalette.AsteroidCity1;
+% cmap = linspecer(4);
 figure1=figure(202998);clf;hold on
 lw = 0.8; ms = 4;
 plot(x,y1,'LineWidth',lw,'Color',cmap(1,:),'Marker','o','MarkerEdgeColor',cmap(1,:),'MarkerFaceColor',[1,1,1],'MarkerSize',ms); 
diff --git a/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/Copy_of_FIGURE_WAVELENGTH.m b/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/Copy_of_FIGURE_WAVELENGTH.m
index 21f881f..986fcdb 100644
--- a/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/Copy_of_FIGURE_WAVELENGTH.m
+++ b/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/Copy_of_FIGURE_WAVELENGTH.m
@@ -6,7 +6,7 @@ db = DBHandler("dataBase", "labor_highspeed", "type", database_type);
 
 pam_levels = [4, 6, 8];      % three tiles
 bitrate_set = 360e9;
-fiberL = 10;
+fiberL = 2;
 
 fields = [
     db.getTableFieldNames('power_state_info');
@@ -96,7 +96,7 @@ end
 %% ============================================================
 %  PLOT — 1×3 (PAM-4, PAM-6, PAM-8)
 % ============================================================
-fig = figure(9110); clf;
+fig = figure(9112); clf;
 tiledlayout(1,3,'TileSpacing','compact','Padding','compact');
 
 lw = 1.8;
@@ -223,19 +223,19 @@ ylabel('');
 
 end
 
-pos = 1e3.*[2.7770    1.2017    1.4000    0.3200];
-set(fig, 'Position', pos);
+% pos = 1e3.*[2.7770    1.2017    1.4000    0.3200];
+% set(fig, 'Position', pos);
 
 %% === EXPORT ===
-outfile = 'C:\Users\Silas\Documents\latex\JLT_400G copy\media\matlab2tikz\wavelength_analysis.tikz';
-matlab2tikz(outfile, ...
-    'width','\fwidth', ...
-    'height','\fheight', ...
-    'showInfo',false, ...
-    'extraAxisOptions',{ ...
-        'legend style={font=\footnotesize}', ...
-        'legend columns=1' ...
-    });
+% outfile = 'C:\Users\Silas\Documents\latex\JLT_400G copy\media\matlab2tikz\wavelength_analysis.tikz';
+% matlab2tikz(outfile, ...
+%     'width','\fwidth', ...
+%     'height','\fheight', ...
+%     'showInfo',false, ...
+%     'extraAxisOptions',{ ...
+%         'legend style={font=\footnotesize}', ...
+%         'legend columns=1' ...
+%     });
 
 
 
diff --git a/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/FIGURE_introduction.m b/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/FIGURE_introduction.m
index f9af8de..147c7d9 100644
--- a/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/FIGURE_introduction.m
+++ b/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/FIGURE_introduction.m
@@ -7,113 +7,133 @@ data = readtable(tablename,"Delimiter",';','DecimalSeparator',',');
 %   PLOT
 % ============================================================
 %% 1. DATA EXTRACTION & SETUP
-% Extract columns from the table
+%% 1. DATA EXTRACTION & SETUP
 raw_M       = data.M;
 raw_baud    = data.BaudRate;
 raw_net     = data.NetRate;
-raw_names   = string(data.ZoteroCode); 
-raw_band    = string(data.Band); % Spalte 'Band' als String extrahieren
+raw_codes   = string(data.ZoteroCode); 
+raw_names   = string(data.Name);
+raw_band    = string(data.Band);
 
-% Filter: Keep only PAM 2, 4, 6, 8 and rows where Rates are not NaN
+% Filter Valid Data
 target_M = [2, 4, 6, 8];
 validIdx = ismember(raw_M, target_M) & ~isnan(raw_baud) & ~isnan(raw_net);
 
-% Apply filter to create the working variables
 Mvals   = raw_M(validIdx);
 baud    = raw_baud(validIdx);
 netrate = raw_net(validIdx);
+codes   = raw_codes(validIdx);
 names   = raw_names(validIdx);
-bands   = raw_band(validIdx);    % Gefilterte Bands
+bands   = raw_band(validIdx);
 
-% Setup Lists for Formatting
 pam_list = target_M;
+colors = flip(cbrewer2('SET1',4));
 
-% Visual Settings
-colors = lines(4); % 4 Farben für PAM-2,4,6,8
-
-% 2. PLOT
+%% 2. PLOT (For Visual Check only)
 figure; hold on;
-ms = 40;      % Marker size (etwas größer für bessere Sichtbarkeit)
-lw = 1.5;     % Line width für Fit
+ms = 20;     
+lw = 0.5;   
 
-for k = 1:length(pam_list)   % Loop durch PAM-Formate
+for k = 1:length(pam_list) 
     M = pam_list(k);
-    
-    % Logical mask for current PAM
     idxPam = (Mvals == M);
     
-    % Extract for this PAM
     x = baud(idxPam);
     y = netrate(idxPam);
+    b = bands(idxPam);
     n = names(idxPam);
-    b = bands(idxPam); % Bands für diese PAM-Gruppe
-    
-    % Get color for this PAM format
     col = colors(k,:);
     
-    % ----- LEGEND FLAG -----
     firstLegend = true;
     
-    % ---- PLOT ALL POINTS (Marker basierend auf Band) ----
     for i = 1:sum(idxPam)
-        
-        % === MARKER LOGIC ===
-        currentBand = b(i);
-        if strcmpi(currentBand, 'O')
-            marker = 'o'; % Kreis für O-Band
-        elseif strcmpi(currentBand, 'C')
-            marker = 'd'; % Diamond für C-Band
+
+        % Marker Logic
+        ms = 20;   
+        if strcmpi(b(i), 'O')
+            marker = 'o'; 
+        elseif strcmpi(b(i), 'C')
+            marker = 'd'; 
         else
-            marker = 's'; % Fallback (Quadrat), falls andere Bands auftauchen
+            marker = 's'; 
+        end
+        if strcmpi(n(i), 'THIS WORK')
+            marker = 'pentagram'; 
+            ms = 100;
         end
         
-        % Plotting
+        % Plot Scatter
         if firstLegend
-            h = scatter(x(i), y(i), ms, ...
-                'Marker', marker, ...
-                'MarkerEdgeColor', col, ...
-                'MarkerFaceColor', col, ...
+            scatter(x(i), y(i), ms, 'Marker', marker, ...
+                'MarkerEdgeColor', col, 'MarkerFaceColor', col, ...
                 'DisplayName', sprintf('PAM-%d', M));
             firstLegend = false;
         else
-            h = scatter(x(i), y(i), ms, ...
-                'Marker', marker, ...
-                'MarkerEdgeColor', col, ...
-                'MarkerFaceColor', col, ...
+            scatter(x(i), y(i), ms, 'Marker', marker, ...
+                'MarkerEdgeColor', col, 'MarkerFaceColor', col, ...
                 'HandleVisibility','off');
         end
-        
-        % ====== CUSTOM DATATIP CONTENT ======
-        dt = h.DataTipTemplate;
-        dt.DataTipRows(1).Label = 'Baud rate';
-        dt.DataTipRows(2).Label = 'Net rate';
-        dt.DataTipRows(end+1) = dataTipTextRow('Band', b(i)); % Zeigt Band an
-        dt.DataTipRows(end+1) = dataTipTextRow('Paper', n(i));
     end
     
-    % ---- Fit (PAM-specific) ----
+    % Fit lines
     if length(x) >= 3
         [p, S, mu] = polyfit(x, y, 2); 
         xfit = linspace(min(x), max(x), 200);
         yfit = polyval(p, xfit, S, mu);
-        
-        plot(xfit, yfit, ':', ...
-            'LineWidth', lw, ...
-            'Color', col, ...
-            'HandleVisibility', 'off'); 
+        plot(xfit, yfit, '-', 'LineWidth', lw, 'Color', col, 'HandleVisibility', 'off'); 
     end
 end
 
 grid on; box on;
 xlabel('Baud rate [GBd]');
 ylabel('Net rate [Gb/s]');
-title('Transmission Records (Circle=O, Diamond=C)');
-legend('Location','northwest');
-set(gca,'FontSize',11);
+% title('Check Command Window for TikZ Code');
+% legend('Location','northwest');
 
+%% 3. GENERATE TIKZ ANNOTATION CODE
+% This prints the manual \draw commands to the console
+
+%% GENERATE TIKZ ANNOTATION CODE
+% This prints the manual \draw commands to the console
+
+%% GENERATE TIKZ ANNOTATION CODE (Colored Borders + Tiny Font)
+%% GENERATE TIKZ ANNOTATION CODE (No Arrow, Close Text)
+fprintf('\n\n%% ===========================================================\n');
+fprintf('%% COPY THE FOLLOWING LINES INTO YOUR .TEX FILE \n');
+fprintf('%% (Paste them just before \\end{axis})\n');
+fprintf('%% ===========================================================\n\n');
+
+for i = 1:length(baud)
+    bx = baud(i);
+    by = netrate(i);
+    key = codes(i);
+    M_val = Mvals(i);
+    
+    % --- PLACEMENT LOGIC ---
+    if M_val == 8
+        % PAM-8: Place Top-Left
+        % 'south east' anchor means the text's bottom-right corner touches the coordinate
+        % shift moves it slightly up and left to clear the marker
+        anchorStr = 'south east';
+        shiftStr  = 'shift={(-3pt, 3pt)}'; 
+    else
+        % Others: Place Bottom-Right
+        % 'north west' anchor means the text's top-left corner touches the coordinate
+        % shift moves it slightly down and right
+        anchorStr = 'north west';
+        shiftStr  = 'shift={(3pt, -3pt)}';
+    end
+    
+    % --- PRINT COMMAND ---
+    % Uses \node directly at the coordinate (axis cs:...)
+    fprintf('\\node[anchor=%s, %s, font=\\tiny, fill=white, inner sep=1pt] at (axis cs:%.2f, %.2f) {\\cite{%s}};\n', ...
+        anchorStr, shiftStr, bx, by, key);
+end
+fprintf('\n')
 
 %% === EXPORT ===
-outfile = 'C:\Users\Silas\Documents\latex\JLT_400G copy\media\matlab2tikz\highspeedresults.tikz';
+outfile = 'C:\Users\Silas\Documents\latex\JLT_400G_submission\media\matlab2tikz\highspeedresults_test.tikz';
+
 matlab2tikz(outfile, ...
     'width','\fwidth', ...
     'height','\fheight', ...
diff --git a/projects/ML_based_MLSE/analyze_filter_length.m b/projects/ML_based_MLSE/analyze_filter_length.m
index bcf2a91..1b31d6b 100644
--- a/projects/ML_based_MLSE/analyze_filter_length.m
+++ b/projects/ML_based_MLSE/analyze_filter_length.m
@@ -25,10 +25,10 @@ h = abs([0.3 0.9 0.3]); h = h/norm(h);
 
 % h = [1	-1.67085330039878	1.17918163282514	-0.805210559745616	0.571564213123367	-0.296337147529674	0.00649773445209780	0.0854177610195952	-0.0576009020965258	0.0520994427061551	-0.0624586034913656	0.0553280962699552	-0.00705582559925755	-0.0336399056707792	0.0706903719452810	-0.0334124287931977	0.0131699455037966	0.0587431373842994	-0.0515902976066452	0.00647904355473619	0.0137506750904990	-0.0547974515885928	0.00994735499340592	-0.0135513582534086	-0.00463322575007739	0.0277311946101940];
 % h = h/norm(h);
-symbols_filt = Symbols.filter(1,h);
+symbols_filt = Symbols.filter(h,1);
 symbols_noi = symbols_filt;
 symbols_noi.signal = awgn(symbols_filt.signal,SNR_dB,'measured');
-symbols_noi.spectrum;
+symbols_noi.spectrum();
 
 % --- Generate all parameter pairs
 [O,D] = ndgrid(order_range, delta_range);
@@ -42,7 +42,7 @@ ce_vec  = nan(size(pairs,1),1);
 ce_training = nan(size(pairs,1),training_len);
 
 % --- Parallel loop over parameter pairs
-parfor k = 1:size(pairs,1)
+for k = 1:size(pairs,1)
     order_k = pairs(k,1);
     delta_k = pairs(k,2);
 
@@ -55,7 +55,7 @@ parfor k = 1:size(pairs,1)
     try
         ml = ML_MLSE("epochs_tr",training_len,"epochs_dd",1,"len_tr",2^15, ...
             "mu_dd",0.1,"mu_tr",0.1,"order",order_k,"sps",1, ...
-            "traceback_depth",128,"L",2,"delta",delta_k,"adaptive_mu",0);
+            "traceback_depth",128,"L",3,"delta",delta_k,"adaptive_mu",0);
 
         [y_ml,y_ref] = ml.process(symbols_noi,Symbols);
         ref_bits = PAMmapper(M,0).demap(y_ref);