From 01a38814554d2074153519b39e3b7dddd8f02075 Mon Sep 17 00:00:00 2001 From: "silas (home)" Date: Wed, 7 Jan 2026 08:30:14 +0100 Subject: [PATCH] few changes in WDM stuff cleaned up the base system... but its not yet a good minimal example... --- projects/IMDD_base_system/imdd_it.m | 532 +++++++++--------- projects/IMDD_base_system/imdd_model.m | 215 +++---- .../submit_handle.m | 0 projects/WDM/WDM_auswertung.m | 333 +---------- projects/WDM/WDM_model_10km_queue.m | 102 +++- projects/WDM/plot_BER_vs_ROP.m | 199 +++++++ projects/WDM/plot_FEC_violin.m | 213 +++++++ 7 files changed, 841 insertions(+), 753 deletions(-) rename projects/{Messung_Zürich => IMDD_base_system}/submit_handle.m (100%) create mode 100644 projects/WDM/plot_BER_vs_ROP.m create mode 100644 projects/WDM/plot_FEC_violin.m diff --git a/projects/IMDD_base_system/imdd_it.m b/projects/IMDD_base_system/imdd_it.m index 748708f..c9d9cf4 100644 --- a/projects/IMDD_base_system/imdd_it.m +++ b/projects/IMDD_base_system/imdd_it.m @@ -1,117 +1,133 @@ -% basePath = 'C:\Users\Silas\Documents\MATLAB\Datensätze\sioe_labor\'; -% db = DBHandler("pathToDB",[basePath,'silas_labor.db']); + if 1 uloops = struct; uloops.precomp = [1]; - uloops.db_precode = [0]; - uloops.bitrate = [224].*1e9; %[300,330,360,390,420,450,480] [224,336,360,390,420,448] for MPI + uloops.bitrate = [300].*1e9; %[300,330,360,390,420,450,480] [224,336,360,390,420,448] for MPI % uloops.laser_wavelength = [1293,1297.5,1302,1306.5,1310,1313.4,1318,1322.7,1327.4]; - uloops.laser_wavelength = [1310]; + uloops.laser_wavelength = [1293]; uloops.M = [4]; - uloops.link_length = [1]; % 1,2,3,5,6,8,10 - % uloops.interference_attenuation = [0,3,6,9,12,15,18,21,24,27,30,45]; + uloops.link_length = [0:2:10]; % 1,2,3,5,6,8,10 + uloops.alpha = [0.5,-0.5,0]; + wh = DataStorage(uloops); wh.addStorage("ber"); - % wh = submit_simulations(wh,"parallel",0,"simulation_mode",0); - wh = submit_handle(@imdd_model,wh,"serial",0); + wh = submit_handle(@imdd_model,wh,"parallel",1); end -wh_ana = wh_master; - -cols = cbrewer2('Paired',8); - -figure() - -for precomp = [0,1] - wavelength=uloops.laser_wavelength; - for m = [6] - - baudrate = wh_ana.parameter.bitrate.values; - - - %VNLE - precode = 1; - a = wh_ana.getStoValue('ber',precomp, precode, baudrate , wavelength, m, uloops.link_length); - ber_vnle_pc = cellfun(@(x) x.vnle_pf_package{1,1}.ber_vnle, a); - %MLSE - ber_mlse_pc = cellfun(@(x) x.vnle_pf_package{1,1}.ber_mlse, a); - %DB - ber_dbtgt_pc = cellfun(@(x) x.dbtgt_package{1,1}.ber, a); - - precode = 0; - a = wh_ana.getStoValue('ber',precomp, precode, baudrate , wavelength, m, uloops.link_length); - ber_vnle = cellfun(@(x) x.vnle_pf_package{1,1}.ber_vnle, a); - %MLSE - ber_mlse = cellfun(@(x) x.vnle_pf_package{1,1}.ber_mlse, a); - %DB - ber_dbtgt = cellfun(@(x) x.dbtgt_package{1,1}.ber, a); - - if precomp - legndname1 = ['Pre-Emphasis']; - else - legndname1 = ['No Pre-Emphasis']; - end - - - baudrate = floor( uloops.bitrate.*1e-9 ./log2(m) ) .* 2.5 .* 1e9; - - subplot(1,3,1) - hold on - title(sprintf('%d km | %d nm | PAM %d',uloops.link_length,wavelength,m)); - title(sprintf('PAM %d',m)); - plot(baudrate.*1e-9,ber_vnle,'DisplayName',['Pre-Emphasis: ', num2str(precomp), '| Diff.-Code: 0'],'Color',cols(1+precomp,:),'LineStyle','-','HandleVisibility','on'); - plot(baudrate.*1e-9,ber_vnle_pc,'DisplayName',['Pre-Emphasis: ', num2str(precomp), '| Diff.-Code: 1'],'Color',cols(1+precomp,:),'LineStyle','-.','HandleVisibility','on'); - xticks(baudrate.*1e-9); - set(gca, 'YScale', 'log'); - ylim([5e-5 0.4]); - xlim([min(baudrate(2).*1e-9), max(baudrate.*1e-9) ]); - yline([3.8e-3, 2e-2],'HandleVisibility','off'); - legend - beautifyBERplot() - xlabel('Bit Rate in Gbps'); - ylabel('BER'); - - - - subplot(1,3,2) - hold on - % title(sprintf('%d km | %d nm | PAM %d',uloops.link_length,wavelength,m)); - title(sprintf('PAM %d',m)); - plot(baudrate.*1e-9,ber_dbtgt,'DisplayName',['Pre-Emphasis: ', num2str(precomp), '| Diff.-Code: 0'],'Color',cols(5+precomp,:),'LineStyle','-','HandleVisibility','on'); - plot(baudrate.*1e-9,ber_dbtgt_pc,'DisplayName',['Pre-Emphasis: ', num2str(precomp), '| Diff.-Code: 1'],'Color',cols(5+precomp,:),'LineStyle','-.','HandleVisibility','on'); - xticks(baudrate.*1e-9); - set(gca, 'YScale', 'log'); - ylim([5e-5 0.4]); - xlim([min(baudrate(2).*1e-9), max(baudrate.*1e-9) ]); - yline([3.8e-3, 2e-2],'HandleVisibility','off'); - legend - beautifyBERplot() - xlabel('Bit Rate in Gbps'); - ylabel('BER'); - - - subplot(1,3,3) - hold on - % title(sprintf('%d km | %d nm | PAM %d',uloops.link_length,wavelength,m)); - title(sprintf('PAM %d',m)); - plot(baudrate.*1e-9,ber_mlse,'DisplayName',['Pre-Emphasis: ', num2str(precomp), '| Diff.-Code: 0'],'Color',cols(3+precomp,:),'LineStyle','-','HandleVisibility','on'); - plot(baudrate.*1e-9,ber_mlse_pc,'DisplayName',['Pre-Emphasis: ', num2str(precomp), '| Diff.-Code: 1'],'Color',cols(3+precomp,:),'LineStyle','-.','HandleVisibility','on'); - xticks(baudrate.*1e-9); - set(gca, 'YScale', 'log'); - ylim([5e-5 0.4]); - xlim([min(baudrate(2).*1e-9), max(baudrate.*1e-9) ]); - yline([3.8e-3, 2e-2],'HandleVisibility','off'); - legend - beautifyBERplot() - xlabel('Bit Rate in Gbps'); - ylabel('BER'); - end +%% +figure +hold on +for alpha = uloops.alpha + a=wh.getStoValue('ber',1, [300].*1e9 , 1293, 4, uloops.link_length,alpha); + ffe = cellfun(@(x) x.ffe_results.metrics.BER, a); + plot(uloops.link_length,ffe,'DisplayName',sprintf('Alpha: %d',alpha),'LineStyle','-','HandleVisibility','on'); end -% -% + +set(gca, 'YScale', 'log'); +ylim([5e-5 0.4]); +yline([3.8e-3, 2e-2],'HandleVisibility','off'); +legend +beautifyBERplot() +ylabel('BER'); + + +% +% wh_ana = wh_master; +% +% cols = cbrewer2('Paired',8); +% +% figure() +% +% for precomp = [0,1] +% wavelength=uloops.laser_wavelength; +% for m = [6] +% +% baudrate = wh_ana.parameter.bitrate.values; +% +% +% %VNLE +% precode = 1; +% a = wh_ana.getStoValue('ber',precomp, precode, baudrate , wavelength, m, uloops.link_length); +% ber_vnle_pc = cellfun(@(x) x.vnle_pf_package{1,1}.ber_vnle, a); +% %MLSE +% ber_mlse_pc = cellfun(@(x) x.vnle_pf_package{1,1}.ber_mlse, a); +% %DB +% ber_dbtgt_pc = cellfun(@(x) x.dbtgt_package{1,1}.ber, a); +% +% precode = 0; +% a = wh_ana.getStoValue('ber',precomp, precode, baudrate , wavelength, m, uloops.link_length); +% ber_vnle = cellfun(@(x) x.vnle_pf_package{1,1}.ber_vnle, a); +% %MLSE +% ber_mlse = cellfun(@(x) x.vnle_pf_package{1,1}.ber_mlse, a); +% %DB +% ber_dbtgt = cellfun(@(x) x.dbtgt_package{1,1}.ber, a); +% +% if precomp +% legndname1 = ['Pre-Emphasis']; +% else +% legndname1 = ['No Pre-Emphasis']; +% end +% +% +% baudrate = floor( uloops.bitrate.*1e-9 ./log2(m) ) .* 2.5 .* 1e9; +% +% subplot(1,3,1) +% hold on +% title(sprintf('%d km | %d nm | PAM %d',uloops.link_length,wavelength,m)); +% title(sprintf('PAM %d',m)); +% plot(baudrate.*1e-9,ber_vnle,'DisplayName',['Pre-Emphasis: ', num2str(precomp), '| Diff.-Code: 0'],'Color',cols(1+precomp,:),'LineStyle','-','HandleVisibility','on'); +% plot(baudrate.*1e-9,ber_vnle_pc,'DisplayName',['Pre-Emphasis: ', num2str(precomp), '| Diff.-Code: 1'],'Color',cols(1+precomp,:),'LineStyle','-.','HandleVisibility','on'); +% xticks(baudrate.*1e-9); +% set(gca, 'YScale', 'log'); +% ylim([5e-5 0.4]); +% xlim([min(baudrate(2).*1e-9), max(baudrate.*1e-9) ]); +% yline([3.8e-3, 2e-2],'HandleVisibility','off'); +% legend +% beautifyBERplot() +% xlabel('Bit Rate in Gbps'); +% ylabel('BER'); +% +% +% +% subplot(1,3,2) +% hold on +% % title(sprintf('%d km | %d nm | PAM %d',uloops.link_length,wavelength,m)); +% title(sprintf('PAM %d',m)); +% plot(baudrate.*1e-9,ber_dbtgt,'DisplayName',['Pre-Emphasis: ', num2str(precomp), '| Diff.-Code: 0'],'Color',cols(5+precomp,:),'LineStyle','-','HandleVisibility','on'); +% plot(baudrate.*1e-9,ber_dbtgt_pc,'DisplayName',['Pre-Emphasis: ', num2str(precomp), '| Diff.-Code: 1'],'Color',cols(5+precomp,:),'LineStyle','-.','HandleVisibility','on'); +% xticks(baudrate.*1e-9); +% set(gca, 'YScale', 'log'); +% ylim([5e-5 0.4]); +% xlim([min(baudrate(2).*1e-9), max(baudrate.*1e-9) ]); +% yline([3.8e-3, 2e-2],'HandleVisibility','off'); +% legend +% beautifyBERplot() +% xlabel('Bit Rate in Gbps'); +% ylabel('BER'); +% +% +% subplot(1,3,3) +% hold on +% % title(sprintf('%d km | %d nm | PAM %d',uloops.link_length,wavelength,m)); +% title(sprintf('PAM %d',m)); +% plot(baudrate.*1e-9,ber_mlse,'DisplayName',['Pre-Emphasis: ', num2str(precomp), '| Diff.-Code: 0'],'Color',cols(3+precomp,:),'LineStyle','-','HandleVisibility','on'); +% plot(baudrate.*1e-9,ber_mlse_pc,'DisplayName',['Pre-Emphasis: ', num2str(precomp), '| Diff.-Code: 1'],'Color',cols(3+precomp,:),'LineStyle','-.','HandleVisibility','on'); +% xticks(baudrate.*1e-9); +% set(gca, 'YScale', 'log'); +% ylim([5e-5 0.4]); +% xlim([min(baudrate(2).*1e-9), max(baudrate.*1e-9) ]); +% yline([3.8e-3, 2e-2],'HandleVisibility','off'); +% legend +% beautifyBERplot() +% xlabel('Bit Rate in Gbps'); +% ylabel('BER'); +% end +% end +% % +% % % cols = linspecer(7);%cbrewer2('Set2',10); % % @@ -163,163 +179,163 @@ end % end - -m = 6; -ir = [2,2.5,3]; -cols = linspecer(6); -baudrate_gather = []; -for i = 1:3 - m = uloops.M(i); - - %%% GET VNLE VALS - precode = 0; - precomp = 1; - a = wh_master.getStoValue('ber',precomp, precode, uloops.bitrate , uloops.laser_wavelength, m, uloops.link_length); - ber_vnle(i,:) = cellfun(@(x) x.vnle_pf_package{1,1}.ber_vnle, a); - - %%% GET DB VALS - precode = 1; - precomp = 0; - a = wh_master.getStoValue('ber',precomp, precode, uloops.bitrate , uloops.laser_wavelength, m, uloops.link_length); - ber_db(i,:) = cellfun(@(x) x.dbtgt_package{1,1}.ber, a); - - %%% GET MLSE VALS - precode = 0; - precomp = 0; - a = wh_master.getStoValue('ber',precomp, precode, uloops.bitrate , uloops.laser_wavelength, m, uloops.link_length); - ber_mlse(i,:) = cellfun(@(x) x.vnle_pf_package{1,1}.ber_mlse, a); - - inf_rate_pam(i,:) = cellfun(@(x) x.vnle_pf_package{1,1}.air, a); - inf_rate_pam(i,:) = inf_rate_pam(i,:)./log2(m); - - bitrate = floor( uloops.bitrate.*1e-9 ./log2(m) ) .* ir(i) .* 1e9; - baudrate = floor( uloops.bitrate.*1e-9 ./log2(m) ) .* 1e9; - baudrate_gather = union(baudrate_gather,baudrate); - baudrate_ticks = 100:20:240; - bitrate_ticks = 300:30:480; - - tp = TransmissionPerformance; - netRatesVNLE = tp.calculateNetRate(bitrate, 'NGMI', inf_rate_pam(i,:), 'BER', ber_vnle(i,:)); - - %%% NGMI - figure(12) - hold on - title(sprintf('Performance at 1310 nm')); - plot(baudrate.*1e-9,inf_rate_pam(i,:),'DisplayName',sprintf('NGMI; PAM %d',m),'Color',cols(i,:),'LineStyle','-'); - xlabel('Baud rate in GBd'); - ylabel('NGMI') - beautifyBERplot() - xticks(baudrate_ticks); - xlim([min(baudrate_ticks) max(baudrate_ticks)]); - - - %%% AIR - figure(14) - hold on - title(sprintf('Performance at 1310 nm')); - plot(baudrate.*1e-9,inf_rate_pam(i,:).*bitrate.*1e-9,'DisplayName',sprintf('AIR; PAM %d',m),'Color',cols(i,:),'LineStyle','-'); - xlabel('Baud rate in GBd'); - ylabel('AIR'); - beautifyBERplot() - xticks(baudrate_ticks); - xlim([min(baudrate_ticks) max(baudrate_ticks)]); - - %%% RATES - figure(16) - hold on - title(sprintf('Performance at 1310 nm')); - if i == 1 - hv = 'on'; - else - hv = 'off'; - end - plot(baudrate*1e-9,netRatesVNLE.SDHD.NetRate.*1e-9,'DisplayName',sprintf('SD+HD FEC',m),'Color',cols(i,:),'LineStyle','-','HandleVisibility',hv,'Marker','o'); - plot(baudrate.*1e-9,netRatesVNLE.HD.NetRate.*1e-9,'DisplayName',sprintf('HD FEC',m),'Color',cols(i,:),'LineStyle',':','HandleVisibility',hv,'Marker','diamond'); - plot(baudrate.*1e-9,netRatesVNLE.KP4_hamming.NetRate*1e-9,'DisplayName',sprintf('KP4+Hamming',m),'Color',cols(i,:),'LineStyle','-.','HandleVisibility',hv,'Marker','square'); - xticks(baudrate_ticks); - xlim([min(baudrate_ticks) max(baudrate_ticks)]); - xlabel('Baud rate in GBd'); - ylabel('Net Bitrate in Gbps') - beautifyBERplot() - ylim([250 410]) - - %%% CODE OVERHEAD IN % - figure(18) - hold on - title(sprintf('Performance at 1310 nm')); - plot(baudrate*1e-9,100.*(1-netRatesVNLE.SDHD.CodeRate)./netRatesVNLE.SDHD.CodeRate,'DisplayName',sprintf('SD+HD FEC',m),'Color',cols(i,:),'LineStyle','-','HandleVisibility',hv,'Marker','o'); - plot(baudrate.*1e-9,100.*(1-netRatesVNLE.HD.CodeRate)./netRatesVNLE.HD.CodeRate,'DisplayName',sprintf('HD FEC',m),'Color',cols(i,:),'LineStyle',':','HandleVisibility',hv,'Marker','diamond'); - plot(baudrate.*1e-9,100.*(1-netRatesVNLE.KP4_hamming.CodeRate)./netRatesVNLE.KP4_hamming.CodeRate,'DisplayName',sprintf('KP4+Hamming',m),'Color',cols(i,:),'LineStyle','-.','HandleVisibility',hv,'Marker','square'); - xticks(baudrate_ticks); - xlim([min(baudrate_ticks) max(baudrate_ticks)]); - xlabel('Baud rate in GBd'); - ylabel('FEC Overhead in %') - beautifyBERplot() - - %%% CLASSIC BER - figure(22) - subplot(1,4,i) - hold on - plot(baudrate*1e-9,ber_vnle(i,:),'DisplayName',sprintf('Tx precomp + VNLE',m),'Color',cols(i,:),'LineStyle','-','HandleVisibility','on','Marker','o'); - plot(baudrate*1e-9,ber_mlse(i,:),'DisplayName',sprintf('VNLE + PF + MLSE',m),'Color',cols(i,:),'LineStyle',':','HandleVisibility','on','Marker','square'); - plot(baudrate*1e-9,ber_db(i,:),'DisplayName',sprintf('Diff. Code + DB tgt.',m),'Color',cols(i,:),'LineStyle','--','HandleVisibility','on','Marker','diamond'); - yline(4.85e-3,'HandleVisibility','off'); - yline(2e-2,'HandleVisibility','off'); - xticks(baudrate*1e-9); - xlim([min(baudrate*1e-9) max(baudrate*1e-9)]); - ylim([1e-4 0.3]); - xlabel('Baudrate in GBd'); - if i == 1 - ylabel('BER') - end - beautifyBERplot() - set(gca, 'YScale', 'log'); -legend - subplot(1,4,4) - hold on - if m == 4 - - plot(bitrate*1e-9,ber_db(i,:),'DisplayName',sprintf('Diff. Code + DB tgt.'),'Color',cols(i,:),'LineStyle','--','HandleVisibility','on','Marker','diamond'); - - elseif m == 6 - - plot(bitrate*1e-9,ber_vnle(i,:),'DisplayName',sprintf('VNLE + PF + MLSE'),'Color',cols(i,:),'LineStyle','-','HandleVisibility','on','Marker','o'); - % plot(bitrate*1e-9,ber_mlse(i,:),'DisplayName',sprintf('MLSE',m),'Color',cols(i,:),'LineStyle',':','HandleVisibility',hv,'Marker','square'); - - elseif m ==8 - - plot(bitrate*1e-9,ber_vnle(i,:),'DisplayName',sprintf('Tx precomp + VNLE'),'Color',cols(i,:),'LineStyle','-','HandleVisibility','on','Marker','o'); - % plot(bitrate*1e-9,ber_mlse(i,:),'DisplayName',sprintf('MLSE',m),'Color',cols(i,:),'LineStyle',':','HandleVisibility',hv,'Marker','square'); - - end - yline(4.85e-3,'HandleVisibility','off'); - yline(2e-2,'HandleVisibility','off'); - xticks(bitrate_ticks); - xlim([min(bitrate_ticks) max(bitrate_ticks)]); - ylim([1e-4 0.3]); - xlabel('Gross Bitrate in Gbps'); - % ylabel('BER') - beautifyBERplot() - set(gca, 'YScale', 'log'); - - -end - - - -figure() -title(sprintf('%d km | 1310 nm | PAM %d | VNLE',uloops.link_length,uloops.M)); -hold on -line([min(uloops.bitrate.*1e-9) max(uloops.bitrate.*1e-9)],[min(uloops.bitrate.*1e-9) max(uloops.bitrate.*1e-9)],'Color',[.7,.7,.7],'Marker','none','Handlevisibility','off'); -plot(uloops.bitrate.*1e-9,cellfun(@min, inf_rate_vnle).*uloops.bitrate./log2(uloops.M).*1e-9,'DisplayName',sprintf('AIR'),'Color',cols(1,:),'LineStyle',':'); -plot(uloops.bitrate.*1e-9,netRatesVNLE.SDHD.NetRate.*1e-9,'DisplayName',sprintf('SD+HD'),'Color',cols(2,:),'LineStyle',':'); -plot(uloops.bitrate.*1e-9,netRatesVNLE.HD.NetRate.*1e-9,'DisplayName',sprintf('HD'),'Color',cols(3,:),'LineStyle',':'); -plot(uloops.bitrate.*1e-9,netRatesVNLE.KP4_hamming.NetRate.*1e-9,'DisplayName',sprintf('KP4+Hamming'),'Color',cols(4,:),'LineStyle',':'); -beautifyBERplot() -xlim([min(uloops.bitrate.*1e-9) max(uloops.bitrate.*1e-9)]) -xlabel('Gross Bitrate in Gbps'); -ylabel('Net Bitrate in Gbps'); -legend +% +% m = 6; +% ir = [2,2.5,3]; +% cols = linspecer(6); +% baudrate_gather = []; +% for i = 1:3 +% m = uloops.M(i); +% +% %%% GET VNLE VALS +% precode = 0; +% precomp = 1; +% a = wh_master.getStoValue('ber',precomp, precode, uloops.bitrate , uloops.laser_wavelength, m, uloops.link_length); +% ber_vnle(i,:) = cellfun(@(x) x.vnle_pf_package{1,1}.ber_vnle, a); +% +% %%% GET DB VALS +% precode = 1; +% precomp = 0; +% a = wh_master.getStoValue('ber',precomp, precode, uloops.bitrate , uloops.laser_wavelength, m, uloops.link_length); +% ber_db(i,:) = cellfun(@(x) x.dbtgt_package{1,1}.ber, a); +% +% %%% GET MLSE VALS +% precode = 0; +% precomp = 0; +% a = wh_master.getStoValue('ber',precomp, precode, uloops.bitrate , uloops.laser_wavelength, m, uloops.link_length); +% ber_mlse(i,:) = cellfun(@(x) x.vnle_pf_package{1,1}.ber_mlse, a); +% +% inf_rate_pam(i,:) = cellfun(@(x) x.vnle_pf_package{1,1}.air, a); +% inf_rate_pam(i,:) = inf_rate_pam(i,:)./log2(m); +% +% bitrate = floor( uloops.bitrate.*1e-9 ./log2(m) ) .* ir(i) .* 1e9; +% baudrate = floor( uloops.bitrate.*1e-9 ./log2(m) ) .* 1e9; +% baudrate_gather = union(baudrate_gather,baudrate); +% baudrate_ticks = 100:20:240; +% bitrate_ticks = 300:30:480; +% +% tp = TransmissionPerformance; +% netRatesVNLE = tp.calculateNetRate(bitrate, 'NGMI', inf_rate_pam(i,:), 'BER', ber_vnle(i,:)); +% +% %%% NGMI +% figure(12) +% hold on +% title(sprintf('Performance at 1310 nm')); +% plot(baudrate.*1e-9,inf_rate_pam(i,:),'DisplayName',sprintf('NGMI; PAM %d',m),'Color',cols(i,:),'LineStyle','-'); +% xlabel('Baud rate in GBd'); +% ylabel('NGMI') +% beautifyBERplot() +% xticks(baudrate_ticks); +% xlim([min(baudrate_ticks) max(baudrate_ticks)]); +% +% +% %%% AIR +% figure(14) +% hold on +% title(sprintf('Performance at 1310 nm')); +% plot(baudrate.*1e-9,inf_rate_pam(i,:).*bitrate.*1e-9,'DisplayName',sprintf('AIR; PAM %d',m),'Color',cols(i,:),'LineStyle','-'); +% xlabel('Baud rate in GBd'); +% ylabel('AIR'); +% beautifyBERplot() +% xticks(baudrate_ticks); +% xlim([min(baudrate_ticks) max(baudrate_ticks)]); +% +% %%% RATES +% figure(16) +% hold on +% title(sprintf('Performance at 1310 nm')); +% if i == 1 +% hv = 'on'; +% else +% hv = 'off'; +% end +% plot(baudrate*1e-9,netRatesVNLE.SDHD.NetRate.*1e-9,'DisplayName',sprintf('SD+HD FEC',m),'Color',cols(i,:),'LineStyle','-','HandleVisibility',hv,'Marker','o'); +% plot(baudrate.*1e-9,netRatesVNLE.HD.NetRate.*1e-9,'DisplayName',sprintf('HD FEC',m),'Color',cols(i,:),'LineStyle',':','HandleVisibility',hv,'Marker','diamond'); +% plot(baudrate.*1e-9,netRatesVNLE.KP4_hamming.NetRate*1e-9,'DisplayName',sprintf('KP4+Hamming',m),'Color',cols(i,:),'LineStyle','-.','HandleVisibility',hv,'Marker','square'); +% xticks(baudrate_ticks); +% xlim([min(baudrate_ticks) max(baudrate_ticks)]); +% xlabel('Baud rate in GBd'); +% ylabel('Net Bitrate in Gbps') +% beautifyBERplot() +% ylim([250 410]) +% +% %%% CODE OVERHEAD IN % +% figure(18) +% hold on +% title(sprintf('Performance at 1310 nm')); +% plot(baudrate*1e-9,100.*(1-netRatesVNLE.SDHD.CodeRate)./netRatesVNLE.SDHD.CodeRate,'DisplayName',sprintf('SD+HD FEC',m),'Color',cols(i,:),'LineStyle','-','HandleVisibility',hv,'Marker','o'); +% plot(baudrate.*1e-9,100.*(1-netRatesVNLE.HD.CodeRate)./netRatesVNLE.HD.CodeRate,'DisplayName',sprintf('HD FEC',m),'Color',cols(i,:),'LineStyle',':','HandleVisibility',hv,'Marker','diamond'); +% plot(baudrate.*1e-9,100.*(1-netRatesVNLE.KP4_hamming.CodeRate)./netRatesVNLE.KP4_hamming.CodeRate,'DisplayName',sprintf('KP4+Hamming',m),'Color',cols(i,:),'LineStyle','-.','HandleVisibility',hv,'Marker','square'); +% xticks(baudrate_ticks); +% xlim([min(baudrate_ticks) max(baudrate_ticks)]); +% xlabel('Baud rate in GBd'); +% ylabel('FEC Overhead in %') +% beautifyBERplot() +% +% %%% CLASSIC BER +% figure(22) +% subplot(1,4,i) +% hold on +% plot(baudrate*1e-9,ber_vnle(i,:),'DisplayName',sprintf('Tx precomp + VNLE',m),'Color',cols(i,:),'LineStyle','-','HandleVisibility','on','Marker','o'); +% plot(baudrate*1e-9,ber_mlse(i,:),'DisplayName',sprintf('VNLE + PF + MLSE',m),'Color',cols(i,:),'LineStyle',':','HandleVisibility','on','Marker','square'); +% plot(baudrate*1e-9,ber_db(i,:),'DisplayName',sprintf('Diff. Code + DB tgt.',m),'Color',cols(i,:),'LineStyle','--','HandleVisibility','on','Marker','diamond'); +% yline(4.85e-3,'HandleVisibility','off'); +% yline(2e-2,'HandleVisibility','off'); +% xticks(baudrate*1e-9); +% xlim([min(baudrate*1e-9) max(baudrate*1e-9)]); +% ylim([1e-4 0.3]); +% xlabel('Baudrate in GBd'); +% if i == 1 +% ylabel('BER') +% end +% beautifyBERplot() +% set(gca, 'YScale', 'log'); +% legend +% subplot(1,4,4) +% hold on +% if m == 4 +% +% plot(bitrate*1e-9,ber_db(i,:),'DisplayName',sprintf('Diff. Code + DB tgt.'),'Color',cols(i,:),'LineStyle','--','HandleVisibility','on','Marker','diamond'); +% +% elseif m == 6 +% +% plot(bitrate*1e-9,ber_vnle(i,:),'DisplayName',sprintf('VNLE + PF + MLSE'),'Color',cols(i,:),'LineStyle','-','HandleVisibility','on','Marker','o'); +% % plot(bitrate*1e-9,ber_mlse(i,:),'DisplayName',sprintf('MLSE',m),'Color',cols(i,:),'LineStyle',':','HandleVisibility',hv,'Marker','square'); +% +% elseif m ==8 +% +% plot(bitrate*1e-9,ber_vnle(i,:),'DisplayName',sprintf('Tx precomp + VNLE'),'Color',cols(i,:),'LineStyle','-','HandleVisibility','on','Marker','o'); +% % plot(bitrate*1e-9,ber_mlse(i,:),'DisplayName',sprintf('MLSE',m),'Color',cols(i,:),'LineStyle',':','HandleVisibility',hv,'Marker','square'); +% +% end +% yline(4.85e-3,'HandleVisibility','off'); +% yline(2e-2,'HandleVisibility','off'); +% xticks(bitrate_ticks); +% xlim([min(bitrate_ticks) max(bitrate_ticks)]); +% ylim([1e-4 0.3]); +% xlabel('Gross Bitrate in Gbps'); +% % ylabel('BER') +% beautifyBERplot() +% set(gca, 'YScale', 'log'); +% +% +% end +% +% +% +% figure() +% title(sprintf('%d km | 1310 nm | PAM %d | VNLE',uloops.link_length,uloops.M)); +% hold on +% line([min(uloops.bitrate.*1e-9) max(uloops.bitrate.*1e-9)],[min(uloops.bitrate.*1e-9) max(uloops.bitrate.*1e-9)],'Color',[.7,.7,.7],'Marker','none','Handlevisibility','off'); +% plot(uloops.bitrate.*1e-9,cellfun(@min, inf_rate_vnle).*uloops.bitrate./log2(uloops.M).*1e-9,'DisplayName',sprintf('AIR'),'Color',cols(1,:),'LineStyle',':'); +% plot(uloops.bitrate.*1e-9,netRatesVNLE.SDHD.NetRate.*1e-9,'DisplayName',sprintf('SD+HD'),'Color',cols(2,:),'LineStyle',':'); +% plot(uloops.bitrate.*1e-9,netRatesVNLE.HD.NetRate.*1e-9,'DisplayName',sprintf('HD'),'Color',cols(3,:),'LineStyle',':'); +% plot(uloops.bitrate.*1e-9,netRatesVNLE.KP4_hamming.NetRate.*1e-9,'DisplayName',sprintf('KP4+Hamming'),'Color',cols(4,:),'LineStyle',':'); +% beautifyBERplot() +% xlim([min(uloops.bitrate.*1e-9) max(uloops.bitrate.*1e-9)]) +% xlabel('Gross Bitrate in Gbps'); +% ylabel('Net Bitrate in Gbps'); +% legend % plot(uloops.bitrate.*1e-9,ber_vnle,'DisplayName',sprintf('NGMI MLSE; %d km',len),'Color',cols(1,:),'LineStyle','-'); diff --git a/projects/IMDD_base_system/imdd_model.m b/projects/IMDD_base_system/imdd_model.m index ac6bd0c..06b1812 100644 --- a/projects/IMDD_base_system/imdd_model.m +++ b/projects/IMDD_base_system/imdd_model.m @@ -19,19 +19,13 @@ fdac = 256e9; fadc = 256e9; random_key = 1; -interference_attenuation = 0; -is_mpi = 1; - -precomp = 0; -db_precode = 0; - -db_encode = 0; - rcalpha = 0.05; kover = 16; + vbias_rel = 0.5; -u_pi = 2.9; +u_pi = 3; vbias = -vbias_rel*u_pi; + laser_wavelength = 1293; laser_linewidth = 0; tx_bw_nyquist = 0.8; @@ -40,7 +34,7 @@ tx_bw_nyquist = 0.8; link_length = 1; % RX -rop = -5; +rop = -8; rx_bw_nyquist = 0.8; vnle_order1 = 50; @@ -68,7 +62,7 @@ mu_dfe = 0.0004; dfe_ = sum(dfe_order)>0; -doub_mode = db_mode.no_db; +duob_mode = db_mode.no_db; %%% change specific parameter if given in varargin % Parse optional input arguments @@ -90,43 +84,6 @@ if ~isempty(varargin) end end -if doub_mode ~= db_mode.db_encoded - if precomp == 0 && db_precode == 1 - doub_mode = db_mode.db_precoded; - - db_precode = 1; % preceded data (in my measurement set, this corresponds to low precomp too!) - discard_precode = 0; % - emulate_precode = 0; - legendentry = 'low precomp; precoded'; - disp('low precomp; precoded') - elseif precomp == 1 && db_precode == 1 - doub_mode = db_mode.db_emulate; - - db_precode = 0; % preceded data (in my measurement set, this corresponds to low precomp too!) - discard_precode = 0; % - emulate_precode = 1; - legendentry = 'high precomp; precoded'; - disp('high precomp; precoded') - elseif precomp == 0 && db_precode == 0 - doub_mode = db_mode.db_discard; - - db_precode = 1; % preceded data (in my measurement set, this corresponds to low precomp too!) - discard_precode = 1; % - emulate_precode = 0; - legendentry = 'no precomp; not precoded'; - disp('no precomp; not precoded') - elseif precomp == 1 && db_precode == 0 - doub_mode = db_mode.no_db; - - db_precode = 0; % preceded data (in my measurement set, this corresponds to low precomp too!) - discard_precode = 0; % - emulate_precode = 0; - legendentry = 'high precomp; not precoded'; - disp('high precomp; not precoded') - end -else - -end fsym_ = floor( bitrate*1e-9./log2(M) ).*1e9; @@ -137,28 +94,19 @@ end f_nyquist = fsym/2; - - -Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",16,"alpha",rcalpha); -rcalpha = 1; Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"alpha",rcalpha); -db_precode = 0; -db_encode = 0; -apply_pulsef = 1; [Digi_sig,Symbols,Tx_bits] = PAMsource(... "fsym",fsym,"M",M,"order",18,"useprbs",0,... "fs_out",fdac,... "applyclipping",0,"clipfactor",1.5,... "applypulseform",apply_pulsef,"pulseformer",Pform,... "randkey",random_key,... - "db_precode",db_precode,"db_encode",db_encode,... + 'duobinary_mode',duob_mode,... "mrds_code",0,"mrds_blocklength",512).process(); Digi_sig.spectrum("displayname",'Digi Spectrum','fignum',10,'normalizeTo0dB',1); - - %%%%% AWG % El_sig = M8199A("kover",kover).process(Digi_sig); El_sig = AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",0,"kover",kover,"bit_resolution",12,"upsampling_method","samplehold","precomp_sinc_rolloff",1).process(Digi_sig); @@ -172,13 +120,17 @@ El_sig = Filter('filtdegree',4,"f_cutoff",tx_bwl,"fs",fdac*kover,"filterType",fi % El_sig.spectrum("displayname",'Digi Spectrum','fignum',100,'normalizeTo0dB',1); %%%%% Electrical Driver Amplifier %%%%%% -El_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",3).process(El_sig); +% El_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",3).process(El_sig); El_sig = El_sig.normalize("mode","oneone"); +scaling = 0.6*(u_pi/2-abs(vbias-u_pi/2)); +El_sig = El_sig .* scaling; %%%%% MODULATE E/O CONVERSION %%%%%% -[Opt_sig] = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs,"lambda",laser_wavelength,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",random_key+1).process(El_sig); +[Opt_sig] = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs,"lambda",laser_wavelength,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",random_key+1,"alpha",alpha).process(El_sig); -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.spectrum("displayname",'Opt Spectrum','fignum',10,'normalizeTo0dB',1); + +Opt_sig = Fiber("fsimu",Opt_sig.fs,"fiber_length",link_length,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.07).process(Opt_sig); %%%%%% ROP %%%%%% Rx_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",rop).process(Opt_sig); @@ -202,110 +154,75 @@ Scpe_sig = Scope("fsimu",fdac*kover,"fadc",fadc,... "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(Rx_sig); -Scpe_cell{1} = Scpe_sig; - - -if db_precode - Symbols_precoded = Symbols; -end - output = struct(); -vnle_pf_package = {}; -vnle_dfe_package = {}; -dbtgt_package = {}; + +%%%%%% Sample to 2x fsym %%%%%% +Scpe_sig = Scpe_sig.resample("fs_out",2*fsym); +Scpe_sig.signal = Scpe_sig.signal(1:2*length(Symbols)); + +%%%%%% Sync Rx signal with reference %%%%%% +[Scpe_sig,~] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym,"debug_plots",1); + +Scpe_sig = Filter('filtdegree',4,"f_cutoff",Symbols.fs.*0.5,"fs",Scpe_sig.fs,"filterType",filtertypes.gaussian,"active",true).process(Scpe_sig); + +Scpe_sig = Scpe_sig - mean(Scpe_sig.signal); + +%%% EQUALIZING -proc_occ = min(1,length(Scpe_cell)); -for occ = 1%:proc_occ +% -------------------- 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); - Scpe_sig = Scpe_cell{occ}; +output.ffe_results = ffe(eq_ffe,M,Scpe_sig,Symbols,Tx_bits, ... + "precode_mode",duob_mode,'showAnalysis',0,"postFFE",[], ... + "eth_style_symbol_mapping",0); - %%%%%% Sample to 2x fsym %%%%%% - Scpe_sig = Scpe_sig.resample("fs_out",2*fsym); +output.ffe_results.metrics.print - %%%%%% Sync Rx signal with reference %%%%%% - [Scpe_sig,~] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym); +% -------------------- DFE -------------------- +eq_dfe = EQ("Ne",ffe_order,"Nb",[2,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); - Scpe_sig = Filter('filtdegree',4,"f_cutoff",Symbols.fs.*0.5,"fs",Scpe_sig.fs,"filterType",filtertypes.gaussian,"active",true).process(Scpe_sig); +output.dfe_results = ffe(eq_dfe,M,Scpe_sig,Symbols,Tx_bits, ... + "precode_mode",duob_mode,'showAnalysis',0,"postFFE",[], ... + "eth_style_symbol_mapping",0); - Scpe_sig = Scpe_sig - mean(Scpe_sig.signal); - % - % Pform = Pulseformer("fsym",Scpe_sig.fs,"fdac",2*fsym,"pulse","rrc","pulselength",16,"alpha",rcalpha,"matched",0); - % - % Scpe_sig_matched = Pform.process(Scpe_sig); - % - % Scpe_sig.spectrum("normalizeTo0dB",0,"fignum",336,"displayname","scope "); - % Scpe_sig_matched.spectrum("normalizeTo0dB",0,"fignum",336,"displayname","matched"); - - %%% EQUALIZING +output.dfe_results.metrics.print("description",'DFE'); - % eq_mlse = FFE_DCremoval("epochs_tr",5,"epochs_dd",5,"len_tr",len_tr,"mu_dd",mu_ffe(1),"mu_tr",0,"order",ffe_order(1),"sps",2,"decide",0,"dc_buffer_len",1,"mu_dc",0.05); - % eq_mlse = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",len_tr,"mu_dd",mu_ffe(1),"mu_tr",0,"order",ffe_order(1),"sps",2,"decide",0); - % eq_mlse = FFE_DCremoval("epochs_tr",5,"epochs_dd",5,"len_tr",len_tr,"mu_dd",mu_ffe(1),"mu_tr",0,"order",ffe_order(1),"sps",2,"decide",0,"dc_buffer_len",512,"mu_dc",0.05); +% -------------------- VNLE + MLSE -------------------- +pf_ncoeffs = 1; +ffe_order3 = [50, 5, 5]; +eq_v = EQ("Ne",ffe_order3,"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); - mu_ffe = [mu_ffe1 mu_ffe2 mu_ffe3]; - vnle_order=[vnle_order1,vnle_order2,vnle_order3]; +mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); - % %%%%% VNLE + DFE %%%% - if 0 - - eq_vnle_dfe = EQ("Ne",vnle_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_2 = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",2001,"sps",1,"decide",0); - - [result] = vnle(eq_vnle_dfe,M,Scpe_sig,Symbols,Tx_bits,"precode_mode",doub_mode,"showAnalysis",1,"postFFE",[]); - vnle_dfe_package{occ} = result; - - end - %%%%% VNLE + PF + MLSE %%%% - if 1 - - % len_tr = length(Symbols)-1000; - eq_vnle_ = EQ("Ne",vnle_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); - % eq_vnle_ = VNLE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",[0.0004 0.0005 0.0006],"mu_tr",0,"order",vnle_order,"sps",2,"decide",0); - pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); - mlse_ = MLSE_viterbi("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); - - [result] = vnle_postfilter_mlse(eq_vnle_,pf_,mlse_,M,Scpe_sig,Symbols,Tx_bits,"precode_mode",doub_mode,'showAnalysis',1); - vnle_pf_package{occ} = result; - - end +[output.vnle_results, output.mlse_results] = vnle_postfilter_mlse(eq_v, pf_, mlse_, M, Scpe_sig, Symbols, Tx_bits, ... + "precode_mode", duob_mode, 'showAnalysis', 0, "postFFE", [], "eth_style_symbol_mapping", 0); - %%%%% Duobinary Targeting %%%% - if 1 +% -------------------- DB target -------------------- +mlse_db_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",M,'trellis_states',PAMmapper(M,0).levels); +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); +output.dbt_results = duobinary_target(eq_,mlse_db_, M, Scpe_sig, Symbols, Tx_bits, ... + "precode_mode", duob_mode, 'showAnalysis', 0, "postFFE", []); - mlse_db = MLSE_viterbi("DIR",[1,1],"duobinary_output",0,"M",M,"trellis_states",PAMmapper(M,0).levels); - eq_db = EQ("Ne",vnle_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); - - [result] = duobinary_target(eq_db, mlse_db, M, Scpe_sig, Symbols, Tx_bits, "precode_mode", doub_mode,'showAnalysis',0); - dbtgt_package{occ} = result; +output.dbt_results.metrics.print("description",'Duobinary'); - end - - %%%%%% %db signaling => db encoded %%%%% - if 0 - mlse_db_enc = MLSE_viterbi("DIR",[1,1],"duobinary_output",0,"M",M,"trellis_states",PAMmapper(M,0).levels); - eq_db_enc = EQ("Ne",vnle_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); - [result] = duobinary_signaling(eq_db_enc, mlse_db_enc,M, Scpe_sig ,Symbols, Tx_bits); - dbenc_package{occ} = result; - end - - - % autoArrangeFigures; - disp('- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ') - fprintf('\n') - - -end - -output.vnle_dfe_package = vnle_dfe_package; -output.vnle_pf_package = vnle_pf_package; -output.dbtgt_package = dbtgt_package; - -if ~isempty(curFolder) - cd(curFolder); -end +disp('- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ') +fprintf('\n') end \ No newline at end of file diff --git a/projects/Messung_Zürich/submit_handle.m b/projects/IMDD_base_system/submit_handle.m similarity index 100% rename from projects/Messung_Zürich/submit_handle.m rename to projects/IMDD_base_system/submit_handle.m diff --git a/projects/WDM/WDM_auswertung.m b/projects/WDM/WDM_auswertung.m index a6a666c..4c213fa 100644 --- a/projects/WDM/WDM_auswertung.m +++ b/projects/WDM/WDM_auswertung.m @@ -1,325 +1,12 @@ -%% Robust eval/plot script for res.* (handles missing/partial dims) -% res.ffe stored like: output_ffe{ch, rop, realiz, eval_distance} -% BUT: some dimensions may be singleton or not fully filled (e.g. only 0 km, only 1 realization) +filename = "C:\Users\Silas\Documents\MATLAB\Datensätze\FWM_2025\WDM_20260106_131152_n261_14408161_10km_16ch_400ghz_co_alpha0_5.mat"; +filename = "C:\Users\Silas\Documents\MATLAB\Datensätze\FWM_2025\WDM_20260106_131152_n262_14408165_10km_16ch_400ghz_co_alpha-0_5.mat"; +filename = "C:\Users\Silas\Documents\MATLAB\Datensätze\FWM_2025\WDM_20260106_134921_n223_14409338_10km_16ch_400ghz_co_alpha0_0.mat"; +res = load(filename,'res'); +res = res.res; -% -------------------- Basic metadata -------------------- -rop = res.settings.rop(:); -wavelengthplan = res.settings.wavelengthplan(:); -distances = res.eval_dist_km(:); % can be [0], [2 10], etc. +%% +% Routine A: plot BER curves and compute crossings +S = plot_BER_vs_ROP(res, 'fec', 3.8e-3); -N_ch = numel(wavelengthplan); -N_rop = numel(rop); - -% -------------------- Determine dims robustly -------------------- -% Ensure we have 4D addressing even if MATLAB collapses trailing singletons -dims = size(res.ffe); -if numel(dims) < 4 - dims(end+1:4) = 1; -end -N_realiz = dims(3); -N_distances = dims(4); - -% If res.eval_dist_km length differs from dims(4), trust dims(4) -if numel(distances) ~= N_distances - % best-effort fallback - distances = (1:N_distances).'; -end - -% -------------------- Plot settings -------------------- -fec = 3.8e-3; % choose one; you can switch to 2.2e-4 if needed - -qLow = 0.00; -qHigh = 1.00; - -% Color per wavelength/channel -try - cols = ccbrewer2('Set1', N_ch); -catch - cols = linspecer(N_ch); -end - -% Containers for FEC crossings per channel (per eval distance) -Sffe = cell(N_distances, N_ch); -Sdfe = cell(N_distances, N_ch); -Svnle = cell(N_distances, N_ch); -Smlse = cell(N_distances, N_ch); -Sdbt = cell(N_distances, N_ch); - -% -------------------- Choose which eval distances to plot -------------------- -% If you want all: eval_list = 1:N_distances; -% If you want a specific distance (e.g., last): eval_list = N_distances; -eval_list = 1;%1:N_distances; - -% -------------------- BER vs ROP figures (one figure per eval distance) -------------------- -for eval_ptr = eval_list - - figure('Name', sprintf('BER vs ROP @ %s', distLabel(distances, eval_ptr))); - hold on; - - for ch = 1:N_ch - - % ---- Extract cell slices robustly: [N_rop x N_realiz] ---- - ffe_cells = sliceCells4D(res.ffe , ch, eval_ptr, N_rop); - dfe_cells = sliceCells4D(res.dfe , ch, eval_ptr, N_rop); - vnle_cells = sliceCells4D(res.vnle, ch, eval_ptr, N_rop); - mlse_cells = sliceCells4D(res.mlse, ch, eval_ptr, N_rop); - dbt_cells = sliceCells4D(res.dbt , ch, eval_ptr, N_rop); - - % ---- FEC crossings (works with incomplete columns) ---- - [Sffe{eval_ptr,ch}, ~] = fecCrossings(rop, ffe_cells, fec); - [Sdfe{eval_ptr,ch}, ~] = fecCrossings(rop, dfe_cells, fec); - [Svnle{eval_ptr,ch}, ~] = fecCrossings(rop, vnle_cells, fec); - [Smlse{eval_ptr,ch}, ~] = fecCrossings(rop, mlse_cells, fec); - [Sdbt{eval_ptr,ch}, ~] = fecCrossings(rop, dbt_cells, fec); - - % ---- Extract BER matrices using only complete realizations ---- - ffe_mat = extractCompleteBER(ffe_cells); % N_rop x K - dfe_mat = extractCompleteBER(dfe_cells); - vnle_mat = extractCompleteBER(vnle_cells); - mlse_mat = extractCompleteBER(mlse_cells); - dbt_mat = extractCompleteBER(dbt_cells); - - showLegend = (ch == 1); % one legend entry per technique - - % Plot bands + mean (only if non-empty) - if ~isempty(ffe_mat) - plotBandMeanBL(rop, ffe_mat, cols(ch,:), sprintf('FFE @ %dnm', round(wavelengthplan(ch))), qLow, qHigh, '--s', showLegend); - end - % if ~isempty(dfe_mat) - % plotBandMeanBL(rop, dfe_mat, cols(ch,:), sprintf('DFE @ %dnm', round(wavelengthplan(ch))), qLow, qHigh, '-o', showLegend); - % end - % if ~isempty(vnle_mat) - % plotBandMeanBL(rop, vnle_mat, cols(ch,:), sprintf('VNLE @ %dnm', round(wavelengthplan(ch))), qLow, qHigh, '--x', showLegend); - % end - % if ~isempty(mlse_mat) - % plotBandMeanBL(rop, mlse_mat, cols(ch,:), sprintf('VNLE+PF+MLSE @ %dnm', round(wavelengthplan(ch))), qLow, qHigh, '-o', showLegend); - % end - % if ~isempty(dbt_mat) - % plotBandMeanBL(rop, dbt_mat, cols(ch,:), sprintf('DBt+MLSE @ %dnm', round(wavelengthplan(ch))), qLow, qHigh, '--v', showLegend); - % end - - end - - set(gca,'XScale','linear','YScale','log','TickLabelInterpreter','latex','FontSize',11); - yline([3.8e-3, 2.2e-4], 'HandleVisibility','off','LineWidth',1.5); - ylabel('BER'); - xlabel('ROP [dB]'); - title(sprintf('BER vs. ROP @ %s', distLabel(distances, eval_ptr))); - xlim([min(rop) max(rop)]); - ylim([1e-5 0.3]); - grid on; - legend show; -end - - -%% -------------------- VIOLIN of ROP at FEC crossing -------------------- -% Build technique groups; each group is {N_distances x N_ch} cell entries -% You can pick a distance to show (e.g., eval_ptr=1 for 0 km, or last) -eval_ptr_violin = min(N_distances, max(1, N_distances)); % default: last available -% eval_ptr_violin = 1; - -% Choose which technique(s) to show -techNames = {'VNLE'}; -S_groups = {Svnle(eval_ptr_violin,:)}; - -% Only plot if violinplot exists -if exist('violinplot','file') == 2 - figure('Name', sprintf('FEC crossing violin @ %s', distLabel(distances, eval_ptr_violin))); hold on; - - colsTech = linspecer(numel(S_groups)); - - for i = 1:numel(S_groups) - S_cell = S_groups{i}; % 1 x N_ch cell, each cell is 1xK crossings - - % Pad to rectangular: rows = realizations (max K), cols = wavelengths (N_ch) - Kmax = max(cellfun(@numel, S_cell)); - if isempty(Kmax) || Kmax == 0 - continue; - end - S_mat = zeros(Kmax, N_ch); - for ch = 1:N_ch - k = numel(S_cell{ch}); - if k > 0 - S_mat(1:k, ch) = S_cell{ch}(:); - end - end - - catLabels = arrayfun(@(nm) sprintf('%d nm', nm), round(wavelengthplan), 'UniformOutput', false); - violinplot(S_mat, catLabels, ... - 'ViolinColor', colsTech(i,:), ... - 'ViolinAlpha', 0.10, ... - 'MarkerSize', 20, ... - 'ShowMedian', true, ... - 'EdgeColor', colsTech(i,:), ... - 'ShowWhiskers', false, ... - 'ShowData', true, ... - 'ShowBox', false, ... - 'Bandwidth', 0.05); - - % Add a legend proxy - plot(nan, nan, 'o', 'Color', colsTech(i,:), 'DisplayName', techNames{i}); - end - - ylim([-10,-6]) - ylabel('ROP at FEC crossing'); - title(sprintf('ROP to cross BER %.2e @ %s', fec, distLabel(distances, eval_ptr_violin))); - grid on; box on; - legend show; -else - fprintf('violinplot.m not found on path -> skipping violin plot.\n'); -end - - -%% ================= helpers ================= - -function cells2D = sliceCells4D(C4, ch, eval_ptr, N_rop) -% Returns a [N_rop x N_realiz] cell array (with N_realiz inferred) -% Works even if trailing dimensions are singleton. -% -% If C4 is smaller than expected, missing entries are returned as []. - -dims = size(C4); -if numel(dims) < 4, dims(end+1:4) = 1; end -N_realiz = dims(3); - -cells2D = cell(N_rop, N_realiz); - -% bounds (in case stored arrays are smaller) -chMax = dims(1); -ropMax = dims(2); -rMax = dims(3); -eMax = dims(4); - -if ch > chMax || eval_ptr > eMax - return; % all empty -end - -ropUse = min(N_rop, ropMax); -rUse = min(N_realiz, rMax); - -% Extract and place into a consistent sized cell matrix -tmp = squeeze(C4(ch, 1:ropUse, 1:rUse, eval_ptr)); -% squeeze may return vector/empty if rUse==1 etc. Normalize: -tmp = reshape(tmp, ropUse, rUse); - -cells2D(1:ropUse, 1:rUse) = tmp; -end - -function lbl = distLabel(distances, eval_ptr) -if eval_ptr <= numel(distances) - d = distances(eval_ptr); - if isfinite(d) - lbl = sprintf('%.0f km', d); - return; - end -end -lbl = sprintf('eval\\_%d', eval_ptr); -end - -function plotBandMeanBL(x, Y, color, techLabel, qLow, qHigh, lineSpec, showLegend) -% Y: (nPoints x nRealizations) -if isempty(Y), return; end - -% Remove realizations that are entirely zero -badCols = all(Y == 0, 1); -Y(:, badCols) = []; - -if isempty(Y), return; end - -% Avoid log(0) -Y(Y==0) = 1e-8; - -mu = mean(Y, 2, 'omitnan'); -lo = quantile(Y, qLow, 2); -hi = quantile(Y, qHigh, 2); -b = [mu - lo, hi - mu]; - -% boundedline optional; fallback to simple plot if missing -if 0%exist('boundedline','file') == 2 - [hl, hp] = boundedline(x(:), mu(:), b, lineSpec, 'alpha', 'transparency', 0.18); - set(hl, 'Color', color, 'LineWidth', 1.4, 'MarkerSize', 4); - set(hp, 'FaceColor', color, 'HandleVisibility','off'); -else - hl = plot(x(:), mu(:), lineSpec, 'LineWidth', 1.4, 'MarkerSize', 4); - set(hl, 'Color', color); -end - -if showLegend - set(hl, 'DisplayName', techLabel); -else - set(hl, 'HandleVisibility','off'); -end - -if 0%exist('outlinebounds','file') == 2 && exist('boundedline','file') == 2 - ho = outlinebounds(hl, hp); - set(ho, 'linestyle', ':', 'color', color, 'linewidth', 1, 'HandleVisibility','off'); -end -end - -function [S, noCrossingMask] = fecCrossings(rop, cellsNxR, fec) -% cellsNxR: N_rop x R cell array; each cell is struct with .metrics.BER -Y = extractCompleteBER(cellsNxR); % -> N_rop x K -if isempty(Y) - S = []; - noCrossingMask = []; - return; -end - -ok = mean(Y,1,'omitnan') <= 0.1; -Y = Y(:, ok); -if isempty(Y) - S = []; - noCrossingMask = []; - return; -end - -nR = size(Y,2); -S = nan(1,nR); -noCrossingMask = true(1,nR); - -rop = rop(:); -for j = 1:nR - y = Y(:,j); - above = (y > fec); - idx = find(above(1:end-1) & ~above(2:end), 1, 'first'); - if ~isempty(idx) - x1 = rop(idx); y1 = y(idx); - x2 = rop(idx+1); y2 = y(idx+1); - if isfinite(y1) && isfinite(y2) && y2 ~= y1 - t = (fec - y1) / (y2 - y1); - S(j) = x1 + t*(x2 - x1); - noCrossingMask(j) = false; - end - end -end -end - -function Y = extractCompleteBER(cellSlice) -% Keep only those realization columns where ALL N_rop entries are non-empty -if isempty(cellSlice), Y = []; return; end -nR = size(cellSlice,2); -keep = false(1,nR); - -for r = 1:nR - col = cellSlice(:,r); - keep(r) = all(cellfun(@(c) ~isempty(c), col)); -end - -if ~any(keep), Y = []; return; end - -Y = cellfun(@(c) c.metrics.BER, cellSlice(:,keep), 'UniformOutput', true); -end - -function Y = extractCompleteAlphas(cellSlice) -if isempty(cellSlice), Y = []; return; end -nR = size(cellSlice,2); -keep = false(1,nR); - -for r = 1:nR - col = cellSlice(:,r); - keep(r) = all(cellfun(@(c) ~isempty(c), col)); -end - -if ~any(keep), Y = []; return; end -Y = cellfun(@(c) c.metrics.Alpha, cellSlice(:,keep), 'UniformOutput', true); -end +%% Routine B: violin plot (independent) +plot_FEC_violin(res, 'tech','VNLE', 'fec',3.8e-3, 'eval_ptr',3, 'ylim',[-10 -3]); diff --git a/projects/WDM/WDM_model_10km_queue.m b/projects/WDM/WDM_model_10km_queue.m index 24b7e4e..e81cbb2 100644 --- a/projects/WDM/WDM_model_10km_queue.m +++ b/projects/WDM/WDM_model_10km_queue.m @@ -12,6 +12,7 @@ arguments options.rand_key = 1; options.num_realiz = 1; options.fwm_mitigation_technique = "co"; + options.chirpalpha = 0; end %% Add the imdd_simulation framework to the path @@ -59,7 +60,9 @@ else end if ~exist(output_root,'dir'), mkdir(output_root); end -fname = sprintf('WDM_%s_%s_%s_%dkm_%dch_%dghz_%s.mat', char(t), host, jobid, options.fiber_length_km(end), options.num_channels, options.channel_spacing.*1e-9, options.fwm_mitigation_technique); +fname = sprintf('WDM_%s_%s_%s_%dkm_%dch_%dghz_%s_alpha%0.1f', char(t), host, jobid, options.fiber_length_km(end), options.num_channels, options.channel_spacing.*1e-9, options.fwm_mitigation_technique, options.chirpalpha); +fname = strrep(fname,'.','_'); +fname = [fname,'.mat']; % Worker memory logs directory memlog_dir = fullfile(output_root, 'memlogs'); @@ -116,12 +119,12 @@ upsample_pow = 2^nextpow2(upsample_required); s.f_opt = fdac*kover*upsample_pow; s.f_opt_nyq = s.f_opt/2; -s.rop = -12:0.75:-0.75; -% s.rop = -12:1:0; +s.rop = -10:0.5:0; + %% ---------- Intermediate evaluation points (distance dimension) ---------- segment_length = 1; % km (must match fiber loop below) -eval_dist_km = [0]; +eval_dist_km = [2,4,6,8,10]; eval_dist_km = eval_dist_km(eval_dist_km <= link_length); if link_length > 0 @@ -208,7 +211,9 @@ for realiz = 1:s.num_realiz % 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,"alpha",0.8).process(El_sig); + "linewidth",laser_linewidth,"randomkey",s.random_key+l+realiz,"alpha",options.chirpalpha).process(El_sig); + + s.alpha = options.chirpalpha; % El_sig not needed after EML clear El_sig @@ -244,7 +249,7 @@ for realiz = 1:s.num_realiz nSegments = round(nSegments); zdw = 1310; - randomize_D = true; + randomize_D = false; % Guard for 0 km: avoid calling getDispersionVector(0,...) if it doesn't support it if nSegments > 0 @@ -256,7 +261,6 @@ for realiz = 1:s.num_realiz eval_ptr = 1; % =================== Queue throttle (prevents OOM) =================== % Limit how many futures can be outstanding (running + queued + finished-not-yet-fetched). - % Start conservatively; you can raise to 4*NumWorkers if stable. maxInFlight = max(2, p.NumWorkers); % ===================================================================== @@ -300,6 +304,10 @@ for realiz = 1:s.num_realiz end % ----------------------------------------------------- + save_results_per_realization( ... + s, eval_dist_km, ... + output_ffe, output_dfe, output_vnle, output_mlse, output_dbt, ... + output_root, fname); for seg = 1:nSegments @@ -318,7 +326,6 @@ for realiz = 1:s.num_realiz cnt = 0; total = length(s.rop)*N; - fprintf('total of %d jobs to enqueue at 0 km\n', total); for ri = 1:length(s.rop) for l = 1:N @@ -338,8 +345,8 @@ for realiz = 1:s.num_realiz len_tr, mu_dc, mu_ffe, mu_dfe, dfe_order, duob_mode, ... memlog_dir); - fprintf('Enqueued job %d/%d for realiz %d/%d, l=%d/%d, ri=%d/%d at 0 km (inflight=%d/%d)\n', ... - cnt, total, realiz, s.num_realiz, l, N, ri, length(s.rop), numel(F), maxInFlight); + fprintf('Seg: %d - Enqueued job %d/%d for realiz %d/%d, l=%d/%d, ri=%d/%d (inflight=%d/%d)\n', ... + seg, cnt, total, realiz, s.num_realiz, l, N, ri, length(s.rop), numel(F), maxInFlight); meta(end+1) = struct('l',l,'ri',ri,'realiz',realiz,'eval_ptr',eval_ptr); end @@ -356,24 +363,38 @@ for realiz = 1:s.num_realiz end - while ~isempty(F) + drainIterMax = 10; % 10 * 60s = 10 minutes + for drainIter = 1:drainIterMax + + if isempty(F) + break; + end + + % Status line: how many futures are in which states? + st = {F.State}; + nUnavail = sum(strcmp(st,'unavailable')); + nFinished = sum(strcmp(st,'finished')); + nRunning = sum(strcmp(st,'running')); + nQueued = sum(strcmp(st,'queued')); + nFailed = sum(strcmp(st,'failed')); + fprintf('[drain %d/%d] F=%d | finished=%d running=%d queued=%d failed=%d unavailable=%d\n', ... + drainIter, drainIterMax, numel(F), nFinished, nRunning, nQueued, nFailed, nUnavail); + + % Try to fetch at least one result (blocking) [F, meta, output_ffe, output_dfe, output_vnle, output_mlse, output_dbt] = ... collect_done(F, meta, output_ffe, output_dfe, output_vnle, output_mlse, output_dbt, true); + + % If still not empty, wait a bit before the next drain attempt + if ~isempty(F) + pause(60); + end end - %% Save results (per realization) - res = struct(); - res.settings = s; - res.eval_dist_km = eval_dist_km; - res.ffe = output_ffe; - res.dfe = output_dfe; - res.vnle = output_vnle; - res.mlse = output_mlse; - res.dbt = output_dbt; - 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')); + save_results_per_realization( ... + s, eval_dist_km, ... + output_ffe, output_dfe, output_vnle, output_mlse, output_dbt, ... + output_root, fname); % Per-realization large arrays that are no longer needed clear Opt_sig_wdm_fib Symbols Tx_bits Dvec @@ -524,6 +545,16 @@ end first = true; while ~isempty(F) + + idxUn = strcmp({F.State}, 'unavailable'); + if any(idxUn) + warning('collect_done: dropping %d unavailable futures.', sum(idxUn)); + F(idxUn) = []; + meta(idxUn) = []; + if isempty(F), break; end + end + + if first timeout = timeout_first; first = false; @@ -555,6 +586,31 @@ end end +function save_results_per_realization( ... + s, eval_dist_km, ... + output_ffe, output_dfe, output_vnle, output_mlse, output_dbt, ... + output_root, fname) + +% Assemble result struct +res = struct(); +res.settings = s; +res.eval_dist_km = eval_dist_km; +res.ffe = output_ffe; +res.dfe = output_dfe; +res.vnle = output_vnle; +res.mlse = output_mlse; +res.dbt = output_dbt; + +% Save (HDF5-based for large data) +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 + + + %% ========================= Memory logging helpers ========================= function logfile = make_worker_logfile(memlog_dir) diff --git a/projects/WDM/plot_BER_vs_ROP.m b/projects/WDM/plot_BER_vs_ROP.m new file mode 100644 index 0000000..b8ac31f --- /dev/null +++ b/projects/WDM/plot_BER_vs_ROP.m @@ -0,0 +1,199 @@ +%% ========================= +% Routine A: BER vs ROP plots +%% ========================= +function S = plot_BER_vs_ROP(res, varargin) +% S = plot_BER_vs_ROP(res, 'fec', 3.8e-3, 'eval_list', [], 'colormap', 'Spectral') +% +% Creates one BER-vs-ROP figure per evaluated distance. +% Returns S struct with FEC crossings: S.Sffe, S.Sdfe, S.Svnle, S.Smlse, S.Sdbt +% +% res.* assumed: res.ffe{ch,rop,realiz,eval_distance} etc. + +p = inputParser; +p.addParameter('fec', 3.8e-3, @(x)isnumeric(x)&&isscalar(x)); +p.addParameter('eval_list', [], @(x)isnumeric(x)); +p.addParameter('colormap', 'Spectral', @(x)ischar(x) || isstring(x)); +p.parse(varargin{:}); +fec = p.Results.fec; + +% -------------------- Basic metadata -------------------- +rop = res.settings.rop(:); +wavelengthplan = res.settings.wavelengthplan(:); +distances = res.eval_dist_km(:); + +N_ch = numel(wavelengthplan); +N_rop = numel(rop); + +% -------------------- Determine dims robustly -------------------- +dims = size(res.ffe); +if numel(dims) < 4, dims(end+1:4) = 1; end +N_distances = dims(4); + +if numel(distances) ~= N_distances + distances = (1:N_distances).'; +end + +% -------------------- Choose eval distances -------------------- +eval_list = p.Results.eval_list; +if isempty(eval_list) + eval_list = 1:N_distances; +end + +% -------------------- Colors -------------------- +try + cols = cbrewer2(char(p.Results.colormap), N_ch); +catch + cols = linspecer(N_ch); +end + +% -------------------- Crossings containers -------------------- +S = struct(); +S.rop = rop; +S.wavelengthplan = wavelengthplan; +S.distances = distances; +S.fec = fec; + +S.Sffe = cell(N_distances, N_ch); +S.Sdfe = cell(N_distances, N_ch); +S.Svnle = cell(N_distances, N_ch); +S.Smlse = cell(N_distances, N_ch); +S.Sdbt = cell(N_distances, N_ch); + +% -------------------- Plot per distance -------------------- +for eval_ptr = eval_list + + figure('Name', sprintf('BER vs ROP @ %s', distLabel(distances, eval_ptr))); + hold on; + + for ch = 1:N_ch + + % Extract cell slices + ffe_cells = sliceCells4D(res.ffe , ch, eval_ptr, N_rop); + dfe_cells = sliceCells4D(res.dfe , ch, eval_ptr, N_rop); + vnle_cells = sliceCells4D(res.vnle, ch, eval_ptr, N_rop); + mlse_cells = sliceCells4D(res.mlse, ch, eval_ptr, N_rop); + dbt_cells = sliceCells4D(res.dbt , ch, eval_ptr, N_rop); + + % Crossings (only depends on cells) + [S.Sffe{eval_ptr,ch}, ~] = fecCrossings(rop, ffe_cells, fec); + [S.Sdfe{eval_ptr,ch}, ~] = fecCrossings(rop, dfe_cells, fec); + [S.Svnle{eval_ptr,ch}, ~] = fecCrossings(rop, vnle_cells, fec); + [S.Smlse{eval_ptr,ch}, ~] = fecCrossings(rop, mlse_cells, fec); + [S.Sdbt{eval_ptr,ch}, ~] = fecCrossings(rop, dbt_cells, fec); + + % BER matrices (complete realizations only) + ffe_mat = extractCompleteBER(ffe_cells); % N_rop x K + + if ~isempty(ffe_mat) + % Your current style: plot each realization curve (no boundedline) + plot(rop, ffe_mat, 'Color', cols(ch,:), ... + 'DisplayName', sprintf('FFE @ %dnm', round(wavelengthplan(ch)))); + end + end + + set(gca,'XScale','linear','YScale','log', ... + 'TickLabelInterpreter','latex','FontSize',11); + + yline([3.8e-3, 2.2e-4], 'HandleVisibility','off','LineWidth',1.5); + + ylabel('BER'); + xlabel('ROP [dB]'); + title(sprintf('BER vs. ROP @ %s', distLabel(distances, eval_ptr))); + xlim([min(rop) max(rop)]); + ylim([1e-5 0.3]); + grid on; + legend show; + beautifyBERplot; +end +end + + + +%% ========================= +% Shared helpers (unchanged) +%% ========================= + +function cells2D = sliceCells4D(C4, ch, eval_ptr, N_rop) +dims = size(C4); +if numel(dims) < 4, dims(end+1:4) = 1; end +N_realiz = dims(3); + +cells2D = cell(N_rop, N_realiz); + +chMax = dims(1); +ropMax = dims(2); +rMax = dims(3); +eMax = dims(4); + +if ch > chMax || eval_ptr > eMax + return; +end + +ropUse = min(N_rop, ropMax); +rUse = min(N_realiz, rMax); + +tmp = squeeze(C4(ch, 1:ropUse, 1:rUse, eval_ptr)); +tmp = reshape(tmp, ropUse, rUse); + +cells2D(1:ropUse, 1:rUse) = tmp; +end + +function lbl = distLabel(distances, eval_ptr) +if eval_ptr <= numel(distances) + d = distances(eval_ptr); + if isfinite(d) + lbl = sprintf('%.0f km', d); + return; + end +end +lbl = sprintf('eval\\_%d', eval_ptr); +end + +function [S, noCrossingMask] = fecCrossings(rop, cellsNxR, fec) +Y = extractCompleteBER(cellsNxR); +if isempty(Y) + S = []; + noCrossingMask = []; + return; +end + +ok = mean(Y,1,'omitnan') <= 0.1; +Y = Y(:, ok); +if isempty(Y) + S = []; + noCrossingMask = []; + return; +end + +nR = size(Y,2); +S = nan(1,nR); +noCrossingMask = true(1,nR); + +rop = rop(:); +for j = 1:nR + y = Y(:,j); + above = (y > fec); + idx = find(above(1:end-1) & ~above(2:end), 1, 'first'); + if ~isempty(idx) + x1 = rop(idx); y1 = y(idx); + x2 = rop(idx+1); y2 = y(idx+1); + if isfinite(y1) && isfinite(y2) && y2 ~= y1 + t = (fec - y1) / (y2 - y1); + S(j) = x1 + t*(x2 - x1); + noCrossingMask(j) = false; + end + end +end +end + +function Y = extractCompleteBER(cellSlice) +if isempty(cellSlice), Y = []; return; end +nR = size(cellSlice,2); +keep = false(1,nR); +for r = 1:nR + col = cellSlice(:,r); + keep(r) = all(cellfun(@(c) ~isempty(c), col)); +end +if ~any(keep), Y = []; return; end +Y = cellfun(@(c) c.metrics.BER, cellSlice(:,keep), 'UniformOutput', true); +end diff --git a/projects/WDM/plot_FEC_violin.m b/projects/WDM/plot_FEC_violin.m new file mode 100644 index 0000000..5ce9217 --- /dev/null +++ b/projects/WDM/plot_FEC_violin.m @@ -0,0 +1,213 @@ +function plot_FEC_violin(res, varargin) +% plot_fec_violin_from_res Self-contained violin plot (no dependency on other routines) +% +% Usage: +% plot_fec_violin_from_res('WDM_20260106_....mat', 'tech', 'FFE', 'fec', 3.8e-3, 'eval_ptr', 2, 'ylim', [-10 -6]); +% plot_fec_violin_from_res(res, 'tech', 'FFE', 'fec', 3.8e-3, 'eval_ptr', 2); +% +% Notes: +% - Works with partially-filled res.* cell arrays (missing eval distances / missing realizations). +% - Uses only complete realizations per channel (all ROP points present). +% - Requires violinplot.m on path. + +p = inputParser; +p.addParameter('tech', 'FFE', @(x)ischar(x) || isstring(x)); +p.addParameter('fec', 3.8e-3, @(x)isnumeric(x)&&isscalar(x)); +p.addParameter('eval_ptr', [], @(x)isnumeric(x)&&isscalar(x)); +p.addParameter('ylim', [], @(x)isnumeric(x)&&numel(x)==2); +p.addParameter('title_prefix', 'ROP at FEC crossing', @(x)ischar(x)||isstring(x)); +p.parse(varargin{:}); + +if exist('violinplot','file') ~= 2 + error('violinplot.m not found on path. Add it to path first.'); +end + + +% ---------- Metadata ---------- +rop = res.settings.rop(:); +wavelengthplan = res.settings.wavelengthplan(:); +distances = res.eval_dist_km(:); + +N_ch = numel(wavelengthplan); +N_rop = numel(rop); + +% Determine distances dim robustly from the stored cell arrays (prefer ffe) +dims = size(res.ffe); +if numel(dims) < 4, dims(end+1:4) = 1; end +N_distances = dims(4); +if numel(distances) ~= N_distances + distances = (1:N_distances).'; +end + +% Choose eval distance +eval_ptr = p.Results.eval_ptr; +if isempty(eval_ptr) + eval_ptr = N_distances; % default: last available +end +if eval_ptr < 1 || eval_ptr > N_distances + error('eval_ptr=%d out of range. Available: 1..%d', eval_ptr, N_distances); +end + +% Choose technique field +tech = upper(string(p.Results.tech)); +switch tech + case "FFE", C4 = res.ffe; + case "DFE", C4 = res.dfe; + case "VNLE", C4 = res.vnle; + case "MLSE", C4 = res.mlse; + case "DBT", C4 = res.dbt; + otherwise + error('Unknown tech "%s". Use one of: FFE, DFE, VNLE, MLSE, DBT.', tech); +end + +fec = p.Results.fec; + +% ---------- Compute crossings per channel ---------- + +% ---------- Plot ---------- +figure('Name', sprintf('%s violin @ %s', tech, distLabel(distances, eval_ptr))); +hold on; +c = linspecer(N_distances); +for eval_ptr = 1:N_distances + S_cell = cell(1, N_ch); + + for ch = 1:N_ch + cells = sliceCells4D(C4, ch, eval_ptr, N_rop); % N_rop x N_realiz cell + S_cell{ch} = fecCrossings_only(rop, cells, fec); % 1 x K crossings (may be empty) + end + + % Pad to matrix for violinplot: rows=realizations (max K), cols=channels + Kmax = max(cellfun(@numel, S_cell)); + if isempty(Kmax) || Kmax == 0 + warning('No FEC crossings found for %s at eval_ptr=%d (%s).', tech, eval_ptr, distLabel(distances, eval_ptr)); + return; + end + + S_mat = NaN(Kmax, N_ch); + for ch = 1:N_ch + k = numel(S_cell{ch}); + if k > 0 + S_mat(1:k, ch) = S_cell{ch}(:); + end + end + + catLabels = arrayfun(@(nm) sprintf('%d nm', nm), round(wavelengthplan), 'UniformOutput', false); + + if numel(S_mat)>numel(wavelengthplan) + S_mat(isnan(S_mat)) = 0; + v=violinplot(S_mat, catLabels, ... + 'ViolinColor', c(eval_ptr,:), ... + 'ViolinAlpha', 0.10, ... + 'MarkerSize', 20, ... + 'ShowMedian', true, ... + 'EdgeColor', c(eval_ptr,:), ... + 'ShowWhiskers', false, ... + 'ShowData', true, ... + 'ShowBox', false, ... + 'Bandwidth', 0.05); + else + channs = 1:numel(wavelengthplan); + scatter( wavelengthplan ,S_mat,10,c(eval_ptr,:),'Marker','o','LineWidth',2,'DisplayName',distLabel(distances, eval_ptr)); + + channs = 1:numel(wavelengthplan); + scatter( round(wavelengthplan(S_mat==0)),S_mat(S_mat==0),10,c(eval_ptr,:),'Marker','o','HandleVisibility','off'); + end + +end + +ylabel('ROP at FEC crossing'); +title(sprintf('%s: %s | BER %.2e | %s', ... + p.Results.title_prefix, tech, fec, distLabel(distances, eval_ptr))); +grid on; box on; + +if ~isempty(p.Results.ylim) + ylim(p.Results.ylim); +end + +end + + +%% ================= helpers (self-contained) ================= + +function cells2D = sliceCells4D(C4, ch, eval_ptr, N_rop) +dims = size(C4); +if numel(dims) < 4, dims(end+1:4) = 1; end +N_realiz = dims(3); + +cells2D = cell(N_rop, N_realiz); + +chMax = dims(1); +ropMax = dims(2); +rMax = dims(3); +eMax = dims(4); + +if ch > chMax || eval_ptr > eMax + return; +end + +ropUse = min(N_rop, ropMax); +rUse = min(N_realiz, rMax); + +tmp = squeeze(C4(ch, 1:ropUse, 1:rUse, eval_ptr)); +tmp = reshape(tmp, ropUse, rUse); + +cells2D(1:ropUse, 1:rUse) = tmp; +end + +function lbl = distLabel(distances, eval_ptr) +if eval_ptr <= numel(distances) + d = distances(eval_ptr); + if isfinite(d) + lbl = sprintf('%.0f km', d); + return; + end +end +lbl = sprintf('eval\\_%d', eval_ptr); +end + +function S = fecCrossings_only(rop, cellsNxR, fec) +% Returns 1xK crossing positions for all complete realizations (drops bad ones). +Y = extractCompleteBER(cellsNxR); % N_rop x K +if isempty(Y), S = []; return; end + +% Optional quality gate (same as your previous scripts) +ok = mean(Y,1,'omitnan') <= 0.1; +Y = Y(:, ok); +if isempty(Y), S = []; return; end + +nR = size(Y,2); +S = nan(1,nR); + +rop = rop(:); +for j = 1:nR + y = Y(:,j); + above = (y > fec); + idx = find(above(1:end-1) & ~above(2:end), 1, 'first'); + if ~isempty(idx) + x1 = rop(idx); y1 = y(idx); + x2 = rop(idx+1); y2 = y(idx+1); + if isfinite(y1) && isfinite(y2) && y2 ~= y1 + t = (fec - y1) / (y2 - y1); + S(j) = x1 + t*(x2 - x1); + end + end +end + +S = S(isfinite(S)); % keep only successful crossings +end + +function Y = extractCompleteBER(cellSlice) +% Keep only those realization columns where ALL N_rop entries are non-empty +if isempty(cellSlice), Y = []; return; end +nR = size(cellSlice,2); +keep = false(1,nR); + +for r = 1:nR + col = cellSlice(:,r); + keep(r) = all(cellfun(@(c) ~isempty(c), col)); +end + +if ~any(keep), Y = []; return; end + +Y = cellfun(@(c) c.metrics.BER, cellSlice(:,keep), 'UniformOutput', true); +end