diff --git a/projects/WDM/WDM_auswertung.m b/projects/WDM/WDM_auswertung.m index 508291d..60651c1 100644 --- a/projects/WDM/WDM_auswertung.m +++ b/projects/WDM/WDM_auswertung.m @@ -1,68 +1,68 @@ +%saved like this: output_ffe{ch,rop,realiz,eval_distance} -try - rop = res.settings.rop; % 12 points - wavelengthplan = res.settings.wavelengthplan; -catch - wavelengthplan = [1295,1305,1315,1325]; - wavelengthplan = calcWavelengthPlan(16,400e9,1310); - rop = -8.25:0.75:0; -end - -N = length(wavelengthplan); +dims = size(res.ffe); +rop = res.settings.rop; +wavelengthplan = res.settings.wavelengthplan; +N_distances = dims(4); %usually 2km and 10km are evaualted +distances = res.eval_dist_km; %[2,10] +N_ch = length(wavelengthplan); %16 channels 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 = 3.8e-3; -Sffe = cell(1,N); -Svnle = cell(1,N); -Smlse = cell(1,N); -Sdbt = cell(1,N); +Sffe = cell(1,N_ch); +Svnle = cell(1,N_ch); +Smlse = cell(1,N_ch); +Sdbt = cell(1,N_ch); % 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) qHigh = 1; % upper quantile (e.g., 0.99 for old script) -cols = linspecer(N); % one color per wavelength (Ch) -cols = cbrewer2('set1',N); +cols = linspecer(N_ch); % one color per wavelength (Ch) +cols = cbrewer2('RdBu',N_ch); -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),[]); +for eval_ptr = 2 + for ch = 1:N_ch - [Sffe{l}, noX_ffe] = fecCrossings(rop, ffe_cells, fec); + % Slice 12x50 cell arrays + ffe_cells = reshape(squeeze(res.ffe(ch,:,:,eval_ptr)),length(rop),[]); + vnle_cells = reshape(squeeze(res.vnle(ch,:,:,eval_ptr)),length(rop),[]); + mlse_cells = reshape(squeeze(res.mlse(ch,:,:,eval_ptr)),length(rop),[]); + dbt_cells = reshape(squeeze(res.dbt(ch,:,:,eval_ptr)),length(rop),[]); - [Svnle{l}, noX_ffe] = fecCrossings(rop, vnle_cells, fec); + [Sffe{ch}, noX_ffe] = fecCrossings(rop, ffe_cells, fec); - [Smlse{l}, noX_ffe] = fecCrossings(rop, mlse_cells, fec); + [Svnle{ch}, noX_ffe] = fecCrossings(rop, vnle_cells, fec); - [Sdbt{l}, noX_ffe] = fecCrossings(rop, dbt_cells, fec); + [Smlse{ch}, noX_ffe] = fecCrossings(rop, mlse_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 + [Sdbt{ch}, noX_ffe] = fecCrossings(rop, dbt_cells, fec); - showLegend = 1; % one legend entry per technique + % 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 - % 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'); + showLegend = 1; % one legend entry per technique - % plotBandMeanBL(rop, vnle_mat, cols(l,:), sprintf('VNLE @ %.1d nm',round(wavelengthplan(l))), qLow, qHigh, '--x', showLegend); + % Plot shaded band + mean line with boundedline + % 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'); - plotBandMeanBL(rop, mlse_mat, cols(l,:), sprintf('VNLE+PF+MLSE @ %.1d nm',round(wavelengthplan(l))), qLow, qHigh, '-o', showLegend); + % plotBandMeanBL(rop, vnle_mat, cols(ch,:), sprintf('VNLE @ %.1d nm',round(wavelengthplan(l))), qLow, qHigh, '--x', showLegend); - % plotBandMeanBL(rop, dbt_mat, cols(l,:), sprintf('DBt.+MLSE @ %.1d nm',round(wavelengthplan(l))), qLow, qHigh, '--v', showLegend); + % plotBandMeanBL(rop, mlse_mat, cols(ch,:), sprintf('VNLE+PF+MLSE @ %.1d nm',round(wavelengthplan(ch))), qLow, qHigh, '-o', showLegend); - set(gca,'XScale','linear','YScale','log','TickLabelInterpreter','latex','FontSize',11); - yline([3.8e-3, 2.2e-4], 'HandleVisibility','off','LineWidth',1.5); + % 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); + yline([3.8e-3, 2.2e-4], 'HandleVisibility','off','LineWidth',1.5); + + end end ylabel('BER'); @@ -76,23 +76,23 @@ legend show; %% VIOLIN S_cell = Sdbt; -S_cell =Smlse; +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}); + S_mat = NaN(Kmax, N_ch); + for ch = 1:N_ch + k = numel(S_cell{i}{ch}); if k > 0 - S_mat(1:k, l) = S_cell{i}{l}; + S_mat(1:k, ch) = S_cell{i}{ch}; 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, ... @@ -105,9 +105,9 @@ for i = 1:length(S_cell) 'ShowData', true, ... 'ShowBox', false, ... 'Bandwidth', 0.05); - + ylim([floor(min(S_mat,[],'all')), ceil(max(S_mat,[],'all'))]) - % ylim([-8 0]); + % ylim([-8 0]); ylabel('ROP at FEC crossing'); title(sprintf('RROP to cross BER %.2e', fec)); grid on; box on; @@ -115,51 +115,43 @@ for i = 1:length(S_cell) end - - - - - - - - %% ================= helper ================= function plotBandMeanBL(x, Y, color, techLabel, qLow, qHigh, lineSpec, showLegend) - % Y: (nPoints x nRealizations) - % Remove realizations that are entirely zero (like removeZeros behavior) - badCols = all(Y == 0, 1); - Y(:, badCols) = []; - - Y(Y==0) = 1e-8; - % Stats across realizations - mu = mean(Y, 2, 'omitnan'); % mean line - lo = quantile(Y, qLow, 2); % lower bound - hi = quantile(Y, qHigh, 2); % upper bound +% Y: (nPoints x nRealizations) +% Remove realizations that are entirely zero (like removeZeros behavior) +badCols = all(Y == 0, 1); +Y(:, badCols) = []; - % Convert to asymmetric distances required by boundedline: - % b(:,1) = distance to lower side; b(:,2) = distance to upper side - b = [mu - lo, hi - mu]; +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 - % 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 +% Convert to asymmetric distances required by boundedline: +% b(:,1) = distance to lower side; b(:,2) = distance to upper side +b = [mu - lo, hi - mu]; - % Single legend entry per technique (use first wavelength only) - if showLegend - set(hl, 'DisplayName', techLabel); - else - set(hl, 'HandleVisibility','off'); - end +% 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 - % 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 +% 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) @@ -173,52 +165,52 @@ function [S, noCrossingMask, Y_keep] = fecCrossings(rop, cells12xR, fec) % 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 +% 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 +% 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); +% 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); +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'); + % 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 ~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 + 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; +Y_keep = Y; end @@ -226,28 +218,28 @@ 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); +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); +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 diff --git a/projects/WDM/WDM_model_10km.m b/projects/WDM/WDM_model_10km.m index 85cf2d7..d8977eb 100644 --- a/projects/WDM/WDM_model_10km.m +++ b/projects/WDM/WDM_model_10km.m @@ -33,6 +33,7 @@ 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; @@ -44,6 +45,7 @@ if isempty(p) || p.NumWorkers ~= cpus 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'); @@ -106,7 +108,7 @@ Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",1 N = numel(s.wavelengthplan); 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_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. 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); % 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) 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); @@ -181,7 +184,7 @@ for realiz = 1:s.num_realiz Tx_bits = cell(1,N); %% ---------- TX per channel ---------- - parfor l = 1:N + for l = 1:N [Digi_sig,Symbols{l},Tx_bits{l}] = PAMsource( ... "fsym",fsym,"M",s.M,"order",18,"useprbs",0, ... @@ -235,6 +238,8 @@ for realiz = 1:s.num_realiz 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, ... @@ -247,10 +252,9 @@ for realiz = 1:s.num_realiz 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); - % 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 l = 1:N %%%%%% ROP %%%%%% 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; + % -------------------- 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]; @@ -348,6 +365,7 @@ for realiz = 1:s.num_realiz 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; diff --git a/projects/WDM/WDM_model_10km_queue.m b/projects/WDM/WDM_model_10km_queue.m new file mode 100644 index 0000000..7e1832d --- /dev/null +++ b/projects/WDM/WDM_model_10km_queue.m @@ -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