%%% Run parameters % TX 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, 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, 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 = options.fiber_length_km; s.pmd = 0.1; s.gamma = 0.0023; s.M = 4; m = floor(log2(s.M)*10)/10; 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_dc = 0; 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); 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; 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,... "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(); % 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); % 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",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",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 = 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); 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); 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) 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) 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_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",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_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 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} = ffe_results; %VNLE 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} = vnle_results; output_mlse{l,ri,realiz} = mlse_results; % DB tgt. 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} = dbt_results; end end res = struct(); res.settings = s; res.ffe = output_ffe; res.vnle = output_vnle; res.mlse = output_mlse; res.dbt = output_dbt; % Save results 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