Merge branch 'main' of cau-git.rz.uni-kiel.de:nt/mitarbeiter/silas/imdd_simulation

This commit is contained in:
sutef391
2025-12-22 07:29:09 +01:00
164 changed files with 14366 additions and 3330 deletions

View File

@@ -15,9 +15,9 @@ m = floor(log2(M)*10)/10;
fsym = 224e9;
fdac = 2*fsym;
fadc = 2*fsym;
random_key = 2;
s.random_key = 100;
% Laser / Modulator
% Laser / s.Modulator
vbias_rel = 0.5;
u_pi = 3.2;
vbias = -vbias_rel*u_pi;
@@ -47,8 +47,8 @@ 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);
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;
@@ -57,38 +57,38 @@ 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;
s.f_opt = fdac*kover*upsample_pow;
s.f_opt_nyq = s.f_opt/2;
signal_cell = {};
Symbols = {};
Tx_bits = {};
rop = -8.25:0.75:0;
s.rop = -6:0.75:-0.75;
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);
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:num_realiz
for realiz = 1:s.num_realiz
parfor l = 1:N
[Digi_sig,Symbols{l},Tx_bits{l}] = PAMsource(...
"fsym",fsym,"M",M,"order",18,"useprbs",0,...
"fsym",fsym,"M",s.M,"order",18,"useprbs",0,...
"fs_out",fdac,...
"applyclipping",0,"clipfactor",1.5,...
"applypulseform",apply_pulsef,"pulseformer",Pform,...
"randkey",random_key+l+realiz,...
"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",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 = s.M8199B("kover",kover).process(Digi_sig);
% El_sig.spectrum("fignum",101,"displayname",'bla','normalizeTo0dB',0,'lambda0_nm',1310,'useWavelengthAxis',0);
%%%%% Electrical Driver Amplifier %%%%%%
@@ -96,10 +96,10 @@ for realiz = 1:num_realiz
% 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);
%%%%% 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",30).process(Eml_out);
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,...
@@ -119,12 +119,12 @@ for realiz = 1:num_realiz
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
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",gamma,"lambda",zdw,"n_waveplates",10,"SS_dphimax",0.01,...
"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
@@ -133,12 +133,12 @@ for realiz = 1:num_realiz
% 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)
% 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)
parfor ri = 1:length(rop)
for ri = 1:length(s.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_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);
@@ -147,9 +147,9 @@ for realiz = 1:num_realiz
%%%%%% 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 = 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);
% PD_sig.spectrum("fignum",222,"displayname",'bla','normalizeTo0dB',1);
%%%%%% Low-pass RX (PD, El. Connectors and Scope %%%%%%
rx_bwl = 100e9;
@@ -176,7 +176,7 @@ 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},...
ffe_results = ffe(eq_ffe,s.M,Rx_sig,Symbols{l},Tx_bits{l},...
"precode_mode",duob_mode,...
'showAnalysis',0,...
"postFFE",[],...
@@ -194,12 +194,12 @@ for realiz = 1:num_realiz
useviterbi = 0;
if useviterbi
mlse_ = MLSE_viterbi("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels);
mlse_ = MLSE_viterbi("duobinary_output",0,'M',s.M,'trellis_states',PAMmapper(s.M,0).levels);
else
mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels);
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_, M, Rx_sig, Symbols{l},Tx_bits{l}, ...
[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", [],...
@@ -209,18 +209,17 @@ for realiz = 1:num_realiz
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);
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",M,"trellis_states",PAMmapper(M,0).levels);
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_, M, Rx_sig, Symbols{l},Tx_bits{l}, ...
dbt_results = duobinary_target(eq_, mlse_db_, s.M, Rx_sig, Symbols{l},Tx_bits{l}, ...
"precode_mode", duob_mode, ...
'showAnalysis', 0,...
"postFFE", []);
@@ -232,6 +231,7 @@ for realiz = 1:num_realiz
end
res = struct();
res.settings = s;
res.ffe = output_ffe;
res.vnle = output_vnle;
res.mlse = output_mlse;
@@ -240,38 +240,13 @@ for realiz = 1:num_realiz
% 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
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,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, res.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, res.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');
xlabel('ROP')
title('BER vs. ROP');
set(gca, 'XScale', 'linear', ...
'YScale', 'log', ...
'TickLabelInterpreter', 'latex', ...
'FontSize', 11);
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
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.
% 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