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

Merge branch 'main' of cau-git.rz.uni-kiel.de:nt/mitarbeiter/silas/imdd_simulation

Browse Source 
# Conflicts:
#	projects/WDM/WDM_auswertung.m
#	projects/WDM/WDM_model_old.m
This commit is contained in:
Silas Oettinghaus
2026-01-12 12:49:16 +01:00
parent 16846d4eb7 3b44315eea
commit ad14cd738f
26 changed files with 2941 additions and 1070 deletions
Download Patch File Download Diff File Expand all files Collapse all files

293
projects/WDM/WDM_auswertung.m
View File

@@ -1,248 +1,63 @@
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');
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;
% 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
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;
idx = strcmp(T.scheme,"co") & ...
T.alpha == 0.4 & ...
T.date >= datetime(2026,1,9) & ...
T.date <= datetime(2026,1,10,23,59,59);
end
T_sel = T(idx,:);
i = 2;
res = load(fullfile(T_sel.folder(i),T_sel.file(i)),'res');
res = res.res;
%%
% Routine A: plot BER curves and compute crossings
S = plot_BER_vs_ROP(res, 'fec', 3.8e-3);
%% ================= 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
%% Routine B: violin plot (independent)
plot_FEC_violin(res, 'tech','VNLE', 'fec',3.8e-3, 'ylim',[-10 0]);

265
projects/WDM/WDM_model.m
View File

@@ -1,25 +1,96 @@
%%% Run parameters
% TX
% --- FIRST LINE: evaluate settings located beside this script ---
run(fullfile(fileparts(mfilename('fullpath')),'WDM_settings.m'));
num_realiz = 50;
s.wavelengthplan = calcWavelengthPlan(16,400e9,1310);
function WDM_model(options)
arguments
options.num_channels = 16;
options.channel_spacing = 400e9;
options.fiber_length_km = 0;
options.rand_key = 1;
options.num_realiz = 1;
options.fwm_mitigation_technique = "co";
end
%%
% Add the imdd_simulation framework to the path
if ispc
addpath(genpath('C:\Users\Silas\Documents\MATLAB\imdd_simulation'));
else
% Linux path on the cluster
addpath(genpath('/work_beegfs/sutef391/imdd_simulation'));
end
% Quiet the ambiguous CET warning (best is to set TZ in sbatch; see below)
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
% Use a per-job, node-local JobStorageLocation to avoid stale locks on $HOME
% Prefer $TMPDIR if your cluster provides it, else tempdir().
tmpbase = getenv('TMPDIR');
if isempty(tmpbase), tmpbase = tempdir; end
jsl = fullfile(tmpbase, sprintf('matlab_jobstorage_%s_%s', ...
getenv('USER'), getenv('SLURM_JOB_ID')));
if ~exist(jsl,'dir'); mkdir(jsl); end
% Configure the local cluster explicitly and start the pool
c = parcluster('local');
c.NumWorkers = cpus;
c.JobStorageLocation = jsl;
p = gcp('nocreate');
if isempty(p) || p.NumWorkers ~= cpus
if ~isempty(p), delete(p); end
p = parpool(c, cpus); % avoids the “queued” state
end
fprintf('parpool up with %d workers; JobStorage=%s\n', p.NumWorkers, c.JobStorageLocation);
% result filename (timestamp + optional job id)
t = datetime('now','TimeZone','local','Format','yyyyMMdd_HHmmss');
jobid = getenv('SLURM_JOB_ID'); if isempty(jobid), jobid = 'nojid'; end
host = getenv('HOSTNAME'); if isempty(host), host = 'localhost'; end
% Output directory depends on platform
% create/ use folders foroptions.fiber_length_km, options.num_channels, options.channel_spacing.*1e-9, options.fwm_mitigation_technique
foldname = sprintf('%dkm_%dch_%dghz_%s', options.fiber_length_km(end), options.num_channels, options.channel_spacing.*1e-9, options.fwm_mitigation_technique);
if ispc
output_root = fullfile('C:\Users\Silas\Documents\MATLAB\Datensätze\FWM_2025\',foldname,'\');
else
output_root = fullfile('/work_beegfs/sutef391/results_WDM',foldname,'\');
end
if ~exist(output_root,'dir'), mkdir(output_root); end
% Build filename
t = datetime('now','TimeZone','local','Format','yyyyMMdd_HHmmss');
jobid = getenv('SLURM_JOB_ID'); if isempty(jobid), jobid = 'nojid'; end
host = getenv('HOSTNAME'); if isempty(host), host = 'localhost'; end
fname = sprintf('WDM_%s_%s_%s_%dkm_%dch_%dghz_%s.mat', char(t), host, jobid, options.fiber_length_km(end), options.num_channels, options.channel_spacing.*1e-9, options.fwm_mitigation_technique);
%%
s.num_realiz = options.num_realiz;
% s.wavelengthplan = calcWavelengthPlan(16,400e9,1310);
s.wavelengthplan = calcWavelengthPlan(options.num_channels,options.channel_spacing,1310);
% wavelengthplan = [1295,1305,1315,1325];
link_length = 10;
pmd = 0.1;
gamma = 0.0023;
link_length = options.fiber_length_km;
s.pmd = 0.1;
s.gamma = 0.0023;
M = 4;
m = floor(log2(M)*10)/10;
s.M = 4;
m = floor(log2(s.M)*10)/10;
fsym = 112e9;
fdac = 2*fsym;
fadc = 2*fsym;
s.random_key = 100;
fadc = 120000000000;
s.random_key = options.rand_key;
% Laser / s.Modulator
vbias_rel = 0.5;
u_pi = 3.2;
u_pi = 4.6;
vbias = -vbias_rel*u_pi;
laser_linewidth = 0e6;
@@ -38,7 +109,7 @@ mu_dc = 0.005;
mu_ffe = [mu_ffe1 mu_ffe3 mu_ffe3];
mu_dfe = 0.0004;
%DB Stuff
%DB Stuff
db_precode = 0;
db_encode = 0;
duob_mode = db_mode.no_db;
@@ -49,10 +120,10 @@ Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",1
N = numel(s.wavelengthplan);
f_plan = physconst('lightspeed')./(s.wavelengthplan.*1e-9);
margin = 25e12; % some THz left and right
margin = 5e12; % some THz left and right
f_span = (max(f_plan)+margin)-(min(f_plan)-margin);
f_nyq = f_span/2;
kover = 8;
kover = 4;
upsample_required = f_nyq./(fdac*kover/2);
upsample_pow = 2^nextpow2(upsample_required);
upsample_ceil = ceil(upsample_required);
@@ -64,18 +135,35 @@ signal_cell = {};
Symbols = {};
Tx_bits = {};
s.rop = -6:0.75:-0.75;
s.rop = -12:0.75:-0.75;
output_ffe = cell(length(s.wavelengthplan),length(s.rop),s.num_realiz);
output_vnle = cell(length(s.wavelengthplan),length(s.rop),s.num_realiz);
output_mlse = cell(length(s.wavelengthplan),length(s.rop),s.num_realiz);
output_dbt = cell(length(s.wavelengthplan),length(s.rop),s.num_realiz);
s.p_launch = 3;
s.p = options.fwm_mitigation_technique;
switch s.p
case "co"
pol_rot = 100.*ones(1,length(s.wavelengthplan));
d_local = 0;
case "pair"
pol_rot = repmat([100,100,0,0],1,length(s.wavelengthplan)/4);
d_local = 0;
case "alt"
pol_rot = repmat([100,0,100,0],1,length(s.wavelengthplan)/4);
d_local = 0;
case "seg"
pol_rot = 100.*ones(1,length(s.wavelengthplan));
d_local = 3;
end
for realiz = 1:s.num_realiz
parfor l = 1:N
[Digi_sig,Symbols{l},Tx_bits{l}] = PAMsource(...
"fsym",fsym,"M",s.M,"order",18,"useprbs",0,...
"fs_out",fdac,...
@@ -84,93 +172,97 @@ for realiz = 1:s.num_realiz
"randkey",s.random_key+l+realiz,...
"db_precode",db_precode,"db_encode",db_encode,...
"mrds_code",0,"mrds_blocklength",512,"duobinary_mode",duob_mode).process();
% Digi_sig.spectrum("fignum",101,"displayname",'bla','normalizeTo0dB',0,'lambda0_nm',1310,'useWavelengthAxis',0);
Lp_awg = Filter('filtdegree',3,"f_cutoff",100e9,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true);
El_sig = AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",1,"kover",kover,"bit_resolution",12,"upsampling_method","samplehold","precomp_sinc_rolloff",0,"H_lpf",Lp_awg,"dac_max",0.6,"dac_min",-0.6).process(Digi_sig);
Lp_awg = Filter('filtdegree',3,"f_cutoff",56e9,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true);
El_sig = AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",1,"kover",kover,"bit_resolution",6,"upsampling_method","samplehold","precomp_sinc_rolloff",0,"H_lpf",Lp_awg,"dac_max",0.6,"dac_min",-0.6).process(Digi_sig);
% El_sig = s.M8199B("kover",kover).process(Digi_sig);
% El_sig.spectrum("fignum",101,"displayname",'bla','normalizeTo0dB',0,'lambda0_nm',1310,'useWavelengthAxis',0);
%%%%% Electrical Driver Amplifier %%%%%%
El_sig = El_sig.normalize("mode","oneone");
scaling = 0.6*(u_pi/2-abs(vbias-u_pi/2)); % scale to 60% of available modulator curve
El_sig = El_sig .* scaling;
% El_sig = El_sig.setPower(1,"dBm");
% figure;histogram(El_sig.signal);
%%%%% s.MODULATE E/O CONVERSION %%%%%
Eml_out = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs,"lambda",s.wavelengthplan(l),"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",s.random_key+l+realiz).process(El_sig);
signal_cell{l} = Polarization_Controller("mode","rot_power","desired_power",100).process(Eml_out);
signal_cell{l} = Polarization_Controller("mode","rot_power","desired_power",pol_rot(l)).process(Eml_out);
end
Opt_sig_wdm = Optical_Multiplex("fs_in",fdac*kover,"fs_out",upsample_pow*fdac*kover,...
"lambda_center",1310,"random_key",0,"filtype",1,"B",200e9).process(signal_cell);
Opt_sig_wdm = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",3+10*log10(N)).process(Opt_sig_wdm);
"lambda_center",1310,"random_key",0,"filtype",1,"B",120e9).process(signal_cell);
Opt_sig_wdm = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",s.p_launch+10*log10(N)).process(Opt_sig_wdm);
% Opt_sig_wdm.spectrum("fignum",101,"displayname",'bla','normalizeTo0dB',0,'lambda0_nm',1310,'useWavelengthAxis',0);
% Opt_sig_wdm.spectrum("fignum",101,"displayname",'bla','normalizeTo0dB',1,'max_num_lines',2);
%%%%%% Fiber %%%%%%
Opt_sig_wdm_fib=Opt_sig_wdm;
segment_length = 1; % km
nSegments = link_length/segment_length;
zdw = 1310;
D_local = 0; %if ~=0, simulation uses "segmented fiber with d+,d-)
d_local = d_local; %if ~=0, simulation uses "segmented fiber with d+,d-)
randomize_D = true;
Dvec = getDispersionVector(nSegments, D_local, zdw, randomize_D, s.random_key+realiz);
Dvec = getDispersionVector(nSegments, d_local, zdw, randomize_D, s.random_key+realiz);
propdist = 0;
for seg = 1:nSegments
Opt_sig_wdm_fib = DP_Fiber("L",segment_length,"D",Dvec(s),"Dpmd",pmd,"Ds",0.07,...
"beat_len",10,"corr_len",100,"dz",1,"manakov",0,...
"gamma",s.gamma,"lambda",zdw,"n_waveplates",10,"SS_dphimax",0.01,...
"SS_dzmax",50,"SS_dzmin",10,"X_alpha",0.3,"X_beta",0,"rng",1).process(Opt_sig_wdm_fib);
Opt_sig_wdm_fib = DP_Fiber("L",segment_length,"D",Dvec(seg),"Dpmd",s.pmd,"Ds",0.07,...
"beat_len",10,"corr_len",100,"dz",1,"manakov",0,...
"gamma",s.gamma,"lambda",zdw,"n_waveplates",10,"SS_dphimax",0.01,...
"SS_dzmax",50,"SS_dzmin",10,"X_alpha",0.3,"X_beta",0,"rng",1).process(Opt_sig_wdm_fib);
propdist = segment_length;
end
Opt_sig_wdm_fib.spectrum("fignum",realiz,"displayname",'bla','lambda0_nm',1310,'useWavelengthAxis',0);
% Opt_sig_wdm_fib.move_it_spectrum("fignum",100212,"displayname",'bla');
% Opt_sig = Fiber("fsimu",Opt_sig.fs,"fiber_length",s.link_length/1000,"alpha",0.3,"D",0,"lambda0",1310,"s.gamma",0,"Dslope",0.07).process(Opt_sig)
%%%%%% Demux after 2 km %%%%%%
Opt_sig_wdm_demux = Optical_Demultiplex("attenuation",0,"B",200e9,"filtype",1,"fs_out",fdac*kover,"fs_in",fdac*kover*upsample_pow,"lambda_center",1310).process(Opt_sig_wdm_fib);
for ri = 1:length(s.rop)
%%%%%% ROP %%%%%%
Opt_sig_wdm_rx = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",s.rop(ri)+10*log10(N)).process(Opt_sig_wdm_fib);
Opt_sig_wdm_demux = Optical_Demultiplex("attenuation",0,"B",200e9,"filtype",1,"fs_out",Opt_sig_wdm_rx.fs/upsample_pow,"fs_in",Opt_sig_wdm_rx.fs,"lambda_center",1310).process(Opt_sig_wdm_rx);
PD_cell = {};
for l = 1:N
parfor l = 1:N
%%%%%% ROP %%%%%%
Opt_sig_wdm_rx = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",s.rop(ri)).process(Opt_sig_wdm_demux{l}); % rop+10*log10(N)
%%%%%% PD Square Law %%%%%%
assert(fdac*kover==Opt_sig_wdm_demux{l}.fs,'Sampling Frequencies do not match! Check previous steps');
PD_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20,"nep",1.8e-11,"randomkey",s.random_key+l+realiz).process(Opt_sig_wdm_demux{l});
assert(fdac*kover==Opt_sig_wdm_rx.fs,'Sampling Frequencies do not match! Check previous steps');
PD_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20,"nep",1.8e-11,"randomkey",s.random_key+l+realiz).process(Opt_sig_wdm_rx);
% PD_sig.spectrum("fignum",222,"displayname",'bla','normalizeTo0dB',1);
%%%%%% Low-pass RX (PD, El. Connectors and Scope %%%%%%
rx_bwl = 100e9;
PD_sig = Filter('filtdegree',4,"f_cutoff",rx_bwl,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(PD_sig);
% %%%%%% Low-pass Scope %%%%%%
Lp_scpe = Filter('filtdegree',4,"f_cutoff",110e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true);
Lp_scpe = Filter('filtdegree',4,"f_cutoff",80e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true);
%%%%%% Scope %%%%%%
Scpe_sig = Scope("fsimu",fdac*kover,"fadc",fadc,...
"delay",0,"fixed_delay",0,"filtertype",filtertypes.butterworth,...
"samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,...
"adcresolution",8,"quantbuffer",0.1,'block_dc',1,'lpf_active',1,'H_lpf',Lp_scpe).process(PD_sig);
"adcresolution",8,"quantbuffer",0.1,'block_dc',1,'lpf_active',0,'H_lpf',Lp_scpe).process(PD_sig);
Scpe_sig_2sps = Scpe_sig.resample("fs_out",2*fsym);
% Scpe_sig.spectrum("fignum",222,"displayname",'bla','normalizeTo0dB',1);
[~, Scpe_cell, ~, found_sync] = Scpe_sig_2sps.tsynch("reference", Symbols{l}, "fs_ref", fsym, "debug_plots", 1);
Rx_sig = Scpe_cell{1};
Rx_sig = Rx_sig.normalize("mode","rms");
% FFE
@@ -227,7 +319,7 @@ for realiz = 1:s.num_realiz
output_dbt{l,ri,realiz} = dbt_results;
end
end
res = struct();
@@ -237,6 +329,10 @@ for realiz = 1:s.num_realiz
res.mlse = output_mlse;
res.dbt = output_dbt;
%%%%%% Demux after final (10) km %%%%%%
% Save results
save(fullfile(output_root, fname), 'res', '-v7.3');
fprintf('Saved results to: %s\n', fullfile(output_root, fname));
@@ -244,36 +340,7 @@ for realiz = 1:s.num_realiz
end
function dispersion_vector = getDispersionVector(N, D, ref_zdw, randomize_ZDW, randomkey)
% s.MATLAB version of the Python generator shown above.
% Returns an N×1 vector (ps/(nm·km)).
%
% D is the nominal dispersion magnitude. For D>0 the link is segmented with
% alternating sign (+D, -D, +D, …). For D==0 it is flat (0) except for
% ZDW randomization. The ZDW detuning is ~N(0, 2 nm) around 1310 nm and is
% converted to dispersion via 0.09 ps/(nm·km) per nm.
% constants (matching the Python code)
meanLambda_nm = 1310; % center wavelength
sigma_nm = 2; % ZDW sigma
Dslope = 0.07; % ps/(nm·km) per nm detuning
% random ZDW-induced dispersion offset
if randomize_ZDW
rng(randomkey, 'twister');
rand_zdws_nm = meanLambda_nm + sigma_nm .* randn(N,1);
rand_D = (rand_zdws_nm - ref_zdw) .* Dslope; % ps/(nm·km)
else
rand_D = zeros(N,1);
end
% nominal segmented pattern (match Python intent; keep length N)
if D > 0
base = (-1) .^ ((0:N-1).'); % +1,-1,+1,-1,...
else % D == 0 (or anything else)
base = ones(N,1);
end
dispersion_vector = base .* D + rand_D; % ps/(nm·km)
end

379
projects/WDM/WDM_model_10km.m Normal file
View File

@@ -0,0 +1,379 @@
function WDM_model_10km(options)
%%% Run parameters
% TX
arguments
options.num_channels = 16;
options.channel_spacing = 400e9;
options.fiber_length_km = 0;
options.rand_key = 1;
options.num_realiz = 1;
options.fwm_mitigation_technique = "co";
end
%% Add the imdd_simulation framework to the path
if ispc
addpath(genpath('C:\Users\Silas\Documents\MATLAB\imdd_simulation'));
else
% Linux path on the cluster
addpath(genpath('/work_beegfs/sutef391/imdd_simulation'));
end
% Quiet the ambiguous CET warning (best is to set TZ in sbatch)
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
% Use a per-job, node-local JobStorageLocation to avoid stale locks on $HOME
tmpbase = getenv('TMPDIR');
if isempty(tmpbase), tmpbase = tempdir; end
jsl = fullfile(tmpbase, sprintf('matlab_jobstorage_%s_%s', ...
getenv('USER'), getenv('SLURM_JOB_ID')));
if ~exist(jsl,'dir'); mkdir(jsl); end
if 0
% Configure the local cluster explicitly and start the pool
c = parcluster('local');
c.NumWorkers = cpus;
c.JobStorageLocation = jsl;
p = gcp('nocreate');
if isempty(p) || p.NumWorkers ~= cpus
if ~isempty(p), delete(p); end
p = parpool(c, cpus);
end
fprintf('parpool up with %d workers; JobStorage=%s\n', p.NumWorkers, c.JobStorageLocation);
end
%% result filename (timestamp + optional job id)
t = datetime('now','TimeZone','local','Format','yyyyMMdd_HHmmss');
jobid = getenv('SLURM_JOB_ID'); if isempty(jobid), jobid = 'nojid'; end
host = getenv('HOSTNAME'); if isempty(host), host = 'localhost'; end
% Output directory depends on platform
foldname = sprintf('%dkm_%dch_%dghz_%s', options.fiber_length_km(end), options.num_channels, options.channel_spacing.*1e-9, options.fwm_mitigation_technique);
if ispc
output_root = fullfile('C:\Users\Silas\Documents\MATLAB\Datensätze\FWM_2025\',foldname,'\');
else
output_root = fullfile('/work_beegfs/sutef391/results_WDM',foldname,'\');
end
if ~exist(output_root,'dir'), mkdir(output_root); end
fname = sprintf('WDM_%s_%s_%s_%dkm_%dch_%dghz_%s.mat', char(t), host, jobid, options.fiber_length_km(end), options.num_channels, options.channel_spacing.*1e-9, options.fwm_mitigation_technique);
%% Settings
s.num_realiz = options.num_realiz;
s.wavelengthplan = calcWavelengthPlan(options.num_channels, options.channel_spacing, 1310);
link_length = options.fiber_length_km;
s.pmd = 0;%0.1;
s.gamma = 0;%0.0023;
s.M = 4;
fsym = 112e9;
fdac = 2*fsym;
fadc = 120000000000;
s.random_key = options.rand_key;
% Laser / s.Modulator
vbias_rel = 0.5;
u_pi = 4.6;
vbias = -vbias_rel*u_pi;
laser_linewidth = 0e6;
% EQ SETTINGS
vnle_order1 = 50;
vnle_order2 = 3;
vnle_order3 = 3;
vnle_order = [vnle_order1,vnle_order2,vnle_order3];
dfe_order = [0 0 0];
len_tr = 4096*2;
mu_ffe1 = 0.0001;
mu_ffe2 = 0.0008;
mu_ffe3 = 0.001;
mu_dc = 0.005;
mu_ffe = [mu_ffe1 mu_ffe3 mu_ffe3];
mu_dfe = 0.0004;
% DB Stuff
db_precode = 0;
db_encode = 0;
duob_mode = db_mode.no_db;
apply_pulsef = 0;
rcalpha = 0.05;
Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",16,"alpha",rcalpha);
N = numel(s.wavelengthplan);
f_plan = physconst('lightspeed')./(s.wavelengthplan.*1e-9);
margin = 2e12; % some THz left and right
f_span = (max(f_plan)+margin)-(min(f_plan)-margin);
f_nyq = f_span/2;
kover = 4;
upsample_required = f_nyq./(fdac*kover/2);
upsample_pow = 2^nextpow2(upsample_required);
upsample_ceil = ceil(upsample_required); %#ok<NASGU>
s.f_opt = fdac*kover*upsample_pow;
s.f_opt_nyq = s.f_opt/2;
signal_cell = {};
Symbols = {};
Tx_bits = {};
s.rop = -12;%:0.75:-0;
%% ---------- Intermediate evaluation points (distance dimension) ----------
% Evaluate BER at these intermediate distances (km), plus always include the final link_length if > 0.
segment_length = 1; % km (must match fiber loop below)
eval_dist_km = [10];
eval_dist_km = eval_dist_km(eval_dist_km <= link_length);
% Always include final distance (if > 0) and avoid duplicates
if link_length > 0
if isempty(eval_dist_km) || eval_dist_km(end) ~= link_length
eval_dist_km = unique([eval_dist_km link_length], 'stable');
end
end
% Convert to segment indices; require integer multiples of segment_length
eval_seg = eval_dist_km ./ segment_length;
if any(abs(eval_seg - round(eval_seg)) > 1e-12)
error('eval_dist_km must be integer multiples of segment_length=%g km.', segment_length);
end
eval_seg = round(eval_seg);
nEval = numel(eval_seg);
% -------------------------------------------------------------------------
%% Preallocate outputs (add eval distance as 4th dimension)
output_ffe = cell(length(s.wavelengthplan), length(s.rop), s.num_realiz, nEval);
output_dfe = cell(length(s.wavelengthplan), length(s.rop), s.num_realiz, nEval);
output_vnle = cell(length(s.wavelengthplan), length(s.rop), s.num_realiz, nEval);
output_mlse = cell(length(s.wavelengthplan), length(s.rop), s.num_realiz, nEval);
output_dbt = cell(length(s.wavelengthplan), length(s.rop), s.num_realiz, nEval);
s.p_launch = 3;
s.p = options.fwm_mitigation_technique;
switch s.p
case "co"
pol_rot = 100.*ones(1,length(s.wavelengthplan));
d_local = 0;
case "pair"
pol_rot = repmat([100,100,0,0],1,length(s.wavelengthplan)/4);
d_local = 0;
case "alt"
pol_rot = repmat([100,0,100,0],1,length(s.wavelengthplan)/4);
d_local = 0;
case "seg"
pol_rot = 100.*ones(1,length(s.wavelengthplan));
d_local = 3;
otherwise
error('Unknown fwm_mitigation_technique: %s', string(s.p));
end
for realiz = 1:s.num_realiz
% Reset per-realization storage (so each realiz writes only its slice)
signal_cell = cell(1,N);
Symbols = cell(1,N);
Tx_bits = cell(1,N);
%% ---------- TX per channel ----------
for l = 1:N
[Digi_sig,Symbols{l},Tx_bits{l}] = PAMsource( ...
"fsym",fsym,"M",s.M,"order",18,"useprbs",0, ...
"fs_out",fdac, ...
"applyclipping",0,"clipfactor",1.5, ...
"applypulseform",apply_pulsef,"pulseformer",Pform, ...
"randkey",s.random_key+l+realiz, ...
"db_precode",db_precode,"db_encode",db_encode, ...
"mrds_code",0,"mrds_blocklength",512,"duobinary_mode",duob_mode ...
).process();
Lp_awg = Filter('filtdegree',3,"f_cutoff",56e9,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true);
El_sig = AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",1,"kover",kover, ...
"bit_resolution",6,"upsampling_method","samplehold","precomp_sinc_rolloff",0, ...
"H_lpf",Lp_awg,"dac_max",0.6,"dac_min",-0.6).process(Digi_sig);
% Electrical Driver Amplifier
El_sig = El_sig.normalize("mode","oneone");
scaling = 0.6*(u_pi/2-abs(vbias-u_pi/2));
El_sig = El_sig .* scaling;
% E/O Conversion
Eml_out = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs, ...
"lambda",s.wavelengthplan(l),"bias",vbias,"u_pi",u_pi, ...
"linewidth",laser_linewidth,"randomkey",s.random_key+l+realiz,"alpha",0.8).process(El_sig);
signal_cell{l} = Polarization_Controller("mode","rot_power","desired_power",pol_rot(l)).process(Eml_out);
end
%% ---------- WDM mux + launch ----------
Opt_sig_wdm = Optical_Multiplex("fs_in",fdac*kover,"fs_out",upsample_pow*fdac*kover, ...
"lambda_center",1310,"random_key",0,"filtype",1,"B",120e9).process(signal_cell);
Opt_sig_wdm = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power", ...
"amplification_db",s.p_launch+10*log10(N)).process(Opt_sig_wdm);
%% ---------- Fiber propagation ----------
Opt_sig_wdm_fib = Opt_sig_wdm;
nSegments = link_length/segment_length;
if abs(nSegments - round(nSegments)) > 1e-12
error('fiber_length_km=%g must be an integer multiple of segment_length=%g km.', link_length, segment_length);
end
nSegments = round(nSegments);
zdw = 1310;
randomize_D = true;
Dvec = getDispersionVector(nSegments, d_local, zdw, randomize_D, s.random_key+realiz);
eval_ptr = 1; % points into eval_seg
for seg = 1:nSegments
fprintf('Realiz %d/%d: Segment %d/%d \n', realiz, s.num_realiz, seg, nSegments);
Opt_sig_wdm_fib = DP_Fiber("L",segment_length,"D",Dvec(seg),"Dpmd",s.pmd,"Ds",0.07, ...
"beat_len",10,"corr_len",100,"dz",1,"manakov",0, ...
"gamma",s.gamma,"lambda",zdw,"n_waveplates",10,"SS_dphimax",0.01, ...
"SS_dzmax",50,"SS_dzmin",10,"X_alpha",0.3,"X_beta",0,"rng",1).process(Opt_sig_wdm_fib);
% -------- Evaluate at intermediate distance (if scheduled) --------
if eval_ptr <= nEval && seg == eval_seg(eval_ptr)
%%%%%% Demux at this distance %%%%%%
Opt_sig_wdm_demux = Optical_Demultiplex("attenuation",0,"B",200e9,"filtype",1, ...
"fs_out",fdac*kover,"fs_in",fdac*kover*upsample_pow,"lambda_center",1310).process(Opt_sig_wdm_fib);
for ri = 1:length(s.rop)
for l = 1:N
%%%%%% ROP %%%%%%
Opt_sig_wdm_rx = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power", ...
"amplification_db",s.rop(ri)).process(Opt_sig_wdm_demux{l});
%%%%%% PD Square Law %%%%%%
assert(fdac*kover==Opt_sig_wdm_rx.fs,'Sampling Frequencies do not match! Check previous steps');
PD_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20, ...
"nep",1.8e-11,"randomkey",s.random_key+l+realiz).process(Opt_sig_wdm_rx);
%%%%%% Low-pass RX (PD, El. Connectors and Scope) %%%%%%
rx_bwl = 100e9;
PD_sig = Filter('filtdegree',4,"f_cutoff",rx_bwl,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(PD_sig);
%%%%%% Low-pass Scope %%%%%%
Lp_scpe = Filter('filtdegree',4,"f_cutoff",80e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true);
%%%%%% Scope %%%%%%
Scpe_sig = Scope("fsimu",fdac*kover,"fadc",fadc, ...
"delay",0,"fixed_delay",0,"filtertype",filtertypes.butterworth, ...
"samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0, ...
"adcresolution",8,"quantbuffer",0.1,'block_dc',1,'lpf_active',0,'H_lpf',Lp_scpe).process(PD_sig);
Scpe_sig_2sps = Scpe_sig.resample("fs_out",2*fsym);
[~, Scpe_cell, ~, found_sync] = Scpe_sig_2sps.tsynch("reference", Symbols{l}, "fs_ref", fsym, "debug_plots", 1); %#ok<ASGLU>
Rx_sig = Scpe_cell{1};
Rx_sig = Rx_sig.normalize("mode","rms");
% -------------------- FFE --------------------
ffe_order = [50, 0, 0];
eq_ffe = EQ("Ne",ffe_order,"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",0);
ffe_results = ffe(eq_ffe,s.M,Rx_sig,Symbols{l},Tx_bits{l}, ...
"precode_mode",duob_mode, ...
'showAnalysis',0, ...
"postFFE",[], ...
"eth_style_symbol_mapping",0);
output_ffe{l,ri,realiz,eval_ptr} = ffe_results;
% -------------------- DFE --------------------
dfe_order = [50, 0, 0];
eq_dfe = EQ("Ne",dfe_order,"Nb",[2,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",0);
dfe_results = ffe(eq_dfe,s.M,Rx_sig,Symbols{l},Tx_bits{l}, ...
"precode_mode",duob_mode, ...
'showAnalysis',0, ...
"postFFE",[], ...
"eth_style_symbol_mapping",0);
output_dfe{l,ri,realiz,eval_ptr} = dfe_results;
% -------------------- VNLE + MLSE --------------------
pf_ncoeffs = 1;
ffe_order = [50, 5, 5];
eq_ = EQ("Ne",ffe_order,"Nb",dfe_order,"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);
pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1);
useviterbi = 0;
if useviterbi
mlse_ = MLSE_viterbi("duobinary_output",0,'M',s.M,'trellis_states',PAMmapper(s.M,0).levels);
else
mlse_ = MLSE("duobinary_output",0,'M',s.M,'trellis_states',PAMmapper(s.M,0).levels);
end
[vnle_results, mlse_results] = vnle_postfilter_mlse(eq_, pf_, mlse_, s.M, Rx_sig, Symbols{l},Tx_bits{l}, ...
"precode_mode", duob_mode, ...
'showAnalysis', 0, ...
"postFFE", [], ...
"eth_style_symbol_mapping", 0);
output_vnle{l,ri,realiz,eval_ptr} = vnle_results;
output_mlse{l,ri,realiz,eval_ptr} = mlse_results;
% -------------------- DB target --------------------
useviterbi = 0;
if useviterbi
mlse_db_ = MLSE_viterbi("duobinary_output",0,'M',s.M,'trellis_states',PAMmapper(s.M,0).levels);
else
mlse_db_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",s.M,'trellis_states',PAMmapper(s.M,0).levels);
end
ffe_order = [50, 5, 5];
eq_ = EQ("Ne",ffe_order,"Nb",dfe_order,"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_, s.M, Rx_sig, Symbols{l},Tx_bits{l}, ...
"precode_mode", duob_mode, ...
'showAnalysis', 0, ...
"postFFE", []);
output_dbt{l,ri,realiz,eval_ptr} = dbt_results;
end
end
eval_ptr = eval_ptr + 1;
end
% ----------------------------------------------------------------
end
%% Save results (per realization)
res = struct();
res.settings = s;
res.eval_dist_km = eval_dist_km;
res.ffe = output_ffe;
res.dfe = output_dfe;
res.vnle = output_vnle;
res.mlse = output_mlse;
res.dbt = output_dbt;
save(fullfile(output_root, fname), 'res', '-v7.3');
fprintf('Saved results to: %s\n', fullfile(output_root, fname));
disp(datetime('now','TimeZone','local','Format','yyyyMs.Mdd_HHmmss'));
end
end

900
projects/WDM/WDM_model_10km_queue.m Normal file
View File

@@ -0,0 +1,900 @@
function WDM_model_10km_queue(options)
%WDM_model_10km_queue_memopt
% Queue/pipeline version with:
% (a) reduced variable lifetime / fewer unnecessary copies (clear large temporaries early)
% (b) worker-side memory monitoring (RSS logging to per-worker files)
%%% Run parameters
arguments
options.num_channels = 8;
options.channel_spacing = 400e9;
options.fiber_length_km = 10;
options.rand_key = 1;
options.num_realiz = 1;
options.fwm_mitigation_technique = "co";
options.chirpalpha = 0;
end
%% Add the imdd_simulation framework to the path
if ispc
addpath(genpath('C:\Users\Silas\Documents\MATLAB\imdd_simulation'));
else
addpath(genpath('/work_beegfs/sutef391/imdd_simulation'));
end
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'))-1; end
% Use a per-job, node-local JobStorageLocation to avoid stale locks on $HOME
tmpbase = getenv('TMPDIR');
if isempty(tmpbase), tmpbase = tempdir; end
jsl = fullfile(tmpbase, sprintf('matlab_jobstorage_%s_%s', ...
getenv('USER'), getenv('SLURM_JOB_ID')));
if ~exist(jsl,'dir'); mkdir(jsl); end
% Start pool once
c = parcluster('local');
c.NumWorkers = cpus;
c.JobStorageLocation = jsl;
p = gcp('nocreate');
if isempty(p) || p.NumWorkers ~= cpus
if ~isempty(p), delete(p); end
p = parpool(c, cpus,"IdleTimeout",300);
end
fprintf('parpool up with %d workers; JobStorage=%s\n', p.NumWorkers, c.JobStorageLocation);
%% result filename (timestamp + optional job id)
t = datetime('now','TimeZone','local','Format','yyyyMMdd_HHmmss');
jobid = getenv('SLURM_JOB_ID'); if isempty(jobid), jobid = 'nojid'; end
host = getenv('HOSTNAME'); if isempty(host), host = 'localhost'; end
foldname = sprintf('%dkm_%dch_%dghz_%s', options.fiber_length_km(end), options.num_channels, options.channel_spacing.*1e-9, options.fwm_mitigation_technique);
if ispc
output_root = fullfile('C:\Users\Silas\Documents\MATLAB\Datensätze\FWM_2025\',foldname,'\');
else
output_root = fullfile('/work_beegfs/sutef391/results_WDM',foldname,'\');
end
if ~exist(output_root,'dir'), mkdir(output_root); end
fname = sprintf('WDM_%s_%s_%s_%dkm_%dch_%dghz_%s_alpha%0.1f', char(t), host, jobid, options.fiber_length_km(end), options.num_channels, options.channel_spacing.*1e-9, options.fwm_mitigation_technique, options.chirpalpha);
fname = strrep(fname,'.','_');
fname = [fname,'.mat'];
% Worker memory logs directory
memlog_dir = fullfile(output_root, 'memlogs');
if ~exist(memlog_dir,'dir'), mkdir(memlog_dir); end
%% Settings
s.num_realiz = options.num_realiz;
s.wavelengthplan = calcWavelengthPlan(options.num_channels, options.channel_spacing, 1310);
link_length = options.fiber_length_km;
s.pmd = 0;%0.1;
s.gamma = 0;%0.0023;
s.M = 4;
fsym = 112e9;
fdac = 2*fsym;
fadc = 120000000000;
s.random_key = options.rand_key;
% Laser / s.Modulator
vbias_rel = 0.5;
u_pi = 4.6;
vbias = -vbias_rel*u_pi;
laser_linewidth = 0e6;
% EQ SETTINGS
dfe_order = [0 0 0];
len_tr = 4096*2;
mu_ffe1 = 0.0001;
mu_ffe2 = 0.0008;
mu_ffe3 = 0.001;
mu_dc = 0.005;
mu_ffe = [mu_ffe1 mu_ffe3 mu_ffe3];
mu_dfe = 0.0004;
% DB Stuff
db_precode = 0;
db_encode = 0;
duob_mode = db_mode.no_db;
apply_pulsef = 0;
rcalpha = 0.05;
Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",16,"alpha",rcalpha);
N = numel(s.wavelengthplan);
f_plan = physconst('lightspeed')./(s.wavelengthplan.*1e-9);
margin = 2e12; % some THz left and right
f_span = (max(f_plan)+margin)-(min(f_plan)-margin);
f_nyq = f_span/2;
kover = 4;
upsample_required = f_nyq./(fdac*kover/2);
upsample_pow = 2^nextpow2(upsample_required);
s.f_opt = fdac*kover*upsample_pow;
s.f_opt_nyq = s.f_opt/2;
s.rop = -10:1:0;
%% ---------- Intermediate evaluation points (distance dimension) ----------
segment_length = 1; % km (must match fiber loop below)
eval_dist_km = [2,4,6,8,10];
eval_dist_km = eval_dist_km(eval_dist_km <= link_length);
if link_length > 0
if isempty(eval_dist_km) || eval_dist_km(end) ~= link_length
eval_dist_km = unique([eval_dist_km link_length], 'stable');
end
end
eval_seg = eval_dist_km ./ segment_length;
if any(abs(eval_seg - round(eval_seg)) > 1e-12)
error('eval_dist_km must be integer multiples of segment_length=%g km.', segment_length);
end
eval_seg = round(eval_seg);
nEval = numel(eval_seg);
% -------------------------------------------------------------------------
%% Preallocate outputs (eval distance as 4th dimension)
output_ffe = cell(N, length(s.rop), s.num_realiz, nEval);
output_dfe = cell(N, length(s.rop), s.num_realiz, nEval);
output_vnle = cell(N, length(s.rop), s.num_realiz, nEval);
output_mlse = cell(N, length(s.rop), s.num_realiz, nEval);
output_dbt = cell(N, length(s.rop), s.num_realiz, nEval);
s.p_launch = 3;
s.p = options.fwm_mitigation_technique;
switch s.p
case "co"
pol_rot = 100.*ones(1,N);
d_local = 0;
case "pair"
pol_rot = repmat([100,100,0,0],1,N/4);
d_local = 0;
case "alt"
pol_rot = repmat([100,0,100,0],1,N/4);
d_local = 0;
case "seg"
pol_rot = 100.*ones(1,N);
d_local = 3;
otherwise
error('Unknown fwm_mitigation_technique: %s', string(s.p));
end
for realiz = 1:s.num_realiz
% Per-realization TX storage (needed later for DSP: Symbols/Tx_bits)
signal_cell = cell(1,N);
Symbols = cell(1,N);
Tx_bits = cell(1,N);
% -------- Job queue containers --------
F = parallel.FevalFuture.empty(0,1);
meta = struct('l',{},'ri',{},'realiz',{},'eval_ptr',{});
% -------------------------------------
%% ---------- TX per channel ----------
for l = 1:N
[Digi_sig,Symbols{l},Tx_bits{l}] = PAMsource( ...
"fsym",fsym,"M",s.M,"order",15,"useprbs",0, ...
"fs_out",fdac, ...
"applyclipping",0,"clipfactor",1.5, ...
"applypulseform",apply_pulsef,"pulseformer",Pform, ...
"randkey",s.random_key+l+realiz, ...
"db_precode",db_precode,"db_encode",db_encode, ...
"mrds_code",0,"mrds_blocklength",512,"duobinary_mode",duob_mode ...
).process();
Lp_awg = Filter('filtdegree',3,"f_cutoff",56e9,"fs",fdac*kover, ...
"filterType",filtertypes.gaussian,"active",true);
El_sig = AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",1,"kover",kover, ...
"bit_resolution",6,"upsampling_method","samplehold","precomp_sinc_rolloff",0, ...
"H_lpf",Lp_awg,"dac_max",0.6,"dac_min",-0.6).process(Digi_sig);
% Digi_sig not needed after AWG
clear Digi_sig
% Electrical Driver Amplifier
El_sig = El_sig.normalize("mode","oneone");
scaling = 0.6*(u_pi/2-abs(vbias-u_pi/2));
El_sig = El_sig .* scaling;
% E/O Conversion
Eml_out = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs, ...
"lambda",s.wavelengthplan(l),"bias",vbias,"u_pi",u_pi, ...
"linewidth",laser_linewidth,"randomkey",s.random_key+l+realiz,"alpha",options.chirpalpha).process(El_sig);
s.alpha = options.chirpalpha;
% El_sig not needed after EML
clear El_sig
signal_cell{l} = Polarization_Controller("mode","rot_power","desired_power",pol_rot(l)).process(Eml_out);
% Eml_out not needed after pol controller
clear Eml_out Lp_awg
end
disp('Signal generated for all channels.');
%% ---------- WDM mux + launch ----------
Opt_sig_wdm = Optical_Multiplex("fs_in",fdac*kover,"fs_out",upsample_pow*fdac*kover, ...
"lambda_center",1310,"random_key",0,"filtype",1,"B",120e9).process(signal_cell);
% IMPORTANT: signal_cell is not needed anymore after multiplex (Symbols/Tx_bits remain)
clear signal_cell
Opt_sig_wdm = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power", ...
"amplification_db",s.p_launch+10*log10(N)).process(Opt_sig_wdm);
%% ---------- Fiber propagation ----------
Opt_sig_wdm_fib = Opt_sig_wdm;
% Opt_sig_wdm no longer needed as separate handle/object
clear Opt_sig_wdm
nSegments = link_length/segment_length;
if abs(nSegments - round(nSegments)) > 1e-12
error('fiber_length_km=%g must be an integer multiple of segment_length=%g km.', link_length, segment_length);
end
nSegments = round(nSegments);
zdw = 1310;
randomize_D = false;
% Guard for 0 km: avoid calling getDispersionVector(0,...) if it doesn't support it
if nSegments > 0
Dvec = getDispersionVector(nSegments, d_local, zdw, randomize_D, s.random_key+realiz);
else
Dvec = [];
end
eval_ptr = 1;
% =================== Queue throttle (prevents OOM) ===================
% Limit how many futures can be outstanding (running + queued + finished-not-yet-fetched).
maxInFlight = max(2, p.NumWorkers);
% =====================================================================
% -------- Evaluate at 0 km (BTB) if requested --------
if eval_ptr <= nEval && eval_seg(eval_ptr) == 0
disp('0 km before demux.');
Opt_sig_wdm_demux = Optical_Demultiplex("attenuation",0,"B",200e9,"filtype",1, ...
"fs_out",fdac*kover,"fs_in",fdac*kover*upsample_pow,"lambda_center",1310).process(Opt_sig_wdm_fib);
cnt = 0;
total = length(s.rop)*N;
fprintf('total of %d jobs to enqueue at 0 km\n', total);
for ri = 1:length(s.rop)
for l = 1:N
cnt = cnt + 1;
% ---- Throttle before enqueueing more futures ----
[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);
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, ...
'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);
end
end
clear Opt_sig_wdm_demux
eval_ptr = eval_ptr + 1;
end
% -----------------------------------------------------
for seg = 1:nSegments
fprintf('Realiz %d/%d: Segment %d/%d \n', realiz, s.num_realiz, seg, nSegments);
Opt_sig_wdm_fib = DP_Fiber("L",segment_length,"D",Dvec(seg),"Dpmd",s.pmd,"Ds",0.07, ...
"beat_len",10,"corr_len",100,"dz",1,"manakov",0, ...
"gamma",s.gamma,"lambda",zdw,"n_waveplates",10,"SS_dphimax",0.01, ...
"SS_dzmax",50,"SS_dzmin",10,"X_alpha",0.3,"X_beta",0,"rng",1).process(Opt_sig_wdm_fib);
% -------- Evaluate at intermediate distance (enqueue jobs) --------
if eval_ptr <= nEval && seg == eval_seg(eval_ptr)
Opt_sig_wdm_demux = Optical_Demultiplex("attenuation",0,"B",200e9,"filtype",1, ...
"fs_out",fdac*kover,"fs_in",fdac*kover*upsample_pow,"lambda_center",1310).process(Opt_sig_wdm_fib);
cnt = 0;
total = length(s.rop)*N;
for ri = 1:length(s.rop)
for l = 1:N
cnt = cnt + 1;
% ---- Throttle before enqueueing more futures ----
[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);
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, ...
'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);
end
end
clear Opt_sig_wdm_demux
eval_ptr = eval_ptr + 1;
end
% ----------------------------------------------------------------
% Non-blocking harvest (your existing line)
[F, meta, output_ffe, output_dfe, output_vnle, output_mlse, output_dbt] = ...
collect_done(F, meta, output_ffe, output_dfe, output_vnle, output_mlse, output_dbt, false);
end
drainIterMax = 10; % 10 * 60s = 10 minutes
for drainIter = 1:drainIterMax
if isempty(F)
break;
end
% Status line: how many futures are in which states?
st = {F.State};
nUnavail = sum(strcmp(st,'unavailable'));
nFinished = sum(strcmp(st,'finished'));
nRunning = sum(strcmp(st,'running'));
nQueued = sum(strcmp(st,'queued'));
nFailed = sum(strcmp(st,'failed'));
fprintf('[drain %d/%d] F=%d | finished=%d running=%d queued=%d failed=%d unavailable=%d\n', ...
drainIter, drainIterMax, numel(F), nFinished, nRunning, nQueued, nFailed, nUnavail);
% Try to fetch at least one result (blocking)
[F, meta, output_ffe, output_dfe, output_vnle, output_mlse, output_dbt] = ...
collect_done(F, meta, output_ffe, output_dfe, output_vnle, output_mlse, output_dbt, true);
% If still not empty, wait a bit before the next drain attempt
if ~isempty(F)
pause(60);
end
end
save_results_per_realization( ...
s, eval_dist_km, ...
output_ffe, output_dfe, output_vnle, output_mlse, output_dbt, ...
output_root, fname);
% Per-realization large arrays that are no longer needed
clear Opt_sig_wdm_fib Symbols Tx_bits Dvec
end
end % end main
%% ========================= Local functions =========================
function [ffe_results, dfe_results, vnle_results, mlse_results, dbt_results] = rx_job( ...
Opt_sig_chan, rop_db, 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)
% NOTE: keep plotting OFF in workers
debug_plots = 0;
% Create per-worker logfile (avoid contention)
logfile = make_worker_logfile(memlog_dir);
log_mem(logfile, 'job_start', l, ri, realiz, eval_ptr, rop_db, Opt_sig_chan);
%%%%%% ROP %%%%%%
Opt_sig_wdm_rx = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power", ...
"amplification_db", rop_db).process(Opt_sig_chan);
% Opt_sig_chan no longer needed after amp
clear Opt_sig_chan
log_mem(logfile, 'after_amp', l, ri, realiz, eval_ptr, rop_db, Opt_sig_wdm_rx);
%%%%%% PD Square Law %%%%%%
assert(fdac*kover==Opt_sig_wdm_rx.fs,'Sampling Frequencies do not match! Check previous steps');
PD_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20, ...
"nep",1.8e-11,"randomkey",s.random_key + l + realiz).process(Opt_sig_wdm_rx);
% Opt_sig_wdm_rx not needed after PD
clear Opt_sig_wdm_rx
log_mem(logfile, 'after_pd', l, ri, realiz, eval_ptr, rop_db, PD_sig);
%%%%%% Low-pass RX (PD, El. Connectors and Scope) %%%%%%
rx_bwl = 100e9;
PD_sig = Filter('filtdegree',4,"f_cutoff",rx_bwl,"fs",fdac*kover, ...
"filterType",filtertypes.butterworth,"active",true).process(PD_sig);
log_mem(logfile, 'after_rx_lpf', l, ri, realiz, eval_ptr, rop_db, PD_sig);
%%%%%% Low-pass Scope %%%%%%
Lp_scpe = Filter('filtdegree',4,"f_cutoff",80e9,"fs",fadc, ...
"filterType",filtertypes.butterworth,"active",true);
%%%%%% Scope %%%%%%
Scpe_sig = Scope("fsimu",fdac*kover,"fadc",fadc, ...
"delay",0,"fixed_delay",0,"filtertype",filtertypes.butterworth, ...
"samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0, ...
"adcresolution",8,"quantbuffer",0.1,'block_dc',1,'lpf_active',0,'H_lpf',Lp_scpe).process(PD_sig);
% PD_sig no longer needed after scope
clear PD_sig Lp_scpe
log_mem(logfile, 'after_scope', l, ri, realiz, eval_ptr, rop_db, Scpe_sig);
Scpe_sig_2sps = Scpe_sig.resample("fs_out",2*fsym);
% Scpe_sig no longer needed after resample
clear Scpe_sig
log_mem(logfile, 'after_resample', l, ri, realiz, eval_ptr, rop_db, Scpe_sig_2sps);
[~, Scpe_cell, ~, ~] = Scpe_sig_2sps.tsynch("reference", Symbols_l, "fs_ref", fsym, "debug_plots", debug_plots);
% Scpe_sig_2sps no longer needed after sync
clear Scpe_sig_2sps
log_mem(logfile, 'after_sync', l, ri, realiz, eval_ptr, rop_db);
Rx_sig = Scpe_cell{1}.normalize("mode","rms");
% Scpe_cell no longer needed
clear Scpe_cell
log_mem(logfile, 'after_rxsig', l, ri, realiz, eval_ptr, rop_db, Rx_sig);
% -------------------- FFE --------------------
ffe_order = [50, 0, 0];
eq_ffe = EQ("Ne",ffe_order,"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",0);
ffe_results = ffe(eq_ffe,s.M,Rx_sig,Symbols_l,Tx_bits_l, ...
"precode_mode",duob_mode,'showAnalysis',0,"postFFE",[], ...
"eth_style_symbol_mapping",0);
clear eq_ffe
log_mem(logfile, 'after_ffe', l, ri, realiz, eval_ptr, rop_db);
% -------------------- DFE --------------------
eq_dfe = EQ("Ne",ffe_order,"Nb",[2,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",0);
dfe_results = ffe(eq_dfe,s.M,Rx_sig,Symbols_l,Tx_bits_l, ...
"precode_mode",duob_mode,'showAnalysis',0,"postFFE",[], ...
"eth_style_symbol_mapping",0);
clear eq_dfe
log_mem(logfile, 'after_dfe', l, ri, realiz, eval_ptr, rop_db);
% -------------------- VNLE + MLSE --------------------
pf_ncoeffs = 1;
ffe_order3 = [50, 5, 5];
eq_v = EQ("Ne",ffe_order3,"Nb",dfe_order, ...
"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);
pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1);
mlse_ = MLSE("duobinary_output",0,'M',s.M,'trellis_states',PAMmapper(s.M,0).levels);
[vnle_results, mlse_results] = vnle_postfilter_mlse(eq_v, pf_, mlse_, s.M, Rx_sig, Symbols_l, Tx_bits_l, ...
"precode_mode", duob_mode, 'showAnalysis', 0, "postFFE", [], "eth_style_symbol_mapping", 0);
clear eq_v pf_ mlse_
log_mem(logfile, 'after_vnle_mlse', l, ri, realiz, eval_ptr, rop_db);
% -------------------- DB target --------------------
mlse_db_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",s.M,'trellis_states',PAMmapper(s.M,0).levels);
ffe_order = [50, 5, 5];
eq_ = EQ("Ne",ffe_order,"Nb",dfe_order,"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_, s.M, Rx_sig, Symbols_l, Tx_bits_l, ...
"precode_mode", duob_mode, 'showAnalysis', 0, "postFFE", []);
clear mlse_db_
log_mem(logfile, 'job_end', l, ri, realiz, eval_ptr, rop_db);
% Rx_sig no longer needed
clear Rx_sig Symbols_l Tx_bits_l
end
function [F, meta, output_ffe, output_dfe, output_vnle, output_mlse, output_dbt] = ...
collect_done(F, meta, output_ffe, output_dfe, output_vnle, output_mlse, output_dbt, block)
if nargin < 9, block = false; end
if isempty(F), return; end
if block
timeout_first = Inf;
else
timeout_first = 0;
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;
else
timeout = 0;
end
try
[k, ffe_r, dfe_r, vnle_r, mlse_r, dbt_r] = fetchNext(F, timeout);
catch
break;
end
if isempty(k)
break;
end
m = meta(k);
output_ffe{m.l, m.ri, m.realiz, m.eval_ptr} = ffe_r;
output_dfe{m.l, m.ri, m.realiz, m.eval_ptr} = dfe_r;
output_vnle{m.l, m.ri, m.realiz, m.eval_ptr} = vnle_r;
output_mlse{m.l, m.ri, m.realiz, m.eval_ptr} = mlse_r;
output_dbt{m.l, m.ri, m.realiz, m.eval_ptr} = dbt_r;
F(k) = [];
meta(k) = [];
end
end
function save_results_per_realization( ...
s, eval_dist_km, ...
output_ffe, output_dfe, output_vnle, output_mlse, output_dbt, ...
output_root, fname)
% Assemble result struct
res = struct();
res.settings = s;
res.eval_dist_km = eval_dist_km;
res.ffe = output_ffe;
res.dfe = output_dfe;
res.vnle = output_vnle;
res.mlse = output_mlse;
res.dbt = output_dbt;
% Save (HDF5-based for large data)
save(fullfile(output_root, fname), 'res', '-v7.3');
fprintf('Saved results to: %s\n', fullfile(output_root, fname));
disp(datetime('now','TimeZone','local','Format','yyyyMs.Mdd_HHmmss'));
end
%% ========================= Memory logging helpers =========================
function logfile = make_worker_logfile(memlog_dir)
% One file per worker process to avoid write contention.
pid = get_pid_safe();
t = getCurrentTask();
if isempty(t)
tid = -1;
else
tid = t.ID;
end
logfile = fullfile(memlog_dir, sprintf('memlog_pid%d_task%d.txt', pid, tid));
end
function log_mem(logfile, tag, l, ri, realiz, eval_ptr, rop_db, varargin)
% Append one line with RSS/VmSize plus optional "largest variable" info.
% Uses /proc on Linux where available.
ts = datetime('now','TimeZone','local','Format','yyyy-MM-dd HH:mm:ss.SSS');
[rssMB, vmsMB] = proc_mem_mb();
% Optional: include size of a specific variable/object if provided
extra = "";
if ~isempty(varargin)
try
x = varargin{1}; %#ok<NASGU>
w = whos('x');
extra = sprintf(' | x_bytes=%d', w.bytes);
catch
extra = " | x_bytes=NA";
end
end
line = sprintf('%s | %s | l=%d ri=%d realiz=%d eval=%d rop=%.3f | RSS=%.1fMB Vm=%.1fMB%s\n', ...
char(ts), tag, l, ri, realiz, eval_ptr, rop_db, rssMB, vmsMB, extra);
fid = fopen(logfile, 'a');
if fid ~= -1
fwrite(fid, line);
fclose(fid);
end
end
function [rssMB, vmsMB] = proc_mem_mb()
% RSS/VmSize from /proc (Linux). Falls back to NaN if unavailable.
rssMB = NaN; vmsMB = NaN;
if isunix
try
txt = fileread('/proc/self/status');
rssMB = parse_kb_field(txt, 'VmRSS:') / 1024;
vmsMB = parse_kb_field(txt, 'VmSize:') / 1024;
return;
catch
% fall through
end
end
% Fallback (Windows): memory() sometimes works
if ispc
try
m = memory;
% MemUsedMATLAB is bytes
rssMB = double(m.MemUsedMATLAB) / 1024^2;
vmsMB = NaN;
catch
end
end
end
function kb = parse_kb_field(txt, key)
kb = NaN;
idx = strfind(txt, key);
if isempty(idx), return; end
i = idx(1) + length(key);
% read until end of line
j = find(txt(i:end)==newline, 1, 'first') + i - 2;
val = strtrim(txt(i:j));
% format: "123456 kB"
parts = split(val);
kb = str2double(parts{1});
end
function pid = get_pid_safe()
pid = -1;
try
pid = feature('getpid');
catch
% no-op
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
% 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 for at least one to finish (or time out quickly)
try
wait(F, 'finished', 1);
catch
pause(0.1);
end
% ---- 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));
% 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);
end
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
try, st{k} = Fs(k).State; catch, st{k} = ''; end
try, idbad(k) = (Fs(k).ID < 0); catch, idbad(k) = false; end
end
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
% ---- Harvest done futures ----
for jj = 1:numel(idxDone)
ii = idxDone(jj);
m = meta(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
% finished-with-error
try
err = F(ii).Error;
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
output_ffe{m.l, m.ri, m.realiz, m.eval_ptr} = ffe_r;
output_dfe{m.l, m.ri, m.realiz, m.eval_ptr} = dfe_r;
output_vnle{m.l, m.ri, m.realiz, m.eval_ptr} = vnle_r;
output_mlse{m.l, m.ri, m.realiz, m.eval_ptr} = mlse_r;
output_dbt{m.l, m.ri, m.realiz, m.eval_ptr} = dbt_r;
end
% 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

49
projects/WDM/WDM_model_16x200.m → projects/WDM/WDM_model_old.m
View File

@@ -3,23 +3,23 @@
% --- FIRST LINE: evaluate settings located beside this script ---
run(fullfile(fileparts(mfilename('fullpath')),'WDM_settings.m'));
s.num_realiz = 2;
num_realiz = 50;
s.wavelengthplan = calcWavelengthPlan(16,400e9,1310);
% wavelengthplan = [1295,1305,1315,1325];
link_length = 1;
s.pmd = 0.1;
s.gamma = 0.0023;
link_length = 10;
pmd = 0.1;
gamma = 0.0023;
s.M = 4;
m = floor(log2(s.M)*10)/10;
M = 4;
m = floor(log2(M)*10)/10;
fsym = 112e9;
fdac = 2*fsym;
fadc = 120000000000;
fadc = 2*fsym;
s.random_key = 100;
% Laser / s.Modulator
vbias_rel = 0.5;
u_pi = 4.6;
u_pi = 3.2;
vbias = -vbias_rel*u_pi;
laser_linewidth = 0e6;
@@ -96,14 +96,13 @@ for realiz = 1:s.num_realiz
"mrds_code",0,"mrds_blocklength",512,"duobinary_mode",duob_mode).process();
% Digi_sig.spectrum("fignum",101,"displayname",'bla','normalizeTo0dB',0,'lambda0_nm',1310,'useWavelengthAxis',0);
Lp_awg = Filter('filtdegree',3,"f_cutoff",56e9,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true);
El_sig = AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",1,"kover",kover,"bit_resolution",6,"upsampling_method","samplehold","precomp_sinc_rolloff",0,"H_lpf",Lp_awg,"dac_max",0.6,"dac_min",-0.6).process(Digi_sig);
Lp_awg = Filter('filtdegree',3,"f_cutoff",100e9,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true);
El_sig = AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",1,"kover",kover,"bit_resolution",12,"upsampling_method","samplehold","precomp_sinc_rolloff",0,"H_lpf",Lp_awg,"dac_max",0.6,"dac_min",-0.6).process(Digi_sig);
% El_sig = s.M8199B("kover",kover).process(Digi_sig);
% El_sig.spectrum("fignum",101,"displayname",'bla','normalizeTo0dB',0,'lambda0_nm',1310,'useWavelengthAxis',0);
%%%%% Electrical Driver Amplifier %%%%%%
El_sig = El_sig.normalize("mode","oneone");
El_sig = El_sig .* u_pi .* 0.5;
% El_sig = El_sig.setPower(1,"dBm");
% figure;histogram(El_sig.signal);
@@ -114,7 +113,7 @@ for realiz = 1:s.num_realiz
end
Opt_sig_wdm = Optical_Multiplex("fs_in",fdac*kover,"fs_out",upsample_pow*fdac*kover,...
"lambda_center",1310,"random_key",0,"filtype",1,"B",120e9).process(signal_cell);
"lambda_center",1310,"random_key",0,"filtype",1,"B",200e9).process(signal_cell);
Opt_sig_wdm = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",3+10*log10(N)).process(Opt_sig_wdm);
@@ -145,19 +144,21 @@ for realiz = 1:s.num_realiz
% Opt_sig_wdm_fib.move_it_spectrum("fignum",100212,"displayname",'bla');
% Opt_sig = Fiber("fsimu",Opt_sig.fs,"fiber_length",s.link_length/1000,"alpha",0.3,"D",0,"lambda0",1310,"s.gamma",0,"Dslope",0.07).process(Opt_sig)
Opt_sig_wdm_demux = Optical_Demultiplex("attenuation",0,"B",200e9,"filtype",1,"fs_out",fdac*kover,"fs_in",fdac*kover*upsample_pow,"lambda_center",1310).process(Opt_sig_wdm_fib);
for ri = 1:length(s.rop)
%%%%%% ROP %%%%%%
Opt_sig_wdm_rx = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",s.rop(ri)+10*log10(N)).process(Opt_sig_wdm_fib);
Opt_sig_wdm_demux = Optical_Demultiplex("attenuation",0,"B",200e9,"filtype",1,"fs_out",Opt_sig_wdm_rx.fs/upsample_pow,"fs_in",Opt_sig_wdm_rx.fs,"lambda_center",1310).process(Opt_sig_wdm_rx);
PD_cell = {};
for l = 1:N
%%%%%% ROP %%%%%%
Opt_sig_wdm_rx = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",s.rop(ri)).process(Opt_sig_wdm_demux{l}); % rop+10*log10(N)
%%%%%% PD Square Law %%%%%%
assert(fdac*kover==Opt_sig_wdm_rx.fs,'Sampling Frequencies do not match! Check previous steps');
PD_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20,"nep",1.8e-11,"randomkey",s.random_key+l+realiz).process(Opt_sig_wdm_rx);
assert(fdac*kover==Opt_sig_wdm_demux{l}.fs,'Sampling Frequencies do not match! Check previous steps');
PD_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20,"nep",1.8e-11,"randomkey",s.random_key+l+realiz).process(Opt_sig_wdm_demux{l});
% PD_sig.spectrum("fignum",222,"displayname",'bla','normalizeTo0dB',1);
%%%%%% Low-pass RX (PD, El. Connectors and Scope %%%%%%
@@ -165,13 +166,13 @@ for realiz = 1:s.num_realiz
PD_sig = Filter('filtdegree',4,"f_cutoff",rx_bwl,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(PD_sig);
% %%%%%% Low-pass Scope %%%%%%
Lp_scpe = Filter('filtdegree',4,"f_cutoff",80e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true);
Lp_scpe = Filter('filtdegree',4,"f_cutoff",110e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true);
%%%%%% Scope %%%%%%
Scpe_sig = Scope("fsimu",fdac*kover,"fadc",fadc,...
"delay",0,"fixed_delay",0,"filtertype",filtertypes.butterworth,...
"samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,...
"adcresolution",8,"quantbuffer",0.1,'block_dc',1,'lpf_active',0,'H_lpf',Lp_scpe).process(PD_sig);
"adcresolution",8,"quantbuffer",0.1,'block_dc',1,'lpf_active',1,'H_lpf',Lp_scpe).process(PD_sig);
Scpe_sig_2sps = Scpe_sig.resample("fs_out",2*fsym);
% Scpe_sig.spectrum("fignum",222,"displayname",'bla','normalizeTo0dB',1);

57
projects/WDM/WDM_settings.m
View File

@@ -2,60 +2,3 @@
% Add the imdd_simulation framework to the path
if ispc
addpath(genpath('C:\Users\Silas\Documents\MATLAB\imdd_simulation'));
else
% Linux path on the cluster
addpath(genpath('/work_beegfs/sutef391/imdd_simulation'));
end
% Quiet the ambiguous CET warning (best is to set TZ in sbatch; see below)
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
% Use a per-job, node-local JobStorageLocation to avoid stale locks on $HOME
% Prefer $TMPDIR if your cluster provides it, else tempdir().
tmpbase = getenv('TMPDIR');
if isempty(tmpbase), tmpbase = tempdir; end
jsl = fullfile(tmpbase, sprintf('matlab_jobstorage_%s_%s', ...
getenv('USER'), getenv('SLURM_JOB_ID')));
if ~exist(jsl,'dir'); mkdir(jsl); end
% Configure the local cluster explicitly and start the pool
c = parcluster('local');
c.NumWorkers = cpus;
c.JobStorageLocation = jsl;
p = gcp('nocreate');
if isempty(p) || p.NumWorkers ~= cpus
if ~isempty(p), delete(p); end
p = parpool(c, cpus); % avoids the “queued” state
end
fprintf('parpool up with %d workers; JobStorage=%s\n', p.NumWorkers, c.JobStorageLocation);
% result filename (timestamp + optional job id)
t = datetime('now','TimeZone','local','Format','yyyyMMdd_HHmmss');
jobid = getenv('SLURM_JOB_ID'); if isempty(jobid), jobid = 'nojid'; end
host = getenv('HOSTNAME'); if isempty(host), host = 'localhost'; end
% Output directory depends on platform
if ispc
output_root = fullfile('C:\Users\Silas\Documents\MATLAB\Datensätze\FWM_2025\');
else
output_root = '/work_beegfs/sutef391/results_WDM';
end
if ~exist(output_root,'dir'), mkdir(output_root); end
% Build filename
t = datetime('now','TimeZone','local','Format','yyyyMMdd_HHmmss');
jobid = getenv('SLURM_JOB_ID'); if isempty(jobid), jobid = 'nojid'; end
host = getenv('HOSTNAME'); if isempty(host), host = 'localhost'; end
fname = sprintf('WDM_%s_%s_%s.mat', char(t), host, jobid);

32
projects/WDM/getDispersionVector.m Normal file
View File

@@ -0,0 +1,32 @@
function dispersion_vector = getDispersionVector(N, D, ref_zdw, randomize_ZDW, randomkey)
% s.MATLAB version of the Python generator shown above.
% Returns an N×1 vector (ps/(nm·km)).
%
% D is the nominal dispersion magnitude. For D>0 the link is segmented with
% alternating sign (+D, -D, +D, …). For D==0 it is flat (0) except for
% ZDW randomization. The ZDW detuning is ~N(0, 2 nm) around 1310 nm and is
% converted to dispersion via 0.09 ps/(nm·km) per nm.
% constants (matching the Python code)
meanLambda_nm = 1310; % center wavelength
sigma_nm = 2; % ZDW sigma
Dslope = 0.07; % ps/(nm·km) per nm detuning
% random ZDW-induced dispersion offset
if randomize_ZDW
rng(randomkey, 'twister');
rand_zdws_nm = meanLambda_nm + sigma_nm .* randn(N,1);
rand_D = (rand_zdws_nm - ref_zdw) .* Dslope; % ps/(nm·km)
else
rand_D = zeros(N,1);
end
% nominal segmented pattern (match Python intent; keep length N)
if D > 0
base = (-1) .^ ((0:N-1).'); % +1,-1,+1,-1,...
else % D == 0 (or anything else)
base = ones(N,1);
end
dispersion_vector = base .* D + rand_D; % ps/(nm·km)
end

198
projects/WDM/plot_BER_vs_ROP.m Normal file
View File

@@ -0,0 +1,198 @@
%% =========================
% Routine A: BER vs ROP plots
%% =========================
function S = plot_BER_vs_ROP(res, varargin)
% S = plot_BER_vs_ROP(res, 'fec', 3.8e-3, 'eval_list', [], 'colormap', 'Spectral')
%
% Creates one BER-vs-ROP figure per evaluated distance.
% Returns S struct with FEC crossings: S.Sffe, S.Sdfe, S.Svnle, S.Smlse, S.Sdbt
%
% res.* assumed: res.ffe{ch,rop,realiz,eval_distance} etc.
p = inputParser;
p.addParameter('fec', 3.8e-3, @(x)isnumeric(x)&&isscalar(x));
p.addParameter('eval_list', [], @(x)isnumeric(x));
p.addParameter('colormap', 'Spectral', @(x)ischar(x) || isstring(x));
p.parse(varargin{:});
fec = p.Results.fec;
% -------------------- Basic metadata --------------------
rop = res.settings.rop(:);
wavelengthplan = res.settings.wavelengthplan(:);
distances = res.eval_dist_km(:);
N_ch = numel(wavelengthplan);
N_rop = numel(rop);
% -------------------- Determine dims robustly --------------------
dims = size(res.ffe);
if numel(dims) < 4, dims(end+1:4) = 1; end
N_distances = dims(4);
if numel(distances) ~= N_distances
distances = (1:N_distances).';
end
% -------------------- Choose eval distances --------------------
eval_list = p.Results.eval_list;
if isempty(eval_list)
eval_list = 1:N_distances;
end
% -------------------- Colors --------------------
try
cols = cbrewer2(char(p.Results.colormap), N_ch);
catch
cols = linspecer(N_ch);
end
% -------------------- Crossings containers --------------------
S = struct();
S.rop = rop;
S.wavelengthplan = wavelengthplan;
S.distances = distances;
S.fec = fec;
S.Sffe = cell(N_distances, N_ch);
S.Sdfe = cell(N_distances, N_ch);
S.Svnle = cell(N_distances, N_ch);
S.Smlse = cell(N_distances, N_ch);
S.Sdbt = cell(N_distances, N_ch);
% -------------------- Plot per distance --------------------
for eval_ptr = eval_list
figure('Name', sprintf('BER vs ROP @ %s', distLabel(distances, eval_ptr)));
hold on;
for ch = 1:N_ch
% Extract cell slices
ffe_cells = sliceCells4D(res.ffe , ch, eval_ptr, N_rop);
dfe_cells = sliceCells4D(res.dfe , ch, eval_ptr, N_rop);
vnle_cells = sliceCells4D(res.vnle, ch, eval_ptr, N_rop);
mlse_cells = sliceCells4D(res.mlse, ch, eval_ptr, N_rop);
dbt_cells = sliceCells4D(res.dbt , ch, eval_ptr, N_rop);
% Crossings (only depends on cells)
[S.Sffe{eval_ptr,ch}, ~] = fecCrossings(rop, ffe_cells, fec);
[S.Sdfe{eval_ptr,ch}, ~] = fecCrossings(rop, dfe_cells, fec);
[S.Svnle{eval_ptr,ch}, ~] = fecCrossings(rop, vnle_cells, fec);
[S.Smlse{eval_ptr,ch}, ~] = fecCrossings(rop, mlse_cells, fec);
[S.Sdbt{eval_ptr,ch}, ~] = fecCrossings(rop, dbt_cells, fec);
% BER matrices (complete realizations only)
ffe_mat = extractCompleteBER(ffe_cells); % N_rop x K
if ~isempty(ffe_mat)
plot(rop, ffe_mat, 'Color', cols(ch,:), ...
'DisplayName', sprintf('FFE @ %dnm', round(wavelengthplan(ch))));
end
end
set(gca,'XScale','linear','YScale','log', ...
'TickLabelInterpreter','latex','FontSize',11);
yline([3.8e-3, 2.2e-4], 'HandleVisibility','off','LineWidth',1.5);
ylabel('BER');
xlabel('ROP [dB]');
title(sprintf('BER vs. ROP @ %s', distLabel(distances, eval_ptr)));
xlim([min(rop) max(rop)]);
ylim([1e-5 0.3]);
grid on;
legend show;
beautifyBERplot;
end
end
%% =========================
% Shared helpers (unchanged)
%% =========================
function cells2D = sliceCells4D(C4, ch, eval_ptr, N_rop)
dims = size(C4);
if numel(dims) < 4, dims(end+1:4) = 1; end
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
return;
end
ropUse = min(N_rop, ropMax);
rUse = min(N_realiz, rMax);
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);
end
function [S, noCrossingMask] = fecCrossings(rop, cellsNxR, fec)
Y = extractCompleteBER(cellsNxR);
if isempty(Y)
S = [];
noCrossingMask = [];
return;
end
ok = mean(Y,1,'omitnan') <= 0.1;
Y = Y(:, ok);
if isempty(Y)
S = [];
noCrossingMask = [];
return;
end
nR = size(Y,2);
S = nan(1,nR);
noCrossingMask = true(1,nR);
rop = rop(:);
for j = 1:nR
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);
S(j) = x1 + t*(x2 - x1);
noCrossingMask(j) = false;
end
end
end
end
function Y = extractCompleteBER(cellSlice)
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

219
projects/WDM/plot_FEC_violin.m Normal file
View File

@@ -0,0 +1,219 @@
function plot_FEC_violin(res, varargin)
% 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(res, 'tech','VNLE', 'fec',3.8e-3, 'eval_ptr',[], 'ylim',[-10 -6], 'bandwidth',0.15);
% -------- args --------
p = inputParser;
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;
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
rop = res.settings.rop(:);
wavelengthplan = res.settings.wavelengthplan(:);
distances = res.eval_dist_km(:);
N_ch = numel(wavelengthplan);
N_rop = numel(rop);
dims = size(C4);
if numel(dims) < 4, dims(end+1:4) = 1; end
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
% -------- 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
% -------- plot --------
figure('Name', sprintf('%s FEC violin @ %s', tech, distLabel(distances, eval_ptr)));
hold on;
% 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');
grid on; box on;
if ~isempty(opt.ylim)
ylim(opt.ylim);
end
end
% ================= helpers =================
function cells2D = sliceCells4D(C4, ch, eval_ptr, N_rop)
dims = size(C4);
if numel(dims) < 4, dims(end+1:4) = 1; end
N_realiz = dims(3);
cells2D = cell(N_rop, N_realiz);
if ch > dims(1) || eval_ptr > dims(4)
return;
end
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 [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
% optional quality gate (keep your style)
ok = mean(Y,1,'omitnan') <= 0.1;
Y = Y(:,ok);
K_total = size(Y,2);
if K_total == 0
S = [];
K_cross = 0;
return;
end
rop = rop(:);
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);
S(j) = x1 + t*(x2 - x1);
end
end
end
K_cross = sum(isfinite(S));
S = S(isfinite(S));
end
function Y = extractCompleteBER(cellSlice)
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 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