more and more

Classes/04_DSP/Equalizer/FFE.m

@@ -313,7 +313,7 @@ classdef FFE < handle
             end
             obj.mu_optimization_iter = 0;
             obj.mu_optimization = bayesopt(@(p)obj.muObjective(p,x,d),vars, ...
-                "MaxObjectiveEvaluations",10, ...
+                "MaxObjectiveEvaluations",20, ...
                 "AcquisitionFunctionName","expected-improvement-plus", ...
                 "IsObjectiveDeterministic",false, ...
                 "Verbose",0, ...
@@ -341,7 +341,7 @@ classdef FFE < handle
             end
         end
 
-        function objective = muObjective(obj,params,x,d)
+        function objective_db = muObjective(obj,params,x,d)
             old_debug = obj.save_debug;
             old_mu_dc = obj.mu_dc;
             obj.save_debug = 1;

Functions/EQ_blocks/preprocessSignal.m

@@ -33,7 +33,7 @@ switch preprocessMode
         error('preprocessSignal:InvalidMode', 'Unsupported preprocessing mode "%s".', preprocessMode);
 end
 
-[Scpe_sig, Scpe_cell] = Scpe_sig.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", options.debug_plots);
+[Scpe_sig, Scpe_cell, inverted] = Scpe_sig.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", options.debug_plots);
  % Scpe_sig = Scpe_cell{1};
 
 % Apply Gaussian filter

Functions/EQ_recipes/mpi_recipe_dev.m

@@ -21,18 +21,14 @@ function output = mpi_recipe_dev(Scpe_sig_raw, Symbols, Tx_bits, options)
         "debug_plots", options.debug_plots);
     
     
-    mu_dc = 0; % 1e-5
+    mu_dc = 0.005;%1e-5;
-    dc_buffer_len = 0;
-    ffe_buffer_len = 0;
-    smoothing_buffer_length = options.userParameters.smoothing_length;
-    smoothing_buffer_update = 1;
 
     eq_settings = { ...
         "epochs_tr", 5, ...
-        "epochs_dd", 5, ...
+        "epochs_dd", 3, ...
         "len_tr", 4096*2, ...
-        "mu_dd", 1e-1, ...
+        "mu_dd",3.005e-03, ...
-        "mu_tr", 0.4, ...
+        "mu_tr",1.063e-01, ...
         "order", 25, ...
         "sps", 2, ...
         "decide", 0, ...
@@ -43,36 +39,54 @@ function output = mpi_recipe_dev(Scpe_sig_raw, Symbols, Tx_bits, options)
 
     eq_ffe = FFE(eq_settings{:});
     
-    % showLevelScatter(Scpe_sig_raw, Symbols, ...
+    showLevelScatter(Scpe_sig_raw, Symbols, ...
-    %     "fsym", options.fsym, ...
+        "fsym", options.fsym, ...
-    %     "fignum", options.dataTable.run_id, ...
+        "fignum", 400, ...
-    %     "normalize", true);
+        "normalize", true);
 
-    %%
+    % Scpe_sig_raw.spectrum("normalizeTo0dB",1,"fft_length",4096*4,"fignum",401);
-    % tic
-    % ffe_results = ffe(eq_ffe, options.M, Scpe_sig, Symbols, Tx_bits, ...
-    %     "precode_mode", options.duob_mode, ...
-    %     'showAnalysis', options.debug_plots, ...
-    %     "postFFE", [], ...
-    %     "eth_style_symbol_mapping", 0);
-    % toc
-    % ffe_results.metrics.print("description",'FFE');
-    % output.ffe_package = ffe_results;
 
-    %%
+    %options.dataTable.sir;
-    eq_ffe_dcr = FFE_DCremoval_adaptive_mu(eq_settings{:}, ...
+
-        "dc_buffer_len",dc_buffer_len, ...
+    %% NORMAL FFE
-        "ffe_buffer_len",ffe_buffer_len,...
+    if  1
-        "smoothing_buffer_length",smoothing_buffer_length,...
+
-        "smoothing_buffer_update",smoothing_buffer_update);
+        ffe_results = ffe(eq_ffe, options.M, Scpe_sig, Symbols, Tx_bits, ...
-    tic
-    ffe_results_dcr = ffe(eq_ffe_dcr, options.M, Scpe_sig, Symbols, Tx_bits, ...
             "precode_mode", options.duob_mode, ...
             'showAnalysis', options.debug_plots, ...
             "postFFE", [], ...
             "eth_style_symbol_mapping", 0);
-    toc
+        
+        ffe_results.metrics.print("description",sprintf('Normal FFE; SIR %d dB',options.dataTable.sir));
+        output.ffe_package = ffe_results;
+
+    end
+
+    %% MPI Reduction 
+    if 1
+
+        % dc_buffer_len = 0;
+        % ffe_buffer_len = 0;
+        % smoothing_buffer_length = options.userParameters.smoothing_length;
+        % smoothing_buffer_update = 1;
+        
+        % eq_ffe_dcr = FFE_DCremoval_adaptive_mu(eq_settings{:}, ...
+        %     "dc_buffer_len",dc_buffer_len, ...
+        %     "ffe_buffer_len",ffe_buffer_len,...
+        %     "smoothing_buffer_length",smoothing_buffer_length,...
+        %     "smoothing_buffer_update",smoothing_buffer_update);
+
+        eq_ = FFE_DCremoval("epochs_tr",5,"epochs_dd",3,"len_tr",4096*2,"mu_dd",...
+            0.0002,"mu_tr",0,"order",25,"sps",2,"decide",0,...
+            "mu_dc",0.005,"dc_buffer_len",1);
+        
+        ffe_results_dcr = ffe(eq_, options.M, Scpe_sig, Symbols, Tx_bits, ...
+            "precode_mode", options.duob_mode, ...
+            'showAnalysis', options.debug_plots, ...
+            "postFFE", [], ...
+            "eth_style_symbol_mapping", 0);
+        
         ffe_results_dcr.metrics.print("description",'FFE DCR');
         output.ffe_dcr_package = ffe_results_dcr;
-
+    end
 end

Functions/EQ_visuals/showEQNoisePSD.m

@@ -23,7 +23,8 @@ end
     end
 
     % Ensure the figure is ready before calling spectrum
-    eq_noise.spectrum("displayname", options.displayname, "fignum", fig.Number, "normalizeTo0dB", 1,"color",options.color);
+    eq_noise = eq_noise - mean(eq_noise.signal);
+    eq_noise.spectrum("displayname", options.displayname, "fignum", fig.Number, "normalizeTo0dB", 0,"color",options.color,"fft_length",4096*2);
 
     if ~isnan(options.postfilter_taps)
         % Hold on to the figure for further plotting

Functions/EQ_visuals/showLevelHistogram.m

@@ -38,6 +38,8 @@ end
 
     eq_signal = max(min(eq_signal,3),-3);
 
+    eq_signal = eq_signal .* -1;
+
     %%% histogram 
     clf
     if isempty(ref_symbol_uncoded)
@@ -47,6 +49,15 @@ end
         received_sd = NaN(numel(constellation),numel(ref_symbols));
         lvlcol = linspecer(numel(constellation));
 
+        lvlcol = cbrewer2('RdBu',numel(constellation)+4);
+
+        % remove 4 colors from the middle to avoid the bright central colors
+        mid = ceil(size(lvlcol,1)/2);
+        rm_idx = mid + (-1:2); % two before mid and two after (2x2 removal centered)
+        rm_idx = max(1,min(size(lvlcol,1),rm_idx));
+        lvlcol(rm_idx,:) = [];
+        lvlcol = lvlcol(1:numel(constellation),:);
+
         for lvl = 1:numel(constellation)
             class_mask = ref_symbols == constellation(lvl);
             received_sd(lvl,class_mask) = eq_signal(class_mask);
@@ -71,7 +82,15 @@ end
         % DB lobe is drawn separately so the DB-level structure stays visible.
         db_constellation = unique(ref_symbols);
         classes = unique(ref_symbol_uncoded);
-        lvlcol = linspecer(numel(classes));
+        % lvlcol = linspecer(numel(classes));
+        lvlcol = cbrewer2('RdBu',numel(classes)+4);
+
+        % remove 4 colors from the middle to avoid the bright central colors
+        mid = ceil(size(lvlcol,1)/2);
+        rm_idx = mid + (-1:2); % two before mid and two after (2x2 removal centered)
+        rm_idx = max(1,min(size(lvlcol,1),rm_idx));
+        lvlcol(rm_idx,:) = [];
+        lvlcol = lvlcol(1:numel(sir_group_labels),:);
 
         for db_lvl = 1:numel(db_constellation)
             db_mask = ref_symbols == db_constellation(db_lvl);

Functions/EQ_visuals/showLevelScatter.m

@@ -21,15 +21,19 @@ arguments
     options.clear (1,1) logical = true
     options.windowLength (1,1) double {mustBePositive, mustBeInteger} = 500
     options.xLimits double = []
-    options.yLimits (1,2) double = [-3 3]
+    options.yLimits (1,2) double = [-2.5 2.5]
     options.showStdAnnotations (1,1) logical = true
+    options.scatterAlpha (1,1) double {mustBeGreaterThanOrEqual(options.scatterAlpha, 0), mustBeLessThanOrEqual(options.scatterAlpha, 1)} = 0.2
+    options.scatterSize (1,1) double {mustBePositive} = 1
+    options.avgLineMaxPoints (1,1) double {mustBePositive, mustBeInteger} = 1000
+    options.avgLineSmoothWindow (1,1) double {mustBePositive, mustBeInteger} = 1
 end
 
 fsym = resolveSymbolRate(rxInput, refSymbols, options);
 [symbols_for_lvl, avg_for_lvl, xAxisUs, info] = prepareLevelScatterData(rxInput, refSymbols, fsym, options);
 
 if options.showPlot
-    plotLevelScatter(symbols_for_lvl, avg_for_lvl, refSymbols, xAxisUs, options);
+    info.plot = plotLevelScatter(symbols_for_lvl, avg_for_lvl, refSymbols, xAxisUs, options);
 end
 end
 
@@ -73,6 +77,7 @@ info.shifts = [];
 info.fsym = fsym;
 info.shiftFs = fsym;
 info.varianceByLevel = [];
+info.plot = struct();
 end
 
 function [symbols_for_lvl, avg_for_lvl, info] = prepareSignalInput(rxSignal, refSymbols, fsym, options, info)
@@ -250,50 +255,95 @@ else
 end
 end
 
-function plotLevelScatter(symbols_for_lvl, avg_for_lvl, refSymbols, xAxisUs, options)
+function plotInfo = plotLevelScatter(symbols_for_lvl, avg_for_lvl, refSymbols, xAxisUs, options)
+plotInfo = struct();
 if isempty(symbols_for_lvl)
     return
 end
 
 if isnan(options.fignum)
-    figure;
+    figHandle = figure;
 else
-    figure(options.fignum);
+    figHandle = figure(options.fignum);
     if options.clear
         clf;
     end
 end
 
-hold on
+ax = gca;
+hold(ax, "on");
 numLevels = size(symbols_for_lvl, 1);
-cols = cbrewer2("Paired", 2*numLevels);
+cols = cbrewer2("RdBu", numLevels);
+% cols = linspecer(numLevels);
+xLimits = getXLimits(options.xLimits, xAxisUs, symbols_for_lvl);
 
-for levelIdx = 1:numLevels
+if 1
-    scatter(xAxisUs, symbols_for_lvl(levelIdx, :), 10, ".", ...
+    for levelIdx = numLevels:-1:1
-        "MarkerEdgeColor", cols((2*levelIdx)-1, :), ...
+        scatter(ax, xAxisUs, symbols_for_lvl(levelIdx, :), options.scatterSize, ".", ...
-        "MarkerEdgeAlpha", 0.5);
+            "MarkerEdgeColor", cols(levelIdx, :), ...
+            "MarkerEdgeAlpha", options.scatterAlpha);
+    end
 end
 
 for levelIdx = 1:numLevels
-    plot(xAxisUs, avg_for_lvl(levelIdx, :), ...
+    [xReduced, yReduced] = reduceAverageLineForPlot( ...
-        "LineWidth", 1, ...
+        xAxisUs, avg_for_lvl(levelIdx, :), ...
-        "Color", cols(2*levelIdx, :));
+        options.avgLineMaxPoints, options.avgLineSmoothWindow);
+    plot(ax, xReduced, yReduced, ...
+        "LineWidth", 2, ...
+        "Color", cols(levelIdx, :));
 end
 
-refValues = numericSignal(refSymbols);
+% refValues = unique(numericSignal(refSymbols));
-yline(unique(refValues), "HandleVisibility", "off");
+% for refIdx = 1:numel(refValues)
+%     yline(ax, refValues(refIdx), "--", ...
+%         "Color", [0.35 0.35 0.35], ...
+%         "LineWidth", 0.75, ...
+%         "HandleVisibility", "off");
+% end
+
 if options.showStdAnnotations
     annotateLevelStd(symbols_for_lvl, avg_for_lvl, xAxisUs, options.xLimits);
 end
 
-xlabel('Time in $\mu$s');
+xlabel(ax, 'Time in $\mu$s');
-ylabel('Normalized Amplitude');
+ylabel(ax, 'Normalized Amplitude');
-xlim(getXLimits(options.xLimits, xAxisUs, symbols_for_lvl));
+xlim(ax, xLimits);
-ylim(options.yLimits);
+ylim(ax, options.yLimits);
-grid on
+grid(ax, "on");
+figure(figHandle);
 drawnow;
 end
 
+function [xReduced, yReduced] = reduceAverageLineForPlot(xAxisUs, yTrace, maxPoints, smoothWindow)
+valid = isfinite(xAxisUs) & isfinite(yTrace);
+xValid = xAxisUs(valid);
+yValid = yTrace(valid);
+
+if isempty(xValid)
+    xReduced = [];
+    yReduced = [];
+    return
+end
+
+[xValid, uniqueIdx] = unique(xValid, "stable");
+yValid = yValid(uniqueIdx);
+
+if smoothWindow > 1
+    yValid = smoothdata(yValid, "movmean", smoothWindow);
+end
+
+numOut = min(maxPoints, numel(xValid));
+if numOut >= numel(xValid)
+    xReduced = xValid;
+    yReduced = yValid;
+    return
+end
+
+xReduced = linspace(xValid(1), xValid(end), numOut);
+yReduced = interp1(xValid, yValid, xReduced, "pchip");
+end
+
 function xLimits = getXLimits(configuredLimits, xAxisUs, symbols_for_lvl)
 if ~isempty(configuredLimits)
     xLimits = configuredLimits;
@@ -308,12 +358,12 @@ end
 
 lastFilledColumn = find(filledColumns, 1, "last");
 xMax = xAxisUs(lastFilledColumn);
-xLimits = [0, 1.05*xMax];
+xLimits = [0, 1*xMax];
 end
 
 function annotateLevelStd(symbols_for_lvl, avg_for_lvl, xAxisUs, configuredXLimits)
 xLimits = getXLimits(configuredXLimits, xAxisUs, symbols_for_lvl);
-xText = xLimits(1) + 0.96*diff(xLimits);
+xText = xLimits(1) + 0.2*diff(xLimits);
 
 for levelIdx = 1:size(symbols_for_lvl, 1)
     levelSamples = symbols_for_lvl(levelIdx, :);
@@ -328,9 +378,9 @@ for levelIdx = 1:size(symbols_for_lvl, 1)
         yText = median(levelSamples(~isnan(levelSamples)), 'omitnan');
     end
 
-    label = ['std = ', sprintf('%.3f', levelStd)];
+    label = ['$\sigma^2_',sprintf('%d', levelIdx),'$ = ', sprintf('%.3f', levelStd)];
     text(xText, yText, label, ...
-        'Interpreter', 'none', ...
+        'Interpreter', 'latex', ...
         'HorizontalAlignment', 'right', ...
         'VerticalAlignment', 'middle', ...
         'FontSize', 10, ...

Functions/Job_Processing/runBatch.m

@@ -25,6 +25,10 @@ function batchResults = runBatch(workerFcn, jobs, options)
     jobs = normalizeJobs(jobs);
     nJobs = numel(jobs);
     batchResults = cell(1, nJobs);
+    jobDurations = nan(1, nJobs);
+    lastJobDuration = nan;
+    batchStart = tic;
+    effectiveWorkerCount = 1;
     h = [];
 
     if nJobs == 0
@@ -34,29 +38,46 @@ function batchResults = runBatch(workerFcn, jobs, options)
     if options.waitbar
         h = waitbar(0, char(options.waitbarMessage));
         cleanupWaitbar = onCleanup(@() closeWaitbar(h));
+        updateWaitbar(options.waitbar, h, 0, nJobs, jobDurations, ...
+            lastJobDuration, 0, effectiveWorkerCount);
     end
 
     switch mode
         case processingMode.parallel
             pool = setupParallelPool(options.numWorkers, options.idleTimeout, options.cancelExistingQueue);
+            effectiveWorkerCount = min(pool.NumWorkers, nJobs);
             futures = parallel.FevalFuture.empty(nJobs, 0);
 
             for jobIdx = 1:nJobs
                 fprintf('[%s] Submitted.\n', jobs(jobIdx).label);
-                futures(jobIdx) = parfeval(pool, workerFcn, 1, jobs(jobIdx).args{:});
+                futures(jobIdx) = parfeval(pool, @runTimedJob, 3, workerFcn, jobs(jobIdx).args);
             end
 
             consumedIdx = false(nJobs, 1);
             for completedCount = 1:nJobs
                 try
-                    [jobIdx, jobResult] = fetchNext(futures);
+                    [jobIdx, jobResult, jobRuntime, jobError] = fetchNext(futures);
                     consumedIdx(jobIdx) = true;
+                    jobDurations(jobIdx) = jobRuntime;
+                    lastJobDuration = jobRuntime;
+
+                    if isempty(jobError)
                         batchResults{jobIdx} = jobResult;
-                    fprintf('[%s] Completed (%d/%d).\n', jobs(jobIdx).label, completedCount, nJobs);
+                        fprintf('[%s] Completed (%d/%d, runtime %s).\n', ...
+                            jobs(jobIdx).label, completedCount, nJobs, formatDuration(jobRuntime));
 
                         if ~isempty(options.resultHandler)
                             options.resultHandler(jobResult, jobs(jobIdx), jobIdx);
                         end
+                    else
+                        batchResults{jobIdx} = jobError;
+
+                        if ~isempty(options.errorHandler)
+                            options.errorHandler(jobError, jobs(jobIdx), jobIdx);
+                        else
+                            defaultErrorHandler(jobError, jobs(jobIdx), jobIdx);
+                        end
+                    end
                 catch fetchErr
                     errorJobIdx = findErroredFuture(futures, consumedIdx);
                     if isempty(errorJobIdx)
@@ -74,21 +95,28 @@ function batchResults = runBatch(workerFcn, jobs, options)
                     end
                 end
 
-                updateWaitbar(options.waitbar, h, completedCount, nJobs);
+                updateWaitbar(options.waitbar, h, completedCount, nJobs, ...
+                    jobDurations, lastJobDuration, toc(batchStart), effectiveWorkerCount);
             end
 
         case processingMode.serial
             for jobIdx = 1:nJobs
+                jobStart = tic;
                 try
                     fprintf('[%s] Running.\n', jobs(jobIdx).label);
                     jobResult = workerFcn(jobs(jobIdx).args{:});
+                    jobDurations(jobIdx) = toc(jobStart);
+                    lastJobDuration = jobDurations(jobIdx);
                     batchResults{jobIdx} = jobResult;
-                    fprintf('[%s] Completed (%d/%d).\n', jobs(jobIdx).label, jobIdx, nJobs);
+                    fprintf('[%s] Completed (%d/%d, runtime %s).\n', ...
+                        jobs(jobIdx).label, jobIdx, nJobs, formatDuration(jobDurations(jobIdx)));
 
                     if ~isempty(options.resultHandler)
                         options.resultHandler(jobResult, jobs(jobIdx), jobIdx);
                     end
                 catch ME
+                    jobDurations(jobIdx) = toc(jobStart);
+                    lastJobDuration = jobDurations(jobIdx);
                     batchResults{jobIdx} = ME;
 
                     if ~isempty(options.errorHandler)
@@ -98,7 +126,8 @@ function batchResults = runBatch(workerFcn, jobs, options)
                     end
                 end
 
-                updateWaitbar(options.waitbar, h, jobIdx, nJobs);
+                updateWaitbar(options.waitbar, h, jobIdx, nJobs, ...
+                    jobDurations, lastJobDuration, toc(batchStart), effectiveWorkerCount);
             end
 
         otherwise
@@ -149,16 +178,69 @@ function jobs = normalizeJobs(jobs)
     end
 end
 
-function updateWaitbar(useWaitbar, h, completedCount, totalCount)
+function [jobResult, jobRuntime, jobError] = runTimedJob(workerFcn, jobArgs)
+    jobStart = tic;
+    jobError = [];
+
+    try
+        jobResult = workerFcn(jobArgs{:});
+    catch ME
+        jobResult = [];
+        jobError = ME;
+    end
+
+    jobRuntime = toc(jobStart);
+end
+
+function updateWaitbar(useWaitbar, h, completedCount, totalCount, jobDurations, lastJobDuration, elapsedSeconds, effectiveWorkerCount)
     if ~useWaitbar || ~isgraphics(h)
         return
     end
 
-    waitbar(completedCount / totalCount, h, ...
+    completedDurations = jobDurations(~isnan(jobDurations));
-        sprintf('Completed %d/%d jobs', completedCount, totalCount));
+    if isempty(completedDurations)
+        averageJobTimeText = 'calculating...';
+        lastJobTimeText = 'calculating...';
+        remainingTimeText = 'calculating...';
+    else
+        averageJobSeconds = mean(completedDurations);
+        remainingJobs = totalCount - completedCount;
+        remainingSeconds = averageJobSeconds * remainingJobs / max(1, effectiveWorkerCount);
+        averageJobTimeText = formatDuration(averageJobSeconds);
+        lastJobTimeText = formatDuration(lastJobDuration);
+        remainingTimeText = formatDuration(remainingSeconds);
+    end
+
+    message = sprintf(['Completed %d/%d jobs\n' ...
+        'Elapsed total: %s\n' ...
+        'Last job time: %s\n' ...
+        'Avg. time/job: %s\n' ...
+        'Estimated remaining: %s'], ...
+        completedCount, totalCount, formatDuration(elapsedSeconds), ...
+        lastJobTimeText, averageJobTimeText, remainingTimeText);
+
+    waitbar(completedCount / totalCount, h, message);
     drawnow;
 end
 
+function text = formatDuration(seconds)
+    if isempty(seconds) || isnan(seconds) || isinf(seconds)
+        text = 'unknown';
+        return
+    end
+
+    seconds = max(0, seconds);
+    hours = floor(seconds / 3600);
+    minutes = floor(mod(seconds, 3600) / 60);
+    wholeSeconds = floor(mod(seconds, 60));
+
+    if hours > 0
+        text = sprintf('%d:%02d:%02d', hours, minutes, wholeSeconds);
+    else
+        text = sprintf('%02d:%02d', minutes, wholeSeconds);
+    end
+end
+
 function closeWaitbar(h)
     if isgraphics(h)
         delete(h);

Functions/beautifyBERplot.m

@@ -19,66 +19,85 @@ end
 ax = gca;
 
 % ---------------------------------------------------------
-% Extract *only* real data lines from ax.Children
+% Extract real data objects from ax.Children
 % ---------------------------------------------------------
 children = ax.Children;
 isLine   = arrayfun(@(h) isa(h,'matlab.graphics.chart.primitive.Line'), children);
-lines    = children(isLine);
+isScatter = arrayfun(@(h) isa(h,'matlab.graphics.chart.primitive.Scatter'), children);
 
 % Remove helper lines (e.g. yline, fit overlays)
-lines = lines(~strcmp({lines.Tag},'ConstantLine'));
+isConstantLine = arrayfun(@(h) isprop(h,'Tag') && strcmp(h.Tag,'ConstantLine'), children);
 
-n = numel(lines);
+dataObjects = children((isLine | isScatter) & ~isConstantLine);
-if n == 0
+n = numel(dataObjects);
-    return
-end
 
 % ---------------------------------------------------------
 % Fixed color palette (deterministic across figures)
 % ---------------------------------------------------------
-try
+if n > 0
+    try
         cmap = linspecer(n);
-catch
+    catch
-    cmap = lines(n).Color; %#ok<NASGU>
         cmap = linespecer_fallback(n);
+    end
 end
 
 markers = {'o','s','o','o','^','v','d','>'};
 lw = 0.8;
-ms = 4;
+ms = 2;
 
 % ---------------------------------------------------------
-% Apply line + marker styling
+% Apply line/scatter + marker styling
 % ---------------------------------------------------------
 for i = 1:n
-    ln = lines(n-i+1);   % preserve plotting order
+    dataObj = dataObjects(n-i+1);   % preserve plotting order
-    ln.LineWidth = lw;
+
+    if isa(dataObj,'matlab.graphics.chart.primitive.Line')
+        dataObj.LineWidth = lw;
 
         if options.setcolors
-        ln.Color = cmap(i,:);
+            dataObj.Color = cmap(i,:);
         end
 
         if options.setmarkers
-        if strcmp(ln.Marker,'none')
+            if strcmp(dataObj.Marker,'none')
-            ln.Marker = markers{mod(i-1,numel(markers))+1};
+                dataObj.Marker = markers{mod(i-1,numel(markers))+1};
             end
-        ln.MarkerSize = ms;
+            dataObj.MarkerSize = ms;
             if options.setcolors
-            ln.MarkerEdgeColor = cmap(i,:);
+                dataObj.MarkerEdgeColor = cmap(i,:);
+            end
+            dataObj.MarkerFaceColor = [1 1 1];
+        end
+    elseif isa(dataObj,'matlab.graphics.chart.primitive.Scatter')
+        if options.setcolors
+            dataObj.CData = cmap(i,:);
+            dataObj.MarkerEdgeColor = cmap(i,:);
+            markerFaceColor = dataObj.MarkerFaceColor;
+            hasNoFace = ischar(markerFaceColor) && strcmp(markerFaceColor,'none');
+            if ~hasNoFace
+                dataObj.MarkerFaceColor = cmap(i,:);
+            end
+        end
+
+        if options.setmarkers
+            if strcmp(dataObj.Marker,'none')
+                dataObj.Marker = markers{mod(i-1,numel(markers))+1};
+            end
+            dataObj.SizeData = ms^2;
         end
-        ln.MarkerFaceColor = [1 1 1];
     end
 end
 
 % ---------------------------------------------------------
 % Optional smoothing / fitting overlay
 % ---------------------------------------------------------
-if options.polyfit
+if options.polyfit && n > 0
     hold on
     for i = 1:n
-        ln = lines(n-i+1);
+        dataObj = dataObjects(n-i+1);
-        x = ln.XData(:);
+        x = dataObj.XData(:);
-        y = ln.YData(:);
+        y = dataObj.YData(:);
         valid = isfinite(x) & isfinite(y);
 
         if sum(valid) < options.polyorder+1
@@ -112,7 +131,15 @@ if options.polyfit
                 yf = polyval(p,xf);
         end
 
-        lightcol = ln.Color + 0.4*(1-ln.Color);
+        if isa(dataObj,'matlab.graphics.chart.primitive.Line')
+            basecol = dataObj.Color;
+        else
+            basecol = dataObj.CData(1,:);
+            if numel(basecol) ~= 3
+                basecol = cmap(i,:);
+            end
+        end
+        lightcol = basecol + 0.4*(1-basecol);
         lightcol(lightcol>1)=1;
 
         plot(xf,yf,'-','Color',lightcol,...

Functions/mat2tikz_improved.m

@@ -3,9 +3,13 @@ arguments
     % Default to the path in your example if no argument is provided
     filename (1,1) string = 'C:\Users\Silas\Documents\Dissertation\00_Examples\tikz\textfig.tikz';
     options.cleanfigure = false;
+    options.cleanTargetResolution = 600;
+    options.cleanScalePrecision (1,1) double = 1;
 end
 if options.cleanfigure
-    cleanfigure;
+    cleanfigure( ...
+        'targetResolution', options.cleanTargetResolution, ...
+        'scalePrecision', options.cleanScalePrecision);
 end
 matlab2tikz(char(filename), ...
     'width', '\fwidth', ...

projects/Diss/MPI_revisit/analytic_mpi_evaluation_with_levels.m

@@ -0,0 +1,261 @@
+%% MPI level-dependent variance: analytic + Monte-Carlo PAM-4
+
+% clear; clc;
+
+%% Parameters
+df       = 1e6;                    % Laser linewidth [Hz]
+SIR_dB   = 20;                     % Interference attenuation [dB], field attenuation
+alpha    = 10^(-SIR_dB/20);        % Field attenuation [linear]
+n_fiber  = 1.467;
+c        = physconst('lightspeed');
+
+L        = linspace(0,250,50);     % Interference delay [m]
+tau      = n_fiber./c.*L;          % Delay [s]
+
+tau_c    = 1/(pi*df);
+L_c      = (c/n_fiber)*tau_c;
+
+%% PAM-4 optical power levels
+
+M = 4;
+ER_dB = 10;
+r_ER = 10^(ER_dB/10);              % P3/P0
+
+% Equally spaced optical powers with finite ER
+P0 = 1/(r_ER - 1);
+P_levels = P0 + (0:M-1).'/(M-1);   % optical powers
+Pavg = mean(P_levels);
+
+level_axis = (0:M-1).';
+
+%% Analytic level-dependent variance
+
+coh_factor = (1 - exp(-2*pi*df.*tau)).^2;
+
+var_levels_analytic = zeros(M, length(tau));
+
+for m = 1:M
+    var_levels_analytic(m,:) = ...
+        2 * alpha^2 * P_levels(m) * Pavg .* coh_factor;
+end
+var_levels_saturation = 2 * alpha^2 * P_levels * Pavg;
+
+sigma_levels_analytic = sqrt(var_levels_analytic);
+
+slope_vs_index_analytic = zeros(1, length(tau));
+
+for k = 1:length(tau)
+    p_k = polyfit(level_axis, var_levels_analytic(:,k), 1);
+    slope_vs_index_analytic(k) = p_k(1);
+end
+
+%% Monte-Carlo simulation
+
+fs                  = 100e9;             % Sampling rate [Hz]
+Tsim                = 100e-6;             % Simulation duration [s]
+N                   = round(Tsim*fs);
+max_delay_samples   = round(max(tau)*fs);
+phase_noise_std     = sqrt(2*pi*df/fs);
+
+num_realizations    = 50;
+
+use_rms_normalize   = false;             % false: comparable to analytic model
+
+mc_var_levels   = zeros(num_realizations, M, length(L));
+mc_sigma_levels = zeros(num_realizations, M, length(L));
+mc_slope        = zeros(num_realizations, length(L));
+
+parfor rr = 1:num_realizations
+
+    % Random PAM symbols
+    sym_idx = randi(M, N + max_delay_samples, 1);   % column
+    P_seq   = P_levels(sym_idx);                    % column
+    A_seq   = sqrt(P_seq);                          % column
+
+    % Laser phase noise random walk
+    dphi = phase_noise_std * randn(N + max_delay_samples, 1);  % column
+    phi  = cumsum(dphi);                                      % column
+
+    P_direct   = P_seq(max_delay_samples+1 : max_delay_samples+N);
+    A_direct   = A_seq(max_delay_samples+1 : max_delay_samples+N);
+    phi_direct = phi(max_delay_samples+1 : max_delay_samples+N);
+
+    varlev_k = zeros(M, length(L));
+    slope_k  = zeros(1, length(L));
+
+    for kk = 1:length(tau)
+
+        nd = round(tau(kk)*fs);
+
+        A_delayed   = A_seq(max_delay_samples+1-nd : max_delay_samples+N-nd);
+        A_delayed = sqrt(Pavg) * ones(N,1);
+        phi_delayed = phi(max_delay_samples+1-nd : max_delay_samples+N-nd);
+
+        onlybeating = 0;
+        if ~onlybeating
+            E = A_direct .* exp(1j*phi_direct) + ...
+                alpha .* A_delayed .* exp(1j*phi_delayed);
+
+            responsivity = 1;
+            I = abs(E).^2 .* responsivity;
+        else
+            I_B = 2 * alpha .* A_direct .* A_delayed .* ...
+                cos(phi_direct - phi_delayed);
+
+            I = I_B;
+        end
+
+        if use_rms_normalize
+            I = I - mean(I);
+            I = I ./ rms(I);
+        end
+
+        % Remove deterministic level means before variance estimation
+        for m = 1:M
+            mask = sym_idx(max_delay_samples+1 : max_delay_samples+N).' == m;
+
+            I_m = I(mask);
+            I_m = I_m - mean(I_m, 'omitnan');
+
+            varlev_k(m,kk) = var(I_m, 0, 'omitnan');
+        end
+
+        p_k = polyfit(level_axis, varlev_k(:,kk), 1);
+        slope_k(kk) = p_k(1);
+    end
+
+    mc_var_levels(rr,:,:) = varlev_k;
+    mc_sigma_levels(rr,:,:) = sqrt(varlev_k);
+    mc_slope(rr,:)        = slope_k;
+end
+
+avg_mc_var_levels   = squeeze(mean(mc_var_levels, 1));
+std_mc_var_levels   = squeeze(std(mc_var_levels, 0, 1));
+avg_mc_sigma_levels = squeeze(mean(mc_sigma_levels, 1));
+std_mc_sigma_levels = squeeze(std(mc_sigma_levels, 0, 1));
+
+avg_mc_slope = mean(mc_slope, 1);
+std_mc_slope = std(mc_slope, 0, 1);
+
+%% Plot 1: variance per PAM level
+
+script_dir = fileparts(mfilename('fullpath'));
+if isempty(script_dir)
+    repo_root = pwd;
+else
+    repo_root = fileparts(fileparts(fileparts(script_dir)));
+end
+addpath(genpath(fullfile(repo_root, 'Libs', 'boundedlines')));
+
+cols = flip(cbrewer2('RdBu',4));
+x_norm = L./L_c;
+bound_alpha = 0.32;
+
+figure; hold on;
+
+for m = 1:M
+    showLegendEntry = m == 1;
+
+    hAnalytic = plot(x_norm, var_levels_analytic(m,:), ...
+        'LineWidth', 2, ...
+        'Color', cols(m,:), ...
+        'DisplayName', 'Analytic');
+    if ~showLegendEntry
+        set(hAnalytic, 'HandleVisibility', 'off');
+    end
+
+    [hl, hp] = boundedline(x_norm, avg_mc_var_levels(m,:), std_mc_var_levels(m,:), '-o', ...
+        'alpha', ...
+        'transparency', bound_alpha, ...
+        'Color', cols(m,:), ...
+        'LineWidth', 1.2);
+    set(hl, 'DisplayName', 'Simulation', 'Marker','o', 'MarkerFaceColor', cols(m,:), 'MarkerSize',2);
+    if ~showLegendEntry
+        set(hl, 'HandleVisibility', 'off');
+    end
+    set(hp, 'HandleVisibility', 'off');
+
+    hSat = yline(var_levels_saturation(m), '-.k', ...
+        'LineWidth', 1.2, ...
+        'DisplayName', 'Saturation $= 2\alpha^2 A^2 A_m^2$');
+    if ~showLegendEntry
+        set(hSat, 'HandleVisibility', 'off');
+    end
+
+    text(x_norm(end) * 0.98, var_levels_saturation(m), sprintf('$m = %d$', m), ...
+        'Interpreter', 'latex', ...
+        'HorizontalAlignment', 'right', ...
+        'VerticalAlignment', 'bottom', ...
+        'Color', 'k', ...
+        'FontSize', 10);
+end
+
+grid on;
+xlabel('$L/L_c$', 'Interpreter', 'latex');
+ylabel('$\sigma_m^2$', 'Interpreter', 'latex');
+title(sprintf('MPI-induced PAM-4 level variance, SIR = %d dB, ER = %g dB', SIR_dB, ER_dB));
+legend('Location', 'northwest', 'Interpreter', 'latex');
+
+%% Plot 2: standard deviation per PAM level
+
+figure; hold on;
+
+for m = 1:M
+    showLegendEntry = m == 1;
+
+    hAnalytic = plot(x_norm, sigma_levels_analytic(m,:), ...
+        'LineWidth', 2, ...
+        'Color', cols(m,:), ...
+        'DisplayName', 'Analytic');
+    if ~showLegendEntry
+        set(hAnalytic, 'HandleVisibility', 'off');
+    end
+
+    [hl, hp] = boundedline(x_norm, avg_mc_sigma_levels(m,:), std_mc_sigma_levels(m,:), '-o', ...
+        'alpha', ...
+        'transparency', bound_alpha, ...
+        'Color', cols(m,:), ...
+        'LineWidth', 1.2);
+    set(hl, 'DisplayName', 'Simulation', 'Marker','o', 'MarkerFaceColor', cols(m,:), 'MarkerSize',2);
+    if ~showLegendEntry
+        set(hl, 'HandleVisibility', 'off');
+    end
+    set(hp, 'HandleVisibility', 'off');
+
+    text(x_norm(end) * 0.98, sigma_levels_analytic(m,end), sprintf('$m = %d$', m), ...
+        'Interpreter', 'latex', ...
+        'HorizontalAlignment', 'right', ...
+        'VerticalAlignment', 'bottom', ...
+        'Color', cols(m,:), ...
+        'FontSize', 10);
+end
+
+grid on;
+xlabel('$L/L_c$', 'Interpreter', 'latex');
+ylabel('$\sigma_m$', 'Interpreter', 'latex');
+title(sprintf('MPI-induced PAM-4 level standard deviation, SIR = %d dB, ER = %g dB', SIR_dB, ER_dB));
+legend('Location', 'northwest', 'Interpreter', 'latex');
+
+%% Plot 3: slope of level variance
+
+figure; hold on;
+
+plot(x_norm, slope_vs_index_analytic, ...
+    'LineWidth', 2, ...
+    'DisplayName', 'Analytic');
+bound_alpha =0.1;
+[hl, hp] = boundedline(x_norm, avg_mc_slope, std_mc_slope, '-o', ...
+    'alpha', ...
+    'transparency', bound_alpha, ...
+    'LineWidth', 1.2);
+set(hl, 'DisplayName', 'Monte Carlo');
+set(hp, 'HandleVisibility', 'off');
+
+grid on;
+xlabel('$L/L_c$', 'Interpreter', 'latex');
+ylabel('Slope of level variance');
+title('MPI-induced variance slope');
+legend('Location', 'southeast');
+
+%%
+% mat2tikz_improved("C:\Users\Silas\Documents\6971e0b65b380ca6d71c837f\02_IMDD_System\tikz\mpi\analytical_mpi_variance_per_level_v2.tikz");

projects/Diss/MPI_revisit/calculate_and_store_std.m

@@ -0,0 +1,135 @@
+function calculate_and_store_std()
+   
+
+    db = DBHandler("type","mysql","dataBase",'labor');
+    savePath = 'W:\labdata\ECOC Silas\ecoc_2025\';
+    
+    fp = QueryFilter();
+    M = 4;
+    % fp.where('Runs', 'pam_level','EQUALS', M);
+    % fp.where('Runs', 'symbolrate','EQUALS', 112e9);
+    fp.where('Runs', 'fiber_length','EQUALS', 0);
+    % fp.where('Runs', 'is_mpi','EQUALS', 1);
+    % fp.where('Runs','run_id','GREATER_EQUAL',3153);
+    fp.where('Runs', 'db_mode',                   'EQUALS', '"no_db"');
+    
+    [dataTable, ~] = db.queryDB(fp, db.getTableFieldNames('Runs'));
+    
+    [~, uniqueIdx] = unique(dataTable.run_id);
+    dataTable = dataTable(uniqueIdx,:);
+    
+    nRuns = height(dataTable);
+    
+    std_levels  = nan(M, nRuns);
+    mean_levels = nan(M, nRuns);
+    p           = nan(2, nRuns);   % p(1,:) slope, p(2,:) offset
+    
+    parfor i = 1:nRuns
+    
+        dataTable_ = dataTable(i,:);
+        run_id = dataTable_.run_id;
+
+        if ~rawLevelStatsAreNaN(dataTable_)
+            fprintf('Skipped run_id %d (%d/%d): raw level stats already stored\n', run_id, i, nRuns);
+            continue
+        end
+    
+        fsym = dataTable_.symbolrate;
+        % Optional, only if needed later
+        % duob_mode = db_mode.(strrep(char(dataTable_.db_mode),'"',''));
+    
+        % Load TX symbols
+        try
+            Symbols_tmp = load(fullfile(savePath, char(dataTable_.tx_symbols_path)));
+            Symbols = Symbols_tmp.Symbols;
+    
+    
+            % Load RX raw signal
+            rxPath = char(dataTable_.rx_raw_path);
+            Scpe_tmp = load(fullfile(savePath, rxPath));
+            Scpe_sig_raw = Scpe_tmp.Scpe_sig_raw;
+    
+            [separated_levels_sig, ~] = showLevelScatter(Scpe_sig_raw, Symbols, ...
+                "fsym", fsym, ...
+                "syncFs", 2*fsym, ...
+                "fignum", 400, ...
+                "normalize", true, ...
+                "debug_plots", false, ...
+                "showPlot", false);
+    
+            separated_levels_sig = flip(separated_levels_sig,1);
+    
+            std_i  = std(separated_levels_sig, 0, 2, 'omitnan');
+            mean_i = mean(separated_levels_sig, 2, 'omitnan');
+    
+            db_local = DBHandler("type","mysql","dataBase",'labor');
+            updateRawLevelStats(db_local, run_id, std_i, mean_i);
+    
+            fprintf('Added to db: run_id %d (%d/%d)\n', run_id, i, nRuns);
+        catch
+            fprintf('Failed run_id %d (%d/%d)\n', run_id, i, nRuns);
+            continue
+        end
+    
+    end
+    
+    
+    
+    %%
+    
+    for i = 1:nRuns
+        % Fit variance over empirical level means
+        p(:,i) = polyfit(mean_levels(:,i), std_levels(:,i), 1).';
+    end
+    
+    var_slope  = p(1,:).';
+    var_offset = p(2,:).';
+
+end
+
+function updateRawLevelStats(db, run_id, std_i, mean_i)
+    stdJson = escapeSqlString(jsonencode(std_i(:).'));
+    meanJson = escapeSqlString(jsonencode(mean_i(:).'));
+    
+    query = sprintf([ ...
+        'UPDATE Runs ' ...
+        'SET std_rawlevels = ''%s'', mean_rawlevels = ''%s'' ' ...
+        'WHERE run_id = %d'], ...
+        stdJson, meanJson, run_id);
+    
+    db.executeSQL(query);
+end
+    
+function escaped = escapeSqlString(value)
+    escaped = strrep(char(value), '''', '''''');
+end
+
+function answer = rawLevelStatsAreNaN(runRow)
+    answer = isDbNaN(runRow.std_rawlevels) || isDbNaN(runRow.mean_rawlevels);
+end
+
+function answer = isDbNaN(value)
+    if iscell(value)
+        if isempty(value)
+            answer = true;
+            return
+        end
+        value = value{1};
+    end
+
+    if isempty(value)
+        answer = true;
+    elseif isnumeric(value)
+        answer = isempty(value) || all(isnan(value), 'all');
+    elseif isstring(value) || ischar(value)
+        value = strtrim(string(value));
+        answer = value == "" || strcmpi(value, "NaN") || strcmpi(value, "null");
+    else
+        try
+            answer = any(ismissing(value), 'all');
+        catch
+            answer = false;
+        end
+    end
+end
+

projects/Diss/MPI_revisit/investigate_database.m

@@ -0,0 +1,155 @@
+
+
+dsp_options.database_type = "mysql";
+dsp_options.dataBase = "labor";
+dsp_options.server = "192.168.178.192";
+dsp_options.port = 3306;
+dsp_options.user = "silas";
+dsp_options.password = "silas";
+
+db = DBHandler("dataBase", dsp_options.dataBase, "type", dsp_options.database_type, ...
+    "server", dsp_options.server, "port", dsp_options.port, ...
+    "user", dsp_options.user, "password", dsp_options.password);
+
+fp = QueryFilter();
+% fp.where('Runs', 'symbolrate','EQUALS', 112e9);
+[dataTable,~] = db.queryDB(fp, db.getTableFieldNames('Runs'));
+
+%% TOTAL RUNS
+
+totalRuns = height(dataTable);
+
+%% Difference between "opti" and not opti (before and after run_id = 2661)
+
+fp = QueryFilter();
+fp.where('Runs','run_id','LESS_THAN',3153);
+fp.where('Runs', 'db_mode','EQUALS', '"no_db"');
+% fp.where('Runs', 'interference_path_length','EQUALS', 300);
+[dataTable_1,~] = db.queryDB(fp, db.getTableFieldNames('Runs'));
+
+fp = QueryFilter();
+fp.where('Runs','run_id','GREATER_EQUAL',3153);
+% fp.where('Runs', 'interference_path_length','EQUALS', 300);
+[dataTable_2,~] = db.queryDB(fp, db.getTableFieldNames('Runs'));
+
+%% prepare x/y and plotting without forcing datetime priority for x
+% simplified: assume columns exist and no datetime special cases
+xCol = 'symbolrate';
+yCol = 'v_awg';
+
+% get x and y (fallback to NaN vectors of appropriate length)
+if ismember(xCol, dataTable_1.Properties.VariableNames)
+    x1 = dataTable_1.(xCol);
+else
+    x1 = NaN(height(dataTable_1),1);
+end
+if ismember(xCol, dataTable_2.Properties.VariableNames)
+    x2 = dataTable_2.(xCol);
+else
+    x2 = NaN(height(dataTable_2),1);
+end
+
+if ismember(yCol, dataTable_1.Properties.VariableNames)
+    y1 = dataTable_1.(yCol);
+else
+    y1 = NaN(height(dataTable_1),1);
+end
+if ismember(yCol, dataTable_2.Properties.VariableNames)
+    y2 = dataTable_2.(yCol);
+else
+    y2 = NaN(height(dataTable_2),1);
+end
+
+% convert text to numeric where needed
+toNum = @(v) double(string(v));
+if iscell(x1) || isstring(x1) || ischar(x1), x1 = toNum(x1); end
+if iscell(x2) || isstring(x2) || ischar(x2), x2 = toNum(x2); end
+if iscell(y1) || isstring(y1) || ischar(y1), y1 = toNum(y1); end
+if iscell(y2) || isstring(y2) || ischar(y2), y2 = toNum(y2); end
+
+% simple scatter plots side-by-side
+figure;
+subplot(1,2,1);
+scatter(x1, y1, 36, 'b', 'filled');
+xlabel(xCol); ylabel(yCol); title('dataTable\_1'); grid on;
+
+subplot(1,2,2);
+scatter(x2, y2, 36, 'r', 'filled');
+xlabel(xCol); ylabel(yCol); title('dataTable\_2'); grid on;
+
+linkaxes(findall(gcf,'Type','axes'),'y');
+
+%% Count interference path lengths by PAM and symbolrate
+
+pamCol = 'pam_level';
+iplCol = 'interference_path_length';
+rateCol = 'symbolrate';
+pamLevels = [2, 4, 6, 8];
+
+% Current analysis uses only the post-opti/post-threshold runs.
+% Use [dataTable_1; dataTable_2] here if the pre-threshold runs should be included.
+allData = dataTable_2;
+
+requiredCols = {pamCol, iplCol, rateCol};
+missingCols = setdiff(requiredCols, allData.Properties.VariableNames);
+if ~isempty(missingCols)
+    error('Missing required column(s): %s', strjoin(missingCols, ', '));
+end
+
+pam = double(string(allData.(pamCol)));
+ipl = double(string(allData.(iplCol)));
+symbolrate = double(string(allData.(rateCol)));
+
+validRows = ~isnan(pam) & ~isnan(ipl) & ~isnan(symbolrate);
+iplValues = unique(ipl(validRows));
+symbolrateValues = unique(symbolrate(validRows));
+
+if isempty(iplValues) || isempty(symbolrateValues)
+    error('No valid rows found for PAM/interference_path_length/symbolrate plotting.');
+end
+
+iplCategories = categorical(string(iplValues), string(iplValues), string(iplValues));
+symbolrateLabels = compose('%.0f GBd', symbolrateValues*1e-9);
+
+figure('Units', 'normalized', 'Position', [0.05 0.1 0.9 0.7]);
+tiledlayout(1, 5, 'TileSpacing', 'compact', 'Padding', 'compact');
+
+for k = 1:numel(pamLevels)
+    ax = nexttile;
+    pamMask = validRows & pam == pamLevels(k);
+    counts = zeros(numel(iplValues), numel(symbolrateValues));
+    for rateIdx = 1:numel(symbolrateValues)
+        for iplIdx = 1:numel(iplValues)
+            counts(iplIdx, rateIdx) = sum(pamMask ...
+                & ipl == iplValues(iplIdx) ...
+                & symbolrate == symbolrateValues(rateIdx));
+        end
+    end
+
+    bar(ax, iplCategories, counts, 'stacked');
+    xlabel(ax, 'interference\_path\_length');
+    ylabel(ax, 'count');
+    title(ax, sprintf('PAM %d', pamLevels(k)));
+    xtickangle(ax, 45);
+    grid(ax, 'on');
+end
+
+axAll = nexttile;
+countsAll = zeros(numel(iplValues), numel(symbolrateValues));
+for rateIdx = 1:numel(symbolrateValues)
+    for iplIdx = 1:numel(iplValues)
+        countsAll(iplIdx, rateIdx) = sum(validRows ...
+            & ipl == iplValues(iplIdx) ...
+            & symbolrate == symbolrateValues(rateIdx));
+    end
+end
+
+bar(axAll, iplCategories, countsAll, 'stacked');
+xlabel(axAll, 'interference\_path\_length');
+ylabel(axAll, 'count');
+title(axAll, 'All PAMs');
+xtickangle(axAll, 45);
+grid(axAll, 'on');
+
+lgd = legend(axAll, symbolrateLabels, 'Location', 'eastoutside');
+title(lgd, 'symbolrate');

projects/Diss/MPI_revisit/investigate_mpi_algorithms.m

@@ -3,9 +3,9 @@ dsp_options = struct();
 dsp_options.mode = "run_id";
 dsp_options.recipe = @mpi_recipe_dev;
 dsp_options.append_to_db = false;
-dsp_options.max_occurences = 3;
 dsp_options.start_occurence = 1;
-dsp_options.debug_plots = true;
+dsp_options.max_occurences = 10;
+dsp_options.debug_plots = false;
 
 dsp_options.database_type = "mysql";
 
@@ -24,21 +24,21 @@ db = DBHandler("dataBase", [dsp_options.dataBase],...
 
 %%
 fp = QueryFilter();
+fp.where('Runs','run_id','GREATER_EQUAL',3153);
 fp.where('Runs', 'symbolrate',                'EQUALS', 112e9);
 fp.where('Runs', 'fiber_length',              'EQUALS', 0);
 fp.where('Runs', 'interference_path_length',  'EQUALS', 1000);
-fp.where('Runs', 'sir',                       'EQUALS', 23);
+fp.where('Runs', 'sir',                       'LESS_EQUAL', 30);
-% fp.where('Runs', 'db_mode',                   'EQUALS', '"no_db"');
+fp.where('Runs', 'db_mode',                   'EQUALS', '"no_db"');
 % fp.where('Runs', 'is_mpi',                    'EQUALS', 1);
 fp.where('Runs', 'pam_level',                 'EQUALS', 4);
-fp.where('Runs', 'wavelength',                'EQUALS', 1310);
+% fp.where('Runs', 'v_bias',                    'EQUALS', 2.65);
-fp.where('Runs', 'v_bias',                    'EQUALS', 2.65);
 
 [dataTable,~] = db.queryDB(fp, db.getTableFieldNames('Runs'));
 
-[~, uniqueSirRows] = unique(dataTable.sir, "stable");
+% [~, uniqueSirRows] = unique(dataTable.sir, "stable");
-dataTable = dataTable(uniqueSirRows, :);
+% dataTable = dataTable(uniqueSirRows, :);
-% sort rows by sir (ascending) and extract run_ids
+
 [~, sortIdx] = sort(dataTable.sir, 'descend');
 dataTable = dataTable(sortIdx, :);
 run_ids = dataTable.run_id;
@@ -50,27 +50,55 @@ end
 
 %% === Warehouse setup ===
 dsp_options.userParameters = struct();
-dsp_options.userParameters.smoothing_length = logspace(1,5.5,5);
+% dsp_options.userParameters.smoothing_length = [0,linspace(100,1000,10),2500,5000,10000];
-
 wh = DataStorage(dsp_options.userParameters);
 
+
+%%
+n_realizations = (dsp_options.max_occurences - dsp_options.start_occurence + 1);
+n_userparams = prod(wh.dim);
+n_run_ids = numel(run_ids);
+parallel_jobs = n_userparams * n_run_ids;
+queried_jobs = n_realizations * n_userparams * n_run_ids;
+
+fprintf("-> [ %d run_id(s) × %d userParam combination(s) = %d parallel job(s) ] × %d realizations = %d total jobs \n", ...
+    n_run_ids, n_userparams, parallel_jobs, n_realizations, queried_jobs);
+
 %% === Run ===
 % submitJobs returns the raw per-job results and the filled Warehouse. For a
 % single run_id and remove_dc = [0, 1], results is a 2-by-1 cell array.
-[results, wh] = submitJobs(run_ids, dsp_options, processingMode.parallel, ...
+[results, wh] = submitJobs(run_ids, dsp_options, processingMode.serial, ...
     "wh", wh, ...
     "waitbar", true);
 
 %% Analyze
 
 storageNames = fieldnames(wh.sto);
-x_vars = dsp_options.userParameters.smoothing_length;
+x_vars = run_ids;%dsp_options.userParameters.smoothing_length;
+x_vars = reshape(x_vars,1,[]);
 
 figure(2026);hold on
 for i = 1:numel(run_ids)
     disp(run_ids(i))
-    result = wh.getStoValue(storageNames{1},x_vars);
+    result = wh.getStoValue(storageNames{1}, x_vars);
-    ber = cellfun(@(packageCell) packageCell{1}.metrics.BER, result);
+    % Each cell in result is a cell array (packageCell) containing multiple packages.
-    plot(x_vars,ber);
+    % Extract BERs from all packages and, for example, average them per x_var.
-    beautifyBERplot("logscale",true,"setcolors",true,"setmarkers",true);
+    ber_all = cellfun(@(packageCell) cellfun(@(pkg) pkg.metrics.BER, packageCell), result, 'UniformOutput', false);
+    % Convert to numeric matrix (rows = packages, cols = x_vars) by padding if needed
+    maxPackages = max(cellfun(@numel, ber_all));
+    ber_mat = nan(maxPackages, numel(ber_all));
+    for k = 1:numel(ber_all)
+        ber_mat(1:numel(ber_all{k}), k) = ber_all{k};
+    end
+    % Choose aggregation: mean across packages (ignore NaNs)
+    ber = mean(ber_mat, 1, 'omitnan');
+    
+    
+    x = dataTable.sir;
+    y = ber;
+    
+    plot(x,ber);
+    scatter(x,ber_mat)
+    beautifyBERplot("logscale",true,"setcolors",false,"setmarkers",true);
 end
+

projects/Diss/MPI_revisit/investigate_results.m

@@ -0,0 +1,86 @@
+% === DSP settings ===
+dsp_options.database_type = "mysql";
+dsp_options.dataBase = "labor";
+dsp_options.server = "192.168.178.192";
+dsp_options.port = 3306;
+dsp_options.user = "silas";
+dsp_options.password = "silas";
+
+db = DBHandler("dataBase", dsp_options.dataBase, "type", dsp_options.database_type, ...
+    "server", dsp_options.server, "port", dsp_options.port, ...
+    "user", dsp_options.user, "password", dsp_options.password);
+
+%%
+fp = QueryFilter();
+fp.where('Runs','run_id','GREATER_EQUAL',3153);
+% fp.where('Runs', 'symbolrate',                'EQUALS', 112e9);
+fp.where('Runs', 'fiber_length',              'EQUALS', 0);
+% fp.where('Runs', 'interference_path_length',  'EQUALS', 300);
+% fp.where('Runs', 'sir',                       'EQUALS', 23);
+fp.where('Runs', 'db_mode',                   'EQUALS', '"no_db"');
+% fp.where('Runs', 'is_mpi',                    'EQUALS', 1);
+% fp.where('Runs', 'pam_level',                 'EQUALS', 4);
+% fp.where('Runs', 'v_bias',                    'EQUALS', 2.65);
+
+[dataTable,~] = db.queryDB(fp, db.getTableFieldNames('mpi_superview'));
+dataTable_save = dataTable;
+
+%%
+
+% No compensation method
+dataTable = dataTable(dataTable.dc_buffer_len ==            1, :);
+dataTable = dataTable(dataTable.ffe_buffer_len ==           1, :);
+dataTable = dataTable(dataTable.smoothing_buffer_len ==     0, :);
+dataTable = dataTable(dataTable.smoothing_buffer_update ==  0, :);
+dataTable = dataTable(dataTable.DCmu ==                     0, :);
+
+% BER over SIR
+plotTable = sortrows(dataTable, 'sir');
+validBer = ~isnan(plotTable.sir) & ~isnan(plotTable.BER) & plotTable.BER > 0 & plotTable.BER < 0.1;
+
+figure;
+hold on;
+scatter(plotTable.sir(validBer), plotTable.BER(validBer), 36, ...
+    'filled', ...
+    'DisplayName', 'BER');
+
+xlabel('SIR [dB]');
+ylabel('BER');
+title(sprintf('BER over SIR (%d filtered rows)', height(dataTable)));
+legend('Location', 'best');
+beautifyBERplot("setcolors", true, "setmarkers", true,"logscale",true,"polyfit",false);
+
+%%
+
+
+% Any Value over SIR
+varname = 'BER';
+groupvar = 'pam_level';
+requiredCols = {'sir', varname, groupvar};
+missingCols = setdiff(requiredCols, dataTable.Properties.VariableNames);
+if ~isempty(missingCols)
+    error('Missing required column(s) for grouped plot: %s', strjoin(missingCols, ', '));
+end
+
+plotTable = sortrows(dataTable, 'sir');
+validline = ~isnan(plotTable.sir) ...
+    & ~isnan(plotTable.(varname)) ...
+    & ~isnan(plotTable.(groupvar));% & plotTable.BER > 0 & plotTable.BER < 0.1;
+groupValues = unique(plotTable.(groupvar)(validline));
+
+figure;
+hold on;
+for groupIdx = 1:numel(groupValues)
+    currentGroup = groupValues(groupIdx);
+    groupRows = validline & plotTable.(groupvar) == currentGroup;
+
+    scatter(plotTable.sir(groupRows),  plotTable.(varname)(groupRows), 36, ...
+        'filled', ...
+        'DisplayName', sprintf('PAM %g', currentGroup));
+end
+
+xlabel('SIR [dB]');
+ylabel(varname);
+title(sprintf('%s over SIR grouped by %s (%d filtered rows)', varname, groupvar, height(dataTable)));
+legend('Location', 'best');
+beautifyBERplot("setcolors", true, "setmarkers", true,"logscale",true,"polyfit",false);

projects/Diss/MPI_revisit/plot_mpi_timesignal.m

@@ -0,0 +1,453 @@
+
+% Diesen PLot habe ich in der Diss als Tabelle
+
+% Create a compact 3-column by 4-row MPI time-signal figure table in LaTeX.
+% 
+% Use these plots from 04_Experimental_Evaluation/tikz/mpi/:
+% - columns: SIR = 20 dB, 30 dB, 50 dB
+% - rows: path length = 1 m, 20 m, 300 m, 1000 m
+% - file pattern: mpi_timesignal_<pathlen>m_<sir>_db.tikz
+% 
+% Use a tabular layout with thin black table gridlines:
+% - \arrayrulewidth = 0.2pt
+% - \arrayrulecolor{black}
+% - vertical and horizontal rules around all cells
+% - rotated row labels for 1 m, 20 m, 300 m, 1000 m
+% - center the rotated row labels vertically against the plot cells
+% - use a dedicated left-column macro so the SIR 20 dB plots are slightly wider:
+% - left column fwidth = 0.255\textwidth
+% - inner columns fwidth = 0.240\textwidth
+% - all fheight = 0.102\textheight
+% 
+% For the TikZ plot files used in the table:
+% - leftmost/SIR 20 dB plots keep the y-label:
+% ylabel={Norm. Amplitude}
+% - SIR 30 dB and 50 dB plots have no y-label:
+% ylabel={}
+% - SIR 30 dB and 50 dB plots hide y ticks:
+% ytick=\empty
+% - axis labels and tick labels use scriptsize:
+% every axis/.append style={font=\scriptsize}
+% xlabel style={font=\color{white!15!black}\scriptsize}
+% ylabel style={font=\color{white!15!black}\scriptsize}
+% - sigma annotation node boxes use tiny:
+% \node[font=\tiny, ...]
+%     - node boxes:
+%     inner sep=0.5pt
+%     at (axis cs:10,...)
+%     sigma labels use normal math subscripts, e.g. $\sigma^2_2$, not $\sigma^2\_2$
+% - all plotted line widths and grid style widths were normalized as needed:
+% addplot line width=1.0pt
+% - grid should be visually off:
+% remove xmajorgrids and ymajorgrids
+% keep style definitions available:
+% grid style={line width=0.4pt, solid, color=black!20}
+% minor grid style={line width=0.2pt, solid, color=black!10}
+% 
+% db = DBHandler("type","mysql","dataBase",'labor');
+% savePath = 'W:\labdata\ECOC Silas\ecoc_2025\';
+
+
+for sirval = [20,30,50]
+    for pathlen = [1,20,300,1000]
+    
+        fp = QueryFilter();
+        % fp.where('Runs', 'run_id','EQUALS', 987);
+        M = 4;
+        fp.where('Runs', 'pam_level','EQUALS', M);     
+        % fp.where('Runs', 'symbolrate','NOT_EQUAL', 112e9);
+        fp.where('Runs', 'fiber_length','EQUALS', 0);   
+        fp.where('Runs', 'is_mpi','EQUALS', 1); 
+        
+        fp.where('Runs', 'interference_path_length','EQUALS', pathlen);      
+        % fp.where('Runs', 'loop_id','GREATER_THAN', 11);  
+    
+        fp.where('Runs', 'sir','GREATER_EQUAL',sirval);  
+        fp.where('Runs','run_id','GREATER_EQUAL',3153);
+        
+        % Sort results by SIR value ascending and reorder corresponding rows
+        [dataTable,sql_query] = db.queryDB(fp, db.getTableFieldNames('Runs'));
+        [~, idx] = sort(dataTable.sir, 'ascend', 'MissingPlacement', 'last');
+        dataTable = dataTable(idx, :);
+        
+        dataTable = dataTable(1,:);
+        dataTable_ = dataTable;
+        [~, uniqueIdx] = unique(dataTable.run_id);  % Get unique run_id indices
+        
+        % dataTable.SIR = -7 - round(dataTable.power_mpi_interference);
+        
+        fsym = dataTable.symbolrate;
+        M = double(dataTable_.pam_level);
+        duob_mode = db_mode.(strrep(char(dataTable_.db_mode),'"',''));
+        
+        Tx_signal = load([savePath, char(dataTable_.tx_signal_path),'.mat']);
+        Tx_signal = Tx_signal.Digi_sig;
+        
+        Tx_bits     = load([savePath, char(dataTable_.tx_bits_path)]);
+        Tx_bits     = Tx_bits.Bits;
+        Symbols_mapped = PAMmapper(M,0).map(Tx_bits);
+        Symbols_mapped.fs = dataTable_.symbolrate;
+        
+        Symbols  = load([savePath, char(dataTable_.tx_symbols_path)]);
+        Symbols  = Symbols.Symbols;
+        
+        Scpe_sig_raw = load([savePath, char(dataTable_.rx_raw_path(1))]);
+        Scpe_sig_raw = Scpe_sig_raw.Scpe_sig_raw;
+        Scpe_sig_raw = Scpe_sig_raw.normalize("mode","rms");
+        % Scpe_sig_raw.plot("displayname",['Scope Signal (Run ID: ',num2str(dataTable_.run_id)],"fignum",dataTable_.run_id,"clear",0);
+        % 
+        % Scpe_sig_raw.spectrum("normalizeTo0dB",1);
+        
+        
+        disp(num2str(dataTable_.interference_path_length));
+        
+        % Scpe_sig_raw.eye(dataTable.symbolrate,dataTable.pam_level);
+        
+        [separated_levels_sig, avg_sig, plotInfo] = showLevelScatter(Scpe_sig_raw, Symbols, ...
+            "fsym", fsym, ...
+            "syncFs", 2*fsym, ...
+            "fignum", sirval+pathlen, ...
+            "normalize", true, ...
+            "debug_plots", false, ...
+            "showPlot", true, ...
+            "showStdAnnotations", true, ...
+            "yLimits", [-2.5 2.5], ...
+            "xLimits",[0 25],...
+            "scatterAlpha", 0.25, ...
+            "scatterSize", 1, ...
+            "avgLineMaxPoints", 500, ...
+            "avgLineSmoothWindow", 5);
+            
+        title(sprintf("%d m, %d db", pathlen,sirval));
+        
+        exportDir = "C:/Users/Silas/Documents/6971e0b65b380ca6d71c837f/04_Experimental_Evaluation/tikz/mpi";
+        exportName = string(sprintf("mpi_timesignal_v2_%dm_%d_db", pathlen,sirval));
+        pngFile = fullfile(exportDir, exportName + "-1.png");
+        tikzFile = fullfile(exportDir, exportName + ".tikz");
+        pngTikzPath = "04_Experimental_Evaluation/tikz/mpi/" + exportName + "-1.png";
+        
+        figure(sirval + pathlen);
+        exportNakedPlotWithTikz(gca, pngFile, tikzFile, pngTikzPath, 150);
+    
+    end
+end
+%%
+
+
+
+% std_levels = std(separated_levels_sig,0,2,'omitnan');
+% mean_levels = mean(separated_levels_sig,2,'omitnan');
+% std_tot = std(Scpe_sig_raw.signal,0,1,'omitnan');
+% 
+% var_levels = std_levels.^2;
+% p = polyfit(mean_levels, var_levels, 1);
+% 
+% var_slope  = p(1);
+% var_offset = p(2);c
+
+function exportNakedPlotWithTikz(sourceAx, pngFile, tikzFile, pngTikzPath, resolutionDpi)
+outputDir = fileparts(pngFile);
+if ~exist(outputDir, "dir")
+    mkdir(outputDir);
+end
+
+exportScatterLayerPng(sourceAx, pngFile, resolutionDpi);
+
+linePlots = collectTikzLinePlots(sourceAx);
+textAnnotations = collectTikzTextAnnotations(sourceAx);
+writeTikzPngAxis(tikzFile, pngTikzPath, sourceAx, linePlots, textAnnotations);
+end
+
+function exportScatterLayerPng(sourceAx, pngFile, resolutionDpi)
+sourceFig = ancestor(sourceAx, "figure");
+xLimits = sourceAx.XLim;
+yLimits = sourceAx.YLim;
+scatterHandles = findGraphicsObjectsByType(sourceAx, "scatter");
+
+exportFig = figure( ...
+    "Visible", "off", ...
+    "Color", "white", ...
+    "Units", sourceFig.Units, ...
+    "Position", sourceFig.Position);
+cleanupFigure = onCleanup(@() close(exportFig));
+
+exportAx = copyobj(sourceAx, exportFig);
+exportAx.Units = "normalized";
+exportAx.Position = [0 0 1 1];
+exportAx.XLim = xLimits;
+exportAx.YLim = yLimits;
+exportAx.XLimMode = "manual";
+exportAx.YLimMode = "manual";
+stripAxesToScatterOnly(exportAx);
+hideAxesInkForRasterExport(exportAx);
+
+fprintf("Exporting scatter PNG with XLim=[%g %g], YLim=[%g %g], scatter objects=%d: %s\n", ...
+    xLimits(1), xLimits(2), yLimits(1), yLimits(2), numel(scatterHandles), pngFile);
+drawnow;
+exportFig.PaperPositionMode = "auto";
+print(exportFig, char(pngFile), "-dpng", sprintf("-r%d", resolutionDpi));
+end
+
+function stripAxesToScatterOnly(ax)
+plotObjects = findall(ax);
+for objIdx = 1:numel(plotObjects)
+    obj = plotObjects(objIdx);
+    if obj == ax || ~isvalid(obj)
+        continue
+    end
+
+    objType = string(get(obj, "Type"));
+    if lower(objType) ~= "scatter"
+        delete(obj);
+    end
+end
+end
+
+function hideAxesInkForRasterExport(ax)
+ax.Visible = "on";
+ax.Color = "white";
+ax.Box = "off";
+ax.XGrid = "off";
+ax.YGrid = "off";
+ax.XMinorGrid = "off";
+ax.YMinorGrid = "off";
+ax.XTick = [];
+ax.YTick = [];
+ax.XColor = "white";
+ax.YColor = "white";
+ax.Title.String = "";
+ax.XLabel.String = "";
+ax.YLabel.String = "";
+end
+
+function handles = findGraphicsObjectsByType(parentHandle, objectType)
+allHandles = findall(parentHandle);
+matches = false(size(allHandles));
+for handleIdx = 1:numel(allHandles)
+    currentHandle = allHandles(handleIdx);
+    if ~isvalid(currentHandle)
+        continue
+    end
+
+    matches(handleIdx) = lower(string(get(currentHandle, "Type"))) == lower(string(objectType));
+end
+
+handles = allHandles(matches);
+end
+
+function linePlots = collectTikzLinePlots(sourceAx)
+lineHandles = findGraphicsObjectsByType(sourceAx, "line");
+lineHandles = flipud(lineHandles(:));
+
+linePlots = repmat(struct( ...
+    "xData", [], ...
+    "yData", [], ...
+    "color", [], ...
+    "lineWidth", []), numel(lineHandles), 1);
+validLineCount = 0;
+
+for lineIdx = 1:numel(lineHandles)
+    lineHandle = lineHandles(lineIdx);
+    xData = lineHandle.XData(:);
+    yData = lineHandle.YData(:);
+    validSamples = isfinite(xData) & isfinite(yData);
+
+    if ~any(validSamples)
+        continue
+    end
+
+    validLineCount = validLineCount + 1;
+    linePlots(validLineCount).xData = xData(validSamples);
+    linePlots(validLineCount).yData = yData(validSamples);
+    linePlots(validLineCount).color = lineHandle.Color;
+    linePlots(validLineCount).lineWidth = lineHandle.LineWidth;
+end
+
+linePlots = linePlots(1:validLineCount);
+end
+
+function textAnnotations = collectTikzTextAnnotations(sourceAx)
+textHandles = findGraphicsObjectsByType(sourceAx, "text");
+textHandles = flipud(textHandles(:));
+xLimits = sourceAx.XLim;
+yLimits = sourceAx.YLim;
+
+textAnnotations = repmat(struct( ...
+    "x", [], ...
+    "y", [], ...
+    "text", "", ...
+    "anchor", "center"), numel(textHandles), 1);
+validTextCount = 0;
+
+for textIdx = 1:numel(textHandles)
+    textHandle = textHandles(textIdx);
+    labelText = string(textHandle.String);
+    if strlength(strtrim(labelText)) == 0
+        continue
+    end
+
+    textPosition = textHandle.Position;
+    if textPosition(1) < xLimits(1) || textPosition(1) > xLimits(2) || ...
+            textPosition(2) < yLimits(1) || textPosition(2) > yLimits(2)
+        continue
+    end
+
+    validTextCount = validTextCount + 1;
+    textAnnotations(validTextCount).x = textPosition(1);
+    textAnnotations(validTextCount).y = textPosition(2);
+    textAnnotations(validTextCount).text = labelText;
+    textAnnotations(validTextCount).anchor = matlabTextAlignmentToTikzAnchor( ...
+        textHandle.HorizontalAlignment, textHandle.VerticalAlignment);
+end
+
+textAnnotations = textAnnotations(1:validTextCount);
+end
+
+function writeTikzPngAxis(tikzFile, pngTikzPath, sourceAx, linePlots, textAnnotations)
+outputDir = fileparts(tikzFile);
+if ~exist(outputDir, "dir")
+    mkdir(outputDir);
+end
+
+xLimits = sourceAx.XLim;
+yLimits = sourceAx.YLim;
+fid = fopen(tikzFile, "w");
+if fid < 0
+    error("plot_mpi_timesignal:TikzOpenFailed", ...
+        "Could not open TikZ export file: %s", tikzFile);
+end
+cleanupFile = onCleanup(@() fclose(fid));
+
+axisOptions = {
+    "              every axis/.append style={font=\scriptsize},"
+    "width=\fwidth,"
+    "height=\fheight,"
+    "at={(0\fwidth,0\fheight)},"
+    "scale only axis,"
+    "axis on top,"
+    "xmin=" + formatPgfNumber(xLimits(1)) + ","
+    "xmax=" + formatPgfNumber(xLimits(2)) + ","
+    "xlabel style={font=\color{white!15!black}\small},"
+    "xlabel={Time in $\mu$s},"
+    "ymin=" + formatPgfNumber(yLimits(1)) + ","
+    "ymax=" + formatPgfNumber(yLimits(2)) + ","
+    "ylabel style={font=\color{white!15!black}\small},"
+    "ylabel={Normalized Amplitude},"
+    "axis background/.style={fill=white},"
+    "xmajorgrids,"
+    "ymajorgrids,"
+    "grid style={line width=0.4pt, dotted, color=black!20},"
+    "scaled ticks=false,"
+    "tick label style={/pgf/number format/fixed, /pgf/number format/1000 sep={}},"
+    "legend columns=1"
+    };
+
+fprintf(fid, "%% This file was generated by plot_mpi_timesignal.m.\n");
+fprintf(fid, "%%\n");
+for lineIdx = 1:numel(linePlots)
+    fprintf(fid, "%s\n", formatTikzColorDefinition(lineIdx, linePlots(lineIdx).color));
+end
+if ~isempty(linePlots)
+    fprintf(fid, "%%\n");
+end
+fprintf(fid, "%s\n\n", "\begin{tikzpicture}");
+fprintf(fid, "%s\n", "\begin{axis}[%");
+for optionIdx = 1:numel(axisOptions)
+    fprintf(fid, "%s\n", axisOptions{optionIdx});
+end
+fprintf(fid, "%s\n", "]");
+
+graphicsLine = sprintf( ...
+    "\\addplot [forget plot] graphics [xmin=%s, xmax=%s, ymin=%s, ymax=%s] {%s};", ...
+    formatPgfNumber(xLimits(1)), ...
+    formatPgfNumber(xLimits(2)), ...
+    formatPgfNumber(yLimits(1)), ...
+    formatPgfNumber(yLimits(2)), ...
+    char(strrep(pngTikzPath, "\", "/")));
+fprintf(fid, "%s\n\n", graphicsLine);
+writeTikzLinePlots(fid, linePlots);
+writeTikzTextAnnotations(fid, textAnnotations);
+fprintf(fid, "%s\n\n", "\end{axis}");
+fprintf(fid, "%s", "\end{tikzpicture}%");
+end
+
+function valueText = formatPgfNumber(value)
+valueText = string(sprintf("%.15g", value));
+end
+
+function colorDefinition = formatTikzColorDefinition(colorIdx, rgbColor)
+colorDefinition = sprintf( ...
+    "\\definecolor{mycolor%d}{rgb}{%.5f,%.5f,%.5f}%%", ...
+    colorIdx, rgbColor(1), rgbColor(2), rgbColor(3));
+end
+
+function writeTikzLinePlots(fid, linePlots)
+for lineIdx = 1:numel(linePlots)
+    fprintf(fid, "\\addplot [color=mycolor%d, line width=%.1fpt, forget plot]\n", ...
+        lineIdx, linePlots(lineIdx).lineWidth);
+    fprintf(fid, "  table[row sep=crcr]{%%\n");
+
+    for sampleIdx = 1:numel(linePlots(lineIdx).xData)
+        fprintf(fid, "%s\t%s\\\\\n", ...
+            formatPgfNumber(linePlots(lineIdx).xData(sampleIdx)), ...
+            formatPgfNumber(linePlots(lineIdx).yData(sampleIdx)));
+    end
+
+    fprintf(fid, "};\n\n");
+end
+end
+
+function writeTikzTextAnnotations(fid, textAnnotations)
+for textIdx = 1:numel(textAnnotations)
+    fprintf(fid, "\\node[font=\\scriptsize, anchor=%s, fill=white, draw=black!40, rounded corners=2pt, inner sep=2pt]\n", ...
+        textAnnotations(textIdx).anchor);
+    fprintf(fid, "  at (axis cs:%s,%s) {%s};\n\n", ...
+        formatPgfNumber(textAnnotations(textIdx).x), ...
+        formatPgfNumber(textAnnotations(textIdx).y), ...
+        escapeTikzText(textAnnotations(textIdx).text));
+end
+end
+
+function anchor = matlabTextAlignmentToTikzAnchor(horizontalAlignment, verticalAlignment)
+horizontalAlignment = string(horizontalAlignment);
+verticalAlignment = string(verticalAlignment);
+
+if horizontalAlignment == "left"
+    horizontalAnchor = "west";
+elseif horizontalAlignment == "right"
+    horizontalAnchor = "east";
+else
+    horizontalAnchor = "";
+end
+
+if verticalAlignment == "top"
+    verticalAnchor = "north";
+elseif verticalAlignment == "bottom"
+    verticalAnchor = "south";
+else
+    verticalAnchor = "";
+end
+
+if strlength(horizontalAnchor) > 0 && strlength(verticalAnchor) > 0
+    anchor = verticalAnchor + " " + horizontalAnchor;
+elseif strlength(horizontalAnchor) > 0
+    anchor = horizontalAnchor;
+elseif strlength(verticalAnchor) > 0
+    anchor = verticalAnchor;
+else
+    anchor = "center";
+end
+end
+
+function textOut = escapeTikzText(textIn)
+textOut = char(textIn);
+% textOut = strrep(textOut, "\", "\textbackslash{}");
+textOut = strrep(textOut, "{", "\{");
+textOut = strrep(textOut, "}", "\}");
+textOut = strrep(textOut, "%", "\%");
+textOut = strrep(textOut, "&", "\&");
+textOut = strrep(textOut, "#", "\#");
+textOut = strrep(textOut, "_", "\_");
+% textOut = strrep(textOut, "$", "\$");
+end
+

projects/Diss/MPI_revisit/verify_mpi_noise_power.m

@@ -0,0 +1,663 @@
+%% Investigate stored std values per PAM level
+% Assumes mpi_superview contains:
+% run_id, pam_level, sir, interference_path_length, std_rawlevels, mean_rawlevels
+% std_rawlevels / mean_rawlevels stored as JSON arrays, e.g. [0.1,0.2,0.3,0.4]
+
+dsp_options.database_type = "mysql";
+dsp_options.dataBase = "labor";
+dsp_options.server = "192.168.178.192";
+dsp_options.port = 3306;
+dsp_options.user = "silas";
+dsp_options.password = "silas";
+
+db = DBHandler("dataBase", dsp_options.dataBase, "type", dsp_options.database_type, ...
+    "server", dsp_options.server, "port", dsp_options.port, ...
+    "user", dsp_options.user, "password", dsp_options.password);
+
+
+
+
+
+%% Query
+fp = QueryFilter();
+% fp.where('Runs','run_id','GREATER_EQUAL',3153);
+fp.where('Runs', 'fiber_length', 'EQUALS', 0);
+fp.where('Runs', 'db_mode', 'EQUALS', '"no_db"');
+fp.where('Runs', 'interference_path_length', 'NOT_EQUAL', 0);
+fp.where('Runs', 'symbolrate',                'EQUALS', 112e9);
+fp.where('Runs', 'pam_level','EQUALS', 4);
+[dataTable,~] = db.queryDB(fp, db.getTableFieldNames('Runs'));
+
+
+% Basic filtering
+
+% Keep only rows with available level statistics
+valid = hasStoredLevelStats(dataTable);
+
+dataTable = dataTable(valid,:);
+
+% % Optional analysis filters
+% dataTable = dataTable(dataTable.sir < 25, :);
+% dataTable = dataTable(dataTable.pam_level == 4, :);
+
+nRuns = height(dataTable);
+M = 8;
+
+std_levels  = nan(M,nRuns);
+mean_levels = nan(M,nRuns);
+var_levels  = nan(M,nRuns);
+std_slope   = nan(nRuns,1);
+var_slope   = nan(nRuns,1);
+var_offset  = nan(nRuns,1);
+gamma       = nan(nRuns,1);
+
+% Decode JSON arrays and compute variance slopes
+
+for i = 1:nRuns
+
+    std_i  = jsondecode(char(dataTable.std_rawlevels(i)));
+    mean_i = jsondecode(char(dataTable.mean_rawlevels(i)));
+
+    std_i  = std_i(:);
+    mean_i = mean_i(:);
+
+    % Keep only valid entries
+    valid_i = isfinite(std_i) & isfinite(mean_i);
+
+    if ~sum(valid_i)
+        continue
+    end
+    std_i  = std_i(valid_i);
+    mean_i = mean_i(valid_i);
+
+    % Sort by measured level mean
+    [mean_i, idx] = sort(mean_i, 'ascend');
+    std_i = std_i(idx);
+
+    var_i = std_i.^2;
+
+    % Store
+    std_levels(1:numel(std_i),i)  = std_i;
+    mean_levels(1:numel(mean_i),i) = mean_i;
+    var_levels(1:numel(var_i),i)  = var_i;
+
+    % 1) Fit variance over measured raw level mean
+    p_k = polyfit(mean_i, var_i, 1);
+
+    % 2) Normalized x-axis slope, better cross-PAM comparable
+    x_i = (mean_i - min(mean_i)) ./ (max(mean_i) - min(mean_i));
+    p_std = polyfit(x_i, std_i, 1);
+    std_slope(i) = p_std(1);
+
+    var_slope(i)  = p_k(1);
+    var_offset(i) = p_k(2);
+
+    gamma(i) = std_i(end) / std_i(1);
+end
+
+% Add calculated values to table
+dataTable.var_slope_calc  = var_slope;
+dataTable.std_slope_calc  = std_slope;
+dataTable.var_offset_calc = var_offset;
+dataTable.gamma_calc      = gamma;
+
+dataTable.sir = -7 - dataTable.power_mpi_interference;
+
+
+
+
+%% Plot 1: variance per level at fixed path length, grouped by PAM level
+
+% NICHT GENUTZT
+path_length_filter = -1;
+idx_path = dataTable.interference_path_length > path_length_filter;
+pam_levels = unique(dataTable.pam_level(idx_path));
+cols = linspecer(numel(pam_levels));
+figure; hold on;
+sir_values_plot = dataTable.sir(idx_path & isfinite(dataTable.sir));
+if isempty(sir_values_plot)
+    sir_limits = [0, 1];
+else
+    sir_limits = [min(sir_values_plot), max(sir_values_plot)];
+end
+if diff(sir_limits) == 0
+    sir_limits = sir_limits + [-0.5, 0.5];
+end
+sir_cmap = cbrewer2('RdBu',256);%parula(256);
+colormap(sir_cmap);
+caxis(sir_limits);
+mkr = ['+','o','*','square'];
+for g = 1:numel(pam_levels)
+
+    pamLevel = pam_levels(g);
+    idx_g = idx_path & dataTable.pam_level == pamLevel;
+
+    mean_g = mean_levels(:,idx_g);
+    var_g = std_levels(:,idx_g);
+    run_sir_g = dataTable.sir(idx_g);
+    sir_g = repmat(dataTable.sir(idx_g).', M, 1);
+    run_id_g = repmat(dataTable.run_id(idx_g).', M, 1);
+    path_length_g = repmat(dataTable.interference_path_length(idx_g).', M, 1);
+    pam_level_g = repmat(dataTable.pam_level(idx_g).', M, 1);
+
+    % Scatter all runs
+    hScatter = scatter(mean_g(:), var_g(:), 5, sir_g(:), ...
+        'filled', ...
+        'MarkerFaceAlpha', 0.55, ...
+        'MarkerEdgeAlpha', 0.35, ...
+        'HandleVisibility','off');
+
+    for runIdx = 1:5:size(mean_g, 2)
+        if isfinite(run_sir_g(runIdx))
+            sirColorIdx = round(interp1(sir_limits, [1, size(sir_cmap,1)], run_sir_g(runIdx), 'linear', 'extrap'));
+            sirColorIdx = min(max(sirColorIdx, 1), size(sir_cmap,1));
+            lineColor = sir_cmap(sirColorIdx,:);
+        else
+            lineColor = [0.5, 0.5, 0.5];
+        end
+
+        hLine = plot(mean_g(:,runIdx), var_g(:,runIdx), ...
+            'LineWidth', 0.1, ...
+            'LineStyle', '-', ...
+            'Marker', 'none', ...
+            'HandleVisibility', 'off');
+        try
+            hLine.Color = [lineColor, 0.25];
+        catch
+            hLine.Color = lineColor;
+        end
+    end
+
+    hScatter.DataTipTemplate.DataTipRows(1).Label = 'mean level';
+    hScatter.DataTipTemplate.DataTipRows(2).Label = 'variance';
+    hScatter.DataTipTemplate.DataTipRows(3).Label = 'SIR';
+    hScatter.DataTipTemplate.DataTipRows(end+1) = dataTipTextRow('run id', run_id_g(:));
+    hScatter.DataTipTemplate.DataTipRows(end+1) = dataTipTextRow('interference length', path_length_g(:));
+    hScatter.DataTipTemplate.DataTipRows(end+1) = dataTipTextRow('pam_level', pam_level_g(:));
+
+    % Group mean curve
+    mean_level_g = mean(mean_g, 2, 'omitnan');
+    mean_var_g = mean(var_g, 2, 'omitnan');
+
+    plot(mean_level_g, mean_var_g, '-o', ...
+        'Color', [0,0,0], ...
+        'LineWidth', 2, ...
+        'MarkerFaceColor', cols(g,:), ...
+        'MarkerEdgeColor', cols(g,:), ...
+        'DisplayName', sprintf('PAM %.0f', pamLevel), 'MarkerSize', 5,'Marker','o');
+end
+
+grid on;
+xlabel('Mean raw level');
+ylabel('Stored level variance');
+title(sprintf('Level variance at %.0f m path length, grouped by PAM level', path_length_filter));
+legend('Location','best');
+cb = colorbar;
+ylabel(cb, 'SIR / dB');
+
+%% Plot 1b: variance per level at fixed path length, grouped by SIR
+
+sir_filter_max_slope = 60;
+idx_path_sir = idx_path ...
+    & isfinite(dataTable.sir) ...
+    & dataTable.sir < sir_filter_max_slope;
+
+sir_group_edges = [-inf,17.5 20, 22.5, 25, 30, 40, inf];
+sir_group_labels = {'SIR < 17.5 dB','17.5 < SIR < 20 dB', '20 < SIR < 22.5 dB','22.5 < SIR < 25 dB', '25 < SIR < 30 dB', '30 < SIR < 40 dB', sprintf('40 < SIR < %.0f dB', sir_filter_max_slope)};
+sir_group_cols =  cbrewer2('RdBu',numel(sir_group_labels));
+pam_marker_styles = {'o','square','diamond','^','v','>','<','pentagram','hexagram'};
+
+figure; hold on;
+pam_levels_plot = unique(dataTable.pam_level(idx_path_sir));
+max_slope_lines = numel(sir_group_labels) * numel(pam_levels_plot);
+slope_annotation_lines = cell(max_slope_lines, 1);
+slope_annotation_idx = 0;
+fit_line_x = cell(max_slope_lines, 1);
+fit_line_y = cell(max_slope_lines, 1);
+fit_line_color = cell(max_slope_lines, 1);
+fit_line_label = cell(max_slope_lines, 1);
+mean_point_x = cell(max_slope_lines, 1);
+mean_point_y = cell(max_slope_lines, 1);
+mean_point_marker = cell(max_slope_lines, 1);
+fit_text_x = nan(max_slope_lines, 1);
+fit_text_y = nan(max_slope_lines, 1);
+fit_text_label = cell(max_slope_lines, 1);
+fit_line_idx = 0;
+
+for sirGroupIdx = numel(sir_group_labels):-1:1
+    idx_sir_group = idx_path_sir ...
+        & dataTable.sir > sir_group_edges(sirGroupIdx) ...
+        & dataTable.sir < sir_group_edges(sirGroupIdx+1);
+
+    if ~any(idx_sir_group)
+        continue
+    end
+
+    pam_levels_sir = unique(dataTable.pam_level(idx_sir_group));
+    for pamIdx = 1:numel(pam_levels_sir)
+        pamLevel = pam_levels_sir(pamIdx);
+        idx_sir_pam = idx_sir_group & dataTable.pam_level == pamLevel;
+        markerStyle = pam_marker_styles{mod(pamIdx-1,numel(pam_marker_styles))+1};
+
+        mean_g = mean_levels(:,idx_sir_pam);
+        var_g = std_levels(:,idx_sir_pam);
+        sir_g = repmat(dataTable.sir(idx_sir_pam).', M, 1);
+        run_id_g = repmat(dataTable.run_id(idx_sir_pam).', M, 1);
+        path_length_g = repmat(dataTable.interference_path_length(idx_sir_pam).', M, 1);
+        pam_level_g = repmat(dataTable.pam_level(idx_sir_pam).', M, 1);
+
+        hScatter = scatter(mean_g(:), var_g(:), 5, sir_group_cols(sirGroupIdx,:), markerStyle, ...
+            'filled', ...
+            'MarkerFaceAlpha', 0.25, ...
+            'MarkerEdgeAlpha', 0.20, ...
+            'HandleVisibility','off');
+
+        hScatter.DataTipTemplate.DataTipRows(1).Label = 'mean level';
+        hScatter.DataTipTemplate.DataTipRows(2).Label = 'std';
+        hScatter.DataTipTemplate.DataTipRows(end+1) = dataTipTextRow('SIR', sir_g(:));
+        hScatter.DataTipTemplate.DataTipRows(end+1) = dataTipTextRow('run id', run_id_g(:));
+        hScatter.DataTipTemplate.DataTipRows(end+1) = dataTipTextRow('interference length', path_length_g(:));
+        hScatter.DataTipTemplate.DataTipRows(end+1) = dataTipTextRow('pam_level', pam_level_g(:));
+
+        mean_level_g = mean(mean_g, 2, 'omitnan');
+        mean_std_g = mean(var_g, 2, 'omitnan');
+        valid_mean_curve = isfinite(mean_level_g) & isfinite(mean_std_g);
+
+        if nnz(valid_mean_curve) >= 2
+            p_mean_slope = polyfit(mean_level_g(valid_mean_curve), mean_std_g(valid_mean_curve), 1);
+            fit_x_g = linspace(min(mean_level_g(valid_mean_curve)), max(mean_level_g(valid_mean_curve)), 100);
+            fit_y_g = polyval(p_mean_slope, fit_x_g);
+            slope_text_g = strrep(sprintf('%.2f', p_mean_slope(1)), '.', ',');
+            label_x_g = fit_x_g(end);
+            label_y_g = fit_y_g(end);
+        else
+            fit_x_g = [];
+            fit_y_g = [];
+            slope_text_g = 'NaN';
+            label_x_g = NaN;
+            label_y_g = NaN;
+        end
+        slope_annotation_idx = slope_annotation_idx + 1;
+        slope_annotation_lines{slope_annotation_idx} = sprintf('%s, PAM %.0f: %s', sir_group_labels{sirGroupIdx}, pamLevel, slope_text_g);
+
+        fit_line_idx = fit_line_idx + 1;
+        mean_point_x{fit_line_idx} = mean_level_g;
+        mean_point_y{fit_line_idx} = mean_std_g;
+        mean_point_marker{fit_line_idx} = markerStyle;
+        fit_line_x{fit_line_idx} = fit_x_g;
+        fit_line_y{fit_line_idx} = fit_y_g;
+        fit_line_color{fit_line_idx} = sir_group_cols(sirGroupIdx,:);
+        fit_line_label{fit_line_idx} = sprintf('%s, PAM %.0f', sir_group_labels{sirGroupIdx}, pamLevel);
+        fit_text_x(fit_line_idx) = label_x_g;
+        fit_text_y(fit_line_idx) = label_y_g;
+        fit_text_label{fit_line_idx} = sprintf('  m=%s', slope_text_g);
+    end
+end
+
+for lineIdx = 1:fit_line_idx
+    plot(mean_point_x{lineIdx}, mean_point_y{lineIdx}, ...
+        'LineStyle', 'none', ...
+        'Color', fit_line_color{lineIdx}, ...
+        'Marker', 'x', ...
+        'MarkerFaceColor', fit_line_color{lineIdx}, ...
+        'MarkerEdgeColor', fit_line_color{lineIdx}, ...
+        'HandleVisibility', 'off', ...
+        'MarkerSize', 3);
+
+    plot(fit_line_x{lineIdx}, fit_line_y{lineIdx}, '-', ...
+        'Color', fit_line_color{lineIdx}, ...
+        'LineWidth', 0.8, ...
+        'DisplayName', fit_line_label{lineIdx});
+
+    text(fit_text_x(lineIdx), fit_text_y(lineIdx), fit_text_label{lineIdx}, ...
+        'Color', fit_line_color{lineIdx}, ...
+        'FontSize', 8, ...
+        'VerticalAlignment', 'middle', ...
+        'HorizontalAlignment', 'left', ...
+        'Interpreter', 'none');
+end
+
+grid on;
+xlabel('Mean raw level');
+ylabel('Stored level std');
+title(sprintf('Mean level-std curves at %.0f m path length, grouped by SIR and PAM level', path_length_filter));
+legend('Location','eastoutside');
+slope_annotation_lines = slope_annotation_lines(1:slope_annotation_idx);
+
+
+%% Plot 2: slope vs interference path length
+
+figure; hold on;
+valid = dataTable.interference_path_length > 0;
+scatter(dataTable.interference_path_length(valid), dataTable.var_slope_calc(valid), ...
+    45, dataTable.sir(valid), 'filled');
+
+grid on;
+xlabel('Interference path length / m');
+ylabel('Slope of level variance');
+title('Level-variance slope vs. interference path length');
+cb = colorbar;
+ylabel(cb, 'SIR / dB');
+
+%% Plot 3:
+
+figure; hold on;
+
+valid = dataTable.interference_path_length > 50;
+
+sir_group_edges = [-inf,20, 30, inf];
+for sirGroupIdx = 1:numel(sir_group_edges)-1
+    idx_sir_group = validSlopeSir ...
+        & dataTable.sir > sir_group_edges(sirGroupIdx) ...
+        & dataTable.sir < sir_group_edges(sirGroupIdx+1)...
+         & dataTable.interference_path_length > 50;
+
+    scatter(dataTable.sir(idx_sir_group), ...
+        dataTable.var_slope_calc(idx_sir_group), ...
+        8, 'filled', ...
+        'MarkerFaceColor', sir_group_cols(sirGroupIdx,:), ...
+        'MarkerEdgeColor', sir_group_cols(sirGroupIdx,:), ...
+        'MarkerFaceAlpha', 1, ...
+        'MarkerEdgeAlpha', 1, ...
+        'HandleVisibility','off');
+end
+
+set(gca, 'YScale', 'log');
+grid on;
+xlabel('SIR');
+ylabel('Slope of level variance');
+
+
+%% Plot 4:
+
+sir_filter_max_slope = 60;
+validSlopeSir = isfinite(dataTable.interference_path_length) ...
+    & isfinite(dataTable.var_slope_calc) ...
+    & isfinite(dataTable.sir) ...
+    & dataTable.sir < sir_filter_max_slope ...
+    & dataTable.interference_path_length > 0 ...
+    & dataTable.var_slope_calc > 0;
+
+path_lengths_slope = unique(dataTable.interference_path_length(validSlopeSir));
+
+sir_group_edges = [-inf,20, 30, inf];
+sir_group_labels = {'SIR < 20 dB','20 < SIR < 30 dB', sprintf('30 < SIR < %.0f dB', sir_filter_max_slope)};
+tmp_cols = cbrewer2('RdBu', numel(sir_group_labels)+4);
+% remove 4 colors from the middle to avoid the bright central colors
+mid = ceil(size(tmp_cols,1)/2);
+rm_idx = mid + (-1:2); % two before mid and two after (2x2 removal centered)
+rm_idx = max(1,min(size(tmp_cols,1),rm_idx));
+tmp_cols(rm_idx,:) = [];
+sir_group_cols = tmp_cols(1:numel(sir_group_labels),:);
+
+
+figure(6);clf; hold on;
+
+% Plot raw values first. Use one single-color scatter object per SIR group
+% so matlab2tikz does not need unsupported RGB scatter CData.
+for sirGroupIdx = 1:numel(sir_group_labels)
+    idx_sir_group = validSlopeSir ...
+        & dataTable.sir > sir_group_edges(sirGroupIdx) ...
+        & dataTable.sir < sir_group_edges(sirGroupIdx+1);
+
+    scatter(dataTable.interference_path_length(idx_sir_group), ...
+        dataTable.var_slope_calc(idx_sir_group), ...
+        8, 'filled', ...
+        'MarkerFaceColor', sir_group_cols(sirGroupIdx,:), ...
+        'MarkerEdgeColor', sir_group_cols(sirGroupIdx,:), ...
+        'MarkerFaceAlpha', 0.35, ...
+        'MarkerEdgeAlpha', 0.35, ...
+        'HandleVisibility','off');
+end
+
+% Plot grouped averages on top of the raw points.
+for sirGroupIdx = 1:numel(sir_group_labels)
+    idx_sir_group = validSlopeSir ...
+        & dataTable.sir > sir_group_edges(sirGroupIdx) ...
+        & dataTable.sir < sir_group_edges(sirGroupIdx+1);
+
+    mean_slope = nan(numel(path_lengths_slope),1);
+
+    for g = 1:numel(path_lengths_slope)
+        idx_g = idx_sir_group ...
+            & dataTable.interference_path_length == path_lengths_slope(g);
+        slope_g = dataTable.var_slope_calc(idx_g);
+        slope_g = slope_g(isfinite(slope_g));
+        if numel(slope_g) >= 3
+            slope_g = slope_g(~isoutlier(slope_g, 'median'));
+        end
+        mean_slope(g) = mean(slope_g, 'omitnan');
+    end
+
+    valid_mean_slope = isfinite(mean_slope) & mean_slope > 0;
+    plot(path_lengths_slope(valid_mean_slope), mean_slope(valid_mean_slope), ...
+        '-o', ...
+        'Color', sir_group_cols(sirGroupIdx,:), ...
+        'MarkerFaceColor', sir_group_cols(sirGroupIdx,:), ...
+        'MarkerEdgeColor', sir_group_cols(sirGroupIdx,:), ...
+        'LineWidth', 1, ...
+        'DisplayName', sir_group_labels{sirGroupIdx});
+end
+
+xlim([10,1000])
+ylim([0.006,0.1])
+set(gca, 'XScale', 'log');
+set(gca, 'YScale', 'log');
+grid on;
+xlabel('Interference path length / m');
+ylabel('Mean slope of level variance');
+title(sprintf('Grouped level-variance slope vs. interference path length, SIR < %.0f dB', sir_filter_max_slope));
+xticks(path_lengths_slope);
+xticklabels(compose('%.0f', path_lengths_slope));
+legend('Location','northwest');
+beautifyBERplot("setcolors",false)
+% mat2tikz_improved('C:/Users/Silas/Documents/6971e0b65b380ca6d71c837f/04_Experimental_Evaluation/tikz/mpi/interference_path_len_vs_slope_logdomain.tikz');
+
+
+%% Optional: slope vs SIR sanity check
+
+figure; hold on;
+sir_cmap = cbrewer2('RdBu',256);
+
+validSlopeCategory = isfinite(dataTable.interference_path_length) ...
+    & isfinite(dataTable.std_slope_calc) ...
+    & isfinite(dataTable.sir);
+path_lengths_plot = unique(dataTable.interference_path_length(validSlopeCategory));
+sir_values_plot = dataTable.sir(validSlopeCategory);
+if isempty(sir_values_plot)
+    sir_limits = [0, 1];
+else
+    sir_limits = [min(sir_values_plot), max(sir_values_plot)];
+end
+if diff(sir_limits) == 0
+    sir_limits = sir_limits + [-0.5, 0.5];
+end
+colormap(sir_cmap);
+caxis(sir_limits);
+has_symbolrate = ismember('symbolrate', dataTable.Properties.VariableNames);
+
+for g = 1:numel(path_lengths_plot)
+    idx_g = validSlopeCategory ...
+        & dataTable.interference_path_length == path_lengths_plot(g);
+    n_g = sum(idx_g);
+    sir_g = dataTable.sir(idx_g);
+    slope_g = dataTable.var_slope_calc(idx_g);
+    run_id_g = dataTable.run_id(idx_g);
+    pam_level_g = dataTable.pam_level(idx_g);
+    if has_symbolrate
+        symbolrate_g = dataTable.symbolrate(idx_g);
+    else
+        symbolrate_g = nan(n_g,1);
+    end
+
+    if n_g == 1 || range(sir_g) == 0
+        x_g = g * ones(n_g,1);
+    else
+        x_g = g + 0.9 * ((sir_g - min(sir_g)) ./ range(sir_g) - 0.5);
+    end
+
+    hScatter = scatter(x_g, slope_g, ...
+        5, sir_g, 'filled', ...
+        'DisplayName', sprintf('%.0f m', path_lengths_plot(g)));
+    hScatter.DataTipTemplate.DataTipRows(1).Label = 'interference path category';
+    hScatter.DataTipTemplate.DataTipRows(2).Label = 'slope';
+    hScatter.DataTipTemplate.DataTipRows(3).Label = 'SIR';
+    hScatter.DataTipTemplate.DataTipRows(end+1) = dataTipTextRow('run id', run_id_g);
+    hScatter.DataTipTemplate.DataTipRows(end+1) = dataTipTextRow('PAM format', pam_level_g);
+    hScatter.DataTipTemplate.DataTipRows(end+1) = dataTipTextRow('symbol rate', symbolrate_g);
+end
+
+grid on;
+xlabel('Interference path length / m');
+ylabel('Slope of level variance');
+title('Level-variance slope grouped by interference path length');
+xticks(1:numel(path_lengths_plot));
+xticklabels(compose('%.0f', path_lengths_plot));
+cb = colorbar;
+ylabel(cb, 'SIR / dB');
+
+%% Optional: 3D slope vs SIR by interference path length
+
+figure(7);clf; hold on;
+sir_cmap_3d = cbrewer2('RdBu',256);
+
+validSlope3D = isfinite(dataTable.interference_path_length) ...
+    & isfinite(dataTable.var_slope_calc) ...
+    & isfinite(dataTable.sir);
+path_lengths_3d = unique(dataTable.interference_path_length(validSlope3D));
+path_categories_3d = 0:numel(path_lengths_3d)-1;
+sir_values_3d = dataTable.sir(validSlope3D);
+if isempty(sir_values_3d)
+    sir_limits_3d = [0, 1];
+else
+    sir_limits_3d = [min(sir_values_3d), max(sir_values_3d)];
+end
+if diff(sir_limits_3d) == 0
+    sir_limits_3d = sir_limits_3d + [-0.5, 0.5];
+end
+colormap(gca, sir_cmap_3d);
+clim(sir_limits_3d);
+has_symbolrate = ismember('symbolrate', dataTable.Properties.VariableNames);
+
+for g = 1:numel(path_lengths_3d)
+    idx_g = validSlope3D ...
+        & dataTable.interference_path_length == path_lengths_3d(g);
+    n_g = sum(idx_g);
+    sir_g = dataTable.sir(idx_g);
+    slope_g = dataTable.var_slope_calc(idx_g);
+    path_length_g = dataTable.interference_path_length(idx_g);
+    run_id_g = dataTable.run_id(idx_g);
+    pam_level_g = dataTable.pam_level(idx_g);
+    if has_symbolrate
+        symbolrate_g = dataTable.symbolrate(idx_g);
+    else
+        symbolrate_g = nan(n_g,1);
+    end
+
+    path_category_g = path_categories_3d(g) * ones(n_g,1);
+    hScatter = scatter3(path_category_g, sir_g, slope_g, ...
+        8, sir_g, 'filled', ...
+        'DisplayName', sprintf('%.0f m', path_lengths_3d(g)));
+    hScatter.DataTipTemplate.DataTipRows(1).Label = 'interference path category';
+    hScatter.DataTipTemplate.DataTipRows(2).Label = 'SIR';
+    hScatter.DataTipTemplate.DataTipRows(3).Label = 'slope';
+    hScatter.DataTipTemplate.DataTipRows(end+1) = dataTipTextRow('interference path length', path_length_g);
+    hScatter.DataTipTemplate.DataTipRows(end+1) = dataTipTextRow('run id', run_id_g);
+    hScatter.DataTipTemplate.DataTipRows(end+1) = dataTipTextRow('PAM format', pam_level_g);
+    hScatter.DataTipTemplate.DataTipRows(end+1) = dataTipTextRow('symbol rate', symbolrate_g);
+end
+
+set(gca, 'XScale', 'lin');
+set(gca, 'YScale', 'lin');
+set(gca, 'ZScale', 'log');
+grid on;
+box on;
+xlabel('Interference path length / m');
+ylabel('SIR / dB');
+zlabel('Slope');
+title('Level-variance slope vs. SIR by interference path length');
+xticks(path_categories_3d);
+xticklabels(compose('%.0f', path_lengths_3d));
+view([-35 20]);
+cb = colorbar;
+ylabel(cb, 'SIR / dB');
+
+
+%%
+
+
+
+sir_all = -7 - dataTable.power_mpi_interference;
+% for i = 1:height(dataTable)
+%     sigma_all(i,:)   = jsondecode(char(dataTable.std_rawlevels(i)))';
+% end
+slope_var = dataTable.var_slope_calc;
+slope_std = dataTable.std_slope_calc;
+pl = dataTable.interference_path_length;
+col = cbrewer2('RdBu',256);
+% map path lengths to colormap indices
+pl_unique = unique(pl(isfinite(pl)));
+if isempty(pl_unique)
+    pl_idx = ones(size(pl));
+else
+    % normalize pl to [1, size(col,1)]
+    pl_norm = (pl - min(pl_unique)) ./ max(1, (max(pl_unique)-min(pl_unique)));
+    pl_idx = round(1 + pl_norm * (size(col,1)-1));
+    pl_idx(~isfinite(pl_idx)) = 1;
+end
+
+figure(); hold on;
+% scatter slope_var colored by path length
+for i = 1:numel(sir_all)
+    scatter(sir_all(i), slope_var(i), 30, col(pl_idx(i),:), '.', 'MarkerEdgeColor', col(pl_idx(i),:));
+end
+% scatter slope_std colored by path length with different marker edge brightness
+for i = 1:numel(sir_all)
+    scatter(sir_all(i), slope_std(i), 30, col(pl_idx(i),:), '.', 'MarkerEdgeColor', col(pl_idx(i),:));
+end
+cb = colorbar;
+colormap(col);
+caxis([min(pl_unique), max(pl_unique)]);
+ylabel(cb, 'Interference path length / m');
+
+function valid = hasStoredLevelStats(dataTable)
+valid = ~cellfun(@isempty, cellstr(string(dataTable.std_rawlevels))) & ...
+        ~cellfun(@isempty, cellstr(string(dataTable.mean_rawlevels)));
+end
+
+function [stdLevels, meanLevels] = decodeStoredLevelStats(dataTable)
+nRuns = height(dataTable);
+numLevels = 0;
+stdCells = cell(nRuns, 1);
+meanCells = cell(nRuns, 1);
+
+for i = 1:nRuns
+    std_i = jsondecode(char(dataTable.std_rawlevels(i)));
+    mean_i = jsondecode(char(dataTable.mean_rawlevels(i)));
+
+    std_i = std_i(:);
+    mean_i = mean_i(:);
+
+    valid_i = isfinite(std_i) & isfinite(mean_i);
+    std_i = std_i(valid_i);
+    mean_i = mean_i(valid_i);
+
+    [mean_i, idx] = sort(mean_i, 'ascend');
+    std_i = std_i(idx);
+
+    stdCells{i} = std_i;
+    meanCells{i} = mean_i;
+    numLevels = max(numLevels, numel(std_i));
+end
+
+stdLevels = nan(numLevels, nRuns);
+meanLevels = nan(numLevels, nRuns);
+
+for i = 1:nRuns
+    std_i = stdCells{i};
+    mean_i = meanCells{i};
+    stdLevels(1:numel(std_i), i) = std_i;
+    meanLevels(1:numel(mean_i), i) = mean_i;
+end
+end
+

projects/ECOC_2025_MPI/load_signal_standalone.m

@@ -7,19 +7,19 @@ fp = QueryFilter();
 % fp.where('Runs', 'run_id','EQUALS', 987);
 M = 4;
 fp.where('Runs', 'pam_level','EQUALS', M);     
-fp.where('Runs', 'symbolrate','EQUALS', 112e9);
+% fp.where('Runs', 'symbolrate','NOT_EQUAL', 112e9);
 fp.where('Runs', 'fiber_length','EQUALS', 0);   
 fp.where('Runs', 'is_mpi','EQUALS', 1);                
 fp.where('Runs', 'interference_path_length','EQUALS', 1000);      
 % fp.where('Runs', 'loop_id','GREATER_THAN', 11);   
-fp.where('Runs', 'sir','EQUALS',20);  
+fp.where('Runs', 'sir','LESS_EQUAL',20);  
-
+% fp.where('Runs','run_id','GREATER_EQUAL',3153);
 
 [dataTable,sql_query] = db.queryDB(fp, db.getTableFieldNames('Runs'));
 
 [~, uniqueIdx] = unique(dataTable.run_id);  % Get unique run_id indices
 
-dataTable.SIR = -7 - round(dataTable.power_mpi_interference);
+% dataTable.SIR = -7 - round(dataTable.power_mpi_interference);
 
 %%
 for i = 1:size(dataTable,1)
@@ -41,16 +41,29 @@ for i = 1:size(dataTable,1)
     Symbols  = Symbols.Symbols;
 
     Scpe_sig_raw = load([savePath, char(dataTable_.rx_raw_path(1))]);
-    Scpe_sig_raw = Scpe_sig_raw.Scpe_sig_raw.normalize("mode","rms");
+    Scpe_sig_raw = Scpe_sig_raw.Scpe_sig_raw;
+    % Scpe_sig_raw = Scpe_sig_raw.normalize("mode","rms");
     Scpe_sig_raw.plot("displayname",['Scope Signal (Run ID: ',num2str(dataTable_.run_id)],"fignum",dataTable_.run_id,"clear",0);
 
+    Scpe_sig_raw.spectrum("normalizeTo0dB",1);
+
 
     disp(num2str(dataTable_.interference_path_length));
-    [sep_sig, avg_sig] = showLevelScatter(Scpe_sig_raw, Symbols, ...
+    [separated_levels_sig, avg_sig] = showLevelScatter(Scpe_sig_raw, Symbols, ...
         "fsym", fsym, ...
         "syncFs", 2*fsym, ...
         "fignum", 400, ...
         "normalize", true, ...
-        "debug_plots", false);
+        "debug_plots", false,"showPlot",true,"showStdAnnotations",false,"yLimits",[-3 3]);
-    var(sep_sig,0,2,'omitnan')
+
+    std_levels = std(separated_levels_sig,0,2,'omitnan');
+    mean_levels = mean(separated_levels_sig,2,'omitnan');
+    std_tot = std(Scpe_sig_raw.signal,0,1,'omitnan');
+
+    var_levels = std_levels.^2;
+    p = polyfit(mean_levels, var_levels, 1);
+
+    var_slope  = p(1);
+    var_offset = p(2);c
+
 end

projects/ECOC_2025_MPI/plots_from_database/bias_vs_ber.m

@@ -4,43 +4,68 @@
 % database_name = 'ecoc2025.db';
 % database = DBHandler("pathToDB", [basePath, database_name],"type",'mysql');
 
-database = DBHandler("dataBase",'labor',"type","mysql","user","silas","password","silas","server","192.168.178.192");
+% database = DBHandler("dataBase",'labor',"type","mysql","user","silas","password","silas","server","192.168.178.192");
+% 
+% filterParams = database.tables;
+% filterParams.Configurations = struct( ...
+%     'symbolrate',      112e9, ... %[224,336,360,390,420,448]
+%     'fiber_length',     0, ...
+%     'db_mode',        '"no_db"', ...
+%     'interference_attenuation', [], ...
+%     'interference_path_length', [], ...
+%     'is_mpi',           1, ...
+%     'pam_level',        4, ...
+%     'wavelength',       1310, ...
+%     'precomp_amp',      -64, ...
+%     'signal_attenuation', '0', ...
+%     'v_awg',            [], ...
+%     'v_bias',           [] ...
+%     );
+% 
+% % filterParams.EqualizerParameters.diff_precode = int32(db_mode.db_encoded);
+% % filterParams.EqualizerParameters.equalizer_structure = int32(equalizer_structure.vnle);
+% 
+% selectedFields = {'Configurations.run_id' 'Runs.date_of_run' 'Runs.rx_raw_path'	'Configurations.bitrate' 'Configurations.v_bias' 'Configurations.v_awg' 'Configurations.precomp_amp'	'Configurations.symbolrate'	'Configurations.pam_level'...
+%     'Configurations.db_mode'	'Configurations.rop_attenuation'	'Configurations.is_mpi'	'Configurations.interference_attenuation' 'Configurations.signal_attenuation'	...
+%     'EqualizerParameters.equalizer_structure'	'EqualizerParameters.diff_precode' 'EqualizerParameters.eq_id'	'Measurements.power_pd_in' ...
+%     'Measurements.power_mpi_interference'	'Measurements.power_mpi_signal'	'Results.BER' 'Results.BER_precoded'	'Results.SNR'	'Results.GMI'	'Results.Alpha' 'Results.date_of_processing'};
+% 
+% [dataTable,sql_query] = database.queryDB(filterParams, selectedFields);
 
-filterParams = database.tables;
-filterParams.Configurations = struct( ...
-    'symbolrate',      112e9, ... %[224,336,360,390,420,448]
-    'fiber_length',     0, ...
-    'db_mode',        '"no_db"', ...
-    'interference_attenuation', [], ...
-    'interference_path_length', [], ...
-    'is_mpi',           1, ...
-    'pam_level',        4, ...
-    'wavelength',       1310, ...
-    'precomp_amp',      -64, ...
-    'signal_attenuation', '0', ...
-    'v_awg',            [], ...
-    'v_bias',           [] ...
-    );
 
-% filterParams.EqualizerParameters.diff_precode = int32(db_mode.db_encoded);
+%%
-% filterParams.EqualizerParameters.equalizer_structure = int32(equalizer_structure.vnle);
 
-selectedFields = {'Configurations.run_id' 'Runs.date_of_run' 'Runs.rx_raw_path'	'Configurations.bitrate' 'Configurations.v_bias' 'Configurations.v_awg' 'Configurations.precomp_amp'	'Configurations.symbolrate'	'Configurations.pam_level'...
-    'Configurations.db_mode'	'Configurations.rop_attenuation'	'Configurations.is_mpi'	'Configurations.interference_attenuation' 'Configurations.signal_attenuation'	...
-    'EqualizerParameters.equalizer_structure'	'EqualizerParameters.diff_precode' 'EqualizerParameters.eq_id'	'Measurements.power_pd_in' ...
-    'Measurements.power_mpi_interference'	'Measurements.power_mpi_signal'	'Results.BER' 'Results.BER_precoded'	'Results.SNR'	'Results.GMI'	'Results.Alpha' 'Results.date_of_processing'};
 
-[dataTable,sql_query] = database.queryDB(filterParams, selectedFields);
+dsp_options = []; db = DBHandler("dataBase","labor","type","mysql","server","192.168.178.192","port",3306,"user","silas","password","silas");
+
+fp = QueryFilter();
+fp.where('Runs', 'symbolrate',                'EQUALS', 112e9);
+fp.where('Runs', 'fiber_length',              'EQUALS', 0);
+fp.where('Runs', 'interference_path_length',  'EQUALS', 1000);
+% fp.where('Runs', 'sir',                       'EQUALS', 23);
+% fp.where('Runs', 'db_mode',                   'EQUALS', '"no_db"');
+% fp.where('Runs', 'is_mpi',                    'EQUALS', 1);
+fp.where('Runs', 'pam_level',                 'EQUALS', 4);
+fp.where('Runs', 'wavelength',                'EQUALS', 1310);
+% fp.where('Runs', 'v_bias',                    'EQUALS', 2.65);
+
+[dataTable,~] = db.queryDB(fp, db.getTableFieldNames('mpi_superview'));
+
+%%
 
 dataTable = cleanUpTable(dataTable);
 
 % Filter by time
-startTime = datetime('2025-04-11 13:00:00', 'InputFormat', 'yyyy-MM-dd HH:mm:ss');
+
-stopTime  = datetime('2025-04-11 14:00:00', 'InputFormat', 'yyyy-MM-dd HH:mm:ss');
+dataTable.date_of_run = datetime(dataTable.date_of_run, 'InputFormat', 'yyyy-MM-dd HH:mm:ss');
-dataTable.date_of_run = datetime(dataTable.date_of_run, 'InputFormat', 'yyyy-MM-dd HH:mm:ss.SSSSSS');
+dataTable.date_of_processing = datetime(dataTable.date_of_processing, 'InputFormat', 'yyyy-MM-dd HH:mm:ss'); %.SSSSSS
-dataTable.date_of_processing = datetime(dataTable.date_of_processing, 'InputFormat', 'yyyy-MM-dd HH:mm:ss.SSSSSS');
+
-dataTable = dataTable(dataTable.date_of_processing > startTime, :);
+% startTime = datetime('2025-04-11 13:00:00', 'InputFormat', 'yyyy-MM-dd HH:mm:ss');
-dataTable = dataTable(dataTable.date_of_processing < stopTime, :);
+% stopTime  = datetime('2025-04-11 14:00:00', 'InputFormat', 'yyyy-MM-dd HH:mm:ss');
+% dataTable = dataTable(dataTable.date_of_processing > startTime, :);
+% dataTable = dataTable(dataTable.date_of_processing < stopTime, :);
+
+%%
 
 % Group by smth
 y_var = 'SNR';
@@ -52,6 +77,7 @@ dataTableGrpd = groupIt(fixedVars,dataTable);
 
 plotRealizations = 1;
 
+%%
 % Create a new figure
 mkr = 'x';
 

projects/ECOC_2025_MPI/plots_from_database/plot_bias_awg_vs_ber.m

@@ -46,13 +46,17 @@ fn = [db.getTableFieldNames('mpi_superview')];
 
 dataTable = cleanUpTable(dataTable);
 
-% Filter by time
+
-startTime = datetime('2025-04-11 13:00:00', 'InputFormat', 'yyyy-MM-dd HH:mm:ss');
+%% Filter by time
-stopTime  = datetime('2025-04-11 14:00:00', 'InputFormat', 'yyyy-MM-dd HH:mm:ss');
+
-dataTable.date_of_run = datetime(dataTable.date_of_run, 'InputFormat', 'yyyy-MM-dd HH:mm:ss.SSSSSS');
+dataTable.date_of_run = datetime(dataTable.date_of_run, 'InputFormat', 'yyyy-MM-dd HH:mm:ss');
-dataTable.date_of_processing = datetime(dataTable.date_of_processing, 'InputFormat', 'yyyy-MM-dd HH:mm:ss.SSSSSS');
+dataTable.date_of_processing = datetime(dataTable.date_of_processing, 'InputFormat', 'yyyy-MM-dd HH:mm:ss');
-dataTable = dataTable(dataTable.date_of_processing > startTime, :);
+
-dataTable = dataTable(dataTable.date_of_processing < stopTime, :);
+%%
+% startTime = datetime('2025-04-11 13:00:00', 'InputFormat', 'yyyy-MM-dd HH:mm:ss');
+% stopTime  = datetime('2025-04-11 14:00:00', 'InputFormat', 'yyyy-MM-dd HH:mm:ss');
+% dataTable = dataTable(dataTable.date_of_processing > startTime, :);
+% dataTable = dataTable(dataTable.date_of_processing < stopTime, :);
 
 % Group by smth
 y_var = 'BER';
@@ -62,6 +66,8 @@ loop_var = 'eq_id';
 fixedVars = {'eq_id',loop_var};
 dataTableGrpd = groupIt(fixedVars,dataTable);
 
+
+%%
 plotRealizations = 1;
 
 % Create a new figure