Merge branch 'main' of https://cau-git.rz.uni-kiel.de/nt/mitarbeiter/silas/imdd_simulation

# Conflicts:
#	Classes/00_signals/Signal.m
#	Functions/EQ_structures/duobinary_signaling.m
#	Functions/EQ_structures/vnle_postfilter_mlse.m
#	Functions/EQ_visuals/showLevelScatter.m

Classes/00_signals/Opticalsignal.m

@@ -5,6 +5,7 @@ classdef Opticalsignal < Signal
     properties
         nase
         lambda
+        polrot
 
     end
 
@@ -18,6 +19,7 @@ classdef Opticalsignal < Signal
                 options.logbook
                 options.lambda
                 options.nase
+                options.polrot
             end
 
             obj = obj@Signal(signal);

Classes/00_signals/Signal.m

@@ -108,6 +108,7 @@ classdef Signal
                 options.logbook
                 options.nase
                 options.lambda
+                options.polrot
             end
 
             fn = fieldnames(options);
@@ -122,7 +123,7 @@ classdef Signal
 
                 %convert to optical
 
-                o_sig = Opticalsignal(obj.signal,"fs",obj.fs,"lambda",options.lambda,"logbook",obj.logbook,"nase",options.nase);
+                o_sig = Opticalsignal(obj.signal,"fs",obj.fs,"lambda",options.lambda,"logbook",obj.logbook,"nase",options.nase,"polrot",options.polrot);
 
 
             elseif isa(obj,'Informationsignal')
@@ -141,7 +142,7 @@ classdef Signal
 
             arguments
                 obj
-                options.fignum
+                options.fignum = randi(1000)
                 options.displayname = [];
                 options.timeframe = 0;
                 options.clear = 0;
@@ -187,7 +188,7 @@ classdef Signal
                 ylabel('Amplitude');
             end
 
-           
+
             % Add legend if not already present
             if isempty(get(gca, 'Legend'))
                 legend;
@@ -269,6 +270,7 @@ classdef Signal
             SignalPower = [obj.power];
             Nase        = [0];
             SignalCopy  = obj.signal;
+            SignalCopy  = [];
             ModifierName = class(CallingModifier);
 
 
@@ -342,83 +344,151 @@ classdef Signal
 
             arguments
                 obj
-                options.fignum
+                options.fignum = 2025
                 options.displayname = "";
                 options.color = [];
+                options.linestyle = '-';
                 options.normalizeToNyquist = 0;
+                options.addDCoffset = 0;
+                options.normalizeToDC = 0;
                 options.normalizeTo0dB = 0;
                 options.max_num_lines = []; % Leave empty or omit to disable line rotation
                 options.fft_length = [];
+                % --- NEW options ---
+                options.useWavelengthAxis (1,1) logical = false  % plot x-axis in wavelength
+                options.lambda0_nm (1,1) double = 1310           % center wavelength [nm]
             end
 
             if isempty(options.fft_length)
                 options.fft_length = 2^(nextpow2(length(obj.signal))-9);
             end
+
             
+
             if options.normalizeToNyquist == 0
-                [p_lin,w] = pwelch(obj.signal,hanning(options.fft_length),options.fft_length/2,options.fft_length,obj.fs,"centered","power","mean");
-                w = w.*1e-9;
+                [p_lin,f_Hz] = pwelch(obj.signal, hanning(options.fft_length), ...
+                    options.fft_length/2, options.fft_length, ...
+                    obj.fs, "centered", "power", "mean");
+                f_GHz = f_Hz*1e-9; % keep frequency vector for frequency axis
             else
-                [p_lin,w] = pwelch(obj.signal,hanning(options.fft_length),options.fft_length/2,options.fft_length,"centered","power","mean");
+                [p_lin,f_rad] = pwelch(obj.signal, hanning(options.fft_length), ...
+                    options.fft_length/2, options.fft_length, ...
+                    "centered", "power", "mean");
+                % In normalized mode, pwelch returns rad/sample centered on 0.
+                % We'll keep f_rad for the x-axis in that mode.
             end
 
+            % p_lin = movmean(p_lin,4);
+
             if options.normalizeTo0dB
-                p_lin = p_lin./ max(p_lin);
+                p_lin = p_lin ./ max(p_lin);
                 p_dbm = 10*log10(p_lin); % normalized to 0 dB
-                ylab = "normalized to 0 dB";
+                ylab = "Normalized PSD";
             else
                 p_dbm = 10*log10(p_lin);
                 ylab = "Power (dB/Hz)";
             end
 
+            % --- If requested, build wavelength axis from frequency offset ---
+            if options.useWavelengthAxis && options.normalizeToNyquist == 0
+                c = physconst('LightSpeed');         % [m/s]
+                lambda0_m = options.lambda0_nm*1e-9; % center wavelength [m]
+                f_c = c / lambda0_m;                 % carrier frequency [Hz]
+
+                % exact mapping
+                f_abs    = f_c + f_Hz;               % absolute frequency [Hz]
+                lambda_m = c ./ f_abs;               % wavelength [m]
+                lambda_nm = lambda_m * 1e9;          % wavelength [nm]
+
+                % assign axis
+                x_vec   = lambda_nm(:);
+                x_label = "Wavelength [nm]";
+
+                % Sort to ensure axis is ascending
+                [x_vec, sortIdx] = sort(x_vec, 'ascend');
+                p_dbm = p_dbm(sortIdx, :);
+            else
+                % Frequency or normalized axes
+                if options.normalizeToNyquist == 0
+                    x_vec   = f_GHz;
+                    x_label = "Frequency in GHz";
+                else
+                    x_vec   = f_rad;     % normalized frequency in rad/sample
+                    x_label = "Normalized Frequency";
+                end
+            end
+
             figure(options.fignum);
             ax = gca;
             hold on
 
-            if isempty(options.color)
-                hLine = plot(w,p_dbm,'DisplayName',options.displayname,'LineWidth',1);
-            else
-                hLine = plot(w,p_dbm,'DisplayName',options.displayname,'LineWidth',1,'Color',options.color);
+            p_dbm = p_dbm+options.addDCoffset;
+
+            if options.normalizeToDC
+                [~,min_idx]=min(abs(f_GHz));
+                pow_at_dc = p_dbm(min_idx);
+                p_dbm = p_dbm-pow_at_dc;
+            end
+            % p_dbm = movmean(p_dbm,10);
+
+            for s = 1:min(size(p_dbm))
+                if isempty(options.color)
+                    plot(x_vec, p_dbm(:,s), 'DisplayName', options.displayname, 'LineWidth', 1);
+                else
+                    plot(x_vec, p_dbm(:,s), 'DisplayName', options.displayname, 'LineWidth', 1, 'Color', options.color,'LineStyle',options.linestyle);
+                end
             end
 
-            % If user wants to limit the number of lines, check and remove old lines
+            % Limit number of lines if requested
             if ~isempty(options.max_num_lines) && options.max_num_lines > 0
                 allLines = findall(ax, 'Type', 'Line');
                 if length(allLines) > options.max_num_lines
-                    % Sort lines by creation order. Usually, the oldest lines appear first in allLines.
-                    % If needed, you can sort by UserData or other criteria.
                     numToRemove = length(allLines) - options.max_num_lines;
                     delete(allLines(1:numToRemove));
                 end
             end
 
-            if options.normalizeToNyquist == 0
-                xlabel("Frequency in GHz");
-                edgetick = 2^(nextpow2(obj.fs*1e-9));
-                xticks(-edgetick:16:edgetick);
-                xlim([100*round(min(w)/100,1)-10, 100*round(max(w)/100,1)+10])
-                xlim([-128 128]);%256GSa/s
+            % Axis labels and limits
+            xlabel(x_label);
+
+            if options.useWavelengthAxis && options.normalizeToNyquist == 0
+                xlim([min(x_vec) max(x_vec)]);
             else
-                xlabel("Normalized Frequency");
-                xlim([-pi, pi]);
+                if options.normalizeToNyquist == 0
+                    % Keep your existing freq handling (you can fine-tune as needed)
+                    % xlim([-128 128]); % example for 256 GSa/s if desired
+                    xlim([min(x_vec) max(x_vec)]);
+                else
+                    xlim([-pi, pi]);
+                end
             end
 
             ylabel(ylab);
 
-            try
-            ylim([max(min(floor(min(p_dbm))-3, ax.YLim(1)),-40), min(max(ceil(max(p_dbm))+3, ax.YLim(2)),10)]);
-            catch
-            ylim([floor(min(p_dbm))-3, ceil(max(p_dbm))+3]);
+            % --- Y-Axis scaling (auto with margin) ---
+            y_min = min(p_dbm(:));
+            y_max = max(p_dbm(:));
+
+            % Add 5% dynamic range margin on both sides
+            y_range = y_max - y_min;
+            if y_range == 0
+                y_range = 10; % fallback if flat
             end
-            ylim([floor(min(p_dbm))-3, ceil(max(p_dbm))+3]);
-            yticks(-200:10:10);
-            grid on; grid minor;
-            legend
-            % legend('Interpreter','none');
+
+            y_margin = 0.05 * y_range;
+            ylim([y_min - y_margin, y_max + y_margin]);
+
+            % Set ticks automatically, avoid overpopulation
+            try
+            yticks(round(linspace(y_min, y_max, min(10, max(4, ceil(y_range/10))))));
+            end
+            grid on;
+
 
         end
 
 
+
         function move_it_spectrum(obj,options)
             arguments
                 obj
@@ -546,7 +616,7 @@ classdef Signal
                 options.unit power_notation = power_notation.dBm
             end
 
-            pow = mean(abs(obj.signal).^2);
+            pow = sum(mean(abs(obj.signal).^2));
 
             switch options.unit
                 case power_notation.dBm
@@ -661,9 +731,9 @@ classdef Signal
                 options.fs_ref = 0;
                 options.debug_plots = 0;
             end
-     
-            
-            
+
+
+
             S = {};
             inverted = -1;
             sequenceFound = 0;
@@ -693,10 +763,10 @@ classdef Signal
                 return
             end
 
-            if mean(w) > 10 
-                warning('No sync possible, check code of findpeaks! Look at correlation')
+            pkpos = sort(pkpos);
+
+            if mean(w) > 10 || mean(p) > 10
                 return
-                
             else
                 sequenceFound = 1;
             end
@@ -707,33 +777,33 @@ classdef Signal
                 findpeaks(abs(co./max(co)),'MinPeakDistance',length(b)/2,'MinPeakHeight',0.2,'NPeaks',maxpeaknum,'SortStr','descend')
             end
 
-            
+
 
             shifts = lags(pkpos);
             sequenceStarts = shifts;
-            % shifts = shifts(shifts>=0);
+            shifts = shifts(shifts>=0);
 
             if numel(shifts) > 0
-    
+
                 %Cut occurences of ref signal from signal (only positive shifts)
                 if all(sign(co(pkpos))==-1)
                     inverted = 1;
                 end
-               
+
                 for c = shifts
                     sig = obj.delay(-c,'mode','samples');
-                    sig.signal = sig.signal(1:length(b)) .* -inverted;
+                    sig.signal = sig.signal(1:length(b));% .* -inverted;
                     S{end+1,1} = sig;
                 end
-    
-                %return/keep the sinal with the highest correlation (only within positive shifts)
-                [~,idx]=max(pks);
-                obj.signal = S{idx}.signal;
-                %put signal with highest corr. to first index in S array
-                swap = S{1};
-                S{1} = S{idx};
-                S{idx} = swap;
-    
+
+                % %return/keep the sinal with the highest correlation (only within positive shifts)
+                % [~,idx]=max(pks);
+                % obj.signal = S{idx}.signal;
+                % %put signal with highest corr. to first index in S array
+                % swap = S{1};
+                % S{1} = S{idx};
+                % S{idx} = swap;
+
                 for c = 1:numel(shifts)
                     S{c}.logbook = [];
                 end
@@ -762,10 +832,6 @@ classdef Signal
         end
 
 
-
-
-
-
         %%
         function obj = filter(obj,a,b)
 
@@ -912,11 +978,11 @@ classdef Signal
             elseif mode == 1
                 % generate eye diagram using histogram
 
-                maxA = max(sig(100:end-100))*2;
-                minA = min(sig(100:end-100))*2;
+                maxA = max(sig(100:end-100))*1.3;
+                minA = min(sig(100:end-100))*1.3;
 
-                % maxA = 0.0015;
-                % minA = 0;
+                maxA = 0.12;
+                minA = -0.08;
 
                 difference= maxA-minA;
 
@@ -927,7 +993,7 @@ classdef Signal
                     hist_data(:,n)=flip(nn.'); %without flip, the eye is upside down :-(
                 end
 
-                ax = gca;
+
 
                 plot_data = 20*log10(hist_data);
                 plot_data(plot_data==-Inf) = 0;
@@ -944,27 +1010,29 @@ classdef Signal
                 if isa(obj,'Opticalsignal')
                     title(['Optical Eye ',options.displayname])
                     ylabel("Power in mW");
-                    y_tickstring = string(linspace(maxA.*1e3,minA.*1e3,16));
+                    y_tickstring = string(linspace(maxA.*1e3,minA.*1e3,6));
                     min_ = min(abs(obj.signal(100:end-100)).^2);
                     max_ = abs(max(obj.signal(100:end-100)).^2);
                 elseif isa(obj,'Electricalsignal')
                     title(['Electrical Eye ',options.displayname])
                     ylabel("Voltage in V");
-                    y_tickstring = string(linspace(maxA,minA,16));
+                    y_tickstring = string(linspace(maxA,minA,6));
                     min_ = min(obj.signal(100:end-100));
                     max_ = abs(max(obj.signal(100:end-100)));
                 else
                     title(['Digital Eye ',options.displayname])
                     ylabel("Digital Signal Amplitude");
-                    y_tickstring = string(linspace(maxA,minA,16));
+                    y_tickstring = string(linspace(maxA,minA,6));
                     min_ = min(obj.signal(100:end-100));
                     max_ = abs(max(obj.signal(100:end-100)));
                 end
                 xlabel('Time in ps')
 
+
+
                 % add information
 
-                if 1
+                if 0
 
                     pwr_dbm = round(obj.power,3);
                     pwr_lin = obj.power("unit",power_notation.W);
@@ -1073,18 +1141,20 @@ classdef Signal
 
                     end
 
-                    yticks(linspace(0,histpoints,16));
-                    y_tickstring = sprintfc('%.2f', y_tickstring);
-                    yticklabels(y_tickstring);
+                    grid off
 
-                    xticks(linspace(0,histpoints_horizontal,8))
-                    x_tickstring = sprintfc('%.2f', linspace(0, 2/fsym, 8) .* 1e12);
-                    xticklabels(x_tickstring);
                 end
 
 
+                yticks(linspace(0,histpoints,6));
+                y_tickstring = sprintfc('%.2f', y_tickstring);
+                yticklabels(y_tickstring);
 
+                xticks(linspace(0,histpoints_horizontal,6))
+                x_tickstring = sprintfc('%.2f', linspace(0, 2/fsym, 8) .* 1e12);
+                xticklabels(x_tickstring);
 
+                % 
             end
 
             % disp('h');

Classes/01_transmit/ChannelFreqResp.m

@@ -139,6 +139,8 @@ classdef ChannelFreqResp < handle
             fnew = linspace(0,fstarget/2,length(Target)/2+1);
             fnew = fnew(2:end-1);
 
+
+
             % Old frequency axis (should be much coarser)
             idx_old  = find((obj.faxis > 0) .* (obj.faxis < fstarget/2)); %positions of all Frequencies smaller than fs/2
             int_fold = obj.faxis(idx_old); %old frequencies from 0 to fs/2
@@ -149,6 +151,10 @@ classdef ChannelFreqResp < handle
             % interpolate the frequency response that had a coarse frequency resolution (e.g. 256 bins) to the current frequency resolution (e.g. 21843 bins)
             iH  = interp1(int_fold, real(H_inv(idx_old)) ,fnew, 'linear') + 1i*interp1(int_fold, imag(H_inv(idx_old)) ,fnew, 'linear');
 
+            if 0
+                figure(7);hold on;plot(fnew,20*log10(abs(iH)))
+            end
+
             % set all NaN values to the fist/ last non-NaN value
             nH    = find(~isnan(iH),1,'first');
             iH(1:nH)=iH(nH);
@@ -171,6 +177,8 @@ classdef ChannelFreqResp < handle
             % five frequencies -> should be the vaue at f=0=DC component?
             iH  = iH./mean(abs(iH)); %why 1:5??
 
+            dc = 20*log10(abs(mean(abs(iH(1:5)))));
+
             % set maximum amplification
             % set als values higher than hmax to hmax and keep the
             % phase information by multiplication with respective
@@ -206,7 +214,7 @@ classdef ChannelFreqResp < handle
 
         function plot(obj)
 
-            figure(55);
+            figure();
             clf;
 
             Havg = obj.H;
@@ -225,7 +233,7 @@ classdef ChannelFreqResp < handle
             xlim([0.2 .5*max(obj.faxis)*1e-9]);
             grid on;
 
-            figure(56);
+            figure();
             clf;
 
 
@@ -246,7 +254,7 @@ classdef ChannelFreqResp < handle
             xlim([0.2 .5*max(obj.faxis)*1e-9]); grid on;
 
             %%% plot for publication
-            figure(1234);
+            figure(101);
             hold all;
             box on;
             title('Magnitude Freq. Response');
@@ -269,6 +277,25 @@ classdef ChannelFreqResp < handle
             legend('Interpreter','latex')
             grid on;
 
+            % 
+            figure(100); hold on;
+            Havg = obj.H;
+            Havg = Havg./max(Havg);
+            Hall = obj.H_all;
+
+            for i = 1:size(Hall,1)
+                Hall(i,:) = Hall(i,:)./max(Hall(i,:));
+            end
+
+            %1)
+            col = cbrewer2('Paired',8);
+            hold all;box on;title('Magnitude Freq. Response');
+            plot(obj.faxis/1e9, 20*log10(abs(Hall)),'linewidth',0.1,'LineStyle','-','Color',col(1,:),'HandleVisibility','off') ;
+            xlim([0.2 .5*max(obj.faxis)*1e-9]);
+            plot(obj.faxis/1e9, 20*log10(abs(Havg)),'LineWidth',2,'Color',col(2,:));
+            grid on;
+
+
 
         end
 

Classes/01_transmit/M8199B.m

@@ -17,7 +17,7 @@ classdef M8199B < AWG
 
             fdac            = 256e9;
             
-            Lp_awg      = Filter('filtdegree',4,"f_cutoff",75e9,"fs",fdac*options.kover,"filterType",filtertypes.butterworth,"active",true);
+            Lp_awg      = Filter('filtdegree',3,"f_cutoff",75e9,"fs",fdac*options.kover,"filterType",filtertypes.gaussian,"active",true);
 
             obj = obj@AWG("fdac",fdac,"dac_min",dac_min,"dac_max",dac_max,"lpf_active",1,"H_lpf",Lp_awg,"kover",options.kover,...
                         "bit_resolution",5.5,"normalize2dac",1,"upsampling_method","samplehold");

Classes/01_transmit/PAMmapper.m

@@ -29,7 +29,7 @@ classdef PAMmapper
 
             obj.levels = obj.get_levels();
 
-            obj.scaling = rms(obj.get_levels());
+            obj.scaling = obj.get_scaling;%rms(obj.get_levels());
 
             obj.eth_style = options.eth_style;
 
@@ -45,7 +45,7 @@ classdef PAMmapper
                 signal_in  = signal_in.logbookentry(lbdesc,obj);
                 out = signal_in;
             else
-                out = signal_in;
+                out = obj.map_(signal_in);
             end
 
         end
@@ -283,6 +283,21 @@ classdef PAMmapper
             end
         end
 
+        function scaling = get_scaling(obj)
+            switch obj.M
+                case 2
+                    scaling = 1;
+                case 4
+                    scaling = sqrt(5);
+                case 6
+                    scaling = sqrt(10);
+                case 8
+                    scaling = sqrt(21);
+                case 16
+                    scaling = 1;
+            end
+        end
+
         function [data_out] = demap_(obj,data_in)
             
             data_in= data_in';
@@ -328,8 +343,9 @@ classdef PAMmapper
 
                     if ~obj.eth_style
 
-                        data_in = data_in/(sqrt(mean(abs(data_in).^2)));
+                        % data_in = data_in/(sqrt(mean(abs(data_in).^2)));
                         data_in = data_in*sqrt(10);
+                        %data_in = data_in*sqrt((5^2 + 3^2 + 1^2 + 5^2 + 3^2 + 1^2)/6);
 
                         if size(data_in,2) > 1
                             data_in = data_in.';
@@ -449,7 +465,7 @@ classdef PAMmapper
         function [Signal_out] = quantize(obj,Signal_in)
 
             constellation = obj.get_levels();
-            constellation = constellation ./ rms(constellation);
+            constellation = constellation ./ obj.scaling;
 
             issignalclass = 0;
             if isa(Signal_in,'Signal')
@@ -480,7 +496,9 @@ classdef PAMmapper
         end
 
         function bitmap = showBitMapping(obj)
-            bitmap = obj.demap([obj.levels ./ obj.scaling]');
+
+            bitmap = obj.demap((obj.levels ./ obj.scaling)');
+
         end
 
     end

Classes/01_transmit/PAMsource.m

@@ -136,6 +136,13 @@ classdef PAMsource
 
             %%%%%% Duobinary %%%%%%%%%%%
 
+            if obj.db_precode
+                obj.duobinary_mode = db_mode.db_precoded;
+            end
+            if obj.db_encode
+                obj.duobinary_mode = db_mode.db_encoded;
+            end
+
             switch obj.duobinary_mode
 
                 case db_mode.no_db
@@ -164,7 +171,7 @@ classdef PAMsource
             %%%%% Pulse-forming %%%%%%
             if obj.applypulseform
                 %%% MY CODE
-                if 0
+                if 1
                     digi_sig = obj.pulseformer.process(symbols);
 
                 else

Classes/02_etc/Amplifier.m

@@ -108,7 +108,7 @@ classdef Amplifier
             if obj.gain_mode == gain_mode.output_power
 
                 %get linear gain for output power mode
-                pow_in = mean(abs(xin.^2)) ; % lin input power
+                pow_in = sum(mean(abs(xin.^2)),2) ; % lin input power
                 pow_out = 10^(obj.amplification_db/10 - 3) ; % dBm to lin
                 a_lin = sqrt(pow_out/pow_in) ;
 

Classes/02_optical/EML.m

@@ -70,7 +70,7 @@ classdef EML
             [signalclass_in.signal,obj] = obj.process_(signalclass_in.signal); 
 
              % cast the inform. signal to electrical signal 
-            signalclass_in = Opticalsignal(signalclass_in,"fs",obj.fsimu,"logbook",signalclass_in.logbook,"lambda",obj.lambda*1e-9,"nase",0);
+            signalclass_in = Opticalsignal(signalclass_in,"fs",obj.fsimu,"logbook",signalclass_in.logbook,"lambda",obj.lambda*1e-9,"nase",0,"polrot",0);
 
              % append to logbook
              lbdesc = [num2str(obj.lambda),' nm Laser with ',num2str(obj.power),' dBm P_out. Linew.=',num2str(obj.linewidth*1e-6),' MHz. Modulation mode: ',char(obj.mode) ];

Classes/02_optical/Fiber.m

@@ -1,189 +1,171 @@
 classdef Fiber
-    %FIBER Summary of this class goes here
-    %   Detailed explanation goes here
+    % Fiber: Simulate optical fiber signal propagation using
+    % the split-step Fourier method (SSFM).
 
     properties
+        % Simulation sampling frequency [Hz]
         fsimu
+        % Fiber length [m]
         fiber_length
+        % Attenuation coefficient [dB/km]
         alpha
+        % Dispersion parameter D [s/m^2]
         D
+        % Dispersion slope [s/m^3]
         Dslope
+        % Reference wavelength [m]
         lambda0
+        % Nonlinear coefficient [1/(W·m)]
         gamma
+        % Maximum allowed nonlinear phase change per step [rad]
         dphimax
 
+        % Derived parameters:
+        % Second-order dispersion coefficient β2 [s^2/m]
         b2
+        % Third-order dispersion coefficient β3 [s^3/m]
         b3
+        % Linear attenuation constant [1/m]
         alpha_lin
+        % Frequency-domain linear operator per unit length
         linstep
     end
 
     methods
         function obj = Fiber(options)
-            %FIBER Construct an instance of this class
-            %   Detailed explanation goes here
+            % Constructor: initialize fiber physical and simulation parameters.
             arguments
-                options.fsimu
-                options.fiber_length = 0
-                options.alpha = 0.2
-                options.D = 17
-                options.Dslope = 0.06
-                options.lambda0 = 1550
-                options.gamma = 0
-                options.dphimax = 5e-3
+                options.fsimu            % Sampling frequency [Hz]
+                options.fiber_length = 0 % Fiber length [km]
+                options.alpha = 0.2      % Attenuation [dB/km]
+                options.D = 17           % Dispersion D parameter [ps/(nm·km)]
+                options.Dslope = 0.06    % Dispersion slope [ps/(nm^2·km)]
+                options.lambda0 = 1550   % Reference wavelength [nm]
+                options.gamma = 0        % Nonlinear coefficient [1/(W·km)]
+                options.dphimax = 5e-3   % Max phase step [rad]
             end
 
-            obj.fsimu           = options.fsimu;
-            obj.fiber_length    = options.fiber_length*1000; %km
-            obj.alpha           = options.alpha;
-            obj.D               = options.D*1e-6;
-            obj.Dslope          = options.Dslope*1e3;
-            obj.lambda0         = options.lambda0*1e-9;
-            obj.gamma           = options.gamma;
-            obj.dphimax         = options.dphimax;
-
-
-
+            % Assign inputs to object properties, converting units to SI.
+            obj.fsimu        = options.fsimu;
+            obj.fiber_length = options.fiber_length * 1e3;   % km → m
+            obj.alpha        = options.alpha;
+            obj.D            = options.D * 1e-6;            % ps/(nm·km) → s/(m·m)
+            obj.Dslope       = options.Dslope * 1e3;        % ps/(nm^2·km) → s/m^2
+            obj.lambda0      = options.lambda0 * 1e-9;      % nm → m
+            obj.gamma        = options.gamma;               % 1/(W·km) (assumed → 1/(W·m) internally)
+            obj.dphimax      = options.dphimax;
         end
 
-        function signalclass_out = process(obj,signalclass_in)
+        function signalclass_out = process(obj, signalclass_in)
+            % process: Apply fiber propagation to input signal class.
+            %   Calls the internal SSFM routine and logs the operation.
 
-            % actual processing of the signal (steps 1. - 3.)
-            signalclass_in = obj.process_(signalclass_in);
+            % Run internal split-step propagation
+            signalclass = obj.process_(signalclass_in);
 
-            % append to logbook
-            lbdesc = 'Fiber ';
-            signalclass_in  = signalclass_in.logbookentry(lbdesc);
-
-            % write to output
-            signalclass_out = signalclass_in;
+            % Add entry to logbook for tracking
+            signalclass = signalclass.logbookentry('Fiber ');
 
+            % Return processed signal class
+            signalclass_out = signalclass;
         end
 
-        function opt_sig = process_(obj,opt_sig)
-            %METHOD1 Summary of this method goes here
-            %   Detailed explanation goes here
+        function opt_sig = process_(obj, opt_sig)
+            % process_: Internal routine for one-step fiber propagation.
+            %   Computes linear and (optionally) nonlinear effects.
 
-            signal = opt_sig.signal;
+            % Extract time-domain field and wavelength
+            signal        = opt_sig.signal;
             lambda_signal = opt_sig.lambda;
 
-            obj.D = obj.D + (lambda_signal-obj.lambda0)*obj.Dslope;
-            
-            obj.b2 = -obj.D*lambda_signal^2/(2*pi*Constant.LightSpeed);
-            obj.b3 = ((lambda_signal.^2/(2*pi*Constant.LightSpeed)).^2*obj.Dslope);
-            obj.alpha_lin = obj.alpha/10*log(10)/1000;
-            
+            % Update dispersion parameter D for current wavelength
+            obj.D = obj.D + (lambda_signal - obj.lambda0) * obj.Dslope;
+
+            % Compute dispersion coefficients (β2, β3)
+            obj.b2 = -obj.D * lambda_signal^2 / (2 * pi * Constant.LightSpeed);
+            obj.b3 = (lambda_signal^2 / (2 * pi * Constant.LightSpeed))^2 * obj.Dslope;
+
+            % Convert attenuation from dB/km to linear 1/m
+            obj.alpha_lin = obj.alpha / 10 * log(10) / 1e3;
+
+            % Build frequency axis for FFT operations
             N = length(signal);
-            faxis = linspace(-obj.fsimu/2,obj.fsimu/2,N+1);
-            faxis = ifftshift(faxis(:,1:end-1));
-            faxis = faxis';
-            
-            obj.linstep = -obj.alpha_lin/2 - 2*1j*pi^2*obj.b2*faxis.^2 - 4/3*1j*pi^3*obj.b3*faxis.^3;
+            faxis = linspace(-obj.fsimu/2, obj.fsimu/2, N+1);
+            faxis = faxis(1:end-1);          % drop redundant endpoint
+            faxis = ifftshift(faxis)';       % center zero frequency
 
-            if 0
-                H = exp((obj.linstep)*obj.fiber_length);
-
-                figure(222)
-                hold on
-                plot(faxis.*1e-9,abs(real(Y)),'LineStyle','-','DisplayName','Abs(real) part of complex TF');
-                xlabel("Frequency in GHz")
-                ylabel("$ R|(H(\omega, L))|$")
-            end
+            % Define linear operator per meter in frequency domain
+            obj.linstep = -obj.alpha_lin/2 ...                 % half-step loss
+                          - 1j*2*pi^2*obj.b2 .* faxis.^2 ... % second-order dispersion
+                          - 1j*(4/3)*pi^3*obj.b3 .* faxis.^3; % third-order dispersion
 
+            % Choose linear-only or nonlinear SSFM based on gamma
             if obj.gamma ~= 0
+                % Full adaptive SSFM with nonlinear Schrödinger solver
                 opt_out = obj.NLSE(signal);
             else
-                opt_out = ifft( fft(signal) .* exp(obj.linstep*obj.fiber_length) ); % only one linear step
+                % Single-step linear propagation in frequency domain
+                opt_out = ifft( fft(signal) .* exp(obj.linstep * obj.fiber_length) );
             end
 
-            %TODO: attenuate nase ... 
-
+            % Update output field in signal class
             opt_sig.signal = opt_out;
-
         end
 
-        function [yout] = NLSE(obj, xin)
-
-            maxPow = obj.gamma.*max(abs(xin).^2);
-
-            Leff = obj.dphimax / maxPow ;
-
-            dz = Leff;
+        function yout = NLSE(obj, xin)
+            % NLSE: Solve nonlinear Schrödinger equation by adaptive split-step
+            %   xin: input time-domain field
+            %   Returns yout: output time-domain field after propagation
 
+            % Initial estimate of effective length per step
+            maxPow = obj.gamma * max(abs(xin).^2);
+            Leff   = obj.dphimax / maxPow;
+            dz     = Leff;
             z_prop = 0;
 
-            yout = fft(xin);
-
-            yout = ((yout).*exp(obj.linstep*dz/2));
+            % Apply initial half-step linear operator
+            yout = fft(xin) .* exp(obj.linstep * dz/2);
 
+            % Loop until full fiber length is reached
             while true
-
+                % Inverse FFT to time domain for nonlinear phase shift
                 yout = ifft(yout);
 
-                Leff = dz;
-
+                % Nonlinear phase rotation per segment
                 power = abs(yout).^2;
-                Hnl = exp( -1j*obj.gamma*power*Leff);
-                yout = yout .* Hnl;
+                Hnl   = exp(-1j * obj.gamma * power * dz);
+                yout  = yout .* Hnl;
 
+                % Accumulate propagation distance
                 z_prop = z_prop + dz;
 
-                maxPow = obj.gamma*max(abs(yout).^2);
-                Leff = obj.dphimax/maxPow;
-
-                dz_new = Leff;
+                % Compute adaptive step for next interval
+                maxPow = obj.gamma * max(abs(yout).^2);
+                dz_new = obj.dphimax / maxPow;
 
+                % If remaining length shorter than new step, finish loop
                 if z_prop + dz_new > obj.fiber_length
                     dz_new = obj.fiber_length - z_prop;
-                    break
+                    break;
                 end
 
-                yout = fft(yout);
-
-                yout = ((yout).*exp(obj.linstep*(dz/2+dz_new/2)));
-
-                dz = dz_new;
-
+                % Half-step linear operator bridging segments
+                yout = fft(yout) .* exp(obj.linstep * ((dz/2) + (dz_new/2)));
+                dz   = dz_new;
             end
 
-            yout = fft(yout);
-
-            yout = ((yout).*exp(obj.linstep*(dz/2+dz_new/2)));
-
+            % Final propagation segment: combine half-steps and nonlinear
+            yout = fft(yout) .* exp(obj.linstep * ((dz/2) + (dz_new/2)));
             yout = ifft(yout);
 
-            Leff = dz;
-
+            % Last nonlinear phase shift
             power = abs(yout).^2;
+            Hnl   = exp(-1j * obj.gamma * power * dz_new);
+            yout  = yout .* Hnl;
 
-            Hnl = exp( -1j*obj.gamma*power*Leff);
-
-            yout = yout .* Hnl;
-
-            yout = fft(yout);
-
-            yout = ((yout).*exp(obj.linstep*(dz/2)));
-
+            % Final half-step linear operator to complete SSFM
+            yout = fft(yout) .* exp(obj.linstep * (dz_new/2));
             yout = ifft(yout);
-
-
         end
-
-
     end
 end
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-

Classes/04_DSP/Coding/Duobinary.m

@@ -53,6 +53,7 @@ classdef Duobinary
             %     bk(k+1) =mod(data(k)-bk(k),M);
             % end
 
+            % THIS WAS USED!
             for k = 2:numel(data)
                 bk(k) = mod(data(k)-bk(k-1),M);
             end

Classes/04_DSP/Equalizer/FFE.m

@@ -3,24 +3,34 @@ classdef FFE < handle
     % 1) Training mode (stable performance when you use NLMS)
     % 2) Decision directed mode
 
-    % Eq     = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",0);
+    %LMS: mu in order of 0.001 for acceptable convergence speed
+    %NLMS: mu in order of 0.01 for acceptable convergence speed
+    %RLS: mu is lambda -> 0.99 -> 1 (has a strong dependency on this! use a loop to find out best values)
+
+    %  FFE("epochs_tr",5,"epochs_dd",2,"len_tr",2^13,"mu_dd",mu_dd,"mu_tr",mu_tr,"order",25,"sps",2,"decide",0, "adaption",adaption_method(adaption),"dd_mode",use_dd_mode);
 
     properties
         sps % usually 2
         order
         e
+        e_tr
         error
+        debug_struct
 
         len_tr
         mu_tr
         epochs_tr
 
-        mu_dd
+        adaption_technique % nlms, lms, rls
+        dd_mode % 1 or 0 to set DD-mode on or off
+        mu_dd %weight update in dd mode
         epochs_dd
 
-        constellation
+        P % covariance matrix of rls
 
-        decide
+        constellation % symbol constellation
+
+        decide %wether to return the (hard) decisions or the result after FFE (soft)
     end
 
     methods
@@ -34,6 +44,8 @@ classdef FFE < handle
                 options.mu_tr = 0;
                 options.epochs_tr = 5;
 
+                options.adaption_technique adaption_method = adaption_method.lms;
+                options.dd_mode = 1; 
                 options.mu_dd = 1e-5;
                 options.epochs_dd = 5;
 
@@ -59,10 +71,15 @@ classdef FFE < handle
 
             obj.constellation = unique(D.signal);
 
+
+            delta   = 0.05;
+            obj.P = (1/delta) * eye(obj.order);
+
             % Training Mode
             training = 1;
             showviz = 0;
             obj.equalize(X.signal, D.signal,obj.mu_tr,obj.epochs_tr,obj.len_tr,training,showviz);
+            obj.e_tr = obj.e;
 
             % Decision Directed Mode
             n = X.length;
@@ -71,7 +88,7 @@ classdef FFE < handle
             [signal,decision]=obj.equalize(X.signal, D.signal,obj.mu_dd,obj.epochs_dd,n,training,showviz);
 
             % Output Signal
-            if obj.decide 
+            if obj.decide
                 X.signal = decision;
             else
                 X.signal = signal;
@@ -83,24 +100,25 @@ classdef FFE < handle
 
             Noi = X;
             Noi = X - D;
-            
+
 
         end
 
-        function [y,d_hat] = equalize(obj,x,d,mio,epochs,N,training,showviz)
+        function [y,d_hat] = equalize(obj,x,d,mu,epochs,N,training,showviz)
 
             arguments
-                obj 
-                x 
-                d 
-                mio 
-                epochs 
-                N 
-                training 
+                obj
+                x
+                d
+                mu
+                epochs
+                N
+                training
                 showviz
             end
 
             x = [zeros(floor(obj.order/2),1); x; zeros(obj.order,1)];
+            lambda  = mu;
 
             if training
                 mask = ones(obj.order,1);
@@ -111,8 +129,17 @@ classdef FFE < handle
             end
 
             mask = ones(obj.order,1);
+            maincursor_pos=ceil(length(obj.e)/2);
+            always_ideal_decision = 0;
+            save_debug = 0;
+            grad =0;
+            weight = 0;
+            update = 0;
+
+            if mu == 0 || (~obj.dd_mode && ~training)
+                epochs = 1;
+            end
 
-            
             for epoch = 1 : epochs
                 symbol = 0;
                 for sample = 1 : obj.sps : N
@@ -124,30 +151,75 @@ classdef FFE < handle
                     y(symbol,1) = (obj.e.*mask).' * U;    % Calculating output of LMS  __ * |
 
                     if training
-                        d_hat(symbol,1) = d(symbol);  
+                        d_hat(symbol,1) = d(symbol);
                     else
-                        [~,symbol_idx]  = min(abs(y(symbol) - obj.constellation)); % decision for closest constellation point
-                        d_hat(symbol,1) = obj.constellation(symbol_idx);
+                        if ~always_ideal_decision
+                            [~,symbol_idx]  = min(abs(y(symbol) - obj.constellation)); % decision for closest constellation point
+                            d_hat(symbol,1) = obj.constellation(symbol_idx);
+                        else
+                            d_hat(symbol,1) = d(symbol);
+                        end
                     end
 
-                    err(symbol) =  y(symbol) - d_hat(symbol);   % Instantaneous error
+                    % err(symbol) =  y(symbol) - d_hat(symbol);   % Instantaneous error
+                    err(symbol) =  d_hat(symbol) - y(symbol);   % Instantaneous error
 
                     true_err(symbol) =  y(symbol) - d(symbol);   % Instantaneous error
 
-                    if mio ~= 0
-                        obj.e = obj.e - (mio * err(symbol) * U) ;    % Weight update rule of LMS
-                    else
-                        normalizationfactor =  (U.' * U);
-                        obj.e = obj.e - err(symbol) * U / normalizationfactor;    % Weight update rule of NLMS
+                    if training || obj.dd_mode 
+                        switch obj.adaption_technique
+    
+                            case adaption_method.lms
+    
+                                % mu used as update weight (suggestion: 0.001)
+                                weight = mu;
+                                grad = err(symbol) * U;
+                                update = grad * weight;
+                                obj.e = obj.e + update;
+    
+                            case adaption_method.nlms
+    
+                                % mu used as update weight (suggestion: 0.01-0.05; bit higher during tr)
+                                normU = ((U.'*U)) + eps;
+                                weight = mu / normU;
+                                grad = err(symbol) * U;
+                                update = grad * weight;
+                                obj.e = obj.e + update;
+    
+                            case adaption_method.rls
+
+    
+                                % RLS‐Gain:
+                                denom = lambda + U.' * obj.P * U;
+                                k = (obj.P * U) / denom;
+    
+                                % Gewichtsupdate:
+                                update = k * err(symbol);
+                                obj.e = obj.e + update;
+    
+                                % P-Matrix‐Update:
+                                obj.P = (1/lambda) * (obj.P - k * (U.' * obj.P));
+   
+    
+                        end
                     end
 
-                    obj.error(epoch,symbol) = err(symbol) * err(symbol)';  % Instantaneous square error
+
+                    if save_debug
+                        obj.debug_struct.error(epoch,symbol) = err(symbol) * err(symbol)';  % Instantaneous square error
+                        if training
+                            obj.debug_struct.error_tr(epoch,symbol) = err(symbol) * err(symbol)';  % Instantaneous square error
+                            obj.debug_struct.update_tr(epoch,symbol) = update.'*update ./ rms(obj.e);
+                        end
+                        obj.debug_struct.main_cursor(epoch,symbol) = abs(obj.e(maincursor_pos));
+                        obj.debug_struct.mu_nlms(epoch,symbol) = weight;
+                        obj.debug_struct.update_gradient(epoch,symbol) = grad.'*grad;
+                        obj.debug_struct.update(epoch,symbol) = update.'*update ./ rms(obj.e);
+                    end
 
                 end
             end
 
         end
-
     end
-end
-
+end

Classes/04_DSP/Equalizer/FFE_DCremoval.m

@@ -57,6 +57,8 @@ classdef FFE_DCremoval < handle
             obj.e = zeros(obj.order,1);
             obj.error = 0;
 
+            obj.dc_buffer_len = floor(obj.dc_buffer_len);
+
         end
 
         function [X,Noi] = process(obj, X, D)
@@ -155,16 +157,15 @@ classdef FFE_DCremoval < handle
                 end
             end
 
-            if ~training
-                figure(1122)
-                hold on
-                scatter(1:numel(e_dc_save),e_dc_save,1,'.');
-                scatter(1:numel(y),y,1,'.');
-
-            end
+            % if ~training
+            %     figure(1122)
+            %     hold on
+            %     scatter(1:numel(e_dc_save),e_dc_save,1,'.');
+            %     % scatter(1:numel(y),y,1,'.');
+            % 
+            % end
             
         end
         
     end
-end
-
+end

Classes/04_DSP/Equalizer/FFE_DFE.m

@@ -60,7 +60,7 @@ classdef FFE_DFE < handle
 
         end
 
-        function [X] = process(obj, X, D)
+        function [X,Noi] =  process(obj, X, D)
 
             % actual processing of the signal (steps 1. - 3.)
             % 1 normalize RMS
@@ -88,6 +88,9 @@ classdef FFE_DFE < handle
             X.fs = D.fs; %change sampling frequency of outgoing signal from fdac e.g. 2 sps to symbol spaced = fsym
             lbdesc = [num2str(obj.ffe_order),' tap FFE'];
             X  = X.logbookentry(lbdesc); % append to logbook
+
+            Noi = X;
+            Noi = X - D;
             
 
         end
@@ -185,8 +188,6 @@ classdef FFE_DFE < handle
             obj.e = coeff(1:obj.ffe_order);
             obj.b = coeff(obj.ffe_order+1:end);
 
-
-
         end
 
     end

Classes/04_DSP/Equalizer/Postfilter.m

@@ -31,7 +31,7 @@ classdef Postfilter < handle
 
         end
 
-        function signalclass_out = process(obj,signalclass_in,noiseclass_in,options)
+        function [signalclass_out,noiseclass_out] = process(obj,signalclass_in,noiseclass_in,options)
 
             arguments
                 obj 
@@ -55,7 +55,8 @@ classdef Postfilter < handle
             else
                 
             end
-           
+
+            noiseclass_in = noiseclass_in.filter(obj.coefficients,1);
             signalclass_in = signalclass_in.filter(obj.coefficients,1);
 
             % append to logbook
@@ -64,7 +65,7 @@ classdef Postfilter < handle
 
             % write to output
             signalclass_out = signalclass_in;
-
+            noiseclass_out = noiseclass_in;
         end
 
         function showFilter(obj,options)

Classes/04_DSP/Sequence Detection/MLSE.m

@@ -6,6 +6,10 @@ classdef MLSE < handle
         DIR
         trellis_states
         duobinary_output
+        trellis_state_mode
+        trellis_exclusion
+        debug
+        scale_mode
     end
 
     methods (Access=public)
@@ -19,7 +23,10 @@ classdef MLSE < handle
                 options.DIR double = [1];
                 options.trellis_states double = [-3 -1 1 3];
                 options.duobinary_output logical = false;
-
+                options.trellis_state_mode = 2;
+                options.trellis_exclusion = 0;
+                options.scale_mode = 2;
+                options.debug = 0;
             end
 
             %
@@ -29,28 +36,64 @@ classdef MLSE < handle
                     obj.(fn{n}) = options.(fn{n});
                 end
             end
-
-            % do more stuff
-
         end
 
-        function [signalclass_hd,signalclass_sd] = process(obj,signalclass,ref_symbolclass)
+        function [signalclass_hd,LLR,GMI] = process(obj,signalclass,ref_symbolclass)
+
+            arguments
+                obj
+                signalclass
+                ref_symbolclass
+            end
 
             data_in =   signalclass.signal;
+            shape_in = size(data_in);
             data_ref = ref_symbolclass.signal;
 
-            [data_out_hd,data_out_sd] = obj.process_(data_in,data_ref);
+            [data_out_hd,LLR,GMI] = obj.process_(data_in,data_ref);
+            try
+                data_out_hd = reshape(data_out_hd,shape_in(1),shape_in(2));
+            catch
+                warning('output reshaping failed after MLSE');
+            end
 
             signalclass_hd = signalclass;
             signalclass_hd.signal = data_out_hd;
 
-            signalclass_sd = signalclass;
-            signalclass_sd.signal = data_out_sd;
+            % signalclass_sd = signalclass;
+            % signalclass_sd.signal = data_out_sd;
 
         end
 
-        function [VITERBI_ESTIMATION_SYMBOLS,soft_decisions] = process_(obj,data_in,data_ref)
+        function [VITERBI_ESTIMATION_SYMBOLS,LLR_maxlogmap,GMI] = process_(obj,data_in,data_ref)
 
+            
+            arguments
+                obj
+                data_in
+                data_ref
+            end
+
+            debug                   = obj.debug;
+
+            trellis_state_mode      = obj.trellis_state_mode;   % General: States should match the target states of the prev. EQ (EQ's job was to reduce the error between signal and the target)
+                                            % 0 = use provided states (MUST provide the correct states); 
+                                            % 1 = normalize to = 1 rms;
+                                            % 2 = use target symbols; 
+                                            % 3 = use statistical levels
+                                            % 3 analyzes avg of rx signal levels - can help with nonlinear impairments
+            
+            trellis_exclusion       = obj.trellis_exclusion; % PAM-6 only (only if data is NOT precoded!)
+            
+            scale_mode              = obj.scale_mode; % scale_mode: 
+                                            % 0 = no scaling, 
+                                            % 1 = RMS→scale MODEL, 
+                                            % 2 = MMSE/time-corr→scale MODEL, 
+                                            % 3 = RMS→scale DATA,  
+                                            % 4 = MMSE/time-corr→scale DATA
+                
+            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+            %%%%%% PREPARATIONS %%%%%%%%
 
             % remove unnecessary zeros at start of impulse response to keep
             % number of trellis states minimal
@@ -63,399 +106,557 @@ classdef MLSE < handle
                 obj.DIR = [0 obj.DIR];
             end
 
+            % impulse respnse to remove from signal
+            obj.DIR = flip(obj.DIR); %i.e. -0.2676   -0.0478    1.0000
 
-            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-            %%%%%% PREPARATIONS %%%%%%%%
+            tx_bits = PAMmapper(obj.M,0,"eth_style",0).demap(data_ref);
 
-            %%%% Separate the equalized signal into the respective levels based on the actually transmitted level
-            constellation    = unique(data_ref);
-            decisionLevels   = (constellation(1:end-1) + constellation(2:end)) / 2;
-            tx_bits = PAMmapper(numel(constellation),0,"eth_style",1).demap(data_ref);
+            % Trellis States
+            if trellis_state_mode == 1 % Normalize the Trellis states to =1 RMS
 
+                obj.trellis_states = obj.trellis_states ./ rms(obj.trellis_states);
 
-            % impulse respnse i.e. [0.5, 1.0000]
-            obj.DIR = flip(obj.DIR);
+            elseif trellis_state_mode == 2 %simply use the states from the ref signal (should be a robust option)
+
+                obj.trellis_states = reshape(unique(data_ref),1,length(unique(data_ref)));
+
+            elseif trellis_state_mode == 3 %use_statistical_levels
+
+                %%%% Separate the equalized signal into the respective levels based on the actually transmitted level
+                constellation    = unique(data_ref);
+
+                % find actual levels from rx signal
+                symbols_for_lvl = NaN(numel(constellation),length(data_ref));
+                for l = 1:numel(constellation)
+                    level_amplitude = constellation(l);
+                    symbols_for_lvl(l,data_ref==level_amplitude) = data_in(data_ref==level_amplitude);
+                end
+
+                %replace the trellis states
+                avg_levels = mean(symbols_for_lvl,2,'omitnan');
+                obj.trellis_states = sort(avg_levels)';
+
+                %also replace the whole ref signal (PAM-M) levels
+                [~, idx] = ismember(data_ref, unique(data_ref));
+                data_ref = avg_levels(idx);
+
+            end
 
-            % RMS normalization of input data
-            data_in = data_in ./ rms(data_in);
 
             % seems to be the only way to use combvec for a flexible amount
             % of vectors. 'combs' contains all trellis states
             pre_comb_mat = repmat(obj.trellis_states,length(obj.DIR)-1,1);
             pre_comb_cell = mat2cell(pre_comb_mat,ones(1,size(pre_comb_mat,1)),size(pre_comb_mat,2));
             combs = fliplr(combvec(pre_comb_cell{:}).');
+            first_sym = combs(:,1); % das ist das älteste/ trailing Symbol aus der sequenz
+            last_sym = combs(:,end); %hiermit wird entschieden/ das ist das cursor symbol am ende der sequenz
 
-            % Save first and last symbol of each state
-            first_sym = combs(:,1);
-            last_sym = combs(:,end);
             states = sum(combs,2);
+            nStates = length(last_sym);
 
-            % Calculate all possible input symbols for the desired impulse
-            % response. Row number is the index of the previous state,
-            % column number is the index of the next state
-            % noise free received == branch metrics
-            noise_free_received = zeros(length(states),length(states));
-            count_row = 1;
-            count_col = 1;
-            for l1 = 1:length(states)
-                for l2 = 1:length(states)
-                    if sum(combs(l2,2:end) == combs(l1,1:end-1)) == size(combs,2)-1
-                        noise_free_received(count_row,count_col) = sum(combs(l2,:).*obj.DIR(end:-1:2)) + last_sym(l1)*obj.DIR(1);
-                    else
-                        noise_free_received(count_row,count_col) = inf;
+            % % Calculate all possible input symbols for the desired impulse
+            % % response. Row number is the index of the previous state,
+            % % column number is the index of the next state
+            % % noise free received == branch metrics
+            % assumes: last_sym = combs(:,end);  % already defined earlier
+            levels = sort(unique(obj.trellis_states(:)).');
+            edges  = [levels(1) levels(end)];    % edge levels (0 and 5 in PAM6)
+
+            noise_free_received = inf(nStates,nStates);   % rows: to, cols: from
+            edge_edge_mask      = false(nStates,nStates); % rows: to, cols: from
+
+            for from = 1:nStates
+                for to = 1:nStates
+                    % valid transition if shift-register overlap holds
+                    if all(combs(to,2:end) == combs(from,1:end-1))
+                        % noiseless sample for the 'to' state reached from 'from'
+                        noise_free_received(to,from) = ...
+                            dot(combs(to,:), obj.DIR(end:-1:2)) + last_sym(from)*obj.DIR(1);
+
+                        % mark edge→edge candidate (to be excluded only on even→odd steps)
+                        edge_edge_mask(to,from) = ...
+                            (last_sym(from)==edges(1) || last_sym(from)==edges(2)) && ...
+                            (last_sym(to)  ==edges(1) || last_sym(to)  ==edges(2));
                     end
-                    count_row = count_row + 1;
                 end
-                count_col = count_col + 1;
-                count_row = 1;
             end
 
-            % match amplitude levels of input signal to those of the calculated ideal symbols
-            % i.e. match the rms values of data_in to noise_free_received
-            if isreal(data_in)
-                if obj.M == round(obj.M)
-                    data_in = data_in * rms(noise_free_received(noise_free_received ~= inf),'all','omitnan');
-                end
+            h = flip(obj.DIR(:)).';
+            data_in = data_in(:);
+            y_ideal = conv(data_ref(:), h, "same");
+
+            switch scale_mode
+                case 0
+                    g = 1; b = 0;
+                case 1  % RMS: scale model to data
+                    g = rms(data_in)/rms(y_ideal);  b = mean(data_in) - g*mean(y_ideal);
+                case 2  % MMSE/time-corr: scale states to data
+                    [c,lags] = xcorr(data_in(:), y_ideal, 64);
+                    [~,ix]   = max(abs(c));
+                    lag      = lags(ix);
+                    y_ideal  = circshift(y_ideal, lag);
+                    mu_y = mean(data_in(:));
+                    mu_i = mean(y_ideal);
+                    y_c = data_in(:)-mu_y;
+                    yi_c = y_ideal-mu_i;
+                    g = (yi_c'*y_c)/(yi_c'*yi_c);
+                    b = mu_y - g*mu_i;
+                case 3  % RMS flipped: scale data to model
+                    gd = rms(y_ideal)/rms(data_in); bd = mean(y_ideal) - gd*mean(data_in);
+                    data_in = gd*data_in + bd;
+                    g = 1; b = 0;
+                case 4  % MMSE/time-corr flipped: scale data to states
+                    [c,lags] = xcorr(data_in(:), y_ideal(:), 64);
+                    [~,ix]   = max(abs(c));
+                    lag      = lags(ix);
+                    y_ideal  = circshift(y_ideal(:), lag);
+                    mu_y = mean(data_in(:));
+                    mu_i = mean(y_ideal);
+                    y_c  = data_in(:) - mu_y;    % data_in centered
+                    yi_c = y_ideal      - mu_i;  % ideal   centered
+                    g = (y_c' * yi_c) / (y_c' * y_c);
+                    b = mu_i - g * mu_y;
+                    data_in = g * data_in(:) + b;
+                    g = 1; b = 0;
+            end
+
+            % apply (g,b) to model and compute common sigma
+            noise_free_received = g*noise_free_received + b;
+            last_sym            = g*last_sym + b;
+            sigma2              = mean(abs(data_in - (g*y_ideal + b)).^2);
+            inv2s2              = 1/(2*sigma2);
+
+            if debug
+                figure(100); clf; hold on
+                showLevelScatter(data_in, data_ref, "fignum", 100);
+                yline(noise_free_received(:), 'DisplayName','Transition States','Color','red','HandleVisibility','off');
+                yline(obj.trellis_states(:), 'DisplayName','Transition States','Color','green','LineWidth',2,'HandleVisibility','off')
             end
 
             %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-            %%%%% FORWARD PASS %%%%%
+            %%%%% FORWARD PASS (VITERBI -Alpha's) %%%%%
 
             % Initialize the output vector
-            pm      = zeros(length(states),length(states));
-            bm_fw   = zeros(length(states),length(states),length(data_in));
+            pm      = zeros(nStates,nStates);
+            bm_fw   = zeros(nStates,nStates,length(data_in));
+
+            % first start is evaluated without ISI/ wihout the full Impulse response
+            % so simply use the constellation here
+            bm =  -(data_in(1) - last_sym).^2 * inv2s2;
+            pm = pm + bm;
+            [alpha(:,1),pm_survivor_fw_idx(:,1)] = max(pm,[],2);
+            pm = repmat(alpha(:,1).',nStates,1);
+            bm_fw(:,:,1) = pm;
 
             % Forward Recursion (FSM Computation)
-            for n = 1:length(data_in)
+            for n = 2:length(data_in)
 
+                bm =  -(data_in(n) - noise_free_received).^2 * inv2s2;
+
+                % exclude edge→edge transitions only for even->odd steps && PAM-6
+                if mod(n,2) == 0 && obj.M == 6 && trellis_exclusion
+                    bm(edge_edge_mask) = -Inf;
+                end
 
-                bm =  abs(data_in(n) - noise_free_received).^2;
                 pm = pm + bm;
-                [pm_survivor_fw(:,n),pm_survivor_fw_idx(:,n)] = min(pm,[],2);         % choose lowest path metric as new state
-                pm = repmat(pm_survivor_fw(:,n).',length(states),1);               % update pm (chosen state to 2nd dimension -> FROM state)
+                [alpha(:,n),pm_survivor_fw_idx(:,n)] = max(pm,[],2);        % choose lowest path metric as new state (get min distance for all state transitions towards a new state)
+                pm = repmat(alpha(:,n).',nStates,1);                        % update pm (chosen state to 2nd dimension -> FROM state)
 
                 bm_fw(:,:,n) = bm;
 
             end
 
             % we can now get the best path as min
-            best_fw_path = NaN(1,length(data_in)+1);
+            viterbi_path = NaN(1,length(data_in));
 
             % find ideal trellis path by going through the trellis backwards
-            [~,best_fw_path(length(data_in)+1)] = min(pm_survivor_fw(:,length(data_in)));
+            [~,viterbi_path(length(data_in))] = max(alpha(:,length(data_in)));
+            for n = length(data_in):-1:2
+                viterbi_path(n-1) = pm_survivor_fw_idx(viterbi_path(n),n);
+            end
+
+
+            if debug
+                % alpha_ = alpha - min(alpha) + eps;
+                % figure();hold on;
+                % n = 10;
+                % scatter(1:n,obj.trellis_states(repmat([1:numel(obj.trellis_states)]',1,n)),abs(alpha_(:,end-n+1:end)),'Marker','o','LineWidth',1);
+                % scatter(1:n,obj.trellis_states(viterbi_path(end-n+1:end)),500,'Marker','x','LineWidth',1,'MarkerEdgeColor','green');
+                % % scatter(1:n,data_ref(end-n+1:end),500,'Marker','x','LineWidth',1,'MarkerEdgeColor','red');
+                % yticks(obj.trellis_states);
+                % ylim([min(obj.trellis_states)-1 max(obj.trellis_states)+1]);
+
+                %% ----- Build true path from known sequence (data_ref) -----
+                % Memory length used by VA (L = length(obj.DIR)-1 symbols stored in state)
+                L = size(combs,2);                 % each row of combs is the L-tap state vector
+                N = length(data_in);
+
+                % Map each trellis level to an index (1..M)
+                [~, level_to_idx] = ismember(levels, levels); %#ok<ASGLU>  % identity map
+                [ok_ref, ref_idx] = ismember(data_ref(:).', levels);
+                if ~all(ok_ref)
+                    warning('Some data_ref symbols are not in "levels". True-path build may fail.');
+                end
+
+                % Precompute next_state(from_state, u_idx) LUT such that:
+                % combs(next_state,:) == [ combs(from_state,2:end) , levels(u_idx) ]
+                next_state = zeros(size(combs,1), numel(levels), 'uint32');
+                for from = 1:size(combs,1)
+                    prefix = combs(from,2:end);  % what must match in 'to' for a valid transition
+                    for ui = 1:numel(levels)
+                        target = [prefix, levels(ui)];
+                        % find the unique 'to' whose state vector equals target
+                        to = find(all(bsxfun(@eq, combs, target), 2), 1, 'first');
+                        if isempty(to), to = 0; end
+                        next_state(from, ui) = to;
+                    end
+                end
+
+                % Initialize true path at n=1: pick any state with last_sym == data_ref(1)
+                % Prefer one whose suffix matches the first available history if L>1.
+                cand = find(last_sym == data_ref(1));
+                if isempty(cand)
+                    % fallback: choose closest in amplitude (should not happen if levels match)
+                    [~,ix] = min(abs(last_sym - data_ref(1)));
+                    cand = ix;
+                end
+                true_state_path = zeros(1,N,'uint32');
+                true_state_path(1) = cand(1);
+
+                % Propagate forward using the known inputs data_ref(n)
+                for n = 2:N
+                    ui = ref_idx(n);     % index of the actual transmitted level at time n
+                    from = true_state_path(n-1);
+                    if from==0 || ui==0
+                        true_state_path(n) = 0;
+                    else
+                        true_state_path(n) = next_state(from, ui);
+                        if true_state_path(n)==0
+                            % Safety fallback: if no valid transition found (should not happen)
+                            % choose any to-state whose vector matches shift+current symbol
+                            target = [combs(from,2:end), levels(ui)];
+                            to = find(all(bsxfun(@eq, combs, target), 2), 1, 'first');
+                            if isempty(to), to = from; end
+                            true_state_path(n) = to;
+                        end
+                    end
+                end
+
+                %% ----- Collect branch metrics along decoded vs. true path -----
+                bm_decoded = nan(1,N);
+                bm_true    = nan(1,N);
+
+                % n=1 in your code stores pm into bm_fw(:,:,1); real BMs start at n>=2
+                for n = 2:N
+                    % Decoded path: to = viterbi_path(n), from = survivor that fed it
+                    to_d   = viterbi_path(n);
+                    from_d = pm_survivor_fw_idx(to_d, n);
+                    bm_decoded(n) = bm_fw(to_d, from_d, n);
+
+                    % True path: transition true_state_path(n-1) -> true_state_path(n)
+                    to_t   = true_state_path(n);
+                    from_t = true_state_path(n-1);
+                    if to_t>0 && from_t>0
+                        bm_true(n) = bm_fw(to_t, from_t, n);
+                    end
+                end
+
+                % Convert to "cost" for intuitive plotting (your BM is a log-likelihood)
+                cost_dec = -bm_decoded;
+                cost_true= -bm_true;
+
+                %% ----- Plot a short window for clarity -----
+                win = max(2, N-20000):N;  % last 200 samples (adjust as needed)
+                figure('Color','w'); hold on; grid on; box on;
+                plot(win, cost_true(win),  'LineWidth',1.2, 'DisplayName','True path cost (−BM)');
+                plot(win, cost_dec(win),   'LineWidth',1.2, 'DisplayName','Decoded path cost (−BM)');
+                xlabel('Time index n'); ylabel('Branch cost'); title('Branch metrics along true vs decoded path');
+                legend('Location','best');
+
+                %% ----- Optional: overlay symbol levels for the same window -----
+                yyaxis right
+                plot(win, data_in(win), ':', 'LineWidth',0.8, 'DisplayName','y(n)');
+                ylabel('Amplitude');
+
+            end
+
+            VITERBI_ESTIMATION_SYMBOLS(1:length(data_in)) = first_sym(viterbi_path);
 
-            
             %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-            %%%%% BACKWARD PASS %%%%%
+            %%%%% BACKWARD PASS (Beta's) %%%%%
 
             % Initialize the output vector
-            pm                 = zeros(length(states),length(states));
-            pm_survivor_bw     = zeros(length(states),length(data_in)+1);
-            pm_survivor_bw_idx = zeros(length(states),length(data_in)+1);
-            bm_bw   = zeros(length(states),length(states),length(data_in)+1); 
+            pm       = zeros(nStates,nStates);
+            beta     = zeros(nStates,length(data_in));
+            pm_survivor_bw_idx = zeros(nStates,length(data_in));
+            bm_bw   = zeros(nStates,nStates,length(data_in));
 
             % starting with the state that has the lowest sum path
             % metric, follow the stored information about the
             % predecessor
-            for h = length(data_in):-1:1
+            for h = length(data_in)-1:-1:1
 
-                best_fw_path(h) = pm_survivor_fw_idx(best_fw_path(h+1),h);
+                bm =  -(data_in(h+1) - noise_free_received).^2 * inv2s2;
+
+                % exclude edge→edge transitions only for even->odd steps && PAM-6
+                if mod(h+1, 2) == 0 && obj.M == 6 && trellis_exclusion
+                    bm(edge_edge_mask) = -Inf;
+                end
 
-                bm =  abs(data_in(h) - noise_free_received).^2;
                 pm = pm + bm.';
-                [pm_survivor_bw(:,h),pm_survivor_bw_idx(:,h)] = min(pm,[],2);      % choose lowest path metric as new state
-                pm = repmat(pm_survivor_bw(:,h).',length(states),1);               % update pm (chosen state to 2nd dimension -> FROM state)
+                [beta(:,h),pm_survivor_bw_idx(:,h)] = max(pm,[],2);      % choose lowest path metric as new state
+                pm = repmat(beta(:,h).',nStates,1);               % update pm (chosen state to 2nd dimension -> FROM state)
 
                 bm_bw(:,:,h) = bm;
 
             end
 
-            VITERBI_ESTIMATION_IDX(1:length(data_in)) = first_sym(best_fw_path(2:end));
-            VITERBI_ESTIMATION_SYMBOLS(1:length(data_in)) = constellation(best_fw_path(2:end));
-
             %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-            %%%%% FORWARD PASS %%%%%
-            
+            %%%%% FORWARD PASS PAM 2,4,8 (Combine Alpha and Beta to yield LLP's) %%%%%
+
             %calc the log probabilities (llp's)
 
-            for k = 2:length(data_in)
-                
-                fsm     = repmat(pm_survivor_fw(:,k-1)',[length(states),1]);   %pm_survivor_fw size: length(states)xlength(sequence)
-                bm      = bm_fw(:,:,k);                                         %bm_fw size: length(states)xlength(states)xlength(sequence)
-                bsm     = repmat(pm_survivor_bw(:,k)',[length(states),1]);       %pm_survivor_bw size: length(states)xlength(sequence)+1
+            for k = 1:length(data_in)
+
+                if k == 1
+
+                    alpha_     = repmat(alpha(:,k)',[nStates,1])';
+                    beta_      = beta(:,k);
+
+                    LLP(:,k) = max(alpha_ + beta_,[],2);
+
+                else
+
+                    alpha_     = repmat(alpha(:,k-1)',[nStates,1])';
+                    gamma_     = bm_fw(:,:,k)';
+                    beta_      = beta(:,k);
+
+                    LLP(:,k) = max(alpha_ + gamma_,[],1) + beta_';
+                end
 
-                llp(:,k) = min(fsm+bm+bsm,[],2);
-                
             end
 
-            % Compute soft output PAM4 stream from the metric_sym
-            soft_output = zeros(length(data_in),1);      % Expected symbol value per stage
-            symbol_prob = zeros(length(data_in), length(states)); % Store full probability distribution (optional)
-            llp = llp';
-            for n = 1:length(data_in)
-                metrics = llp(n, :);   % A posteriori metric for each PAM4 candidate
-                % For numerical stability, subtract the maximum metric before exponentiating
-                maxMetric = min(metrics);
-                expMetrics = exp(metrics - maxMetric);
-                probs = expMetrics / sum(expMetrics);  % Normalize to get probabilities
-                symbol_prob(n, :) = probs;              % (Optional) store distribution for analysis
-                % Compute the soft output as the expected value of the PAM4 symbols
-                soft_output(n) = sum(probs .* constellation');
+            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+            %%%%% FORWARD PASS PAM2,4,8 %%%%%
+
+            % These are interchangeable... second is chatgpt:
+            nml_LLP = LLP - max(LLP); %subtract highest value for better numerical stability, LLP's are not always close to zero
+            expLLP = exp(nml_LLP);
+            state_prob = expLLP ./ sum(expLLP); % sums to one (or numerically close to one)
+
+            % compute symbol‐posteriors from LLP in the log‐domain:
+            amax = max(LLP,[],1);
+            logZ    = amax + log(sum(exp(LLP - amax), 1));
+            logPstate = LLP - logZ;          % still in log‐domain
+            state_prob    = exp(logPstate);      % exact, sums to 1
+
+
+            if obj.M == 6
+
+                num_bits = 5;
+
+                % all possible transitions (for now 36, including the "edges"
+                % of the QAM 32 constellation)
+                states = PAMmapper(6,0,"eth_style",0).levels;
+                pam6transitions = combvec(states,states)'; % pam6transitions =
+                % [-5    -5;
+                % -3    -5;
+                % -1    -5; ...
+
+                pam6bits = PAMmapper(6,0,"eth_style",0).demap(reshape(pam6transitions',[],1)./sqrt(10));
+
+                pam6bits = reshape(pam6bits',5,[])';
+
+                [ok1, idx_sym_1] = ismember(pam6transitions(:,1), states);
+                [ok2, idx_sym_2] = ismember(pam6transitions(:,2), states);
+                assert(all(ok1)&all(ok2), 'Some transition amplitude not found in trellis_states')
+                pam6ind = [idx_sym_1, idx_sym_2];
+
+                tx_bits_pam6_reshaped = reshape(tx_bits,5,[])'; % N x 5
+
+                numPairs = floor(size(LLP,2)/2);
+                LLR_exact     = zeros(numPairs,5);
+                LLR_maxlogmap = zeros(numPairs,5);
+
+                for k = 1:numPairs
+                    symbol1 = 2*k-1;
+                    symbol2 = 2*k;
+
+                    LLP1 = LLP    (:,symbol1);
+                    LLP2 = LLP    (:,symbol2);
+                    prob1 = state_prob(:,symbol1);
+                    prob2 = state_prob(:,symbol2);
+
+                    % 36 joint‐metrics M = log P(i)*P(j) = L1(i)+L2(j)
+                    Mij = LLP1(pam6ind(:,1)) + LLP2(pam6ind(:,2));
+                    pij = prob1(pam6ind(:,1)) .* prob2(pam6ind(:,2));
+
+                    % now for each of the 5 bits do exact-LLR or max-log
+                    for b = 1:num_bits
+                        idx_sym_1 = pam6bits(:,b)==1;
+                        idx_bit_1 = pam6bits(:,b)==0;
+
+                        %--- exact LLR from probabilities
+                        P1 = sum(pij(idx_sym_1));
+                        P0 = sum(pij(idx_bit_1));
+                        LLR_exact(k,b) = log(P1./P0);
+
+                        %--- max-log:
+                        LLR_maxlogmap(k,b) = max( Mij(idx_sym_1) ) - max( Mij(idx_bit_1) ); % N x num_bits
+                    end
+                end
+
+                MI = zeros(1, num_bits);
+                for k = 1:num_bits
+
+                    idx_bit_1 = (tx_bits_pam6_reshaped(:,k) == 0); %wo sind die 1en
+                    idx_sym_1 = (tx_bits_pam6_reshaped(:,k) == 1); %wo sind die 0en
+
+                    %LLR's for all actually transmitted ones or zeros
+                    llr0 = LLR_exact(idx_bit_1,k);
+                    llr1 = LLR_exact(idx_sym_1,k);
+
+                    % Calculate mutual information for bit position k
+                    I0 = mean(log2(1 + exp(llr0))); % exp(--LLR) = exp(positive) > 1
+                    I1 = mean(log2(1 + exp(-llr1))); % exp(-+LLR) = exp(negative) < 1
+                    MI(k) = 1 - 0.5 * (I0 + I1);
+                end
+
+                GMI = sum(MI);  % Total mutual information per symbol
+                GMI = GMI/2;
+
+            else
+
+                % Number of symbols and bits per symbol
+                num_bits = log2(length(obj.trellis_states)); % 2 bits per symbol
+                % bit_mapping = PAMmapper(length(obj.trellis_states),0,"eth_style",0).demap(first_sym./rms(first_sym));
+                bit_mapping = PAMmapper(length(obj.trellis_states),0,"eth_style",0).showBitMapping;
+                % Initialize LLR storage
+                LLR_maxlogmap = zeros(length(data_in),num_bits);
+                LLR_exact = zeros(length(data_in),num_bits);
+
+                % Compute bit-wise LLRs
+                for bit_idx = 1:num_bits
+
+                    % Find indices where bit is 0 and where it is 1
+                    idx_bit_0 = bit_mapping(:,bit_idx) == 0;
+                    idx_bit_1 = bit_mapping(:,bit_idx) == 1;
+
+                    % Sum over log-probabilities
+                    % Max-Log approximation: using max instead of sum)
+                    LLR_maxlogmap(:,bit_idx) = max(LLP(idx_bit_1,:), [], 1) - max(LLP(idx_bit_0,:), [], 1);
+
+                    % Sum probabilities over states for which the bit is 1 and 0, respectively.
+                    P0 = sum(state_prob(idx_bit_0, :),1);
+                    P1 = sum(state_prob(idx_bit_1, :),1);
+                    LLR_exact(:,bit_idx) = log(P1./P0); % N x num_bits
+
+
+                end
+
+                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+                %%%%% CALC NGMI %%%%%
+
+                MI = zeros(1, num_bits);
+                LLR_exact = LLR_exact;
+                for k = 1:num_bits
+
+                    idx_bit_0 = (tx_bits(:,k) == 0); %wo sind die 1en
+                    idx_bit_1 = (tx_bits(:,k) == 1); %wo sind die 0en
+
+                    %LLR's for all actually transmitted ones or zeros
+                    llr0 = LLR_exact(idx_bit_0,k);
+                    llr1 = LLR_exact(idx_bit_1,k);
+
+                    % Calculate mutual information for bit position k
+                    I0 = mean(log2(1 + exp(llr0))); % exp(--LLR) = exp(positive) > 1
+                    I1 = mean(log2(1 + exp(-llr1))); % exp(-+LLR) = exp(negative) < 1
+                    MI(k) = 1 - 0.5 * (I0 + I1);
+                end
+
+                GMI = sum(MI);  % Total mutual information for 2 symbols
+
             end
 
-            % Number of symbols and bits per symbol
-            num_symbols = constellation;
-            num_bits = 2; % 2 bits per symbol
+            VITERBI_ESTIMATION_SYMBOLS = VITERBI_ESTIMATION_SYMBOLS./rms(VITERBI_ESTIMATION_SYMBOLS);
 
-            bit_mapping = PAMmapper(4,0,"eth_style",1).showBitMapping; % Each row corresponds to the symbol above
+            if debug
+                %%% DEBUG PLOT LIKELIHOOD RATIOS %%%
+                figure(115);clf
+                subplot(2,1,1)
+                for bit = 1:num_bits
+                    hold on;
+                    histogram(LLR_exact(:,bit),1000,"DisplayName",sprintf('Actual LLR of Bit Pos %d',bit),'LineStyle','none','FaceAlpha',0.4);
+                end
+                legend
 
-            % Initialize LLR storage
-            llr = zeros(num_bits, length(data_in));
-
-            % Compute bit-wise LLRs
-            for bit_idx = 1:num_bits
-                % Find indices where bit is 0 and where it is 1
-                idx_bit_0 = find(bit_mapping(:,bit_idx) == 0);
-                idx_bit_1 = find(bit_mapping(:,bit_idx) == 1);
-
-                % Sum over log-probabilities (Max-Log approximation: using min instead of sum)
-                llr(:,bit_idx) = min(llp(:,idx_bit_1), [], 1) - min(llp(:,idx_bit_0), [], 1);
-               
+                subplot(2,1,2)
+                for bit = 1:num_bits
+                    hold on;
+                    histogram(LLR_maxlogmap(:,bit),1000,"DisplayName",sprintf('Max Log LLR of Bit Pos %d',bit),'LineStyle','none','FaceAlpha',0.4);
+                end
+                legend
             end
 
-            % Convert LLR values to a hard-decision bit stream
-            bit_stream = llr < 0;
-            [~,~,ber_llr,~] = calc_ber(bit_stream',tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
-            fprintf('LLR BER   : %.2e \n',ber_llr);
-
-            % directly decide based on lowest LLP index
-            [~,llp_based_state_seq]=min(llp);
-            LLP_EST(1:length(data_in)) = constellation(llp_based_state_seq);
-            rx_bits = PAMmapper(numel(constellation),0,"eth_style",1).demap(LLP_EST');
-            [~,~,ber_llp,~] = calc_ber(rx_bits,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
-            fprintf('LLP BER   : %.2e \n',ber_llp);
-            % [~,~,ber_llp,~] = calc_ber(circshift(rx_bits,1),tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
-            % fprintf('LLP BER +1: %.2e \n',ber_llp);
-            % [~,~,ber_llp,~] = calc_ber(circshift(rx_bits,-1),tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
-            % fprintf('LLP BER -1: %.2e \n',ber_llp);
-            
-            %%%% DECIDE based on Viterbi traceback
-            rx_bits = PAMmapper(numel(constellation),0,"eth_style",1).demap(VITERBI_ESTIMATION_SYMBOLS');
-            [~,~,ber_viterbi,~] = calc_ber(rx_bits,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
-            fprintf('Viterbi BER: %.2e \n',ber_viterbi);
-            % [~,~,ber_viterbi,~] = calc_ber(circshift(rx_bits,1),tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
-            % fprintf('Viterbi BER: %.2e \n',ber_viterbi);
-
-            % directly decide based on the FW path metrics
-            [~,fw_direct_state_seq]=min(pm_survivor_fw);
-            FW_EST(1:length(data_in)) = constellation(fw_direct_state_seq);
-            rx_bits = PAMmapper(numel(constellation),0,"eth_style",1).demap(FW_EST');
-            [~,~,ber_fw,~] = calc_ber(rx_bits,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
-            fprintf('FW BER: %.2e \n',ber_fw);
-
-            % directly decide based on the BW path metrics
-            [~,bw_direct_state_seq]=min(pm_survivor_bw);
-            BW_EST(1:length(data_in)) = constellation(bw_direct_state_seq(2:end));
-            rx_bits = PAMmapper(numel(constellation),0,"eth_style",1).demap(BW_EST');
-            [~,~,ber_bw,~] = calc_ber(rx_bits,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
-            fprintf('BW BER: %.2e \n',ber_bw);
-            % [~,~,ber_viterbi,~] = calc_ber(circshift(rx_bits,1),tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
-            % fprintf('BW BER: %.2e \n',ber_viterbi);
-            % [~,~,ber_viterbi,~] = calc_ber(circshift(rx_bits,-1),tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
-            % fprintf('BW BER: %.2e \n',ber_viterbi);
-
-            PAMmapper(4,0,"eth_style",1).showBitMapping
 
 
             %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-            tx_symbolpos     = zeros(numel(constellation),length(data_ref));
-            est_symbolpos    = zeros(numel(constellation),length(data_ref));
-            for lvl = 1:numel(constellation)
-                tx_symbolpos(lvl,data_ref==constellation(lvl)) = 1;
-                est_symbolpos(lvl,VITERBI_ESTIMATION_IDX==first_sym(lvl)) = 1;
-                est_symbolpos_llm(lvl,llp_based_state_seq==lvl) = 1;
+            %%%%% CHECK BER's %%%%%
+
+            if debug
+                tx_bits = reshape(tx_bits',[],1);
+                fprintf('\n')
+                disp('Start DEBUG MLSE:')
+                % DECIDE based on Viterbi traceback
+                VITERBI_ESTIMATION_SYMBOLS = VITERBI_ESTIMATION_SYMBOLS./rms(VITERBI_ESTIMATION_SYMBOLS);
+                rx_bits = PAMmapper(obj.M,0,"eth_style",0).demap(VITERBI_ESTIMATION_SYMBOLS');
+                rx_bits = reshape(rx_bits',[],1);
+                [~,numErr,ber_viterbi,~] = calc_ber(rx_bits,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
+                fprintf('Viterbi BER: %.2e \n',ber_viterbi);
+                fprintf('Viterbi Errors = %d\n', numErr);
+
+                pairs = reshape(VITERBI_ESTIMATION_SYMBOLS,2,[]).';
+                levels = sort(unique(VITERBI_ESTIMATION_SYMBOLS));
+                isedge = ismember(pairs, [levels(1) levels(end)]);
+                isforbidden = sum(isedge,2)==2;
+                fprintf('Found %d forbidden transitions (even→odd edges).\n', nnz(isforbidden));
+
+
+                % Convert LLR values to a hard-decision bit stream
+                bit_stream = LLR_maxlogmap > 0; %ratio separates lower or higher than =0 -> simply decode for the negative values
+                bit_stream = reshape(bit_stream',[],1);
+                [~,~,ber_llr,~] = calc_ber(bit_stream,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
+                fprintf('LLR MaxLogMAP BER   : %.2e \n',ber_llr);
+
+                % Convert LLR values to a hard-decision bit stream
+                bit_stream = LLR_exact > 0; %ratio separates lower or higher than =0 -> simply decode for the negative value
+                bit_stream = reshape(bit_stream',[],1);
+                [~,~,ber_llr,~] = calc_ber(bit_stream,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
+                fprintf('LLR BER   : %.2e \n',ber_llr);
+
+                % DECIDE based on lowest LLP index and check BER
+                [~,llp_based_state_seq]=max(LLP);
+                LLP_EST(1:length(data_in)) = first_sym(llp_based_state_seq);
+                LLP_EST = LLP_EST./rms(LLP_EST);
+                rx_bits = PAMmapper(obj.M,0,"eth_style",0).demap(LLP_EST');
+                rx_bits = reshape(rx_bits',[],1);
+                [~,~,ber_llp,~] = calc_ber(rx_bits,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
+                fprintf('LLP BER   : %.2e \n',ber_llp);
+
+                % directly decide based on the FW path metrics
+                [~,fw_direct_state_seq]=max(alpha);
+                FW_EST(1:length(data_in)) = first_sym(fw_direct_state_seq);
+                FW_EST = FW_EST./rms(FW_EST);
+                rx_bits = PAMmapper(obj.M,0,"eth_style",0).demap(FW_EST');
+                rx_bits = reshape(rx_bits',[],1);
+                [~,~,ber_fw,~] = calc_ber(rx_bits,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
+                fprintf('FW BER   : %.2e \n',ber_fw);
+                disp('Stop DEBUG MLSE:')
+                fprintf('\n')
             end
 
-            est_symbolpos= logical(est_symbolpos);
-            tx_symbolpos = logical(tx_symbolpos);
-            est_symbolpos_llm = logical(est_symbolpos_llm);
-
-            figure(299);clf;hold on;
-            llp_filt = NaN(length(states),length(data_in));
-            llp_ = llp - min(llp,[],1);
-            for c = 2%:numel(constellation)
-                for c2 = 1:3%numel(constellation)
-
-                    idx = est_symbolpos_llm(c,:);
-                    llp_filt(c2,idx) = (llp(c2,idx));
-
-                    plotidx = 2:200000;
-                    scatter(plotidx,llp_filt(c2,plotidx),1,'o','LineWidth',1,'DisplayName',sprintf('LLP of Symbol %d | %d was transmitted',c2,c));
-
-                end
-            end
-
-            % Assume llp is a 4 x N matrix (4 PAM-4 symbols, N time steps)
-            [numSymbols, numTimeSteps] = size(llp);
-            P_symbols = zeros(numSymbols, numTimeSteps); % To hold the soft output probabilities
-
-            for k = 1:numTimeSteps
-
-                % Compute the exponentials. (Use -llp_shifted if llp's are costs.)
-                exponents = -llp_(:, k);
-
-                % Normalize to form a probability vector.
-                P_symbols(:, k) = exponents / sum(exponents);
-            end
-
-            % Optionally, if you prefer a single soft output value per symbol (an expectation),
-            % define your PAM-4 constellation levels, e.g.:
-            soft_output = sum(P_symbols .* constellation, 1);  % 1 x N vector of soft outputs
-
-
-            %%%% BITWISE LLR's ??? %%%%%
-            figure(100);
-            clf
-            hold on
-            title('LLR between inner and outer bits')
-            bitmapping = PAMmapper(numel(constellation),0).demap(constellation);
-            for bp = 1:size(bitmapping,2)
-                pos_bitone = find(bitmapping(:,bp)==1);
-                pos_bitzero = find(bitmapping(:,bp)~=1);
-
-                llr_bits(bp,:) = min(llp(pos_bitone,:),[],1) - min(llp(pos_bitzero,:),[],1);
-
-                subplot(size(bitmapping,2),1,bp)
-                scatter(1:length(data_ref),llr_bits(bp,:),1,'.');
-            end
-            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-            % From: Log-Likelihood Probabilities LLP --> To: Log-Likelihood Ratios LLR
-            for s = 1:length(states)-1
-                llr(s,:) = llp_(s,:) - llp_(s+1,:); % subtract llp's of successive const. points to get llr's
-            end
-            
-            % Build Sets of LLR's that contain only those values at
-            % timepoint k, where the symbols:
-            % a) were actually transmitted: set "S"
-            % b) were decoded by viterbi: set "SC"
-            S = NaN(length(states)-1,length(data_in));
-            SC = NaN(length(states)-1,length(data_in));
-
-            llp_inf = llp_;
-            llp_inf(llp_==0) = Inf;
-            [~,scnd_idx] = min(llp_inf,[],1) ;
-            for s = 1:length(states)-1
-
-                %%% SET "S"
-                % find all time indices k, where const point s was actually transmitted
-                indice = tx_symbolpos(s,:)==1;
-                S(s,indice) = llr(s,indice);
-                %calc mean of all llr's where symbol s was transmitted
-                K(s,1) = mean(S(s,indice),'omitnan');
-                
-                % find all time indices k, where const. point s+1 was actually transmitted
-                indice_plusone = tx_symbolpos(s+1,:)==1;
-                S(s,indice_plusone) = llr(s,indice_plusone);
-                K(s,2) = mean(S(s,indice_plusone),'omitnan');
-
-                
-                %%% SET "SC"
-                % find all time indices k, where symbol s was decoded 
-                % idx: find positions in time where a symbol s was decoded
-                idx = est_symbolpos_llm(s,:);
-                % idx 2: find positions where the strongest competitor is
-                % s+1, i.e. the second best llp is at s+1
-                idx2 = scnd_idx == s+1;
-                SC(s,idx&idx2) = llr(s,idx&idx2);
-                KC(s,1) = mean(SC(s,idx&idx2),'omitnan');
-
-                % find all time indices k, where symbol s+1 was decoded 
-                % idx: find positions in time where a symbol s+1 was decoded
-                idx = est_symbolpos_llm(s+1,:);
-                idx2 = scnd_idx == s;
-                % idx 2: find positions where the strongest competitor is
-                % s, i.e. the second best llp is at s
-                SC(s,idx&idx2) = llr(s,idx&idx2);
-                KC(s,2) = mean(SC(s,idx&idx2),'omitnan');
-
-                % scale the sets, using the average (?) of the 
-                llrcn(s,:) = SC(s,:) * ( constellation(s+1)-constellation(s) ) ./ ( K(s,2)-K(s,1)) + decisionLevels(s); 
-
-            end
-            
-            figure(111)
-            title("Bla")
-            clf
-            for s = 1:length(states)-1
-                figure(1111)
-                hold on
-                title(sprintf('llrcn = llp %d - llp %d',s, s+1,s, s+1));
-                scatter(1:length(llrcn),llrcn(s,:),1,'.');
-
-
-                % STUFENARTIGE LLR'S UND LLP'S %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-                figure(111)
-                
-                subplot(1,length(states),s)
-                hold on
-                title(sprintf('llr = llp %d - llp %d \n SC=llr(tx sym = %d or %d)',s, s+1,s, s+1));
-                scatter(1:length(llr),llr(s,:),1,'.');
-                scatter(1:length(SC),SC(s,:),1,'.');
-                yline([K(s,1),K(s,2)]);
-                low = min(min(llr));
-                hi = max(max(llr));
-                ylim([low,hi]);
-
-                subplot(1,length(states),length(states))
-                hold on
-                scatter(1:length(llr),llr(s,:),1,'.');
-                ylim([low,hi]);
-
-
-                % HISTOGRAMME %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-                figure(112)
-                subplot(length(states),1,s)
-                hold on
-                title(sprintf('histogram of llr; between const point %d - %d',s, s+1));
-                histogram(llrcn(s,:),10000,'EdgeAlpha',0)
-                % histogram(SC(s,:),10000,'EdgeAlpha',0)
-                xlim([-20, 20]);
-                subplot(length(states),1,length(states))
-                hold on
-                histogram(llrcn(s,:),10000,'EdgeAlpha',0)
-                % histogram(SC(s,:),10000,'EdgeAlpha',0)
-                % xlim([-20, 20]);
-
-            end
-
-            softdecisions = mean(llrcn,1,'omitnan');
-
-            distance = abs(VITERBI_ESTIMATION_SYMBOLS-llrcn);
-            [win_cost,win_idx] = min(distance,[],1);
-
-            for i = 1:length(data_in)
-                soft_decisions(i) = llrcn(win_idx(i),i);
-            end
-
-            soft_decisions = max(llrcn,[],1);
-
-            showLevelHistogram(soft_decisions,data_ref)
-
-            figure()
-            % scatter(1:length(data_in),VITERBI_ESTIMATION_SYMBOLS,1,'.');
-            scatter(1:length(data_in),soft_decisions,1,'.');
-
-
-            MLM_ESTIMATION(1:length(data_in)) = PAMmapper(numel(constellation),0).quantize(soft_decisions)';
-            rx_bits = PAMmapper(numel(constellation),0).demap(MLM_ESTIMATION');
-            [~,~,ber_mlm,~] = calc_ber(rx_bits,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
-            fprintf('MLM BER: %.2e \n',ber_mlm);
-
         end
-    end
 
-    methods (Access=private)
 
     end
 end

Classes/04_DSP/TransmissionPerformance.m

@@ -42,6 +42,7 @@ classdef TransmissionPerformance
             0.8733, 0.8790, 0.8848, 0.8905, 0.8962, ...
             0.9019, 0.9077, 0.9134, 0.9191, 0.9248, ...
             0.9306, 0.9363, 0.9420, 0.9477, 0.9535];
+
         NGMITHRESHOLDS_SDHD    = [0.8116, 0.8167, 0.8241, 0.8317, 0.8401, ...
             0.8459, 0.8512, 0.8574, 0.8685, 0.8746, ...
             0.8829, 0.8892, 0.8958, 0.9022, 0.9090,...
@@ -60,7 +61,7 @@ classdef TransmissionPerformance
 
         %% LUT for KP4-FEC and Inner Code https://grouper.ieee.org/groups/802/3/dj/public/23_03/patra_3dj_01b_2303.pdf
         
-        CODE_RATE_KP4_AND_INNER = [0.885799];
+        CODE_RATE_KP4_AND_INNER = [0.885799]; %https://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=9979198 -> 199.8 / 224 = 0.891 als code rate 
         BERTHRESHOLDS_KP4_AND_INNER = 4.85e-3; 
 
         % https://www.ieee802.org/3/bs/public/14_11/parthasarathy_3bs_01a_1114.pdf
@@ -153,11 +154,21 @@ classdef TransmissionPerformance
                 netrates.SDHD.Threshold = NaN(1, numMeasurements);
             end
             if ~isempty(ber)
+                netrates.STAIR.GrossRate      = NaN(1, numMeasurements);
+                netrates.STAIR.NetRate      = NaN(1, numMeasurements);
+                netrates.STAIR.CodeRate = NaN(1, numMeasurements);
+                netrates.STAIR.Threshold = NaN(1, numMeasurements);
+
                 netrates.HD.GrossRate      = NaN(1, numMeasurements);
                 netrates.HD.NetRate      = NaN(1, numMeasurements);
                 netrates.HD.CodeRate = NaN(1, numMeasurements);
                 netrates.HD.Threshold = NaN(1, numMeasurements);
 
+                netrates.KP4.GrossRate      = NaN(1, numMeasurements);
+                netrates.KP4.NetRate      = NaN(1, numMeasurements);
+                netrates.KP4.CodeRate = NaN(1, numMeasurements);
+                netrates.KP4.Threshold = NaN(1, numMeasurements);
+
                 netrates.KP4_hamming.GrossRate      = NaN(1, numMeasurements);
                 netrates.KP4_hamming.NetRate      = NaN(1, numMeasurements);
                 netrates.KP4_hamming.CodeRate = NaN(1, numMeasurements);
@@ -200,11 +211,47 @@ classdef TransmissionPerformance
                     end
                     if ~isempty(idxBER)
                         codeRate = obj.CODE_RATES_HD(idxBER);
+                        netrates.STAIR.NetRate(i)      = grossRate(i) * codeRate;
+                        netrates.STAIR.GrossRate(i)    = grossRate(i) ;
+                        netrates.STAIR.CodeRate(i)     = codeRate;
+                        netrates.STAIR.Threshold(i)    = obj.BERTHRESHOLDS_HD(idxBER);
+                    end
+
+                    % HD FEC 3,8e-3
+                    idxBER = [];
+                    for j = length(obj.BERTHRESHOLDS_HDFEC):-1:1
+                        if ber(i) <= obj.BERTHRESHOLDS_HDFEC(j)
+                            idxBER = j;
+                            break;
+                        end
+                    end
+                    if ~isempty(idxBER)
+                        codeRate = obj.CODE_RATE_HDFEC(idxBER);
                         netrates.HD.NetRate(i)      = grossRate(i) * codeRate;
                         netrates.HD.GrossRate(i)    = grossRate(i) ;
                         netrates.HD.CodeRate(i)     = codeRate;
-                        netrates.HD.Threshold(i)    = obj.BERTHRESHOLDS_HD(idxBER);
+                        netrates.HD.Threshold(i)    = obj.CODE_RATE_HDFEC(idxBER);
                     end
+
+
+                    % KP4
+                    idxBER = [];
+                    for j = length(obj.BERTHRESHOLDS_KP4):-1:1
+                        if ber(i) <= obj.BERTHRESHOLDS_KP4(j)
+                            idxBER = j;
+                            break;
+                        end
+                    end
+                    if ~isempty(idxBER)
+                        codeRate = obj.CODE_RATE_KP4(idxBER);
+                        netrates.KP4.NetRate(i)      = grossRate(i) * codeRate;
+                        netrates.KP4.GrossRate(i)    = grossRate(i) ;
+                        netrates.KP4.CodeRate(i)     = codeRate;
+                        netrates.KP4.Threshold(i)    = obj.CODE_RATE_KP4(idxBER);
+                    end
+
+
+                    % KP4 + inner Hamming
                     idxBER = [];
                     for j = length(obj.BERTHRESHOLDS_KP4_AND_INNER):-1:1
                         if ber(i) <= obj.BERTHRESHOLDS_KP4_AND_INNER(j)
@@ -220,6 +267,7 @@ classdef TransmissionPerformance
                         netrates.KP4_hamming.Threshold(i)    = obj.BERTHRESHOLDS_KP4_AND_INNER(idxBER);
                     end
 
+                    % O FEC
                     idxBER = [];
                     for j = length(obj.BERTHRESHOLDS_O_FEC):-1:1
                         if ber(i) <= obj.BERTHRESHOLDS_O_FEC(j)

Classes/DataBaseHandler/DBHandler.m

@@ -5,7 +5,7 @@ classdef DBHandler < handle
 
     properties
         conn        % Database connection object
-        pathToDB    % Path to the SQLite database
+        dataBase    % Path to the SQLite database
         tableNames  % Cell array containing names of all tables in the database
         tables = struct();  % Structure containing MATLAB tables for each database table
         distinctValues
@@ -21,7 +21,7 @@ classdef DBHandler < handle
             %   obj = DBHandler('pathToDB', 'path/to/database.db');
 
             arguments
-                options.pathToDB = ""; % Default value for pathToDB if not provided
+                options.dataBase = ""; % Default value for pathToDB if not provided
                 options.type = "mysql";
             end
 
@@ -37,12 +37,15 @@ classdef DBHandler < handle
             try
 
                 if options.type == "sqlite"
-                    obj.conn = sqlite(obj.pathToDB);
+
+                    obj.conn = sqlite(obj.dataBase);
+
                 elseif options.type == "mysql"
-                    datasource = "jdbc:mysql://134.245.243.254:3306/labor";
+
+                    % datasource = "jdbc:mysql://134.245.243.254:3306/labor";
 
                     obj.conn  = database( ...
-                        "labor", ...                            % Database name
+                        string(obj.dataBase), ...                            % Database name
                         "silas", ...                            % Username
                         "silas", ...                            % Password (or getSecret)
                         "Vendor", "MySQL", ...
@@ -60,7 +63,9 @@ classdef DBHandler < handle
                 obj.refresh();
 
             else
+
                 error('DB seems to be corrupt')
+
             end
 
         end
@@ -68,9 +73,11 @@ classdef DBHandler < handle
 
         function obj = refresh(obj)
             % Get table names and the first rows of each table to understand the structure
+            warning off
             obj.getTableNames();
             obj.getTables();
-            obj.getDistinctValues();
+            warning on
+            % obj.getDistinctValues();
         end
 
         function obj = getTableNames(obj)
@@ -83,7 +90,7 @@ classdef DBHandler < handle
                     result = fetch(obj.conn, 'SELECT name FROM sqlite_master WHERE type="table"');
                     obj.tableNames = result.name;
                 end
-                
+
             catch e
                 error('Failed to retrieve table names: %s', e.message);
             end
@@ -169,14 +176,6 @@ classdef DBHandler < handle
         function healthyDB = dbIsHealthy(obj)
             healthyDB = false;
 
-            num_runs = obj.fetch('SELECT COUNT(*) AS total_runs FROM Runs');
-
-            num_configs = obj.fetch('SELECT COUNT(*) AS total_configurations FROM Configurations');
-
-            num_meas = obj.fetch('SELECT COUNT(*) AS total_measurements FROM Measurements');
-
-            assert((num_runs{1,1}==num_configs{1,1})&&(num_configs{1,1}==num_meas{1,1}),'Different num of entries per table');
-
             %Check for any duplicate paths
             duplictae_raw = obj.fetch("SELECT COALESCE(Runs.rx_raw_path,'NaN') AS rx_raw_path,  COUNT(*) AS occurrences FROM Runs GROUP BY rx_raw_path HAVING COUNT(*) > 1");
             duplictae_sync = obj.fetch("SELECT COALESCE(Runs.rx_sync_path,'NaN') AS rx_sync_path, COUNT(*) AS occurrences FROM Runs GROUP BY rx_sync_path HAVING COUNT(*) > 1");
@@ -186,11 +185,11 @@ classdef DBHandler < handle
                     fprintf('Raw Rx Paths: Found %d duplictaes of %s \n',duplictae_raw.occurrences(i),duplictae_raw.rx_raw_path(i));
                 end
             end
+
             healthyDB = true;
         end
 
 
-
         function lastID = appendToTable(obj, tableName, newRow)
             % appendToTable Appends a new row to the specified table
             %
@@ -206,40 +205,50 @@ classdef DBHandler < handle
                 error('Table %s does not exist in the database or has not been fetched.', tableName);
             end
 
-            % Convert newRow to a table if it is a struct
+            % Handle struct preprocessing before table conversion
             if isstruct(newRow)
                 fields = fieldnames(newRow);
-                emptyFields = structfun(@isempty,newRow);
-                if sum(emptyFields) > 0
-                    emptyFieldNames = fields(emptyFields);  % use () to get a cell array
-                    for idx = 1:numel(emptyFieldNames)
-                        newRow.(emptyFieldNames{idx}) = NaN;
+
+                % % Handle empty fields
+                % emptyFields = structfun(@isempty, newRow);
+                % if sum(emptyFields) > 0
+                %     emptyFieldNames = fields(emptyFields);
+                %     for idx = 1:numel(emptyFieldNames)
+                %         newRow.(emptyFieldNames{idx}) = NaN;
+                %     end
+                % end
+
+                % Convert any non-scalar numeric (including [] and vectors) into JSON
+                for i = 1:numel(fields)
+                    name  = fields{i};
+                    value = newRow.(name);
+
+                    if isnumeric(value) && ~isscalar(value)
+                        % jsonencode([]) -> "[]"
+                        % jsonencode([a,b,c]) -> "[a,b,c]"
+                        newRow.(name) = jsonencode(value);
                     end
                 end
+
+                % Convert to table
                 newRow = struct2table(newRow);
             end
 
             % Ensure the new row matches the structure of the existing table
             existingTableStructure = obj.tables.(tableName);
 
-            % Perform data type checks and conversions
+            % Perform remaining data type checks and conversions
             for colName = newRow.Properties.VariableNames
-                % Extract the value and its intended column type
                 value = newRow.(colName{1});
-                existingValue = existingTableStructure.(colName{1});
 
-                % If the value is a class object, convert it to JSON format
-                if isobject(value) && ~isdatetime(value) && ~isa(value,"string")
+                if iscell(value) && ~isempty(value)
+                    if ischar(value{1}) || isstring(value{1})
+                        newRow.(colName{1}) = value{1};
+                    end
+                elseif isobject(value) && ~isdatetime(value) && ~isa(value, "string")
                     newRow.(colName{1}) = string(jsonencode(value));
-
-                    % If the value is a character array, convert it to a string
                 elseif ischar(value)
                     newRow.(colName{1}) = string(value);
-
-                elseif isdatetime(value)
-
-                    %  newRow.(colName{1}) = string(value);
-
                 end
 
             end
@@ -247,7 +256,6 @@ classdef DBHandler < handle
             % Parameters for retry logic
             maxRetries = 50;
             basePause = 0.05; % seconds
-
             attempt = 0;
             success = false;
 
@@ -256,37 +264,41 @@ classdef DBHandler < handle
                     if obj.type == "mysql"
                         newRow_ = obj.convertTableToCellStrings(newRow);
                     end
-                    sqlwrite(obj.conn, tableName, newRow);
-                    success = true; % Write successful
+                    sqlwrite(obj.conn, tableName, newRow,"Catalog",obj.dataBase);
+                    success = true;
                 catch e
                     if contains(e.message, 'database is locked') || contains(e.message, 'cannot rollback transaction')
                         attempt = attempt + 1;
-                        pauseTime = basePause * (1 + rand()); % Add random jitter to reduce collision chance
-                        fprintf('Database locked. Retry %d/%d after %.3f seconds...\n', attempt, maxRetries, pauseTime);
+                        pauseTime = basePause * (1 + rand());
+                        fprintf('Database locked. Retry %d/%d after %.3f seconds...\n', ...
+                            attempt, maxRetries, pauseTime);
                         pause(pauseTime);
                     else
+                        fprintf('Error details:\n');
+                        fprintf('Column values:\n');
+                        disp(newRow);
                         error('Failed to append to the table %s: %s', tableName, e.message);
                     end
                 end
             end
 
             if ~success
-                error('Failed to append to table %s after %d retries due to database lock.', tableName, maxRetries);
+                error('Failed to append to table %s after %d retries due to database lock.', ...
+                    tableName, maxRetries);
             end
 
-            % Retrieve the measurement_id of the newly inserted row for linking other tables
-
+            % Retrieve the ID of the newly inserted row
             if obj.type == "mysql"
-                queue = "SELECT LAST_INSERT_ID()";
+                query = "SELECT LAST_INSERT_ID()";
             elseif obj.type == "sqlite"
-                queue = "SELECT last_insert_rowid()";
+                query = "SELECT last_insert_rowid()";
             end
 
-            result = fetch(obj.conn, queue);
-            lastID = result{1, 1};  % Access the value directly from the table
+            result = fetch(obj.conn, query);
+            lastID = result{1, 1};
         end
 
-        function exists = checkIfRunExists(obj, table2check, column2check, value2check)
+        function [exists, count] = checkIfRunExists(obj, table2check, column2check, value2check)
             % checkIfRunExists Checks if a specific value exists in a specified column of a table
             %
             % Usage:
@@ -311,7 +323,11 @@ classdef DBHandler < handle
             end
 
             % Construct the query to check for the value in the specified column
-            query = sprintf('SELECT COUNT(*) FROM %s WHERE %s = "%s"', table2check, column2check, value2check);
+            if isnumeric(value2check)
+                query = sprintf('SELECT COUNT(*) FROM %s WHERE %s = %d', table2check, column2check, value2check);
+            else
+                query = sprintf('SELECT COUNT(*) FROM %s WHERE %s = "%s"', table2check, column2check, value2check);
+            end
 
             % Execute the query and pass the value2check to avoid SQL injection issues
             try
@@ -325,12 +341,21 @@ classdef DBHandler < handle
             exists = count > 0;
 
             if exists
-                disp(['The value "', value2check, '" already exists in the column "', column2check, '" of the table "', table2check, '".']);
+                % disp(['The value "', num2str(value2check), '" already exists in the column "', column2check, '" of the table "', table2check, '".']);
             else
                 % disp(['The value "', value2check, '" does not exist in the column "', column2check, '" of the table "', table2check, '".']);
             end
         end
 
+        function hashStr = calcHash(~,object)
+            jsonStr = jsonencode(object);
+            md = java.security.MessageDigest.getInstance('MD5');
+            md.update(uint8(jsonStr));
+            hashBytes = typecast(md.digest, 'uint8');
+            hashStr = lower(dec2hex(hashBytes)');
+            hashStr = lower(strtrim(hashStr(:)')); % Convert to a lowercase string
+        end
+
 
         function resultID = addProcessingResult(obj, run_id, resultData, eqParamsData)
             % addProcessingResult Adds a processing result and links it to an EqualizerParameters entry.
@@ -348,18 +373,19 @@ classdef DBHandler < handle
             %   resultID: The result_id of the newly inserted ProcessingResults entry.
 
             % 1. Compute hash for equalizer parameters
-            jsonStr = jsonencode(eqParamsData);
-            md = java.security.MessageDigest.getInstance('MD5');
-            md.update(uint8(jsonStr));
-            hashBytes = typecast(md.digest, 'uint8');
-            hashStr = lower(dec2hex(hashBytes)');
-            hashStr = lower(strtrim(hashStr(:)')); % Convert to a lowercase string
+            % jsonStr = jsonencode(eqParamsData);
+            % md = java.security.MessageDigest.getInstance('MD5');
+            % md.update(uint8(jsonStr));
+            % hashBytes = typecast(md.digest, 'uint8');
+            % hashStr = lower(dec2hex(hashBytes)');
+            % hashStr = lower(strtrim(hashStr(:)')); % Convert to a lowercase string
+            hashStr = obj.calcHash(eqParamsData);
 
             % Add hash to equalizer parameters
-            eqParamsData.config_hash = hashStr;
+            eqParamsData.hash = hashStr;
 
             % 2. Check if an equalizer configuration with the same hash exists
-            queryStr = sprintf('SELECT eq_id FROM EqualizerParameters WHERE config_hash = ''%s''', eqParamsData.config_hash);
+            queryStr = sprintf('SELECT eq_id FROM Equalizer WHERE hash = ''%s''', eqParamsData.hash);
             existingEntry = obj.fetch(queryStr);
 
             if ~isempty(existingEntry)
@@ -367,27 +393,30 @@ classdef DBHandler < handle
                 eq_id = existingEntry{1,1};
             else
                 % Insert the new equalizer configuration and get its eq_id
-                eq_id = obj.appendToTable('EqualizerParameters', eqParamsData);
+                eqParamsData = eqParamsData.toStruct;
+                eq_id = obj.appendToTable('Equalizer', eqParamsData);
             end
 
             % 3. Add the equalizer configuration reference and run_id to resultData
-            resultData.eqParam_id = eq_id;
+            resultData = resultData.toStruct;
+            resultData.eq_id = eq_id;
             resultData.run_id = run_id;
 
             % 4. Compute hash for the processing result
             tempResultData = rmfield(resultData, 'date_of_processing');
-            resultJsonStr = jsonencode(tempResultData);
-            md2 = java.security.MessageDigest.getInstance('MD5');  % Create a new MD5 instance
-            md2.update(uint8(resultJsonStr));
-            resultHashBytes = typecast(md2.digest, 'uint8');
-            resultHashStr = lower(dec2hex(resultHashBytes)');
-            resultHashStr = lower(strtrim(resultHashStr(:)')); % Convert to a lowercase string
+            % resultJsonStr = jsonencode(tempResultData);
+            % md2 = java.security.MessageDigest.getInstance('MD5');  % Create a new MD5 instance
+            % md2.update(uint8(resultJsonStr));
+            % resultHashBytes = typecast(md2.digest, 'uint8');
+            % resultHashStr = lower(dec2hex(resultHashBytes)');
+            % resultHashStr = lower(strtrim(resultHashStr(:)')); % Convert to a lowercase string
+            resultHashStr = obj.calcHash(tempResultData);
 
             % Add the result hash to resultData
-            resultData.result_hash = resultHashStr;
+            resultData.hash = resultHashStr;
 
             % 5. Check if an identical processing result already exists
-            queryStr2 = sprintf('SELECT result_id FROM Results WHERE result_hash = ''%s''', resultData.result_hash);
+            queryStr2 = sprintf('SELECT result_id FROM Results WHERE hash = ''%s''', resultData.hash);
             existingResult = obj.fetch(queryStr2);
 
             if ~isempty(existingResult)
@@ -397,7 +426,66 @@ classdef DBHandler < handle
                 return;
             end
 
-            % 6. Insert the processing result
+            % 6. check if obj.tables.Results matches Metricstruct
+            % Fields to exclude from comparison
+            excludeFields = {'result_id', 'run_id', 'eq_id'};
+
+            % 7. Chack for new fields in Metric struct and append to
+            % database if necessary
+            ms=Metricstruct;
+            metricFields = setdiff(fieldnames(ms), excludeFields);
+            tableFields = setdiff(fieldnames(obj.tables.Results), excludeFields);
+
+            % Check matches and find missing fields
+            matches = all(ismember(metricFields, tableFields));
+
+            if ~matches
+                missingFields = setdiff(metricFields, tableFields);
+
+
+                % If there are missing fields, add them to the SQL table
+                if ~isempty(missingFields)
+                    for i = 1:length(missingFields)
+                        fieldName = missingFields{i};
+
+                        % Determine SQL data type based on MATLAB class
+                        fieldValue = ms.(fieldName);
+                        if isnumeric(fieldValue)
+                            if isinteger(fieldValue)
+                                sqlType = 'INTEGER';
+                            else
+                                sqlType = 'REAL';
+                            end
+                        elseif ischar(fieldValue) || isstring(fieldValue)
+                            sqlType = 'TEXT';
+                        elseif isdatetime(fieldValue)
+                            sqlType = 'DATETIME';
+                        elseif iscell(fieldValue) || isstruct(fieldValue) || islogical(fieldValue)
+                            sqlType = 'TEXT'; % Store as JSON
+                        else
+                            sqlType = 'TEXT'; % Default to TEXT for unknown types
+                        end
+
+                        % Create ALTER TABLE query
+                        queryStr = sprintf('ALTER TABLE Results ADD COLUMN %s %s', fieldName, sqlType);
+
+                        try
+                            % Execute the query
+                            obj.fetch(queryStr);
+                            fprintf('Added field "%s" of type %s to Results table\n', fieldName, sqlType);
+                        catch ME
+                            fprintf('Error adding field "%s": %s\n', fieldName, ME.message);
+                        end
+                    end
+                else
+                    fprintf('No missing fields to add.\n');
+                end
+
+            end
+
+
+
+            % 8. Insert the processing result
             resultID = obj.appendToTable('Results', resultData);
         end
 
@@ -486,17 +574,32 @@ classdef DBHandler < handle
         end
 
 
-
-
         function executeSQL(obj, query)
             % This method executes an SQL statement using MATLAB's execute function.
             execute(obj.conn, query);
         end
 
-        function answer = fetch(obj,query)
-            answer = fetch(obj.conn,query);
+
+        function answer = fetch(obj, query)
+            maxFast = 20; maxSlow = 30;
+            for attempt = 1:maxSlow
+                try
+                    answer = fetch(obj.conn, query);
+                    return
+                catch ME
+                    if attempt < maxFast
+                        pause(0.1)
+                    else
+                        pause(1)
+                    end
+                    lastErr = ME;
+                end
+            end
+            error('Database fetch failed after %d attempts:\n%s', maxSlow, lastErr.getReport())
         end
 
+
+
         function [result,query] = queryDB(obj, filterParams, selectedFields)
             % getPathsWithFlexibleFilter Retrieves values from Runs table with flexible filtering
             % and lets the user select which fields to include in the SELECT statement.
@@ -518,22 +621,6 @@ classdef DBHandler < handle
                 selectedFields = [];
             end
 
-            % Step 1: Prompt the user to input filter parameters if not provided
-            if isempty(filterParams)
-                filterParams = obj.promptFilterParameters();
-            end
-
-            % Step 2: Prompt the user to select fields to include in the SELECT statement
-            if isempty(selectedFields)
-                selectedFields = obj.promptSelectFields();
-            else
-                if iscell(selectedFields)
-
-                elseif isstruct(selectedFields)
-
-                end
-            end
-
             % Step 3: Construct the SQL query based on the inputs
             query = obj.constructSQLQuery(filterParams, selectedFields);
 
@@ -558,7 +645,7 @@ classdef DBHandler < handle
                         cleanedTable.(varNames{i}) = numCol;
                     else
                         % Clean double-quoted SQL literals (e.g., ""no_db"")
-                        cleanedTable.(varNames{i}) = strrep(string(col), '""', '"');
+                        cleanedTable.(varNames{i}) = strrep(string(col), '"', '');
                     end
                 end
             end
@@ -619,132 +706,238 @@ classdef DBHandler < handle
             end
         end
 
-
         function query = constructSQLQuery(obj, filterParams, selectedFields)
             % constructSQLQuery Constructs the SQL query based on filter parameters and selected fields.
-            %
-            % If selectedFields is provided as a struct, it is converted to a cell array.
-            % The conversion takes the field names and creates entries like:
-            %    {'selectedFields.fieldName'} for each field.
 
-            % Input check for selectedFields: if it's a struct, convert it to a cell array.
-            if isstruct(selectedFields)
+            arguments
+                obj
+                filterParams
+                selectedFields
+            end
+
+            % -------- Step 1: Normalize selectedFields to {'Table.field', ...} --------
+            if isempty(selectedFields) || (ischar(selectedFields) && strcmpi(selectedFields, 'all'))
+                selectedFields = obj.getTableFieldNames('Runs'); % default
+            elseif isstruct(selectedFields)
                 newFields = {};
                 tableNames = fieldnames(selectedFields);
                 for t = 1:numel(tableNames)
                     tableStruct = selectedFields.(tableNames{t});
-                    fieldNames = fieldnames(tableStruct);
-                    for f = 1:numel(fieldNames)
-                        if isequal(tableStruct.(fieldNames{f}), 1)
-                            newFields{end+1} = sprintf('%s.%s', tableNames{t}, fieldNames{f});
+                    fns = fieldnames(tableStruct);
+                    for f = 1:numel(fns)
+                        if isequal(tableStruct.(fns{f}), 1)
+                            newFields{end+1} = sprintf('%s.%s', tableNames{t}, fns{f}); %#ok<AGROW>
                         end
                     end
                 end
                 selectedFields = newFields;
             end
 
-            % Construct the SELECT clause dynamically based on user selection.
-            % (Assuming that when provided as a cell array, each entry is of the form
-            % 'TableName.fieldName' or, in our conversion case, 'selectedFields.fieldName'.)
-            selectClause = 'SELECT DISTINCT ';
-            for i = 1:numel(selectedFields)
-                fieldParts = strsplit(selectedFields{i}, '.');
-                % If the field comes from the struct conversion, its first part is 'selectedFields'
-                % and the actual field name is in the second part.
-                if strcmp(fieldParts{1}, 'selectedFields')
-                    tableName = fieldParts{1};  % not used for type checking below
-                    fieldName = fieldParts{2};
-                else
-                    tableName = fieldParts{1};
-                    fieldName = fieldParts{2};
-                end
+            % Parse the table names actually referenced by the SELECT
+            reqTables = unique(cellfun(@(s) extractBefore(s, '.'), selectedFields, ...
+                'UniformOutput', false));
 
-                % Decide on COALESCE depending on the field type.
-                % If the table is known in obj.tables and the field is numeric, use 'NaN'.
-                if isfield(obj.tables, tableName) && isfield(obj.tables.(tableName), fieldName) && isnumeric(obj.tables.(tableName).(fieldName))
-                    selectClause = [selectClause, 'COALESCE(', selectedFields{i}, ', ''NaN'') AS ', fieldName];
-                else
-                    selectClause = [selectClause, 'COALESCE(', selectedFields{i}, ', '''') AS ', fieldName];
-                end
+            % -------- Step 2: Build SELECT with COALESCE wrapper as you already do ----
+            selectClause = obj.generateCoalesceString(selectedFields);
 
-                if i < numel(selectedFields)
-                    selectClause = [selectClause, ', '];
-                else
-                    selectClause = [selectClause, ' '];
+            % -------- Step 3: FROM and minimal JOIN plan ------------------------------
+            % Decide main table: prefer the first explicitly referenced table, else 'Runs'
+            if ~isempty(reqTables)
+                mainTable = reqTables{1};
+            else
+                mainTable = 'Runs';
+            end
+
+            % If WHERE references a table not in reqTables (e.g., Runs.*), make sure it’s present.
+            whereClause = '';
+            if ~isempty(filterParams)
+                whereClause = obj.generateWhereClause(filterParams);
+                % Heuristic: add 'Runs' if WHERE clause mentions 'Runs.'
+                if contains(whereClause, 'Runs.')
+                    reqTables = unique([reqTables; {'Runs'}]); %#ok<AGROW>
                 end
             end
 
+            % Ensure main table is included
+            if ~ismember(mainTable, reqTables)
+                reqTables = unique([mainTable; reqTables]); %#ok<AGROW>
+            end
 
-            % --- Adaptive FROM Clause ---
-            mainTable = 'Runs';
-            fromClause = ['FROM ', mainTable, ' '];
+            % We’ll build joins only for the required tables (minus the main)
+            otherTables = setdiff(reqTables, {mainTable});
 
-            % Prepare join buffers
-            normalJoins = '';
-            equalizerJoin = '';
+            % Keep track of what’s already in the FROM graph (start with main)
+            present = string(mainTable);
+            joins = strings(0,1);
 
-            tableNamesAll = fieldnames(obj.tables);
-            for t = 1:numel(tableNamesAll)
-                tableName = tableNamesAll{t};
+            % Helper lambdas
+            hasField  = @(tbl, fld) isfield(obj.tables.(char(tbl)), char(fld));
+            canJoinBy = @(left, right, key) hasField(left, key) && hasField(right, key);
 
-                if strcmpi(tableName, mainTable) || strcmpi(tableName, 'sqlite_sequence')
-                    continue;
+
+            % A small helper that adds a LEFT JOIN if the right table isn't present yet
+            function addJoinByKey(rightTbl, key)
+                if any(present == string(rightTbl))
+                    return; % already joined
                 end
+                % Prefer to join against an already-present table that has the key
+                anchor = '';
+                for k = 1:numel(present)
+                    if canJoinBy(char(present(k)), rightTbl, key)
+                        anchor = char(present(k));
+                        break;
+                    end
+                end
+                if isempty(anchor)
+                    % No anchor in current graph; if the right table is 'Equalizer' and key is eq_id,
+                    % try to ensure a bridge table with eq_id exists (Results or a dashboard view).
+                    if strcmpi(rightTbl,'Equalizer') && strcmpi(key,'eq_id')
+                        % Bring in one eq_id-capable table if it is requested
+                        bridgeOrder = {'Results','dashboard_old','dashboard_new','dashboard_ungrouped'};
+                        for b = 1:numel(bridgeOrder)
+                            br = bridgeOrder{b};
+                            if ismember(br, reqTables) && ~any(present == string(br)) && hasField(obj.tables.(br),'eq_id')
+                                % Attach bridge by run_id if possible, otherwise leave for eq_id
+                                if any(present == "Runs") && hasField(obj.tables.(br),'run_id') && hasField(obj.tables.('Runs'),'run_id')
+                                    joins(end+1,1) = "LEFT JOIN " + br + " ON Runs.run_id = " + br + ".run_id";
+                                    present(end+1,1) = string(br);
+                                    anchor = br; % we can now anchor Equalizer on eq_id to this
+                                    break;
+                                else
+                                    % Fallback: anchor to main if it shares eq_id
+                                    for k = 1:numel(present)
+                                        pk = char(present(k));
+                                        if canJoinBy(pk, br, 'eq_id')
+                                            joins(end+1,1) = "LEFT JOIN " + br + " ON " + pk + ".eq_id = " + br + ".eq_id";
+                                            present(end+1,1) = string(br);
+                                            anchor = br;
+                                            break;
+                                        end
+                                    end
+                                    if ~isempty(anchor), break; end
+                                end
+                            end
+                        end
+                    end
+                end
+                % Re-scan for an anchor (maybe the bridge helped)
+                if isempty(anchor)
+                    for k = 1:numel(present)
+                        if canJoinBy(char(present(k)), rightTbl, key)
+                            anchor = char(present(k));
+                            break;
+                        end
+                    end
+                end
+                if isempty(anchor)
+                    % As a final fallback, if the main table is Runs and right has run_id, join by run_id
+                    if strcmpi(mainTable,'Runs') && hasField(obj.tables.(rightTbl),'run_id') && hasField(obj.tables.('Runs'),'run_id')
+                        anchor = 'Runs';
+                        key = 'run_id';
+                    end
+                end
+                if isempty(anchor)
+                    % Could not find a path; skip join silently (or throw if you prefer strict)
+                    return;
+                end
+                joins(end+1,1) = "LEFT JOIN " + rightTbl + " ON " + anchor + "." + key + " = " + rightTbl + "." + key;
+                present(end+1,1) = string(rightTbl);
+            end
 
-                if isfield(obj.tables.(tableName), 'run_id')
-                    % Direct join to Runs
-                    normalJoins = [normalJoins, 'LEFT JOIN ', tableName, ' ON ', mainTable, '.run_id = ', tableName, '.run_id '];
-                elseif isfield(obj.tables.(tableName), 'eq_id')
-                    % Equalizer depends on Results, collect this join separately
-                    equalizerJoin = ['LEFT JOIN ', tableName, ' ON Results.eqParam_id = ', tableName, '.eq_id '];
+            % First pass: if WHERE uses Runs.* and mainTable isn’t Runs, ensure Runs is in the graph
+            if contains(string(whereClause), "Runs.") && ~any(present == "Runs")
+                % Try to join Runs to whatever has run_id (mainTable ideally)
+                if hasField(mainTable, 'run_id') && hasField('Runs', 'run_id')
+                    joins(end+1,1) = "LEFT JOIN Runs ON " + string(mainTable) + ".run_id = Runs.run_id";
+                    present(end+1,1) = "Runs";
                 end
             end
 
-            % Combine joins: normal joins first, Equalizer last
-            fromClause = [fromClause, normalJoins, equalizerJoin];
+            % Join the required tables with minimal edges
+            for i = 1:numel(otherTables)
+                tbl = otherTables{i};
+                % Prefer run_id join if possible, else eq_id, else skip
+                if any(present == "Runs") && hasField(tbl,'run_id')
+                    addJoinByKey(tbl, 'run_id');
+                elseif hasField(tbl,'run_id') && hasField(mainTable,'run_id')
+                    addJoinByKey(tbl, 'run_id');
+                elseif hasField(tbl,'eq_id')
+                    addJoinByKey(tbl, 'eq_id');
+                else
+                    % no obvious key; skip
+                end
+            end
 
+            % Build the FROM clause
+            fromClause = "FROM " + string(mainTable) + " " + strjoin(joins, " ");
 
+            % -------- Step 4: WHERE (unchanged logic) ------------------------------
+            if ~isempty(whereClause)
+                query = char(strjoin([selectClause, fromClause, "WHERE " + string(whereClause)], " "));
+            else
+                query = char(strjoin([selectClause, fromClause], " "));
+            end
+        end
 
-
-            % --- WHERE Clause Construction ---
-            baseQuery = [selectClause, ' ', fromClause, 'WHERE '];
+        function whereClause = generateWhereClause(obj, filterParams)
             filterClauses = [];
-            tableNames_ = fieldnames(filterParams);
-            for t = 1:numel(tableNames_)
-                tableName = tableNames_{t};
+
+            % If filterParams is a DbFilterParameter object, get its internal structure
+            if isa(filterParams, 'QueryFilter')
+                filterParams = filterParams.toStruct();
+            end
+
+            % Now proceed with the structure
+            tableNames = fieldnames(filterParams);
+
+            for t = 1:numel(tableNames)
+                tableName = tableNames{t};
                 tableParams = filterParams.(tableName);
                 fieldNames = fieldnames(tableParams);
+
                 for i = 1:numel(fieldNames)
                     fieldName = fieldNames{i};
                     value = tableParams.(fieldName);
                     fullName = sprintf('%s.%s', tableName, fieldName);
 
-                    % Handle various types of values for SQL query construction
-                    if isempty(value)
+                    % Skip empty values
+                    if isempty(value) || (isa(value, 'QueryFilter') && isempty(value.value))
                         continue;
-                    elseif isnumeric(value) && isnan(value)
-                        filterClause = sprintf('%s IS NULL', fullName);
-                    elseif isnumeric(value) && ~isEnumeration(value)
-                        filterClause = sprintf('%s = %f', fullName, value);
-                    elseif islogical(value) || (isnumeric(value) && ismember(value, [0, 1])) && ~isEnumeration(value)
-                        filterClause = sprintf('%s = %d', fullName, value);
-                    elseif ischar(value) || isstring(value)
-                        filterClause = sprintf('%s = ''%s''', fullName, char(value));
-                    elseif isEnumeration(value)
-                        filterClause = sprintf('%s = ''%s''', fullName, value);
-                    else
-                        error('Unsupported data type for field "%s".', fullName);
                     end
-                    filterClauses = [filterClauses, filterClause, ' AND '];
+
+                    % Handle Filter class
+                    if isa(value, 'SqlFilter')
+                        if isnumeric(value.value)
+                            filterClause = sprintf('%s %s %f', ...
+                                fullName, value.operator, value.value);
+                        elseif ischar(value.value) || isstring(value.value)
+                            filterClause = sprintf('%s %s ''%s''', ...
+                                fullName, value.operator, char(value.value));
+                        else
+                            continue; % Skip unsupported types
+                        end
+                        filterClauses = [filterClauses, filterClause, ' AND '];
+                    else
+                        % Handle direct values (legacy support)
+                        if isnumeric(value) && isnan(value)
+                            filterClause = sprintf('%s IS NULL', fullName);
+                        elseif isnumeric(value)
+                            filterClause = sprintf('%s = %f', fullName, value);
+                        elseif ischar(value) || isstring(value)
+                            filterClause = sprintf('%s = ''%s''', fullName, char(value));
+                        else
+                            continue; % Skip unsupported types
+                        end
+                        filterClauses = [filterClauses, filterClause, ' AND '];
+                    end
                 end
             end
 
-            % Remove trailing ' AND ' if any filters were added.
+            % Remove trailing ' AND ' if any filters were added
             if ~isempty(filterClauses)
-                filterClauses = filterClauses(1:end-5);
-                query = [selectClause, ' ', fromClause, 'WHERE ', filterClauses];
+                whereClause = filterClauses(1:end-5);
             else
-                query = [selectClause, ' ', fromClause];
+                whereClause = '';
             end
         end
 
@@ -980,5 +1173,59 @@ classdef DBHandler < handle
         end
 
 
+        function fieldNames = getTableFieldNames(obj, tableName)
+            % Returns all field names for a given table as a cell array in the format {'TableName.fieldName'}
+            if isfield(obj.tables, tableName)
+                % Get raw field names
+                rawFields = fieldnames(obj.tables.(tableName));
+
+                % Create cell array with table name prefix
+                fieldNames = cellfun(@(x) [tableName, '.', x], ...
+                    rawFields, ...
+                    'UniformOutput', false);
+            else
+                error('Table "%s" not found in obj.tables.', tableName);
+            end
+        end
+
+
+        function coalesceStr = generateCoalesceString(obj, selectedFields)
+            % Generates a COALESCE string for selected fields
+            % Input:
+            %   selectedFields: cell array of strings in format {'Table.field'}
+            %                  e.g., {'Runs.run_id', 'Runs.bitrate'}
+            % Output:
+            %   coalesceStr: string with COALESCE statements
+
+            arguments
+                obj
+                selectedFields cell
+            end
+
+            % Initialize cell array to store each COALESCE statement
+            coalesceStatements = cell(length(selectedFields), 1);
+
+            % Generate COALESCE statement for each field
+            for i = 1:length(selectedFields)
+                % Split table and field name
+                parts = strsplit(selectedFields{i}, '.');
+                if length(parts) ~= 2
+                    error('Field name must be in format "Table.field": %s', selectedFields{i});
+                end
+                tableName = parts{1};
+                fieldName = parts{2};
+
+                % Generate COALESCE statement
+                coalesceStatements{i} = sprintf('COALESCE(%s.%s, ''NaN'') AS %s ', ...
+                    tableName, fieldName, fieldName);
+            end
+
+            % Join with comma, newline and MATLAB string continuation
+            coalesceStr = strjoin(coalesceStatements, [', ' sprintf('\n    ')]);
+            % Add initial newline and spacing for formatting
+            coalesceStr = [sprintf('SELECT DISTINCT  \n    ') coalesceStr];
+        end
+
+
     end
 end

Classes/DataBaseHandler/copyStylingFrom.m

@@ -1,76 +0,0 @@
-function copyStylingFrom(figNumSource, figNumTgt)
-    % Get handles to the source and target figures
-    sourceFig = figure(figNumSource);
-    targetFig = figure(figNumTgt);
-    
-    % Get axes of source and target figures
-    sourceAxes = findall(sourceFig, 'type', 'axes');
-    targetAxes = findall(targetFig, 'type', 'axes');
-    
-    % Ensure the number of axes match
-    if length(sourceAxes) ~= length(targetAxes)
-        error('Number of axes in source and target figures must be the same.');
-    end
-
-    % Loop through each pair of axes and copy styling properties
-    for i = 1:length(sourceAxes)
-        copyAxesProperties(sourceAxes(i), targetAxes(i));
-    end
-    
-    % Apply general figure properties if desired
-    targetFig.Color = sourceFig.Color; % Background color
-end
-
-function copyAxesProperties(sourceAx, targetAx)
-    % List of properties to copy from source to target axes
-    propsToCopy = {'XColor', 'YColor', 'ZColor', 'FontSize', 'FontName', ...
-                   'GridColor', 'GridLineStyle', 'MinorGridColor', 'Box', ...
-                   'XGrid', 'YGrid', 'ZGrid', 'XMinorGrid', 'YMinorGrid', 'ZMinorGrid', ...
-                   'LineWidth', 'TitleFontSizeMultiplier', 'LabelFontSizeMultiplier'};
-    
-    % Copy properties from source to target
-    for i = 1:length(propsToCopy)
-        try
-            targetAx.(propsToCopy{i}) = sourceAx.(propsToCopy{i});
-        catch
-            % Skip property if it doesn't exist or can't be copied
-        end
-    end
-    
-    % Copy axis labels and titles
-    targetAx.Title.String = sourceAx.Title.String;
-    targetAx.XLabel.String = sourceAx.XLabel.String;
-    targetAx.YLabel.String = sourceAx.YLabel.String;
-    targetAx.ZLabel.String = sourceAx.ZLabel.String;
-    
-    % Copy children elements like lines, patches, etc.
-    sourceChildren = allchild(sourceAx);
-    targetChildren = allchild(targetAx);
-    
-    % Ensure the number of children elements match
-    if length(sourceChildren) ~= length(targetChildren)
-        warning('Number of elements in source and target axes differ. Styling may not be applied completely.');
-    end
-
-    % Copy properties of children (like lines, patches, etc.), except colors and legends
-    for i = 1:min(length(sourceChildren), length(targetChildren))
-        copyObjectProperties(sourceChildren(i), targetChildren(i));
-    end
-end
-
-function copyObjectProperties(sourceObj, targetObj)
-    % List of common properties to copy for plot elements (lines, patches, etc.)
-    propsToCopy = {'LineStyle', 'LineWidth', 'Marker', 'MarkerSize', ...
-                   'MarkerEdgeColor', 'MarkerFaceColor', 'DisplayName'};
-    
-    % Copy properties from source to target, excluding colors
-    for i = 1:length(propsToCopy)
-        try
-            if ~contains(propsToCopy{i}, 'Color') % Skip color properties
-                targetObj.(propsToCopy{i}) = sourceObj.(propsToCopy{i});
-            end
-        catch
-            % Skip property if it doesn't exist or can't be copied
-        end
-    end
-end

Classes/Warehouse_class/classes/DataStorage.m

@@ -51,8 +51,9 @@ classdef DataStorage < handle
         function save(obj,path)
             try
                 save(path,"obj");
-            catch
-
+            catch e
+                disp(e.message)
+                disp('Provide save path')
             end
         end
 
@@ -136,7 +137,21 @@ classdef DataStorage < handle
                     if ~isempty(tmp)
 
                         if isa(tmp,'double')
-                            value(i) = tmp ;
+                            try
+                                value(i,:) = tmp ;
+                            catch
+                                a = size(value,2);
+                                b = size(tmp,2);
+                                if a > b
+                                    value(i,:) =[tmp,NaN(1,a-b)] ;
+                                elseif a < b
+                                    value(i,:) = tmp(1:size(value,2)) ;
+                                else
+                                    error('unknwon case...')
+                                end
+
+                                
+                            end
                         elseif isa(tmp,'Signal') || isa(tmp,'struct') ||  isa(tmp,'Exfo_laser') ||  isa(tmp,'DC_supply')
                             if i == 1
                                 value = {};
@@ -150,7 +165,7 @@ classdef DataStorage < handle
                                 end
                                 value{i} = tmp{1} ;
                             else
-                                value{i} = tmp{1} ;
+                                value{i} = tmp ;
                             end
 
                         else

Classes/Warehouse_class/functions/fwm_plots/CompleteRoutine.m

@@ -1,251 +0,0 @@
-% Script, that shows the data management routine :-)
-
-loadExistingWareHouse = 0;
-
-if loadExistingWareHouse
-
-    [file, path] = uigetfile();
-    wh = load([path filesep file]);
-    wh = wh.wh;
-    wh.showInfo;
-
-else
-
-    % 1) Define all your parameters, best practice directly constructs a
-    % structure
-
-    params = struct;
-
-    params.l               = [2,10];
-
-    params.dispersion      = [0];
-
-    params.sgm             = [0];
-
-%     params.pol             =  ["YXYXYXYX","YXXYYXXY","YYYYYYYY"];
-    params.pol             =  ["alternated","paired","copolarized"];
-
-    params.p_in            = [3];
-
-    params.p_out           = [-12,-11,-10,-9,-8,-7,-6,-5,-4,-3,-2];
-
-    params.pmd             = [0.1];
-
-    params.gamma           = [0.0023];
-
-    params.realization     = [1:20];
-
-    params.numchannels      = [16];
-
-    params.center_wavelength =  floor([getSweepWavelengths(35, 50e9, 1310)] .* 1000) ./ 1000 ;
-    params.center_wavelength =  [1285 1287 1290 1292 1295];
-    params.center_wavelength = 1310;
-
-    params.channelspacing = [400e9];
-
-    params.random_zdw = [0];
-
-    %wh = warehouse :-)
-    wh = DataStorage(params);
-
-    wh.showInfo;
-
-    wh.addStorage("ber");
-
-    wh.addStorage("totalBer");
-
-end
-
-
-
-%2) Simulate a bunch of data - TO BE IMPLEMENTED HERE - for now use scripts
-%from Sebastian
-
-%3) Once the simulation folder is around, specifiy path and analyze dirs
-
-path = uigetdir('C:\Users\Silas\Documents\MATLAB\Raw_Cluster_Simulations');
-
-allMat = getAllFilesInFolder(path,'.mat');
-allErr = getAllFilesInFolder(path,'.err');
-
-%allMat = dir([path filesep '*.mat']);
-
-%allErr = dir([path filesep '*.err']);
-
-if numel(allMat) == 0
-    warning('You defined an empty folder. Could not locate any .mat file.')
-else
-    fprintf('%-20s', 'Err Files:'); fprintf('%-12s', num2str(numel(allErr))); fprintf('\n');
-    fprintf('%-20s', 'Mat Files:'); fprintf('%-12s', num2str(numel(allMat))); fprintf('\n');
-    fprintf('%-20s', 'Missing Mat Files:'); fprintf('%-12s', num2str(numel(allErr)-numel(allMat))); fprintf('\n');
-end
-
-%4) Now load that data
-
-f = waitbar(0,'Please wait...');
-cnt = 0;
-
-for num = 1:numel(allMat)
-
-     fileName    = allMat(num).name;
-     fileFolder  = allMat(num).path;
-     fileExt = allMat(num).ext;
-% 
-    matFile     = load([fileFolder filesep fileName fileExt]);
-    matFile     = matFile.loop_data;
-
-    % ____________________________________
-    % FIND THE DATAPOINT CURRENTLY LOADED
-    zdw = 1310;
-
-    channelplan = "symmetric";
-
-    channelspacing = str2double(strrep(regexp(fileName,'(_chsp)+([\d]*)','match'),'_chsp','')).*1e9;
-
-    numchannels = str2double(strrep(regexp(fileName,'(ch)+(_)+([\d]*)','match'),'ch_',''));
-
-    center_wavelength = str2double(insertAfter(strrep(regexp(fileName,'(lambda)+([\d]*)','match'),'lambda',''),4,'.'));
-
-    center_wavelength = floor(center_wavelength * 1000) / 1000;
-
-    if center_wavelength == 2192
-        continue
-    end
-    
-    center_wavelength = 1310;
-
-    random_zdw = str2double(strrep(regexp(fileName,'(rzwd)+([\d])','match'),'rzwd',''));
-
-    l =     str2double(strrep(regexp(fileName,'([L])+(_)+([\d]*)','match'),'L_',''));
-
-    d =     str2double(strrep(regexp(fileName,'([D])+(_)+([\d]*)','match'),'D_',''));
-
-    if d == 0
-        sgm = false;
-    else
-        sgm = true;
-    end
-
-
-    if numel(regexp(fileName,'(YYYY)','match')) > 1
-        pol = "copolarized";
-    elseif numel(regexp(fileName,'(YXXY)','match')) > 1
-        pol = "paired";
-    elseif numel(regexp(fileName,'(YXYX)','match')) > 1
-        pol = "alternated";
-    else
-        pol = "copolarized";
-    end
-
-    p_in =  str2double(strrep(regexp(fileName,'(pow_)+([-,\d]{1})','match'),'pow_',''));
-
-    pmd = 0.1;
-
-    gamma = 0.0023;
-
-    realiz = str2double(strrep(regexp(fileName,'(r)+([-,\d]{1,3})','match'),'r',''));
-
-    
-
-    % ____________________________________
-    % Get the information you want from current file
-    rop=[];
-    ber = [];
-    for pow = 2:12
-
-        module_number = '';
-        for p = 1:11 %11 because there are 11 ROP branches in model
-
-            % get ROP
-            if p == 1
-                p_out =  matFile.dp_optatten_para.atten;
-            else
-                p_out =  matFile.("dp_optatten__"+(p)+"_para").atten;
-            end
-
-            p_out = round(p_out-10*log10(numel(matFile.config.parameters.common.wavelengthPlan)));
-
-            for c = 1:numel(matFile.config.parameters.common.wavelengthPlan)
-
-                ber(c) = matFile.("prms_compare_wdm__"+(p+1)+"_out"){1,c}.ber;
-
-            end
-
-            totalBer = matFile.("prms_compare_wdm__"+(p+1)+"_out"){1,end}.totalBer;
-
-            if totalBer > 0.2 && channelspacing == 400e9 &&  pol == "alternated"
-                disp("stopping here");
-                pause;
-            end
-           
-
-            % ____________________________________
-            % Add value to warehouse at the correct position
-
-    
-
-            wh.addValueToStorage(ber,'ber',l,d,sgm,pol,p_in,p_out,pmd,gamma,realiz,numchannels, center_wavelength,channelspacing,random_zdw);
-            wh.getStoValue('ber',l,d,sgm,pol,p_in,p_out,pmd,gamma,realiz,numchannels, center_wavelength,channelspacing,random_zdw);
-            wh.addValueToStorage(totalBer,'totalBer',l,d,sgm,pol,p_in,p_out,pmd,gamma,realiz,numchannels,center_wavelength,channelspacing,random_zdw);
-
-        end
-
-    end
-    waitbar(num/numel(allMat),f,'Loading your data');
-end
-
-close(f)
-
-
-
-
-
-% 4) Hey! the warehouse is here and (hopefully) filled with data  :-)
-
-% Create a save dialog
-defaultDir = 'C:\Users\Silas\Documents\MATLAB\Raw_Cluster_Simulations\';
-defaultExt = '*.mat';
-[filename, pathname] = uiputfile(fullfile(defaultDir, defaultExt),'', 'wh.mat');
-
-% Check if the user pressed Cancel
-if isequal(filename, 0) || isequal(pathname, 0)
-    disp('Save operation canceled.');
-else
-    % Save the variable to the selected file
-    save(fullfile(pathname, filename), 'wh');
-    disp(['Variable "wh" saved to: ', fullfile(pathname, filename)]);
-end
-
-
-
-
-function matFileStructArray = getAllFilesInFolder(folderPath,extension)
-    % Get a list of all files in the current folder
-    currentFolderFiles = dir(fullfile(folderPath, '*'));
-
-    % Exclude '.' and '..' directories
-    currentFolderFiles = currentFolderFiles(~ismember({currentFolderFiles.name}, {'.', '..'}));
-
-    % Initialize the structure array for .mat files
-    matFileStructArray = struct('path', {}, 'name', {}, 'ext', {});
-
-    % Loop over each file in the current folder
-    for i = 1:length(currentFolderFiles)
-        currentFile = currentFolderFiles(i);
-
-        % Check if the current item is a file and has a .mat extension
-        if ~currentFile.isdir && endsWith(currentFile.name, extension, 'IgnoreCase', true)
-            % If it's a .mat file, add it to the structure array
-            [matFileStructArray(end + 1).path,matFileStructArray(end+1).name, matFileStructArray(end+1).ext] = fileparts(fullfile(folderPath, currentFile.name));
-        elseif currentFile.isdir
-            % If it's a directory, recursively call the function
-            subfolderPath = fullfile(folderPath, currentFile.name);
-            subfolderMatFiles = getAllFilesInFolder(subfolderPath,extension);
-
-            % Add .mat files from the subfolder to the structure array
-            matFileStructArray = [matFileStructArray, subfolderMatFiles];
-        end
-    end
-end
-
-

Classes/Warehouse_class/functions/fwm_plots/CompleteRoutine_DifferentChannels_Fig3.m

@@ -1,250 +0,0 @@
-% Script, that shows the data management routine :-)
-
-loadExistingWareHouse = 0;
-
-if loadExistingWareHouse
-
-    [file, path] = uigetfile();
-    wh = load([path filesep file]);
-    wh = wh.wh;
-    wh.showInfo;
-
-else
-
-    % 1) Define all your parameters, best practice directly constructs a
-    % structure
-
-    params = struct;
-
-    params.l               = [2, 10];
-
-    params.dispersion      = [0, 3];
-
-    params.sgm             = [0, 1];
-
-%     params.pol             =  ["YXYXYXYX","YXXYYXXY","YYYYYYYY"];
-    params.pol             =  ["alternated","paired","copolarized"];
-
-    params.p_in            = [3];
-
-    params.p_out           = [-12,-11,-10,-9,-8,-7,-6,-5,-4,-3,-2];
-
-    params.pmd             = [0.1];
-
-    params.gamma           = [0.0023];
-
-    params.realization     = [1:20];
-
-    params.numchannels      = [1,2,4,8,16];
-
-    params.center_wavelength =  floor([getSweepWavelengths(35, 50e9, 1310)] .* 1000) ./ 1000 ;
-    params.center_wavelength =  [1285 1287 1290 1292 1295];
-    params.center_wavelength = 1310;
-
-    params.channelspacing = [400e9];
-
-    params.random_zdw = [0,1];
-
-    %wh = warehouse :-)
-    wh = DataStorage(params);
-
-    wh.showInfo;
-
-    wh.addStorage("ber");
-
-    wh.addStorage("totalBer");
-
-end
-
-
-
-%2) Simulate a bunch of data - TO BE IMPLEMENTED HERE - for now use scripts
-%from Sebastian
-
-%3) Once the simulation folder is around, specifiy path and analyze dirs
-
-path = uigetdir('C:\Users\Silas\Documents\MATLAB\Raw_Cluster_Simulations');
-
-allMat = getAllFilesInFolder(path,'.mat');
-allErr = getAllFilesInFolder(path,'.err');
-
-%allMat = dir([path filesep '*.mat']);
-
-%allErr = dir([path filesep '*.err']);
-
-if numel(allMat) == 0
-    warning('You defined an empty folder. Could not locate any .mat file.')
-else
-    fprintf('%-20s', 'Err Files:'); fprintf('%-12s', num2str(numel(allErr))); fprintf('\n');
-    fprintf('%-20s', 'Mat Files:'); fprintf('%-12s', num2str(numel(allMat))); fprintf('\n');
-    fprintf('%-20s', 'Missing Mat Files:'); fprintf('%-12s', num2str(numel(allErr)-numel(allMat))); fprintf('\n');
-end
-
-%4) Now load that data
-
-f = waitbar(0,'Please wait...');
-cnt = 0;
-
-for num = 1:numel(allMat)
-
-     fileName    = allMat(num).name;
-     fileFolder  = allMat(num).path;
-     fileExt = allMat(num).ext;
-% 
-    matFile     = load([fileFolder filesep fileName fileExt]);
-    matFile     = matFile.loop_data;
-
-    % ____________________________________
-    % FIND THE DATAPOINT CURRENTLY LOADED
-    zdw = 1310;
-
-    channelplan = "symmetric";
-
-    channelspacing = str2double(strrep(regexp(fileName,'(_chsp)+([\d]*)','match'),'_chsp','')).*1e9;
-
-    numchannels = str2double(strrep(regexp(fileName,'(ch)+(_)+([\d]*)','match'),'ch_',''));
-
-    center_wavelength = str2double(insertAfter(strrep(regexp(fileName,'(lambda)+([\d]*)','match'),'lambda',''),4,'.'));
-
-    center_wavelength = floor(center_wavelength * 1000) / 1000;
-
-    if center_wavelength == 2192
-        continue
-    end
-    
-    center_wavelength = 1310;
-
-    random_zdw = str2double(strrep(regexp(fileName,'(rzwd)+([\d])','match'),'rzwd',''));
-
-    l =     str2double(strrep(regexp(fileName,'([L])+(_)+([\d]*)','match'),'L_',''));
-
-    d =     str2double(strrep(regexp(fileName,'([D])+(_)+([\d]*)','match'),'D_',''));
-
-    if d == 0
-        sgm = false;
-    else
-        sgm = true;
-    end
-
-
-    if numel(regexp(fileName,'(YYYY)','match')) > 1
-        pol = "copolarized";
-    elseif numel(regexp(fileName,'(YXXY)','match')) > 1
-        pol = "paired";
-    elseif numel(regexp(fileName,'(YXYX)','match')) > 1
-        pol = "alternated";
-    else
-        pol = "copolarized";
-    end
-
-    p_in =  str2double(strrep(regexp(fileName,'(pow_)+([-,\d]{1})','match'),'pow_',''));
-
-    pmd = 0.1;
-
-    gamma = 0.0023;
-
-    realiz = str2double(strrep(regexp(fileName,'(r)+([-,\d]{1,3})','match'),'r',''));
-
-    
-
-    % ____________________________________
-    % Get the information you want from current file
-    rop=[];
-    ber = [];
-    for pow = 2:12
-
-        module_number = '';
-        for p = 1:11 %11 because there are 11 ROP branches in model
-
-            % get ROP
-            if p == 1
-                p_out =  matFile.dp_optatten_para.atten;
-            else
-                p_out =  matFile.("dp_optatten__"+(p)+"_para").atten;
-            end
-
-            p_out = round(p_out-10*log10(numel(matFile.config.parameters.common.wavelengthPlan)));
-
-            for c = 1:numel(matFile.config.parameters.common.wavelengthPlan)
-
-                ber(c) = matFile.("prms_compare_wdm__"+(p+1)+"_out"){1,c}.ber;
-
-            end
-
-            totalBer = matFile.("prms_compare_wdm__"+(p+1)+"_out"){1,end}.totalBer;
-
-            if totalBer > 0.2 && channelspacing == 400e9 &&  pol == "alternated"
-                disp("stopping here");
-                pause;
-            end
-
-            % ____________________________________
-            % Add value to warehouse at the correct position
-
-    
-
-            wh.addValueToStorage(ber,'ber',l,d,sgm,pol,p_in,p_out,pmd,gamma,realiz,numchannels, center_wavelength,channelspacing,random_zdw);
-            wh.getStoValue('ber',l,d,sgm,pol,p_in,p_out,pmd,gamma,realiz,numchannels, center_wavelength,channelspacing,random_zdw);
-            wh.addValueToStorage(totalBer,'totalBer',l,d,sgm,pol,p_in,p_out,pmd,gamma,realiz,numchannels,center_wavelength,channelspacing,random_zdw);
-
-        end
-
-    end
-    waitbar(num/numel(allMat),f,'Loading your data');
-end
-
-close(f)
-
-
-
-
-
-% 4) Hey! the warehouse is here and (hopefully) filled with data  :-)
-
-% Create a save dialog
-defaultDir = 'C:\Users\Silas\Documents\MATLAB\Raw_Cluster_Simulations\';
-defaultExt = '*.mat';
-[filename, pathname] = uiputfile(fullfile(defaultDir, defaultExt),'', 'wh.mat');
-
-% Check if the user pressed Cancel
-if isequal(filename, 0) || isequal(pathname, 0)
-    disp('Save operation canceled.');
-else
-    % Save the variable to the selected file
-    save(fullfile(pathname, filename), 'wh');
-    disp(['Variable "wh" saved to: ', fullfile(pathname, filename)]);
-end
-
-
-
-
-function matFileStructArray = getAllFilesInFolder(folderPath,extension)
-    % Get a list of all files in the current folder
-    currentFolderFiles = dir(fullfile(folderPath, '*'));
-
-    % Exclude '.' and '..' directories
-    currentFolderFiles = currentFolderFiles(~ismember({currentFolderFiles.name}, {'.', '..'}));
-
-    % Initialize the structure array for .mat files
-    matFileStructArray = struct('path', {}, 'name', {}, 'ext', {});
-
-    % Loop over each file in the current folder
-    for i = 1:length(currentFolderFiles)
-        currentFile = currentFolderFiles(i);
-
-        % Check if the current item is a file and has a .mat extension
-        if ~currentFile.isdir && endsWith(currentFile.name, extension, 'IgnoreCase', true)
-            % If it's a .mat file, add it to the structure array
-            [matFileStructArray(end + 1).path,matFileStructArray(end+1).name, matFileStructArray(end+1).ext] = fileparts(fullfile(folderPath, currentFile.name));
-        elseif currentFile.isdir
-            % If it's a directory, recursively call the function
-            subfolderPath = fullfile(folderPath, currentFile.name);
-            subfolderMatFiles = getAllFilesInFolder(subfolderPath,extension);
-
-            % Add .mat files from the subfolder to the structure array
-            matFileStructArray = [matFileStructArray, subfolderMatFiles];
-        end
-    end
-end
-
-

Classes/Warehouse_class/functions/fwm_plots/automate_JLT_plots.m

@@ -1,394 +0,0 @@
-
-
-%automate plots
-[file, path] = uigetfile("C:\Users\Silas\Documents\MATLAB\Datensätze\Raw_Cluster_Simulations\session_februar_24\wh_mi_nacht.mat");
-wh = load([path filesep file]);
-wh = wh.wh;
-
-plotJob = struct();
-width = 350;
-height = 200;
-plotJob.Position = [100 100 width 100+height];
-cols = cbrewer2("paired",12);
-plotJob.color = cols(1,:);
-plotJob.l = 2;
-plotJob.ch = 16;
-plotJob.d = 0;
-plotJob.sgm = 0;
-plotJob.pol = "copolarized";
-plotJob.p_in = 3;
-plotJob.gamma = 0;
-plotJob.pmd = 0;
-plotJob.channelspacing = 400e9;
-plotJob.randzdw = 0;
-
-
-plotJob.plot_ber_curve = 1;
-plotJob.plot_3dber_curve = 0;
-plotJob.plot_violin = 0;
-plotJob.plot_wavelength_sweep = 0;
-plotJob.plot_wavelength_sweep_failure_rate = 0;
-
-
-plotJob.dataStatArg = 'Lineplot with quartiles';
-plotJob.plotTypeArg = 'Lines';
-plotJob.displayname = 'a';
-plotJob.title = 'title';
-plotJob.figName = '16 Chann__';
-plotJob.xAxisLabel = 'ROP per Channel in dBm';
-plotJob.yAxisLabel = 'BER';
-
-% createbercurves(wh,plotJob)
-createviolinplots(wh,plotJob);
-% createsweepplots(wh,plotJob);
-
-
-%% 1
-function createbercurves(wh,plotJob)
-width = 1650;
-height = 400;
-s = 100;
-e = 100;
-
-cols = cbrewer2("paired",12);
-numRows = 2;
-numCols = 4;
-
-plotJob.figName = '16 Chann_200G';
-plotJob.channelspacing = 400e9;
-plotJob.ch = 16;
-Len = [2,2,2,2,10,10,10,10];
-Pol = ["copolarized","alternated","paired","copolarized","copolarized","alternated","paired","copolarized"];
-Title = ["Co Polarized","Alternating Pol. Interl.","Paired Pol. Interl.","Link Segmentation",];
-D = [0,0,0,3,0,0,0,3];
-Sgm = [0,0,0,1,0,0,0,1];
-
-colidx = [4,8,6];
-P_launch = [0,3,6];
-
-fig = figure('Name',plotJob.figName);
-fig.Position = plotJob.Position;
-fig.Units = "centimeters";
-fig.Position = [0 0 18 7];
-t = tiledlayout(numRows,numCols,'TileSpacing','compact','Padding','compact');
-for idx = 1:(numRows * numCols)
-    % Create subplot
-    % sp = subplot(numRows, numCols, idx);
-    nexttile;
-    plotJob.l = Len(idx);
-    plotJob.pol = Pol(idx);
-    plotJob.d = D(idx);
-    plotJob.sgm = Sgm(idx);
-    plotJob.randzdw = 1;
-
-    for i = 1:3
-
-        plotJob.p_in = P_launch(i);
-        plotJob.color = cols(colidx(i),:);
-        hold on
-        plotCurve(wh, plotJob);
-    
-    end
-    %
-    if idx ~= 1 && idx ~= 5 % For example, hide y-axis for subplot 1
-        set(gca, 'YTickLabel',[]); % Hide y-axis ticks and labels
-        set(gca,'YGrid','on');
-        set(gca, 'YLabel', []);
-    end
-    if idx ~= 5 && idx ~= 6 && idx ~= 7 && idx ~= 8
-        set(gca, 'XLabel', []);
-        set(gca, 'XTickLabel', []);
-
-    end
-    grid on
-
-    g = gca;
-    pos = g.Position;
-    if idx <= 4
-        title(Title(idx),'FontSize',8);
-        %              a = annotation('textbox', pos-[0.0020   -0.1434    0.0947    0.3121], 'String', "FEC: 3.8e-3","LineStyle","none","FontSize",8,'LineWidth',0.1,'Interpreter','latex',"FontUnits","points",'Vert','middle','FitBoxToText','on');
-        %              a = annotation('textbox', pos, 'String', "2 km","LineStyle","none","FontSize",8,'LineWidth',0.1,'Interpreter','latex',"FontUnits","points",'Vert','bottom','FitBoxToText','on');
-    else
-        %              a = annotation('textbox', pos, 'String', "FEC: 3.8e-3","LineStyle","none","FontSize",8,'LineWidth',0.1,'Interpreter','latex',"FontUnits","points",'Vert','middle','FitBoxToText','on');
-        %              a = annotation('textbox', pos, 'String', "10 km","LineStyle","none","FontSize",8,'LineWidth',0.1,'Interpreter','latex',"FontUnits","points",'Vert','bottom','FitBoxToText','on');
-    end
-
-end
-
-% Create textbox
-annotation(fig,'textbox',...
-    [0.0696078431372547 0.246851385390432 0.0656862745098043 0.0453400503778337],...
-    'String','10 km',...
-    'LineStyle','none',...
-    'Interpreter','latex',...
-    'FontSize',8,...
-    'FitBoxToText','off');
-
-% Create textbox
-annotation(fig,'textbox',...
-    [0.288235294117647 0.239294710327459 0.0656862745098043 0.0453400503778338],...
-    'String','10 km',...
-    'LineStyle','none',...
-    'Interpreter','latex',...
-    'FontSize',8,...
-    'FitBoxToText','off');
-
-% Create textbox
-annotation(fig,'textbox',...
-    [0.516666666666666 0.241813602015117 0.0656862745098042 0.0453400503778338],...
-    'String','10 km',...
-    'LineStyle','none',...
-    'Interpreter','latex',...
-    'FontSize',8,...
-    'FitBoxToText','off');
-
-% Create textbox
-annotation(fig,'textbox',...
-    [0.742156862745097 0.236775818639802 0.0656862745098042 0.0453400503778339],...
-    'String','10 km',...
-    'LineStyle','none',...
-    'Interpreter','latex',...
-    'FontSize',8,...
-    'FitBoxToText','off');
-
-% Create textbox
-annotation(fig,'textbox',...
-    [0.071996471804854 0.578899159967702 0.0656862745098039 0.0453400503778341],...
-    'String',{'2 km'},...
-    'LineStyle','none',...
-    'Interpreter','latex',...
-    'FontSize',8,...
-    'FitBoxToText','off');
-
-% Create textbox
-annotation(fig,'textbox',...
-    [0.29596893566113 0.576253721089996 0.0656862745098041 0.0453400503778341],...
-    'String',{'2 km'},...
-    'LineStyle','none',...
-    'Interpreter','latex',...
-    'FontSize',8,...
-    'FitBoxToText','off');
-
-% Create textbox
-annotation(fig,'textbox',...
-    [0.524883914268912 0.580785315705112 0.0656862745098041 0.0453400503778341],...
-    'String',{'2 km'},...
-    'LineStyle','none',...
-    'Interpreter','latex',...
-    'FontSize',8,...
-    'FitBoxToText','off');
-
-% Create textbox
-annotation(fig,'textbox',...
-    [0.753758338909501 0.581291504465311 0.0656862745098037 0.0453400503778341],...
-    'String',{'2 km'},...
-    'LineStyle','none',...
-    'Interpreter','latex',...
-    'FontSize',8,...
-    'FitBoxToText','off');
-
-a=sgtitle(['N=',num2str(plotJob.ch),'; $\Delta f_{\mathrm{ch}}$= ',num2str(plotJob.channelspacing*1e-9),' GHz'],'FontSIze',10);
-a.Interpreter = "latex";
-
-lgd = legend('$P_{\mathrm{in}}=0$ dBm','$P_{\mathrm{in}}=3$ dBm','$P_{\mathrm{in}}=6$ dBm','Interpreter','latex');
-lgd.NumColumns = 3;
-lgd.Layout.Tile = 'south';
-
-copygraphics(t,'BackgroundColor','none');
-end
-
-%% 2
-function createviolinplots(wh,plotJob)
-width = 350;
-height = 200;
-s = 100;
-e = 100;
-
-cols = cbrewer2("paired",12);
-numRows = 1;
-numCols = 4;
-
-
-plotJob.ch = 16;
-plotJob.p_in = 3;
-plotJob.randzdw = 0;
-
-Pol = ["copolarized","copolarized","alternated","paired",];
-Title = ["Co Pol.","Link Segmentation","Paired Pol. Interl.","Alternating Pol. Interl."];
-D = [0,3,0,0];
-Sgm = [0,1,0,0];
-
-colidx = [2];
-Len = [2];
-
-plotJob.figName = ['_Violin',num2str(plotJob.ch),' Channels; ',num2str(plotJob.channelspacing*1e-9),' GHz; ',num2str(plotJob.p_in),' dBm; randomized: ', num2str(plotJob.randzdw)];
-fig = figure('Name',plotJob.figName);
-fig.Position = plotJob.Position;
-fig.Units = "centimeters";
-fig.Position = [0 0 18 7];
-
-t = tiledlayout(numRows,numCols,'TileSpacing','compact','Padding','compact');
-for idx = 1:(numRows * numCols)
-    % Create subplot
-    %subplot(numRows, numCols, idx);
-    nexttile;
-
-    plotJob.pol = Pol(idx);
-    plotJob.d = D(idx);
-    plotJob.sgm = Sgm(idx);
-
-    for i = 1
-
-        plotJob.color = cols(colidx(i),:);
-        plotJob.l = Len(i);
-        hold on
-        plotViolin(wh, plotJob);
-
-    end
-
-    if idx ~= 1 % For example, hide y-axis for subplot 1
-        %set(gca, 'YTickLabel',[]); % Hide y-axis ticks and labels
-        set(gca, 'YGrid','on');
-        set(gca, 'YLabel', []);
-    end
-
-    if idx <= 4
-        title(Title(idx));
-    end
-
-
-end
-
-% Create textbox
-annotation(fig,'textbox',...
-    [0.300019607843137 0.816120906801009 0.108803921568628 0.0906801007556676],...
-    'String','$P_{\mathrm{in}}=3$ dBm',...
-    'LineStyle','none',...
-    'Interpreter','latex',...
-    'FontSize',8,...
-    'FitBoxToText','off');
-
-% Create textbox
-annotation(fig,'textbox',...
-    [0.530411764705882 0.81612090680101 0.108803921568628 0.0906801007556676],...
-    'String','$P_{\mathrm{in}}=3$ dBm',...
-    'LineStyle','none',...
-    'Interpreter','latex',...
-    'FontSize',8,...
-    'FitBoxToText','off');
-
-% Create textbox
-annotation(fig,'textbox',...
-    [0.755901960784313 0.816120906801011 0.108803921568628 0.0906801007556676],...
-    'String','$P_{\mathrm{in}}=3$ dBm',...
-    'LineStyle','none',...
-    'Interpreter','latex',...
-    'FontSize',8,...
-    'FitBoxToText','off');
-
-% Create textbox
-annotation(fig,'textbox',...
-    [0.0696274509803918 0.584382871536529 0.108803921568628 0.0906801007556676],...
-    'String','$P_{\mathrm{in}}=3$ dBm',...
-    'LineStyle','none',...
-    'Interpreter','latex',...
-    'FontSize',8,...
-    'FitBoxToText','off');
-
-copygraphics(t,'BackgroundColor','none');
-%     lgd = legend('$P_{\mathrm{in}}=0$ dBm','$P_{\mathrm{in}}=3$ dBm','$P_{\mathrm{in}}=6$ dBm','Interpreter','latex');
-%     lgd.NumColumns = 3;
-%     lgd.Layout.Tile = 'south';
-
-end
-
-%% 3
-function createsweepplots(wh,plotJob)
-
-width = 650;
-height = 200;
-s = 100;
-e = 100;
-plotJob.Position = [0 0 width e+height];
-
-cols = cbrewer2("paired",12);
-numRows = 1;
-numCols = 4;
-
-
-plotJob.channelspacing = 200e9;
-plotJob.ch = 16;
-plotJob.randzdw = 1;
-plotJob.l = 10;
-
-plotJob.p_in = 3;
-
-Pol = ["copolarized","alternated","paired","copolarized"];
-Title = ["Co Polarized","Alternating Pol. Interl.","Paired Pol. Interl.","Link Segmentation",];
-D = [0,0,0,3];
-Sgm = [0,0,0,1];
-Channelspacing = [200e9, 200e9];
-PlotTypeArg = ["--","-"];
-colidx = [6,8,2,4];
-Len = [2,10];
-
-plotJob.figName = [num2str(plotJob.ch),num2str(plotJob.channelspacing*1e-9),num2str(plotJob.p_in),'...'];
-plotJob.figName = "10km 400ghz";
-
-
-
-
-fig = findall(groot, 'Type', 'figure', 'Name', plotJob.figName);
-
-if isvalid(fig)
-    figure(fig)
-    fig = get(fig);
-    AxesMain =  fig.CurrentAxes;
-    hold on
-else
-    fig = figure('name',char(plotJob.figName));
-    AxesMain    = gca;
-    hold on
-end
-
-fig.Position = plotJob.Position;
-fig.Units = "centimeters";
-fig.Position = [0 0 18 7];
-
-for j = 1
-
-    plotJob.channelspacing = Channelspacing(j);
-    plotJob.plotTypeArg = PlotTypeArg(j);
-
-    for idx = 1:4
-
-        plotJob.color = cols(colidx(idx),:);
-        plotJob.pol = Pol(idx);
-        plotJob.d = D(idx);
-        plotJob.sgm = Sgm(idx);
-        plotJob.displayname = [char(plotJob.pol)];
-        hold on
-        plotBerVsZdwFailureRate(wh, plotJob);
-
-    end
-end
-legend('Location', 'southoutside', 'Orientation', 'horizontal');
-
-
-%plot channel positions
-hold on
-chpos = calcWavelengthPlan(plotJob.ch, plotJob.channelspacing, 1310);
-xline(chpos,'LineWidth',2,'Alpha',0.4,'HandleVisibility','off');
-
-chpos = calcWavelengthPlan(plotJob.ch, plotJob.channelspacing, chpos(4));
-xline(chpos,'LineWidth',2,'LineStyle','--','Alpha',0.1,'HandleVisibility','off');
-
-title(['N=',num2str(plotJob.ch),' $\Delta f_{\mathrm{ch}}$= ',num2str(plotJob.channelspacing*1e-9),' GHz'],'FontSize',10,'Interpreter','latex');
-
-
-
-
-end
-
-
-

Classes/Warehouse_class/functions/fwm_plots/dispersion_only.m

@@ -1,83 +0,0 @@
-%automate plots
-[file, path] = uigetfile("C:\Users\Silas\Documents\MATLAB\Datensätze\Raw_Cluster_Simulations\");
-wh = load([path filesep file]);
-wh = wh.wh;
-
-plotJob = struct();
-width = 350;
-height = 200;
-plotJob.Position = [100 100 width 100+height];
-cols = cbrewer2("paired",12);
-plotJob.color = cols(1,:);
-plotJob.l = 1;
-plotJob.ch = 1;
-
-plotJob.sgm = 1;
-plotJob.pol = "copolarized";
-plotJob.p_in = 3;
-plotJob.gamma = 0.0023;
-plotJob.pmd = 0.1;
-plotJob.channelspacing = 400e9;
-plotJob.randzdw = 0;
-
-
-plotJob.plot_ber_curve = 1;
-plotJob.plot_3dber_curve = 0;
-plotJob.plot_violin = 0;
-plotJob.plot_wavelength_sweep = 0;
-plotJob.plot_wavelength_sweep_failure_rate = 0;
-
-
-plotJob.dataStatArg = 'Lineplot with quartiles';
-plotJob.plotTypeArg = 'Lines';
-plotJob.displayname = 'a';
-plotJob.title = 'title';
-plotJob.figName = '1 Chann__';
-plotJob.xAxisLabel = 'ROP per Channel in dBm';
-plotJob.yAxisLabel = 'BER';
-
-
-plotJob.d = 0;
-
-xAxis = wh.parameter.p_out.values;
-D = wh.parameter.dispersion.values;
-
-figure()
-ber_ = [];
-for d_ = 0:39
-    if d_ == 0
-        plotJob.sgm = 0;
-        ber_(d_+1,:) = wh.getStoValue('ber',plotJob.l,d_,plotJob.sgm,string(plotJob.pol),plotJob.p_in,xAxis,plotJob.pmd,plotJob.gamma,1,plotJob.ch,1310,plotJob.channelspacing,plotJob.randzdw)';
-    else
-        plotJob.sgm = 1;
-        ber_(d_+1,:) = wh.getStoValue('ber',plotJob.l,d_,plotJob.sgm,string(plotJob.pol),plotJob.p_in,xAxis,plotJob.pmd,plotJob.gamma,1,plotJob.ch,1310,plotJob.channelspacing,plotJob.randzdw)';
-    end
-    hold on
-    plot(xAxis,ber_(d_+1,:))
-    set(gca,'yscale','log');
-end   
-yline(3.8e-3);
-
-
-hdfec = 3.8e-3.*ones(size(xAxis));
-for i = 1:size(ber_,1)
-    ber_series = ber_(i,:);
-    a = InterX([hdfec(:)';xAxis],[ber_series;xAxis]);
-    cross(i) = a(2);
-end
-
-col = cbrewer2('Paired',8);
-figure()
-plot(D,cross,'Marker','o','MarkerSize',5,'MarkerEdgeColor',[1,1,1],'MarkerFaceColor',col(2,:),'Color',col(1,:),'LineWidth',1);
-grid minor
-xlabel('Accumulated Dispersion')
-ylabel('Required ROP to reach FEC limit in dB')
-line([D(16),D(16)],[-10,cross(16)],'linestyle','--')
-line([0,D(16)],[cross(16),cross(16)],'linestyle','--')
-
-line([D(29),D(29)],[-10,cross(29)],'linestyle','--')
-line([0,D(29)],[cross(29),cross(29)],'linestyle','--')
-
-
-
-

Classes/Warehouse_class/functions/fwm_plots/dispersion_validation_miniskript.m

@@ -1,52 +0,0 @@
-
-wh = load("C:\Users\Silas\Documents\MATLAB\Datensätze\Raw_Cluster_Simulations\session_dispersion_validation\wh_with_variation.mat");
-wh = wh.wh;
-
-lambda = 1295;
-
-figure(3)
-plot(xAxis,getber(wh,lambda),'DisplayName',['w:', num2str(lambda), 'variation']);
-yline(3.8e-3,'HandleVisibility','off');
-legend
-set(gca,'yscale','log');
-grid(gca,'on');
-grid(gca,'minor');
-grid minor
-fontsize(gca,8,"points")
-fig.Units = "centimeters";
-fig.Position = [2 2 8.5 7];
-set(gca,'TickLabelInterpreter','latex')
-ylim([1e-5,0.5]);
-xlim([min(xAxis),-3]);
-
-wh = load("C:\Users\Silas\Documents\MATLAB\Datensätze\Raw_Cluster_Simulations\session_dispersion_validation\wh_no_variation.mat");
-wh = wh.wh;
-
-
-figure(3)
-hold on
-plot(xAxis,getber(wh,lambda),'DisplayName',['w:', num2str(lambda), 'no variation']);
-yline(3.8e-3,'HandleVisibility','off');
-legend
-set(gca,'yscale','log');
-grid(gca,'on');
-grid(gca,'minor');
-grid minor
-fontsize(gca,8,"points")
-fig.Units = "centimeters";
-fig.Position = [2 2 8.5 7];
-set(gca,'TickLabelInterpreter','latex')
-ylim([1e-5,0.5]);
-xlim([min(xAxis),-3]);
-
-
-function ber = getber(wh,lambda)
-    realization = wh.parameter.realization.values(1:end);
-    xAxis = wh.parameter.p_out.values;
-    ber = [];
-    for xl = 1:numel(xAxis)
-        p_out = xAxis(xl);
-        temp = wh.getStoValue('ber',10,0,0,"copolarized",3,p_out,0.1,0.0023,realization,1,lambda,400e9,1);
-        ber(xl) = mean(temp,'all');
-    end
-end

Classes/Warehouse_class/functions/fwm_plots/generatePlots.m

@@ -1,61 +0,0 @@
-function generatePlots(wh,plotJob)
-
-% 0) Test for valid query:
-p_out = wh.parameter.p_out.values(1);
-realization = 9;
-
-
-if 1 %~isempty(wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.ch,1310))
-    % test violin
-
-    baseName = plotJob.figName;
-
-    width = 350;
-    height = 200;
-    s = 100;
-    e = 100;
-
-    if plotJob.plot_ber_curve
-        plotJob.Position = [100 100 width e+height];
-        plotJob.figName = [baseName, ' zdwvsber'];
-        plotCurve(wh, plotJob);
-    end
-
-    if plotJob.plot_3dber_curve
-        plotJob.Position = [100 100 width e+height];
-        plotJob.figName = [baseName, ' zdwvsber'];
-        plot3dCurve(wh, plotJob);
-    end
-
-    if plotJob.plot_wavelength_sweep
-        plotJob.Position = [100 100 width e+height];
-        plotJob.figName = [baseName, ' zdwvsber'];
-        plotBerVsZDW(wh, plotJob);
-    end
-
-    if plotJob.plot_wavelength_sweep_failure_rate
-        plotJob.Position = [100 100 width e+height];
-        plotJob.figName = [baseName, ' zdwvsber'];
-        plotBerVsZdwFailureRate(wh, plotJob);
-    end
-
-   if plotJob.plot_violin
-        plotJob.Position = [s+width 100 width  e+height];
-        plotJob.figName = [baseName, ' violin'];
-        plotViolin(wh, plotJob);
-   end
-
-
-    if 0
-    %2) plotHistogram
-        plotJob.Position = [s+2*width 100 width e+height];
-        plotJob.figName = [baseName, ' FEC crossing'];
-        plotHistogram(wh,plotJob)
-    end
-
-
-else
-    warndlg('The requested Datapoint is not available... This can occur for some edgecase constellations... ')
-end
-
-end

Classes/Warehouse_class/functions/fwm_plots/plot3dCurve.m

@@ -1,248 +0,0 @@
-function plotCurve(wh,plotJob)
-%PLOTCURVE Summary of this function goes here
-%   Detailed explanation goes here
-
-fig = findall(groot, 'Type', 'figure', 'Name', plotJob.figName);
-
-if isvalid(fig)
-    figure(fig)
-    fig = get(fig);
-    AxesMain =  fig.CurrentAxes;
-    hold on
-else
-    fig = figure('name',char(plotJob.figName));
-    AxesMain    = gca;
-    hold on
-end
-
-col = plotJob.color;
-
-% we want to fetch all realizations
-if plotJob.pmd == 0
-    realization = 499;
-    % realization = 0:7;
-else
-    realization = wh.parameter.realization.values(1:end);
-end
-realization = wh.parameter.realization.values(1:end);
-% get all xAxis values
-xAxis = wh.parameter.p_out.values;
-
-% Fetch Data from Warehouse
-for xl = 1:numel(xAxis)
-    p_out = xAxis(xl);
-    if string(plotJob.dataStatArg) == "Worst"
-        temp = wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.ch,1310,plotJob.channelspacing,plotJob.randzdw);
-        temp = removeZeros(temp);
-        ber(xl) = max(temp,[],'all');
-        linew = 1.0;
-        markersz = 3;
-        linestyle = '-';
-    elseif string(plotJob.dataStatArg) == "AVG"
-        temp = wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.ch,1310,plotJob.channelspacing,plotJob.randzdw);
-        temp = removeZeros(temp);
-        ber(xl) = mean(temp,'all');
-        linew = 1.0;
-        markersz = 3;
-        linestyle = '-';
-    elseif string(plotJob.dataStatArg) == "Best"
-        temp = wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.ch,1310,plotJob.channelspacing,plotJob.randzdw);
-        temp = removeZeros(temp);
-        ber(xl) = min(temp,[],'all');
-        linew = 1.0;
-        markersz = 3;
-        linestyle = '-';
-
-    elseif string(plotJob.dataStatArg) == "All Channels; mean(PMD Realizations)"
-        temp = wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.ch,1310,plotJob.channelspacing,plotJob.randzdw);
-        temp = removeZeros(temp);
-
-        ber(:,xl) = mean(temp,1,"omitnan").';
-
-        linew = 1;
-        markersz = 3;
-        linestyle = '--';
-
-    elseif string(plotJob.dataStatArg) == "Lineplot with quartiles"
-
-        temp = wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.ch,1310,plotJob.channelspacing,plotJob.randzdw);
-        temp = removeZeros(temp);
-        ber(xl) = mean(temp,"all","omitnan").';
-
-
-        upperq(xl) = quantile(temp,0.9,"all");
-        lowerq(xl) = quantile(temp,0.1,"all");
-        
-        if lowerq(xl) == 0
-            lowerq(xl) = lowerq(xl-1);
-        end
-        %         upperq(xl) =  0.5*std(tmp,1,'all','omitnan');
-        %         lowerq(xl) =  0.5*std(tmp,1,'all','omitnan');
-
-        %         upperq(xl) = max(dataNoNans);
-        %         lowerq(xl) = min(dataNoNans);
-
-        %                 upperq(xl) = mean(tmp,"all","omitnan") + 1.96 * (std(tmp,1,'all','omitnan')/sqrt(numel(tmp)));
-        %                 lowerq(xl) = mean(tmp,"all","omitnan") - 1.96 * (std(tmp,1,'all','omitnan')/sqrt(numel(tmp)));
-
-        linew = 1;
-        markersz = 1;
-        linestyle = '-';
-
-    elseif string(plotJob.dataStatArg) == "All Channels ;All PMD Realizations"
-
-        tmp = reshape(wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.ch,1310,plotJob.channelspacing,plotJob.randzdw).',[],1);
-        ber(1:size(tmp,1),xl) = tmp;
-
-        linew = 0.3;
-        markersz = 2;
-        linestyle = ':';
-    end
-end
-
-xAxis = xAxis;
-
-% Plot Data
-if string(plotJob.plotTypeArg) == "Scatter"
-    for rlz = 1:size(ber,1)
-
-        if rlz < size(ber,1)
-            scatter(xAxis,ber(rlz,:),markersz,'MarkerEdgeColor',col,'MarkerFaceColor',col,'Parent',  AxesMain(1),'HandleVisibility','off');
-        else
-            scatter(xAxis,ber(rlz,:),markersz,'MarkerEdgeColor',col,'MarkerFaceColor',col,'Parent',  AxesMain(1),'DisplayName',[plotJob.dataStatArg,' Ch.: ', plotJob.displayname]);
-        end
-
-    end
-
-elseif string(plotJob.plotTypeArg) == "Lines"
-
-    %         if ~anynan(ber)
-    %             [xAxis,ber] = interpCurve(xAxis, ber);
-    %         end
-
-    for rlz = 1:size(ber,1)
-% 
-        if (string(plotJob.dataStatArg) == "All Channels ;All PMD Realizations")||(string(plotJob.dataStatArg) == "All Channels; mean(PMD Realizations)")
-            ch = mod(rlz,plotJob.ch);
-            if ch == 0; ch = plotJob.ch; end
-        else
-            ch = plotJob.dataStatArg;
-        end
-
-        if rlz <= size(ber,1)
-            
-
-            s = plot3(xAxis,repmat(ch,1,numel(xAxis)),ber(rlz,:),linestyle,'Marker',"o",'MarkerSize',markersz,'MarkerFaceColor',plotJob.color,'LineWidth',linew,'Color',col,'Parent',  AxesMain,'HandleVisibility','off');
-            
-            s.DataTipTemplate.Interpreter = "latex";
-            s.DataTipTemplate.DataTipRows(1).Label = "Ch: ";
-            s.DataTipTemplate.DataTipRows(1).Value = repmat(ch,size(ber));
-            s.DataTipTemplate.DataTipRows(1);
-            s.DataTipTemplate.DataTipRows(2) = [];
-
-            
-        else
-
-            if string(plotJob.dataStatArg) == "Lineplot with quartiles"
-                [hl,hp] = boundedline(xAxis,ber(rlz,:),([(ber(rlz,:)-lowerq(rlz,:));upperq(rlz,:)-ber(rlz,:)]'),'-*', 'alpha','Color',col,'transparency', 0.2);
-                hp.LineWidth = 1.2;
-
-                ho = outlinebounds(hl,hp);
-                set(ho, 'linestyle', ':', 'color', col);
-            else
-
-                s = plot(xAxis,ber(rlz,:),linestyle,'Marker',"o",'MarkerFaceColor',col,'MarkerSize',markersz,'LineWidth',linew,'Color',col,'Parent',  AxesMain,'DisplayName',[plotJob.displayname]);
-                
-                s.DataTipTemplate.Interpreter = "latex";
-                s.DataTipTemplate.DataTipRows(1).Label = "Ch: ";
-                s.DataTipTemplate.DataTipRows(1).Value = repmat(string(ch),size(ber));
-                s.DataTipTemplate.DataTipRows(1)
-                s.DataTipTemplate.DataTipRows(2) = [];
-            end
-
-        end
-
-    end
-end
-
-% Draw FEC Threshold Line
-%get x data of first children:
-%get all linear Values
-if 0
-    linear   = wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,wh.parameter.p_out.values,0,0,499,"symmetric");
-    linear   = mean(linear,2);
-    lincurve = findall(AxesMain, 'Type', 'line','DisplayName','Linear Baseline');
-    %
-    if isempty(lincurve)
-        xdata = AxesMain.Children(1).XData;
-        hdfec = 3.8e-3.*ones(size(xdata));
-        plot(xdata,linear,'-','Marker',"o",'MarkerSize',3,'LineWidth',4,'Color',[0.6400    0.6400    0.6400],'MarkerFaceColor',[0.6400    0.6400    0.6400],'Parent', AxesMain,'DisplayName','Linear Baseline');
-        %h = get(gca,'Children');
-        %set(gca,'Children',[h(2) h(1)])
-    end
-end
-
-feccurve = findall(AxesMain, 'Type', 'line','DisplayName','FEC $3.8*10^{-3}$');
-%
-if isempty(feccurve)
-    xdata = xAxis;
-    hdfec = 3.8e-3.*ones(size(xdata));
-    for ch = 1:plotJob.ch
-        plot3(xdata,repmat(ch,1,numel(xAxis)),hdfec,':','MarkerSize',4,'Color','black','MarkerFaceColor','black','LineWidth',0.5,'Parent', AxesMain,'DisplayName','FEC $3.8*10^{-3}$','HandleVisibility','off');
-    end
-    %h = get(gca,'Children');
-    %set(gca,'Children',[h(2) h(1)])
-end
-
-% Figure Settings
-%title(AxesMain,plotJob.title,"Interpreter","none");
-
-xlabel(AxesMain,plotJob.xAxisLabel,"Interpreter","none");
-
-ylabel(AxesMain,plotJob.yAxisLabel,"Interpreter","none");
-
-set(AxesMain,'zscale','log');
-
-grid(AxesMain,'on');
-
-grid(AxesMain,'minor');
-
-grid minor
-
-view(AxesMain,[42.0619302949062 23.4176470588235]);
-%legend(AxesMain);
-
-fontsize(AxesMain,8,"points")
-fontname(AxesMain,"Arial")
-
-fig.Position = plotJob.Position;
-fig.Units = "centimeters";
-fig.Position = [2 2 8.5 7];
-
-set(AxesMain,'TickLabelInterpreter','none')
-
-set(AxesMain.Legend,'Interpreter','none')
-% set(gcf,'Units','centimeters')
-% set(gcf,'Position',[2 2 9 4.5])
-
-zlim([1e-4,0.3]);
-
-xlim([min(xAxis),-3]);
-
-
-annotation('textbox', [0.125, 0.32, 0.1, 0.1], 'String', "FEC 3.8e-3","LineStyle","none","FontSize",8,"FontUnits","points")
-
-
-hold off
-
-
-end
-
-
-function vec = removeZeros(vec)
-    % Find rows that contain only zeros
-    rows_to_remove = all(vec == 0, 2);
-    
-    % Remove rows with only zeros
-    vec(rows_to_remove, :) = [];
-end

Classes/Warehouse_class/functions/fwm_plots/plotBerVsZDW.m

@@ -1,300 +0,0 @@
-function plotBerVsZDW(wh,plotJob)
-
-    fig = findall(groot, 'Type', 'figure', 'Name', plotJob.figName);
-    
-    if isvalid(fig)
-        figure(fig)
-        fig = get(fig);
-        AxesMain =  fig.CurrentAxes;
-        hold on
-    else
-        fig = figure('name',char(plotJob.figName));
-        AxesMain    = gca;
-        hold on
-    end
-    
-    
-    col = plotJob.color;
-    
-    % we want to fetch all realizations
-    realization = wh.parameter.realization.values(1:end);
-    
-    % get all xAxis values
-    xAxis = wh.parameter.p_out.values;
-    
-    % get all center wavelengths
-    wavelengths = wh.parameter.center_wavelength.values;
-
-%     totber = NaN(numel(realization),1,numel(xAxis),numel(wavelengths));
-%     ber = zeros(numel(realization),plotJob.ch,numel(xAxis),numel(wavelengths));
-    
-    %get BER values for query
-    for w = 2:numel(wavelengths)
-        
-        for xl = 1:numel(xAxis)
-    
-            c_wavelen = wavelengths(w);
-            p_out = xAxis(xl);
-    
-            % dim1 : realiz; dim2: channels, dim3: rop, dim4, c_wavelength
-            temp  = removeZeros(wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.ch,c_wavelen,plotJob.channelspacing,plotJob.randzdw));
-            ber(1:size(temp,1),:,xl,w)     = temp;
-            
-        end
-    end
-    
-    hdfec = 3.8e-3.*ones(size(xAxis));
-    zdw_ = [];
-    zdw_chann = [];
-    zdw_tot = [];
-    cf_tot = [];
-    cf_ = [];
-    S = [];
-    Stot = [];
-    S_chann = [];
-    cf_chann = [];
-    
-    cnt  = 0;
-    %get fec thresholds
-    %linear = squeeze(linear);
-    for c_wavelen = 1:size(ber,4)
-        for realiz = 1:size(ber,1)
-            for chann = 1:size(ber,2)
-            
-                %finde Schnittpunkt zwischen FEC und BER Kurve
-                temp_ber = squeeze(ber(realiz,chann,:,c_wavelen)).';
-                if ~all(temp_ber == 0)
-                   %nur wenn nicht alles nullen sind
-                    crossing_ch = InterX([hdfec;xAxis],[temp_ber;xAxis]);
-                else
-                    continue
-                end
-
-                %Req. FEC Ergebnis einsortieren
-                if ~isempty(crossing_ch)
-                    if crossing_ch(2) == 0
-                        print("d")
-                    end
-                    S(realiz,chann,c_wavelen) = crossing_ch(2);
-                else
-                    S(realiz,chann,c_wavelen) = -1;
-                    cnt = cnt +1;
-                end
-        
-        
-            end
-        end
-    end
-    
-    temp_max = -inf;
-    for i = 1:plotJob.ch
-        hold on
-        %S:: 1.dim: realiz; 2.dim: channel; 3.dim: center wavelength
-        %squeeze a channel:
-        temp_data = squeeze(S(:,i,:));
-
-        %remove realizations that have no entry (only zero)
-        temp_data = removeZeros(temp_data);
-
-        %replace zeros with NAN (e.g. for the wavelengths that have missing realizations)
-        temp_data(temp_data==0) = NaN;
-
-        %plot required ROP for channel and all realizations that cross the
-        %FEC limit
-        scatter(wavelengths,temp_data ,5,plotJob.color,'Marker','.');
-
-%         %plot mean per channel
-%         temp_mean = mean(temp_data,'omitnan');
-%         hold on
-%         plot(wavelengths,temp_mean,'Marker','*');
-
-        %get max overall value 
-        temp_max = max(temp_max,max(temp_data));
-    end
-
-
-    
-    %plot mean overall
-    mean_overall = squeeze(mean(S,2));
-    mean_overall(mean_overall==0) = NaN;
-    %mean_overall(mean_overall==-1) = NaN;
-    mean_overall=mean(mean_overall,1,'omitnan');
-    plot(wavelengths,mean_overall,'Color',plotJob.color);
-
-    %plot max overall
-    scatter(wavelengths,temp_max ,35,plotJob.color,'Marker','v');
-
-
-    %plot channel positions
-    hold on
-    chpos = calcWavelengthPlan(plotJob.ch, 400e9, 1310);
-    xline(chpos,'LineWidth',2,'Alpha',0.2);
-    chpos = calcWavelengthPlan(plotJob.ch, 400e9, chpos(4));
-    xline(chpos,'LineWidth',2,'Alpha',0.2);
-
-    fig.Position = plotJob.Position;
-
-    ylabel('Penalty in dB');
-    xlabel('Wavelength in nm');
-
-    xlim([min(wavelengths),max(wavelengths) ]);
-
-    grid minor;
-    set(gca, 'color', 'none');
-    legend = [];
-    
-    fontsize(AxesMain,8,"points")
-    
-    fig.Position = plotJob.Position;
-    fig.Units = "centimeters";
-    fig.Position = [2 2 8.5 7];
-    
-    set(AxesMain,'TickLabelInterpreter','latex')
-    
-    set(AxesMain.Legend,'Interpreter','latex')
-    
-    
-    
-    
-    
-    
-    % 
-    % 
-    % 
-    % 
-    % 
-    % distinct_cf = unique(cf_chann);
-    % 
-    % for i = 1:length(distinct_cf)
-    %     indices = find(cf_chann==distinct_cf(i));
-    %     cf(i) = distinct_cf(i);
-    %     worst_fec_cross(i) = max(S_chann(indices));
-    %     avg_fec_cross(i) = mean(S_chann(indices));
-    % end
-    % 
-    % avg_fec_cross = smooth(avg_fec_cross,5);
-    % 
-    % figure(224)
-    % hold on
-    % scatter(cf_,S,10.*abs(S-mean(S)).*ones(size(S)),'DisplayName',['AVG'],'MarkerEdgeColor',col,'MarkerFaceColor',col,'Marker','.');
-    % hold on
-    % scatter(cf(2:end),worst_fec_cross(2:end),15,'DisplayName',['Worst'],'MarkerEdgeColor',col,'MarkerFaceColor',col,'Marker','.','HandleVisibility','off');
-    % plot(cf(2:end),avg_fec_cross(2:end),'DisplayName',['AVG'],'LineWidth',1,'LineStyle','-','Color',col,'Marker','none','MarkerFaceColor',col,'MarkerSize',2);
-    % plot(cf(2:end),worst_fec_cross(2:end),'DisplayName',['AVG'],'LineWidth',0.5,'LineStyle','-','Color',col,'Marker','none','MarkerFaceColor',col,'MarkerSize',2);
-    % ghzgrid = hz2nm(nm2hz(1310)+[0:1:20].*400e9);
-    % set(gca,'xtick',sort(ghzgrid))
-    % xlim([min(cf(cf~=0)), 1310.1]);
-    % 
-    % 
-    % % With matlab internal errorbar function...
-    % figure(221)
-    % %plot(cf,avg_fec_cross,'DisplayName',['AVG'],'LineWidth',1,'Color',col,'Marker','none','MarkerFaceColor',col,'MarkerSize',2);
-    % hold on
-    % %plot(cf_tot,min(S_chann(1:length(cf_tot),:),[],2),'DisplayName',['AVG'],'LineWidth',1,'LineStyle',':','Color',col,'Marker','^','MarkerFaceColor',col,'MarkerSize',2);
-    % plot(cf,worst_fec_cross,'DisplayName',['AVG'],'LineWidth',2,'LineStyle','-','Color',col,'Marker','none','MarkerFaceColor',col,'MarkerSize',2);
-    % ghzgrid = hz2nm(nm2hz(1310)+[0:1:20].*400e9);
-    % set(gca,'xtick',sort(ghzgrid))
-    % xlim([1302, 1310.1]);
-    % xline(ghzgrid,'LineStyle',':','Color',[.7 .7 .7]);
-    
-    
-    %with 
-    % Stot = movmean(Stot,5);
-    % figure(222)
-    % [hl,hp] = boundedline(cf_tot,Stot,[(Stot'-min(S_chann(1:length(cf_tot),:),[],2)),(max(S_chann(1:length(cf_tot),:),[],2)-Stot')], 'alpha','Color',col,'transparency', 0.05);
-    % ho = outlinebounds(hl,hp);
-    % set(ho, 'linestyle', ':', 'color', col, 'marker', '.','linewidth',0.5);
-    % hold on
-    % 
-    % ghzgrid = hz2nm(nm2hz(1310)+[0:1:12].*400e9);
-    % xline(ghzgrid);
-    
-    
-    
-    
-    %plot the total ber
-    
-    % %figure(22);
-    % hold on;
-    % b= movmean(Stot,3);
-    % plot(AxesMain,cf_tot,b,'LineWidth',2,'DisplayName',[plotJob.dataStatArg,' Ch.: ', plotJob.displayname],'Color',col,'Marker','o');
-    % hold on
-    
-    
-    
-    
-    
-    %scatter(AxesMain,zdw_,S,'Marker','+','MarkerEdgeColor',col,'MarkerFaceAlpha',0.4,'MarkerEdgeAlpha',0.4,'LineWidth',0.5,'Parent',  AxesMain(1),'DisplayName',[plotJob.dataStatArg,' Ch.: ', plotJob.displayname]);
-    %ylim(AxesMain,[-9.3 -7]);
-    
-    % for i = 1:numel(chp)
-    %     hold on
-    %     xline(AxesMain,chp(i),'Color',colr(i,:),'DisplayName',['CH: ', num2str(i)],'LineWidth',1.5);
-    %     hold off
-    % end
-    
-    
-    
-    
-    % 
-    % a = movmean(sortrows([zdw_; S]'),10,'Endpoints','discard');
-    % 
-    % sorted = sortrows([zdw_; S]');
-    % figure(2)
-    % scatter(sorted(:,1),sorted(:,2))
-    % 
-    % ber_sorted = sort(S);
-    % mean(ber_sorted);
-    % std(ber_sorted);
-    % z1 = [];
-    % penalty = mean(ber_sorted):0.01:mean(ber_sorted)+2;
-    % for i = 1:length(penalty)
-    %     l = penalty(i);
-    %     if i == 1
-    %         z1 = [z1 sum(ber_sorted(1,:)<l ) / length(ber_sorted)  ];
-    %     else
-    %         z1 = [z1 sum(ber_sorted(1,:)<l & ber_sorted(1,:)>l-0.01) / length(ber_sorted)   ];
-    %     end
-    %     
-    % end
-    % 
-    % penalty_higherthan = 0.5;
-    % probability = sum(z1(find(penalty>mean(ber_sorted)+penalty_higherthan)));
-    % disp(['A penalty of more than 1dB has a probability of: ', num2str(probability)]);
-    % % 
-    % stem(AxesMain,penalty,z1,"filled",'Marker','o','MarkerSize',2,'Color',col);
-    % 
-    % 
-    % [f1,x1]=ecdf(S(end,:));
-    % %figure(23);plot(AxesMain,x1,f1,'r','LineWidth',3, 'Color',col);
-    % 
-    % %plot(AxesMain,a(:,1),a(:,2),'Color',col+1,'Parent',  AxesMain(1));
-    % 
-    % % histogram(AxesMain,S,1000,'EdgeColor','none','FaceAlpha',0.4);
-    % 
-    % 
-    % 
-    % % 
-    % %scatter(AxesMain,zdw_,S,'Marker','+','MarkerEdgeColor',col,'MarkerFaceAlpha',0.4,'MarkerEdgeAlpha',0.4,'LineWidth',0.5,'Parent',  AxesMain(1),'DisplayName',[plotJob.dataStatArg,' Ch.: ', plotJob.displayname]);
-    % hold on
-    % %scatter(AxesMain,zdw_chann(:,1),mean(S_chann,2),'Marker','diamond','MarkerEdgeColor',col,'MarkerFaceAlpha',0.6,'LineWidth',7,'MarkerFaceColor',col,'Parent',  AxesMain(1),'DisplayName',[plotJob.dataStatArg,' Ch.: ', plotJob.displayname]);
-    % hold off
-    % % 
-    % % for rlz = 1:size(S,2)
-    % %     zdwval      = zdw_(rlz);
-    % %     feccrossing = S(rlz);
-    % %     scatter(AxesMain,zdwval,feccrossing,10,'MarkerEdgeColor',col,'MarkerFaceColor',col,'Parent',  AxesMain(1),'DisplayName',[plotJob.dataStatArg,' Ch.: ', plotJob.displayname]);
-    % % end
-    % 
-    % xline([1302.03471732576,1304.30061112288,1306.57440520904,1308.85614097425,1311.14586009814,1313.44360455251,1315.74941660391,1318.06333881618]);
-
-
-end
-
-function vec = removeZeros(vec)
-    % Find rows that contain only zeros
-    rows_to_remove = all(vec == 0, 2);
-    
-    % Remove rows with only zeros
-    vec(rows_to_remove, :) = [];
-end

Classes/Warehouse_class/functions/fwm_plots/plotBerVsZdwFailureRate.m

@@ -1,157 +0,0 @@
-function plotBerVsZdwFailureRate(wh,plotJob)
-
-    fig = findall(groot, 'Type', 'figure', 'Name', plotJob.figName);
-    
-    if isvalid(fig)
-        figure(fig)
-        fig = get(fig);
-        AxesMain =  fig.CurrentAxes;
-        hold on
-    else
-        fig = figure('name',char(plotJob.figName));
-        AxesMain    = gca;
-        hold on
-    end
-    
-    
-    col = plotJob.color;
-    
-    % we want to fetch all realizations
-    realization = wh.parameter.realization.values(1:end);
-    
-    % get all xAxis values
-    xAxis = wh.parameter.p_out.values;
-    
-    % get all center wavelengths
-    wavelengths = wh.parameter.center_wavelength.values;
-    %wavelengths = wavelengths(2:end);
-%     totber = NaN(numel(realization),1,numel(xAxis),numel(wavelengths));
-%     ber = zeros(numel(realization),plotJob.ch,numel(xAxis),numel(wavelengths));
-    
-    %get BER values for query
-    for w = 1:numel(wavelengths)
-        
-        for xl = 1:numel(xAxis)
-    
-            c_wavelen = wavelengths(w);
-            p_out = xAxis(xl);
-    
-            % dim1 : realiz; dim2: channels, dim3: rop, dim4, c_wavelength
-            temp  = removeZeros(wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.ch,c_wavelen,plotJob.channelspacing,plotJob.randzdw));
-            ber(1:size(temp,1),:,xl,w)     = temp;
-            
-        end
-    end
-    
-    hdfec = 3.8e-3.*ones(size(xAxis));
-    zdw_ = [];
-    zdw_chann = [];
-    zdw_tot = [];
-    cf_tot = [];
-    cf_ = [];
-    S = [];
-    Stot = [];
-    S_chann = [];
-    cf_chann = [];
-    
-    cnt  = 0;
-    %get fec thresholds
-    %linear = squeeze(linear);
-    for c_wavelen = 1:size(ber,4)
-        for realiz = 1:size(ber,1)
-            for chann = 1:size(ber,2)
-            
-                %finde Schnittpunkt zwischen FEC und BER Kurve
-                temp_ber = squeeze(ber(realiz,chann,:,c_wavelen)).';
-                if ~all(temp_ber == 0)
-                   %nur wenn nicht alles nullen sind
-                    crossing_ch = InterX([hdfec;xAxis],[temp_ber;xAxis]);
-                else
-                    continue
-                end
-
-                %Req. FEC Ergebnis einsortieren
-                if ~isempty(crossing_ch)
-                    if crossing_ch(2) == 0
-                        print("d")
-                    end
-                    S(realiz,chann,c_wavelen) = crossing_ch(2);
-                else
-                    S(realiz,chann,c_wavelen) = -1;
-                    cnt = cnt +1;
-                end
-        
-        
-            end
-        end
-    end
-    
-    temp_max = -inf;
-    sum_FEC_not_crossed=[];
-    sum_FEC_crossed=[];
-    
-    threshold = plotJob.p_in - 10;
-
-    for i = 1:plotJob.ch
-        hold on
-        %S:: 1.dim: realiz; 2.dim: channel; 3.dim: center wavelength
-        %squeeze a channel:
-        temp_data = squeeze(S(:,i,:));
-
-        %remove realizations that have no entry (only zero)
-        temp_data = removeZeros(temp_data);
-
-        %replace zeros with NAN (e.g. for the wavelengths that have missing realizations)
-        temp_data(temp_data==0) = NaN;
-        
-        %for current channel
-        FEC_crossed     = temp_data < threshold & ~isnan(temp_data);
-        FEC_not_crossed = temp_data >= threshold & ~isnan(temp_data);
-
-        %sum over channels for overall picture
-        sum_FEC_not_crossed(i,:) =  sum(FEC_not_crossed);
-        sum_FEC_crossed(i,:)     =  sum(FEC_crossed);
-
-        failure_rate_channelwise(i,:)     = sum_FEC_not_crossed(i,:)./ ( sum_FEC_crossed(i,:) + sum_FEC_not_crossed(i,:));
-
-    end
-
-    failure_rate_total = sum(sum_FEC_not_crossed,1) ./ ( sum(sum_FEC_crossed,1) + sum(sum_FEC_not_crossed,1) );
-    % plot failure rate (nbetween 0 and 1)
-
-    plot(wavelengths,failure_rate_total,'Color',plotJob.color,'LineWidth',1,'LineStyle',plotJob.plotTypeArg,'Marker','x','MarkerSize',5,'MarkerFaceColor',plotJob.color,'DisplayName',plotJob.displayname);
-
-    %plot max overall
-    %scatter(wavelengths,failure_rate_channelwise ,35,plotJob.color,'Marker','.');
-
-    ylabel('Failure Rate of Link');
-    xlabel('Wavelength in nm');
-
-    xlim([min(wavelengths),max(wavelengths) ]);
-    ylim([0,1]);
-
-    grid minor;
-    set(gca, 'color', 'none');
-    legend = [];
-    
-%     fontsize(AxesMain,8,"points")
-    
-    fig.Position = plotJob.Position;
-    fig.Units = "centimeters";
-    fig.Position = [0 0 12 5 7];
-    
-    try
-    set(AxesMain,'TickLabelInterpreter','latex')
-    
-    set(AxesMain.Legend,'Interpreter','latex')
-    end
-    
-end
-
-function vec = removeZeros(vec)
-    % Find rows that contain only zeros
-    rows_to_remove = all(vec == 0, 2);
-    
-    % Remove rows with only zeros
-    vec(rows_to_remove, :) = [];
-end

Classes/Warehouse_class/functions/fwm_plots/plotChannelSpacingAna.m

@@ -1,51 +0,0 @@
-function plotChannelSpacingAna(wh,plotJob)
-
-xAxis = wh.parameter.p_out.values;
-
-
-realization = wh.parameter.realization.values(1:end);
-
-channelsp = wh.parameter.channelspacing.values(1:end);
-channelsp = [200 400].*1e9;
-for ch = 1:2
-    channspacing = channelsp(ch);
-    for xl = 1:numel(xAxis)
-        p_out = xAxis(xl);
-    
-        curber = wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.numchannels,channspacing);
-        curzdw = wh.getStoValue('zdw',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.numchannels,channspacing);
-    
-        ber(1:size(curber,1),1:size(curber,2),xl) = curber;
-        zdw(1:size(curber,1),1,xl)  = curzdw;
-    
-    end
-    
-    ber = squeeze(mean(ber,1));
-    zdw = squeeze(mean(zdw,1));
-    
-    hdfec = 3.8e-3.*ones(size(xAxis));
-    S = [];
-    wavelength={};
-    
-    zdw_ = [];
-    zdw_chann = [];
-    zdw_tot = [];
-    S = [];
-    S_chann = [];
-    wl = round([1302.03471732576,1304.30061112288,1306.57440520904,1308.85614097425,1311.14586009814,1313.44360455251,1315.74941660391,1318.06333881618],2);
-    wl = 1:16;
-    wl = [1.2930    1.2953    1.2975    1.2998    1.3020    1.3043    1.3066    1.3089    1.3111    1.3134    1.3157    1.3181    1.3204    1.3227    1.3251    1.3274];
-    
-    
-    a = InterX([hdfec(:)';xAxis],[mean(ber,1);xAxis]);
-    if ~isempty(a)
-        s(ch) = a(2);
-    else
-        s(ch) = NaN;
-    end
-end
-
-figure(2224)
-hold on
-plot(channelsp,s,'LineWidth',1,'Color',plotJob.color,'Marker','o');
-

Classes/Warehouse_class/functions/fwm_plots/plotCurve.m

@@ -1,278 +0,0 @@
-function plotCurve(wh,plotJob)
-%PLOTCURVE Summary of this function goes here
-%   Detailed explanation goes here
-
-fig = findall(groot, 'Type', 'figure', 'Name', plotJob.figName);
-
-if isvalid(fig)
-    figure(fig)
-    fig = get(fig);
-    AxesMain =  fig.CurrentAxes;
-    hold on
-else
-    fig = figure('name',char(plotJob.figName));
-    AxesMain    = gca;
-    hold on
-end
-
-col = plotJob.color;
-
-% we want to fetch all realizations
-if plotJob.pmd == 0
-    realization = 499;
-    % realization = 0:7;
-else
-    realization = wh.parameter.realization.values(1:end);
-end
-realization = wh.parameter.realization.values(1:end);
-% get all xAxis values
-xAxis = wh.parameter.p_out.values;
-
-% Fetch Data from Warehouse
-for xl = 1:numel(xAxis)
-    p_out = xAxis(xl);
-    if string(plotJob.dataStatArg) == "Worst"
-        temp = wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.ch,1310,plotJob.channelspacing,plotJob.randzdw);
-        temp = removeZeros(temp);
-        ber(xl) = max(temp,[],'all');
-        linew = 1.0;
-        markersz = 3;
-        linestyle = '-';
-        
-    elseif string(plotJob.dataStatArg) == "AVG"
-        temp = wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.ch,1310,plotJob.channelspacing,plotJob.randzdw);
-        temp = removeZeros(temp);
-        ber(xl) = mean(temp,'all');
-        linew = 1.0;
-        markersz = 3;
-        linestyle = '--';
-
-    elseif string(plotJob.dataStatArg) == "Best"
-        temp = wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.ch,1310,plotJob.channelspacing,plotJob.randzdw);
-        temp = removeZeros(temp);
-        ber(xl) = min(temp,[],'all');
-        linew = 1.0;
-        markersz = 3;
-        linestyle = '-';
-
-    elseif string(plotJob.dataStatArg) == "All Channels; mean(PMD Realizations)"
-        temp = wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.ch,1310,plotJob.channelspacing,plotJob.randzdw);
-        temp = removeZeros(temp);
-
-        ber(:,xl) = mean(temp,1,"omitnan").';
-
-        linew = 1;
-        markersz = 1;
-        linestyle = '-';
-
-    elseif string(plotJob.dataStatArg) == "Lineplot with quartiles"
-
-        temp = wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.ch,1310,plotJob.channelspacing,plotJob.randzdw);
-        temp = removeZeros(temp);
-        ber(xl) = mean(temp,"all","omitnan").';
-
-
-        upperq(xl) = quantile(temp,0.99,"all");
-        lowerq(xl) = quantile(temp,0.04,"all");
-        
-        if lowerq(xl) == 0
-            lowerq(xl) = lowerq(xl-1);
-        end
-        %         upperq(xl) =  0.5*std(tmp,1,'all','omitnan');
-        %         lowerq(xl) =  0.5*std(tmp,1,'all','omitnan');
-
-        %         upperq(xl) = max(dataNoNans);
-        %         lowerq(xl) = min(dataNoNans);
-
-        %                 upperq(xl) = mean(tmp,"all","omitnan") + 1.96 * (std(tmp,1,'all','omitnan')/sqrt(numel(tmp)));
-        %                 lowerq(xl) = mean(tmp,"all","omitnan") - 1.96 * (std(tmp,1,'all','omitnan')/sqrt(numel(tmp)));
-
-        linew = 1;
-        markersz = 1;
-        linestyle = '-';
-
-    elseif string(plotJob.dataStatArg) ==  "All Channels ;All PMD Realizations"
-
-        raw_fetch = wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.ch,1310,plotJob.channelspacing,plotJob.randzdw).';
-        tmp = reshape(raw_fetch,[],1);
-        ber(1:size(tmp,1),xl) = tmp;
-
-        ber_per_chann(:,:,xl) = raw_fetch;
-
-        linew = 0.3;
-        markersz = 2;
-        linestyle = ':';
-        
-    end
-end
-
-
-%routine to remove total outliers (here those wehere the rop curve has a mean BER greater than 0.1)
-% ber_per_chann_clean = NaN(size(ber_per_chann));
-% for ch = 1:size(ber_per_chann,1)
-%     bla = squeeze(ber_per_chann(ch,:,:));
-%     ber_per_chann(ch,find(mean(bla,2)>0.25),:) = NaN;
-%     cleaned = rmoutliers(bla,"mean",'ThresholdFactor',2);
-%     
-%     ber_per_chann_clean(ch,1:size(cleaned,1),1:size(cleaned,2)) = cleaned;
-%    
-% end
-% 
-% ber = [];
-% for rop = 1:size(ber_per_chann,3)
-%     temp = squeeze(ber_per_chann_clean(:,:,rop));
-%     ber(rop) = mean(temp,"all","omitnan").';
-%     upperq(rop) = quantile(temp,0.9,"all");
-%     lowerq(rop) = quantile(temp,0.1,"all");
-% end
-
-
-
-
-xAxis = xAxis;
-
-% Plot Data
-if string(plotJob.plotTypeArg) == "Scatter"
-    for rlz = 1:size(ber,1)
-
-        if rlz < size(ber,1)
-            scatter(xAxis,ber(rlz,:),markersz,'MarkerEdgeColor',col,'MarkerFaceColor',col,'Parent',  AxesMain(1),'HandleVisibility','off');
-        else
-            scatter(xAxis,ber(rlz,:),markersz,'MarkerEdgeColor',col,'MarkerFaceColor',col,'Parent',  AxesMain(1),'DisplayName',[plotJob.dataStatArg,' Ch.: ', plotJob.displayname]);
-        end
-
-    end
-
-elseif string(plotJob.plotTypeArg) == "Lines"
-
-    %         if ~anynan(ber)
-    %             [xAxis,ber] = interpCurve(xAxis, ber);
-    %         end
-
-    for rlz = 1:size(ber,1)
-% 
-        if (string(plotJob.dataStatArg) == "All Channels ;All PMD Realizations")||(string(plotJob.dataStatArg) == "All Channels; mean(PMD Realizations)")
-            ch = mod(rlz,plotJob.ch);
-            if ch == 0; ch = plotJob.ch; end
-        else
-            ch = plotJob.dataStatArg;
-        end
-
-        if rlz < size(ber,1)
-            
-
-            s = plot(xAxis,ber(rlz,:),linestyle,'Marker',"none",'MarkerSize',markersz,'LineWidth',linew,'Color',col,'Parent',  AxesMain,'HandleVisibility','off');
-
-            s.DataTipTemplate.Interpreter = "latex";
-            s.DataTipTemplate.DataTipRows(1).Label = "Ch: ";
-            s.DataTipTemplate.DataTipRows(1).Value = repmat(ch,size(ber));
-            s.DataTipTemplate.DataTipRows(1);
-            s.DataTipTemplate.DataTipRows(2) = [];
-
-
-        else
-
-            if string(plotJob.dataStatArg) == "Lineplot with quartiles"
-                [hl,hp] = boundedline(xAxis,ber(rlz,:),([(ber(rlz,:)-lowerq(rlz,:));upperq(rlz,:)-ber(rlz,:)]'),'-o','alpha','Color',col,'transparency', 0.1,'linewidth',0.7);
-                hl.MarkerFaceColor = col;
-                hl.MarkerSize = 3;
-                set(hp,'HandleVisibility','off');
-                %hp.LineWidth = 1.2;
-
-                ho = outlinebounds(hl,hp);
-                set(ho, 'linestyle', ':', 'color', col,'Linewidth',0.5);
-                set(ho,'HandleVisibility','off');
-            else
-
-                s = plot(xAxis,ber(rlz,:),linestyle,'Marker',"o",'MarkerFaceColor',col,'MarkerSize',markersz,'LineWidth',linew,'Color',col,'Parent',  AxesMain,'DisplayName',[plotJob.displayname]);
-                
-                s.DataTipTemplate.Interpreter = "latex";
-                s.DataTipTemplate.DataTipRows(1).Label = "Ch: ";
-                s.DataTipTemplate.DataTipRows(1).Value = repmat(string(ch),size(ber));
-                s.DataTipTemplate.DataTipRows(1)
-                s.DataTipTemplate.DataTipRows(2) = [];
-            end
-
-        end
-
-    end
-end
-
-% Draw FEC Threshold Line
-%get x data of first children:
-%get all linear Values
-if 0
-    linear   = wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,wh.parameter.p_out.values,0,0,499,"symmetric");
-    linear   = mean(linear,2);
-    lincurve = findall(AxesMain, 'Type', 'line','DisplayName','Linear Baseline');
-    %
-    if isempty(lincurve)
-        xdata = AxesMain.Children(1).XData;
-        hdfec = 3.8e-3.*ones(size(xdata));
-        plot(xdata,linear,'-','Marker',"o",'MarkerSize',3,'LineWidth',4,'Color',[0.6400    0.6400    0.6400],'MarkerFaceColor',[0.6400    0.6400    0.6400],'Parent', AxesMain,'DisplayName','Linear Baseline');
-        %h = get(gca,'Children');
-        %set(gca,'Children',[h(2) h(1)])
-    end
-end
-
-feccurve = findall(AxesMain, 'Type', 'line','DisplayName','FEC $3.8*10^{-3}$');
-%
-if isempty(feccurve)
-    xdata = AxesMain.Children(1).XData;
-    hdfec = 3.8e-3.*ones(size(xdata));
-    plot(xdata,hdfec,'--','MarkerSize',4,'Color','black','MarkerFaceColor','black','LineWidth',1,'Parent', AxesMain,'DisplayName','FEC $3.8*10^{-3}$','HandleVisibility','off');
-    %h = get(gca,'Children');
-    %set(gca,'Children',[h(2) h(1)])
-end
-
-% Figure Settings
-%title(AxesMain,plotJob.title,"Interpreter","none");
-
-xlabel(AxesMain,plotJob.xAxisLabel,"Interpreter","latex");
-
-ylabel(AxesMain,plotJob.yAxisLabel,"Interpreter","latex");
-
-set(AxesMain,'yscale','log');
-
-grid(AxesMain,'on');
-
-grid(AxesMain,'minor');
-
-grid minor
-
-%legend(AxesMain);
-
-fontsize(AxesMain,8,"points")
-% fontname(AxesMain,"Arial")
-
-fig.Position = plotJob.Position;
-fig.Units = "centimeters";
-fig.Position = [2 2 8.5 7];
-
-set(AxesMain,'TickLabelInterpreter','latex')
-
-set(AxesMain.Legend,'Interpreter','latex')
-% set(gcf,'Units','centimeters')
-% set(gcf,'Position',[2 2 9 4.5])
-
-ylim([1e-5,0.3]);
-
-xlim([min(xAxis),-3]);
-
-
-% annotation('textbox', [0.125, 0.32, 0.1, 0.1], 'String', "FEC 3.8e-3","LineStyle","none","FontSize",8,"FontUnits","points")
-
-
-hold off
-
-
-end
-
-
-function vec = removeZeros(vec)
-    % Find rows that contain only zeros
-    rows_to_remove = all(vec == 0, 2);
-    
-    % Remove rows with only zeros
-    vec(rows_to_remove, :) = [];
-end

Classes/Warehouse_class/functions/fwm_plots/plotHistogram.m

@@ -1,151 +0,0 @@
-function plotHistogram(wh,plotJob)
-
-
-fig = findall(groot, 'Type', 'figure', 'Name', plotJob.figName);
-
-if isvalid(fig)
-    figure(fig)
-    fig = get(fig);
-    AxesMain =  fig.CurrentAxes;
-    hold on
-else
-    fig = figure('name',char(plotJob.figName));
-    AxesMain    = gca;
-    hold on
-end
-
-col = plotJob.color;
-
-% we want to fetch all realizations
-realization = wh.parameter.realization.values(1:end);
-
-% get all xAxis values
-xAxis = wh.parameter.p_out.values;
-
-
-% Fetch Data from Warehouse
-for xl = 1:numel(xAxis)
-    p_out = xAxis(xl);
-    if string(plotJob.dataStatArg) == "Worst"
-        ber(xl) = max(wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,string(plotJob.channelplan),string(plotJob.dsp)),[],'all');
-        linew = 2.0;
-        markersz = 3;
-        linestyle = '-';
-    elseif string(plotJob.dataStatArg) == "AVG"
-        ber(xl) = mean(wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,string(plotJob.channelplan),string(plotJob.dsp)),'all');
-        linew = 2.0;
-        markersz = 3;
-        linestyle = ':';
-    elseif string(plotJob.dataStatArg) == "Best"
-        ber(xl) = min(wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,string(plotJob.channelplan),string(plotJob.dsp)),[],'all');
-        linew = 2.0;
-        markersz = 3;
-        linestyle = ':';
-
-    elseif string(plotJob.dataStatArg) == "All Channels; mean(PMD Realizations)"
-        tmp = wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,string(plotJob.channelplan),string(plotJob.dsp));
-        
-        if numel(tmp(tmp==0)) ~= 0
-            disp('Removed all zero values!');
-            tmp(tmp==0) = NaN;
-        end
-        ber(:,xl) = mean(tmp,1,"omitnan").';
-
-        linew = 0.7;
-        markersz = 2;
-        linestyle = '-';
-
-    elseif string(plotJob.dataStatArg) == "All Channels ;All PMD Realizations"
-
-
-        
-        tmp = reshape(wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,string(plotJob.channelplan),string(plotJob.dsp)).',[],1);
-        ber(1:size(tmp,1),xl) = tmp;
-
-        linew = 0.3;
-        markersz = 2;
-        linestyle = ':';
-    end
-end
-
-disp('Removed all zero values!');
-ber(ber==0) = NaN;
-if ~anynan(ber)
-    [xAxis,ber] = interpCurve(xAxis, ber);
-end
-
-% plot FEC Crossing as histogram
-
-hdfec = 3.8e-3.*ones(size(xAxis));
-S = [];
-for i = 1:size(ber,1)
-    a = InterX([hdfec;xAxis],[ber(i,:);xAxis]);
-    if ~isempty(a)
-        S(:,i) = a;
-    end
-end
-
-
-%% SUB 1
-AxesMain = subplot(2,1,1);
-
-hold on
-
-if ~isempty(S)
-    histogram(S(end,:),300,'Normalization','probability','FaceColor',col,'EdgeColor',col,'Parent',AxesMain,'DisplayName',[plotJob.dataStatArg,' Ch.: ', plotJob.displayname],'FaceAlpha',0.4,'EdgeAlpha',0.4);
-end
-
-xlim([-3 ,9 ]);
-ylim([0 .10]);
-
-% Figure Settings
-title(AxesMain,['$P_{in}:$ ',num2str(plotJob.p_in-9), 'dBm; L : ', num2str(plotJob.l),'Km'],'Interpreter','latex');
-
-xlabel(AxesMain,plotJob.xAxisLabel,'Interpreter','latex');
-
-ylabel('PDF')
-
-grid(AxesMain,'on');
-
-grid(AxesMain,'minor');
-
-%legend(AxesMain,'Interpreter','latex');
-
-fontsize(AxesMain,24,"pixels")
-
-hold off
-
-
-%% SUB 2
-AxesMain = subplot(2,1,2);
-
-if ~isempty(S)
-hold on
-
-[f1,x1]=ecdf(S(end,:));
-    plot(x1,f1,'r','LineWidth',3, 'Color',col,'DisplayName',[plotJob.dataStatArg,' Ch.: ', plotJob.displayname]);
-end
-xlim([-3 ,9 ]);
-ylim([0 1]);
-% Figure Settings
-title(AxesMain,['$P_{in}:$ ',num2str(plotJob.p_in-9), 'dBm; L : ', num2str(plotJob.l),'Km'],'Interpreter','latex');
-
-xlabel(AxesMain,plotJob.xAxisLabel,'Interpreter','latex');
-
-ylabel('CDF')
-
-grid(AxesMain,'on');
-
-grid(AxesMain,'minor');
-
-fontsize(AxesMain,24,"pixels")
-
-%legend(AxesMain,'Interpreter','latex');
-
-fig.Position = plotJob.Position;
-
-hold off
-
-
-
-end

Classes/Warehouse_class/functions/fwm_plots/plotViolin.m

@@ -1,206 +0,0 @@
-function plotViolin(wh,plotJob)
-
-fig = findall(groot, 'Type', 'figure', 'Name', plotJob.figName);
-
-if isvalid(fig)
-    figure(fig)
-    fig = get(fig);
-    AxesMain =  fig.CurrentAxes;
-    hold on
-else
-    fig = figure('name',char(plotJob.figName));
-    AxesMain    = gca;
-    hold on
-end
-
-%% Violin
-col = plotJob.color;
-
-% we want to fetch all realizations
-if plotJob.pmd == 0
-    realization = 1;
-else
-    realization = wh.parameter.realization.values(1:end);
-end
-
-%realization = 0:8;
-
-% get all xAxis values
-xAxis = wh.parameter.p_out.values;
-
-% ber = NaN(500,16,10);
-% zdw = NaN(500,1,10);
-
-%get BER values for query
-for xl = 1:numel(xAxis)
-    p_out = xAxis(xl);
-
-    curber = wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.ch,1310,plotJob.channelspacing,plotJob.randzdw);
-    curber= removeZeros(curber);
-    ber(1:size(curber,1),1:size(curber,2),xl) = curber;
-
-end
-
-% to remove outliers set the percentile range
-% a = squeeze(mean(ber,2));
-% out = isoutlier(mean(a,2),"percentiles",[0 100]);
-% ber = ber(~out,:,:);
-% disp(sum(out));
-
-hdfec = 3.8e-3.*ones(size(xAxis));
-S = [];
-wavelength={};
-
-
-S = [];
-S_chann = [];
-
-wl = calcWavelengthPlan(plotJob.ch,plotJob.channelspacing,1310);
-%get fec thresholds
-% [C,ia,ib] =intersect(linx,xAxis);
-% linear = squeeze(linear);
-
-%ber(ber==0) = NaN;
-S_chann_no_crossing = zeros(1,plotJob.ch);
-for chann = 1:size(ber,2)
-
-    for realiz = 1:size(ber,1)
-        
-        ber_series = squeeze(ber(realiz,chann,:)).';
-        if mean(ber_series) > 0.1
-            continue
-        end
-
-        a = InterX([hdfec(:)';xAxis],[ber_series;xAxis]);
-
-        if ~isempty(a)
-            S_chann(realiz,chann) = a(2);
-%             if a(2) > -7 && string(plotJob.pol) == "copolarized"
-%                 continue
-%             end
-            S(end+1) = a(2);
-            wavelength{end+1} = num2str(wl(chann));
-        else
-            S(end+1) = 0;
-            wavelength{end+1} = num2str(wl(chann));
-            S_chann_no_crossing(realiz,chann) = 1;
-            S_chann(realiz,chann) = -1;
-        end
-    
-    end
-end
-
-threshold = -6;
-FEC_crossed     = sum(S_chann < threshold & ~isnan(S_chann),1);
-FEC_not_crossed = sum(S_chann >= threshold & ~isnan(S_chann),1);
-failure_rate = FEC_not_crossed ./ (FEC_crossed + FEC_not_crossed) ;
-
-
-S_chann(S_chann==0) = NaN;
-
-total_avg = mean(S_chann,"all","omitnan");
-
-%figure(2024)
-%C = flip(cbrewer2('Spectral',8));
-if numel(S) <= numel(wl)
-    vs = scatter(1:numel(S),S,50,'o','MarkerEdgeColor','black','MarkerFaceColor',plotJob.color,'LineWidth',1,'HandleVisibility','off');
-    %vs = scatter(1,mean(S),50,'o','MarkerEdgeColor','black','MarkerFaceColor',plotJob.color,'LineWidth',1);
-
-else
-
-    vs = violinplot(S,wavelength,...
-        'ViolinColor',plotJob.color,...
-        'ViolinAlpha',0.1,...
-        'MarkerSize',1,...
-        'ShowMedian',false,...
-        'EdgeColor',plotJob.color,...
-        'ShowWhiskers',false,...
-        'ShowData',false,...
-        'ShowBox',false,...
-        'Bandwidth',0.051 ...
-        );
-
-    
-    hold on
-
-    partly_failed = boolean(ceil(failure_rate));
-    
-    avg = mean(S_chann,1,"omitnan");
-    
-    notfailed = ~partly_failed .* avg;
-    notfailed(notfailed==0) = NaN;
-    scatter(1:size(S_chann,2),notfailed,10,'Marker','x','MarkerEdgeColor','black','MarkerFaceColor',plotJob.color,'LineWidth',0.5,'HandleVisibility','off');
-
-
-    hold on
-    partly_failed = partly_failed.*avg;
-    partly_failed(partly_failed==0) = NaN;
-    s=scatter(1:numel(failure_rate),partly_failed,10,'Marker','x','LineWidth',0.5,'HandleVisibility','off','MarkerEdgeColor','red');
-    s.DataTipTemplate.Interpreter = "latex";
-    s.DataTipTemplate.DataTipRows(1).Label = "Fail Rate: ";
-    s.DataTipTemplate.DataTipRows(1).Value = failure_rate;
-    s.DataTipTemplate.DataTipRows(2) = [];
-    hold off
-end
-
-% ax = gca;
-% ax.XTicks
-
-hold on
-yline(total_avg,'LineWidth',1,'LineStyle','--','DisplayName','System Avg.')
-fig.Position = plotJob.Position;
-
-xticklabels(1:16);
-ylabel('Penalty in dB');
-xlabel('Channel Number');
-ylim([-9.3,-3]);
-xlim([0,plotJob.ch+1]);
-grid minor;
-set(gca, 'color', 'none');
-legend = [];
-
-fontsize(AxesMain,8,"points")
-
-fig.Position = plotJob.Position;
-fig.Units = "centimeters";
-fig.Position = [2 2 8.5 7];
-
-set(AxesMain,'TickLabelInterpreter','latex')
-
-set(AxesMain.Legend,'Interpreter','latex')
-
-
-
-if 0
-ber_sorted = sort(S);
-
-z1 = [];
-penalty = mean(ber_sorted):0.01:mean(ber_sorted)+2;
-
-for i = 1:length(penalty)
-    l = penalty(i);
-    if i == 1
-        z1 = [z1 sum(ber_sorted(1,:)<l ) / length(ber_sorted)  ];
-    else
-        z1 = [z1 sum(ber_sorted(1,:)<l & ber_sorted(1,:)>l-0.01) / length(ber_sorted)   ];
-    end
-    
-end
-
-penalty_higherthan = 0.5;
-probability = sum(z1(find(penalty>mean(ber_sorted)+penalty_higherthan)));
-disp(['A penalty of more than 0.5 dB has a probability of: ', num2str(probability)]);
-end
-
-end
-
-
-function vec = removeZeros(vec)
-    % Find rows that contain only zeros
-    rows_to_remove = all(vec == 0, 2);
-    
-    % Remove rows with only zeros
-    vec(rows_to_remove, :) = [];
-
-
-end

Classes/Warehouse_class/functions/fwm_plots/plot_ber_distribution.m

@@ -1,58 +0,0 @@
-
-%automate plots
-% [file, path] = uigetfile("C:\Users\Silas\Documents\MATLAB\Raw_Cluster_Simulations\session_januar24\wh_complete_at_1310.mat");
-% wh = load([path filesep file]);
-% wh = wh.wh;
-
-plotJob = struct();
-
-plotJob.l = 10; 
-plotJob.ch = 16; 
-plotJob.d = 3; 
-plotJob.sgm = 1; 
-plotJob.pol = "copolarized"; 
-plotJob.p_in = 3; 
-plotJob.gamma = 0.0023;
-plotJob.pmd = 0.1; 
-plotJob.channelspacing = 400e9; 
-plotJob.randzdw = 0; 
-
-ber_per_chann = [];
-% get all xAxis values
-xAxis = wh.parameter.p_out.values;
-realization = wh.parameter.realization.values(1:end);
-% Fetch Data from Warehouse
-for xl = 1:numel(xAxis)
-    p_out = xAxis(xl);
-    raw_fetch = wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.ch,1310,plotJob.channelspacing,plotJob.randzdw).';
-       
-    ber_per_chann(:,:,xl) = raw_fetch;
-end
-%%
-
-figure(2023);
-
-for ch = [1,floor(plotJob.ch/2),ceil(plotJob.ch/2)+1,plotJob.ch] %1:15:size(ber_per_chann,1)
-    
-    for p = 5%1:size(ber_per_chann,3)
-
-        % Extract data for the current row
-        row_data = squeeze(ber_per_chann(ch,:,p));
-        [f, xi] = ksdensity(row_data);
-          % Identify the peak point
-        [max_density, max_index] = max(f);
-        peak_x = xi(max_index);
-
-    end
-    % Create a histogram plot for the current row with a unique color
-    plot(xi, f, 'LineWidth', 2, 'DisplayName', ['Ch. ', num2str(ch)],'LineStyle','--');
-    %histogram(row_data,100, 'DisplayName', ['Channel ', num2str(ch)],  'EdgeColor', 'none');
-    
-    hold on; % Hold the plot for the next iteration
-    text(peak_x, max_density, ['Ch ', num2str(ch)], 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'left');
-end
-
-
-legend show
-
-%%