%%% Run parameters % TX % --- FIRST LINE: evaluate settings located beside this script --- run(fullfile(fileparts(mfilename('fullpath')),'WDM_settings.m')); M = 4; m = floor(log2(M)*10)/10; fsym = 224e9; fdac = 2*fsym; fadc = 2*fsym; link_length = 0; wavelengthplan = calcWavelengthPlan(16,400e9,1310); wavelengthplan = [1295,1305,1315,1325]; random_key = 2; % Laser / Modulator vbias_rel = 0.5; u_pi = 3.2; 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(wavelengthplan); f_plan = physconst('lightspeed')./(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 = 8; upsample_required = f_nyq./(fdac*kover/2); upsample_pow = 2^nextpow2(upsample_required); upsample_ceil = ceil(upsample_required); f_opt = fdac*kover*upsample_pow; f_opt_nyq = f_opt/2; signal_cell = {}; Symbols = {}; Tx_bits = {}; rop = linspace(-11,0,8); num_realiz = 1; gmi_vnle_bitwise = NaN(length(wavelengthplan),length(rop),num_realiz); snr_vnle= NaN(length(wavelengthplan),length(rop),num_realiz); ber_vnle= NaN(length(wavelengthplan),length(rop),num_realiz); output = cell(length(wavelengthplan),length(rop),num_realiz); for realiz = 1:num_realiz for l = 1:N [Digi_sig,Symbols{l},Tx_bits{l}] = PAMsource(... "fsym",fsym,"M",M,"order",18,"useprbs",0,... "fs_out",fdac,... "applyclipping",0,"clipfactor",1.5,... "applypulseform",apply_pulsef,"pulseformer",Pform,... "randkey",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); % El_sig = 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.setPower(1,"dBm"); % figure;histogram(El_sig.signal); %%%%% MODULATE E/O CONVERSION %%%%% Eml_out = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs,"lambda",wavelengthplan(l),"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",random_key+l+realiz).process(El_sig); signal_cell{l} = Polarization_Controller("mode","rot_power","desired_power",30).process(Eml_out); end Opt_sig_wdm = Optical_Multiplex("fs_in",signal_cell{l}.fs,"fs_out",upsample_pow*Eml_out.fs,... "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); % 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; nSegments = 2; zdw = 1310; D_local = 0; %if ~=0, simulation uses "segmented fiber with d+,d-) randomize_D = true; Dvec = getDispersionVector(nSegments, D_local, zdw, randomize_D, random_key+realiz); for s = 1:nSegments Opt_sig_wdm_fib = DP_Fiber("L",link_length/nSegments,"D",Dvec(s),"Dpmd",0.1,"Ds",0.06,... "beat_len",10,"corr_len",100,"dz",1,"manakov",0,... "gamma",0.0023,"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",link_length/1000,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.07).process(Opt_sig) parfor ri = 1:length(rop) %%%%%% ROP %%%%%% Opt_sig_wdm_rx = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",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 %%%%%% 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",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 %%%%%% 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); %%%%%% 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); 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_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,M,Rx_sig,Symbols{l},Tx_bits{l},... "precode_mode",duob_mode,... 'showAnalysis',0,... "postFFE",[],... "eth_style_symbol_mapping",0); output{l,ri,realiz} = ffe_results; % % FFE or VNLE % eq_ = EQ("Ne",[vnle_order1,vnle_order2,vnle_order3],"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.00,"FFEmu",0,"plotfinal",0,"ideal_dfe",0); % % [eq_signal_sd, eq_noise] = eq_.process(Rx_sig, Symbols{l}); % showEQNoisePSD(eq_noise, "fignum",1273876,"displayname",'noise after EQ'); % [mi_gomez] = calc_air(eq_signal_sd, Symbols{l}, "skip_front", 100, "skip_end", 100); % [gmi_vnle_bitwise(l,ri,realiz)] = calc_ngmi(eq_signal_sd,Symbols{l}); % % [gmi_bitwise_2] = calc_gmi_bitwise(eq_signal_sd,Symbols{l}); % snr_vnle(l,ri,realiz) = calc_snr(Symbols{l}, eq_signal_sd-Symbols{l}); % % % eq_signal_sd.plot("displayname",'bla','fignum',199); % % eq_signal_sd.eye(fsym,M,"fignum",103837); % % eq_signal_hd = PAMmapper(M, 0).quantize(eq_signal_sd); % rx_bits = PAMmapper(M,0,"eth_style",0).demap(eq_signal_hd); % [~,tot_err,ber_vnle(l,ri,realiz),a] = calc_ber(rx_bits.signal,Tx_bits{l}.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); % burst_vnle = count_error_bursts(a, 10)./tot_err; % % % showLevelConfusionMatrix(eq_signal_hd,Symbols{l},"M",M,"fignum",200,"displayname",'bla'); % % showLevelScatter(eq_signal_sd,Symbols{l},"displayname",'VNLE Out','f_sym',fsym,'fignum',201); % % show2Dconstellation(eq_signal_sd,Symbols{l},"displayname",'VNLE Out','fignum',2241); % % fprintf('CH %d :BER VNLE: %.2e \n',l,ber_vnle(l,ri,realiz)); % fprintf('CH %d :NGMI VNLE: %.2f \n',l,gmi_vnle_bitwise(l,ri,realiz)./m); end end end figure();hold on; cols = linspecer(N); for l = 1:N % plot(rop,mean(squeeze(ber_vnle(l,:,:)),2,'omitnan'),'Marker','*','DisplayName',sprintf('Ch: %d',wavelengthplan(l))) plot(rop,squeeze(ber_vnle(l,:,:)),'Marker','*','DisplayName',sprintf('Ch: %d',wavelengthplan(l)),'Color',cols(l,:),'HandleVisibility','on') end yline([3.8e-3,2.2e-4],'HandleVisibility','off'); ylabel('BER'); xlabel('ROP') title('BER vs. ROP'); set(gca, 'XScale', 'linear', ... 'YScale', 'log', ... 'TickLabelInterpreter', 'latex', ... 'FontSize', 11); res = struct(); % placeholder res.gmi_vnle_bitwise = gmi_vnle_bitwise; % placeholder res.snr_vnle = snr_vnle; res.ber_vnle = ber_vnle; % --------------------------------------------- % 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); % Save results save(fullfile(output_root, fname), 'res', '-v7.3'); fprintf('Saved results to: %s\n', fullfile(output_root, fname)); % --- save as PNG --- outname = fullfile(output_root, 'BER_vs_ROP.png'); % saves to current folder print(gcf, outname, '-dpng', '-r300'); % 300 dpi fprintf('Saved figure to %s\n', outname); function dispersion_vector = getDispersionVector(N, D, ref_zdw, randomize_ZDW, randomkey) % 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.09; % 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