diff --git a/projects/WDM/WDM_model.m b/projects/WDM/WDM_model.m index 37a39fa..10e6c69 100644 --- a/projects/WDM/WDM_model.m +++ b/projects/WDM/WDM_model.m @@ -3,21 +3,26 @@ % --- FIRST LINE: evaluate settings located beside this script --- run(fullfile(fileparts(mfilename('fullpath')),'WDM_settings.m')); +num_realiz = 1; +% wavelengthplan = calcWavelengthPlan(16,400e9,1310); +wavelengthplan = [1295,1305,1315,1325]; +link_length = 2; +pmd = 0.1; +gamma = 0.0023; + + 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; +laser_linewidth = 0e6; % EQ SETTINGS vnle_order1 = 50; @@ -43,7 +48,6 @@ 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 @@ -63,14 +67,14 @@ 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); +output_ffe = cell(length(wavelengthplan),length(rop),num_realiz); +output_vnle = cell(length(wavelengthplan),length(rop),num_realiz); +output_mlse = cell(length(wavelengthplan),length(rop),num_realiz); +output_dbt = cell(length(wavelengthplan),length(rop),num_realiz); for realiz = 1:num_realiz + for l = 1:N [Digi_sig,Symbols{l},Tx_bits{l}] = PAMsource(... @@ -118,9 +122,9 @@ for realiz = 1:num_realiz 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,... + Opt_sig_wdm_fib = DP_Fiber("L",link_length/nSegments,"D",Dvec(s),"Dpmd",pmd,"Ds",0.07,... "beat_len",10,"corr_len",100,"dz",1,"manakov",0,... - "gamma",0.0023,"lambda",zdw,"n_waveplates",10,"SS_dphimax",0.01,... + "gamma",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 @@ -169,56 +173,85 @@ for realiz = 1:num_realiz - + % 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,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; + output_ffe{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); - + %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',M,'trellis_states',PAMmapper(M,0).levels); + else + mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); + end + + [vnle_results, mlse_results] = vnle_postfilter_mlse(eq_, pf_, mlse_, 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',M,'trellis_states',PAMmapper(M,0).levels); + else + mlse_db_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",M,"trellis_states",PAMmapper(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_, 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.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)); + 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') + % plot(rop,cellfun(@(c) c.metrics.BER, output_ffe(l,:), 'UniformOutput', true),'Marker','*','DisplayName',sprintf('Ch: %d',wavelengthplan(l)),'Color',cols(l,:),'HandleVisibility','on','LineStyle',':'); + plot(rop,cellfun(@(c) c.metrics.BER, output_vnle(l,:), 'UniformOutput', true),'Marker','x','DisplayName',sprintf('Ch: %d',wavelengthplan(l)),'Color',cols(l,:),'HandleVisibility','on','LineStyle','--') + plot(rop,cellfun(@(c) c.metrics.BER, output_mlse(l,:), 'UniformOutput', true),'Marker','o','DisplayName',sprintf('Ch: %d',wavelengthplan(l)),'Color',cols(l,:),'HandleVisibility','on','LineStyle','-') end yline([3.8e-3,2.2e-4],'HandleVisibility','off'); ylabel('BER'); @@ -228,38 +261,8 @@ 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)); +xlim([min(rop) max(rop)]) +ylim([1e-5 0.3]) % --- save as PNG --- outname = fullfile(output_root, 'BER_vs_ROP.png'); % saves to current folder @@ -279,7 +282,7 @@ function dispersion_vector = getDispersionVector(N, D, ref_zdw, randomize_ZDW, r % constants (matching the Python code) meanLambda_nm = 1310; % center wavelength sigma_nm = 2; % ZDW sigma - Dslope = 0.09; % ps/(nm·km) per nm detuning + Dslope = 0.07; % ps/(nm·km) per nm detuning % random ZDW-induced dispersion offset if randomize_ZDW diff --git a/projects/WDM/WDM_settings.m b/projects/WDM/WDM_settings.m index a0bfe85..68a04a3 100644 --- a/projects/WDM/WDM_settings.m +++ b/projects/WDM/WDM_settings.m @@ -52,4 +52,29 @@ 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); \ No newline at end of file +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 +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); \ No newline at end of file