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

new par function

Browse Source
This commit is contained in:
silas (home)
2025-12-31 12:46:56 +01:00
parent 45c9fc2b14
commit 2354d7ac44
3 changed files with 590 additions and 151 deletions
Download Patch File Download Diff File Expand all files Collapse all files

188
projects/WDM/WDM_auswertung.m
View File

@@ -1,45 +1,44 @@
%saved like this: output_ffe{ch,rop,realiz,eval_distance}
try dims = size(res.ffe);
rop = res.settings.rop; % 12 points rop = res.settings.rop;
wavelengthplan = res.settings.wavelengthplan; wavelengthplan = res.settings.wavelengthplan;
catch N_distances = dims(4); %usually 2km and 10km are evaualted
wavelengthplan = [1295,1305,1315,1325]; distances = res.eval_dist_km; %[2,10]
wavelengthplan = calcWavelengthPlan(16,400e9,1310); N_ch = length(wavelengthplan); %16 channels
rop = -8.25:0.75:0;
end
N = length(wavelengthplan);
figure(); hold on; figure(); hold on;
cols = cbrewer2('set2',N); % one color per wavelength (Ch) cols = cbrewer2('set2',N_ch); % one color per wavelength (Ch)
fec = 2.2e-4; fec = 2.2e-4;
fec = 3.8e-3; fec = 3.8e-3;
Sffe = cell(1,N); Sffe = cell(1,N_ch);
Svnle = cell(1,N); Svnle = cell(1,N_ch);
Smlse = cell(1,N); Smlse = cell(1,N_ch);
Sdbt = cell(1,N); Sdbt = cell(1,N_ch);
% Choose your quantile band. For your old style, use 0.04/0.99: % Choose your quantile band. For your old style, use 0.04/0.99:
qLow = 0.0; % lower quantile (e.g., 0.04 for old script) qLow = 0.0; % lower quantile (e.g., 0.04 for old script)
qHigh = 1; % upper quantile (e.g., 0.99 for old script) qHigh = 1; % upper quantile (e.g., 0.99 for old script)
cols = linspecer(N); % one color per wavelength (Ch) cols = linspecer(N_ch); % one color per wavelength (Ch)
cols = cbrewer2('set1',N); cols = cbrewer2('RdBu',N_ch);
for eval_ptr = 2
for ch = 1:N_ch
for l = 1:N
% Slice 12x50 cell arrays % Slice 12x50 cell arrays
ffe_cells = reshape(squeeze(res.ffe(l,:,:)),length(rop),[]); ffe_cells = reshape(squeeze(res.ffe(ch,:,:,eval_ptr)),length(rop),[]);
vnle_cells = reshape(squeeze(res.vnle(l,:,:)),length(rop),[]); vnle_cells = reshape(squeeze(res.vnle(ch,:,:,eval_ptr)),length(rop),[]);
mlse_cells = reshape(squeeze(res.mlse(l,:,:)),length(rop),[]); mlse_cells = reshape(squeeze(res.mlse(ch,:,:,eval_ptr)),length(rop),[]);
dbt_cells = reshape(squeeze(res.dbt(l,:,:)),length(rop),[]); dbt_cells = reshape(squeeze(res.dbt(ch,:,:,eval_ptr)),length(rop),[]);
[Sffe{l}, noX_ffe] = fecCrossings(rop, ffe_cells, fec); [Sffe{ch}, noX_ffe] = fecCrossings(rop, ffe_cells, fec);
[Svnle{l}, noX_ffe] = fecCrossings(rop, vnle_cells, fec); [Svnle{ch}, noX_ffe] = fecCrossings(rop, vnle_cells, fec);
[Smlse{l}, noX_ffe] = fecCrossings(rop, mlse_cells, fec); [Smlse{ch}, noX_ffe] = fecCrossings(rop, mlse_cells, fec);
[Sdbt{l}, noX_ffe] = fecCrossings(rop, dbt_cells, fec); [Sdbt{ch}, noX_ffe] = fecCrossings(rop, dbt_cells, fec);
% Extract BER matrices using only complete realizations (12/12 ROP filled) % Extract BER matrices using only complete realizations (12/12 ROP filled)
ffe_mat = extractCompleteBER(ffe_cells); % 12 x K_ffe ffe_mat = extractCompleteBER(ffe_cells); % 12 x K_ffe
@@ -51,18 +50,19 @@ for l = 1:N
showLegend = 1; % one legend entry per technique showLegend = 1; % one legend entry per technique
% Plot shaded band + mean line with boundedline % Plot shaded band + mean line with boundedline
% plotBandMeanBL(rop, ffe_mat, cols(l,:), sprintf('FFE @ %.1d nm',round(wavelengthplan(l))), qLow, qHigh, '--s', showLegend); % plotBandMeanBL(rop, ffe_mat, cols(ch,:), sprintf('FFE @ %.1d nm',round(wavelengthplan(l))), qLow, qHigh, '--s', showLegend);
% scatter(Sffe,fec.*ones(size(Sffe)),20,'v','MarkerFaceColor','black'); % 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, vnle_mat, cols(ch,:), 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, mlse_mat, cols(ch,:), sprintf('VNLE+PF+MLSE @ %.1d nm',round(wavelengthplan(ch))), qLow, qHigh, '-o', showLegend);
% plotBandMeanBL(rop, dbt_mat, cols(l,:), sprintf('DBt.+MLSE @ %.1d nm',round(wavelengthplan(l))), qLow, qHigh, '--v', showLegend); % plotBandMeanBL(rop, dbt_mat, cols(ch,:), sprintf('DBt.+MLSE @ %.1d nm',round(wavelengthplan(l))), qLow, qHigh, '--v', showLegend);
set(gca,'XScale','linear','YScale','log','TickLabelInterpreter','latex','FontSize',11); set(gca,'XScale','linear','YScale','log','TickLabelInterpreter','latex','FontSize',11);
yline([3.8e-3, 2.2e-4], 'HandleVisibility','off','LineWidth',1.5); yline([3.8e-3, 2.2e-4], 'HandleVisibility','off','LineWidth',1.5);
end
end end
ylabel('BER'); ylabel('BER');
@@ -76,18 +76,18 @@ legend show;
%% VIOLIN %% VIOLIN
S_cell = Sdbt; S_cell = Sdbt;
S_cell =Smlse; S_cell = Smlse;
S_cell = {Svnle,Smlse,Sdbt}; S_cell = {Svnle,Smlse,Sdbt};
S_cell = {Svnle}; S_cell = {Svnle};
figure(5); hold on; figure(5); hold on;
for i = 1:length(S_cell) for i = 1:length(S_cell)
% Pad to rectangular matrix: rows = realizations, cols = wavelengths % Pad to rectangular matrix: rows = realizations, cols = wavelengths
Kmax = max(cellfun(@numel, S_cell{i})); Kmax = max(cellfun(@numel, S_cell{i}));
S_mat = NaN(Kmax, N); S_mat = NaN(Kmax, N_ch);
for l = 1:N for ch = 1:N_ch
k = numel(S_cell{i}{l}); k = numel(S_cell{i}{ch});
if k > 0 if k > 0
S_mat(1:k, l) = S_cell{i}{l}; S_mat(1:k, ch) = S_cell{i}{ch};
end end
end end
@@ -115,51 +115,43 @@ for i = 1:length(S_cell)
end end
%% ================= helper ================= %% ================= helper =================
function plotBandMeanBL(x, Y, color, techLabel, qLow, qHigh, lineSpec, showLegend) function plotBandMeanBL(x, Y, color, techLabel, qLow, qHigh, lineSpec, showLegend)
% Y: (nPoints x nRealizations) % Y: (nPoints x nRealizations)
% Remove realizations that are entirely zero (like removeZeros behavior) % Remove realizations that are entirely zero (like removeZeros behavior)
badCols = all(Y == 0, 1); badCols = all(Y == 0, 1);
Y(:, badCols) = []; Y(:, badCols) = [];
Y(Y==0) = 1e-8; Y(Y==0) = 1e-8;
% Stats across realizations % Stats across realizations
mu = mean(Y, 2, 'omitnan'); % mean line mu = mean(Y, 2, 'omitnan'); % mean line
lo = quantile(Y, qLow, 2); % lower bound lo = quantile(Y, qLow, 2); % lower bound
hi = quantile(Y, qHigh, 2); % upper bound hi = quantile(Y, qHigh, 2); % upper bound
% Convert to asymmetric distances required by boundedline: % Convert to asymmetric distances required by boundedline:
% b(:,1) = distance to lower side; b(:,2) = distance to upper side % b(:,1) = distance to lower side; b(:,2) = distance to upper side
b = [mu - lo, hi - mu]; b = [mu - lo, hi - mu];
% Call boundedline with alpha shading % Call boundedline with alpha shading
[hl, hp] = boundedline(x(:), mu(:), b, lineSpec, 'alpha', ... [hl, hp] = boundedline(x(:), mu(:), b, lineSpec, 'alpha', ...
'transparency', 0.18); 'transparency', 0.18);
% Color styling % Color styling
set(hl, 'Color', color, 'LineWidth', 1.4, 'MarkerSize', 4); set(hl, 'Color', color, 'LineWidth', 1.4, 'MarkerSize', 4);
set(hp, 'FaceColor', color, 'HandleVisibility','off'); % patch hidden in legend set(hp, 'FaceColor', color, 'HandleVisibility','off'); % patch hidden in legend
% Single legend entry per technique (use first wavelength only) % Single legend entry per technique (use first wavelength only)
if showLegend if showLegend
set(hl, 'DisplayName', techLabel); set(hl, 'DisplayName', techLabel);
else else
set(hl, 'HandleVisibility','off'); set(hl, 'HandleVisibility','off');
end end
% Optional: outline the bounds if outlinebounds is available % Optional: outline the bounds if outlinebounds is available
if exist('outlinebounds','file') == 2 if exist('outlinebounds','file') == 2
ho = outlinebounds(hl, hp); ho = outlinebounds(hl, hp);
set(ho, 'linestyle', ':', 'color', color, 'linewidth', 1, ... set(ho, 'linestyle', ':', 'color', color, 'linewidth', 1, ...
'HandleVisibility','off'); 'HandleVisibility','off');
end end
end end
function [S, noCrossingMask, Y_keep] = fecCrossings(rop, cells12xR, fec) function [S, noCrossingMask, Y_keep] = fecCrossings(rop, cells12xR, fec)
@@ -173,32 +165,32 @@ function [S, noCrossingMask, Y_keep] = fecCrossings(rop, cells12xR, fec)
% noCrossingMask 1xK logical mask: true if no crossing for that realization % noCrossingMask 1xK logical mask: true if no crossing for that realization
% Y_keep 12xK numeric BER matrix used for the crossing detection % Y_keep 12xK numeric BER matrix used for the crossing detection
% 1) keep only complete realization columns % 1) keep only complete realization columns
Y = extractCompleteBER(cells12xR); % -> 12 x K Y = extractCompleteBER(cells12xR); % -> 12 x K
if isempty(Y) if isempty(Y)
S = []; S = [];
noCrossingMask = []; noCrossingMask = [];
Y_keep = Y; Y_keep = Y;
return; return;
end end
% 2) optionally drop realizations with mean BER > 0.1 % 2) optionally drop realizations with mean BER > 0.1
ok = mean(Y,1,'omitnan') <= 0.1; ok = mean(Y,1,'omitnan') <= 0.1;
Y = Y(:, ok); Y = Y(:, ok);
if isempty(Y) if isempty(Y)
S = []; S = [];
noCrossingMask = []; noCrossingMask = [];
Y_keep = Y; Y_keep = Y;
return; return;
end end
% 3) find crossings per realization % 3) find crossings per realization
nR = size(Y,2); nR = size(Y,2);
S = nan(1,nR); S = nan(1,nR);
noCrossingMask = true(1,nR); noCrossingMask = true(1,nR);
rop = rop(:); % ensure column rop = rop(:); % ensure column
for j = 1:nR for j = 1:nR
y = Y(:,j); y = Y(:,j);
% sign change from >fec to <=fec (first time it drops below FEC) % sign change from >fec to <=fec (first time it drops below FEC)
@@ -216,9 +208,9 @@ function [S, noCrossingMask, Y_keep] = fecCrossings(rop, cells12xR, fec)
noCrossingMask(j) = false; noCrossingMask(j) = false;
end end
end end
end end
Y_keep = Y; Y_keep = Y;
end end
@@ -226,28 +218,28 @@ end
function Y = extractCompleteBER(cellSlice) function Y = extractCompleteBER(cellSlice)
% cellSlice: 12xR cell array; each cell should be a struct with .metrics.BER % 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. % Keep only those realization columns where ALL 12 ROP entries are valid.
if isempty(cellSlice), Y = []; return; end if isempty(cellSlice), Y = []; return; end
nR = size(cellSlice,2); nR = size(cellSlice,2);
keep = false(1,nR); keep = false(1,nR);
for r = 1:nR for r = 1:nR
col = cellSlice(:,r); col = cellSlice(:,r);
keep(r) = all(cellfun(@(c) ~isempty(c) , col)); keep(r) = all(cellfun(@(c) ~isempty(c) , col));
end end
if ~any(keep), Y = []; return; end if ~any(keep), Y = []; return; end
Y = cellfun(@(c) c.metrics.BER, cellSlice(:,keep), 'UniformOutput', true); Y = cellfun(@(c) c.metrics.BER, cellSlice(:,keep), 'UniformOutput', true);
end end
function Y = extractCompleteAlphas(cellSlice) function Y = extractCompleteAlphas(cellSlice)
% cellSlice: 12xR cell array; each cell should be a struct with .metrics.BER % 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. % Keep only those realization columns where ALL 12 ROP entries are valid.
if isempty(cellSlice), Y = []; return; end if isempty(cellSlice), Y = []; return; end
nR = size(cellSlice,2); nR = size(cellSlice,2);
keep = false(1,nR); keep = false(1,nR);
for r = 1:nR for r = 1:nR
col = cellSlice(:,r); col = cellSlice(:,r);
keep(r) = all(cellfun(@(c) ~isempty(c) , col)); keep(r) = all(cellfun(@(c) ~isempty(c) , col));
end end
if ~any(keep), Y = []; return; end if ~any(keep), Y = []; return; end
Y = cellfun(@(c) c.metrics.Alpha, cellSlice(:,keep), 'UniformOutput', true); Y = cellfun(@(c) c.metrics.Alpha, cellSlice(:,keep), 'UniformOutput', true);
end end

30
projects/WDM/WDM_model_10km.m
View File

@@ -33,6 +33,7 @@ jsl = fullfile(tmpbase, sprintf('matlab_jobstorage_%s_%s', ...
getenv('USER'), getenv('SLURM_JOB_ID'))); getenv('USER'), getenv('SLURM_JOB_ID')));
if ~exist(jsl,'dir'); mkdir(jsl); end if ~exist(jsl,'dir'); mkdir(jsl); end
if 0
% Configure the local cluster explicitly and start the pool % Configure the local cluster explicitly and start the pool
c = parcluster('local'); c = parcluster('local');
c.NumWorkers = cpus; c.NumWorkers = cpus;
@@ -44,6 +45,7 @@ if isempty(p) || p.NumWorkers ~= cpus
p = parpool(c, cpus); p = parpool(c, cpus);
end end
fprintf('parpool up with %d workers; JobStorage=%s\n', p.NumWorkers, c.JobStorageLocation); fprintf('parpool up with %d workers; JobStorage=%s\n', p.NumWorkers, c.JobStorageLocation);
end
%% result filename (timestamp + optional job id) %% result filename (timestamp + optional job id)
t = datetime('now','TimeZone','local','Format','yyyyMMdd_HHmmss'); t = datetime('now','TimeZone','local','Format','yyyyMMdd_HHmmss');
@@ -106,7 +108,7 @@ Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",1
N = numel(s.wavelengthplan); N = numel(s.wavelengthplan);
f_plan = physconst('lightspeed')./(s.wavelengthplan.*1e-9); f_plan = physconst('lightspeed')./(s.wavelengthplan.*1e-9);
margin = 5e12; % some THz left and right margin = 20e12; % some THz left and right
f_span = (max(f_plan)+margin)-(min(f_plan)-margin); f_span = (max(f_plan)+margin)-(min(f_plan)-margin);
f_nyq = f_span/2; f_nyq = f_span/2;
@@ -128,7 +130,7 @@ s.rop = -12:0.75:-0.75;
% Evaluate BER at these intermediate distances (km), plus always include the final link_length if > 0. % 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) segment_length = 1; % km (must match fiber loop below)
eval_dist_km = [2 4 6 8 10]; eval_dist_km = [2 10];
eval_dist_km = eval_dist_km(eval_dist_km <= link_length); eval_dist_km = eval_dist_km(eval_dist_km <= link_length);
% Always include final distance (if > 0) and avoid duplicates % Always include final distance (if > 0) and avoid duplicates
@@ -149,6 +151,7 @@ nEval = numel(eval_seg);
%% Preallocate outputs (add eval distance as 4th dimension) %% Preallocate outputs (add eval distance as 4th dimension)
output_ffe = cell(length(s.wavelengthplan), length(s.rop), s.num_realiz, nEval); 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_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_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); output_dbt = cell(length(s.wavelengthplan), length(s.rop), s.num_realiz, nEval);
@@ -181,7 +184,7 @@ for realiz = 1:s.num_realiz
Tx_bits = cell(1,N); Tx_bits = cell(1,N);
%% ---------- TX per channel ---------- %% ---------- TX per channel ----------
parfor l = 1:N for l = 1:N
[Digi_sig,Symbols{l},Tx_bits{l}] = PAMsource( ... [Digi_sig,Symbols{l},Tx_bits{l}] = PAMsource( ...
"fsym",fsym,"M",s.M,"order",18,"useprbs",0, ... "fsym",fsym,"M",s.M,"order",18,"useprbs",0, ...
@@ -235,6 +238,8 @@ for realiz = 1:s.num_realiz
for seg = 1:nSegments 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, ... 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, ... "beat_len",10,"corr_len",100,"dz",1,"manakov",0, ...
"gamma",s.gamma,"lambda",zdw,"n_waveplates",10,"SS_dphimax",0.01, ... "gamma",s.gamma,"lambda",zdw,"n_waveplates",10,"SS_dphimax",0.01, ...
@@ -247,10 +252,9 @@ for realiz = 1:s.num_realiz
Opt_sig_wdm_demux = Optical_Demultiplex("attenuation",0,"B",200e9,"filtype",1, ... 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); "fs_out",fdac*kover,"fs_in",fdac*kover*upsample_pow,"lambda_center",1310).process(Opt_sig_wdm_fib);
% NOTE: For correctness and to keep changes minimal, we keep the exact DSP chain.
for ri = 1:length(s.rop)
parfor l = 1:N for ri = 1:length(s.rop)
for l = 1:N
%%%%%% ROP %%%%%% %%%%%% ROP %%%%%%
Opt_sig_wdm_rx = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power", ... Opt_sig_wdm_rx = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power", ...
@@ -293,6 +297,19 @@ for realiz = 1:s.num_realiz
output_ffe{l,ri,realiz,eval_ptr} = ffe_results; 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 -------------------- % -------------------- VNLE + MLSE --------------------
pf_ncoeffs = 1; pf_ncoeffs = 1;
ffe_order = [50, 5, 5]; ffe_order = [50, 5, 5];
@@ -348,6 +365,7 @@ for realiz = 1:s.num_realiz
res.settings = s; res.settings = s;
res.eval_dist_km = eval_dist_km; res.eval_dist_km = eval_dist_km;
res.ffe = output_ffe; res.ffe = output_ffe;
res.dfe = output_dfe;
res.vnle = output_vnle; res.vnle = output_vnle;
res.mlse = output_mlse; res.mlse = output_mlse;
res.dbt = output_dbt; res.dbt = output_dbt;

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

@@ -0,0 +1,429 @@
function WDM_model_10km_queue(options)
%WDM_model_10km_queue
% Job-queue / pipeline version of your script:
% - Main thread runs fiber propagation segment-by-segment
% - At each eval distance: demux once and enqueue per-(l,ri) DSP jobs via parfeval
% - While continuing propagation: periodically harvest completed jobs (non-blocking)
% - End of realization: drain remaining jobs, then save
%
% Save this as a NEW file, e.g. "WDM_model_10km_queue.m"
%%% Run parameters
arguments
options.num_channels = 16;
options.channel_spacing = 400e9;
options.fiber_length_km = 2;
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
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')); 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);
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.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.1;
s.gamma = 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 = 25e12; % 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 = -12:0.75:-0.75;
%% ---------- Intermediate evaluation points (distance dimension) ----------
segment_length = 1; % km (must match fiber loop below)
eval_dist_km = [2 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
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",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);
El_sig = El_sig.normalize("mode","oneone");
scaling = 0.6*(u_pi/2-abs(vbias-u_pi/2));
El_sig = El_sig .* scaling;
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",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;
for seg = 1:nSegments
fprintf('Realiz %d/%d: Segment %d/%d \n', realiz, s.num_realiz, seg, nSegments);
if 1
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);
end
% -------- 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);
for ri = 1:length(s.rop)
for l = 1:N
% Enqueue one job per (l,ri)
F(end+1,1) = parfeval(p, @rx_job, 5, ...
Opt_sig_wdm_demux{l}, ... % channel signal at this eval point
s.rop(ri), ... % rop [dB]
Symbols{l}, Tx_bits{l}, ...% reference + bits
s, l, realiz, ... % settings + indices
fdac, kover, fadc, fsym, ...
len_tr, mu_dc, mu_ffe, mu_dfe, dfe_order, duob_mode);
meta(end+1) = struct('l',l,'ri',ri,'realiz',realiz,'eval_ptr',eval_ptr);
end
end
eval_ptr = eval_ptr + 1;
end
% ----------------------------------------------------------------
% Non-blocking harvest: store finished jobs (if any) while we continue propagation
[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
% Drain remaining jobs for this realization (blocking until queue is empty)
while ~isempty(F)
[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);
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 % 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, realiz, ...
fdac, kover, fadc, fsym, len_tr, mu_dc, mu_ffe, mu_dfe, dfe_order, duob_mode)
% NOTE: keep plotting OFF in workers
debug_plots = 0;
%%%%%% ROP %%%%%%
Opt_sig_wdm_rx = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power", ...
"amplification_db", rop_db).process(Opt_sig_chan);
%%%%%% 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, ~, ~] = Scpe_sig_2sps.tsynch("reference", Symbols_l, "fs_ref", fsym, "debug_plots", debug_plots);
Rx_sig = Scpe_cell{1}.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);
% -------------------- DFE --------------------
dfe_order2 = [50, 0, 0];
eq_dfe = EQ("Ne",dfe_order2,"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);
% -------------------- 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);
% -------------------- DB target --------------------
mlse_db_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",s.M,'trellis_states',PAMmapper(s.M,0).levels);
dbt_results = duobinary_target(eq_v, mlse_db_, s.M, Rx_sig, Symbols_l, Tx_bits_l, ...
"precode_mode", duob_mode, 'showAnalysis', 0, "postFFE", []);
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
% Non-blocking mode: try to fetch as many finished as possible with timeout=0
% Blocking mode: wait for at least one result (timeout=Inf), then keep draining with timeout=0
if block
timeout_first = Inf;
else
timeout_first = 0;
end
% Fetch first (maybe blocking), then drain the rest non-blocking
first = true;
while ~isempty(F)
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
% fetchNext throws if timeout=0 and nothing finished in some releases; treat as "nothing ready"
break;
end
if isempty(k)
break; % nothing ready
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;
% Remove completed future & meta entry (important: indices change afterwards)
F(k) = [];
meta(k) = [];
end
end