diff --git a/Classes/04_DSP/Equalizer/FFE.m b/Classes/04_DSP/Equalizer/FFE.m index 18442d6..849539e 100644 --- a/Classes/04_DSP/Equalizer/FFE.m +++ b/Classes/04_DSP/Equalizer/FFE.m @@ -15,7 +15,7 @@ classdef FFE < handle e e_tr error - debug_struct + len_tr mu_tr @@ -31,6 +31,13 @@ classdef FFE < handle constellation % symbol constellation decide %wether to return the (hard) decisions or the result after FFE (soft) + + save_debug = 0; + debug_struct + + optmize_mus = 0; + mu_optimization + mu_optimization_iter = 0; end methods @@ -51,6 +58,9 @@ classdef FFE < handle options.decide = false; + options.save_debug = 0; + options.optmize_mus = 0; + end fn = fieldnames(options); @@ -75,6 +85,12 @@ classdef FFE < handle delta = 0.05; obj.P = (1/delta) * eye(obj.order); + if obj.optmize_mus + obj.optimizeMus(X.signal,D.signal); + obj.e = zeros(obj.order,1); + obj.P = (1/delta) * eye(obj.order); + end + % Training Mode training = 1; showviz = 0; @@ -131,7 +147,6 @@ classdef FFE < handle mask = ones(obj.order,1); maincursor_pos=ceil(length(obj.e)/2); always_ideal_decision = 0; - save_debug = 0; grad =0; weight = 0; update = 0; @@ -205,21 +220,88 @@ classdef FFE < handle end - if save_debug + if obj.save_debug obj.debug_struct.error(epoch,symbol) = err(symbol) * err(symbol)'; % Instantaneous square error + if training obj.debug_struct.error_tr(epoch,symbol) = err(symbol) * err(symbol)'; % Instantaneous square error obj.debug_struct.update_tr(epoch,symbol) = update.'*update ./ rms(obj.e); + else + % if symbol == length(d) + % M = numel(unique(d)); + % bits_ref = PAMmapper(M,0).demap(d); + % bits_eq = PAMmapper(M,0).demap(d_hat); + % [~,~,obj.debug_struct.ber(epoch)] = calc_ber(bits_eq,bits_ref); + % end end - obj.debug_struct.main_cursor(epoch,symbol) = abs(obj.e(maincursor_pos)); - obj.debug_struct.mu_nlms(epoch,symbol) = weight; - obj.debug_struct.update_gradient(epoch,symbol) = grad.'*grad; - obj.debug_struct.update(epoch,symbol) = update.'*update ./ rms(obj.e); + + % obj.debug_struct.main_cursor(epoch,symbol) = abs(obj.e(maincursor_pos)); + % obj.debug_struct.mu_nlms(epoch,symbol) = weight; + % obj.debug_struct.update_gradient(epoch,symbol) = grad.'*grad; + % obj.debug_struct.update(epoch,symbol) = update.'*update ./ rms(obj.e); + end end end end + + function optimizeMus(obj,x,d) + switch obj.adaption_technique + case adaption_method.lms + mu_range = [1e-5, 1e-2]; + case adaption_method.nlms + mu_range = [1e-3, 5e-1]; + case adaption_method.rls + mu_range = [0.98, 0.99999]; + end + + mu_tr_var = optimizableVariable("mu_tr",mu_range,"Transform","log"); + vars = mu_tr_var; + if obj.dd_mode + vars = [vars, optimizableVariable("mu_dd",mu_range,"Transform","log")]; + end + obj.mu_optimization_iter = 0; + obj.mu_optimization = bayesopt(@(p)obj.muObjective(p,x,d),vars, ... + "MaxObjectiveEvaluations",10, ... + "AcquisitionFunctionName","expected-improvement-plus", ... + "IsObjectiveDeterministic",false, ... + "Verbose",0, ... + "PlotFcn",[]); + obj.mu_tr = obj.mu_optimization.XAtMinObjective.mu_tr; + if obj.dd_mode + obj.mu_dd = obj.mu_optimization.XAtMinObjective.mu_dd; + end + fprintf("\nFFE mu opt done: mu_tr=%9.3e, mu_dd=%9.3e, metric=%9.3e\n", ... + obj.mu_tr,obj.mu_dd,obj.mu_optimization.MinObjective); + end + + function objective = muObjective(obj,params,x,d) + old_debug = obj.save_debug; + obj.save_debug = 1; + obj.e = zeros(obj.order,1); + obj.P = (1/0.05) * eye(obj.order); + obj.debug_struct = struct(); + obj.equalize(x,d,params.mu_tr,obj.epochs_tr,obj.len_tr,1,0); + if obj.dd_mode + obj.equalize(x,d,params.mu_dd,obj.epochs_dd,numel(x),0,0); + objective = mean(obj.debug_struct.error(end,:),"omitnan"); + else + objective = mean(obj.debug_struct.error_tr(end,:),"omitnan"); + end + if ~isfinite(objective) + objective = inf; + end + obj.mu_optimization_iter = obj.mu_optimization_iter + 1; + if obj.dd_mode + fprintf("\rFFE mu opt %02d: mu_tr=%9.3e, mu_dd=%9.3e, metric=%9.3e", ... + obj.mu_optimization_iter,params.mu_tr,params.mu_dd,objective); + else + fprintf("\rFFE mu opt %02d: mu_tr=%9.3e, metric=%9.3e", ... + obj.mu_optimization_iter,params.mu_tr,objective); + end + obj.save_debug = old_debug; + end end -end \ No newline at end of file +end diff --git a/Classes/DataBaseHandler/DBHandler.m b/Classes/DataBaseHandler/DBHandler.m index 4ba8138..fad15e4 100644 --- a/Classes/DataBaseHandler/DBHandler.m +++ b/Classes/DataBaseHandler/DBHandler.m @@ -54,7 +54,7 @@ classdef DBHandler < handle options.password, ... % Password (or getSecret) "Vendor", "MySQL", ... "Server", options.server, ... - "PortNumber", 3306, ... + "PortNumber", options.port, ... "JDBCDriverLocation", "C:\Users\Silas\Documents\mysql-connector-j-9.3.0\mysql-connector-j-9.3.0.jar"); end diff --git a/Functions/EQ_structures/dsp_runid.m b/Functions/EQ_structures/dsp_runid.m index ea67304..39526d2 100644 --- a/Functions/EQ_structures/dsp_runid.m +++ b/Functions/EQ_structures/dsp_runid.m @@ -7,6 +7,10 @@ arguments options.parameters = struct(); options.database_type options.dataBase + options.server = "134.245.243.254"; + options.port = 3306; + options.user = "silas"; + options.password = "silas"; options.load_file_path = struct(); options.storage_path options.mode @@ -18,7 +22,8 @@ try if options.mode == "load_run_id" || options.append_to_db database = DBHandler("dataBase", [options.dataBase], "type", options.database_type, ... - 'user', "silas", 'password', "silas", 'server', "134.245.243.254"); + "user", options.user, "password", options.password, ... + "server", options.server, "port", options.port); end if options.mode == "load_run_id" @@ -68,7 +73,6 @@ function output = initializeOutput() end function dspInput = loadDspInputFromFiles(run_id, options) - %#ok Tx_bits = load(options.load_file_path.tx_bits_path); Symbols = load(options.load_file_path.tx_symbols_path); Scpe_sig_raw = load(options.load_file_path.rx_raw_path); diff --git a/Functions/EQ_structures/ffe.m b/Functions/EQ_structures/ffe.m index 05acd96..a58e215 100644 --- a/Functions/EQ_structures/ffe.m +++ b/Functions/EQ_structures/ffe.m @@ -82,6 +82,9 @@ try eq_.b2 = []; eq_.b3 = []; end +try + eq_.mu_optimization = []; +end ffe_results.config = Equalizerstruct(); ffe_results.config.eq = jsonencode(eq_); diff --git a/Functions/EQ_structures/vnle_postfilter_mlse.m b/Functions/EQ_structures/vnle_postfilter_mlse.m index 0e6aea5..803ec70 100644 --- a/Functions/EQ_structures/vnle_postfilter_mlse.m +++ b/Functions/EQ_structures/vnle_postfilter_mlse.m @@ -274,32 +274,32 @@ showEQNoisePSD(eq_noise, "fignum", 338, "displayname", 'Residual Noise after VNL for t = 1:4 pf_.ncoeff = t; [~,~] = pf_.process(eq_signal_sd, eq_noise); - showEQNoisePSD(eq_noise, "fignum", 3388, "displayname", 'Residual Noise after VNLE', 'postfilter_taps', pf_.coefficients); + showEQNoisePSD(eq_noise, "fignum", 339, "displayname", 'Residual Noise after VNLE', 'postfilter_taps', pf_.coefficients); end -tx_symbols.spectrum("displayname",'Equalized Signal','fignum',1234,'normalizeTo0dB',1); +tx_symbols.spectrum("displayname",'Equalized Signal','fignum',340,'normalizeTo0dB',1); if ~isempty(postFFE) - showEQcoefficients('n1', postFFE.e, "displayname", 'Coefficients', 'fignum', 338); + showEQcoefficients('n1', postFFE.e, "displayname", 'Coefficients', 'fignum', 341); end % showEQcoefficients('n1', eq_.e1,'n2', eq_.e2,'n3', eq_.e3, "displayname", 'Coefficients', 'fignum', 339); showEQfilter(eq_.e, eq_signal_sd.fs.*2); figure(340); clf; -eq_signal_sd.eye(eq_signal_sd.fs, M, "fignum", 340); +eq_signal_sd.eye(eq_signal_sd.fs, M, "fignum", 342); figure(341); clf; -showLevelHistogram(eq_signal_sd, tx_symbols, "fignum", 341); +showLevelHistogram(eq_signal_sd, tx_symbols, "fignum", 343); warning off -showLevelScatter(eq_signal_sd, tx_symbols, "fignum", 400); +showLevelScatter(eq_signal_sd, tx_symbols, "fignum", 344); % showLevelScatter(rx_signal.resample("fs_out", tx_symbols.fs), tx_symbols, "fignum", 401); drawnow; warning on -whitened_noise.spectrum("displayname",'after postfilter','fignum',342); -eq_noise.spectrum("displayname",'before postfilter','fignum',342); +whitened_noise.spectrum("displayname",'after postfilter','fignum',345); +eq_noise.spectrum("displayname",'before postfilter','fignum',345); end \ No newline at end of file diff --git a/Functions/Job_Processing/submitJobs.m b/Functions/Job_Processing/submitJobs.m index 26fe083..a11f057 100644 --- a/Functions/Job_Processing/submitJobs.m +++ b/Functions/Job_Processing/submitJobs.m @@ -36,17 +36,10 @@ for r = 1:nRunIds for k = 1:nJobsPerRunId jobCounter = jobCounter + 1; optionalVars = buildOptionalVars(k, submit_options.wh); + jobOptions = dsp_options; + jobOptions.parameters = optionalVars; - jobs(jobCounter).args = { ... - run_ids(r), ... - "database_type", dsp_options.database_type, ... - "dataBase", dsp_options.dataBase, ... - "append_to_db", dsp_options.append_to_db, ... - "load_file_path", dsp_options.load_file_path, ... - "max_occurences", dsp_options.max_occurences, ... - "storage_path", dsp_options.storage_path, ... - "mode", dsp_options.mode, ... - "parameters", optionalVars}; + jobs(jobCounter).args = [{run_ids(r)}, structToNameValue(jobOptions)]; jobs(jobCounter).label = sprintf('RunID %d, Job %d', run_ids(r), k); jobs(jobCounter).meta.runIndex = r; jobs(jobCounter).meta.jobIndex = k; @@ -91,6 +84,16 @@ end end end + function nameValue = structToNameValue(options) + names = fieldnames(options); + nameValue = cell(1, 2*numel(names)); + + for i = 1:numel(names) + nameValue{2*i - 1} = string(names{i}); + nameValue{2*i} = options.(names{i}); + end + end + function storeResult(val, job, ~) if ~isempty(wh) jobIndex = job.meta.jobIndex; diff --git a/Theory/Signal_Processing/EQ_adaption.m b/Theory/Signal_Processing/EQ_adaption.m new file mode 100644 index 0000000..34eebe6 --- /dev/null +++ b/Theory/Signal_Processing/EQ_adaption.m @@ -0,0 +1,342 @@ +% minimal example IM/DD + +M = 4; +fsym = 180e9; + +apply_pulsef = 1; +fdac = 256e9; +fadc = 256e9; +random_key = 1; + +rcalpha = 0.05; +kover = 16; + +duob_mode = db_mode.no_db; + +vbias_rel = 0.5; +u_pi = 3; +vbias = -vbias_rel*u_pi; + +laser_wavelength = 1293; +laser_linewidth = 0; +tx_bw_nyquist = 0.8; + +% Channel +link_length = 1; + +% RX +rop = -7; +rx_bw_nyquist = 0.8; + +vnle_order1 = 50; +vnle_order2 = 7; +vnle_order3 = 7; + +vnle_order=[vnle_order1,vnle_order2,vnle_order3]; +dfe_order = [0 0 0]; + +pf_ncoeffs = 1; + +alpha = 0; + +len_tr = 4096*2; + +mu_ffe1 = 0.0001; +mu_ffe2 = 0.0008; +mu_ffe3 = 0.001; +mu_dc = 0.005; +% mu_dc = 0; + +mu_ffe = [mu_ffe1 mu_ffe3 mu_ffe3]; +mu_dfe = 0.0004; + +n_stat_runs = 1; +base_random_key = random_key; +output.eq_stats = struct(); + +for stat_run = 1:n_stat_runs +random_key = base_random_key + stat_run - 1; +Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"alpha",rcalpha); + +[Digi_sig,Symbols,Tx_bits] = PAMsource(... + "fsym",fsym,"M",M,"order",16,"useprbs",0,... + "fs_out",fdac,... + "applyclipping",0,"clipfactor",1.5,... + "applypulseform",apply_pulsef,"pulseformer",Pform,... + "randkey",random_key,... + 'duobinary_mode',duob_mode,... + "mrds_code",0,"mrds_blocklength",512).process(); + +%%%%% 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); +El_sig.spectrum("displayname",'Digi Spectrum','fignum',1,'normalizeTo0dB',1); +xlim([0,130]); +ylim([-30,5]); +% El_sig = El_sig.setPower(0,"dBm"); + +%%%%% Electrical Driver Amplifier %%%%%% +% 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,"alpha",alpha).process(El_sig); + +Opt_sig.spectrum("displayname",'Opt Spectrum','fignum',10,'normalizeTo0dB',1); + +% Opt_sig.eye(fsym,M,"displayname",'eye adter modulator','fignum',2026); + +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); + +%%%%%% PD Square Law %%%%%% +Rx_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20,"nep",1.8e-11).process(Rx_sig); + +%%%%%% Low-pass RX (PD, El. Connectors and Scope %%%%%% +rx_bwl = 80e9; +Rx_sig = Filter('filtdegree',4,"f_cutoff",rx_bwl,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(Rx_sig); + +% %%%%%% Low-pass Scope %%%%%% +Lp_scpe = Filter('filtdegree',4,"f_cutoff",110e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true); + +%%%%%% Scope %%%%%% +Scpe_sig = Scope("fsimu",fdac*kover,"fadc",fadc,... + "delay",0,"fixed_delay",0,"filtertype",filtertypes.butterworth,... + "samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,... + "adcresolution",8,"quantbuffer",0.1,'block_dc',1,'lpf_active',1,'H_lpf',Lp_scpe).process(Rx_sig); +%% + +% 1) matched filter +% pulse is symmetric, hence we can use pulsef firectly as matched filter. +% It feels off (bit I think correct) that the fsym is now the output freq.!! +% -> output 2 sps to omit timing recovery!? +Pform = Pulseformer("fsym",fsym,"fdac",2*fsym,"pulse","rrc","pulselength",16,"alpha",rcalpha,"matched",1); +Scpe_sig = Pform.process(Scpe_sig); +Scpe_sig.spectrum("displayname",'Signal after matched filter','fignum',1,'normalizeTo0dB',1); +% + +% %% +% %%%%%% 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",0); +Scpe_sig.spectrum("displayname",'Opt Spectrum','fignum',11,'normalizeTo0dB',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); +Scpe_sig.signal = Scpe_sig.signal(1:2*length(Symbols)); + +%% -------------------- FFE -------------------- +ffe_options = struct( ... + "epochs_tr", 3, ... + "epochs_dd", 3, ... + "len_tr", 4096, ... + "order", 40, ... + "sps", 2, ... + "decide", 0, ... + "adaption", adaption_method.rls, ... + "dd_mode", 1); + +eq_nlms = FFE("epochs_tr",ffe_options.epochs_tr,"epochs_dd",ffe_options.epochs_dd, ... + "len_tr",ffe_options.len_tr,"mu_dd",0.001,"mu_tr",0.001, ... + "order",ffe_options.order,"sps",ffe_options.sps,"decide",ffe_options.decide, ... + "adaption",ffe_options.adaption,"dd_mode",ffe_options.dd_mode,"save_debug",1,"optmize_mus",1); + +output.ffe_results = ffe(eq_nlms,M,Scpe_sig,Symbols,Tx_bits, ... + "precode_mode",duob_mode,'showAnalysis',0,"postFFE",[], ... + "eth_style_symbol_mapping",0); + +output.ffe_results.metrics.print("description","FFE"); + +%% + +% -------------------- FFE -------------------- +%% + +for i = 3:-1:1 + + mu_adaptation = adaption_method(i); + + switch mu_adaptation + case adaption_method.lms + mu_search_range = [1e-5, 1e-2]; + + case adaption_method.nlms + mu_search_range = [1e-3, 5e-1]; + + case adaption_method.rls + mu_search_range = [0.98, 0.99999]; % here mu is lambda + end + mu_max_evaluations = 10; + + ffe_options = struct( ... + "epochs_tr", 3, ... + "epochs_dd", 2, ... + "len_tr", 4096*2, ... + "order", 20, ... + "sps", 2, ... + "decide", 0, ... + "adaption", mu_adaptation, ... + "dd_mode", 1); + + method_name = char(mu_adaptation); + rng(random_key, "twister"); + mu_dd_var = optimizableVariable("mu_dd", mu_search_range, "Transform", "log"); + mu_tr_var = optimizableVariable("mu_tr", mu_search_range, "Transform", "log"); + global ffe_mu_trace + ffe_mu_trace = table('Size',[0 5], ... + 'VariableTypes',{'double','double','double','double','double'}, ... + 'VariableNames',{'mu_dd','mu_tr','BER','SNR','mse_last_epoch'}); + objective_fn = @(params) ffe_mu_objective(params, ffe_options, M, Scpe_sig, Symbols, Tx_bits, duob_mode); + + output.eq_stats.(method_name).mu_optimization{stat_run} = bayesopt(objective_fn, [mu_dd_var, mu_tr_var], ... + "MaxObjectiveEvaluations", mu_max_evaluations, ... + "AcquisitionFunctionName", "expected-improvement-plus", ... + "IsObjectiveDeterministic", false, ... + "Verbose", 1, ... + "PlotFcn", []); + + output.eq_stats.(method_name).mu_observed(stat_run,:) = output.eq_stats.(method_name).mu_optimization{stat_run}.XAtMinObjective; + output.eq_stats.(method_name).mu_estimated(stat_run,:) = output.eq_stats.(method_name).mu_optimization{stat_run}.XAtMinEstimatedObjective; + output.eq_stats.(method_name).optimization_trace{stat_run} = ffe_mu_trace; + + best_mu_dd = output.eq_stats.(method_name).mu_observed.mu_dd(stat_run); + best_mu_tr = output.eq_stats.(method_name).mu_observed.mu_tr(stat_run); + output.eq_stats.(method_name).mu_used(stat_run,:) = table(best_mu_dd,best_mu_tr, ... + 'VariableNames',{'mu_dd','mu_tr'}); + + eq_nlms = FFE("epochs_tr",ffe_options.epochs_tr,"epochs_dd",ffe_options.epochs_dd, ... + "len_tr",ffe_options.len_tr,"mu_dd",best_mu_dd,"mu_tr",best_mu_tr, ... + "order",ffe_options.order,"sps",ffe_options.sps,"decide",ffe_options.decide, ... + "adaption",ffe_options.adaption,"dd_mode",ffe_options.dd_mode,"save_debug",1); + + output.ffe_results = ffe(eq_nlms,M,Scpe_sig,Symbols,Tx_bits, ... + "precode_mode",duob_mode,'showAnalysis',0,"postFFE",[], ... + "eth_style_symbol_mapping",0); + + output.ffe_results.metrics.print("description",sprintf('%s',mu_adaptation)); + + output.eq_stats.(method_name).error_tr_first_epoch(stat_run,:) = eq_nlms.debug_struct.error_tr(1,:); + + % figure(6) + % hold on + % for epoch = 1:ffe_options.epochs_tr + % plot(1:size(eq_nlms.debug_struct.error_tr,2),movmean(eq_nlms.debug_struct.error_tr(epoch,:),10)); + % end + +end + +end +%% + +figure(5); +clf; +hold on; +method_names = fieldnames(output.eq_stats); +cols = cbrewer2('Paired',6); +cols = cols(2:2:end,:); +for i = 1:numel(method_names) + method_name = method_names{i}; + first_epoch = output.eq_stats.(method_name).error_tr_first_epoch; + + + for j = 1:size(first_epoch,1) + plot(movmean(output.eq_stats.(method_name).error_tr_first_epoch(j,:),10), ... + "DisplayName",method_name,"Color",[0.7,0.7,0.7],"LineWidth",0.1,'HandleVisibility','off'); + end + +end + +for i = 1:numel(method_names) + method_name = method_names{i}; + first_epoch = output.eq_stats.(method_name).error_tr_first_epoch; + + output.eq_stats.(method_name).error_tr_first_epoch_mean = mean(first_epoch,1,"omitnan"); + output.eq_stats.(method_name).error_tr_first_epoch_std = std(first_epoch,0,1,"omitnan"); + plot(movmean(output.eq_stats.(method_name).error_tr_first_epoch_mean,50), ... + "DisplayName",method_name,"Color",cols(i,:)); +end + + +grid on; +xlabel("Training symbol index","Interpreter","tex"); +ylabel("Mean squared error, first epoch","Interpreter","tex"); +legend("Location","best","Interpreter","tex"); + +figure(6); +clf; +tiledlayout(numel(method_names),3,"TileSpacing","compact","Padding","compact"); +metric_names = {'BER','SNR','mse_last_epoch'}; +metric_labels = {'BER','SNR','mean abs(error)^2, last epoch'}; +for i = 1:numel(method_names) + method_name = method_names{i}; + trace = vertcat(output.eq_stats.(method_name).optimization_trace{:}); + observed_mu = output.eq_stats.(method_name).mu_observed; + estimated_mu = output.eq_stats.(method_name).mu_estimated; + for metric_idx = 1:numel(metric_names) + nexttile; + metric = trace.(metric_names{metric_idx}); + if metric_idx ~= 2 + metric = log10(max(metric,realmin)); + end + scatter(trace.mu_tr,trace.mu_dd,25,metric,"filled"); + hold on; + plot(observed_mu.mu_tr,observed_mu.mu_dd,"kp","MarkerSize",11, ... + "MarkerFaceColor","y","LineWidth",1.2,"DisplayName","Observed best"); + plot(estimated_mu.mu_tr,estimated_mu.mu_dd,"kx","MarkerSize",10, ... + "LineWidth",1.6,"DisplayName","Estimated best"); + set(gca,"XScale","log","YScale","log"); + grid on; + xlabel('\mu_{tr}',"Interpreter","tex"); + ylabel('\mu_{dd}',"Interpreter","tex"); + title(sprintf("%s: %s",method_name,metric_labels{metric_idx}),"Interpreter","tex"); + cb = colorbar; + if metric_idx ~= 2 + cb.Label.String = sprintf("log10(%s)",metric_labels{metric_idx}); + else + cb.Label.String = metric_labels{metric_idx}; + end + cb.Label.Interpreter = "tex"; + legend("Location","best","Interpreter","tex"); + end +end + +function objective = ffe_mu_objective(params, ffe_options, M, rx_signal, tx_symbols, tx_bits, duob_mode) + global ffe_mu_trace + try + eq_candidate = FFE("epochs_tr",ffe_options.epochs_tr,"epochs_dd",ffe_options.epochs_dd, ... + "len_tr",ffe_options.len_tr,"mu_dd",params.mu_dd,"mu_tr",params.mu_tr, ... + "order",ffe_options.order,"sps",ffe_options.sps,"decide",ffe_options.decide, ... + "adaption",ffe_options.adaption,"dd_mode",ffe_options.dd_mode,"save_debug",1); + + candidate_results = ffe(eq_candidate, M, rx_signal, tx_symbols, tx_bits, ... + "precode_mode",duob_mode,'showAnalysis',0,"postFFE",[], ... + "eth_style_symbol_mapping",0); + + objective = candidate_results.metrics.BER; + if ~isfinite(objective) + objective = 0.5; + end + mse_last_epoch = mean(eq_candidate.debug_struct.error_tr(end,:),"omitnan"); + ffe_mu_trace = [ffe_mu_trace; table(params.mu_dd,params.mu_tr,objective, ... + candidate_results.metrics.SNR,mse_last_epoch, ... + 'VariableNames',{'mu_dd','mu_tr','BER','SNR','mse_last_epoch'})]; + + fprintf("mu_dd=%9.3e, mu_tr=%9.3e -> BER=%9.3e, SNR=%6.2f dB, MSE=%9.3e\n", ... + params.mu_dd, params.mu_tr, objective, candidate_results.metrics.SNR, mse_last_epoch); + catch ME + objective = 0.5; + ffe_mu_trace = [ffe_mu_trace; table(params.mu_dd,params.mu_tr,objective,NaN,NaN, ... + 'VariableNames',{'mu_dd','mu_tr','BER','SNR','mse_last_epoch'})]; + fprintf("mu_dd=%9.3e, mu_tr=%9.3e -> failed: %s\n", ... + params.mu_dd, params.mu_tr, ME.message); + end +end diff --git a/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/FIGURE_Spectra.m b/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/FIGURE_Spectra.m index e819327..7762768 100644 --- a/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/FIGURE_Spectra.m +++ b/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/FIGURE_Spectra.m @@ -1,6 +1,17 @@ dsp_options.storage_path = 'Z:\2024\sioe_labor\'; dsp_options.max_occurences = 1; -database = DBHandler("dataBase", 'labor_highspeed', "type", 'mysql' ); +% database = DBHandler("dataBase", 'labor_highspeed', "type", 'mysql' ); +dsp_options.database_type = "mysql"; +dsp_options.dataBase = 'labor_highspeed';%'C:\Users\Silas\Documents\MATLAB\Datensätze\sioe_labor\silas_labor_newdsp_newstructure.db'; +dsp_options.server = "192.168.178.192";% "134.245.243.254"; +dsp_options.user = "silas"; +dsp_options.password = "silas"; +dsp_options.storage_path = 'W:\labdata\sioe_labor\'; +database = DBHandler("dataBase", [dsp_options.dataBase],... + "type", dsp_options.database_type,... + "server", dsp_options.server,... + "user", dsp_options.user, "password", dsp_options.password); + rates = [420e9]; cols = cbrewer2('BuPu',25); @@ -10,10 +21,10 @@ cols = cbrewer2('Set1',6); fignum = 200; fig=figure(fignum);clf; -for dbmode = 1%length(rates) +for dbmode = 0%length(rates) - if 0 + if 1 rcalpha = 0.05; fsym = rates/2; pulsef = 1; @@ -41,7 +52,7 @@ for dbmode = 1%length(rates) maxamp = -37; precomp_est = ChannelFreqResp("Nacq",2048,"Navg",100,"Ncp",63,'f_ref',Digi_sig.fs); - precomp_path = "C:\Users\Silas\Documents\MATLAB\Datensätze\sioe_labor\precomp"; + precomp_path = "W:\labdata\sioe_labor\precomp"; precomp_fn = "lab_high_speed"; Digi_sig_pre = precomp_est.precomp(Digi_sig,'maxampdb',maxamp,'loadPath',precomp_path,'fileName',precomp_fn); diff --git a/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/generate_spectrum_plots.m b/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/generate_spectrum_plots.m index f90cf1a..9fd31ed 100644 --- a/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/generate_spectrum_plots.m +++ b/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/generate_spectrum_plots.m @@ -1,7 +1,16 @@ dsp_options.storage_path = 'Z:\2024\sioe_labor\'; dsp_options.max_occurences = 1; -database = DBHandler("dataBase", 'labor_highspeed', "type", 'mysql' ); - +% database = DBHandler("dataBase", 'labor_highspeed', "type", 'mysql' ); + dsp_options.database_type = "mysql"; + dsp_options.dataBase = 'labor_highspeed';%'C:\Users\Silas\Documents\MATLAB\Datensätze\sioe_labor\silas_labor_newdsp_newstructure.db'; + dsp_options.server = "192.168.178.192";% "134.245.243.254"; + dsp_options.user = "silas"; + dsp_options.password = "silas"; + dsp_options.storage_path = 'W:\labdata\sioe_labor\'; + database = DBHandler("dataBase", [dsp_options.dataBase],... + "type", dsp_options.database_type,... + "server", dsp_options.server,... + "user", dsp_options.user, "password", dsp_options.password); rate = 390e9; %% 1 - PAM 4 with preemphasis @@ -14,7 +23,7 @@ fp.where('Runs', 'wavelength','EQUALS', 1310); fp.where('Runs', 'db_mode','EQUALS', 0); fp.where('Runs', 'rop_attenuation','EQUAL', 0); -[dataTable,~] = db.queryDB(fp, database.getTableFieldNames('Runs')); +[dataTable,~] = database.queryDB(fp, database.getTableFieldNames('Runs')); dataTable = queryRunid(dataTable.run_id, database); fsym = dataTable.symbolrate; @@ -46,7 +55,7 @@ fp.where('Runs', 'wavelength','EQUALS', 1310); fp.where('Runs', 'db_mode','EQUALS', 1); fp.where('Runs', 'rop_attenuation','EQUAL', 0); -[dataTable,~] = db.queryDB(fp, database.getTableFieldNames('Runs')); +[dataTable,~] = database.queryDB(fp, database.getTableFieldNames('Runs')); dataTable = queryRunid(dataTable.run_id, database); fsym = dataTable.symbolrate; diff --git a/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/run_dsp_from_db.m b/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/run_dsp_from_db.m index 730c410..e11693f 100644 --- a/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/run_dsp_from_db.m +++ b/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/run_dsp_from_db.m @@ -12,11 +12,14 @@ if dsp_options.mode == "load_run_id" dsp_options.database_type = "mysql"; dsp_options.dataBase = 'labor_highspeed';%'C:\Users\Silas\Documents\MATLAB\Datensätze\sioe_labor\silas_labor_newdsp_newstructure.db'; + dsp_options.server = "192.168.178.192";% "134.245.243.254"; + dsp_options.user = "silas"; + dsp_options.password = "silas"; dsp_options.storage_path = 'W:\labdata\sioe_labor\'; db = DBHandler("dataBase", [dsp_options.dataBase],... "type", dsp_options.database_type,... - "server","192.168.178.192",...% "server","134.245.243.254",... - "user","silas","password","silas"); + "server", dsp_options.server,... + "user", dsp_options.user, "password", dsp_options.password); elseif experiment == "mpi_ecoc_2025" @@ -379,4 +382,4 @@ ylabel('AIR in Gbps'); set(gca, 'XTick', xticks_vals(1:2:end), 'XTickLabel', xtick_labels(1:2:end)); grid on; legend('Location','best'); -% xlim([1, 256]) \ No newline at end of file +% xlim([1, 256]) diff --git a/projects/IMDD_base_system/minimal_example.m b/projects/IMDD_base_system/minimal_example.m index 919cd7e..4a2d867 100644 --- a/projects/IMDD_base_system/minimal_example.m +++ b/projects/IMDD_base_system/minimal_example.m @@ -1,7 +1,7 @@ % minimal example IM/DD M = 4; -fsym = 112e9; +fsym = 180e9; apply_pulsef = 1; fdac = 256e9; @@ -63,8 +63,8 @@ Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"al "mrds_code",0,"mrds_blocklength",512).process(); %%%%% 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); +El_sig = M8199B("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); El_sig.spectrum("displayname",'Digi Spectrum','fignum',1,'normalizeTo0dB',1); xlim([0,130]); ylim([-30,5]); @@ -130,24 +130,24 @@ Scpe_sig = Scpe_sig - mean(Scpe_sig.signal); Scpe_sig.signal = Scpe_sig.signal(1:2*length(Symbols)); %% - -% -------------------- FFE -------------------- -ffe_order = [50, 0, 0]; -eq_ = 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); - -% eq_ = FFE("epochs_tr",4,"epochs_dd",5,"len_tr",4096,"mu_dd",0.01,"mu_tr",0.01,"order",50,"sps",2,"decide",0, "adaption",adaption_method.nlms,"dd_mode",1); -eq_ = FFE_DFE("epochs_tr",5,"epochs_dd",5,"len_tr",512,"ffe_mu_dd",1e-4,"dfe_mu_dd",5e-4,"ffe_mu_tr",0,"dfe_mu_tr",0,"ffe_order",99,"dfe_order",99,"sps",2,"decide",0); - - -output.ffe_results = ffe(eq_,M,Scpe_sig,Symbols,Tx_bits, ... - "precode_mode",duob_mode,'showAnalysis',1,"postFFE",[], ... - "eth_style_symbol_mapping",0); - -output.ffe_results.metrics.print("description",'DFE'); - +if 0 + % -------------------- FFE -------------------- + ffe_order = [50, 0, 0]; + eq_ = 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); + + % eq_ = FFE("epochs_tr",4,"epochs_dd",5,"len_tr",4096,"mu_dd",0.01,"mu_tr",0.01,"order",50,"sps",2,"decide",0, "adaption",adaption_method.nlms,"dd_mode",1); + eq_ = FFE_DFE("epochs_tr",5,"epochs_dd",5,"len_tr",512,"ffe_mu_dd",1e-4,"dfe_mu_dd",5e-4,"ffe_mu_tr",0,"dfe_mu_tr",0,"ffe_order",99,"dfe_order",99,"sps",2,"decide",0); + + + output.ffe_results = ffe(eq_,M,Scpe_sig,Symbols,Tx_bits, ... + "precode_mode",duob_mode,'showAnalysis',1,"postFFE",[], ... + "eth_style_symbol_mapping",0); + + output.ffe_results.metrics.print("description",'DFE'); +end %% % -------------------- VNLE + MLSE -------------------- diff --git a/projects/Messung_Zürich/dsp_all_files.m b/projects/Messung_Zürich/dsp_all_files.m index 9ec98ab..1b9c1f4 100644 --- a/projects/Messung_Zürich/dsp_all_files.m +++ b/projects/Messung_Zürich/dsp_all_files.m @@ -1,5 +1,5 @@ -M = 6; +M = 4; if M == 2 file_codes = {'20250221T035221' '20250221T035354' '20250221T035824' '20250221T035931' '20250221T040035' '20250221T040132' '20250221T040226' '20250221T040523' '20250221T040646' '20250221T040723' '20250221T040843' '20250221T041011' '20250221T041101' '20250221T041244'}; elseif M == 4 @@ -10,7 +10,7 @@ elseif M ==8 file_codes = {'20250221T004926' '20250221T023534' '20250221T024256' '20250221T024629' '20250221T024929' '20250221T025305' '20250221T025505' '20250221T025856' '20250221T030122' '20250221T030311' '20250221T030513'}; end -if 0 +if 1 uloops = struct; uloops.filecode = file_codes; uloops.vnle_order1 = [50]; @@ -19,7 +19,7 @@ if 0 wh = DataStorage(uloops); wh.addStorage("ber"); - wh = submit_handle(@dsp_ief_file,wh,"parallel",1); + wh = submit_handle(@dsp_ief_file,wh,"parallel",0); end if 0 diff --git a/projects/Messung_Zürich/dsp_ief_file.m b/projects/Messung_Zürich/dsp_ief_file.m index 3fd0754..fc29b8a 100644 --- a/projects/Messung_Zürich/dsp_ief_file.m +++ b/projects/Messung_Zürich/dsp_ief_file.m @@ -36,9 +36,11 @@ function [results] = dsp_ief_file(varargin) end pathToMeasurement = "C:\Users\Silas\Documents\MATLAB\Datensätze\IEF_Polariton_2025\36_IMDD_Kiel\Data\20250221"; + pathToMeasurement = "W:\labdata\Plasmonic O-band Modulator (ETH Zurich)\Data\20250221"; filename = findFileByCode(pathToMeasurement, filecode{1}); pathToTimingRecov = "C:\Users\Silas\Documents\MATLAB\Datensätze\IEF_Polariton_2025\36_IMDD_Kiel\Data\TR_ZIP"; + pathToTimingRecov = "W:\labdata\Plasmonic O-band Modulator (ETH Zurich)\Data\TR_ZIP"; filename_2 = findFileByCode(pathToTimingRecov, filecode{1}); % Extract parameters from the filename using an updated regex @@ -166,7 +168,7 @@ function [results] = dsp_ief_file(varargin) %%%%% VNLE + PF + MLSE %%%% if 1 - eq_post = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",1001,"sps",1,"decide",0); + eq_post = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",1.933e-04,"order",1001,"sps",1,"decide",0,"optmize_mus",1,"dd_mode",0); pf_ = Postfilter("ncoeff",1,"useBurg",1); mlse_ = MLSE("duobinary_output",0,'M',config.M,'trellis_states',PAMmapper(config.M,0).levels); mlse_ = MLSE_viterbi("duobinary_output",0,'M',config.M,'trellis_states',PAMmapper(config.M,0).levels);