MPI plotting stuff

Classes/00_signals/Signal.m

@@ -528,10 +528,7 @@ classdef Signal
             y_margin = 0.05 * y_range;
             ylim([y_min - y_margin, y_max + y_margin]);
 
-            % Set ticks automatically, avoid overpopulation
+            yticks(-200:10:200);
-            try
-            yticks(round(linspace(y_min, y_max, min(10, max(4, ceil(y_range/10))))));
-            end
             grid on;
             
             % Add legend if not already present

Functions/EQ_recipes/mpi_recipe_dev.m

@@ -15,12 +15,13 @@ end
 
 %% Execution toggles
 
-run_conventional_ffe        = 0;
+run_conventional_ffe        = 1;
-run_hpf                     = 1;
+run_ma                      = 0;
-run_a2_tracked_levels       = 0;
+run_hpf                     = 0;
-run_a2_residual             = 0;
+run_a2_tracked_levels       = 1;
-run_a1                      = 0;
+run_a2_residual             = 1;
-run_tracking_adaptive       = 0;
+run_a1                      = 1;
+run_tracking_adaptive       = 1;
 plot_output_signals         = 0;
 
 %% Shared fixed EQ settings
@@ -75,12 +76,68 @@ if run_conventional_ffe
         "plain_ffe", "baseline", block_update);
     output.(char(storageName)) = ffe_results;
 
+    eq_noise_ffe = equalized_signal - Symbols;
+
+
     if plot_output_signals
         plotEqSignals(equalized_signal,Symbols,options,400,-1);
     end
 end
 
-%% Conventional FFE
+
+%% Moving Average Filter -> Conventional FFE
+if run_ma
+    eq_ffe = FFE_plain( ...
+        "sps", eq_sps, ...
+        "order", eq_order, ...
+        "decide", false, ...
+        "adaption_technique", eq_adaption, ...
+        "len_tr", eq_len_tr, ...
+        "epochs_tr", eq_epochs_tr, ...
+        "mu_tr", eq_mu_tr, ...
+        "dd_mode", true, ...
+        "epochs_dd", eq_epochs_dd, ...
+        "mu_dd", 0.0012, ...
+        "optmize_mus", false, ...
+        "plot_mu_optimization", options.debug_plots, ...
+        "save_debug", eq_save_debug);
+
+    %%% APPLY MA filter
+    movavg_window_length = options.userParameters.bufferlen;
+    movavg_window_symmetry = "causal"; % "causal" or "noncausal"
+
+    switch movavg_window_symmetry
+        case "causal"
+            dc_estimate = movmean(Scpe_sig.signal, ...
+                [movavg_window_length-1,0],"Endpoints","shrink");
+        case "noncausal"
+            dc_estimate = movmean(Scpe_sig.signal, ...
+                movavg_window_length,"Endpoints","shrink");
+        otherwise
+            error("mpi_recipe_dev:InvalidMovAvgSymmetry", ...
+                "movavg_window_symmetry must be 'causal' or 'noncausal'.");
+    end
+
+    Scpe_sig_filt = Scpe_sig;
+    Scpe_sig_filt.signal = Scpe_sig.signal - dc_estimate;
+    clear dc_estimate;
+
+    storageName = "MovAvg_plus_conventional_FFE";
+    [ffe_results,equalized_signal] = runFfe(eq_ffe, "Moving-average + conventional FFE", ...
+        Scpe_sig_filt, Symbols, Tx_bits, options);
+    clear Scpe_sig_filt;
+    ffe_results = attachMpiReductionConfig(ffe_results, eq_ffe, storageName, ...
+        "MovAvg_plus_conventional_FFE", char(movavg_window_symmetry), block_update);
+    output.(char(storageName)) = ffe_results;
+
+    if plot_output_signals
+        plotEqSignals(equalized_signal,Symbols,options,400,-1);
+    end
+end
+
+
+
+%% High Pass Filter -> Conventional FFE
 if run_hpf
     eq_ffe = FFE_plain( ...
         "sps", eq_sps, ...
@@ -99,11 +156,12 @@ if run_hpf
 
     %%% APPLY HPF in MHZ order
     f = Filter("f_cutoff",options.userParameters.hpf,"filterType","butterworth","lowpass",0,"filtdegree",4,"fs",options.fsym*2);
-    Scpe_sig = f.process(Scpe_sig);
+    Scpe_sig_filt = f.process(Scpe_sig);
 
     storageName = "HPF_plus_conventional_FFE";
     [ffe_results,equalized_signal] = runFfe(eq_ffe, "Conventional FFE", ...
-        Scpe_sig, Symbols, Tx_bits, options);
+        Scpe_sig_filt, Symbols, Tx_bits, options);
+    clear Scpe_sig_filt;
     ffe_results = attachMpiReductionConfig(ffe_results, eq_ffe, storageName, ...
         "HPF_plus_conventional_FFE", "hpf", block_update);
     output.(char(storageName)) = ffe_results;
@@ -191,7 +249,7 @@ if run_a1
         "epochs_dd", eq_epochs_dd, ...
         "mu_dd", 0.012, ...
         "dc_smoothing_a1", 1, ...
-        "dc_avg_bufferlength_a1", 2048, ...
+        "dc_avg_bufferlength_a1", 2048, ... %was 2048
         "dc_avg_update_blocklength_a1", block_update, ...
         "save_debug", eq_save_debug);
 
@@ -238,7 +296,19 @@ if run_tracking_adaptive
         "dc_tracking", "adaptive", block_update);
     output.(char(storageName)) = ffe_results;
 
+    % fignum = 106;
+    % 
+    % eq_noise = equalized_signal - Symbols;
+    % dn = sprintf("FFE DCT; SIR: %d dB",options.dataTable.sir);
+    % showEQNoisePSD(eq_noise, "fignum", fignum, "displayname", dn,"colormode","diverging");
+    % ylim([-70 -30]);
+    % 
+    % dn = sprintf("FFE only; SIR: %d dB",options.dataTable.sir);
+    % showEQNoisePSD(eq_noise_ffe, "fignum", fignum+1, "displayname", dn,"colormode","diverging");
+    % ylim([-70 -30]);
+
     if plot_output_signals
+
         plotEqSignals(equalized_signal,Symbols,options,450,-1);
     end
 end

Functions/EQ_visuals/showEQNoisePSD.m

@@ -5,8 +5,14 @@ arguments
     options.fignum (1,1) double = NaN % Default to NaN if not provided
     options.displayname (1,:) char = '' % Default to an empty string if not provided
     options.color = [];
+    options.colormode (1,1) string {mustBeMember(options.colormode, ["diverging", "qualitative", "paired"])} = "diverging";
 end
 
+    persistent nextGroupIdx
+    if isempty(nextGroupIdx)
+        nextGroupIdx = 0;
+    end
+
     % Determine the figure number to use or create a new figure
     if isnan(options.fignum)
         fig = figure; % Create a new figure and get its handle
@@ -15,16 +21,16 @@ end
     end
     hold on
     ax = gca;
-    % N = numel(ax.Children);
+    linesBefore = findall(ax, 'Type', 'Line');
-    if isempty(options.color)
+
-        N = sum(arrayfun(@(x) strcmp(x.LineStyle, '-'), ax.Children));
+    useAutomaticColor = isempty(options.color);
-        cmap = linspecer(8);
+    if useAutomaticColor
-        options.color = cmap(mod(N, size(cmap, 1)) + 1, :);
+        options.color = selectEQNoisePSDColor(ax, options.colormode);
     end
 
     % Ensure the figure is ready before calling spectrum
     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);
+    eq_noise.spectrum("displayname", options.displayname, "fignum", fig.Number, "normalizeTo0dB", 0,"color",options.color,"fft_length",4096);
 
     if ~isnan(options.postfilter_taps)
         % Hold on to the figure for further plotting
@@ -44,7 +50,88 @@ end
         legend('show');
     end
 
+    if useAutomaticColor
+        nextGroupIdx = nextGroupIdx + 1;
+        newLines = getNewLines(ax, linesBefore);
+        tagEQNoisePSDLines(newLines, nextGroupIdx);
+        recolorEQNoisePSDLines(ax, options.colormode);
+    end
+
     xlim([-eq_noise.fs/2* 1e-9 eq_noise.fs/2* 1e-9]);
     % ylim([-15, 0]);
 
 end
+
+function color = selectEQNoisePSDColor(ax, colormode)
+if colormode == "qualitative"
+    N = sum(arrayfun(@(x) strcmp(x.LineStyle, '-'), ax.Children));
+    cmap = linspecer(8);
+    color = cmap(mod(N, size(cmap, 1)) + 1, :);
+else
+    color = [0 0 0];
+end
+end
+
+function newLines = getNewLines(ax, linesBefore)
+linesAfter = findall(ax, 'Type', 'Line');
+newLines = linesAfter(~ismember(linesAfter, linesBefore));
+end
+
+function tagEQNoisePSDLines(lines, groupIdx)
+for lineIdx = 1:numel(lines)
+    setappdata(lines(lineIdx), 'showEQNoisePSDGroupIdx', groupIdx);
+end
+end
+
+function recolorEQNoisePSDLines(ax, colormode)
+if colormode == "qualitative"
+    return
+end
+
+lines = findall(ax, 'Type', 'Line');
+groupIdx = NaN(size(lines));
+
+for lineIdx = 1:numel(lines)
+    if isappdata(lines(lineIdx), 'showEQNoisePSDGroupIdx')
+        groupIdx(lineIdx) = getappdata(lines(lineIdx), 'showEQNoisePSDGroupIdx');
+    end
+end
+
+lines = lines(~isnan(groupIdx));
+groupIdx = groupIdx(~isnan(groupIdx));
+
+if isempty(lines)
+    return
+end
+
+groups = unique(groupIdx, 'stable');
+groups = sort(groups);
+cols = getEQNoisePSDColormap(colormode, numel(groups));
+
+for groupColorIdx = 1:numel(groups)
+    sameGroup = groupIdx == groups(groupColorIdx);
+    set(lines(sameGroup), 'Color', cols(groupColorIdx, :));
+end
+end
+
+function cols = getEQNoisePSDColormap(colormode, numColors)
+switch colormode
+    case "paired"
+        cols = cbrewer2('paired', numColors);
+    otherwise
+        cols = getDivergingWithoutBrightCenter(numColors);
+end
+end
+
+function cols = getDivergingWithoutBrightCenter(numColors)
+numBaseColors = max(6, numColors);
+centerExcludeFraction = 0.3;
+
+baseCols = cbrewer2('RdYlBu', numBaseColors);
+centerIdx = (numBaseColors + 1) / 2;
+excludeHalfWidth = centerExcludeFraction * numBaseColors / 2;
+keepIdx = find(abs((1:numBaseColors) - centerIdx) > excludeHalfWidth);
+sampleIdx = round(linspace(1, numel(keepIdx), numColors));
+
+cols = baseCols(keepIdx(sampleIdx), :);
+end

projects/Diss/MPI_revisit/algorithms/PLOT_mpi_reduction_db_delay_regimes_pam_comparison.m

@@ -0,0 +1,514 @@
+%% BER over SIR by MPI delay/coherence regime and PAM level
+% 1) gather block_update = 1 data from the database
+% 2) clean data and assign path-length regimes
+% 3) plot PAM4/PAM6/PAM8 BER-over-SIR curves in the same figure
+
+clear; clc;
+
+%% 1) Gather data
+
+studyName = "block_update_sweep";
+selectedBlockUpdate = 1;
+selectedPamLevels = [4 6 8];
+
+algorithmSelection = table( ...
+    ["dc_tracking"; ...
+     "plain_ffe"; ...
+     "a2_tracked_levels"; ...
+     "a2_residual"; ...
+     "a1_moving_average"], ...
+    [true; true; false; false; false], ...
+    'VariableNames', ["algorithm", "enabled"]);
+selectedAlgorithms = algorithmSelection.algorithm(algorithmSelection.enabled);
+
+useBoundedLines = true;
+usePolyfit = true;
+polyfitOrderMax = 4;
+boundaryPolyfitOrderMax = 4;
+maxBerForPlot = 0.1;
+scatterKeepFraction = 0.2;
+boundMode = "fitStd";          % "fitStd", "directStd", "fitCI", or "directCI"
+confidenceLevel = 0.95;
+regimeNames = ["0-1 m", "10-100 m", "300 m", "1000 m"];
+
+db = DBHandler( ...
+    "dataBase", "labor", ...
+    "type", "mysql", ...
+    "server", "192.168.178.192", ...
+    "user", "silas", ...
+    "password", "silas");
+db.refresh();
+
+fp = QueryFilter();
+fp.where('Runs', 'fiber_length', 'EQUALS', 0);
+fp.where('Runs', 'db_mode', 'EQUALS', '"no_db"');
+fp.where('Runs', 'v_bias', 'EQUALS', 2.65);
+fp.where('MpiReductionResults', 'study_name', 'EQUALS', char(studyName));
+fp.where('MpiReductionResults', 'block_update', 'EQUALS', selectedBlockUpdate);
+
+selectedFields = { ...
+    'Runs.run_id', ...
+    'Runs.sir', ...
+    'Runs.pam_level', ...
+    'Runs.symbolrate', ...
+    'Runs.interference_path_length', ...
+    'Runs.power_mpi_interference', ...
+    'Runs.power_pd_in', ...
+    'Runs.v_bias', ...
+    'Runs.loop_id', ...
+    'MpiReductionResults.occurrence_idx', ...
+    'MpiReductionResults.storage_name', ...
+    'MpiReductionResults.algorithm', ...
+    'MpiReductionResults.algorithm_variant', ...
+    'MpiReductionResults.eq_class', ...
+    'MpiReductionResults.block_update', ...
+    'MpiReductionResults.BER', ...
+    'MpiReductionResults.BER_precoded', ...
+    'MpiReductionResults.SNR', ...
+    'MpiReductionResults.GMI', ...
+    'MpiReductionResults.AIR'};
+selectedFields = selectedFields(:);
+
+[rawData, query] = db.queryDB(fp, selectedFields);
+disp(query);
+fprintf("Fetched %d MPI reduction rows.\n", height(rawData));
+
+%% 2) Clean data and assign regimes
+
+data = rawData;
+numericFields = ["run_id", "sir", "pam_level", "symbolrate", ...
+    "interference_path_length", "occurrence_idx", "block_update", ...
+    "power_mpi_interference", "BER", "BER_precoded", "SNR", "GMI", "AIR"];
+for fieldIdx = 1:numel(numericFields)
+    fieldName = numericFields(fieldIdx);
+    if ismember(fieldName, string(data.Properties.VariableNames))
+        data.(char(fieldName)) = numericColumn(data.(char(fieldName)));
+    end
+end
+
+stringFields = ["storage_name", "algorithm", "algorithm_variant", "eq_class"];
+for fieldIdx = 1:numel(stringFields)
+    fieldName = stringFields(fieldIdx);
+    if ismember(fieldName, string(data.Properties.VariableNames))
+        data.(char(fieldName)) = stringColumn(data.(char(fieldName)));
+    end
+end
+
+data(data.run_id == 3866, :) = [];
+data(data.run_id == 3865, :) = [];
+data(data.run_id == 3796, :) = [];
+data(data.run_id == 3797, :) = [];
+data(data.run_id == 3798, :) = [];
+data(data.run_id == 4002, :) = [];
+data(data.run_id == 4200, :) = [];
+data(data.run_id == 4199, :) = [];
+data(data.loop_id == 91, :) = [];
+data(data.loop_id == 59, :) = [];
+
+data = data(isfinite(data.BER) & data.BER > 0 & data.BER < maxBerForPlot, :);
+data.sir_exact = -7 - data.power_mpi_interference;
+
+data.path_regime = strings(height(data), 1);
+data.path_regime(data.interference_path_length == 0 | ...
+    data.interference_path_length == 1) = "0-1 m";
+data.path_regime(data.interference_path_length >= 10 & ...
+    data.interference_path_length <= 100) = "10-100 m";
+data.path_regime(data.interference_path_length == 300) = "300 m";
+data.path_regime(data.interference_path_length == 1000) = "1000 m";
+data = data(data.path_regime ~= "", :);
+
+data = data(ismember(data.pam_level, selectedPamLevels), :);
+data = data(ismember(data.algorithm, selectedAlgorithms), :);
+selectedAlgorithms = selectedAlgorithms(ismember(selectedAlgorithms, unique(data.algorithm, "stable")));
+
+if isempty(data)
+    warning("PLOT_mpi_reduction_db_delay_regimes_pam_comparison:NoRows", ...
+        "No rows remain after BER/study/block/PAM/algorithm/regime filtering.");
+    return
+end
+
+data.clean_keep = true(height(data), 1);
+groupId = findgroups(data.path_regime, data.pam_level, data.algorithm, data.sir);
+for curGroup = unique(groupId(isfinite(groupId))).'
+    rowMask = groupId == curGroup;
+    berValues = data.BER(rowMask);
+    if nnz(rowMask) > 3
+        data.clean_keep(rowMask) = ~isoutlier(berValues);
+    end
+end
+cleanData = data(data.clean_keep, :);
+
+groupVars = ["path_regime", "pam_level", "algorithm", "sir"];
+summaryTable = groupsummary(cleanData, groupVars, ...
+    {"mean", "min", "max"}, "BER");
+sirExactTable = groupsummary(cleanData, groupVars, "median", "sir_exact");
+summaryTable = sortrows(summaryTable, groupVars);
+sirExactTable = sortrows(sirExactTable, groupVars);
+summaryTable.sir_exact = sirExactTable.median_sir_exact;
+summaryTable = addBerIntervalBounds(summaryTable, cleanData, groupVars, confidenceLevel);
+
+fprintf("Cleaned to %d rows across %d regime/PAM/algorithm/SIR groups.\n", ...
+    height(cleanData), height(summaryTable));
+disp(groupcounts(cleanData, ["path_regime", "pam_level", "algorithm"]));
+
+%% 3) Plot PAM curves in one delay-regime comparison figure
+
+
+
+for regimeIdx = 1:numel(regimeNames)
+    regimeName = regimeNames(regimeIdx);
+    figure(); hold on;
+
+    for algIdx = 1:numel(selectedAlgorithms)
+        algorithmName = selectedAlgorithms(algIdx);
+        algColor = algorithmColor(algorithmName);
+        displayName = algorithmDisplayName(algorithmName);
+
+        for pamIdx = 1:numel(selectedPamLevels)
+            pamLevel = selectedPamLevels(pamIdx);
+            [marker, lineStyle] = pamStyle(pamLevel);
+
+            rawMask = cleanData.path_regime == regimeName & ...
+                cleanData.pam_level == pamLevel & ...
+                cleanData.algorithm == algorithmName;
+            curveMask = summaryTable.path_regime == regimeName & ...
+                summaryTable.pam_level == pamLevel & ...
+                summaryTable.algorithm == algorithmName;
+
+            if ~any(curveMask)
+                continue
+            end
+
+            scatterRows = downsampleRows(cleanData(rawMask, :), scatterKeepFraction);
+            % scatter(scatterRows.sir_exact, scatterRows.BER, ...
+            %     14, ...
+            %     "Marker", marker, ...
+            %     "MarkerEdgeColor", algColor, ...
+            %     "MarkerFaceColor", algColor, ...
+            %     "MarkerEdgeAlpha", 0.35, ...
+            %     "MarkerFaceAlpha", 0.35, ...
+            %     "HandleVisibility", "off");
+
+            sirValues = summaryTable.sir_exact(curveMask).';
+            meanBer = summaryTable.mean_BER(curveMask).';
+            [boundCenterBer, boundLowerBer, boundUpperBer] = ...
+                selectBerBounds(summaryTable(curveMask, :), boundMode);
+            valid = isfinite(sirValues) & isfinite(meanBer) & meanBer > 0;
+
+            if useBoundedLines && exist("boundedline", "file") && any(valid)
+                [xBand, centerBand, yBounds] = berIntervalBounds( ...
+                    sirValues, boundCenterBer, boundLowerBer, boundUpperBer, ...
+                    boundMode, boundaryPolyfitOrderMax);
+
+                [hl, hp] = boundedline(xBand, centerBand, yBounds, ...
+                    'alpha', 'transparency', 0.08, ...
+                    'cmap', algColor, ...
+                    'nan', 'fill', ...
+                    'orientation', 'vert');
+                set(hl, "LineStyle", "none", "Marker", "none", "HandleVisibility", "off");
+                set(hp, "LineStyle", "-", "HandleVisibility", "off", "Marker", "none");
+            end
+
+            plot(sirValues(valid), meanBer(valid), ...
+                "LineStyle", lineStyle, ...
+                "Marker", marker, ...
+                "MarkerSize", 3.5, ...
+                "LineWidth", 1, ...
+                "Color", algColor, ...
+                "MarkerFaceColor", "w", ...
+                "MarkerEdgeColor", algColor, ...
+                "DisplayName", sprintf("%s, PAM %.0f", displayName, pamLevel));
+
+            if usePolyfit
+                fitMask = valid & meanBer > 0;
+                if nnz(fitMask) >= 2
+                    fitOrder = min(polyfitOrderMax, nnz(fitMask) - 1);
+                    fitCoeff = polyfit(sirValues(fitMask), log10(meanBer(fitMask)), fitOrder);
+                    xFit = linspace(min(sirValues(fitMask)), max(sirValues(fitMask)), 300);
+                    yFit = 10 .^ polyval(fitCoeff, xFit);
+
+                    % plot(xFit, yFit, ...
+                    %     "LineStyle", "--", ...
+                    %     "LineWidth", 1.1, ...
+                    %     "Color", algColor, ...
+                    %     "HandleVisibility", "off");
+                end
+            end
+        end
+    end
+
+    yline(2.2e-4, "LineWidth", 1, "LineStyle", "--", "HandleVisibility", "off");
+    yline(3.8e-3, "LineWidth", 1, "LineStyle", "--", "HandleVisibility", "off");
+    yline(2e-2, "LineWidth", 1, "LineStyle", "--", "HandleVisibility", "off");
+
+    title(sprintf("Interference path regime %s", regimeName));
+    xlabel("SIR (dB)");
+    ylabel("BER");
+    set(gca, "YScale", "log");
+    ylim([9e-5, maxBerForPlot]);
+    xlim([15, 45]);
+    grid on;
+    box on;
+    legend("Location", "northeast", "Interpreter", "none");
+
+    if exist("beautifyBERplot", "file")
+        beautifyBERplot("logscale", true, "setcolors", false, ...
+            "setmarkers", false, "changemarkers", false);
+    end
+end
+
+% mat2tikz_improved("C:/Users/Silas/Documents/6971e0b65b380ca6d71c837f/04_Experimental_Evaluation/tikz/mpi/ber_vs_sir_delay_regimes_pam_comparison.tikz","cleanfigure",1);
+
+%% Local helpers
+
+function values = numericColumn(values)
+    if iscell(values)
+        values = string(values);
+    end
+    if isstring(values) || ischar(values)
+        values = str2double(values);
+    end
+    values = double(values);
+end
+
+function values = stringColumn(values)
+    if iscell(values)
+        values = string(values);
+    elseif ischar(values)
+        values = string(values);
+    end
+    values = strip(string(values));
+end
+
+function rows = downsampleRows(rows, keepFraction)
+    if keepFraction >= 1 || height(rows) <= 1
+        return
+    end
+
+    keepEvery = max(1, round(1 / keepFraction));
+    keepIdx = 1:keepEvery:height(rows);
+    rows = rows(keepIdx, :);
+end
+
+function summaryTable = addBerIntervalBounds(summaryTable, cleanData, groupVars, confidenceLevel)
+    nGroups = height(summaryTable);
+    ciCenter = NaN(nGroups, 1);
+    ciLower = NaN(nGroups, 1);
+    ciUpper = NaN(nGroups, 1);
+    stdCenter = NaN(nGroups, 1);
+    stdLower = NaN(nGroups, 1);
+    stdUpper = NaN(nGroups, 1);
+
+    for groupIdx = 1:nGroups
+        rowMask = true(height(cleanData), 1);
+        for varIdx = 1:numel(groupVars)
+            varName = groupVars(varIdx);
+            rowMask = rowMask & cleanData.(char(varName)) == summaryTable.(char(varName))(groupIdx);
+        end
+
+        [ciCenter(groupIdx), ciLower(groupIdx), ciUpper(groupIdx)] = ...
+            logBerMeanConfidenceInterval(cleanData.BER(rowMask), confidenceLevel);
+        [stdCenter(groupIdx), stdLower(groupIdx), stdUpper(groupIdx)] = ...
+            logBerMeanStdInterval(cleanData.BER(rowMask));
+    end
+
+    summaryTable.ci_center_BER = ciCenter;
+    summaryTable.ci_lower_BER = ciLower;
+    summaryTable.ci_upper_BER = ciUpper;
+    summaryTable.std_center_BER = stdCenter;
+    summaryTable.std_lower_BER = stdLower;
+    summaryTable.std_upper_BER = stdUpper;
+end
+
+function [centerBer, lowerBer, upperBer] = logBerMeanConfidenceInterval(berValues, confidenceLevel)
+    berValues = berValues(isfinite(berValues) & berValues > 0);
+    if isempty(berValues)
+        centerBer = NaN;
+        lowerBer = NaN;
+        upperBer = NaN;
+        return
+    end
+
+    logBer = log10(berValues(:));
+    centerLog = mean(logBer, "omitnan");
+
+    if numel(logBer) < 2
+        centerBer = 10 .^ centerLog;
+        lowerBer = centerBer;
+        upperBer = centerBer;
+        return
+    end
+
+    alpha = 1 - confidenceLevel;
+    if exist("tinv", "file")
+        tCritical = tinv(1 - alpha/2, numel(logBer) - 1);
+    else
+        tCritical = 1.96;
+    end
+
+    halfWidth = tCritical * std(logBer, 0, "omitnan") / sqrt(numel(logBer));
+    centerBer = 10 .^ centerLog;
+    lowerBer = 10 .^ (centerLog - halfWidth);
+    upperBer = 10 .^ (centerLog + halfWidth);
+end
+
+function [centerBer, lowerBer, upperBer] = logBerMeanStdInterval(berValues)
+    berValues = berValues(isfinite(berValues) & berValues > 0);
+    if isempty(berValues)
+        centerBer = NaN;
+        lowerBer = NaN;
+        upperBer = NaN;
+        return
+    end
+
+    logBer = log10(berValues(:));
+    centerLog = mean(logBer, "omitnan");
+    stdLog = std(logBer, 0, "omitnan");
+
+    centerBer = 10 .^ centerLog;
+    lowerBer = 10 .^ (centerLog - stdLog);
+    upperBer = 10 .^ (centerLog + stdLog);
+end
+
+function [centerBer, lowerBer, upperBer] = selectBerBounds(curveSummary, boundMode)
+    if boundMode == "fitCI" || boundMode == "directCI"
+        centerBer = curveSummary.ci_center_BER.';
+        lowerBer = curveSummary.ci_lower_BER.';
+        upperBer = curveSummary.ci_upper_BER.';
+    else
+        centerBer = curveSummary.std_center_BER.';
+        lowerBer = curveSummary.std_lower_BER.';
+        upperBer = curveSummary.std_upper_BER.';
+    end
+end
+
+function [xBand, centerBand, yBounds] = berIntervalBounds(sirValues, centerBer, lowerBer, upperBer, boundMode, maxOrder)
+    if boundMode == "directCI" || boundMode == "directStd"
+        [xBand, centerBand, yBounds] = directBerBounds(sirValues, centerBer, lowerBer, upperBer);
+        return
+    end
+
+    [xBand, centerBand, yBounds] = fittedBerBounds(sirValues, centerBer, lowerBer, upperBer, maxOrder);
+end
+
+function [xBand, centerBand, yBounds] = directBerBounds(sirValues, centerBer, lowerBer, upperBer)
+    valid = isfinite(sirValues) & isfinite(centerBer) & isfinite(lowerBer) & ...
+        isfinite(upperBer) & centerBer > 0 & lowerBer > 0 & upperBer > 0;
+
+    xBand = sirValues(valid).';
+    centerBand = centerBer(valid).';
+    lowerBand = lowerBer(valid).';
+    upperBand = upperBer(valid).';
+
+    [xBand, orderIdx] = sort(xBand(:));
+    centerBand = centerBand(orderIdx);
+    lowerBand = lowerBand(orderIdx);
+    upperBand = upperBand(orderIdx);
+
+    lowerTmp = min(lowerBand, upperBand);
+    upperBand = max(lowerBand, upperBand);
+    lowerBand = lowerTmp;
+
+    centerBand = min(max(centerBand, lowerBand), upperBand);
+    yBounds = [max(centerBand - lowerBand, 0), max(upperBand - centerBand, 0)];
+end
+
+function [xBand, centerBand, yBounds] = fittedBerBounds(sirValues, meanBer, minBer, maxBer, maxOrder)
+    valid = isfinite(sirValues) & isfinite(meanBer) & isfinite(minBer) & ...
+        isfinite(maxBer) & meanBer > 0 & minBer > 0 & maxBer > 0;
+
+    x = sirValues(valid);
+    yMean = meanBer(valid);
+    yMin = minBer(valid);
+    yMax = maxBer(valid);
+
+    if numel(x) < 2
+        xBand = x(:);
+        centerBand = yMean(:);
+        yBounds = [max(yMean(:) - yMin(:), 0), max(yMax(:) - yMean(:), 0)];
+        return
+    end
+
+    [x, orderIdx] = sort(x(:));
+    yMean = yMean(orderIdx);
+    yMin = yMin(orderIdx);
+    yMax = yMax(orderIdx);
+
+    xBand = linspace(min(x), max(x), 300).';
+    fitOrder = min(maxOrder, numel(unique(x)) - 1);
+
+    if fitOrder < 1
+        centerBand = interp1(x, yMean, xBand, "linear", "extrap");
+        lowerBand = interp1(x, yMin, xBand, "linear", "extrap");
+        upperBand = interp1(x, yMax, xBand, "linear", "extrap");
+    else
+        centerBand = fitLogBer(x, yMean, xBand, fitOrder);
+        lowerBand = fitLogBer(x, yMin, xBand, fitOrder);
+        upperBand = fitLogBer(x, yMax, xBand, fitOrder);
+    end
+
+    lowerTmp = min(lowerBand, upperBand);
+    upperBand = max(lowerBand, upperBand);
+    lowerBand = lowerTmp;
+
+    centerBand = min(max(centerBand, lowerBand), upperBand);
+    yBounds = [max(centerBand - lowerBand, 0), max(upperBand - centerBand, 0)];
+end
+
+function yFit = fitLogBer(x, y, xFit, fitOrder)
+    coeff = polyfit(x, log10(y), fitOrder);
+    yFit = 10 .^ polyval(coeff, xFit);
+end
+
+function label = algorithmDisplayName(algorithmName)
+    algorithmName = string(algorithmName);
+    switch algorithmName
+        case "plain_ffe"
+            label = "FFE only";
+        case "a2_tracked_levels"
+            label = "ACT";
+        case "a2_residual"
+            label = "L-DCA";
+        case "a1_moving_average"
+            label = "DCA";
+        case "dc_tracking"
+            label = "DCT";
+        otherwise
+            label = algorithmName;
+    end
+end
+
+function color = algorithmColor(algorithmName)
+    algorithmName = string(algorithmName);
+    switch algorithmName
+        case "plain_ffe"
+            color = [0.3467 0.5360 0.6907];
+        case "a2_tracked_levels"
+            color = [0.9153 0.2816 0.2878];
+        case "a2_residual"
+            color = [0.4416 0.7490 0.4322];
+        case "a1_moving_average"
+            color = [1.0000 0.5984 0.2000];
+        case "dc_tracking"
+            color = [0.6769 0.4447 0.7114];
+        otherwise
+            color = [0 0 0];
+    end
+end
+
+function [marker, lineStyle] = pamStyle(pamLevel)
+    switch pamLevel
+        case 4
+            lineStyle = "-";
+            marker = "o";
+        case 6
+            lineStyle = "--";
+            marker = "square";
+        case 8
+            lineStyle = ":";
+            marker = "diamond";
+        otherwise
+            lineStyle = "-.";
+            marker = "^";
+    end
+end

projects/Diss/MPI_revisit/investigate_mpi_algorithms.m

@@ -3,7 +3,7 @@ dsp_options = struct();
 dsp_options.mode = "run_id";
 dsp_options.recipe = @mpi_recipe_dev;
 dsp_options.append_to_db = false;
-dsp_options.start_occurence = 5;
+dsp_options.start_occurence = 1;
 dsp_options.max_occurences = 15;
 dsp_options.debug_plots = false;
 
@@ -51,7 +51,7 @@ for i = 1
     fp.where('Runs', 'symbolrate',                'EQUALS', B); % 72 96 112
     fp.where('Runs', 'fiber_length',              'EQUALS', 0);
     fp.where('Runs', 'interference_path_length',  'EQUALS', 1000);
-    fp.where('Runs', 'sir',                       'LESS_EQUAL', 20);
+    fp.where('Runs', 'sir',                       'LESS_EQUAL', 50);
     fp.where('Runs', 'db_mode',                   'EQUALS', '"no_db"');
     % fp.where('Runs', 'is_mpi',                    'EQUALS', 1);
     fp.where('Runs', 'pam_level',                 'EQUALS', M);
@@ -62,14 +62,20 @@ for i = 1
     [~, sortIdx] = sort(dataTable.sir, 'descend');
     dataTable = dataTable(sortIdx, :);
 
-    dataTable = dataTable(1,:);
+    % desired_sir = [15,22,26,30,38];
+    % desired_sir = [22,30,38];
+    % % Filter dataTable to only keep rows with sir in desired_sir
+    % isDesired = ismember(dataTable.sir, desired_sir);
+    % dataTable = dataTable(isDesired, :);
+    % dataTable = dataTable(1,:);
 
     run_ids = dataTable.run_id;
 
     %% === Warehouse setup ===
     dsp_options.userParameters = struct();
     dsp_options.userParameters.block_update = 1;%[1,2,4,8,16,32,64,112,224,448,448*2,1024,2048,4096,8192,16384];%%logspace(-3.8,-1,22);%[linspace(2,4096,22)];
-    dsp_options.userParameters.hpf = [2e6:2e6:10e6,20e6:10e6:100e6];
+    % dsp_options.userParameters.hpf = [2e6:0.5e6:10e6];
+    % dsp_options.userParameters.bufferlen = [112,224,448,448*2,1024,2048,4096,8192,16384];
     wh = DataStorage(dsp_options.userParameters);
     
     %%
@@ -85,7 +91,7 @@ for i = 1
     %% === 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);
 
@@ -109,9 +115,9 @@ storageNames = fieldnames(wh.sto);
 
 x_base = dataTable.sir(:).';
 x_base = dsp_options.userParameters.block_update;
-x_base = dsp_options.userParameters.hpf;
+x_base = dsp_options.userParameters.bufferlen;
-
+% x_base = dsp_options.userParameters.hpf;
-figure(2026); clf; hold on
+figure(2029); hold on
 for storage_idx = 1:numel(storageNames)
     storageName = storageNames{storage_idx};
     result = wh.sto.(storageName);

projects/Diss/MPI_revisit/parallelization_analysis/PLOT_mpi_reduction_db_parallelization_vs_sir.m

@@ -8,7 +8,7 @@ clear; clc;
 
 studyName = "block_update_sweep";
 pathLengthToPlot = 1000;
-selectedPamLevels = 6;      % set [] to use all PAM levels in the query result
+selectedPamLevels = 4;      % set [] to use all PAM levels in the query result
 selectedAlgorithms = [];    % set [] to use all algorithms in the query result
 
 useBoundedLines = true;