try rop = res.settings.rop; % 12 points wavelengthplan = res.settings.wavelengthplan; catch wavelengthplan = [1295,1305,1315,1325]; wavelengthplan = calcWavelengthPlan(16,400e9,1310); rop = -8.25:0.75:0; end N = length(wavelengthplan); figure(); hold on; cols = cbrewer2('set2',N); % one color per wavelength (Ch) fec = 2.2e-4; fec = 3.8e-3; Sffe = cell(1,N); Svnle = cell(1,N); Smlse = cell(1,N); Sdbt = cell(1,N); % Choose your quantile band. For your old style, use 0.04/0.99: qLow = 0.0; % lower quantile (e.g.,s 0.04 for old script) qHigh = 1; % upper quantile (e.g., 0.99 for old script) cols = linspecer(N); % one color per wavelength (Ch) cols = cbrewer2('set1',N); for l = 1:N % Slice 12x50 cell arrays ffe_cells = reshape(squeeze(res.ffe(l,:,:)),length(rop),[]); vnle_cells = reshape(squeeze(res.vnle(l,:,:)),length(rop),[]); mlse_cells = reshape(squeeze(res.mlse(l,:,:)),length(rop),[]); dbt_cells = reshape(squeeze(res.dbt(l,:,:)),length(rop),[]); [Sffe{l}, noX_ffe] = fecCrossings(rop, ffe_cells, fec); [Svnle{l}, noX_ffe] = fecCrossings(rop, vnle_cells, fec); [Smlse{l}, noX_ffe] = fecCrossings(rop, mlse_cells, fec); [Sdbt{l}, noX_ffe] = fecCrossings(rop, dbt_cells, fec); % Extract BER matrices using only complete realizations (12/12 ROP filled) ffe_mat = extractCompleteBER(ffe_cells); % 12 x K_ffe vnle_mat = extractCompleteBER(vnle_cells); % 12 x K_vnle mlse_mat = extractCompleteBER(mlse_cells); % 12 x K_mlse mlse_alpha_mat = extractCompleteAlphas(mlse_cells); % 12 x K_mlse dbt_mat = extractCompleteBER(dbt_cells); % 12 x K_dbt showLegend = 1; % one legend entry per technique % Plot shaded band + mean line with boundedline plotBandMeanBL(rop, ffe_mat, cols(l,:), sprintf('FFE @ %.1d nm',round(wavelengthplan(l))), qLow, qHigh, '--s', showLegend); % scatter(Sffe,fec.*ones(size(Sffe)),20,'v','MarkerFaceColor','black'); % plotBandMeanBL(rop, vnle_mat, cols(l,:), sprintf('VNLE @ %.1d nm',round(wavelengthplan(l))), qLow, qHigh, '--x', showLegend); % plotBandMeanBL(rop, mlse_mat, cols(l,:), sprintf('VNLE+PF+MLSE @ %.1d nm',round(wavelengthplan(l))), qLow, qHigh, '-o', showLegend); % plotBandMeanBL(rop, dbt_mat, cols(l,:), sprintf('DBt.+MLSE @ %.1d nm',round(wavelengthplan(l))), qLow, qHigh, '--v', showLegend); set(gca,'XScale','linear','YScale','log','TickLabelInterpreter','latex','FontSize',11); yline([3.8e-3, 2.2e-4], 'HandleVisibility','off','LineWidth',1.5); end ylabel('BER'); xlabel('ROP'); title('BER vs. ROP'); xlim([min(rop) max(rop)]); ylim([1e-5 0.3]); grid on; legend show; %% S_cell = Sdbt; S_cell =Smlse; S_cell = {Svnle,Smlse,Sdbt}; S_cell = {Svnle}; figure(5); hold on; for i = 1:length(S_cell) % Pad to rectangular matrix: rows = realizations, cols = wavelengths Kmax = max(cellfun(@numel, S_cell{i})); S_mat = NaN(Kmax, N); for l = 1:N k = numel(S_cell{i}{l}); if k > 0 S_mat(1:k, l) = S_cell{i}{l}; end end % --- Violin plot over wavelengths (columns) --- cols=linspecer(3); catLabels = arrayfun(@(nm) sprintf('%d nm', nm), wavelengthplan, 'UniformOutput', false); vs = violinplot(S_mat, catLabels, ... 'ViolinColor', cols(i,:), ... 'ViolinAlpha', 0.10, ... 'MarkerSize', 20, ... 'ShowMedian', true, ... 'EdgeColor', cols(i,:), ... 'ShowWhiskers', false, ... 'ShowData', true, ... 'ShowBox', false, ... 'Bandwidth', 0.05); ylim([floor(min(S_mat,[],'all')), ceil(max(S_mat,[],'all'))]) % ylim([-8 0]); ylabel('ROP at FEC crossing'); title(sprintf('RROP to cross BER %.2e', fec)); grid on; box on; end %% ================= helper ================= function plotBandMeanBL(x, Y, color, techLabel, qLow, qHigh, lineSpec, showLegend) % Y: (nPoints x nRealizations) % Remove realizations that are entirely zero (like removeZeros behavior) badCols = all(Y == 0, 1); Y(:, badCols) = []; Y(Y==0) = 1e-8; % Stats across realizations mu = mean(Y, 2, 'omitnan'); % mean line lo = quantile(Y, qLow, 2); % lower bound hi = quantile(Y, qHigh, 2); % upper bound % Convert to asymmetric distances required by boundedline: % b(:,1) = distance to lower side; b(:,2) = distance to upper side b = [mu - lo, hi - mu]; % Call boundedline with alpha shading [hl, hp] = boundedline(x(:), mu(:), b, lineSpec, 'alpha', ... 'transparency', 0.18); % Color styling set(hl, 'Color', color, 'LineWidth', 1.4, 'MarkerSize', 4); set(hp, 'FaceColor', color, 'HandleVisibility','off'); % patch hidden in legend % Single legend entry per technique (use first wavelength only) if showLegend set(hl, 'DisplayName', techLabel); else set(hl, 'HandleVisibility','off'); end % Optional: outline the bounds if outlinebounds is available if exist('outlinebounds','file') == 2 ho = outlinebounds(hl, hp); set(ho, 'linestyle', ':', 'color', color, 'linewidth', 1, ... 'HandleVisibility','off'); end end function [S, noCrossingMask, Y_keep] = fecCrossings(rop, cells12xR, fec) % cells12xR: 12xR cell array (one wavelength + scheme slice) % each cell must be a struct with .metrics.BER % rop: 12x1 numeric vector of ROP points % fec: scalar FEC threshold (e.g., 3.8e-3) % % Outputs: % S 1xK vector of crossing ROP per kept realization (NaN if none) % noCrossingMask 1xK logical mask: true if no crossing for that realization % Y_keep 12xK numeric BER matrix used for the crossing detection % 1) keep only complete realization columns Y = extractCompleteBER(cells12xR); % -> 12 x K if isempty(Y) S = []; noCrossingMask = []; Y_keep = Y; return; end % 2) optionally drop realizations with mean BER > 0.1 ok = mean(Y,1,'omitnan') <= 0.1; Y = Y(:, ok); if isempty(Y) S = []; noCrossingMask = []; Y_keep = Y; return; end % 3) find crossings per realization nR = size(Y,2); S = nan(1,nR); noCrossingMask = true(1,nR); rop = rop(:); % ensure column for j = 1:nR y = Y(:,j); % sign change from >fec to <=fec (first time it drops below FEC) above = (y > fec); idx = find(above(1:end-1) & ~above(2:end), 1, 'first'); if ~isempty(idx) % linear interpolation between (x1,y1) and (x2,y2) x1 = rop(idx); y1 = y(idx); x2 = rop(idx+1); y2 = y(idx+1); if isfinite(y1) && isfinite(y2) && y2 ~= y1 t = (fec - y1) / (y2 - y1); S(j) = x1 + t*(x2 - x1); noCrossingMask(j) = false; end end end Y_keep = Y; end function Y = extractCompleteBER(cellSlice) % cellSlice: 12xR cell array; each cell should be a struct with .metrics.BER % Keep only those realization columns where ALL 12 ROP entries are valid. if isempty(cellSlice), Y = []; return; end nR = size(cellSlice,2); keep = false(1,nR); for r = 1:nR col = cellSlice(:,r); keep(r) = all(cellfun(@(c) ~isempty(c) , col)); end if ~any(keep), Y = []; return; end Y = cellfun(@(c) c.metrics.BER, cellSlice(:,keep), 'UniformOutput', true); end function Y = extractCompleteAlphas(cellSlice) % cellSlice: 12xR cell array; each cell should be a struct with .metrics.BER % Keep only those realization columns where ALL 12 ROP entries are valid. if isempty(cellSlice), Y = []; return; end nR = size(cellSlice,2); keep = false(1,nR); for r = 1:nR col = cellSlice(:,r); keep(r) = all(cellfun(@(c) ~isempty(c) , col)); end if ~any(keep), Y = []; return; end Y = cellfun(@(c) c.metrics.Alpha, cellSlice(:,keep), 'UniformOutput', true); end