diff --git a/projects/WDM/WDM_model.m b/projects/WDM/WDM_model.m index 77247d0..cdc8a74 100644 --- a/projects/WDM/WDM_model.m +++ b/projects/WDM/WDM_model.m @@ -4,7 +4,7 @@ run(fullfile(fileparts(mfilename('fullpath')),'WDM_settings.m')); num_realiz = 50; -wavelengthplan = calcWavelengthPlan(16,400e9,1310); +s.wavelengthplan = calcWavelengthPlan(16,400e9,1310); % wavelengthplan = [1295,1305,1315,1325]; link_length = 10; pmd = 0.1; @@ -12,7 +12,7 @@ gamma = 0.0023; M = 4; m = floor(log2(M)*10)/10; -fsym = 224e9; +fsym = 112e9; fdac = 2*fsym; fadc = 2*fsym; s.random_key = 100; @@ -86,7 +86,7 @@ for realiz = 1:s.num_realiz "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); + 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 = s.M8199B("kover",kover).process(Digi_sig); % El_sig.spectrum("fignum",101,"displayname",'bla','normalizeTo0dB',0,'lambda0_nm',1310,'useWavelengthAxis',0); diff --git a/projects/WDM/WDM_model_16x200.m b/projects/WDM/WDM_model_16x200.m new file mode 100644 index 0000000..6478aff --- /dev/null +++ b/projects/WDM/WDM_model_16x200.m @@ -0,0 +1,283 @@ +%%% Run parameters +% TX +% --- FIRST LINE: evaluate settings located beside this script --- +run(fullfile(fileparts(mfilename('fullpath')),'WDM_settings.m')); + +s.num_realiz = 2; +s.wavelengthplan = calcWavelengthPlan(16,400e9,1310); +% wavelengthplan = [1295,1305,1315,1325]; +link_length = 0; +pmd = 0.1; +gamma = 0.0023; + +s.M = 4; +m = floor(log2(s.M)*10)/10; +fsym = 112e9; +fdac = 2*fsym; +fadc = 120000000000; +s.random_key = 100; + +% 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 = 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); + +s.f_opt = fdac*kover*upsample_pow; +s.f_opt_nyq = s.f_opt/2; + +signal_cell = {}; +Symbols = {}; +Tx_bits = {}; + +s.rop = -6: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); + +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"); + El_sig = El_sig .* u_pi .* 0.5; + % 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",100).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",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; + + segment_length = 1; % km + nSegments = link_length/segment_length; + 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, s.random_key+realiz); + for seg = 1:nSegments + + Opt_sig_wdm_fib = DP_Fiber("L",segment_length,"D",Dvec(s),"Dpmd",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) + + for ri = 1:length(s.rop) + + %%%%%% ROP %%%%%% + Opt_sig_wdm_rx = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",s.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",s.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",80e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true); + + disp(fadc); + disp(fdac*kover) + + %%%%%% 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 + 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 + +function dispersion_vector = getDispersionVector(N, D, ref_zdw, randomize_ZDW, randomkey) +% s.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.07; % 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 diff --git a/projects/WDM/WDM_settings.m b/projects/WDM/WDM_settings.m index 68a04a3..c64071f 100644 --- a/projects/WDM/WDM_settings.m +++ b/projects/WDM/WDM_settings.m @@ -1,23 +1,6 @@ - - - - - - - - - - - - - - - - - % Add the imdd_simulation framework to the path if ispc @@ -57,8 +40,6 @@ fprintf('parpool up with %d workers; JobStorage=%s\n', p.NumWorkers, c.JobStorag - - % 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