From 00f1c557c0e8b1d4b7441f21c8fa616be612b13d Mon Sep 17 00:00:00 2001 From: Silas Labor Zizou Date: Fri, 11 Apr 2025 15:27:32 +0200 Subject: [PATCH] MPI measurements --- projects/ECOC_2025/b2b_measurement.m | 0 projects/ECOC_2025/dsp_run_id.m | 235 ++++++++++++ projects/ECOC_2025/dsp_standalone.m | 10 + projects/ECOC_2025/measure_run.m | 207 +++++++++++ projects/ECOC_2025/measure_run_script.m | 172 +++++++++ .../plots_from_database/bias_vs_ber.m | 269 ++++++++++++++ .../plot_bias_awg_vs_ber.m | 271 ++++++++++++++ .../plots_from_database/plot_mpi_trial.m | 342 ++++++++++++++++++ projects/ECOC_2025/submit_dsp.m | 54 +++ 9 files changed, 1560 insertions(+) delete mode 100644 projects/ECOC_2025/b2b_measurement.m create mode 100644 projects/ECOC_2025/dsp_run_id.m create mode 100644 projects/ECOC_2025/dsp_standalone.m create mode 100644 projects/ECOC_2025/measure_run.m create mode 100644 projects/ECOC_2025/measure_run_script.m create mode 100644 projects/ECOC_2025/plots_from_database/bias_vs_ber.m create mode 100644 projects/ECOC_2025/plots_from_database/plot_bias_awg_vs_ber.m create mode 100644 projects/ECOC_2025/plots_from_database/plot_mpi_trial.m create mode 100644 projects/ECOC_2025/submit_dsp.m diff --git a/projects/ECOC_2025/b2b_measurement.m b/projects/ECOC_2025/b2b_measurement.m deleted file mode 100644 index e69de29..0000000 diff --git a/projects/ECOC_2025/dsp_run_id.m b/projects/ECOC_2025/dsp_run_id.m new file mode 100644 index 0000000..fb9f689 --- /dev/null +++ b/projects/ECOC_2025/dsp_run_id.m @@ -0,0 +1,235 @@ +function [output] = dsp_run_id(run_id,options) + +arguments + run_id + options.append_to_db = 0; + options.max_occurences = 4; + options.parameters = struct(); + + options.database_path + options.database_name + + options.storage_path +end + +database = DBHandler("pathToDB",[options.database_path,options.database_name]); + +filterParams = database.tables; +filterParams.Configurations = struct('run_id', run_id); + +selectedFields = {'Runs.run_id','Runs.tx_bits_path','Runs.tx_signal_path','Runs.tx_symbols_path','Runs.rx_sync_path','Runs.rx_raw_path',... + 'Configurations.db_mode','Configurations.pam_level','Configurations.bitrate','Configurations.symbolrate','Configurations.fiber_length','Configurations.wavelength','Configurations.precomp_amp','Measurements.power_rop','Configurations.v_bias',... + 'Configurations.interference_attenuation'}; + +[dataTable,sql_query] = database.queryDB(filterParams, selectedFields); +[~, uniqueIdx] = unique(dataTable.run_id); % Get unique run_id indices +dataTable = dataTable(uniqueIdx,:); % Extract unique configurations for each run_id + +fsym = dataTable.symbolrate; +M = double(dataTable.pam_level); +duob_mode = db_mode.(strrep(char(dataTable.db_mode),'"','')); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +len_tr = 4096*2; + +vnle_order1 = 50; +vnle_order2 = 5; +vnle_order3 = 5; + +dfe_order = [0 0 0]; + +pf_ncoeffs = 1; + +mu_ffe1 = 0.0001; +mu_ffe2 = 0.0008; +mu_ffe3 = 0.001; +mu_dfe = 0.0004; +mu_dc = 0.00; + +mu_ffe = [mu_ffe1 mu_ffe2 mu_ffe3]; +vnle_order=[vnle_order1,vnle_order2,vnle_order3]; + +% Overwrite default parameters if given in options.parameters +paramStruct = options.parameters; +if ~isempty(paramStruct) + paramNames = fieldnames(paramStruct); + for i = 1:numel(paramNames) + thisName = paramNames{i}; + thisValue = paramStruct.(thisName); + eval([thisName ' = thisValue;']); + end +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +eq_ = 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); +pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); +mlse_ = MLSE_viterbi("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); +mlse_db_ = MLSE_viterbi("DIR",[1,1],"duobinary_output",0,"M",M,"trellis_states",PAMmapper(M,0).levels); +eq_post = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",2001,"sps",1,"decide",0); + +output = struct(); +vnle_pf_package = {}; +vnle_dfe_package = {}; +dbtgt_package = {}; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +Tx_signal = load([options.storage_path, char(dataTable.tx_signal_path)]); +Tx_signal = Tx_signal.Digi_sig; + +Tx_bits = load([options.storage_path, char(dataTable.tx_bits_path)]); +Tx_bits = Tx_bits.Bits; +Symbols_mapped = PAMmapper(M,0).map(Tx_bits); +Symbols_mapped.fs = dataTable.symbolrate; + +Symbols = load([options.storage_path, char(dataTable.tx_symbols_path)]); +Symbols = Symbols.Symbols; + +found_sync = 0; + +try + Scpe_load = load([options.storage_path, char(dataTable.rx_sync_path)]); + Scpe_cell = Scpe_load.S; + [~,~,found_sync] = Scpe_cell{2}.tsynch("reference",Symbols,"fs_ref",dataTable.symbolrate,"debug_plots",0); +end + +if ~found_sync + Scpe_sig_raw = load([options.storage_path, char(dataTable.rx_raw_path(1))]); + Scpe_sig_raw = Scpe_sig_raw.Scpe_sig_raw; + Scpe_sig_resampled = Scpe_sig_raw.resample("fs_in",Scpe_sig_raw.fs,"fs_out",2*fsym); + [~,Scpe_cell,found_sync] =Scpe_sig_resampled.tsynch("reference",Symbols,"fs_ref",dataTable.symbolrate,"debug_plots",1); +end + +if ~found_sync + if length(Symbols_mapped.signal) == sum(Symbols_mapped.signal == Symbols.signal) + warning('Could not synchronize the received signal with the stored symbols!') + else + [~,Scpe_cell,found_sync] =Scpe_sig_raw.tsynch("reference",Symbols_mapped,"fs_ref",dataTable.symbolrate,"debug_plots",0); + end + if ~found_sync + warning('Could not synchronize the received signal with the stored symbols!') + end +end + +record_realizations = min(options.max_occurences,length(Scpe_cell)); +for occ = 1:record_realizations + + Scpe_sig = Scpe_cell{occ}; + + %%%%%% Sample to 2x fsym %%%%%% + Scpe_sig = Scpe_sig.resample("fs_out",2*fsym); + + %%%%%% Sync Rx signal with reference %%%%%% + [Scpe_sig,~] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym,"debug_plots",0); + + Scpe_sig = Filter('filtdegree',4,"f_cutoff",Symbols.fs.*0.6,"fs",Scpe_sig.fs,"filterType",filtertypes.gaussian,"active",true).process(Scpe_sig); + + Scpe_sig = Scpe_sig - mean(Scpe_sig.signal); + + if duob_mode ~= db_mode.db_encoded + + % %%%%% 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_post = 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",duob_mode,"showAnalysis",1,"postFFE",[]); + vnle_dfe_package{occ} = result; + + end + + %%%%% VNLE + PF + MLSE %%%% + if 1 + + [result] = vnle_postfilter_mlse(eq_,pf_,mlse_,M,Scpe_sig,Symbols,Tx_bits,"precode_mode",duob_mode,'showAnalysis',1,"postFFE",[],"eth_style_symbol_mapping",0); + vnle_pf_package{occ} = result; + + if options.append_to_db + database.addProcessingResult(run_id,result.resultsMLSE, result.equalizerConfigMLSE); + database.addProcessingResult(run_id,result.resultsVNLE, result.equalizerConfigVNLE); + end + end + + %%%%% Duobinary Targeting %%%% + if 1 + [result] = duobinary_target(eq_, mlse_db_, M, Scpe_sig, Symbols, Tx_bits, "precode_mode", duob_mode,'showAnalysis',0,"postFFE",[]); + dbtgt_package{occ} = result; + + if options.append_to_db + database.addProcessingResult(run_id, result.resultsDBtgt, result.equalizerConfigDBtgt); + end + end + + % fprintf("BER VNLE: %.2e | %.2e; BER MLSE: %.2e | %.2e; BER DB tgt: %.2e | %.2e \n",vnle_pf_package{occ}.resultsVNLE.BER,vnle_pf_package{occ}.resultsVNLE.BER_precoded ,vnle_pf_package{occ}.resultsMLSE.BER,vnle_pf_package{occ}.resultsMLSE.BER_precoded,dbtgt_package{occ}.resultsDBtgt.BER,dbtgt_package{occ}.resultsDBtgt.BER_precoded) + % % fprintf("BER VNLE: %.2e | %.2e; BER MLSE: %.2e | %.2e \n",vnle_pf_package{occ}.resultsVNLE.BER,vnle_pf_package{occ}.resultsVNLE.BER_precoded ,vnle_pf_package{occ}.resultsMLSE.BER,vnle_pf_package{occ}.resultsMLSE.BER_precoded); + + occ = 1; % or whatever your loop index is + + % Extract VNLE results for readability + vnle = vnle_pf_package{occ}.resultsVNLE; + mlse = vnle_pf_package{occ}.resultsMLSE; + dbtgt = dbtgt_package{occ}.resultsDBtgt; + + % Print header + fprintf("==== EQUALIZATION RUN-ID %d | PAM-%d | %.2f GBd ====\n\n", run_id, M, Symbols.fs.*1e-9); + + % VNLE Results + fprintf(">> VNLE Results:\n"); + fprintf(" BER %.2e\n", vnle.BER); + fprintf(" BER (pre-code) %.2e\n", vnle.BER_precoded); + fprintf(" SNR: %.2f dB\n", vnle.SNR); + fprintf(" GMI: %.4f\n", vnle.GMI); + fprintf(" Linerate: %.2f Gbps\n", Symbols.fs .* floor(log2(M)*10)/10 .*1e-9); + fprintf(" AIR: %.2f Gbps\n", vnle.AIR.*1e-9); + fprintf("\n"); + + % MLSE Results + fprintf(">> MLSE Results:\n"); + fprintf(" BER : %.2e\n", mlse.BER); + fprintf(" BER (pre-code): %.2e\n", mlse.BER_precoded); + fprintf(" Channel Alpha : %.2f\n", mlse.Alpha); + fprintf("\n"); + + % DB Target Results + fprintf(">> DB Target Results:\n"); + fprintf(" BER: %.2e\n", dbtgt.BER); + fprintf(" BER (pre-code): %.2e\n", dbtgt.BER_precoded); + fprintf("\n"); + + + else + + %%%%%% %db signaling => db encoded %%%%% + if 1 + 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; + + if options.append_to_db + database.addProcessingResult(run_id, result.resultsDBsignaling, result.equalizerConfigDBsignaling); + end + end + + fprintf("BER DB: %.2e \n",dbenc_package{occ}.resultsDBsignaling.BER); + + + end + + + % autoArrangeFigures; + disp('- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ') + fprintf('\n') + + +end + + +output.dataTable = dataTable; +output.vnle_pf_package = vnle_pf_package; +output.dbtgt_package = dbtgt_package; + +end \ No newline at end of file diff --git a/projects/ECOC_2025/dsp_standalone.m b/projects/ECOC_2025/dsp_standalone.m new file mode 100644 index 0000000..371ac40 --- /dev/null +++ b/projects/ECOC_2025/dsp_standalone.m @@ -0,0 +1,10 @@ + + + +basePath = 'C:\Users\sioe\Documents\MATLAB\imdd_simulation\projects\ECOC_2025\'; +database_name='ecoc2025.db'; +savePath = 'Z:\2025\ECOC Silas\ecoc_2025\'; +num_occ = 1; +run_par = false; +run_id = 562; +[out, ~] = submit_dsp(run_id, basePath, database_name, savePath,"parallel",run_par,'max_occurences',num_occ); diff --git a/projects/ECOC_2025/measure_run.m b/projects/ECOC_2025/measure_run.m new file mode 100644 index 0000000..1a9a189 --- /dev/null +++ b/projects/ECOC_2025/measure_run.m @@ -0,0 +1,207 @@ + +function measure_run(options) + +arguments + options.parallel (1,1) logical = false + options.max_occurences = 1; + options.config = []; + options.only_record = 0; +end + + savePath = 'Z:\2025\ECOC Silas\ecoc_2025\'; + + current_folder = 'precomp_optimization_friday'; + fullFolderPath = fullfile(savePath, current_folder); + if ~exist(fullFolderPath, 'dir') + mkdir(fullFolderPath); + end + + basePath = 'C:\Users\sioe\Documents\MATLAB\imdd_simulation\projects\ECOC_2025\'; + database_name='ecoc2025.db'; + db = DBHandler("pathToDB",[basePath,database_name]); + + precomp_path = "C:\Users\sioe\Documents\MATLAB\imdd_simulation\projects\ECOC_2025\"; + precomp_fn = "precomp_1km_1mm_cable_70ghz_pd_shf_t850_2p8_bias_shot2"; + precomp_mode = 2; %0=do nothing ; 1= measure; 2=precomp active + run_lab_automation = 1; + + %% Configuration + + conf = db.tables.Configurations; + referenceFields = fieldnames(conf); + + conf.unique_elab_id = "20250408-842b0a3078172b48dd032795226fbe683190afc4"; + conf.symbolrate = 112e9; + conf.bitrate = conf.symbolrate .* floor(log2(6)*10)/10; + conf.pam_level = 4; + conf.db_mode = db_mode.no_db; + conf.pulsef_alpha = 0.3; + conf.v_bias = 2.8; + conf.v_awg = 0.8; + conf.precomp_amp = -64; + conf.wavelength = 1310; + conf.laser_power = 10; + conf.fiber_length = 0; + conf.pd_in_desired = 7; + conf.is_mpi = 1; + conf.signal_attenuation = 0; + conf.interference_path_length = 0.001; %m + conf.interference_attenuation = 40; + conf.pam_source = []; + + + % === Optional Override of Conf === + if ~isempty(options.config) + paramStruct = options.config; + paramNames = fieldnames(paramStruct); + for i = 1:numel(paramNames) + thisName = paramNames{i}; + thisValue = paramStruct.(thisName); + if isfield(conf, thisName) + conf.(thisName) = thisValue; + else + warning('Unknown parameter: %s', thisName); + end + end + end + + assert(isequal(fieldnames(conf), referenceFields), 'Fieldnames do not match the reference structure!'); + + currentTime = datetime('now', 'Format', 'yyyyMMdd_HHmmss'); + timeStr = char(currentTime); + filename = sprintf('%s_PAM_%d_R_%d',timeStr,conf.pam_level,conf.symbolrate.*1e-9); + + %% Lab Automation + + if run_lab_automation + + % DC Source + dcs = DC_supply("active",[1,1],"voltage",[conf.v_bias, 5]); + dcs.set("voltage",[conf.v_bias, 5]); + [v_bias_meas,i_bias_meas]=dcs.readVals(); + + % Laser + laser = Exfo_laser("mainframe_channel",1,"safety_mode",0,"connection_id",'10','lab_interface','gpib'); + laser.setWavelength(conf.wavelength); + laser.setPower(conf.laser_power); + laser.enableLaser(); + + % Optical Attenuator + voa = OptAtten("active",[0,2,1,1],"value",[0,conf.pd_in_desired,conf.signal_attenuation,conf.interference_attenuation],"wavelength",[conf.wavelength,conf.wavelength,conf.wavelength,conf.wavelength]); + voa.set('value',[0,conf.pd_in_desired,conf.signal_attenuation,conf.interference_attenuation]); + + % PDFA + pdfa = Thor_PDFA("safety_mode",0,"serialport_number",'COM15'); + + % Construct AWG and Scope Modules %%%%%% + fdac = 256e9; + fadc = 160e9; + Scp = ScopeKeysight("model","DSAZ634A",'autoscale',1,"fadc","GSa_160","channel",[0,0,1,0],"recordLen",2000000,"removeDC",1,"extRef",1); + Awg = AwgKeysight("model","M8199A_ILV","fdac",fdac,"scaletodac",[1,1],"skews",[0,0],"voltages",[0,conf.v_awg]); + A2S = Awg2Scope(Awg,Scp,[0,0,3,0],"waitUntilClick",0); + + end + + + %% Signal generation and transmission + + %%%%% Symbol Generation %%%%%% + Pform = Pulseformer("fsym",conf.symbolrate,"fdac",8*conf.symbolrate,"pulse","rc","pulselength",16,"alpha",conf.pulsef_alpha); + + Pamsource = PAMsource(... + "fsym",conf.symbolrate,"M",conf.pam_level,"order",18,"useprbs",1,... + "fs_out",fdac,... + "applyclipping",0,... + "clipfactor",4,... + "applypulseform",1,... + "pulseformer",Pform,... + "randkey",1,... + "duobinary_mode", conf.db_mode); + + conf.pam_source = Pamsource; + + [Digi_sig,Symbols,Bits] = Pamsource.process(); + + + %%%%% Precompensation Routine %%%%%% + if precomp_mode == 1 % measure channel + precomp_est = ChannelFreqResp("Nacq",2048,"Navg",100,"Ncp",63,'f_ref',fdac); + Digi_sig = precomp_est.buildOFDM(); + elseif precomp_mode == 2 % apply precomp + precomp_est = ChannelFreqResp("Nacq",2048,"Navg",100,"Ncp",63,'f_ref',Digi_sig.fs); + Digi_sig = precomp_est.precomp(Digi_sig,'maxampdb',conf.precomp_amp,'loadPath',precomp_path,'fileName',precomp_fn); + end + + Digi_sig.spectrum("displayname",'Tx Spectrum','fignum',10,'normalizeTo0dB',0); + + %%%%% Resample to DAC rate %%%%%% + Digi_sig = Digi_sig.resample("fs_out",Awg.fdac); + + [~,~,Scpe_sig_raw,~] = A2S.process("signal2",Digi_sig); + + % Awg.upload("signal2",Digi_sig); + % Scpe_sig_raw = Scp.read("channel",[0,0,1,0]); + % Scpe_sig_raw = Scpe_sig_raw{[0,0,1,0]==1}; + + %%%%% Precompensation Routine %%%%%% + if precomp_mode == 1 + % Scpe_sig_resampled = Scpe_sig_raw.resample("fs_in",fadc,"fs_out",2*fsym); + precomp_est.estimate(Scpe_sig_raw,"save",true,"savePath",precomp_path,"fileName",precomp_fn); + precomp_est.plot(); + return; + end + + + Scpe_sig_raw.spectrum("displayname",'Rx Spectrum','fignum',10,'normalizeTo0dB',0); + Scpe_sig_raw.plot("clear",1,"displayname",'Rx Signal','fignum',11); + + %% %%%%% Save Routine %%%%%%%%%%%%%%%%%%%%%%%%% + + tx_bits_path=[filesep,current_folder,filesep,filename,'_bits']; + save([savePath,tx_bits_path],"Bits"); + tx_symbols_path=[filesep,current_folder,filesep,filename,'_symbols']; + save([savePath,tx_symbols_path],"Symbols"); + tx_signal_path=[filesep,current_folder,filesep,filename,'_signal']; + save([savePath,tx_signal_path],"Digi_sig"); + rx_raw_path = [filesep,current_folder,filesep,filename,'_rx_signal_raw']; + save([savePath,rx_raw_path],"Scpe_sig_raw"); + + + % Table 1: Append to Runs + newRun = db.tables.Runs; + newRun.run_id = NaN; + newRun.date_of_run = datetime(currentTime, 'InputFormat', 'yyyyMMdd_HHmmss'); + newRun.tx_bits_path = tx_bits_path; + newRun.tx_symbols_path = tx_symbols_path; + newRun.rx_sync_path = []; + newRun.rx_raw_path = rx_raw_path; + newRun.filename = filename; + newRun.tx_signal_path = tx_signal_path; + run_id = db.appendToTable('Runs', newRun); + + % Table 2: Append to Configurations + conf.configuration_id = NaN; + conf.run_id = run_id; + db.appendToTable('Configurations', conf); + + + % Table 3: Append to Measurements + meas = db.tables.Measurements; + meas.measurement_id = NaN; % Auto-increment, leave empty + meas.run_id = run_id; + meas.power_laser = laser.cur_power; + meas.power_rop = []; + meas.power_pd_in = voa.power_state(2); + meas.power_mpi_interference = voa.power_state(4); + meas.power_mpi_signal = voa.power_state(3); + meas.voa_class = voa; + meas.pdfa_class = pdfa; + meas.laser_class = laser; + assert(isequal(fieldnames(meas), fieldnames(db.tables.Measurements)), 'Fieldnames of "meas" do not match the reference structure!'); + db.appendToTable('Measurements', meas); + + %% %%%%% DSP Routine %%%%%%%%%%%%%%%%%%%%%%%%% + + [out, ~] = submit_dsp(run_id, basePath, database_name, savePath,"parallel",options.parallel,'max_occurences',options.max_occurences); + +end \ No newline at end of file diff --git a/projects/ECOC_2025/measure_run_script.m b/projects/ECOC_2025/measure_run_script.m new file mode 100644 index 0000000..be6ef25 --- /dev/null +++ b/projects/ECOC_2025/measure_run_script.m @@ -0,0 +1,172 @@ + + +savePath = 'Z:\2025\ECOC Silas\ecoc_2025\'; + +current_folder = 'precomp_optimization_friday'; +fullFolderPath = fullfile(savePath, current_folder); +if ~exist(fullFolderPath, 'dir') + mkdir(fullFolderPath); +end + +basePath = 'C:\Users\sioe\Documents\MATLAB\imdd_simulation\projects\ECOC_2025\'; +database_name='ecoc2025.db'; +db = DBHandler("pathToDB",[basePath,database_name]); + +precomp_path = "C:\Users\sioe\Documents\MATLAB\imdd_simulation\projects\ECOC_2025\"; +precomp_fn = "precomp_1km_1mm_cable_70ghz_pd_shf_t850_2p8_bias_shot2"; +precomp_mode = 2; %0=do nothing ; 1= measure; 2=precomp active +run_lab_automation = 1; + +%% Configuration + +conf = confRequest; +referenceFields = fieldnames(db.tables.Configurations); + +assert(isequal(fieldnames(conf), referenceFields), 'Fieldnames do not match the reference structure!'); + +currentTime = datetime('now', 'Format', 'yyyyMMdd_HHmmss'); +timeStr = char(currentTime); +filename = sprintf('%s_PAM_%d_R_%d',timeStr,conf.pam_level,conf.symbolrate.*1e-9); + +%% Lab Automation + +if run_lab_automation + + % DC Source + if confPrev.v_bias ~= conf.v_bias + dcs = DC_supply("active",[1,1],"voltage",[conf.v_bias, 5]); + dcs.set("voltage",[conf.v_bias, 5]); + [v_bias_meas,i_bias_meas]=dcs.readVals(); + end + + % Laser + if confPrev.laser_power ~= conf.laser_power + laser = Exfo_laser("mainframe_channel",1,"safety_mode",0,"connection_id",'10','lab_interface','gpib'); + laser.setWavelength(conf.wavelength); + laser.setPower(conf.laser_power); + laser.enableLaser(); + end + + % Optical Attenuator + voa = OptAtten("active",[0,2,1,1],"value",[0,conf.pd_in_desired,conf.signal_attenuation,conf.interference_attenuation],"wavelength",[conf.wavelength,conf.wavelength,conf.wavelength,conf.wavelength]); + voa.set('value',[0,conf.pd_in_desired,conf.signal_attenuation,conf.interference_attenuation]); + + % PDFA + pdfa = Thor_PDFA("safety_mode",0,"serialport_number",'COM15'); + + % Construct AWG and Scope Modules %%%%%% + fdac = 256e9; + fadc = 160e9; + Scp = ScopeKeysight("model","DSAZ634A",'autoscale',1,"fadc","GSa_160","channel",[0,0,1,0],"recordLen",5000000,"removeDC",1,"extRef",1); + Awg = AwgKeysight("model","M8199A_ILV","fdac",fdac,"scaletodac",[1,1],"skews",[0,0],"voltages",[0,conf.v_awg]); + A2S = Awg2Scope(Awg,Scp,[0,0,3,0],"waitUntilClick",0); + +end + + +%% Signal generation and transmission + +%%%%% Symbol Generation %%%%%% +Pform = Pulseformer("fsym",conf.symbolrate,"fdac",8*conf.symbolrate,"pulse","rc","pulselength",16,"alpha",conf.pulsef_alpha); + +Pamsource = PAMsource(... + "fsym",conf.symbolrate,"M",conf.pam_level,"order",19,"useprbs",1,... + "fs_out",fdac,... + "applyclipping",0,... + "clipfactor",4,... + "applypulseform",1,... + "pulseformer",Pform,... + "randkey",1,... + "duobinary_mode", conf.db_mode); + +conf.pam_source = Pamsource; +upload_required = ~isequal(Pamsource, confPrev.pam_source); +confPrev = conf; + +[Digi_sig,Symbols,Bits] = Pamsource.process(); + + +%%%%% Precompensation Routine %%%%%% +if precomp_mode == 1 % measure channel + precomp_est = ChannelFreqResp("Nacq",2048,"Navg",100,"Ncp",63,'f_ref',fdac); + Digi_sig = precomp_est.buildOFDM(); + upload_required = 1; +elseif precomp_mode == 2 % apply precomp + precomp_est = ChannelFreqResp("Nacq",2048,"Navg",100,"Ncp",63,'f_ref',Digi_sig.fs); + Digi_sig = precomp_est.precomp(Digi_sig,'maxampdb',conf.precomp_amp,'loadPath',precomp_path,'fileName',precomp_fn); +end + +Digi_sig.spectrum("displayname",'Tx Spectrum','fignum',10,'normalizeTo0dB',0); + +%%%%% Resample to DAC rate %%%%%% +Digi_sig = Digi_sig.resample("fs_out",Awg.fdac); + +% [~,~,Scpe_sig_raw,~] = A2S.process("signal2",Digi_sig); +if upload_required || confPrev.v_awg ~= conf.v_awg + Awg.upload("signal2",Digi_sig); +end + +Scpe_sig_raw = Scp.read("channel",[0,0,1,0]); +Scpe_sig_raw = Scpe_sig_raw{[0,0,1,0]==1}; + +%%%%% Precompensation Routine %%%%%% +if precomp_mode == 1 + % Scpe_sig_resampled = Scpe_sig_raw.resample("fs_in",fadc,"fs_out",2*fsym); + precomp_est.estimate(Scpe_sig_raw,"save",true,"savePath",precomp_path,"fileName",precomp_fn); + precomp_est.plot(); + return; +end + + +Scpe_sig_raw.spectrum("displayname",'Rx Spectrum','fignum',10,'normalizeTo0dB',0); +Scpe_sig_raw.plot("clear",1,"displayname",'Rx Signal','fignum',11); + +%% %%%%% Save Routine %%%%%%%%%%%%%%%%%%%%%%%%% + +tx_bits_path=[filesep,current_folder,filesep,filename,'_bits']; +save([savePath,tx_bits_path],"Bits"); +tx_symbols_path=[filesep,current_folder,filesep,filename,'_symbols']; +save([savePath,tx_symbols_path],"Symbols"); +tx_signal_path=[filesep,current_folder,filesep,filename,'_signal']; +save([savePath,tx_signal_path],"Digi_sig"); +rx_raw_path = [filesep,current_folder,filesep,filename,'_rx_signal_raw']; +save([savePath,rx_raw_path],"Scpe_sig_raw"); + + +% Table 1: Append to Runs +newRun = db.tables.Runs; +newRun.run_id = NaN; +newRun.date_of_run = datetime(currentTime, 'InputFormat', 'yyyyMMdd_HHmmss'); +newRun.tx_bits_path = tx_bits_path; +newRun.tx_symbols_path = tx_symbols_path; +newRun.rx_sync_path = []; +newRun.rx_raw_path = rx_raw_path; +newRun.filename = filename; +newRun.tx_signal_path = tx_signal_path; +run_id = db.appendToTable('Runs', newRun); + +% Table 2: Append to Configurations +conf.configuration_id = NaN; +conf.run_id = run_id; +db.appendToTable('Configurations', conf); + + +% Table 3: Append to Measurements +meas = db.tables.Measurements; +meas.measurement_id = NaN; % Auto-increment, leave empty +meas.run_id = run_id; +meas.power_laser = laser.cur_power; +meas.power_rop = []; +meas.power_pd_in = voa.power_state(2); +meas.power_mpi_interference = voa.power_state(4); +meas.power_mpi_signal = voa.power_state(3); +meas.voa_class = voa; +meas.pdfa_class = pdfa; +meas.laser_class = laser; +assert(isequal(fieldnames(meas), fieldnames(db.tables.Measurements)), 'Fieldnames of "meas" do not match the reference structure!'); +db.appendToTable('Measurements', meas); + +%% %%%%% DSP Routine %%%%%%%%%%%%%%%%%%%%%%%%% + +[out, ~] = submit_dsp(run_id, basePath, database_name, savePath,"parallel",parallel_dsp,'max_occurences',max_occurences); + diff --git a/projects/ECOC_2025/plots_from_database/bias_vs_ber.m b/projects/ECOC_2025/plots_from_database/bias_vs_ber.m new file mode 100644 index 0000000..60977a2 --- /dev/null +++ b/projects/ECOC_2025/plots_from_database/bias_vs_ber.m @@ -0,0 +1,269 @@ + + +basePath = 'C:\Users\sioe\Documents\MATLAB\imdd_simulation\projects\ECOC_2025\'; +database_name = 'ecoc2025.db'; +database = DBHandler("pathToDB", [basePath, database_name]); + +filterParams = database.tables; +filterParams.Configurations = struct( ... + 'symbolrate', 112e9, ... %[224,336,360,390,420,448] + 'fiber_length', 0, ... + 'db_mode', '"no_db"', ... + 'interference_attenuation', [], ... + 'interference_path_length', [], ... + 'is_mpi', 1, ... + 'pam_level', 4, ... + 'wavelength', 1310, ... + 'precomp_amp', -64, ... + 'signal_attenuation', '0', ... + 'v_awg', [], ... + 'v_bias', [] ... + ); + +% filterParams.EqualizerParameters.diff_precode = int32(db_mode.db_encoded); +% filterParams.EqualizerParameters.equalizer_structure = int32(equalizer_structure.vnle); + +selectedFields = {'Configurations.run_id' 'Runs.date_of_run' 'Runs.rx_raw_path' 'Configurations.bitrate' 'Configurations.v_bias' 'Configurations.v_awg' 'Configurations.precomp_amp' 'Configurations.symbolrate' 'Configurations.pam_level'... + 'Configurations.db_mode' 'Configurations.rop_attenuation' 'Configurations.is_mpi' 'Configurations.interference_attenuation' 'Configurations.signal_attenuation' ... + 'EqualizerParameters.equalizer_structure' 'EqualizerParameters.diff_precode' 'EqualizerParameters.eq_id' 'Measurements.power_pd_in' ... + 'Measurements.power_mpi_interference' 'Measurements.power_mpi_signal' 'Results.BER' 'Results.BER_precoded' 'Results.SNR' 'Results.GMI' 'Results.Alpha' 'Results.date_of_processing'}; + +[dataTable,sql_query] = database.queryDB(filterParams, selectedFields); + +dataTable = cleanUpTable(dataTable); + +% Filter by time +startTime = datetime('2025-04-11 13:00:00', 'InputFormat', 'yyyy-MM-dd HH:mm:ss'); +stopTime = datetime('2025-04-11 14:00:00', 'InputFormat', 'yyyy-MM-dd HH:mm:ss'); +dataTable.date_of_run = datetime(dataTable.date_of_run, 'InputFormat', 'yyyy-MM-dd HH:mm:ss.SSSSSS'); +dataTable.date_of_processing = datetime(dataTable.date_of_processing, 'InputFormat', 'yyyy-MM-dd HH:mm:ss.SSSSSS'); +dataTable = dataTable(dataTable.date_of_processing > startTime, :); +dataTable = dataTable(dataTable.date_of_processing < stopTime, :); + +% Group by smth +y_var = 'SNR'; +x_var = 'v_bias'; + +loop_var = 'v_awg'; +fixedVars = {'eq_id',loop_var}; +dataTableGrpd = groupIt(fixedVars,dataTable); + +plotRealizations = 1; + +% Create a new figure +mkr = 'x'; + +figure(10); +hold on + +unique_loop_var = unique(dataTable.(loop_var)); +cols = linspecer(8); + +for i = 1:numel(unique_loop_var) + + idx = find(dataTable.(loop_var)== unique_loop_var(i), 1, 'first'); + % dispname = equalizer_structure(dataTable.equalizer_structure(idx)); + dispname = num2str(unique_loop_var(i)); + + % Plot SCATTERS: timestamp vs. interference_attenuation + loop_filt_2 = dataTable.(loop_var)==unique_loop_var(i); + if plotRealizations + + y_values = dataTable.(y_var)(loop_filt_2,:); + x_values = dataTable.(x_var)(loop_filt_2,:); + y_values = double(y_values); + x_values = double(x_values); + + sc = scatter(x_values, y_values, 'LineWidth', 0.5,'Marker',mkr,'MarkerEdgeColor',cols(i,:),'HandleVisibility','on','DisplayName',string(dispname)); + + pair_one = {'Run ID', dataTable.run_id(loop_filt_2,:)}; + pair_two = {'Rate', dataTable.bitrate(loop_filt_2,:).*1e-9}; + pair_three = {'PD in', round(dataTable.power_pd_in(loop_filt_2,:),2)}; + addDatatips(sc, pair_one, pair_two,pair_three); + xticks(unique(x_values)); + + end +end + + +% Label the axes and add a title +legend('Interpreter','latex'); +xlabel(x_var); +ylabel(y_var); +title([x_var,' vs. ',y_var]); +if y_var == 'BER' + yline(3.8e-3,'LineWidth',1,'LineStyle','--','HandleVisibility','off'); + ylim([1e-4 0.5]); +end +% Enable grid for better readability +grid on; + +beautifyBERplot; + + + + +function resultTable = groupIt(fixedVars,dataTable) + + % Group by run_id and eq_id (adjust grouping keys as needed) + + [G, groupKeys] = findgroups(dataTable(:, fixedVars)); + + % Preallocate a cell array for aggregated data. + varNames = dataTable.Properties.VariableNames; + nVars = numel(varNames); + aggData = cell(height(groupKeys), nVars); + groupCount = zeros(height(groupKeys), 1); % To store the size of each group + + % Loop over each group. + for i = 1:height(groupKeys) + idx = (G == i); % Logical index for group i + groupCount(i) = sum(idx); % Count number of rows in this group + % For each variable in the table: + for j = 1:nVars + colData = dataTable.(varNames{j}); + if isnumeric(colData) + % For numeric data, compute the mean. + aggData{i, j} = min(colData(idx)); + else + % For non-numeric data, take the first entry. + if iscell(colData) + aggData{i, j} = colData{find(idx, 1)}; + else + aggData{i, j} = colData(find(idx, 1)); + end + end + end + end + + % Convert the aggregated cell array into a table. + resultTable = cell2table(aggData, 'VariableNames', varNames); + + % Append the group count as a new column. + resultTable.nRows = groupCount; + +end + + +function addDatatips(sc, varargin) +% addDatatips Adds custom data tip rows to a scatter plot. +% +% addDatatips(sc, pair1, pair2, ...) adds one or more custom rows to the +% data tip display of the scatter plot identified by sc. +% +% Each pair should be provided as a 1x2 cell array: {label, value}. +% The value can be a scalar or a vector. If a vector is provided, its length +% must match the number of scatter plot points. +% +% Example: +% sc = scatter(x, y, 'LineWidth', 1.5, 'Marker', 'o'); +% pair_one = {'Attenuation', attenuationVector}; +% addDatatips(sc, pair_one); + + numPoints = numel(sc.XData); + + for k = 1:length(varargin) + pair = varargin{k}; + + if ~iscell(pair) || numel(pair) ~= 2 + error('Each pair must be a 1x2 cell array: {label, value}.'); + end + + label = pair{1}; + value = pair{2}; + + % If value is a vector, ensure its length is either 1 or equal to the number of scatter points. + if isvector(value) && numel(value) ~= 1 && numel(value) ~= numPoints + error('The vector for "%s" must be a scalar or have %d elements matching the scatter data points.', label, numPoints); + end + + % Create a new data tip row using the provided label and vector. + newRow = dataTipTextRow(label, value); + sc.DataTipTemplate.DataTipRows(end+1) = newRow; + end +end + + + +function cleanedTable = cleanUpTable(inputTable) +% cleanUpTable Cleans a MATLAB table where numbers and NaNs are stored as strings or structs. +% +% cleanedTable = cleanUpTable(inputTable) +% +% This function goes through all columns of the input table: +% - Converts strings of numbers to numeric values +% - Converts string 'NaN' and struct NaNs to real NaN +% - Converts date strings to datetime (if possible) +% +% Input: +% inputTable - MATLAB table with mixed types +% +% Output: +% cleanedTable - Cleaned MATLAB table with proper numeric types + + cleanedTable = inputTable; + varNames = cleanedTable.Properties.VariableNames; + + for i = 1:numel(varNames) + col = cleanedTable.(varNames{i}); + + % Case 1: If it's a cell array (likely mixed strings/struct) + if iscell(col) + % Convert struct 'NaN' entries to string 'NaN' + col = cellfun(@(x) convertStructToString(x), col, 'UniformOutput', false); + + % Try to convert string numbers to actual numbers + numericCol = str2double(col); + + if all(isnan(numericCol) == strcmpi(col, 'NaN') | cellfun(@isempty, col)) + % If conversion is successful (NaNs correspond to 'NaN' strings), use it + cleanedTable.(varNames{i}) = numericCol; + else + % Else, try to convert to datetime + try + cleanedTable.(varNames{i}) = datetime(col, 'InputFormat', 'yyyy-MM-dd HH:mm:ss.SSSSSS'); + catch + % If it fails, leave as cell array of strings + cleanedTable.(varNames{i}) = string(col); + end + end + + % Case 2: If it's already a string array + elseif isstring(col) + numericCol = str2double(col); + if all(isnan(numericCol) == strcmpi(col, "NaN")) + cleanedTable.(varNames{i}) = numericCol; + else + % Try convert to datetime + try + cleanedTable.(varNames{i}) = datetime(col, 'InputFormat', 'yyyy-MM-dd HH:mm:ss.SSSSSS'); + catch + % Leave as string + end + end + + % Case 3: If it's already numeric, keep as is + elseif isnumeric(col) + continue; + + % Case 4: If it's datetime, keep as is + elseif isdatetime(col) + continue; + + else + % Catch-all for unexpected types, convert to string + cleanedTable.(varNames{i}) = string(col); + end + end +end + +function out = convertStructToString(x) +% Helper function to convert struct NaN to string 'NaN' + if isstruct(x) + out = "NaN"; + elseif isstring(x) || ischar(x) + out = string(x); + else + out = x; + end +end diff --git a/projects/ECOC_2025/plots_from_database/plot_bias_awg_vs_ber.m b/projects/ECOC_2025/plots_from_database/plot_bias_awg_vs_ber.m new file mode 100644 index 0000000..3603ac1 --- /dev/null +++ b/projects/ECOC_2025/plots_from_database/plot_bias_awg_vs_ber.m @@ -0,0 +1,271 @@ + + +basePath = 'C:\Users\sioe\Documents\MATLAB\imdd_simulation\projects\ECOC_2025\'; +database_name = 'ecoc2025.db'; +database = DBHandler("pathToDB", [basePath, database_name]); + +filterParams = database.tables; +filterParams.Configurations = struct( ... + 'symbolrate', 112e9, ... %[224,336,360,390,420,448] + 'fiber_length', 0, ... + 'db_mode', '"no_db"', ... + 'interference_attenuation', [], ... + 'interference_path_length', [], ... + 'is_mpi', 0, ... + 'pam_level', 4, ... + 'wavelength', 1310, ... + 'precomp_amp', [], ... + 'signal_attenuation', [], ... + 'v_awg', 0.8, ... + 'v_bias', 2.8 ... + ); + +% filterParams.EqualizerParameters.diff_precode = int32(db_mode.db_encoded); +% filterParams.EqualizerParameters.equalizer_structure = int32(equalizer_structure.vnle); + +selectedFields = {'Configurations.run_id' 'Runs.date_of_run' 'Runs.rx_raw_path' 'Configurations.bitrate' 'Configurations.v_bias' 'Configurations.v_awg' 'Configurations.precomp_amp' 'Configurations.symbolrate' 'Configurations.pam_level'... + 'Configurations.db_mode' 'Configurations.rop_attenuation' 'Configurations.is_mpi' 'Configurations.interference_attenuation' 'Configurations.signal_attenuation' ... + 'EqualizerParameters.equalizer_structure' 'EqualizerParameters.diff_precode' 'EqualizerParameters.eq_id' 'Measurements.power_pd_in' ... + 'Measurements.power_mpi_interference' 'Measurements.power_mpi_signal' 'Results.BER' 'Results.BER_precoded' 'Results.SNR' 'Results.GMI' 'Results.Alpha' 'Results.date_of_processing'}; + +[dataTable,sql_query] = database.queryDB(filterParams, selectedFields); + +dataTable = cleanUpTable(dataTable); + +% Filter by time +startTime = datetime('2025-04-11 13:00:00', 'InputFormat', 'yyyy-MM-dd HH:mm:ss'); +stopTime = datetime('2025-04-11 14:00:00', 'InputFormat', 'yyyy-MM-dd HH:mm:ss'); +dataTable.date_of_run = datetime(dataTable.date_of_run, 'InputFormat', 'yyyy-MM-dd HH:mm:ss.SSSSSS'); +dataTable.date_of_processing = datetime(dataTable.date_of_processing, 'InputFormat', 'yyyy-MM-dd HH:mm:ss.SSSSSS'); +dataTable = dataTable(dataTable.date_of_processing > startTime, :); +dataTable = dataTable(dataTable.date_of_processing < stopTime, :); + +% Group by smth +y_var = 'BER'; +x_var = 'interference_attenuation'; + +loop_var = 'eq_id'; +fixedVars = {'eq_id',loop_var}; +dataTableGrpd = groupIt(fixedVars,dataTable); + +plotRealizations = 1; + +% Create a new figure +mkr = '*'; + +figure(10); +hold on + +unique_loop_var = unique(dataTable.(loop_var)); +cols = linspecer(8); + +for i = 1:numel(unique_loop_var) + + idx = find(dataTable.(loop_var)== unique_loop_var(i), 1, 'first'); + + dispname = equalizer_structure(dataTable.equalizer_structure(idx)); + + % dispname = num2str(unique_loop_var(i)); + + % Plot SCATTERS: timestamp vs. interference_attenuation + loop_filt_2 = dataTable.(loop_var)==unique_loop_var(i); + if plotRealizations + + y_values = dataTable.(y_var)(loop_filt_2,:); + x_values = dataTable.(x_var)(loop_filt_2,:); + y_values = double(y_values); + x_values = double(x_values); + + sc = scatter(x_values, y_values, 'LineWidth', 0.5,'Marker',mkr,'MarkerEdgeColor',cols(i,:),'HandleVisibility','on','DisplayName',string(dispname)); + + pair_one = {'Run ID', dataTable.run_id(loop_filt_2,:)}; + pair_two = {'Rate', dataTable.bitrate(loop_filt_2,:).*1e-9}; + pair_three = {'PD in', round(dataTable.power_pd_in(loop_filt_2,:),2)}; + addDatatips(sc, pair_one, pair_two,pair_three); + xticks(unique(x_values)); + + end +end + + +% Label the axes and add a title +legend('Interpreter','latex'); +xlabel(x_var); +ylabel(y_var); +title([x_var,' vs. ',y_var]); +if y_var == 'BER' + yline(3.8e-3,'LineWidth',1,'LineStyle','--','HandleVisibility','off'); + ylim([1e-5 0.1]); +end +% Enable grid for better readability +grid on; + +beautifyBERplot; + + + + +function resultTable = groupIt(fixedVars,dataTable) + + % Group by run_id and eq_id (adjust grouping keys as needed) + + [G, groupKeys] = findgroups(dataTable(:, fixedVars)); + + % Preallocate a cell array for aggregated data. + varNames = dataTable.Properties.VariableNames; + nVars = numel(varNames); + aggData = cell(height(groupKeys), nVars); + groupCount = zeros(height(groupKeys), 1); % To store the size of each group + + % Loop over each group. + for i = 1:height(groupKeys) + idx = (G == i); % Logical index for group i + groupCount(i) = sum(idx); % Count number of rows in this group + % For each variable in the table: + for j = 1:nVars + colData = dataTable.(varNames{j}); + if isnumeric(colData) + % For numeric data, compute the mean. + aggData{i, j} = min(colData(idx)); + else + % For non-numeric data, take the first entry. + if iscell(colData) + aggData{i, j} = colData{find(idx, 1)}; + else + aggData{i, j} = colData(find(idx, 1)); + end + end + end + end + + % Convert the aggregated cell array into a table. + resultTable = cell2table(aggData, 'VariableNames', varNames); + + % Append the group count as a new column. + resultTable.nRows = groupCount; + +end + + +function addDatatips(sc, varargin) +% addDatatips Adds custom data tip rows to a scatter plot. +% +% addDatatips(sc, pair1, pair2, ...) adds one or more custom rows to the +% data tip display of the scatter plot identified by sc. +% +% Each pair should be provided as a 1x2 cell array: {label, value}. +% The value can be a scalar or a vector. If a vector is provided, its length +% must match the number of scatter plot points. +% +% Example: +% sc = scatter(x, y, 'LineWidth', 1.5, 'Marker', 'o'); +% pair_one = {'Attenuation', attenuationVector}; +% addDatatips(sc, pair_one); + + numPoints = numel(sc.XData); + + for k = 1:length(varargin) + pair = varargin{k}; + + if ~iscell(pair) || numel(pair) ~= 2 + error('Each pair must be a 1x2 cell array: {label, value}.'); + end + + label = pair{1}; + value = pair{2}; + + % If value is a vector, ensure its length is either 1 or equal to the number of scatter points. + if isvector(value) && numel(value) ~= 1 && numel(value) ~= numPoints + error('The vector for "%s" must be a scalar or have %d elements matching the scatter data points.', label, numPoints); + end + + % Create a new data tip row using the provided label and vector. + newRow = dataTipTextRow(label, value); + sc.DataTipTemplate.DataTipRows(end+1) = newRow; + end +end + + + +function cleanedTable = cleanUpTable(inputTable) +% cleanUpTable Cleans a MATLAB table where numbers and NaNs are stored as strings or structs. +% +% cleanedTable = cleanUpTable(inputTable) +% +% This function goes through all columns of the input table: +% - Converts strings of numbers to numeric values +% - Converts string 'NaN' and struct NaNs to real NaN +% - Converts date strings to datetime (if possible) +% +% Input: +% inputTable - MATLAB table with mixed types +% +% Output: +% cleanedTable - Cleaned MATLAB table with proper numeric types + + cleanedTable = inputTable; + varNames = cleanedTable.Properties.VariableNames; + + for i = 1:numel(varNames) + col = cleanedTable.(varNames{i}); + + % Case 1: If it's a cell array (likely mixed strings/struct) + if iscell(col) + % Convert struct 'NaN' entries to string 'NaN' + col = cellfun(@(x) convertStructToString(x), col, 'UniformOutput', false); + + % Try to convert string numbers to actual numbers + numericCol = str2double(col); + + if all(isnan(numericCol) == strcmpi(col, 'NaN') | cellfun(@isempty, col)) + % If conversion is successful (NaNs correspond to 'NaN' strings), use it + cleanedTable.(varNames{i}) = numericCol; + else + % Else, try to convert to datetime + try + cleanedTable.(varNames{i}) = datetime(col, 'InputFormat', 'yyyy-MM-dd HH:mm:ss.SSSSSS'); + catch + % If it fails, leave as cell array of strings + cleanedTable.(varNames{i}) = string(col); + end + end + + % Case 2: If it's already a string array + elseif isstring(col) + numericCol = str2double(col); + if all(isnan(numericCol) == strcmpi(col, "NaN")) + cleanedTable.(varNames{i}) = numericCol; + else + % Try convert to datetime + try + cleanedTable.(varNames{i}) = datetime(col, 'InputFormat', 'yyyy-MM-dd HH:mm:ss.SSSSSS'); + catch + % Leave as string + end + end + + % Case 3: If it's already numeric, keep as is + elseif isnumeric(col) + continue; + + % Case 4: If it's datetime, keep as is + elseif isdatetime(col) + continue; + + else + % Catch-all for unexpected types, convert to string + cleanedTable.(varNames{i}) = string(col); + end + end +end + +function out = convertStructToString(x) +% Helper function to convert struct NaN to string 'NaN' + if isstruct(x) + out = "NaN"; + elseif isstring(x) || ischar(x) + out = string(x); + else + out = x; + end +end diff --git a/projects/ECOC_2025/plots_from_database/plot_mpi_trial.m b/projects/ECOC_2025/plots_from_database/plot_mpi_trial.m new file mode 100644 index 0000000..9747fa2 --- /dev/null +++ b/projects/ECOC_2025/plots_from_database/plot_mpi_trial.m @@ -0,0 +1,342 @@ + + +basePath = 'C:\Users\sioe\Documents\MATLAB\imdd_simulation\projects\ECOC_2025\'; +database_name = 'ecoc2025.db'; +database = DBHandler("pathToDB", [basePath, database_name]); + +filterParams = database.tables; +filterParams.Configurations = struct( ... + 'symbolrate', 112e9, ... %[224,336,360,390,420,448] + 'fiber_length', 0, ... + 'db_mode', '"no_db"', ... + 'interference_attenuation', [], ... + 'interference_path_length', 0.06, ... + 'is_mpi', 1, ... + 'pam_level', 4, ... + 'wavelength', 1310, ... + 'precomp_amp', [], ... + 'signal_attenuation', [], ... + 'v_awg', 0.9, ... + 'v_bias', 2.8 ... + ); + +% filterParams.EqualizerParameters.diff_precode = int32(db_mode.db_encoded); +% filterParams.EqualizerParameters.equalizer_structure = int32(equalizer_structure.vnle); + +selectedFields = {'Configurations.run_id' 'Runs.date_of_run' 'Runs.rx_raw_path' 'Configurations.bitrate' 'Configurations.v_bias' 'Configurations.v_awg' 'Configurations.precomp_amp' 'Configurations.symbolrate' 'Configurations.pam_level'... + 'Configurations.db_mode' 'Configurations.rop_attenuation' 'Configurations.is_mpi' 'Configurations.interference_attenuation' 'Configurations.signal_attenuation' ... + 'EqualizerParameters.equalizer_structure' 'EqualizerParameters.diff_precode' 'EqualizerParameters.eq_id' 'Measurements.power_pd_in' ... + 'Measurements.power_mpi_interference' 'Measurements.power_mpi_signal' 'Results.BER' 'Results.BER_precoded' 'Results.SNR' 'Results.GMI' 'Results.Alpha' 'Results.date_of_processing'}; + +[dataTable,sql_query] = database.queryDB(filterParams, selectedFields); +dataTable.SIR = dataTable.power_mpi_signal - dataTable.power_mpi_interference; +dataTable = cleanUpTable(dataTable); + +% Filter by time +startTime = datetime('2025-04-11 14:40:00', 'InputFormat', 'yyyy-MM-dd HH:mm:ss'); +stopTime = datetime('2025-04-11 15:30:00', 'InputFormat', 'yyyy-MM-dd HH:mm:ss'); +dataTable.date_of_run = datetime(dataTable.date_of_run, 'InputFormat', 'yyyy-MM-dd HH:mm:ss.SSSSSS'); +dataTable.date_of_processing = datetime(dataTable.date_of_processing, 'InputFormat', 'yyyy-MM-dd HH:mm:ss.SSSSSS'); +dataTable = dataTable(dataTable.date_of_processing > startTime, :); +dataTable = dataTable(dataTable.date_of_processing < stopTime, :); + +% Group by smth +y_var = 'BER'; +x_var = 'SIR'; + +loop_var = 'eq_id'; +fixedVars = {'eq_id',x_var}; +% dataTableGrpd_mean = groupIt(fixedVars,dataTable); +dataTableGrpd_mean = groupIt(fixedVars, dataTable, @mean); +dataTableGrpd_min = groupIt(fixedVars, dataTable, @min); +dataTableGrpd_max = groupIt(fixedVars, dataTable, @max); + +plotRealizations = 1; +% Create a new figure +mkr = '*'; +figure(); +hold on + +unique_loop_var = unique(dataTable.(loop_var)); +cols = linspecer(numel(unique_loop_var)); % Ensure color count matches + +for i = 1:numel(unique_loop_var) + + % Prepare filtered data for this loop variable + loopValue = unique_loop_var(i); + + loopFiltGrpd = dataTableGrpd_mean.(loop_var) == loopValue; + if ~any(loopFiltGrpd) + continue; % Skip if no data for this loop var + end + + % Extract values + x_values = dataTableGrpd_mean.(x_var)(loopFiltGrpd, :); + y_mean = dataTableGrpd_mean.(y_var)(loopFiltGrpd, :); + y_min = dataTableGrpd_min.(y_var)(loopFiltGrpd, :); + y_max = dataTableGrpd_max.(y_var)(loopFiltGrpd, :); + + % Compute bounds: distance from mean + y_lower = y_mean - y_min; + y_upper = y_max - y_mean; + y_bounds = [y_lower, y_upper]; + + % Display name (optional) + idx = find(dataTable.(loop_var) == loopValue, 1, 'first'); + dispname = equalizer_structure(dataTable.equalizer_structure(idx)); + + % Plot bounded line + [hl, hp] = boundedline(x_values, y_mean, y_bounds, ... + 'alpha', 'transparency', 0.2, ... + 'cmap', cols(i,:), ... + 'nan', 'fill', ... + 'orientation', 'vert'); + + % Style the main line: thinnest, dotted, no marker + set(hl, 'LineWidth', 1, 'LineStyle', '-', 'Marker', 'none', ... + 'Color', cols(i,:), 'DisplayName', string(dispname)); + + % Hide patch (shaded area) from legend + set(hp, 'HandleVisibility', 'off','LineStyle',':','LineWidth',0.5,'Marker','none'); + + % Add invisible scatter for DataTips + sc_fake = scatter(x_values, y_mean, ... + 'Marker', 'o', 'MarkerEdgeColor', 'none', 'MarkerFaceColor', 'none', ... + 'HandleVisibility', 'off', 'PickableParts', 'all'); + + % Add data tips to the invisible scatter + pair_one = {'Run ID', dataTableGrpd_mean.run_id(loopFiltGrpd, :)}; + pair_two = {'Rate', dataTableGrpd_mean.bitrate(loopFiltGrpd, :) * 1e-9}; + pair_three = {'PD in', round(dataTableGrpd_mean.power_pd_in(loopFiltGrpd, :), 2)}; + addDatatips(sc_fake, pair_one, pair_two, pair_three); + + % Optionally: outline bounds for better visibility (optional) + % hnew = outlinebounds(hl, hp); + + % Tick marks (x-axis) + xticks(round(unique(x_values),2)); + + % Optional: scatter realizations + if plotRealizations + loopFiltSingle = dataTable.(loop_var) == loopValue; + x_single = double(dataTable.(x_var)(loopFiltSingle, :)); + y_single = double(dataTable.(y_var)(loopFiltSingle, :)); + + sc = scatter(x_single, y_single, 'Marker', mkr, 'MarkerEdgeColor', cols(i, :), ... + 'LineWidth', 0.5, 'HandleVisibility', 'on', 'DisplayName', string(dispname)); + + pair_one = {'Run ID', dataTable.run_id(loopFiltSingle, :)}; + pair_two = {'Rate', dataTable.bitrate(loopFiltSingle, :) * 1e-9}; + pair_three = {'PD in', round(dataTable.power_pd_in(loopFiltSingle, :), 2)}; + addDatatips(sc, pair_one, pair_two, pair_three); + end +end + +% Label axes and title +legend('Interpreter', 'latex'); +xlabel(x_var); +ylabel(y_var); +title([x_var, ' vs. ', y_var]); + +if y_var == 'BER' + yline(3.8e-3, 'LineWidth', 1, 'LineStyle', '--', 'HandleVisibility', 'off'); + ylim([1e-5, 0.1]); +end + +% Enable grid and beautify +grid on; +beautifyBERplot; + + + + +function resultTable = groupIt(fixedVars, dataTable, aggregationFunction) +% groupIt Groups data in a table based on fixedVars and applies aggregationFunction to numeric data. +% +% resultTable = groupIt(fixedVars, dataTable, aggregationFunction) +% +% Inputs: +% fixedVars - Cell array of variable names to group by +% dataTable - Input MATLAB table +% aggregationFunction - Function handle (e.g., @mean, @min, @max) +% +% Output: +% resultTable - Grouped and aggregated table + + % Group data + [G, groupKeys] = findgroups(dataTable(:, fixedVars)); + + % Prepare aggregation + varNames = dataTable.Properties.VariableNames; + nVars = numel(varNames); + aggData = cell(height(groupKeys), nVars); + groupCount = zeros(height(groupKeys), 1); % Store number of rows in each group + + % Loop over groups + for i = 1:height(groupKeys) + idx = (G == i); % Logical index for group i + groupCount(i) = sum(idx); % Count rows in group + + % Loop over each variable + for j = 1:nVars + colData = dataTable.(varNames{j}); + + if isnumeric(colData) + % Numeric: apply aggregation function (skip empty groups safely) + if any(idx) + aggData{i, j} = aggregationFunction(colData(idx)); + else + aggData{i, j} = NaN; + end + else + % Non-numeric: take first non-empty value + if iscell(colData) + nonEmptyIdx = find(idx & ~cellfun(@isempty, colData), 1); + if ~isempty(nonEmptyIdx) + aggData{i, j} = colData{nonEmptyIdx}; + else + aggData{i, j} = []; + end + else + nonEmptyIdx = find(idx, 1); + if ~isempty(nonEmptyIdx) + aggData{i, j} = colData(nonEmptyIdx); + else + aggData{i, j} = []; + end + end + end + end + end + + % Convert aggregated data to table + resultTable = cell2table(aggData, 'VariableNames', varNames); + + % Add group size as new column + resultTable.nRows = groupCount; +end + + + +function addDatatips(sc, varargin) +% addDatatips Adds custom data tip rows to a scatter plot. +% +% addDatatips(sc, pair1, pair2, ...) adds one or more custom rows to the +% data tip display of the scatter plot identified by sc. +% +% Each pair should be provided as a 1x2 cell array: {label, value}. +% The value can be a scalar or a vector. If a vector is provided, its length +% must match the number of scatter plot points. +% +% Example: +% sc = scatter(x, y, 'LineWidth', 1.5, 'Marker', 'o'); +% pair_one = {'Attenuation', attenuationVector}; +% addDatatips(sc, pair_one); + + numPoints = numel(sc.XData); + + for k = 1:length(varargin) + pair = varargin{k}; + + if ~iscell(pair) || numel(pair) ~= 2 + error('Each pair must be a 1x2 cell array: {label, value}.'); + end + + label = pair{1}; + value = pair{2}; + + % If value is a vector, ensure its length is either 1 or equal to the number of scatter points. + if isvector(value) && numel(value) ~= 1 && numel(value) ~= numPoints + error('The vector for "%s" must be a scalar or have %d elements matching the scatter data points.', label, numPoints); + end + + % Create a new data tip row using the provided label and vector. + newRow = dataTipTextRow(label, value); + sc.DataTipTemplate.DataTipRows(end+1) = newRow; + end +end + + + +function cleanedTable = cleanUpTable(inputTable) +% cleanUpTable Cleans a MATLAB table where numbers and NaNs are stored as strings or structs. +% +% cleanedTable = cleanUpTable(inputTable) +% +% This function goes through all columns of the input table: +% - Converts strings of numbers to numeric values +% - Converts string 'NaN' and struct NaNs to real NaN +% - Converts date strings to datetime (if possible) +% +% Input: +% inputTable - MATLAB table with mixed types +% +% Output: +% cleanedTable - Cleaned MATLAB table with proper numeric types + + cleanedTable = inputTable; + varNames = cleanedTable.Properties.VariableNames; + + for i = 1:numel(varNames) + col = cleanedTable.(varNames{i}); + + % Case 1: If it's a cell array (likely mixed strings/struct) + if iscell(col) + % Convert struct 'NaN' entries to string 'NaN' + col = cellfun(@(x) convertStructToString(x), col, 'UniformOutput', false); + + % Try to convert string numbers to actual numbers + numericCol = str2double(col); + + if all(isnan(numericCol) == strcmpi(col, 'NaN') | cellfun(@isempty, col)) + % If conversion is successful (NaNs correspond to 'NaN' strings), use it + cleanedTable.(varNames{i}) = numericCol; + else + % Else, try to convert to datetime + try + cleanedTable.(varNames{i}) = datetime(col, 'InputFormat', 'yyyy-MM-dd HH:mm:ss.SSSSSS'); + catch + % If it fails, leave as cell array of strings + cleanedTable.(varNames{i}) = string(col); + end + end + + % Case 2: If it's already a string array + elseif isstring(col) + numericCol = str2double(col); + if all(isnan(numericCol) == strcmpi(col, "NaN")) + cleanedTable.(varNames{i}) = numericCol; + else + % Try convert to datetime + try + cleanedTable.(varNames{i}) = datetime(col, 'InputFormat', 'yyyy-MM-dd HH:mm:ss.SSSSSS'); + catch + % Leave as string + end + end + + % Case 3: If it's already numeric, keep as is + elseif isnumeric(col) + continue; + + % Case 4: If it's datetime, keep as is + elseif isdatetime(col) + continue; + + else + % Catch-all for unexpected types, convert to string + cleanedTable.(varNames{i}) = string(col); + end + end +end + +function out = convertStructToString(x) +% Helper function to convert struct NaN to string 'NaN' + if isstruct(x) + out = "NaN"; + elseif isstring(x) || ischar(x) + out = string(x); + else + out = x; + end +end diff --git a/projects/ECOC_2025/submit_dsp.m b/projects/ECOC_2025/submit_dsp.m new file mode 100644 index 0000000..019056d --- /dev/null +++ b/projects/ECOC_2025/submit_dsp.m @@ -0,0 +1,54 @@ +function [output, future] = submit_dsp(run_id, basePath, database_name, savePath, options) +% runDSPInBackground Runs or submits DSP processing based on 'parallel' flag. +% +% [output, future] = runDSPInBackground(..., options) +% - output: result from dsp_run_id (only immediately available in serial mode) +% - future: FevalFuture object if run in parallel, otherwise [] + +arguments + run_id + basePath + database_name + savePath + options.parallel (1,1) logical = true + options.max_occurences = 1; +end + +if options.parallel + % Check if a pool exists + pool = gcp('nocreate'); + if isempty(pool) + parpool; + end + + % Submit the DSP function asynchronously + future = parfeval( ... + @dsp_run_id, 1, ... % One output + run_id, ... + "database_path", basePath, ... + "database_name", database_name, ... + 'storage_path', savePath, ... + 'append_to_db', 1, ... + 'max_occurences', options.max_occurences ... + ); + + output = []; + fprintf('DSP task for run_id %d submitted to the pool.\n', run_id); + + +else + % Run synchronously (debug mode), capture output + fprintf('Running DSP task for run_id %d in main thread (debug mode).\n', run_id); + + output = dsp_run_id( ... + run_id, ... + "database_path", basePath, ... + "database_name", database_name, ... + 'storage_path', savePath, ... + 'append_to_db', 1, ... + 'max_occurences', options.max_occurences ... + ); + + future = []; % No future since it's synchronous +end +end