From 45c9fc2b1433e4c0925c4763dfbb3249eb0f5445 Mon Sep 17 00:00:00 2001 From: "silas (home)" Date: Tue, 23 Dec 2025 12:35:49 +0100 Subject: [PATCH] noew WITH the new code :_) --- projects/WDM/WDM_model_10km.m | 361 ++++++++++++++++++++++++++++++++++ 1 file changed, 361 insertions(+) create mode 100644 projects/WDM/WDM_model_10km.m diff --git a/projects/WDM/WDM_model_10km.m b/projects/WDM/WDM_model_10km.m new file mode 100644 index 0000000..85cf2d7 --- /dev/null +++ b/projects/WDM/WDM_model_10km.m @@ -0,0 +1,361 @@ +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 + +% 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); + +%% 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.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 +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 = 5e12; % 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 + +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.75; + +%% ---------- 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 = [2 4 6 8 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_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 ---------- + parfor 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).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 + + 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); + + % NOTE: For correctness and to keep changes minimal, we keep the exact DSP chain. + for ri = 1:length(s.rop) + + 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}); + + %%%%%% 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 + 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; + + % -------------------- 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.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