silas/imdd_silas
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

chirped WDM analsyis

Browse Source 
minor stuff in FSO link with Magnus
This commit is contained in:
silas (home)
2026-01-12 08:12:42 +01:00
parent 062c5dfbaf
commit dd6094d51e
5 changed files with 389 additions and 184 deletions
Download Patch File Download Diff File Expand all files Collapse all files

38
projects/FSO_transmission/first_analysis.m
View File

@@ -81,9 +81,23 @@ Pform = Pulseformer("fsym",fsym,"fdac",2*fsym,"pulse","rrc","pulselength",16,"al
Rx_matched = Pform.process(Scope_sig);
Rx_matched.spectrum("displayname",'Signal after matched filter','fignum',1);
% timing sync -> at this point we still have no symbol timing recovery, we
% try to do this with 2sps EQ!
[~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_matched.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1);
%%
sys = comm.SymbolSynchronizer('TimingErrorDetector', 'Gardner (non-data-aided)', ...
'SamplesPerSymbol', 2, ...
'DampingFactor', 0.7, ...
'NormalizedLoopBandwidth', 0.01);
Rx_symbolsync = Rx_matched;
[Rx_symbolsync.signal, timing_error] = sys(Rx_matched.signal);
plot(timing_error); % If this is a ramp, you have drift!
%% timing sync -> at this point we still have no symbol timing recovery, we
% % try to do this with 2sps EQ!
[~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_symbolsync.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1);
%% not working..
Rx_synced = Rx_synced_cell{1};
@@ -136,6 +150,22 @@ mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0,"eth_styl
[vnle_results, mlse_results] = vnle_postfilter_mlse(eq_v, pf_, mlse_, M, Rx_synced, Symbols, Bits, ...
"precode_mode", duob_mode, 'showAnalysis', 1, "postFFE", [], "eth_style_symbol_mapping", mapping_style);
mlse_results.metrics.print
mlse_results.metrics.print("description",'MLSE')
fprintf('My EQ: %.1e \n',mlse_results.metrics.BER);
fprintf('Paper: %.1e \n \n',ber_in_paper);
%% -------------------- DB target --------------------
mlse_db_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",M,'trellis_states',PAMmapper(M,0).levels);
eq_ = EQ("Ne",[50, 5, 5],"Nb",[0,0,0],"training_length",len_tr,"training_loops",5,"dd_loops",5, ...
"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1);
dbt_results = duobinary_target(eq_,mlse_db_, M, Rx_synced, Symbols, Bits, ...
"precode_mode", duob_mode, 'showAnalysis', 0, "postFFE", [],"eth_style_symbol_mapping",mapping_style);
dbt_results.metrics.print("description",'Duobinary');
mlse_results.metrics.print
fprintf('My EQ: %.1e \n',dbt_results.metrics.BER);
fprintf('Paper: %.1e \n \n',ber_in_paper);

61
projects/WDM/WDM_auswertung.m
View File

@@ -1,7 +1,58 @@
filename = "C:\Users\Silas\Documents\MATLAB\Datensätze\FWM_2025\WDM_20260106_131152_n261_14408161_10km_16ch_400ghz_co_alpha0_5.mat";
filename = "C:\Users\Silas\Documents\MATLAB\Datensätze\FWM_2025\WDM_20260106_131152_n262_14408165_10km_16ch_400ghz_co_alpha-0_5.mat";
filename = "C:\Users\Silas\Documents\MATLAB\Datensätze\FWM_2025\WDM_20260106_134921_n223_14409338_10km_16ch_400ghz_co_alpha0_0.mat";
res = load(filename,'res');
base = "C:\Users\Silas\Nextcloud\Cluster";
all_files = dir(fullfile(base, "**/*.mat"));
schemes = ["co","pair","alt","seg"];
% Preallocate as table (minimal + convenient)
T = table('Size',[0 10], ...
'VariableTypes', ["string","string","string","datetime","double","double","double","double","double","double"], ...
'VariableNames', ["folder","file","scheme","date","node","jobid","L_km","Nch","df_GHz","alpha"]);
% filename parser
rx = "^WDM_(?<date>\d{8})_(?<time>\d{6})_n(?<node>\d+)_(?<jobid>\d+)_" + ...
"(?<L>\d+)km_(?<Nch>\d+)ch_(?<df>\d+)ghz_(?<scheme>[a-z]+)_alpha(?<alpha>\d+(?:_\d+)?)\.mat$";
for k = 1:numel(all_files)
f = all_files(k);
folder = string(f.folder);
file = string(f.name);
tok = regexp(file, rx, 'names');
if isempty(tok)
continue
end
% scheme from filename is the most reliable
scheme = string(tok.scheme);
% optional: ignore unexpected schemes
if ~any(strcmpi(scheme, schemes)), continue; end
node = str2double(tok.node);
jobid = str2double(tok.jobid);
L_km = str2double(tok.L);
Nch = str2double(tok.Nch);
df_GHz = str2double(tok.df);
date = datetime(strcat(tok.date, tok.time), 'InputFormat','yyyyMMddHHmmss');
% alpha uses "_" as decimal separator in your filenames
alpha = str2double(strrep(tok.alpha, "_", "."));
T(end+1,:) = {folder, file, scheme,date, node, jobid, L_km, Nch, df_GHz, alpha};
end
%%
idx = strcmp(T.scheme,"co") & ...
T.alpha == 0.4 & ...
T.date >= datetime(2026,1,9) & ...
T.date <= datetime(2026,1,10,23,59,59);
T_sel = T(idx,:);
i = 2;
res = load(fullfile(T_sel.folder(i),T_sel.file(i)),'res');
res = res.res;
%%
@@ -9,4 +60,4 @@ res = res.res;
S = plot_BER_vs_ROP(res, 'fec', 3.8e-3);
%% Routine B: violin plot (independent)
plot_FEC_violin(res, 'tech','VNLE', 'fec',3.8e-3, 'eval_ptr',3, 'ylim',[-10 -3]);
plot_FEC_violin(res, 'tech','VNLE', 'fec',3.8e-3, 'ylim',[-10 0]);

205
projects/WDM/WDM_model_10km_queue.m
View File

@@ -6,7 +6,7 @@ function WDM_model_10km_queue(options)
%%% Run parameters
arguments
options.num_channels = 16;
options.num_channels = 8;
options.channel_spacing = 400e9;
options.fiber_length_km = 10;
options.rand_key = 1;
@@ -26,7 +26,7 @@ warning('off','MATLAB:datetime:AmbiguousTimeZone');
%% How many workers?
cpus = str2double(getenv('SLURM_CPUS_PER_TASK'));
if ~isfinite(cpus) || cpus < 1, cpus = max(1, feature('numcores')); end
if ~isfinite(cpus) || cpus < 1, cpus = max(1, feature('numcores'))-1; end
% Use a per-job, node-local JobStorageLocation to avoid stale locks on $HOME
tmpbase = getenv('TMPDIR');
@@ -43,7 +43,7 @@ c.JobStorageLocation = jsl;
p = gcp('nocreate');
if isempty(p) || p.NumWorkers ~= cpus
if ~isempty(p), delete(p); end
p = parpool(c, cpus);
p = parpool(c, cpus,"IdleTimeout",300);
end
fprintf('parpool up with %d workers; JobStorage=%s\n', p.NumWorkers, c.JobStorageLocation);
@@ -119,7 +119,7 @@ upsample_pow = 2^nextpow2(upsample_required);
s.f_opt = fdac*kover*upsample_pow;
s.f_opt_nyq = s.f_opt/2;
s.rop = -10:0.5:0;
s.rop = -10:1:0;
%% ---------- Intermediate evaluation points (distance dimension) ----------
segment_length = 1; % km (must match fiber loop below)
@@ -184,7 +184,7 @@ for realiz = 1:s.num_realiz
for l = 1:N
[Digi_sig,Symbols{l},Tx_bits{l}] = PAMsource( ...
"fsym",fsym,"M",s.M,"order",17,"useprbs",0, ...
"fsym",fsym,"M",s.M,"order",15,"useprbs",0, ...
"fs_out",fdac, ...
"applyclipping",0,"clipfactor",1.5, ...
"applypulseform",apply_pulsef,"pulseformer",Pform, ...
@@ -282,20 +282,21 @@ for realiz = 1:s.num_realiz
[F, meta, output_ffe, output_dfe, output_vnle, output_mlse, output_dbt] = ...
throttle_inflight(F, meta, output_ffe, output_dfe, output_vnle, output_mlse, output_dbt, maxInFlight);
% Enqueue one job per (l,ri)
F(end+1,1) = parfeval(p, @rx_job, 5, ...
m = struct('l',l,'ri',ri,'realiz',realiz,'eval_ptr',eval_ptr);
[F, meta] = enqueue_atomic(F, meta, p, @rx_job, 5, ...
Opt_sig_wdm_demux{l}, ...
s.rop(ri), ...
Symbols{l}, Tx_bits{l}, ...
s, l, ri, realiz, eval_ptr, ...
fdac, kover, fadc, fsym, ...
len_tr, mu_dc, mu_ffe, mu_dfe, dfe_order, duob_mode, ...
memlog_dir);
memlog_dir, ...
'META', m);
fprintf('Enqueued job %d/%d for realiz %d/%d, l=%d/%d, ri=%d/%d at 0 km (inflight=%d/%d)\n', ...
cnt, total, realiz, s.num_realiz, l, N, ri, length(s.rop), numel(F), maxInFlight);
meta(end+1) = struct('l',l,'ri',ri,'realiz',realiz,'eval_ptr',eval_ptr);
end
end
@@ -304,11 +305,6 @@ for realiz = 1:s.num_realiz
end
% -----------------------------------------------------
save_results_per_realization( ...
s, eval_dist_km, ...
output_ffe, output_dfe, output_vnle, output_mlse, output_dbt, ...
output_root, fname);
for seg = 1:nSegments
fprintf('Realiz %d/%d: Segment %d/%d \n', realiz, s.num_realiz, seg, nSegments);
@@ -336,19 +332,21 @@ for realiz = 1:s.num_realiz
[F, meta, output_ffe, output_dfe, output_vnle, output_mlse, output_dbt] = ...
throttle_inflight(F, meta, output_ffe, output_dfe, output_vnle, output_mlse, output_dbt, maxInFlight);
F(end+1,1) = parfeval(p, @rx_job, 5, ...
m = struct('l',l,'ri',ri,'realiz',realiz,'eval_ptr',eval_ptr);
[F, meta] = enqueue_atomic(F, meta, p, @rx_job, 5, ...
Opt_sig_wdm_demux{l}, ...
s.rop(ri), ...
Symbols{l}, Tx_bits{l}, ...
s, l, ri, realiz, eval_ptr, ...
fdac, kover, fadc, fsym, ...
len_tr, mu_dc, mu_ffe, mu_dfe, dfe_order, duob_mode, ...
memlog_dir);
memlog_dir, ...
'META', m);
fprintf('Seg: %d - Enqueued job %d/%d for realiz %d/%d, l=%d/%d, ri=%d/%d (inflight=%d/%d)\n', ...
seg, cnt, total, realiz, s.num_realiz, l, N, ri, length(s.rop), numel(F), maxInFlight);
meta(end+1) = struct('l',l,'ri',ri,'realiz',realiz,'eval_ptr',eval_ptr);
end
end
@@ -546,15 +544,21 @@ end
first = true;
while ~isempty(F)
assert(numel(F)==numel(meta))
idxUn = strcmp({F.State}, 'unavailable');
if any(idxUn)
for ii = find(idxUn)
m = meta(ii);
fprintf('[UNAVAILABLE] realiz=%d eval=%d l=%d ri=%d\n', ...
m.realiz, m.eval_ptr, m.l, m.ri);
end
warning('collect_done: dropping %d unavailable futures.', sum(idxUn));
F(idxUn) = [];
meta(idxUn) = [];
if isempty(F), break; end
end
if first
timeout = timeout_first;
first = false;
@@ -703,65 +707,146 @@ catch
end
end
function [F, meta, output_ffe, output_dfe, output_vnle, output_mlse, output_dbt] = ...
throttle_inflight(F, meta, output_ffe, output_dfe, output_vnle, output_mlse, output_dbt, maxInFlight)
%THROTTLE_INFLIGHT
% Enforces a HARD cap on outstanding futures (running + queued + finished-not-yet-fetched).
% Blocks until numel(F) < maxInFlight, then returns.
% Hard cap on outstanding futures. Robust against:
% - 'unavailable' futures (ID=-1 or State unavailable)
% - finished-with-error futures (State may be 'finished' with non-empty Error)
% - rare F/meta desync
% - race-y property access on FevalFuture arrays (use snapshot Fs)
while numel(F) >= maxInFlight
% Wait until at least one finishes OR fails
% (wait('finished',1) only wakes for finished, so we also handle failures via periodic checks)
% Wait for at least one to finish (or time out quickly)
try
wait(F, 'finished', 1);
catch
% If wait fails for some reason, fall back to a short pause
pause(0.1);
end
% Harvest finished and failed futures (robust snapshot)
n = numel(F);
% ---- Ensure F/meta alignment (best-effort repair) ----
if numel(F) ~= numel(meta)
warning('throttle_inflight: F/meta desync (F=%d meta=%d). Attempting repair.', numel(F), numel(meta));
% Safety: meta must match F length
if numel(meta) ~= n
% Keep them aligned (rare; indicates earlier mismatch)
m = min(numel(meta), n);
F = F(1:m);
% Prefer dropping unusable futures first (likely the unmatched ones)
Fs = F; % snapshot
n = numel(Fs);
st = cell(1,n);
idbad = false(1,n);
for k = 1:n
try, st{k} = Fs(k).State; catch, st{k} = ''; end
try, idbad(k) = (Fs(k).ID < 0); catch, idbad(k) = false; end
end
idxBad = find(strcmp(st,'unavailable') | idbad);
if numel(F) > numel(meta) && ~isempty(idxBad)
need = numel(F) - numel(meta);
idxBad = idxBad(1:min(numel(idxBad), need));
idxBad = idxBad(idxBad >= 1 & idxBad <= numel(F));
if ~isempty(idxBad)
warning('throttle_inflight: dropping %d bad futures to restore alignment.', numel(idxBad));
F(idxBad) = [];
end
end
% Final fallback: truncate both to min length
m = min(numel(F), numel(meta));
F = F(1:m);
meta = meta(1:m);
n = m;
end
states = cell(1,n);
if isempty(F)
break;
end
% ===================== IMPORTANT: snapshot =====================
Fs = F; % stable snapshot for property access
n = numel(Fs);
% ===============================================================
% ---- Compute unavailable indices from snapshot ONLY ----
st = cell(1,n);
idbad = false(1,n);
for k = 1:n
states{k} = F(k).State;
try, st{k} = Fs(k).State; catch, st{k} = ''; end
try, idbad(k) = (Fs(k).ID < 0); catch, idbad(k) = false; end
end
doneFinished = strcmp(states, 'finished');
doneFailed = strcmp(states, 'failed');
done = doneFinished | doneFailed;
idxUn = find(strcmp(st,'unavailable') | idbad);
% Sanitize indices (prevents out-of-range deletions even if something weird happens)
idxUn = idxUn(idxUn >= 1 & idxUn <= numel(F));
if ~isempty(idxUn)
for ii = idxUn(:).'
if ii <= numel(meta)
mm = meta(ii);
fprintf('[UNAVAILABLE@throttle] realiz=%d eval=%d l=%d ri=%d\n', ...
mm.realiz, mm.eval_ptr, mm.l, mm.ri);
else
fprintf('[UNAVAILABLE@throttle] meta_missing for ii=%d\n', ii);
end
end
warning('throttle_inflight: dropping %d unavailable/bad futures.', numel(idxUn));
% Delete in descending order is safest for structs (not strictly necessary, but robust)
idxUn = sort(idxUn, 'descend');
F(idxUn) = [];
meta(idxUn) = [];
continue; % re-check cap after shrinking
end
% ---- Determine "done" (finished or finished-with-error) ----
isFinished = strcmp(st, 'finished');
hasErr = false(1,n);
for k = 1:n
try
hasErr(k) = ~isempty(Fs(k).Error); % use snapshot object
catch
hasErr(k) = false;
end
end
done = isFinished | hasErr;
idxDone = find(done);
idxDone = idxDone(idxDone >= 1 & idxDone <= numel(F)); % sanitize
if isempty(idxDone)
pause(0.05);
continue;
end
% Fetch outputs for finished; mark failed as empty results
% ---- Harvest done futures ----
for jj = 1:numel(idxDone)
ii = idxDone(jj);
m = meta(ii);
if doneFinished(ii)
[ffe_r, dfe_r, vnle_r, mlse_r, dbt_r] = fetchOutputs(F(ii));
if isFinished(ii) && ~hasErr(ii)
try
[ffe_r, dfe_r, vnle_r, mlse_r, dbt_r] = fetchOutputs(F(ii));
catch ME
warning('throttle_inflight: fetchOutputs failed: realiz=%d eval=%d l=%d ri=%d | %s', ...
m.realiz, m.eval_ptr, m.l, m.ri, ME.message);
ffe_r = []; dfe_r = []; vnle_r = []; mlse_r = []; dbt_r = [];
end
else
% Failed
% finished-with-error
try
err = F(ii).Error;
warning('rx_job failed: realiz=%d eval=%d l=%d ri=%d | %s', ...
m.realiz, m.eval_ptr, m.l, m.ri, err.message);
if ~isempty(err)
warning('rx_job error: realiz=%d eval=%d l=%d ri=%d | %s', ...
m.realiz, m.eval_ptr, m.l, m.ri, err.message);
else
warning('rx_job error: realiz=%d eval=%d l=%d ri=%d | (unknown error)', ...
m.realiz, m.eval_ptr, m.l, m.ri);
end
catch
warning('rx_job error: realiz=%d eval=%d l=%d ri=%d | (error unreadable)', ...
m.realiz, m.eval_ptr, m.l, m.ri);
end
ffe_r = []; dfe_r = []; vnle_r = []; mlse_r = []; dbt_r = [];
end
@@ -773,9 +858,43 @@ while numel(F) >= maxInFlight
output_dbt{m.l, m.ri, m.realiz, m.eval_ptr} = dbt_r;
end
% Remove harvested futures/meta entries using numeric indices (safe)
% Remove harvested entries (descending for safety)
idxDone = sort(idxDone, 'descend');
F(idxDone) = [];
meta(idxDone) = [];
end
end
function [F, meta] = enqueue_atomic(F, meta, p, fun, nOut, varargin)
%ENQUEUE_ATOMIC Enqueue a parfeval job and append matching meta atomically.
% Usage:
% m = struct('l',l,'ri',ri,'realiz',realiz,'eval_ptr',eval_ptr);
% [F, meta, f] = enqueue_atomic(F, meta, p, @rx_job, 5, args..., 'META', m);
% Split varargin into args + meta struct
idx = find(strcmp(varargin, 'META'), 1, 'last');
if isempty(idx) || idx == numel(varargin)
error('enqueue_atomic: META struct must be provided as last named argument.');
end
args = varargin(1:idx-1);
m = varargin{idx+1};
% Create future first (local variable), but don't mutate F/meta yet
f = parfeval(p, fun, nOut, args{:});
% If future is unusable, do not append (prevents meta shift)
if f.ID < 0 || strcmp(f.State, 'unavailable')
warning('enqueue_atomic: got unusable future (ID=%d, State=%s). Dropping enqueue: realiz=%d eval=%d l=%d ri=%d', ...
f.ID, string(f.State), m.realiz, m.eval_ptr, m.l, m.ri);
return;
end
% Now append both together (atomic w.r.t. fprintf etc.)
F(end+1,1) = f;
meta(end+1,1) = m;
end

1
projects/WDM/plot_BER_vs_ROP.m
View File

@@ -85,7 +85,6 @@ for eval_ptr = eval_list
ffe_mat = extractCompleteBER(ffe_cells); % N_rop x K
if ~isempty(ffe_mat)
% Your current style: plot each realization curve (no boundedline)
plot(rop, ffe_mat, 'Color', cols(ch,:), ...
'DisplayName', sprintf('FFE @ %dnm', round(wavelengthplan(ch))));
end

266
projects/WDM/plot_FEC_violin.m
View File

@@ -1,29 +1,38 @@
function plot_FEC_violin(res, varargin)
% plot_fec_violin_from_res Self-contained violin plot (no dependency on other routines)
% plot_FEC_violin Violin plot of ROP-at-FEC crossings per wavelength/channel
% Uses the Bechtold community violinplot (MATLAB Central 45134).
%
% Requirements:
% - community violinplot must be on path, either as:
% (a) violinplot.m (shadows MATLAB official), OR
% (b) violinplot_community.m (renamed + function name adjusted)
%
% Usage:
% plot_fec_violin_from_res('WDM_20260106_....mat', 'tech', 'FFE', 'fec', 3.8e-3, 'eval_ptr', 2, 'ylim', [-10 -6]);
% plot_fec_violin_from_res(res, 'tech', 'FFE', 'fec', 3.8e-3, 'eval_ptr', 2);
%
% Notes:
% - Works with partially-filled res.* cell arrays (missing eval distances / missing realizations).
% - Uses only complete realizations per channel (all ROP points present).
% - Requires violinplot.m on path.
% plot_FEC_violin(res, 'tech','VNLE', 'fec',3.8e-3, 'eval_ptr',[], 'ylim',[-10 -6], 'bandwidth',0.15);
% -------- args --------
p = inputParser;
p.addParameter('tech', 'FFE', @(x)ischar(x) || isstring(x));
p.addParameter('fec', 3.8e-3, @(x)isnumeric(x)&&isscalar(x));
p.addParameter('eval_ptr', [], @(x)isnumeric(x)&&isscalar(x));
p.addParameter('ylim', [], @(x)isnumeric(x)&&numel(x)==2);
p.addParameter('title_prefix', 'ROP at FEC crossing', @(x)ischar(x)||isstring(x));
p.addParameter('tech', 'VNLE', @(x)ischar(x)||isstring(x));
p.addParameter('fec', 3.8e-3, @(x)isnumeric(x)&&isscalar(x));
p.addParameter('eval_ptr', [], @(x)isnumeric(x)&&isscalar(x));
p.addParameter('ylim', [], @(x)isnumeric(x)&&numel(x)==2);
p.addParameter('bandwidth', [], @(x)isnumeric(x)&&isscalar(x));
p.parse(varargin{:});
opt = p.Results;
if exist('violinplot','file') ~= 2
error('violinplot.m not found on path. Add it to path first.');
tech = upper(string(opt.tech));
fec = opt.fec;
% -------- pick result field --------
switch tech
case "FFE", C4 = res.ffe;
case "DFE", C4 = res.dfe;
case "VNLE", C4 = res.vnle;
case "MLSE", C4 = res.mlse;
case "DBT", C4 = res.dbt;
otherwise, error('Unknown tech "%s". Use FFE/DFE/VNLE/MLSE/DBT.', tech);
end
% ---------- Metadata ----------
rop = res.settings.rop(:);
wavelengthplan = res.settings.wavelengthplan(:);
distances = res.eval_dist_km(:);
@@ -31,103 +40,97 @@ distances = res.eval_dist_km(:);
N_ch = numel(wavelengthplan);
N_rop = numel(rop);
% Determine distances dim robustly from the stored cell arrays (prefer ffe)
dims = size(res.ffe);
dims = size(C4);
if numel(dims) < 4, dims(end+1:4) = 1; end
N_distances = dims(4);
if numel(distances) ~= N_distances
distances = (1:N_distances).';
N_dist = dims(4);
eval_ptr = opt.eval_ptr;
if isempty(eval_ptr), eval_ptr = N_dist; end
eval_ptr = max(1, min(N_dist, eval_ptr));
% -------- compute crossings per channel + "completeness" --------
S_cell = cell(1, N_ch); % crossings values (finite only)
K_total = zeros(1, N_ch); % number of complete realizations for that channel
K_cross = zeros(1, N_ch); % number of realizations that crossed
for ch = 1:N_ch
cells = sliceCells4D(C4, ch, eval_ptr, N_rop); % N_rop x N_realiz (may contain [])
[S_cell{ch}, K_total(ch), K_cross(ch)] = crossings_and_counts(rop, cells, fec);
end
% Choose eval distance
eval_ptr = p.Results.eval_ptr;
if isempty(eval_ptr)
eval_ptr = N_distances; % default: last available
end
if eval_ptr < 1 || eval_ptr > N_distances
error('eval_ptr=%d out of range. Available: 1..%d', eval_ptr, N_distances);
% -------- build vector+category representation (NO NaNs in xAll) --------
catLabels = arrayfun(@(nm)sprintf('%d nm', round(nm)), wavelengthplan, 'UniformOutput', false);
catsAll = categorical(strings(0,1), catLabels, 'Ordinal', false);
xAll = zeros(0,1);
for ch = 1:N_ch
v = S_cell{ch}(:);
if ~isempty(v)
xAll = [xAll; v];
catsAll = [catsAll; repmat( ...
categorical(string(catLabels{ch}), catLabels, 'Ordinal', false), ...
numel(v), 1)];
end
end
% Choose technique field
tech = upper(string(p.Results.tech));
switch tech
case "FFE", C4 = res.ffe;
case "DFE", C4 = res.dfe;
case "VNLE", C4 = res.vnle;
case "MLSE", C4 = res.mlse;
case "DBT", C4 = res.dbt;
otherwise
error('Unknown tech "%s". Use one of: FFE, DFE, VNLE, MLSE, DBT.', tech);
end
fec = p.Results.fec;
% ---------- Compute crossings per channel ----------
% ---------- Plot ----------
figure('Name', sprintf('%s violin @ %s', tech, distLabel(distances, eval_ptr)));
% -------- plot --------
figure('Name', sprintf('%s FEC violin @ %s', tech, distLabel(distances, eval_ptr)));
hold on;
c = linspecer(N_distances);
for eval_ptr = 1:N_distances
S_cell = cell(1, N_ch);
for ch = 1:N_ch
cells = sliceCells4D(C4, ch, eval_ptr, N_rop); % N_rop x N_realiz cell
S_cell{ch} = fecCrossings_only(rop, cells, fec); % 1 x K crossings (may be empty)
end
% Pad to matrix for violinplot: rows=realizations (max K), cols=channels
Kmax = max(cellfun(@numel, S_cell));
if isempty(Kmax) || Kmax == 0
warning('No FEC crossings found for %s at eval_ptr=%d (%s).', tech, eval_ptr, distLabel(distances, eval_ptr));
return;
end
S_mat = NaN(Kmax, N_ch);
for ch = 1:N_ch
k = numel(S_cell{ch});
if k > 0
S_mat(1:k, ch) = S_cell{ch}(:);
end
end
catLabels = arrayfun(@(nm) sprintf('%d nm', nm), round(wavelengthplan), 'UniformOutput', false);
if numel(S_mat)>numel(wavelengthplan)
S_mat(isnan(S_mat)) = 0;
v=violinplot(S_mat, catLabels, ...
'ViolinColor', c(eval_ptr,:), ...
'ViolinAlpha', 0.10, ...
'MarkerSize', 20, ...
'ShowMedian', true, ...
'EdgeColor', c(eval_ptr,:), ...
'ShowWhiskers', false, ...
'ShowData', true, ...
'ShowBox', false, ...
'Bandwidth', 0.05);
else
channs = 1:numel(wavelengthplan);
scatter( wavelengthplan ,S_mat,10,c(eval_ptr,:),'Marker','o','LineWidth',2,'DisplayName',distLabel(distances, eval_ptr));
channs = 1:numel(wavelengthplan);
scatter( round(wavelengthplan(S_mat==0)),S_mat(S_mat==0),10,c(eval_ptr,:),'Marker','o','HandleVisibility','off');
end
% Violin only if we have any data at all
if ~isempty(xAll)
violinplot_community(xAll, catsAll,'Bandwidth', 0.15, ... % KDE bandwidth (critical!)
'ViolinColor', [0.2 0.4 0.8], ... % or Nx3 for per-channel colors
'ViolinAlpha', 0.15, ...
'EdgeColor', [0.3 0.3 0.3], ...
'MarkerSize', 14, ...
'ShowData', true, ...
'ShowMean', false, ...
'ShowMedian', true, ...
'ShowBox', false, ...
'ShowWhiskers', false);
else
warning('No FEC crossings found at eval_ptr=%d (%s). Plotting only markers.', eval_ptr, distLabel(distances, eval_ptr));
set(gca,'XTick',1:N_ch,'XTickLabel',catLabels);
end
% Marker X positions (violin groups are at 1..N_ch)
xpos = 1:N_ch;
% -------- overlay mean markers (black/red/blue) --------
yBlue = -9;
for ch = 1:N_ch
if K_total(ch) == 0 || K_cross(ch) == 0
% no complete realizations OR none crossed -> blue X at -9
scatter(xpos(ch), yBlue, 80, 'x', 'LineWidth', 2, 'MarkerEdgeColor', [0 0 1], 'HandleVisibility','off');
continue;
end
mu = mean(S_cell{ch}, 'omitnan');
if K_cross(ch) < K_total(ch)
% some complete realizations did not cross -> red X
scatter(xpos(ch), mu, 80, 'x', 'LineWidth', 2, 'MarkerEdgeColor', [1 0 0], 'HandleVisibility','off');
else
% all complete realizations crossed -> black X
scatter(xpos(ch), mu, 80, 'x', 'LineWidth', 2, 'MarkerEdgeColor', [0 0 0], 'HandleVisibility','off');
end
end
title(sprintf('%s | BER %.2e | %s', tech, fec, distLabel(distances, eval_ptr)));
ylabel('ROP at FEC crossing');
title(sprintf('%s: %s | BER %.2e | %s', ...
p.Results.title_prefix, tech, fec, distLabel(distances, eval_ptr)));
grid on; box on;
if ~isempty(p.Results.ylim)
ylim(p.Results.ylim);
if ~isempty(opt.ylim)
ylim(opt.ylim);
end
end
%% ================= helpers (self-contained) =================
% ================= helpers =================
function cells2D = sliceCells4D(C4, ch, eval_ptr, N_rop)
dims = size(C4);
@@ -136,68 +139,64 @@ N_realiz = dims(3);
cells2D = cell(N_rop, N_realiz);
chMax = dims(1);
ropMax = dims(2);
rMax = dims(3);
eMax = dims(4);
if ch > chMax || eval_ptr > eMax
if ch > dims(1) || eval_ptr > dims(4)
return;
end
ropUse = min(N_rop, ropMax);
rUse = min(N_realiz, rMax);
ropUse = min(N_rop, dims(2));
rUse = min(N_realiz, dims(3));
tmp = squeeze(C4(ch, 1:ropUse, 1:rUse, eval_ptr));
tmp = reshape(tmp, ropUse, rUse);
cells2D(1:ropUse, 1:rUse) = tmp;
end
function lbl = distLabel(distances, eval_ptr)
if eval_ptr <= numel(distances)
d = distances(eval_ptr);
if isfinite(d)
lbl = sprintf('%.0f km', d);
return;
end
end
lbl = sprintf('eval\\_%d', eval_ptr);
function [S, K_total, K_cross] = crossings_and_counts(rop, cellsNxR, fec)
% counts "complete" realizations and how many of those crossed
% S returns only finite crossings (no NaN).
Y = extractCompleteBER(cellsNxR); % N_rop x K_total
K_total = size(Y,2);
if K_total == 0
S = [];
K_cross = 0;
return;
end
function S = fecCrossings_only(rop, cellsNxR, fec)
% Returns 1xK crossing positions for all complete realizations (drops bad ones).
Y = extractCompleteBER(cellsNxR); % N_rop x K
if isempty(Y), S = []; return; end
% Optional quality gate (same as your previous scripts)
% optional quality gate (keep your style)
ok = mean(Y,1,'omitnan') <= 0.1;
Y = Y(:, ok);
if isempty(Y), S = []; return; end
Y = Y(:,ok);
K_total = size(Y,2);
nR = size(Y,2);
S = nan(1,nR);
if K_total == 0
S = [];
K_cross = 0;
return;
end
rop = rop(:);
for j = 1:nR
S = nan(1, K_total);
for j = 1:K_total
y = Y(:,j);
above = (y > fec);
idx = find(above(1:end-1) & ~above(2:end), 1, 'first');
if ~isempty(idx)
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);
if isfinite(y1) && isfinite(y2) && (y2 ~= y1)
t = (fec - y1) / (y2 - y1);
S(j) = x1 + t*(x2 - x1);
end
end
end
S = S(isfinite(S)); % keep only successful crossings
K_cross = sum(isfinite(S));
S = S(isfinite(S));
end
function Y = extractCompleteBER(cellSlice)
% Keep only those realization columns where ALL N_rop entries are non-empty
if isempty(cellSlice), Y = []; return; end
nR = size(cellSlice,2);
keep = false(1,nR);
@@ -208,6 +207,13 @@ for r = 1:nR
end
if ~any(keep), Y = []; return; end
Y = cellfun(@(c) c.metrics.BER, cellSlice(:,keep), 'UniformOutput', true);
end
function lbl = distLabel(distances, eval_ptr)
if eval_ptr <= numel(distances) && isfinite(distances(eval_ptr))
lbl = sprintf('%.0f km', distances(eval_ptr));
else
lbl = sprintf('eval_%d', eval_ptr);
end
end