Files
imdd_silas/Theory/GPU/gpu_processing_dpfiber.m
2026-03-25 10:57:48 +01:00

266 lines
8.2 KiB
Matlab
Raw Permalink Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
s.wavelengthplan = calcWavelengthPlan(16, 400e9, 1310);
N = numel(s.wavelengthplan);
link_length = 10;
s.pmd = 0.1;%0.1;
s.gamma = 0.0023;
s.M = 4;
fsym = 112e9;
fdac = 2*fsym;
fadc = 120000000000;
s.random_key = 1;
% Laser / s.Modulator
vbias_rel = 0.5;
u_pi = 4.6;
vbias = -vbias_rel*u_pi;
laser_linewidth = 0e6;
% DB Stuff
duob_mode = db_mode.no_db;
rcalpha = 0.05;
Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",16,"alpha",rcalpha);
s.chirpalpha = 0;
s.p_launch = 3;
s.p = "co";
switch s.p
case "co"
pol_rot = 100.*ones(1,N);
d_local = 0;
case "pair"
pol_rot = repmat([100,100,0,0],1,N/4);
d_local = 0;
case "alt"
pol_rot = repmat([100,0,100,0],1,N/4);
d_local = 0;
case "seg"
pol_rot = 100.*ones(1,N);
d_local = 3;
otherwise
error('Unknown fwm_mitigation_technique: %s', string(s.p));
end
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 = 4;
upsample_required = f_nyq./(fdac*kover/2);
upsample_pow = 2^nextpow2(upsample_required);
s.f_opt = fdac*kover*upsample_pow;
s.f_opt_nyq = s.f_opt/2;
s.rop = -10:1:0;
profile on
%% ---------- TX per channel ----------
for l = 1:N
[Digi_sig,Symbols{l},Tx_bits{l}] = PAMsource( ...
"fsym",fsym,"M",s.M,"order",17,"useprbs",0, ...
"fs_out",fdac, ...
"applyclipping",0,"clipfactor",1.5, ...
"applypulseform",1,"pulseformer",Pform, ...
"randkey",s.random_key+l, ...
"mrds_code",0,"mrds_blocklength",512,"duobinary_mode",duob_mode ...
).process();
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);
% Digi_sig not needed after AWG
clear Digi_sig
% Electrical Driver Amplifier
El_sig = El_sig.normalize("mode","oneone");
scaling = 0.6*(u_pi/2-abs(vbias-u_pi/2));
El_sig = El_sig .* scaling;
% 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,"alpha",s.chirpalpha).process(El_sig);
% El_sig not needed after EML
clear El_sig
signal_cell{l} = Polarization_Controller("mode","rot_power","desired_power",pol_rot(l)).process(Eml_out);
% Eml_out not needed after pol controller
clear Eml_out Lp_awg
end
disp('Signal generated for all channels.');
%% ---------- WDM mux + launch ----------
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.spectrum();
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_fib = Opt_sig_wdm;
segment_length = 1;
nSegments = link_length/segment_length;
if abs(nSegments - round(nSegments)) > 1e-12
error('fiber_length_km=%g must be an integer multiple of segment_length=%g km.', link_length, segment_length);
end
nSegments = round(nSegments);
zdw = 1310;
randomize_D = true;
% Guard for 0 km: avoid calling getDispersionVector(0,...) if it doesn't support it
if nSegments > 0
Dvec = getDispersionVector(nSegments, d_local, zdw, randomize_D, s.random_key);
else
Dvec = [];
end
for seg = 1:nSegments
fprintf('Segment %d/%d \n',seg, 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,"useGPU",true,"useSingle",1).process(Opt_sig_wdm_fib);
end
profile off
profile viewer
Opt_sig_wdm_fib.spectrum();
%% ========== BER Evaluation ==========
fprintf('\n========== BER Evaluation ==========\n');
fprintf('Processing signals through receiver chain...\n');
% Receiver parameters
len_tr = 4096; % Training length
mu_dc = 0.005;
mu_ffe = [0.0001 0.0008 0.001];
mu_dfe = 0.0004;
% Helper function to process through receiver and get BER
function ber = process_receiver(Opt_sig_fib, l, Symbols, Tx_bits, ...
fdac, kover, upsample_pow, fsym, fadc, rop, s, len_tr, mu_dc, mu_ffe, mu_dfe, duob_mode)
% Demux single channel
Opt_sig_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_fib);
% ROP amplifier
Opt_sig_rx = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power", ...
"amplification_db",rop).process(Opt_sig_demux{l});
% Photodiode
PD_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20, ...
"nep",1.8e-11,"randomkey",s.random_key+l).process(Opt_sig_rx);
% Low-pass filter
rx_bwl = 100e9;
PD_sig = Filter('filtdegree',4,"f_cutoff",rx_bwl,"fs",fdac*kover, ...
"filterType",filtertypes.butterworth,"active",true).process(PD_sig);
% Scope
Lp_scpe = Filter('filtdegree',4,"f_cutoff",80e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true);
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);
% Resample to 2 sps
Scpe_sig_2sps = Scpe_sig.resample("fs_out",2*fsym);
% Time sync
[~, Scpe_cell, ~, ~] = Scpe_sig_2sps.tsynch("reference", Symbols{l}, "fs_ref", fsym, "debug_plots", 0);
Rx_sig = Scpe_cell{1};
Rx_sig = Rx_sig.normalize("mode","rms");
% FFE Equalizer
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);
eq_ffe = FFE("epochs_tr",5,"epochs_dd",2,"len_tr",2^13,"mu_dd",6.624e-05,"mu_tr",0.058136,"order",50,"sps",2,"decide",0, "adaption",adaption_method.nlms,"dd_mode",1);
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);
ber = ffe_results.metrics.BER;
end
% Process each ROP value for selected channels using parfor
ber_results = zeros(length(s.rop), N);
% Flatten loop for parfor: iterate over all (ROP, channel) combinations
num_rop = length(s.rop);
rop_vals = s.rop;
ber_flat = zeros(num_rop * N, 1);
parfor idx = 1:(num_rop * N)
% Convert linear index to (ri, l) subscripts
ri = ceil(idx / N);
l = mod(idx - 1, N) + 1;
fprintf('ROP %d dBm, Channel %d/%d\n', rop_vals(ri), l, N);
ber_flat(idx) = process_receiver(Opt_sig_wdm_fib, l, Symbols, Tx_bits, ...
fdac, kover, upsample_pow, fsym, fadc, rop_vals(ri), s, len_tr, mu_dc, mu_ffe, mu_dfe, duob_mode);
end
% Reshape back to [num_rop × N] matrix
ber_results = reshape(ber_flat, [N, num_rop]).';
%% Display BER Results
fprintf('\n========== BER Results ==========\n');
fprintf('ROP [dBm] | ');
for l = 1:N
fprintf('Ch%d | ', l);
end
fprintf('\n');
for ri = 1:length(s.rop)
fprintf('%8d | ', s.rop(ri));
for l = 1:N
fprintf('%.2e | ', ber_results(ri, l));
end
fprintf('\n');
end
% Plot BER vs ROP
figure;
semilogy(s.rop, mean(ber_results, 2), '-o', 'LineWidth', 2);
hold on;
for l = 1:N
semilogy(s.rop, ber_results(:, l), '--', 'LineWidth', 1);
end
hold off;
xlabel('ROP [dBm]');
ylabel('BER');
title('BER vs Received Optical Power (GPU Single Precision)');
legend(['Mean', arrayfun(@(x) sprintf('Ch%d', x), 1:N, 'UniformOutput', false)]);
grid on;
fprintf('\n========== Test Complete ==========\n');