diff --git a/Classes/04_DSP/Equalizer/FFE.m b/Classes/04_DSP/Equalizer/FFE.m index 9ac0ebd..d38c425 100644 --- a/Classes/04_DSP/Equalizer/FFE.m +++ b/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; diff --git a/Functions/EQ_blocks/preprocessSignal.m b/Functions/EQ_blocks/preprocessSignal.m index 2161b2f..679fbb0 100644 --- a/Functions/EQ_blocks/preprocessSignal.m +++ b/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 diff --git a/Functions/EQ_recipes/mpi_recipe_dev.m b/Functions/EQ_recipes/mpi_recipe_dev.m index 5e5abfa..0321960 100644 --- a/Functions/EQ_recipes/mpi_recipe_dev.m +++ b/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 - dc_buffer_len = 0; - ffe_buffer_len = 0; - smoothing_buffer_length = options.userParameters.smoothing_length; - smoothing_buffer_update = 1; + mu_dc = 0.005;%1e-5; eq_settings = { ... "epochs_tr", 5, ... - "epochs_dd", 5, ... + "epochs_dd", 3, ... "len_tr", 4096*2, ... - "mu_dd", 1e-1, ... - "mu_tr", 0.4, ... + "mu_dd",3.005e-03, ... + "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, ... - % "fsym", options.fsym, ... - % "fignum", options.dataTable.run_id, ... - % "normalize", true); + showLevelScatter(Scpe_sig_raw, Symbols, ... + "fsym", options.fsym, ... + "fignum", 400, ... + "normalize", true); - %% - % 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; + % Scpe_sig_raw.spectrum("normalizeTo0dB",1,"fft_length",4096*4,"fignum",401); - %% - 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); - 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_dcr.metrics.print("description",'FFE DCR'); - output.ffe_dcr_package = ffe_results_dcr; + %options.dataTable.sir; + %% NORMAL FFE + if 1 + + 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); + + 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 diff --git a/Functions/EQ_visuals/showEQNoisePSD.m b/Functions/EQ_visuals/showEQNoisePSD.m index bfb0a12..1a16831 100644 --- a/Functions/EQ_visuals/showEQNoisePSD.m +++ b/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 diff --git a/Functions/EQ_visuals/showLevelHistogram.m b/Functions/EQ_visuals/showLevelHistogram.m index 6ef093d..c86828d 100644 --- a/Functions/EQ_visuals/showLevelHistogram.m +++ b/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); diff --git a/Functions/EQ_visuals/showLevelScatter.m b/Functions/EQ_visuals/showLevelScatter.m index 1565487..7b3b62d 100644 --- a/Functions/EQ_visuals/showLevelScatter.m +++ b/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 - scatter(xAxisUs, symbols_for_lvl(levelIdx, :), 10, ".", ... - "MarkerEdgeColor", cols((2*levelIdx)-1, :), ... - "MarkerEdgeAlpha", 0.5); +if 1 + for levelIdx = numLevels:-1:1 + scatter(ax, xAxisUs, symbols_for_lvl(levelIdx, :), options.scatterSize, ".", ... + "MarkerEdgeColor", cols(levelIdx, :), ... + "MarkerEdgeAlpha", options.scatterAlpha); + end end for levelIdx = 1:numLevels - plot(xAxisUs, avg_for_lvl(levelIdx, :), ... - "LineWidth", 1, ... - "Color", cols(2*levelIdx, :)); + [xReduced, yReduced] = reduceAverageLineForPlot( ... + xAxisUs, avg_for_lvl(levelIdx, :), ... + options.avgLineMaxPoints, options.avgLineSmoothWindow); + plot(ax, xReduced, yReduced, ... + "LineWidth", 2, ... + "Color", cols(levelIdx, :)); end -refValues = numericSignal(refSymbols); -yline(unique(refValues), "HandleVisibility", "off"); +% refValues = unique(numericSignal(refSymbols)); +% 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'); -ylabel('Normalized Amplitude'); -xlim(getXLimits(options.xLimits, xAxisUs, symbols_for_lvl)); -ylim(options.yLimits); -grid on +xlabel(ax, 'Time in $\mu$s'); +ylabel(ax, 'Normalized Amplitude'); +xlim(ax, xLimits); +ylim(ax, options.yLimits); +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, ... diff --git a/Functions/Job_Processing/runBatch.m b/Functions/Job_Processing/runBatch.m index b37d675..e091f99 100644 --- a/Functions/Job_Processing/runBatch.m +++ b/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,28 +38,45 @@ 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; - batchResults{jobIdx} = jobResult; - fprintf('[%s] Completed (%d/%d).\n', jobs(jobIdx).label, completedCount, nJobs); + jobDurations(jobIdx) = jobRuntime; + lastJobDuration = jobRuntime; - if ~isempty(options.resultHandler) - options.resultHandler(jobResult, jobs(jobIdx), jobIdx); + if isempty(jobError) + batchResults{jobIdx} = jobResult; + 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); @@ -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, ... - sprintf('Completed %d/%d jobs', completedCount, totalCount)); + completedDurations = jobDurations(~isnan(jobDurations)); + 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); diff --git a/Functions/beautifyBERplot.m b/Functions/beautifyBERplot.m index 8857161..5f84673 100644 --- a/Functions/beautifyBERplot.m +++ b/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); -if n == 0 - return -end +dataObjects = children((isLine | isScatter) & ~isConstantLine); +n = numel(dataObjects); % --------------------------------------------------------- % Fixed color palette (deterministic across figures) % --------------------------------------------------------- -try - cmap = linspecer(n); -catch - cmap = lines(n).Color; %#ok - cmap = linespecer_fallback(n); +if n > 0 + try + cmap = linspecer(n); + catch + 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 - ln.LineWidth = lw; - - if options.setcolors - ln.Color = cmap(i,:); - end + dataObj = dataObjects(n-i+1); % preserve plotting order + + if isa(dataObj,'matlab.graphics.chart.primitive.Line') + dataObj.LineWidth = lw; - if options.setmarkers - if strcmp(ln.Marker,'none') - ln.Marker = markers{mod(i-1,numel(markers))+1}; - end - ln.MarkerSize = ms; if options.setcolors - ln.MarkerEdgeColor = cmap(i,:); + dataObj.Color = cmap(i,:); + end + + if options.setmarkers + if strcmp(dataObj.Marker,'none') + dataObj.Marker = markers{mod(i-1,numel(markers))+1}; + end + dataObj.MarkerSize = ms; + if options.setcolors + 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); - x = ln.XData(:); - y = ln.YData(:); + dataObj = dataObjects(n-i+1); + x = dataObj.XData(:); + 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,... diff --git a/Functions/mat2tikz_improved.m b/Functions/mat2tikz_improved.m index 345f057..0ba1c6d 100644 --- a/Functions/mat2tikz_improved.m +++ b/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', ... @@ -38,4 +42,4 @@ matlab2tikz(char(filename), ... % 'tick label style={/pgf/number format/fixed, /pgf/number format/1000 sep={}}' ... % }); -end \ No newline at end of file +end diff --git a/projects/Diss/MPI_revisit/algorithms/ffe_baseline.fig b/projects/Diss/MPI_revisit/algorithms/ffe_baseline.fig new file mode 100644 index 0000000..b0f2e06 Binary files /dev/null and b/projects/Diss/MPI_revisit/algorithms/ffe_baseline.fig differ diff --git a/projects/Diss/MPI_revisit/algorithms/wh_ffe_baseline.mat b/projects/Diss/MPI_revisit/algorithms/wh_ffe_baseline.mat new file mode 100644 index 0000000..132cedc Binary files /dev/null and b/projects/Diss/MPI_revisit/algorithms/wh_ffe_baseline.mat differ diff --git a/projects/Diss/MPI_revisit/analytic_mpi_evaluation_with_levels.m b/projects/Diss/MPI_revisit/analytic_mpi_evaluation_with_levels.m new file mode 100644 index 0000000..df09e81 --- /dev/null +++ b/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"); \ No newline at end of file diff --git a/projects/Diss/MPI_revisit/calculate_and_store_std.m b/projects/Diss/MPI_revisit/calculate_and_store_std.m new file mode 100644 index 0000000..f86efb6 --- /dev/null +++ b/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 + diff --git a/projects/Diss/MPI_revisit/investigate_database.m b/projects/Diss/MPI_revisit/investigate_database.m new file mode 100644 index 0000000..a39b629 --- /dev/null +++ b/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'); diff --git a/projects/Diss/MPI_revisit/investigate_mpi_algorithms.m b/projects/Diss/MPI_revisit/investigate_mpi_algorithms.m index 04fd076..ca17bf2 100644 --- a/projects/Diss/MPI_revisit/investigate_mpi_algorithms.m +++ b/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', 'db_mode', 'EQUALS', '"no_db"'); +fp.where('Runs', 'sir', 'LESS_EQUAL', 30); +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"); -dataTable = dataTable(uniqueSirRows, :); -% sort rows by sir (ascending) and extract run_ids +% [~, uniqueSirRows] = unique(dataTable.sir, "stable"); +% dataTable = dataTable(uniqueSirRows, :); + [~, 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); - ber = cellfun(@(packageCell) packageCell{1}.metrics.BER, result); - plot(x_vars,ber); - beautifyBERplot("logscale",true,"setcolors",true,"setmarkers",true); + result = wh.getStoValue(storageNames{1}, x_vars); + % Each cell in result is a cell array (packageCell) containing multiple packages. + % Extract BERs from all packages and, for example, average them per x_var. + 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 + diff --git a/projects/Diss/MPI_revisit/investigate_results.m b/projects/Diss/MPI_revisit/investigate_results.m new file mode 100644 index 0000000..d7f01d1 --- /dev/null +++ b/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); diff --git a/projects/Diss/MPI_revisit/plot_mpi_timesignal.m b/projects/Diss/MPI_revisit/plot_mpi_timesignal.m new file mode 100644 index 0000000..ec3462c --- /dev/null +++ b/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_m__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 + diff --git a/projects/Diss/MPI_revisit/verify_mpi_noise_power.m b/projects/Diss/MPI_revisit/verify_mpi_noise_power.m new file mode 100644 index 0000000..44598b5 --- /dev/null +++ b/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 + diff --git a/projects/ECOC_2025_MPI/load_signal_standalone.m b/projects/ECOC_2025_MPI/load_signal_standalone.m index b27238a..bf98730 100644 --- a/projects/ECOC_2025_MPI/load_signal_standalone.m +++ b/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); - var(sep_sig,0,2,'omitnan') + "debug_plots", false,"showPlot",true,"showStdAnnotations",false,"yLimits",[-3 3]); + + 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 diff --git a/projects/ECOC_2025_MPI/plots_from_database/bias_vs_ber.m b/projects/ECOC_2025_MPI/plots_from_database/bias_vs_ber.m index 8941c82..06a8b70 100644 --- a/projects/ECOC_2025_MPI/plots_from_database/bias_vs_ber.m +++ b/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.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, :); -dataTable = dataTable(dataTable.date_of_processing < stopTime, :); + +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 + +% 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 = 'SNR'; @@ -52,6 +77,7 @@ dataTableGrpd = groupIt(fixedVars,dataTable); plotRealizations = 1; +%% % Create a new figure mkr = 'x'; diff --git a/projects/ECOC_2025_MPI/plots_from_database/plot_bias_awg_vs_ber.m b/projects/ECOC_2025_MPI/plots_from_database/plot_bias_awg_vs_ber.m index 30b2657..01b1234 100644 --- a/projects/ECOC_2025_MPI/plots_from_database/plot_bias_awg_vs_ber.m +++ b/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'); -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_processing = datetime(dataTable.date_of_processing, 'InputFormat', 'yyyy-MM-dd HH:mm:ss.SSSSSS'); -dataTable = dataTable(dataTable.date_of_processing > startTime, :); -dataTable = dataTable(dataTable.date_of_processing < stopTime, :); + +%% Filter by time + +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'); + +%% +% 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 diff --git a/workerError.mat b/workerError.mat index a9bdfa3..d0ec048 100644 Binary files a/workerError.mat and b/workerError.mat differ