Faster DP_fiber with GPU processing... Run test\gpu_cpu_comparison.m to see difference on your setup

This commit is contained in:
silas (home)
2026-02-01 13:53:20 +01:00
parent 67689bb70f
commit fba7cbdba2
15 changed files with 1281 additions and 217 deletions

219
test/bayesopt_ffe_tuning.m Normal file
View File

@@ -0,0 +1,219 @@
%% Bayesian Optimization for FFE Parameter Tuning
% This script uses bayesopt to find optimal mu_dd and mu_tr values
% that minimize BER for the FFE equalizer.
clear; clc;
%% Setup - Same as gpu_processing_dpfiber.m
s.wavelengthplan = calcWavelengthPlan(4, 400e9, 1310);
link_length = 10;
s.pmd = 0.1;
s.gamma = 0.0023;
s.M = 4;
fsym = 112e9;
fdac = 2*fsym;
fadc = 120000000000;
s.random_key = 1;
% Laser / Modulator
vbias_rel = 0.5;
u_pi = 4.6;
vbias = -vbias_rel*u_pi;
laser_linewidth = 0e6;
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";
N = numel(s.wavelengthplan);
switch s.p
case "co"
pol_rot = 100.*ones(1,N);
d_local = 0;
end
f_plan = physconst('lightspeed')./(s.wavelengthplan.*1e-9);
margin = 25e12;
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 = -8; % Fixed ROP for optimization
%% Generate TX signals (run once)
fprintf('Generating TX signals...\n');
for l = 1:N
[Digi_sig,Symbols{l},Tx_bits{l}] = PAMsource( ...
"fsym",fsym,"M",s.M,"order",15,"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);
El_sig = El_sig.normalize("mode","oneone");
scaling = 0.6*(u_pi/2-abs(vbias-u_pi/2));
El_sig = El_sig .* scaling;
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);
signal_cell{l} = Polarization_Controller("mode","rot_power","desired_power",pol_rot(l)).process(Eml_out);
end
%% 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 = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power", ...
"amplification_db",s.p_launch+10*log10(N)).process(Opt_sig_wdm);
%% Fiber propagation
segment_length = 1;
nSegments = link_length/segment_length;
nSegments = round(nSegments);
zdw = 1310;
randomize_D = true;
Dvec = getDispersionVector(nSegments, d_local, zdw, randomize_D, s.random_key);
Opt_sig_wdm_fib = Opt_sig_wdm;
fprintf('Running fiber propagation...\n');
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",true).process(Opt_sig_wdm_fib);
end
%% Pre-process to get Rx_sig (do demux once)
fprintf('Pre-processing receiver chain...\n');
l = 1; % Use channel 1 for optimization
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_wdm_fib);
Opt_sig_rx = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power", ...
"amplification_db",s.rop).process(Opt_sig_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).process(Opt_sig_rx);
rx_bwl = 100e9;
PD_sig = Filter('filtdegree',4,"f_cutoff",rx_bwl,"fs",fdac*kover, ...
"filterType",filtertypes.butterworth,"active",true).process(PD_sig);
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);
Scpe_sig_2sps = Scpe_sig.resample("fs_out",2*fsym);
[~, 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");
fprintf('Receiver pre-processing complete. Ready for optimization.\n\n');
%% Define the objective function for bayesopt
function ber = ffe_objective(params, Rx_sig, Symbols_l, Tx_bits_l, M, duob_mode)
mu_dd = params.mu_dd;
mu_tr = params.mu_tr;
try
eq_ffe = FFE("epochs_tr", 5, "epochs_dd", 2, "len_tr", 2^13, ...
"mu_dd", mu_dd, "mu_tr", mu_tr, ...
"order", 50, "sps", 2, "decide", 0, ...
"adaption", adaption_method.nlms, "dd_mode", 1);
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);
ber = ffe_results.metrics.BER;
if ber == 0
ber = 1e-10;
end
if ~isfinite(ber)
ber = 0.5;
end
fprintf(' mu_dd=%.4e, mu_tr=%.4e -> BER=%.4e\n', mu_dd, mu_tr, ber);
catch ME
fprintf(' mu_dd=%.4e, mu_tr=%.4e -> FAILED (%s)\n', mu_dd, mu_tr, ME.message);
ber = 0.5;
end
end
%% Define optimizable variables
mu_dd_var = optimizableVariable('mu_dd', [1e-5, 0.1], 'Transform', 'log');
mu_tr_var = optimizableVariable('mu_tr', [1e-5, 0.1], 'Transform', 'log');
%% Run Bayesian Optimization
fprintf('========== Starting Bayesian Optimization ==========\n');
fprintf('Optimizing mu_dd and mu_tr to minimize BER\n');
fprintf('Search range: mu_dd=[1e-5, 0.1], mu_tr=[1e-5, 0.1]\n\n');
objective_fn = @(params) ffe_objective(params, Rx_sig, Symbols{l}, Tx_bits{l}, s.M, duob_mode);
results = bayesopt(objective_fn, [mu_dd_var, mu_tr_var], ...
'MaxObjectiveEvaluations', 30, ...
'AcquisitionFunctionName', 'expected-improvement-plus', ...
'IsObjectiveDeterministic', false, ...
'ExplorationRatio', 0.5, ...
'Verbose', 1, ...
'PlotFcn', []);
%% Display Results
fprintf('\n========== FFE Optimization Complete ==========\n');
fprintf('Best FFE parameters found:\n');
fprintf(' mu_dd = %.6e\n', results.XAtMinObjective.mu_dd);
fprintf(' mu_tr = %.6e\n', results.XAtMinObjective.mu_tr);
fprintf(' BER = %.6e\n', results.MinObjective);
%% Verify with optimal parameters
fprintf('\nVerifying optimal FFE parameters...\n');
best_mu_dd = results.XAtMinObjective.mu_dd;
best_mu_tr = results.XAtMinObjective.mu_tr;
eq_ffe_best = FFE("epochs_tr", 5, "epochs_dd", 2, "len_tr", 2^13, ...
"mu_dd", best_mu_dd, "mu_tr", best_mu_tr, ...
"order", 50, "sps", 2, "decide", 0, ...
"adaption", adaption_method.nlms, "dd_mode", 1);
ffe_results_best = ffe(eq_ffe_best, s.M, Rx_sig, Symbols{l}, Tx_bits{l}, ...
"precode_mode", duob_mode, ...
'showAnalysis', 1, ...
"postFFE", [], ...
"eth_style_symbol_mapping", 0);
fprintf('\nFinal FFE BER with optimal parameters: %.6e\n', ffe_results_best.metrics.BER);

349
test/gpu_cpu_comparison.m Normal file
View File

@@ -0,0 +1,349 @@
%% GPU vs CPU Comparison Test for DP_Fiber
% This script runs the fiber simulation with and without GPU acceleration
% and compares the numerical results.
% clear; clc;
%% Setup (same as gpu_processing_dpfiber.m but simplified)
s.wavelengthplan = calcWavelengthPlan(4, 400e9, 1310);
link_length = 10;
s.pmd = 0.1;
s.gamma = 0.0023;
s.M = 4;
fsym = 112e9;
fdac = 2*fsym;
fadc = 120000000000;
s.random_key = 1;
% Laser / 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";
N = numel(s.wavelengthplan);
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 = 5e12;
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;
%% TX per channel
for l = 1:N
[Digi_sig,Symbols{l},Tx_bits{l}] = PAMsource( ...
"fsym",fsym,"M",s.M,"order",15,"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);
clear Digi_sig
El_sig = El_sig.normalize("mode","oneone");
scaling = 0.6*(u_pi/2-abs(vbias-u_pi/2));
El_sig = El_sig .* scaling;
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);
clear El_sig
signal_cell{l} = Polarization_Controller("mode","rot_power","desired_power",pol_rot(l)).process(Eml_out);
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 = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power", ...
"amplification_db",s.p_launch+10*log10(N)).process(Opt_sig_wdm);
% Save input for both runs
Opt_sig_input = 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;
if nSegments > 0
Dvec = getDispersionVector(nSegments, d_local, zdw, randomize_D, s.random_key);
else
Dvec = [];
end
%% Run WITHOUT GPU
fprintf('\n========== Running WITHOUT GPU (CPU) ==========\n');
Opt_sig_cpu = Opt_sig_input;
tic;
for seg = 1:nSegments
fprintf('CPU Segment %d/%d \n',seg, nSegments);
Opt_sig_cpu = 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",false).process(Opt_sig_cpu);
end
time_cpu = toc;
fprintf('CPU Time: %.3f seconds\n', time_cpu);
%% Run WITH GPU (Double Precision)
fprintf('\n========== Running WITH GPU (Double Precision) ==========\n');
Opt_sig_gpu_double = Opt_sig_input;
tic;
for seg = 1:nSegments
fprintf('GPU-Double Segment %d/%d \n',seg, nSegments);
Opt_sig_gpu_double = 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",false).process(Opt_sig_gpu_double);
end
time_gpu_double = toc;
fprintf('GPU Double Time: %.3f seconds\n', time_gpu_double);
%% Run WITH GPU (Single Precision)
fprintf('\n========== Running WITH GPU (Single Precision) ==========\n');
Opt_sig_gpu_single = Opt_sig_input;
tic;
for seg = 1:nSegments
fprintf('GPU-Single Segment %d/%d \n',seg, nSegments);
Opt_sig_gpu_single = 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",true).process(Opt_sig_gpu_single);
end
time_gpu_single = toc;
fprintf('GPU Single Time: %.3f seconds\n', time_gpu_single);
%% Compare Results
fprintf('\n========== Numerical Comparison ==========\n');
sig_cpu = Opt_sig_cpu.signal;
sig_gpu_double = Opt_sig_gpu_double.signal;
sig_gpu_single = Opt_sig_gpu_single.signal;
% Check dimensions
fprintf('CPU signal size: [%d x %d]\n', size(sig_cpu,1), size(sig_cpu,2));
fprintf('GPU Double signal size: [%d x %d]\n', size(sig_gpu_double,1), size(sig_gpu_double,2));
fprintf('GPU Single signal size: [%d x %d]\n', size(sig_gpu_single,1), size(sig_gpu_single,2));
% CPU vs GPU Double
fprintf('\n--- CPU vs GPU Double ---\n');
diff_cpu_double = abs(sig_cpu - sig_gpu_double);
max_diff_cpu_double = max(diff_cpu_double(:));
mean_diff_cpu_double = mean(diff_cpu_double(:));
rel_diff_cpu_double = max_diff_cpu_double / max(abs(sig_cpu(:)));
fprintf('Max absolute difference: %.6e\n', max_diff_cpu_double);
fprintf('Mean absolute difference: %.6e\n', mean_diff_cpu_double);
fprintf('Max relative difference: %.6e\n', rel_diff_cpu_double);
% CPU vs GPU Single
fprintf('\n--- CPU vs GPU Single ---\n');
diff_cpu_single = abs(sig_cpu - sig_gpu_single);
max_diff_cpu_single = max(diff_cpu_single(:));
mean_diff_cpu_single = mean(diff_cpu_single(:));
rel_diff_cpu_single = max_diff_cpu_single / max(abs(sig_cpu(:)));
fprintf('Max absolute difference: %.6e\n', max_diff_cpu_single);
fprintf('Mean absolute difference: %.6e\n', mean_diff_cpu_single);
fprintf('Max relative difference: %.6e\n', rel_diff_cpu_single);
% GPU Double vs GPU Single
fprintf('\n--- GPU Double vs GPU Single ---\n');
diff_double_single = abs(sig_gpu_double - sig_gpu_single);
max_diff_double_single = max(diff_double_single(:));
mean_diff_double_single = mean(diff_double_single(:));
rel_diff_double_single = max_diff_double_single / max(abs(sig_gpu_double(:)));
fprintf('Max absolute difference: %.6e\n', max_diff_double_single);
fprintf('Mean absolute difference: %.6e\n', mean_diff_double_single);
fprintf('Max relative difference: %.6e\n', rel_diff_double_single);
% Check tolerances
fprintf('\n--- Tolerance Check ---\n');
tol_double = 1e-10;
tol_single = 1e-5; % Single precision has ~7 significant digits
if max_diff_cpu_double < tol_double
fprintf(' CPU vs GPU Double: EQUIVALENT (diff < %.0e)\n', tol_double);
else
fprintf(' CPU vs GPU Double: DIFFER beyond tolerance (%.0e)\n', tol_double);
end
if max_diff_cpu_single < tol_single
fprintf(' CPU vs GPU Single: ACCEPTABLE (diff < %.0e)\n', tol_single);
else
fprintf(' CPU vs GPU Single: Precision loss detected (diff = %.2e, tol = %.0e)\n', max_diff_cpu_single, tol_single);
end
% Performance comparison
fprintf('\n========== Performance Summary ==========\n');
fprintf('CPU Time: %.3f s\n', time_cpu);
fprintf('GPU Double Time: %.3f s\n', time_gpu_double);
fprintf('GPU Single Time: %.3f s\n', time_gpu_single);
fprintf('\n');
fprintf('Speedup (GPU Double vs CPU): %.2fx\n', time_cpu/time_gpu_double);
fprintf('Speedup (GPU Single vs CPU): %.2fx\n', time_cpu/time_gpu_single);
fprintf('Speedup (GPU Single vs GPU Double): %.2fx\n', time_gpu_double/time_gpu_single);
%% ========== BER Comparison ==========
% Process each fiber output through simplified receiver to check if
% single-precision affects actual BER performance
fprintf('\n========== BER Comparison ==========\n');
fprintf('Processing signals through receiver chain...\n');
% Receiver parameters
rop = -7; % Received optical power [dBm]
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);
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 mode for channel 1
l = 4; % Use first channel for comparison
fprintf('Processing CPU result...\n');
ber_cpu = process_receiver(Opt_sig_cpu, l, Symbols, Tx_bits, ...
fdac, kover, upsample_pow, fsym, fadc, rop, s, len_tr, mu_dc, mu_ffe, mu_dfe, duob_mode);
fprintf('Processing GPU Double result...\n');
ber_gpu_double = process_receiver(Opt_sig_gpu_double, l, Symbols, Tx_bits, ...
fdac, kover, upsample_pow, fsym, fadc, rop, s, len_tr, mu_dc, mu_ffe, mu_dfe, duob_mode);
fprintf('Processing GPU Single result...\n');
ber_gpu_single = process_receiver(Opt_sig_gpu_single, l, Symbols, Tx_bits, ...
fdac, kover, upsample_pow, fsym, fadc, rop, s, len_tr, mu_dc, mu_ffe, mu_dfe, duob_mode);
% Display BER results
fprintf('\n========== BER Results (Channel %d, ROP = %d dBm) ==========\n', l, rop);
fprintf('CPU: BER = %.4e\n', ber_cpu);
fprintf('GPU Double: BER = %.4e\n', ber_gpu_double);
fprintf('GPU Single: BER = %.4e\n', ber_gpu_single);
fprintf('\n--- BER Comparison ---\n');
if ber_cpu == 0 && ber_gpu_double == 0 && ber_gpu_single == 0
fprintf(' All BERs are zero (no errors detected)\n');
else
ber_diff_double = abs(ber_cpu - ber_gpu_double);
ber_diff_single = abs(ber_cpu - ber_gpu_single);
fprintf('|BER_cpu - BER_gpu_double| = %.4e\n', ber_diff_double);
fprintf('|BER_cpu - BER_gpu_single| = %.4e\n', ber_diff_single);
if ber_diff_double < 1e-6 && ber_diff_single < 1e-6
fprintf(' BER differences are negligible\n');
elseif ber_diff_single > ber_diff_double * 10
fprintf(' Single precision shows measurable BER impact\n');
else
fprintf(' BER differences within acceptable range\n');
end
end
fprintf('\n========== Test Complete ==========\n');

View File

@@ -0,0 +1,265 @@
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');

View File

@@ -0,0 +1,70 @@
% 1. Setup Data OUTSIDE the timer
d = gpuDevice;
N = 10000;
fprintf('Preparing data...\n');
A_cpu_double = rand(N, N); % Create double on CPU
A_cpu_single = single(A_cpu_double); % Create single on CPU
% Warmup run (wakes up the GPU from idle state)
A_warm = gpuArray.rand(1000, 1000, 'single');
B_warm = A_warm * A_warm;
wait(d);
fprintf('------------------------------------------------\n');
% TEST 1: Double Precision (The "Slow" way)
% We move data to GPU first so we only measure calculation time
A_gpu_double = gpuArray(A_cpu_double);
wait(d); % Ensure transfer is done before starting timer
fprintf('Running DOUBLE precision test... ');
tic;
B_gpu = A_gpu_double * A_gpu_double;
wait(d); % FORCE MATLAB TO WAIT FOR GPU
time_double = toc;
fprintf('Done.\n');
fprintf('Double Precision Time: %.4f seconds\n', time_double);
% TEST 2: Single Precision (The "Fast" way)
A_gpu_single = gpuArray(A_cpu_single);
wait(d); % Ensure transfer is done
fprintf('Running SINGLE precision test... ');
tic;
B_gpu = A_gpu_single * A_gpu_single;
wait(d); % FORCE MATLAB TO WAIT FOR GPU
time_single = toc;
fprintf('Done.\n');
fprintf('Single Precision Time: %.4f seconds\n', time_single);
% Calculate Speedup
fprintf('------------------------------------------------\n');
fprintf('Speedup Factor using single precision: %.2fx\n', time_double / time_single);
%
% Create 10,000 small matrices (10x10) stacked in a 3D array
A_stack = gpuArray.rand(10, 10, 10000, 'single');
B_stack = gpuArray.rand(10, 10, 10000, 'single');
% BAD: Looping (GPU overhead kills you)
tic;
for i=1:10000
C(:,:,i) = A_stack(:,:,i) * B_stack(:,:,i);
end
wait(d);
loop_time = toc;
% GOOD: Pagefun (Executes all 10,000 mults simultaneously)
tic;
C_stack = pagefun(@mtimes, A_stack, B_stack);
wait(d);
pagefun_time = toc;
fprintf('------------------------------------------------\n');
fprintf('Speedup Factor using pagefun: %.2fx\n', loop_time / pagefun_time);
fprintf('Total VRAM: %.2f GB\n', d.TotalMemory / 1e9);
fprintf('Available VRAM: %.2f GB\n', d.AvailableMemory / 1e9);
fprintf('Usage: %.1f%%\n', 100 * (1 - d.AvailableMemory / d.TotalMemory));