diff --git a/Classes/00_signals/Signal.m b/Classes/00_signals/Signal.m index 97a8f14..613b7a0 100644 --- a/Classes/00_signals/Signal.m +++ b/Classes/00_signals/Signal.m @@ -815,11 +815,18 @@ classdef Signal pkpos = sort(pkpos); - % if mean(w) > 15 || mean(p) > 15 - % return - % else - % sequenceFound = 1; - % end + if isempty(pks) + warning(['Error in findpeaks, ususally the seuqnece is too short. No Peaks detected']); + return + end + + if max(p) < 0.3 || median(w) > 15 + %median(w) > 15 part means “reject if the detected correlation peaks are too broad.” That can be sensible: a true sync peak should often be sharp. + warning(['Error in findpeaks, ususally the seuqnece is too short. max(p) = ',num2str(max(p)),'; median(w)=',num2str(median(w)),'']); + return + end + + sequenceFound = 1; if options.debug_plots figure(121212);clf diff --git a/Classes/01_transmit/PAMmapper.m b/Classes/01_transmit/PAMmapper.m index 83287be..3124edd 100644 --- a/Classes/01_transmit/PAMmapper.m +++ b/Classes/01_transmit/PAMmapper.m @@ -462,6 +462,44 @@ classdef PAMmapper end + function [out, levels] = splitByReferenceLevels(obj, data_in, reference_in, options) + % Split received samples by the transmitted/reference PAM level. + % Unlike separate_pamlevels, this does not decide the RX level. + arguments + obj + data_in + reference_in + options.levels = [] + options.tolerance (1,1) double {mustBeNonnegative} = 0 + end + + data = obj.toNumericVector(data_in); + reference = obj.toNumericVector(reference_in); + + if numel(data) ~= numel(reference) + error("PAMmapper:LengthMismatch", ... + "data_in and reference_in must have the same number of samples."); + end + + if isempty(options.levels) + levels = unique(reference); + else + levels = options.levels(:).'; + end + + out = NaN(numel(levels), numel(reference)); + + for levelIdx = 1:numel(levels) + if options.tolerance == 0 + levelMask = reference == levels(levelIdx); + else + levelMask = abs(reference - levels(levelIdx)) <= options.tolerance; + end + + out(levelIdx, levelMask) = data(levelMask); + end + end + function [Signal_out] = quantize(obj,Signal_in,options) arguments obj @@ -621,5 +659,19 @@ classdef PAMmapper end + methods (Access = private) + + function values = toNumericVector(~, signalLike) + if isa(signalLike, "Signal") + values = signalLike.signal; + else + values = signalLike; + end + + values = values(:).'; + end + + end + end diff --git a/Classes/04_DSP/Equalizer/FFE_DCremoval_adaptive_mu.m b/Classes/04_DSP/Equalizer/FFE_DCremoval_adaptive_mu.m index 911370e..f7cb053 100644 --- a/Classes/04_DSP/Equalizer/FFE_DCremoval_adaptive_mu.m +++ b/Classes/04_DSP/Equalizer/FFE_DCremoval_adaptive_mu.m @@ -1,12 +1,10 @@ classdef FFE_DCremoval_adaptive_mu < handle - % Implementation of plain and simple FFE. - % 1) Training mode (stable performance when you use NLMS) - % 2) Decision directed mode - - % Eq = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",0); + % FFE variant for MPI/DC-removal experiments. + % With dc_buffer_len <= 1, ffe_buffer_len <= 1 and no smoothing, this + % follows FFE.m semantics so MPI-reduction changes can be isolated. properties - sps % usually 2 + sps order e e_tr @@ -16,28 +14,37 @@ classdef FFE_DCremoval_adaptive_mu < handle mu_tr epochs_tr + adaption_technique + dd_mode mu_dd epochs_dd - + dd_len_fraction mu_dc + e_dc + + P + dc_buffer_len - adaptive_mu_mode - ffe_buffer_len - smoothing_buffer_length smoothing_buffer_update constellation decide + + save_debug = 0; + debug_struct + + optmize_mus = 0; + mu_optimization + mu_optimization_iter = 0; end methods function obj = FFE_DCremoval_adaptive_mu(options) arguments(Input) - options.sps = 2; options.order = 15; @@ -45,23 +52,31 @@ classdef FFE_DCremoval_adaptive_mu < handle options.mu_tr = 0; options.epochs_tr = 5; + options.adaption_technique adaption_method = adaption_method.lms; + options.dd_mode = 1; options.mu_dd = 1e-5; options.epochs_dd = 5; + options.dd_len_fraction = 0.25; options.mu_dc = 0.05; options.dc_buffer_len = 1; - - options.ffe_buffer_len = 1; - options.adaptive_mu_mode = 1; - + options.ffe_buffer_len = 1; options.smoothing_buffer_length = 0; options.smoothing_buffer_update = 0; + options.decide = false; + options.save_debug = 0; + options.optmize_mus = 0; end - assert(options.dc_buffer_len>0); + assert(options.dc_buffer_len >= 0); + assert(options.ffe_buffer_len >= 0); + assert(options.smoothing_buffer_length >= 0); + if options.smoothing_buffer_length > 0 + assert(options.smoothing_buffer_update > 0); + end fn = fieldnames(options); for n = 1:numel(fn) @@ -69,283 +84,328 @@ classdef FFE_DCremoval_adaptive_mu < handle end obj.e = zeros(obj.order,1); + obj.e_dc = 0; obj.error = 0; obj.dc_buffer_len = floor(obj.dc_buffer_len); - + obj.ffe_buffer_len = floor(obj.ffe_buffer_len); + obj.smoothing_buffer_length = floor(obj.smoothing_buffer_length); + obj.smoothing_buffer_update = floor(obj.smoothing_buffer_update); end function [X,Noi] = process(obj, X, D) - - % actual processing of the signal (steps 1. - 3.) - % 1 normalize RMS X = X.normalize("mode","rms"); obj.constellation = unique(D.signal); + obj.e_dc = 0; - % if obj.smoothing_buffer_length > 0 - % % Apply A1 filter smoothing - % % Calculate the moving sum with the window size N1 - % moving_sum = movsum(X.signal, [obj.smoothing_buffer_length,0]); - % - % % Initialize the output smoothed signal - % X.signal = X.signal - (1 / obj.smoothing_buffer_length) * moving_sum; - % end + 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.e_dc = 0; + obj.P = (1/delta) * eye(obj.order); + end - % Training Mode training = 1; - obj.equalize(X.signal, D.signal,obj.mu_tr,obj.epochs_tr,obj.len_tr,training); + showviz = 0; + obj.equalize(X.signal,D.signal,obj.mu_tr,obj.epochs_tr,obj.len_tr,training,showviz); obj.e_tr = obj.e; - % Decision Directed Mode - N = X.length; + n = X.length; training = 0; - [signal,decision]=obj.equalize(X.signal, D.signal,obj.mu_dd,obj.epochs_dd,N,training); + if obj.dd_mode + n_dd = obj.ddLength(n); + obj.equalize(X.signal,D.signal,obj.mu_dd,obj.epochs_dd,n_dd,training,showviz); + end + [signal,decision] = obj.applyCurrentTaps(X.signal,n); - % Output Signal if obj.decide X.signal = decision; else X.signal = signal; end - X.fs = D.fs; %change sampling frequency of outgoing signal from fdac e.g. 2 sps to symbol spaced = fsym + + X.fs = D.fs; lbdesc = [num2str(obj.order),' tap FFE']; - X = X.logbookentry(lbdesc); % append to logbook + X = X.logbookentry(lbdesc); + Noi = X; Noi = X - D; - end - function [y,d_hat] = equalize(obj, x, d, mu_lms, epochs, N, training) - % Equalize with adaptive DC-removal, VSS, and parallel-buffered DC updates - % Added: FFE gradient buffering in DD mode (error buffer) with update every obj.dc_buffer_len symbols - + function [y,d_hat] = equalize(obj,x,d,mu,epochs,N,training,showviz) arguments obj x d - mu_lms % LMS step-size (or 0 for NLMS) - epochs % number of training/DD epochs - N % number of samples to process - training % boolean flag: true->training mode, false->DD mode + mu + epochs + N + training + showviz = 0 %#ok end - if isempty(obj.e) - obj.e = zeros(obj.order,1); - end - - % Zero-padding for filter memory x = [zeros(floor(obj.order/2),1); x; zeros(obj.order,1)]; + lambda = mu; - % Initialize storage - numSymbols = ceil(N/obj.sps); - y = zeros(numSymbols,1); - d_hat = zeros(numSymbols,1); - err = NaN(numSymbols,numel(obj.constellation)); - e_dc_save= zeros(numSymbols,1); - - % DC-adaptation parameters - P_err = 0; % running error power - alpha = 0.98; % forgetting factor for error power - err_prev = 0; % previous error sample for VSS correlation - gamma_dc = 1e-6; % meta step-size for DC VSS - mu_min = 1e-6; % lower bound for mu_dc - mu_max = 3e-1; % upper bound for mu_dc - - % DC removal buffer - L = obj.dc_buffer_len; % buffer length - e_dc_buf = NaN(L,1); - e_dc_est = 0; - - % FFE gradient buffer (DD mode only) - L_grad = obj.ffe_buffer_len; % buffer length - if ~training - % each column holds one past gradient of length obj.order - grad_buf = NaN(obj.order, L_grad); + if training + mask = ones(obj.order,1); + else + mask = zeros(obj.order,1); + mask(900:end) = 1; + mask(ceil(length(obj.e)/2)) = 1; end - smth_buffer = zeros(1, obj.smoothing_buffer_length); - smth_mean = 0; - % Main loop - for epoch = 1:epochs - s = 0; - for sample = 1:obj.sps:N - s = s + 1; + mask = ones(obj.order,1); + always_ideal_decision = 0; + grad = 0; + weight = 0; + update = 0; + + if mu == 0 || (~obj.dd_mode && ~training) + epochs = 1; + end + + dc_buffer_enabled = obj.mu_dc ~= 0 && obj.dc_buffer_len > 1; + adaptive_dc_enabled = dc_buffer_enabled && obj.adaptive_mu_mode; + if dc_buffer_enabled + e_dc_buffer = NaN(obj.dc_buffer_len,1); + end + + ffe_buffer_enabled = ~training && obj.ffe_buffer_len > 1 && ... + obj.adaption_technique ~= adaption_method.rls; + if ffe_buffer_enabled + grad_buffer = NaN(obj.order,obj.ffe_buffer_len); + end + + if obj.smoothing_buffer_length > 0 + smoothing_buffer = zeros(1,obj.smoothing_buffer_length); + smoothing_mean = 0; + end + + P_err = 0; + alpha = 0.98; + err_prev = 0; + gamma_dc = 1e-6; + mu_min = 1e-6; + mu_max = 3e-1; + + for epoch = 1 : epochs + symbol = 0; + for sample = 1 : obj.sps : N + symbol = symbol + 1; if obj.smoothing_buffer_length > 0 - smth_buffer = circshift(smth_buffer,1,2); - smth_buffer(1) = x(sample); - if mod(s, obj.smoothing_buffer_update) == 0 - smth_mean = mean(smth_buffer); + smoothing_buffer = circshift(smoothing_buffer,1,2); + smoothing_buffer(1) = x(sample); + if mod(symbol,obj.smoothing_buffer_update) == 0 + smoothing_mean = mean(smoothing_buffer); end - x(sample:sample+obj.sps-1) = x(sample:sample+obj.sps-1)-smth_mean; + x(sample:sample+obj.sps-1) = x(sample:sample+obj.sps-1) - smoothing_mean; end U = x(obj.order+sample-1:-1:sample); - %-- 1) filter output with DC correction - y(s) = e_dc_est + obj.e.'*U; + y(symbol,1) = obj.e_dc + (obj.e.*mask).' * U; - %-- 2) decision if training - [~, idx] = min(abs(d(s) - obj.constellation)); + d_hat(symbol,1) = d(symbol); else - [~, idx] = min(abs(y(s) - obj.constellation)); - end - d_hat(s) = obj.constellation(idx); - - %-- 3) error - e_val = y(s) - d_hat(s); - if epoch == epochs - - err(s,idx) = e_val; - true_err(s,idx) = y(s) - d(s); - end - - %-- 4) tap-weight update: training immediate, DD buffered - if training - % immediate update (LMS or NLMS) - if mu_lms ~= 0 - obj.e = obj.e - mu_lms * e_val * U; + if ~always_ideal_decision + [~,symbol_idx] = min(abs(y(symbol) - obj.constellation)); + d_hat(symbol,1) = obj.constellation(symbol_idx); else - normU = (U.'*U) + eps; - obj.e = obj.e - e_val * U / normU; + d_hat(symbol,1) = d(symbol); end - else - if 0 - % buffer gradient - if mu_lms ~= 0 - grad = e_val * U; - else + end + + err(symbol) = d_hat(symbol) - y(symbol); %#ok + true_err(symbol) = y(symbol) - d(symbol); %#ok + + if training || obj.dd_mode + switch obj.adaption_technique + case adaption_method.lms + weight = mu; + grad = err(symbol) * U; + update = grad * weight; + + case adaption_method.nlms normU = (U.'*U) + eps; - grad = e_val * U / normU; - end - % shift and insert - grad_buf = circshift(grad_buf, 1, 2); - grad_buf(:,1) = grad; - % update once every L symbols - if mod(s, L_grad) == 0 - avg_grad = mean(grad_buf, 2, 'omitnan'); - if mu_lms ~= 0 - obj.e = obj.e - mu_lms * avg_grad; - else - obj.e = obj.e - avg_grad; - end - end + weight = mu / normU; + grad = err(symbol) * U; + update = grad * weight; + + case adaption_method.rls + denom = lambda + U.' * obj.P * U; + k = (obj.P * U) / denom; + update = k * err(symbol); end - end - - - %-- 5) DC adaptation - if obj.mu_dc ~= 0 - - if obj.adaptive_mu_mode - - % VSS for mu_dc - delta_mu = gamma_dc * e_val * err_prev * (U.'*U); - obj.mu_dc = min(max(obj.mu_dc + delta_mu, mu_min), mu_max); - err_prev = e_val; - - % DC buffer update & periodic estimate -- new "P_err" is "e_val^2" - P_err = alpha*P_err + (1-alpha)*e_val^2; - mu_dc_norm = obj.mu_dc / (P_err + eps); + if ffe_buffer_enabled + grad_buffer = circshift(grad_buffer,1,2); + grad_buffer(:,1) = update; + if mod(symbol,obj.ffe_buffer_len) == 0 + obj.e = obj.e + mean(grad_buffer,2,"omitnan"); + end else - - % DC buffer update & periodic estimate - % P_err = alpha*P_err + (1-alpha)*e_val^2; - % mu_dc_norm = obj.mu_dc / (P_err + eps); - - mu_dc_norm = obj.mu_dc; - + obj.e = obj.e + update; end - e_dc_buf = circshift(e_dc_buf, 1); - e_dc_buf(1) = e_dc_est - mu_dc_norm * e_val; - - if mod(s, L) == 0 - e_dc_est = median(e_dc_buf, 'omitnan'); + if obj.adaption_technique == adaption_method.rls + obj.P = (1/lambda) * (obj.P - k * (U.' * obj.P)); end - P_err_save(s) = P_err; - % Pcorr_save(s) = e_val * err_prev; - Ucorr_save(s) = (U.'*U); - mu_dc_save(s) = mu_dc_norm; - e_dc_save(s) = e_dc_est; + if obj.mu_dc ~= 0 + if adaptive_dc_enabled + delta_mu = gamma_dc * err(symbol) * err_prev * (U.'*U); + obj.mu_dc = min(max(obj.mu_dc + delta_mu,mu_min),mu_max); + err_prev = err(symbol); + P_err = alpha*P_err + (1-alpha)*err(symbol)^2; + mu_dc_eff = obj.mu_dc / (P_err + eps); + else + mu_dc_eff = obj.mu_dc; + end + if dc_buffer_enabled + e_dc_buffer = circshift(e_dc_buffer,1); + e_dc_buffer(1) = obj.e_dc + mu_dc_eff * err(symbol); + if mod(symbol,obj.dc_buffer_len) == 0 + obj.e_dc = median(e_dc_buffer,"omitnan"); + end + else + obj.e_dc = obj.e_dc + mu_dc_eff * err(symbol); + end + end end - % store instantaneous squared error - obj.error(epoch, s) = e_val^2; + if obj.save_debug + obj.debug_struct.error(epoch,symbol) = err(symbol) * err(symbol)'; + + if training + obj.debug_struct.error_tr(epoch,symbol) = err(symbol) * err(symbol)'; + obj.debug_struct.update_tr(epoch,symbol) = update.'*update ./ rms(obj.e); + end + end + + % obj.error(epoch,symbol) = err(symbol) * err(symbol)'; end end - - % Optional plotting in DD mode (uncomment if needed) - if 0%~training - - constellation = unique(d); - lvlcol = cbrewer2('Paired', numel(constellation)*2); - lvlcol = lvlcol(2:2:end, :); - - true_err(true_err==0) = NaN; - true_errmoverr = movsum(true_err, 4096, 'omitnan'); - true_errmoverr = true_errmoverr./rms(true_errmoverr); - - moverr = movsum(err, [100,100], 'omitnan'); - moverr = moverr./rms(moverr); - - figure(500); clf - hold on - % 1st subplot: true_errmoverr - % subplot(2,2,1); hold on - % for k = 1:4 - % scatter(1:numSymbols, true_errmoverr(:,k), 1, lvlcol(k,:), '.'); - % end - - % scatter(1:numSymbols, Ucorr_save./rms(Ucorr_save), 1, lvlcol(1,:), '.','DisplayName','Ucorr_save'); - % scatter(1:numSymbols, Pcorr_save./rms(Pcorr_save), 1, lvlcol(1,:), '.','DisplayName','P_corr'); - % scatter(1:numSymbols, P_err_save, 1, lvlcol(1,:), '.','DisplayName','P_err'); - % scatter(1:numSymbols, mu_dc_save, 1, lvlcol(2,:), '.','DisplayName','adapted value of $\mu_{DC}$'); - % scatter(1:numSymbols, sum(moverr,2,'omitnan'), 1, lvlcol(1,:), '.','DisplayName','Mov Error $\hat{d}$ - x over all levels'); - scatter(1:numSymbols, sum(e_dc_save,2,'omitnan'), 1, lvlcol(2,:), '.','DisplayName','Est. Error that is subtracted'); - title('Moving Sum Error'); - hold off - legend - - % 2nd subplot: moverr - subplot(2,2,2); hold on - for k = 1:4 - scatter(1:numSymbols, moverr(:,k), 1, lvlcol(k,:), '.'); - end - title('Moving Sum Error'); - hold off - legend - - % 3rd subplot: err - subplot(2,2,3); hold on - for k = 1:4 - scatter(1:numSymbols, err(:,k), 1, lvlcol(k,:), '.'); - end - title('Error'); - hold off - legend - - % 4th subplot: err + obj.constellation' - subplot(2,2,4); hold on - for k = 1:4 - scatter(1:numSymbols, err(:,k) + obj.constellation(k), 1, lvlcol(k,:), '.'); - end - yline(obj.constellation, '--k'); - title('Error + Constellation'); - hold off - legend - - sgtitle('Error Analysis Subplots'); - - end end + function [y,d_hat] = applyCurrentTaps(obj,x,N) + x = [zeros(floor(obj.order/2),1); x; zeros(obj.order,1)]; + for sample = 1 : obj.sps : N + symbol = (sample - 1) / obj.sps + 1; + U = x(obj.order+sample-1:-1:sample); + y(symbol,1) = obj.e_dc + obj.e.' * U; + [~,symbol_idx] = min(abs(y(symbol) - obj.constellation)); + d_hat(symbol,1) = obj.constellation(symbol_idx); + end + end + function N_dd = ddLength(obj,N) + if isempty(obj.dd_len_fraction) || obj.dd_len_fraction <= 0 || obj.dd_len_fraction >= 1 + N_dd = N; + return + end + + N_dd = floor(N * obj.dd_len_fraction); + N_dd = max(obj.sps,N_dd); + N_dd = min(N,N_dd); + 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_dc_range = [1e-5, 1e-1]; + + 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 + optimize_mu_dc = obj.mu_dc ~= 0; + if optimize_mu_dc + vars = [vars, optimizableVariable("mu_dc",mu_dc_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 + if optimize_mu_dc + obj.mu_dc = obj.mu_optimization.XAtMinObjective.mu_dc; + end + objective_db = 10*log10(obj.mu_optimization.MinObjective); + if obj.dd_mode && optimize_mu_dc + fprintf("\nFFE_DCremoval_adaptive_mu opt done: mu_tr=%9.3e, mu_dd=%9.3e, mu_dc=%9.3e, MSE=%9.3e, MSE_dB=%7.2f dB\n", ... + obj.mu_tr,obj.mu_dd,obj.mu_dc,obj.mu_optimization.MinObjective,objective_db); + elseif obj.dd_mode + fprintf("\nFFE_DCremoval_adaptive_mu opt done: mu_tr=%9.3e, mu_dd=%9.3e, MSE=%9.3e, MSE_dB=%7.2f dB\n", ... + obj.mu_tr,obj.mu_dd,obj.mu_optimization.MinObjective,objective_db); + elseif optimize_mu_dc + fprintf("\nFFE_DCremoval_adaptive_mu opt done: mu_tr=%9.3e, mu_dc=%9.3e, MSE=%9.3e, MSE_dB=%7.2f dB\n", ... + obj.mu_tr,obj.mu_dc,obj.mu_optimization.MinObjective,objective_db); + else + fprintf("\nFFE_DCremoval_adaptive_mu opt done: mu_tr=%9.3e, MSE=%9.3e, MSE_dB=%7.2f dB\n", ... + obj.mu_tr,obj.mu_optimization.MinObjective,objective_db); + end + end + + function objective = muObjective(obj,params,x,d) + old_debug = obj.save_debug; + old_mu_dc = obj.mu_dc; + obj.save_debug = 1; + if isprop(params,"mu_dc") + obj.mu_dc = params.mu_dc; + end + obj.e = zeros(obj.order,1); + obj.e_dc = 0; + 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,obj.ddLength(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 + objective_db = 10*log10(objective); + obj.mu_optimization_iter = obj.mu_optimization_iter + 1; + optimize_mu_dc = isprop(params,"mu_dc"); + if obj.dd_mode && optimize_mu_dc + fprintf("\rFFE_DCremoval_adaptive_mu opt %02d: mu_tr=%9.3e, mu_dd=%9.3e, mu_dc=%9.3e, MSE=%9.3e, MSE_dB=%7.2f dB", ... + obj.mu_optimization_iter,params.mu_tr,params.mu_dd,params.mu_dc,objective,objective_db); + elseif obj.dd_mode + fprintf("\rFFE_DCremoval_adaptive_mu opt %02d: mu_tr=%9.3e, mu_dd=%9.3e, MSE=%9.3e, MSE_dB=%7.2f dB", ... + obj.mu_optimization_iter,params.mu_tr,params.mu_dd,objective,objective_db); + elseif optimize_mu_dc + fprintf("\rFFE_DCremoval_adaptive_mu opt %02d: mu_tr=%9.3e, mu_dc=%9.3e, MSE=%9.3e, MSE_dB=%7.2f dB", ... + obj.mu_optimization_iter,params.mu_tr,params.mu_dc,objective,objective_db); + else + fprintf("\rFFE_DCremoval_adaptive_mu opt %02d: mu_tr=%9.3e, MSE=%9.3e, MSE_dB=%7.2f dB", ... + obj.mu_optimization_iter,params.mu_tr,objective,objective_db); + end + obj.save_debug = old_debug; + obj.mu_dc = old_mu_dc; + end end end - diff --git a/Classes/Warehouse_class/classes/DataStorage.m b/Classes/Warehouse_class/classes/DataStorage.m index 394a12f..49419d3 100644 --- a/Classes/Warehouse_class/classes/DataStorage.m +++ b/Classes/Warehouse_class/classes/DataStorage.m @@ -95,14 +95,14 @@ classdef DataStorage < handle end - function addStorage(obj,varName) - % add a storage - - storage = cell(obj.dim); - - obj.sto.(string(varName)) = storage; - - end + function addStorage(obj,varName) + % add a storage + + storage = cell(obj.getStorageSize()); + + obj.sto.(string(varName)) = storage; + + end function addValueToStorage(obj, valueToStore ,storageVarName, varargin) @@ -273,14 +273,18 @@ classdef DataStorage < handle end %append to index list :-) - fn_=fieldnames(obj.sto); - n_ = fn_{1}; - lin_idx(c,:) = sub2ind(size(obj.sto.(n_)),indices{:}); - % lin_idx(c,:) = eval(['sub2ind(size(obj.sto.',n_,')',str,');']); - - end - - + if isscalar(indices) + lin_idx(c,:) = indices{1}; + else + fn_=fieldnames(obj.sto); + n_ = fn_{1}; + lin_idx(c,:) = sub2ind(size(obj.sto.(n_)),indices{:}); + % lin_idx(c,:) = eval(['sub2ind(size(obj.sto.',n_,')',str,');']); + end + + end + + end % Mapping for single Index @@ -291,24 +295,33 @@ classdef DataStorage < handle - function [phys_indices,param_name] = getPhysIndicesByLinIndex(obj, lin_idx) - % Converts a linear index into the corresponding physical parameter values - % Inputs: - % - lin_idx: The linear index within the storage array + function [phys_indices,param_name] = getPhysIndicesByLinIndex(obj, lin_idx) + % Converts a linear index into the corresponding physical parameter values + % Inputs: + % - lin_idx: The linear index within the storage array % Output: % - phys_indices: A cell array containing the physical parameter values for each dimension - % Initialize output cell array - phys_indices = cell(1, numel(obj.fn)); - - % Convert linear index to subscript indices - [subscripts{1:numel(obj.dim)}] = ind2sub(obj.dim, lin_idx); - - % Map subscripts to physical values for each parameter - for i = 1:numel(obj.fn) - param_name{i} = obj.fn(i); - phys_indices{i} = obj.parameter.(param_name{i}).getPhysForIndex(subscripts{i}); - end + % Initialize output cell array + phys_indices = cell(1, numel(obj.fn)); + param_name = cell(1, numel(obj.fn)); + + if isempty(obj.fn) + return + end + + % Convert linear index to subscript indices + if isscalar(obj.dim) + subscripts = {lin_idx}; + else + [subscripts{1:numel(obj.dim)}] = ind2sub(obj.dim, lin_idx); + end + + % Map subscripts to physical values for each parameter + for i = 1:numel(obj.fn) + param_name{i} = obj.fn(i); + phys_indices{i} = obj.parameter.(param_name{i}).getPhysForIndex(subscripts{i}); + end end function [physStruct, stored_value] = getPhysAndValueByLinIndex(obj, storageVarName, lin_idx) @@ -327,8 +340,14 @@ classdef DataStorage < handle % Initialize an empty structure physStruct = struct(); - % Convert linear index to subscript indices - [subscripts{1:numel(obj.dim)}] = ind2sub(obj.dim, lin_idx); + % Convert linear index to subscript indices + if isempty(obj.fn) + subscripts = {}; + elseif isscalar(obj.dim) + subscripts = {lin_idx}; + else + [subscripts{1:numel(obj.dim)}] = ind2sub(obj.dim, lin_idx); + end % Map subscripts to physical values and parameter names for each dimension for i = 1:numel(obj.fn) @@ -343,17 +362,31 @@ classdef DataStorage < handle stored_value = obj.sto.(storageVarName){lin_idx}; end - function num_elements = getLastLinIndice(obj) - % Returns all possible linear indices for the data structure - % Output: - % - lin_indices: A column vector containing all linear indices for the storage array - - % Calculate the total number of elements in the storage array - num_elements = prod(obj.dim); - end - - - - end + function num_elements = getLastLinIndice(obj) + % Returns all possible linear indices for the data structure + % Output: + % - lin_indices: A column vector containing all linear indices for the storage array + + % Calculate the total number of elements in the storage array + if isempty(obj.dim) + num_elements = 1; + else + num_elements = prod(obj.dim); + end + end + + function storageSize = getStorageSize(obj) + if isempty(obj.dim) + storageSize = [1, 1]; + elseif isscalar(obj.dim) + storageSize = [obj.dim, 1]; + else + storageSize = obj.dim; + end + end + + + + end end diff --git a/Functions/EQ_recipes/dsp_recipe_minimal.m b/Functions/EQ_recipes/dsp_recipe_minimal.m new file mode 100644 index 0000000..fca662b --- /dev/null +++ b/Functions/EQ_recipes/dsp_recipe_minimal.m @@ -0,0 +1,36 @@ +function output = dsp_recipe_minimal(Scpe_sig_raw, Symbols, Tx_bits, options) +%DSP_RECIPE_MINIMAL Minimal example recipe for the DSP job framework. +% This recipe intentionally performs only light preprocessing and records +% summary values. It demonstrates the recipe interface without running a +% full equalizer chain. + + arguments + Scpe_sig_raw + Symbols + Tx_bits + options.fsym + options.M + options.duob_mode + options.dataTable table + options.userParameters struct = struct() + options.debug_plots (1,1) logical = false + end + + Scpe_sig = preprocessSignal(Scpe_sig_raw, Symbols, options.fsym, ... + "mode", "auto", ... + "debug_plots", options.debug_plots); + + eq_ffe = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-1,"mu_tr",0.4,"order",25,... + "sps",2,"decide",0,"optmize_mus",0,"dd_mode",options.userParameters.dd_mode,"adaption_technique","nlms","mu_dc",1.021e-05); + + ffe_results = ffe(eq_ffe, options.M, Scpe_sig, Symbols, Tx_bits, ... + "precode_mode", options.duob_mode, ... + 'showAnalysis', options.debug_plots, ... + "postFFE", [], ... + "eth_style_symbol_mapping", 0); + + ffe_results.config.equalizer_structure = "ffe"; + ffe_results.metrics.print("description",'FFE'); + output.ffe_package = ffe_results; + +end diff --git a/Functions/EQ_blocks/dsp_scope_signal.m b/Functions/EQ_recipes/dsp_scope_signal.m similarity index 100% rename from Functions/EQ_blocks/dsp_scope_signal.m rename to Functions/EQ_recipes/dsp_scope_signal.m diff --git a/Functions/EQ_recipes/mpi_recipe_dev.m b/Functions/EQ_recipes/mpi_recipe_dev.m new file mode 100644 index 0000000..5e5abfa --- /dev/null +++ b/Functions/EQ_recipes/mpi_recipe_dev.m @@ -0,0 +1,78 @@ +function output = mpi_recipe_dev(Scpe_sig_raw, Symbols, Tx_bits, options) +%mpi_recipe_dev Minimal example recipe for the DSP job framework. +% This recipe intentionally performs only light preprocessing and records +% summary values. It demonstrates the recipe interface without running a +% full equalizer chain. + + arguments + Scpe_sig_raw + Symbols + Tx_bits + options.fsym + options.M + options.duob_mode + options.dataTable table + options.userParameters struct = struct() + options.debug_plots (1,1) logical = false + end + + Scpe_sig = preprocessSignal(Scpe_sig_raw, Symbols, options.fsym, ... + "mode", "auto", ... + "debug_plots", options.debug_plots); + + + mu_dc = 0; % 1e-5 + dc_buffer_len = 0; + ffe_buffer_len = 0; + smoothing_buffer_length = options.userParameters.smoothing_length; + smoothing_buffer_update = 1; + + eq_settings = { ... + "epochs_tr", 5, ... + "epochs_dd", 5, ... + "len_tr", 4096*2, ... + "mu_dd", 1e-1, ... + "mu_tr", 0.4, ... + "order", 25, ... + "sps", 2, ... + "decide", 0, ... + "optmize_mus", 0, ... + "dd_mode", 1, ... + "adaption_technique", "nlms", ... + "mu_dc", mu_dc}; + + eq_ffe = FFE(eq_settings{:}); + + % showLevelScatter(Scpe_sig_raw, Symbols, ... + % "fsym", options.fsym, ... + % "fignum", options.dataTable.run_id, ... + % "normalize", true); + + %% + % tic + % ffe_results = ffe(eq_ffe, options.M, Scpe_sig, Symbols, Tx_bits, ... + % "precode_mode", options.duob_mode, ... + % 'showAnalysis', options.debug_plots, ... + % "postFFE", [], ... + % "eth_style_symbol_mapping", 0); + % toc + % ffe_results.metrics.print("description",'FFE'); + % output.ffe_package = ffe_results; + + %% + eq_ffe_dcr = FFE_DCremoval_adaptive_mu(eq_settings{:}, ... + "dc_buffer_len",dc_buffer_len, ... + "ffe_buffer_len",ffe_buffer_len,... + "smoothing_buffer_length",smoothing_buffer_length,... + "smoothing_buffer_update",smoothing_buffer_update); + tic + ffe_results_dcr = ffe(eq_ffe_dcr, options.M, Scpe_sig, Symbols, Tx_bits, ... + "precode_mode", options.duob_mode, ... + 'showAnalysis', options.debug_plots, ... + "postFFE", [], ... + "eth_style_symbol_mapping", 0); + toc + ffe_results_dcr.metrics.print("description",'FFE DCR'); + output.ffe_dcr_package = ffe_results_dcr; + +end diff --git a/Functions/EQ_visuals/showLevelScatter.m b/Functions/EQ_visuals/showLevelScatter.m index f39dcb8..1565487 100644 --- a/Functions/EQ_visuals/showLevelScatter.m +++ b/Functions/EQ_visuals/showLevelScatter.m @@ -1,133 +1,370 @@ -function [symbols_for_lvl,avg_for_lvl] = showLevelScatter(eq_signal,ref_symbols,options) +function [symbols_for_lvl, avg_for_lvl, info] = showLevelScatter(rxInput, refSymbols, options) +%SHOWLEVELSCATTER Plot received samples separated by reference PAM level. +% Supports plain numeric vectors, Signal objects, synchronized scope cell +% arrays, and raw unsynchronized Signal input. Raw Signal input is +% synchronized to refSymbols and stitched before plotting. + arguments - eq_signal - ref_symbols - options.fignum (1,1) double = NaN % Default to NaN if not provided - options.displayname (1,:) char = '' % Default to an empty string if not provided - options.f_sym =1e6; + rxInput + refSymbols + options.fignum (1,1) double = NaN + options.displayname (1,:) char = '' + options.f_sym double = [] + options.fsym double = [] + options.syncFs (1,1) double = 0 + options.shiftFs (1,1) double = 0 + options.shifts double = [] + options.maxOccurences (1,1) double = Inf + options.normalize (1,1) logical = false + options.debug_plots (1,1) logical = false + options.showPlot (1,1) logical = true + options.clear (1,1) logical = true + options.windowLength (1,1) double {mustBePositive, mustBeInteger} = 500 + options.xLimits double = [] + options.yLimits (1,2) double = [-3 3] + options.showStdAnnotations (1,1) logical = true end -plot_shit = 1; +fsym = resolveSymbolRate(rxInput, refSymbols, options); +[symbols_for_lvl, avg_for_lvl, xAxisUs, info] = prepareLevelScatterData(rxInput, refSymbols, fsym, options); -if isa(eq_signal,'Signal') - options.f_sym = eq_signal.fs; - eq_signal = eq_signal.signal; - assert(~isempty(options.f_sym),'No fsym given'); +if options.showPlot + plotLevelScatter(symbols_for_lvl, avg_for_lvl, refSymbols, xAxisUs, options); end -if isa(ref_symbols,'Signal') - ref_symbols = ref_symbols.signal; end +function fsym = resolveSymbolRate(rxInput, refSymbols, options) +if ~isempty(options.fsym) + fsym = options.fsym; +elseif ~isempty(options.f_sym) + fsym = options.f_sym; +elseif isa(refSymbols, "Signal") && ~isempty(refSymbols.fs) + fsym = refSymbols.fs; +elseif isa(rxInput, "Signal") && ~isempty(rxInput.fs) + fsym = rxInput.fs; +elseif iscell(rxInput) && ~isempty(rxInput) && isa(rxInput{1}, "Signal") && ~isempty(rxInput{1}.fs) + fsym = rxInput{1}.fs; +else + fsym = 1e6; +end +end -if plot_shit - % Determine the figure number to use or create a new figure - if isnan(options.fignum) - fig = figure; % Create a new figure and get its handle +function [symbols_for_lvl, avg_for_lvl, xAxisUs, info] = prepareLevelScatterData(rxInput, refSymbols, fsym, options) +refSignal = numericSignal(refSymbols); +info = defaultInfo(fsym); + +if iscell(rxInput) + [symbols_for_lvl, avg_for_lvl, info] = prepareCellInput(rxInput, refSymbols, fsym, options, info); +elseif isa(rxInput, "Signal") + [symbols_for_lvl, avg_for_lvl, info] = prepareSignalInput(rxInput, refSymbols, fsym, options, info); +else + rxSymbols = numericSignal(rxInput); + [symbols_for_lvl, avg_for_lvl] = levelScatterForOneSequence(rxSymbols, refSignal, options.windowLength); +end + +xAxisUs = ((1:size(avg_for_lvl, 2)) / fsym) * 1e6; +end + +function info = defaultInfo(fsym) +info = struct(); +info.found_sync = true; +info.startSamples = 1; +info.shifts = []; +info.fsym = fsym; +info.shiftFs = fsym; +info.varianceByLevel = []; +end + +function [symbols_for_lvl, avg_for_lvl, info] = prepareSignalInput(rxSignal, refSymbols, fsym, options, info) +rxAtSymbolRate = rxSignal.resample("fs_in", rxSignal.fs, "fs_out", fsym); +refSignal = numericSignal(refSymbols); + +if numel(rxAtSymbolRate.signal) == numel(refSignal) + rxSymbols = rxAtSymbolRate.signal; + if options.normalize + rxSymbols = normalizeNumericRms(rxSymbols); + end + + [symbols_for_lvl, avg_for_lvl] = levelScatterForOneSequence(rxSymbols, refSignal, options.windowLength); + info.varianceByLevel = var(symbols_for_lvl, 0, 2, "omitnan"); + return +end + +[scopeCell, shifts, shiftFs, foundSync] = synchronizeRawSignal(rxSignal, refSymbols, fsym, options); +info.found_sync = foundSync; +info.shifts = shifts; +info.shiftFs = shiftFs; + +if isempty(scopeCell) + symbols_for_lvl = []; + avg_for_lvl = []; + warning("showLevelScatter:NoScopeCells", ... + "No synchronized scope signal occurrences available."); + return +end + +[symbols_for_lvl, avg_for_lvl, startSamples] = stitchScopeCells(scopeCell, refSymbols, fsym, shifts, shiftFs, options); +info.startSamples = startSamples; +info.varianceByLevel = var(symbols_for_lvl, 0, 2, "omitnan"); +end + +function [symbols_for_lvl, avg_for_lvl, info] = prepareCellInput(scopeCell, refSymbols, fsym, options, info) +scopeCell = scopeCell(:); +if isempty(scopeCell) + symbols_for_lvl = []; + avg_for_lvl = []; + info.found_sync = false; + warning("showLevelScatter:NoScopeCells", ... + "No synchronized scope signal occurrences available."); + return +end + +shiftFs = options.shiftFs; +if shiftFs <= 0 + shiftFs = fsym; +end + +[symbols_for_lvl, avg_for_lvl, startSamples] = stitchScopeCells(scopeCell, refSymbols, fsym, options.shifts, shiftFs, options); +info.found_sync = true; +info.shifts = options.shifts; +info.shiftFs = shiftFs; +info.startSamples = startSamples; +info.varianceByLevel = var(symbols_for_lvl, 0, 2, "omitnan"); +end + +function [scopeCell, shifts, shiftFs, foundSync] = synchronizeRawSignal(rxSignal, refSymbols, fsym, options) +syncFs = options.syncFs; +if syncFs <= 0 + syncFs = 2*fsym; +end + +syncSignal = rxSignal.resample("fs_in", rxSignal.fs, "fs_out", syncFs); +if options.normalize + syncSignal = syncSignal.normalize("mode", "rms"); +end + +[~, scopeCell, ~, foundSync, shifts] = syncSignal.tsynch( ... + "reference", refSymbols, ... + "fs_ref", fsym, ... + "debug_plots", options.debug_plots); + +if options.shiftFs > 0 + shiftFs = options.shiftFs; +else + shiftFs = syncFs; +end +end + +function [symbols_for_lvl, avg_for_lvl, startSamples] = stitchScopeCells(scopeCell, refSymbols, fsym, shifts, shiftFs, options) +recordOccurrences = min(numel(scopeCell), options.maxOccurences); +scopeCell = scopeCell(1:recordOccurrences); +refSignal = numericSignal(refSymbols); +startSamples = getStartSamples(shifts, recordOccurrences, shiftFs, fsym, numel(refSignal)); + +levelScatter = cell(1, recordOccurrences); +levelAverage = cell(1, recordOccurrences); + +for occurrenceIdx = 1:recordOccurrences + occurrence = scopeCell{occurrenceIdx}; + if isa(occurrence, "Signal") + occurrence = occurrence.resample("fs_out", fsym); + occurrence = occurrence.signal; + end + + [levelScatter{occurrenceIdx}, levelAverage{occurrenceIdx}] = ... + levelScatterForOneSequence(occurrence, refSignal, options.windowLength); +end + +numLevels = size(levelScatter{1}, 1); +traceLength = max(startSamples(:).' + cellfun(@(x) size(x, 2), levelScatter) - 1); +symbols_for_lvl = NaN(numLevels, traceLength); +avg_for_lvl = NaN(numLevels, traceLength); + +for occurrenceIdx = 1:recordOccurrences + writeIdx = startSamples(occurrenceIdx):(startSamples(occurrenceIdx) + size(levelScatter{occurrenceIdx}, 2) - 1); + symbols_for_lvl(:, writeIdx) = levelScatter{occurrenceIdx}; + avg_for_lvl(:, writeIdx) = levelAverage{occurrenceIdx}; +end +end + +function startSamples = getStartSamples(shifts, recordOccurrences, shiftFs, fsym, symbolLength) +if isempty(shifts) + startSamples = ((0:recordOccurrences-1) .* symbolLength) + 1; + return +end + +shifts = shifts(:); +if numel(shifts) ~= recordOccurrences + positiveShifts = shifts(shifts >= 0); + if numel(positiveShifts) >= recordOccurrences + shifts = positiveShifts(1:recordOccurrences); else - fig = figure(options.fignum); % Use the specified figure number + warning("showLevelScatter:ShiftCountMismatch", ... + "Shift count does not match scope cell count. Using sequential stitching."); + startSamples = ((0:recordOccurrences-1) .* symbolLength) + 1; + return + end +end + +startSamples = round((shifts(1:recordOccurrences) - shifts(1)) ./ shiftFs .* fsym) + 1; +end + +function [symbols_for_lvl, avg_for_lvl] = levelScatterForOneSequence(rxSymbols, refSymbols, windowLength) +rxSymbols = numericSignal(rxSymbols); +refSymbols = numericSignal(refSymbols); + +assert(numel(rxSymbols) == numel(refSymbols), ... + 'showLevelScatter:LengthMismatch', ... + 'rxInput and refSymbols must have the same number of samples after resampling/synchronization.'); + +levels = unique(refSymbols); +[symbols_for_lvl, levels] = splitByReferenceLevels(rxSymbols, refSymbols, levels); +avg_for_lvl = NaN(numel(levels), numel(refSymbols)); + +for levelIdx = 1:numel(levels) + levelMask = ~isnan(symbols_for_lvl(levelIdx, :)); + levelSamples = symbols_for_lvl(levelIdx, levelMask); + if isempty(levelSamples) + continue + end + + smoothWindowLength = min(windowLength, numel(levelSamples)); + avg_for_lvl(levelIdx, levelMask) = movmean(levelSamples, smoothWindowLength, 'Endpoints', 'shrink'); + avg_for_lvl(levelIdx, :) = interpolateMissingLevelAverage(avg_for_lvl(levelIdx, :)); +end +end + +function [symbols_for_lvl, levels] = splitByReferenceLevels(rxSymbols, refSymbols, levels) +supportedPamLevels = [2 4 6 8 16]; + +if ismember(numel(levels), supportedPamLevels) + [symbols_for_lvl, levels] = PAMmapper(numel(levels), 0).splitByReferenceLevels( ... + rxSymbols, refSymbols, ... + "levels", levels); +else + symbols_for_lvl = NaN(numel(levels), numel(refSymbols)); + for levelIdx = 1:numel(levels) + levelMask = refSymbols == levels(levelIdx); + symbols_for_lvl(levelIdx, levelMask) = rxSymbols(levelMask); + end +end +end + +function plotLevelScatter(symbols_for_lvl, avg_for_lvl, refSymbols, xAxisUs, options) +if isempty(symbols_for_lvl) + return +end + +if isnan(options.fignum) + figure; +else + figure(options.fignum); + if options.clear clf; end end +hold on +numLevels = size(symbols_for_lvl, 1); +cols = cbrewer2("Paired", 2*numLevels); -rx_symbols = eq_signal; %./ rms(eq_signal); -correct_symbols = ref_symbols; -f_sym = options.f_sym; - -col = cbrewer2('Paired',numel(unique(correct_symbols))*2); -ccnt = -1; - -levels = unique(correct_symbols); -symbols_for_lvl = NaN(numel(levels),length(correct_symbols)); -start = 1; -ende = length(correct_symbols); - -for l = 1:numel(levels) - ccnt = ccnt+2; - - level_amplitude = levels(l); - - symbols_for_lvl(l,correct_symbols==level_amplitude) = rx_symbols(correct_symbols==level_amplitude); - std_lvl(l) = std(symbols_for_lvl(l,:),'omitnan'); - xax_in_sec = ((1:length(correct_symbols)) / f_sym) * 1e6; - - if plot_shit - scatter(xax_in_sec(start:ende),symbols_for_lvl(l,start:ende),10,'.','MarkerFaceAlpha',0.5,'MarkerEdgeAlpha',0.5,'MarkerEdgeColor',col(ccnt,:)); - hold on; - end - +for levelIdx = 1:numLevels + scatter(xAxisUs, symbols_for_lvl(levelIdx, :), 10, ".", ... + "MarkerEdgeColor", cols((2*levelIdx)-1, :), ... + "MarkerEdgeAlpha", 0.5); end -std_lvl = round(std_lvl,2); - -ccnt = 0; -avg_for_lvl = NaN(numel(levels),length(correct_symbols)); -% Add the windowed/ smoothed curves -for l = 1:numel(levels) - ccnt = ccnt+2; - level_amplitude = levels(l); - - L = 500; - movmean = 1/L .* movsum(rx_symbols(correct_symbols==level_amplitude),[L/2,L/2], 'Endpoints', 'fill'); - - avg_for_lvl(l,correct_symbols==level_amplitude) = movmean; - - nanx = isnan(avg_for_lvl(l,:)); - t = 1:numel(avg_for_lvl(l,:)); - avg_for_lvl(l,nanx) = interp1(t(~nanx), avg_for_lvl(l,~nanx), t(nanx)); - - xax_in_sec = ((1:length(correct_symbols)) / f_sym) * 1e6; - % xax_in_sec = 1:length(correct_symbols); - - if plot_shit - plot(xax_in_sec(start:ende),avg_for_lvl(l,start:ende),'Color',col(ccnt,:)); - end - hold on +for levelIdx = 1:numLevels + plot(xAxisUs, avg_for_lvl(levelIdx, :), ... + "LineWidth", 1, ... + "Color", cols(2*levelIdx, :)); end - - -if 0 - annotation(fig,'textbox',... - [0.660523809523809 0.844444444444448 0.133523809523809 0.0603174603174607],... - 'String',['\sigma = ',num2str(std_lvl(4))],... - 'LineWidth',1.8,... - 'LineStyle','none',... - 'FontSize',12,... - 'FitBoxToText','off'); - - % Create textbox - annotation(fig,'textbox',... - [0.667666666666665 0.642857142857147 0.133523809523809 0.0603174603174607],... - 'String',['\sigma = ',num2str(std_lvl(3))],... - 'LineWidth',1.8,... - 'LineStyle','none',... - 'FontSize',12,... - 'FitBoxToText','off'); - - % Create textbox - annotation(fig,'textbox',... - [0.671238095238093 0.442857142857148 0.133523809523809 0.0603174603174608],... - 'String',['\sigma = ',num2str(std_lvl(2))],... - 'LineWidth',1.8,... - 'LineStyle','none',... - 'FontSize',12,... - 'FitBoxToText','off'); - - % Create textbox - annotation(fig,'textbox',... - [0.670047619047616 0.265079365079371 0.133523809523809 0.0603174603174608],... - 'String',['\sigma = ',num2str(std_lvl(1))],... - 'LineWidth',1.8,... - 'LineStyle','none',... - 'FontSize',12,... - 'FitBoxToText','off'); +refValues = numericSignal(refSymbols); +yline(unique(refValues), "HandleVisibility", "off"); +if options.showStdAnnotations + annotateLevelStd(symbols_for_lvl, avg_for_lvl, xAxisUs, options.xLimits); end -if plot_shit -% yline(levels); xlabel('Time in $\mu$s'); ylabel('Normalized Amplitude'); -ylim([-3 3]); +xlim(getXLimits(options.xLimits, xAxisUs, symbols_for_lvl)); +ylim(options.yLimits); +grid on +drawnow; +end + +function xLimits = getXLimits(configuredLimits, xAxisUs, symbols_for_lvl) +if ~isempty(configuredLimits) + xLimits = configuredLimits; + return +end + +filledColumns = any(~isnan(symbols_for_lvl), 1); +if ~any(filledColumns) + xLimits = [xAxisUs(1), xAxisUs(end)]; + return +end + +lastFilledColumn = find(filledColumns, 1, "last"); +xMax = xAxisUs(lastFilledColumn); +xLimits = [0, 1.05*xMax]; +end + +function annotateLevelStd(symbols_for_lvl, avg_for_lvl, xAxisUs, configuredXLimits) +xLimits = getXLimits(configuredXLimits, xAxisUs, symbols_for_lvl); +xText = xLimits(1) + 0.96*diff(xLimits); + +for levelIdx = 1:size(symbols_for_lvl, 1) + levelSamples = symbols_for_lvl(levelIdx, :); + levelStd = std(levelSamples, 0, 2, 'omitnan'); + if isnan(levelStd) + continue + end + + levelAverage = avg_for_lvl(levelIdx, :); + yText = median(levelAverage(~isnan(levelAverage)), 'omitnan'); + if isnan(yText) + yText = median(levelSamples(~isnan(levelSamples)), 'omitnan'); + end + + label = ['std = ', sprintf('%.3f', levelStd)]; + text(xText, yText, label, ... + 'Interpreter', 'none', ... + 'HorizontalAlignment', 'right', ... + 'VerticalAlignment', 'middle', ... + 'FontSize', 10, ... + 'BackgroundColor', 'w', ... + 'Margin', 2, ... + 'EdgeColor', [0.8 0.8 0.8]); end end + +function levelAverage = interpolateMissingLevelAverage(levelAverage) +validSamples = ~isnan(levelAverage); + +if nnz(validSamples) == 0 + return +elseif nnz(validSamples) == 1 + levelAverage(:) = levelAverage(validSamples); + return +end + +t = 1:numel(levelAverage); +levelAverage(~validSamples) = interp1(t(validSamples), levelAverage(validSamples), ... + t(~validSamples), 'linear', 'extrap'); +end + +function values = numericSignal(signalLike) +if isa(signalLike, "Signal") + values = signalLike.signal; +else + values = signalLike; +end + +values = values(:).'; +end + +function values = normalizeNumericRms(values) +values = values ./ sqrt(mean(values.^2, "omitnan")); +end diff --git a/Functions/EQ_visuals/showMpiLevelScatter.m b/Functions/EQ_visuals/showMpiLevelScatter.m new file mode 100644 index 0000000..e999f58 --- /dev/null +++ b/Functions/EQ_visuals/showMpiLevelScatter.m @@ -0,0 +1,40 @@ +function [sep_sig, avg_sig, info] = showMpiLevelScatter(scopeInput, Symbols, options) +%SHOWMPILEVELSCATTER Backward-compatible wrapper around showLevelScatter. +% Prefer showLevelScatter directly for new code. This wrapper keeps older +% MPI call sites working while sharing one plotting/synchronization path. + +arguments + scopeInput + Symbols + options.fsym double = [] + options.syncFs (1,1) double = 0 + options.shiftFs (1,1) double = 0 + options.shifts double = [] + options.maxOccurences (1,1) double = Inf + options.normalize (1,1) logical = true + options.debug_plots (1,1) logical = false + options.fignum (1,1) double = NaN + options.clear (1,1) logical = true + options.showPlot (1,1) logical = true + options.xLimits double = [] + options.yLimits (1,2) double = [-3 3] + options.showStdAnnotations (1,1) logical = true + options.windowLength (1,1) double {mustBePositive, mustBeInteger} = 500 +end + +[sep_sig, avg_sig, info] = showLevelScatter(scopeInput, Symbols, ... + "fsym", options.fsym, ... + "syncFs", options.syncFs, ... + "shiftFs", options.shiftFs, ... + "shifts", options.shifts, ... + "maxOccurences", options.maxOccurences, ... + "normalize", options.normalize, ... + "debug_plots", options.debug_plots, ... + "fignum", options.fignum, ... + "clear", options.clear, ... + "showPlot", options.showPlot, ... + "xLimits", options.xLimits, ... + "yLimits", options.yLimits, ... + "showStdAnnotations", options.showStdAnnotations, ... + "windowLength", options.windowLength); +end diff --git a/Functions/EQ_blocks/dsp_runid.m b/Functions/Job_Processing/dsp_runid.m similarity index 96% rename from Functions/EQ_blocks/dsp_runid.m rename to Functions/Job_Processing/dsp_runid.m index 4561ef9..498eaed 100644 --- a/Functions/EQ_blocks/dsp_runid.m +++ b/Functions/Job_Processing/dsp_runid.m @@ -4,6 +4,7 @@ arguments run_id options.append_to_db = 0; options.max_occurences = 4; + options.start_occurence = 1; options.userParameters = struct(); options.database_type options.dataBase @@ -36,8 +37,8 @@ try dspInput = loadDspInputFromFilePaths(run_id, options); end - options.max_occurences = min(options.max_occurences, length(dspInput.Scpe_cell)); - for r = 1:options.max_occurences + num_occurences = length(dspInput.Scpe_cell); + for r = 1:num_occurences %%%%%%%% CORE EQUALIZATION CALL (Scpe, Symbols, Bits, 'Options') %%%%%%% diff --git a/Functions/Job_Processing/loadAndSyncRunSignals.m b/Functions/Job_Processing/loadAndSyncRunSignals.m new file mode 100644 index 0000000..2c6e633 --- /dev/null +++ b/Functions/Job_Processing/loadAndSyncRunSignals.m @@ -0,0 +1,196 @@ +function [Bits, Symbols, Scpe_cell, found_sync] = loadAndSyncRunSignals(dataTable, options) +%LOADANDSYNCRUNSIGNALS Load and synchronize signal files for one run. +% +% Inputs: +% dataTable - one-row table with run metadata and signal file paths +% options - struct with storage_path, start_occurence and max_occurences +% +% Outputs: +% Bits - transmitted bit reference +% Symbols - transmitted symbol reference +% Scpe_cell - synchronized received signal occurrences +% found_sync - true when a valid synchronization was found + +found_sync = 0; +tempLocalStorage = 1; +Scpe_cell = {}; +loaded_from_cache = false; + +storage_dir = fullfile(prefdir, 'temp_sync_data'); + +if tempLocalStorage == 1 + local_filename = fullfile(storage_dir, sprintf('sync_data_run_%s.mat', num2str(dataTable.run_id))); + if exist(local_filename, 'file') + try + cacheData = load(local_filename, 'Bits', 'Symbols', 'Scpe_cell'); + if isValidSyncCache(cacheData) + Bits = cacheData.Bits; + Symbols = cacheData.Symbols; + Scpe_cell = cacheData.Scpe_cell; + found_sync = 1; + loaded_from_cache = true; + else + warning('loadAndSyncRunSignals:InvalidSyncCache', ... + 'Ignoring incomplete sync cache file "%s".', local_filename); + safeDelete(local_filename); + end + catch + safeDelete(local_filename); + found_sync = 0; + end + end +end + +if ~found_sync + Bits = load(composeStoragePath(options.storage_path, dataTable.tx_bits_path)); + Bits = Bits.Bits; + + M = double(dataTable.pam_level); + fsym = dataTable.symbolrate; + Symbols_mapped = PAMmapper(M,0).map(Bits); + Symbols_mapped.fs = fsym; + + Symbols = load(composeStoragePath(options.storage_path, dataTable.tx_symbols_path)); + Symbols = Symbols.Symbols; + + found_sync = 0; + try + Scpe_load = load(composeStoragePath(options.storage_path, dataTable.rx_sync_path)); + Scpe_cell = Scpe_load.S; + [~,~,~,found_sync] = Scpe_cell{1}.tsynch("reference", Symbols, ... + "fs_ref", fsym, ... + "debug_plots", 0); + catch + % Continue with raw data if pre-synchronized data is unavailable. + end +end + +if ~found_sync + try + Scpe_sig_raw = load(composeStoragePath(options.storage_path, dataTable.rx_raw_path)); + 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", fsym, ... + "debug_plots", 0); + catch + % Continue to mapped-symbol fallback if raw data sync fails. + end +end + +if ~found_sync && exist('Scpe_sig_raw', 'var') + if length(Symbols_mapped.signal) ~= sum(Symbols_mapped.signal == Symbols.signal) + [~, Scpe_cell, ~, found_sync] = Scpe_sig_raw.tsynch("reference", Symbols_mapped, ... + "fs_ref", fsym, ... + "debug_plots", 0); + end +end + +if tempLocalStorage == 1 && found_sync && ~loaded_from_cache + if ~exist(storage_dir, 'dir') + mkdir(storage_dir); + end + + max_local_files = 10; + files = dir(fullfile(storage_dir, 'sync_data_run_*.mat')); + if length(files) >= max_local_files + [~, fileOrder] = sort([files.datenum]); + delete(fullfile(storage_dir, files(fileOrder(1)).name)); + end + + writeSyncCache(local_filename, Bits, Symbols, Scpe_cell); +end + +if found_sync + Scpe_cell = selectSyncedOccurrences(Scpe_cell, options); +else + warning('Could not synchronize the received signal with the stored symbols!'); +end + +end + +function path = composeStoragePath(storagePath, relativePath) +relativePath = firstValue(relativePath); +path = char(string(storagePath) + string(relativePath)); +end + +function value = firstValue(value) +if iscell(value) + value = value{1}; +elseif ~ischar(value) && ~isscalar(value) + value = value(1); +end +end + +function Scpe_cell = selectSyncedOccurrences(Scpe_cell, options) +available_occurences = length(Scpe_cell); +start_occurence = floor(getOption(options, 'start_occurence', 1)); +max_occurences = floor(getOption(options, 'max_occurences', available_occurences)); + +if available_occurences < 1 + warning('loadAndSyncRunSignals:NoSyncedOccurrences', ... + 'Synchronization reported success, but no synced occurrences are available.'); + return +end + +if start_occurence < 1 + error('loadAndSyncRunSignals:InvalidStartOccurrence', ... + 'start_occurence must be >= 1.'); +end + +if max_occurences < 1 + Scpe_cell = Scpe_cell(1:0); + return +end + +if start_occurence > available_occurences + warning('loadAndSyncRunSignals:StartOccurrenceTooHigh', ... + ['Requested start_occurence %d, but only %d synced occurrences are available. ', ... + 'Processing the last occurrence only.'], ... + start_occurence, available_occurences); + Scpe_cell = Scpe_cell(available_occurences); + return +end + +stop_occurence = min(available_occurences, start_occurence + max_occurences - 1); +Scpe_cell = Scpe_cell(start_occurence:stop_occurence); +end + +function value = getOption(options, name, defaultValue) +if isfield(options, name) + value = options.(name); +else + value = defaultValue; +end +end + +function valid = isValidSyncCache(cacheData) +valid = isfield(cacheData, 'Bits') && ... + isfield(cacheData, 'Symbols') && ... + isfield(cacheData, 'Scpe_cell') && ... + iscell(cacheData.Scpe_cell) && ... + ~isempty(cacheData.Scpe_cell); +end + +function writeSyncCache(local_filename, Bits, Symbols, Scpe_cell) +cache_dir = fileparts(local_filename); +temp_filename = [tempname(cache_dir), '.mat']; + +try + save(temp_filename, 'Bits', 'Symbols', 'Scpe_cell'); + movefile(temp_filename, local_filename, 'f'); +catch ME + safeDelete(temp_filename); + rethrow(ME); +end +end + +function safeDelete(filename) +if exist(filename, 'file') + try + delete(filename); + catch + % Another parallel worker may already have removed or replaced it. + end +end +end diff --git a/Functions/Job_Processing/loadAndSyncSignalDataFromDb.m b/Functions/Job_Processing/loadAndSyncSignalDataFromDb.m deleted file mode 100644 index c10f1d7..0000000 --- a/Functions/Job_Processing/loadAndSyncSignalDataFromDb.m +++ /dev/null @@ -1,114 +0,0 @@ -function [Bits, Symbols, Scpe_cell, found_sync] = loadAndSyncSignalDataFromDb(dataTable, options) -% LOADSIGNALDATA Loads and synchronizes signal data from storage -% -% Inputs:d -% dataTable - Table with file paths and configuration -% options - Struct with storage_path and max_occurences -% -% Outputs: -% Symbols_mapped - Mapped symbols from bits -% Symbols - Original symbols -% Scpe_cell - Cell array of synchronized signals -% found_sync - Boolean indicating if synchronization was successful - -found_sync = 0; -tempLocalStorage = 1; - -% Define the fixed storage directory relative to the user's MATLAB preferences directory -storage_dir = fullfile(prefdir, 'temp_sync_data'); - -% Part A: Check and load from local storage if available- -if tempLocalStorage == 1 - local_filename = fullfile(storage_dir, sprintf('sync_data_run_%s.mat', num2str(dataTable.run_id))); - if exist(local_filename, 'file') - % Load from local storage and return - try - load(local_filename, 'Bits', 'Symbols', 'Scpe_cell'); - found_sync = 1; - return - catch - delete(local_filename); - end - end -end - -% If not locally saved, load from storage -if ~found_sync - % Load transmitted bits - Bits = load(fullfile([options.storage_path, char(dataTable.tx_bits_path)])); - Bits = Bits.Bits; - - % Map bits to symbols - M = double(dataTable.pam_level); - fsym = dataTable.symbolrate; - Symbols_mapped = PAMmapper(M,0).map(Bits); - Symbols_mapped.fs = fsym; - - % Load original symbols - Symbols = load(fullfile([options.storage_path, char(dataTable.tx_symbols_path)])); - Symbols = Symbols.Symbols; - - found_sync = 0; - Scpe_cell = {}; - - % Try to load pre-synchronized data - try - Scpe_load = load(fullfile([options.storage_path, char(dataTable.rx_sync_path)])); - Scpe_cell = Scpe_load.S; - [~,~,~,found_sync] = Scpe_cell{2}.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1); - catch - % Continue to next method if this fails - end -end - -% If not found, try with raw data -if ~found_sync - try - 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", fsym, "debug_plots", 1); - catch - % Continue to next method if this fails - end -end - -% Last attempt with mapped symbols -if ~found_sync && exist('Scpe_sig_raw', 'var') - if length(Symbols_mapped.signal) ~= sum(Symbols_mapped.signal == Symbols.signal) - [~, Scpe_cell, ~, found_sync] = Scpe_sig_raw.tsynch("reference", Symbols_mapped, "fs_ref", fsym, "debug_plots", 0); - end -end - -% Part B: Save to local storage if data was loaded and synced -if tempLocalStorage == 1 && found_sync - % Create directory if it doesn't exist - if ~exist(storage_dir, 'dir') - mkdir(storage_dir); - end - -% local_filename = fullfile(storage_dir, sprintf('sync_data_run_%s.mat', num2str(dataTable.run_id))); - - % List existing files and remove oldest if more than N - max_local_files = 10; % Store up to N files - files = dir(fullfile(storage_dir, 'sync_data_run_*.mat')); - if length(files) >= max_local_files - % Sort by date - [~, idx] = sort([files.datenum]); - % Delete oldest file - delete(fullfile(storage_dir, files(idx(1)).name)); - end - - % Save current data - save(local_filename, 'Bits', 'Symbols', 'Scpe_cell'); -end - -% Limit number of occurrences -if found_sync - record_realizations = min(options.max_occurences, length(Scpe_cell)); - Scpe_cell = Scpe_cell(1:record_realizations); -else - warning('Could not synchronize the received signal with the stored symbols!'); -end - -end \ No newline at end of file diff --git a/Functions/Job_Processing/loadDspInputFromFilePaths.m b/Functions/Job_Processing/loadDspInputFromFilePaths.m new file mode 100644 index 0000000..931f3bb --- /dev/null +++ b/Functions/Job_Processing/loadDspInputFromFilePaths.m @@ -0,0 +1,97 @@ +function dspInput = loadDspInputFromFilePaths(run_id, options) +%LOADDSPINPUTFROMFILEPATHS Load explicit signal files and prepare DSP input. + +arguments + run_id + options struct +end + +Tx_bits = load(textScalar(options.load_file_path.tx_bits_path)); +Symbols = load(textScalar(options.load_file_path.tx_symbols_path)); +Scpe_sig_raw = load(textScalar(options.load_file_path.rx_raw_path)); + +Tx_bits = Tx_bits.Bits; +Symbols = Symbols.Symbols; +Scpe_sig_raw = Scpe_sig_raw.Scpe_sig_raw; + +fsym = Symbols.fs; +M = Symbols.logbook.ModifierCopy{1}.M; +duob_mode = Symbols.logbook.ModifierCopy{1}.duobinary_mode; + +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", fsym, ... + "debug_plots", 1); + +if found_sync + Scpe_cell = selectSyncedOccurrences(Scpe_cell, options); +else + warning('Could not synchronize the received signal with the stored symbols!'); +end + +dspInput = struct(); +dspInput.run_id = run_id; +dspInput.dataTable = table(); +dspInput.Tx_bits = Tx_bits; +dspInput.Symbols = Symbols; +dspInput.Scpe_cell = Scpe_cell; +dspInput.found_sync = found_sync; +dspInput.fsym = fsym; +dspInput.M = M; +dspInput.duob_mode = duob_mode; +end + +function value = textScalar(value) +value = firstValue(value); +value = char(string(value)); +end + +function value = firstValue(value) +if iscell(value) + value = value{1}; +elseif ~ischar(value) && ~isscalar(value) + value = value(1); +end +end + +function Scpe_cell = selectSyncedOccurrences(Scpe_cell, options) +available_occurences = length(Scpe_cell); +start_occurence = floor(getOption(options, 'start_occurence', 1)); +max_occurences = floor(getOption(options, 'max_occurences', available_occurences)); + +if available_occurences < 1 + warning('loadDspInputFromFilePaths:NoSyncedOccurrences', ... + 'Synchronization reported success, but no synced occurrences are available.'); + return +end + +if start_occurence < 1 + error('loadDspInputFromFilePaths:InvalidStartOccurrence', ... + 'start_occurence must be >= 1.'); +end + +if max_occurences < 1 + Scpe_cell = Scpe_cell(1:0); + return +end + +if start_occurence > available_occurences + warning('loadDspInputFromFilePaths:StartOccurrenceTooHigh', ... + ['Requested start_occurence %d, but only %d synced occurrences are available. ', ... + 'Processing the last occurrence only.'], ... + start_occurence, available_occurences); + Scpe_cell = Scpe_cell(available_occurences); + return +end + +stop_occurence = min(available_occurences, start_occurence + max_occurences - 1); +Scpe_cell = Scpe_cell(start_occurence:stop_occurence); +end + +function value = getOption(options, name, defaultValue) +if isfield(options, name) + value = options.(name); +else + value = defaultValue; +end +end diff --git a/Functions/Job_Processing/loadDspInputFromRunId.m b/Functions/Job_Processing/loadDspInputFromRunId.m new file mode 100644 index 0000000..ed23272 --- /dev/null +++ b/Functions/Job_Processing/loadDspInputFromRunId.m @@ -0,0 +1,40 @@ +function dspInput = loadDspInputFromRunId(run_id, database, options) +%LOADDSPINPUTFROMRUNID Query run metadata and prepare canonical DSP input. + + arguments + run_id + database + options struct + end + + dataTable = queryRunid(run_id, database); + + % Load signal files referenced by the run metadata, verify/synchronize the + % received signal, optionally cache the sync result, and cap occurrences. + [Tx_bits, Symbols, Scpe_cell, found_sync] = loadAndSyncRunSignals(dataTable, options); + + dspInput = struct(); + dspInput.run_id = run_id; + dspInput.dataTable = dataTable; + dspInput.Tx_bits = Tx_bits; + dspInput.Symbols = Symbols; + dspInput.Scpe_cell = Scpe_cell; + dspInput.found_sync = found_sync; + dspInput.fsym = dataTable.symbolrate; + dspInput.M = double(dataTable.pam_level); + dspInput.duob_mode = parseDbMode(dataTable.db_mode); + +end + +function mode = parseDbMode(rawMode) +if isnumeric(rawMode) + mode = db_mode(rawMode); + return +end + +if iscell(rawMode) + rawMode = rawMode{1}; +end + +mode = db_mode(strrep(string(rawMode), '"', '')); +end diff --git a/Functions/Job_Processing/preprocessRunId.m b/Functions/Job_Processing/preprocessRunId.m deleted file mode 100644 index ba21d00..0000000 --- a/Functions/Job_Processing/preprocessRunId.m +++ /dev/null @@ -1,36 +0,0 @@ -function dspInput = preprocessRunId(run_id, database, options) -%PREPROCESSRUNID Load one run_id and prepare the canonical DSP input. - - arguments - run_id - database - options struct - end - - dataTable = queryRunid(run_id, database); - [Tx_bits, Symbols, Scpe_cell, found_sync] = loadAndSyncSignalDataFromDb(dataTable, options); - - dspInput = struct(); - dspInput.run_id = run_id; - dspInput.dataTable = dataTable; - dspInput.Tx_bits = Tx_bits; - dspInput.Symbols = Symbols; - dspInput.Scpe_cell = Scpe_cell; - dspInput.found_sync = found_sync; - dspInput.fsym = dataTable.symbolrate; - dspInput.M = double(dataTable.pam_level); - dspInput.duob_mode = parseDbMode(dataTable.db_mode); -end - -function mode = parseDbMode(rawMode) - if isnumeric(rawMode) - mode = db_mode(rawMode); - return - end - - if iscell(rawMode) - rawMode = rawMode{1}; - end - - mode = db_mode(strrep(string(rawMode), '"', '')); -end diff --git a/Functions/Job_Processing/runBatch.m b/Functions/Job_Processing/runBatch.m index 70d78d7..b37d675 100644 --- a/Functions/Job_Processing/runBatch.m +++ b/Functions/Job_Processing/runBatch.m @@ -1,6 +1,6 @@ -function results = runBatch(workerFcn, jobs, options) +function batchResults = runBatch(workerFcn, jobs, options) %RUNBATCH Execute a list of jobs in serial or parallel. -% results = runBatch(workerFcn, jobs) executes each jobs(i).args cell +% batchResults = runBatch(workerFcn, jobs) executes each jobs(i).args cell % with workerFcn and returns a cell array aligned with the input jobs. % % Supported job fields: @@ -24,7 +24,7 @@ function results = runBatch(workerFcn, jobs, options) mode = normalizeProcessingMode(options.mode); jobs = normalizeJobs(jobs); nJobs = numel(jobs); - results = cell(1, nJobs); + batchResults = cell(1, nJobs); h = []; if nJobs == 0 @@ -33,7 +33,7 @@ function results = runBatch(workerFcn, jobs, options) if options.waitbar h = waitbar(0, char(options.waitbarMessage)); - cleanupWaitbar = onCleanup(@() closeWaitbar(h)); %#ok + cleanupWaitbar = onCleanup(@() closeWaitbar(h)); end switch mode @@ -41,36 +41,36 @@ function results = runBatch(workerFcn, jobs, options) pool = setupParallelPool(options.numWorkers, options.idleTimeout, options.cancelExistingQueue); futures = parallel.FevalFuture.empty(nJobs, 0); - for idx = 1:nJobs - fprintf('[%s] Submitted.\n', jobs(idx).label); - futures(idx) = parfeval(pool, workerFcn, 1, jobs(idx).args{:}); + for jobIdx = 1:nJobs + fprintf('[%s] Submitted.\n', jobs(jobIdx).label); + futures(jobIdx) = parfeval(pool, workerFcn, 1, jobs(jobIdx).args{:}); end consumedIdx = false(nJobs, 1); for completedCount = 1:nJobs try - [idx, value] = fetchNext(futures); - consumedIdx(idx) = true; - results{idx} = value; - fprintf('[%s] Completed (%d/%d).\n', jobs(idx).label, completedCount, nJobs); + [jobIdx, jobResult] = fetchNext(futures); + consumedIdx(jobIdx) = true; + batchResults{jobIdx} = jobResult; + fprintf('[%s] Completed (%d/%d).\n', jobs(jobIdx).label, completedCount, nJobs); if ~isempty(options.resultHandler) - options.resultHandler(value, jobs(idx), idx); + options.resultHandler(jobResult, jobs(jobIdx), jobIdx); end catch fetchErr - idxErr = findErroredFuture(futures, consumedIdx); - if isempty(idxErr) + errorJobIdx = findErroredFuture(futures, consumedIdx); + if isempty(errorJobIdx) rethrow(fetchErr); end - consumedIdx(idxErr) = true; - errInfo = extractFutureError(futures(idxErr)); - results{idxErr} = errInfo; + consumedIdx(errorJobIdx) = true; + errInfo = extractFutureError(futures(errorJobIdx)); + batchResults{errorJobIdx} = errInfo; if ~isempty(options.errorHandler) - options.errorHandler(errInfo, jobs(idxErr), idxErr); + options.errorHandler(errInfo, jobs(errorJobIdx), errorJobIdx); else - defaultErrorHandler(errInfo, jobs(idxErr), idxErr); + defaultErrorHandler(errInfo, jobs(errorJobIdx), errorJobIdx); end end @@ -78,27 +78,27 @@ function results = runBatch(workerFcn, jobs, options) end case processingMode.serial - for idx = 1:nJobs + for jobIdx = 1:nJobs try - fprintf('[%s] Running.\n', jobs(idx).label); - value = feval(workerFcn, jobs(idx).args{:}); - results{idx} = value; - fprintf('[%s] Completed (%d/%d).\n', jobs(idx).label, idx, nJobs); + fprintf('[%s] Running.\n', jobs(jobIdx).label); + jobResult = workerFcn(jobs(jobIdx).args{:}); + batchResults{jobIdx} = jobResult; + fprintf('[%s] Completed (%d/%d).\n', jobs(jobIdx).label, jobIdx, nJobs); if ~isempty(options.resultHandler) - options.resultHandler(value, jobs(idx), idx); + options.resultHandler(jobResult, jobs(jobIdx), jobIdx); end catch ME - results{idx} = ME; + batchResults{jobIdx} = ME; if ~isempty(options.errorHandler) - options.errorHandler(ME, jobs(idx), idx); + options.errorHandler(ME, jobs(jobIdx), jobIdx); else - defaultErrorHandler(ME, jobs(idx), idx); + defaultErrorHandler(ME, jobs(jobIdx), jobIdx); end end - updateWaitbar(options.waitbar, h, idx, nJobs); + updateWaitbar(options.waitbar, h, jobIdx, nJobs); end otherwise @@ -128,23 +128,23 @@ function jobs = normalizeJobs(jobs) return end - for idx = 1:numel(jobs) - if ~isfield(jobs, 'args') || isempty(jobs(idx).args) - jobs(idx).args = {}; + for jobIdx = 1:numel(jobs) + if ~isfield(jobs, 'args') || isempty(jobs(jobIdx).args) + jobs(jobIdx).args = {}; end - if ~iscell(jobs(idx).args) - error('runBatch:InvalidArgs', 'jobs(%d).args must be a cell array.', idx); + if ~iscell(jobs(jobIdx).args) + error('runBatch:InvalidArgs', 'jobs(%d).args must be a cell array.', jobIdx); end - if ~isfield(jobs, 'label') || isempty(jobs(idx).label) - jobs(idx).label = sprintf('Job %d', idx); + if ~isfield(jobs, 'label') || isempty(jobs(jobIdx).label) + jobs(jobIdx).label = sprintf('Job %d', jobIdx); else - jobs(idx).label = string(jobs(idx).label); + jobs(jobIdx).label = string(jobs(jobIdx).label); end if ~isfield(jobs, 'meta') - jobs(idx).meta = struct(); + jobs(jobIdx).meta = struct(); end end end @@ -189,9 +189,9 @@ function pool = setupParallelPool(numWorkers, idleTimeout, cancelExistingQueue) end end -function idxErr = findErroredFuture(futures, consumedIdx) +function errorJobIdx = findErroredFuture(futures, consumedIdx) readMask = arrayfun(@(future) future.Read, futures).'; - idxErr = find(readMask & ~consumedIdx, 1); + errorJobIdx = find(readMask & ~consumedIdx, 1); end function errInfo = extractFutureError(future) diff --git a/Functions/Job_Processing/submitJobs.m b/Functions/Job_Processing/submitJobs.m index a11f057..de48776 100644 --- a/Functions/Job_Processing/submitJobs.m +++ b/Functions/Job_Processing/submitJobs.m @@ -9,9 +9,11 @@ function [results, wh] = submitJobs(run_ids, dsp_options, submit_mode, submit_op % submit_options : struct with fields % .waitbar (logical) % .wh (DataStorage object) +% .storePackages (string array, optional package whitelist) % % results : cell(nJobsPerRunId, nRunIds) -% wh : updated DataStorage +% wh : updated DataStorage. For multiple run_ids, run_id is added as +% the first storage axis. arguments run_ids int32 = 0 @@ -19,15 +21,19 @@ arguments submit_mode processingMode = processingMode.serial submit_options.waitbar (1,1) logical = true submit_options.wh = DataStorage(struct()) + submit_options.storePackages string = string.empty() end % Normalize run_ids = run_ids(:)'; nRunIds = numel(run_ids); -nJobsPerRunId = submit_options.wh.getLastLinIndice(); +sweepWh = submit_options.wh; +validateNoRunIdSweepParameter(sweepWh); + +nJobsPerRunId = sweepWh.getLastLinIndice(); totalJobs = nRunIds * nJobsPerRunId; -wh = submit_options.wh; +wh = buildStorageWarehouse(sweepWh, run_ids); results = cell(nJobsPerRunId, nRunIds); jobs = repmat(struct('args', {{}}, 'label', "", 'meta', struct()), 1, totalJobs); @@ -35,14 +41,16 @@ jobCounter = 0; for r = 1:nRunIds for k = 1:nJobsPerRunId jobCounter = jobCounter + 1; - optionalVars = buildOptionalVars(k, submit_options.wh); + userParameters = buildUserParameters(k, sweepWh); jobOptions = dsp_options; - jobOptions.parameters = optionalVars; + jobOptions.userParameters = userParameters; 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; + jobs(jobCounter).meta.sweepJobIndex = k; + jobs(jobCounter).meta.storageJobIndex = buildStorageJobIndex(run_ids(r), userParameters, wh); jobs(jobCounter).meta.run_id = run_ids(r); end end @@ -74,16 +82,69 @@ end fprintf('Full report:\n%s\n', getReport(ME,'extended')); end - function optionalVars = buildOptionalVars(jobIndex, wh) - optionalVars = struct(); + function userParameters = buildUserParameters(jobIndex, wh) + userParameters = struct(); if ~isempty(wh.getDimension()) [vals, names] = wh.getPhysIndicesByLinIndex(jobIndex); for pi = 1:numel(names) - optionalVars.(names{pi}) = vals{pi}; + userParameters.(names{pi}) = vals{pi}; end end end + function validateNoRunIdSweepParameter(wh) + if isfield(wh.inputParams, "run_id") + error("submitJobs:RunIdParameterConflict", ... + "Do not define userParameters.run_id. submitJobs manages run_id as a storage axis."); + end + end + + function storageWh = buildStorageWarehouse(sweepWh, runIds) + if isscalar(runIds) + storageWh = sweepWh; + return + end + + storageParameters = struct(); + storageParameters.run_id = runIds; + + sweepParameterNames = fieldnames(sweepWh.inputParams); + for parameterIdx = 1:numel(sweepParameterNames) + parameterName = sweepParameterNames{parameterIdx}; + storageParameters.(parameterName) = sweepWh.inputParams.(parameterName); + end + + storageWh = DataStorage(storageParameters); + end + + function storageJobIndex = buildStorageJobIndex(runId, userParameters, wh) + storageParameterNames = wh.fn; + + if isempty(storageParameterNames) + storageJobIndex = 1; + return + end + + storageSubscripts = cell(1, numel(storageParameterNames)); + for parameterIdx = 1:numel(storageParameterNames) + parameterName = char(storageParameterNames(parameterIdx)); + + if strcmp(parameterName, "run_id") + parameterValue = runId; + else + parameterValue = userParameters.(parameterName); + end + + storageSubscripts{parameterIdx} = wh.getIndexByPhys(parameterName, parameterValue); + end + + if isscalar(storageSubscripts) + storageJobIndex = storageSubscripts{1}; + else + storageJobIndex = sub2ind(wh.getStorageSize(), storageSubscripts{:}); + end + end + function nameValue = structToNameValue(options) names = fieldnames(options); nameValue = cell(1, 2*numel(names)); @@ -95,14 +156,41 @@ end end function storeResult(val, job, ~) - if ~isempty(wh) - jobIndex = job.meta.jobIndex; - wh.addValueToStorageByLinIdx(val.ffe_package, 'ffe_package', jobIndex); - wh.addValueToStorageByLinIdx(val.mlse_package, 'mlse_package', jobIndex); - wh.addValueToStorageByLinIdx(val.vnle_package, 'vnle_package', jobIndex); - wh.addValueToStorageByLinIdx(val.dbtgt_package,'dbtgt_package',jobIndex); - wh.addValueToStorageByLinIdx(val.dbenc_package,'dbenc_package',jobIndex); - wh.addValueToStorageByLinIdx(val.mlmlse_package,'mlmlse_package',jobIndex); + if isempty(wh) || ~isstruct(val) + return + end + + storageJobIndex = job.meta.storageJobIndex; + resultFields = fieldnames(val); + + for resultIdx = 1:numel(resultFields) + storageName = resultFields{resultIdx}; + + if ~shouldStore(storageName) + continue + end + + if isempty(val.(storageName)) + continue + end + + ensureStorage(storageName); + wh.addValueToStorageByLinIdx(val.(storageName), storageName, storageJobIndex); + end + end + + function tf = shouldStore(storageName) + if isempty(submit_options.storePackages) + tf = true; + return + end + + tf = any(string(storageName) == submit_options.storePackages); + end + + function ensureStorage(storageName) + if ~isfield(wh.sto, storageName) + wh.addStorage(storageName); end end diff --git a/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/FIGURE_EYES.m b/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/FIGURE_EYES.m index 0d953af..626862c 100644 --- a/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/FIGURE_EYES.m +++ b/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/FIGURE_EYES.m @@ -31,7 +31,7 @@ M = double(dataTable.pam_level); duob_mode = db_mode(strrep(dataTable.db_mode,'"','')); % Load and Sync signal data from DB -[Tx_bits, Symbols, Scpe_cell, ~] = loadAndSyncSignalDataFromDb(dataTable, dsp_options); +[Tx_bits, Symbols, Scpe_cell, ~] = loadAndSyncRunSignals(dataTable, dsp_options); Scpe_sig_syncd = Scpe_cell{1}; Scpe_sig_syncd.eye(fsym,M,"fignum",rate.*1e-9*M+1,"displayname",' Eye of Signal'); @@ -276,4 +276,4 @@ ax1.Position = pos1; ax2.Box = 'on'; % keep outline but no ticks on the right ax1.Box = 'on'; -ax2.View = [90 -90]; \ No newline at end of file +ax2.View = [90 -90]; 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 7762768..5895cf3 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 @@ -82,7 +82,7 @@ for dbmode = 0%length(rates) M = double(dataTable.pam_level); % Load and Sync signal data from DB - [Tx_bits, Symbols, Scpe_cell, ~] = loadAndSyncSignalDataFromDb(dataTable, dsp_options); + [Tx_bits, Symbols, Scpe_cell, ~] = loadAndSyncRunSignals(dataTable, dsp_options); % Preprocess signal Scpe_sig = preprocessSignal(Scpe_cell{1}, Symbols, fsym); 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 9fd31ed..9e413af 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 @@ -30,7 +30,7 @@ fsym = dataTable.symbolrate; M = double(dataTable.pam_level); % Load and Sync signal data from DB -[Tx_bits, Symbols, Scpe_cell, ~] = loadAndSyncSignalDataFromDb(dataTable, dsp_options); +[Tx_bits, Symbols, Scpe_cell, ~] = loadAndSyncRunSignals(dataTable, dsp_options); % Preprocess signal Scpe_sig = preprocessSignal(Scpe_cell{1}, Symbols, fsym); @@ -62,7 +62,7 @@ fsym = dataTable.symbolrate; M = double(dataTable.pam_level); % Load and Sync signal data from DB -[Tx_bits, Symbols, Scpe_cell, ~] = loadAndSyncSignalDataFromDb(dataTable, dsp_options); +[Tx_bits, Symbols, Scpe_cell, ~] = loadAndSyncRunSignals(dataTable, dsp_options); % Preprocess signal Scpe_sig = preprocessSignal(Scpe_cell{1}, Symbols, fsym); @@ -117,4 +117,4 @@ matlab2tikz(outfile, ... 'legend columns=1' ... }); -end \ No newline at end of file +end 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 ceaa0a2..7ddd071 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 @@ -3,10 +3,10 @@ dsp_options.append_to_db = 0; dsp_options.max_occurences = 1; experiment = "highspeed_2024"; -dsp_options.mode = "load_run_id"; % 'simulate' & 'load_files' +dsp_options.mode = "run_id"; % "run_id", "file_paths", or future "simulate" dsp_options.load_file_path = struct(); -if dsp_options.mode == "load_run_id" +if dsp_options.mode == "run_id" if experiment == "highspeed_2024" @@ -41,7 +41,7 @@ if dsp_options.mode == "load_run_id" end -elseif dsp_options.mode == "load_files" +elseif dsp_options.mode == "file_paths" dsp_options.load_file_path.tx_bits_path = "Z:\2025\ECOC Silas\ecoc_2025\mpi_opti_1000m_pam2 4 6 8\20250417_091513_PAM_4_R_112_bits.mat"'; dsp_options.load_file_path.tx_symbols_path = "Z:\2025\ECOC Silas\ecoc_2025\mpi_opti_1000m_pam2 4 6 8\20250417_091513_PAM_4_R_112_symbols.mat"'; @@ -74,10 +74,10 @@ fp.where('Runs', 'rop_attenuation','EQUAL', 0); [dataTable,~] = db.queryDB(fp, db.getTableFieldNames('Runs')); % === Set LOOPS & Initialize DataStorage === -dsp_options.parameters = struct(); -% dsp_options.parameters.pf_ncoeffs = [1,2];%s[0,logspace(-4,0,10)]; +dsp_options.userParameters = struct(); +dsp_options.userParameters.pf_ncoeffs = [1,2];%s[0,logspace(-4,0,10)]; -wh = DataStorage(dsp_options.parameters); +wh = DataStorage(dsp_options.userParameters); wh.addStorage("ffe_package"); wh.addStorage("mlse_package"); wh.addStorage("vnle_package"); diff --git a/projects/Advanced_DSP_for_400G_IMDD_experiments/a_minimal_example.m b/projects/Advanced_DSP_for_400G_IMDD_experiments/a_minimal_example.m index 7cd7610..bf2bf16 100644 --- a/projects/Advanced_DSP_for_400G_IMDD_experiments/a_minimal_example.m +++ b/projects/Advanced_DSP_for_400G_IMDD_experiments/a_minimal_example.m @@ -34,7 +34,7 @@ for id = run_ids duob_mode = db_mode(strrep(dataTable(1,:).db_mode,'"','')); % Load and Sync signal data from DB - [Tx_bits, Symbols, Scpe_cell, ~] = loadAndSyncSignalDataFromDb(dataTable(1,:), dsp_options); + [Tx_bits, Symbols, Scpe_cell, ~] = loadAndSyncRunSignals(dataTable(1,:), dsp_options); % Preprocess signal Scpe_sig = preprocessSignal(Scpe_cell{1}, Symbols, fsym); @@ -159,4 +159,4 @@ ml_mlse_equalizer = ML_MLSE("epochs_tr",training_epochs,"epochs_dd",1, ... [ml_mlse_results] = ml_mlse(ml_mlse_equalizer, M, Scpe_sig, Symbols, Tx_bits,"precode_mode",duob_mode); -ml_mlse_results.metrics.print("description",'ML pre Eq. + Viterbi') \ No newline at end of file +ml_mlse_results.metrics.print("description",'ML pre Eq. + Viterbi') diff --git a/projects/Advanced_DSP_for_400G_IMDD_experiments/auswertung MPI/dsp_run_id.m b/projects/Advanced_DSP_for_400G_IMDD_experiments/auswertung MPI/dsp_run_id.m index 816505d..9be3a19 100644 --- a/projects/Advanced_DSP_for_400G_IMDD_experiments/auswertung MPI/dsp_run_id.m +++ b/projects/Advanced_DSP_for_400G_IMDD_experiments/auswertung MPI/dsp_run_id.m @@ -4,7 +4,7 @@ arguments run_id options.append_to_db = 0; options.max_occurences = 4; - options.parameters = struct(); + options.userParameters = struct(); end basePath = 'C:\Users\Silas\Documents\MATLAB\Datensätze\sioe_labor\'; @@ -46,8 +46,8 @@ 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; +% Overwrite default parameters if given in options.userParameters +paramStruct = options.userParameters; if ~isempty(paramStruct) paramNames = fieldnames(paramStruct); for i = 1:numel(paramNames) @@ -191,4 +191,4 @@ output.vnle_dfe_package = vnle_dfe_package; output.vnle_pf_package = vnle_pf_package; output.dbtgt_package = dbtgt_package; -end \ No newline at end of file +end diff --git a/projects/Diss/MPI_revisit/investigate_mpi_algorithms.m b/projects/Diss/MPI_revisit/investigate_mpi_algorithms.m new file mode 100644 index 0000000..04fd076 --- /dev/null +++ b/projects/Diss/MPI_revisit/investigate_mpi_algorithms.m @@ -0,0 +1,76 @@ +% === DSP settings === +dsp_options = struct(); +dsp_options.mode = "run_id"; +dsp_options.recipe = @mpi_recipe_dev; +dsp_options.append_to_db = false; +dsp_options.max_occurences = 3; +dsp_options.start_occurence = 1; +dsp_options.debug_plots = true; + +dsp_options.database_type = "mysql"; + +dsp_options.dataBase = "labor"; +dsp_options.storage_path = "W:\labdata\ECOC Silas\ecoc_2025"; + +dsp_options.server = "192.168.178.192"; +dsp_options.port = 3306; +dsp_options.user = "silas"; +dsp_options.password = "silas"; + +db = DBHandler("dataBase", [dsp_options.dataBase],... + "type", dsp_options.database_type,... + "server", dsp_options.server,... + "user", dsp_options.user, "password", dsp_options.password); + +%% +fp = QueryFilter(); +fp.where('Runs', 'symbolrate', 'EQUALS', 112e9); +fp.where('Runs', 'fiber_length', 'EQUALS', 0); +fp.where('Runs', 'interference_path_length', 'EQUALS', 1000); +fp.where('Runs', 'sir', 'EQUALS', 23); +% fp.where('Runs', 'db_mode', 'EQUALS', '"no_db"'); +% fp.where('Runs', 'is_mpi', 'EQUALS', 1); +fp.where('Runs', 'pam_level', 'EQUALS', 4); +fp.where('Runs', 'wavelength', 'EQUALS', 1310); +fp.where('Runs', 'v_bias', 'EQUALS', 2.65); + +[dataTable,~] = db.queryDB(fp, db.getTableFieldNames('Runs')); + +[~, uniqueSirRows] = unique(dataTable.sir, "stable"); +dataTable = dataTable(uniqueSirRows, :); +% sort rows by sir (ascending) and extract run_ids +[~, sortIdx] = sort(dataTable.sir, 'descend'); +dataTable = dataTable(sortIdx, :); +run_ids = dataTable.run_id; + +if isempty(run_ids) + error("run_minimal_recipe:MissingRunIds", ... + "Set run_ids to one or more known run IDs before running this example."); +end + +%% === Warehouse setup === +dsp_options.userParameters = struct(); +dsp_options.userParameters.smoothing_length = logspace(1,5.5,5); + +wh = DataStorage(dsp_options.userParameters); + +%% === Run === +% submitJobs returns the raw per-job results and the filled Warehouse. For a +% single run_id and remove_dc = [0, 1], results is a 2-by-1 cell array. +[results, wh] = submitJobs(run_ids, dsp_options, processingMode.parallel, ... + "wh", wh, ... + "waitbar", true); + +%% Analyze + +storageNames = fieldnames(wh.sto); +x_vars = dsp_options.userParameters.smoothing_length; + +figure(2026);hold on +for i = 1:numel(run_ids) + disp(run_ids(i)) + result = wh.getStoValue(storageNames{1},x_vars); + ber = cellfun(@(packageCell) packageCell{1}.metrics.BER, result); + plot(x_vars,ber); + beautifyBERplot("logscale",true,"setcolors",true,"setmarkers",true); +end diff --git a/projects/ECOC_2025_MPI/auswertung_algorithms/run_offline_dsp.m b/projects/ECOC_2025_MPI/auswertung_algorithms/run_offline_dsp.m index a61927e..5fe9939 100644 --- a/projects/ECOC_2025_MPI/auswertung_algorithms/run_offline_dsp.m +++ b/projects/ECOC_2025_MPI/auswertung_algorithms/run_offline_dsp.m @@ -7,8 +7,8 @@ % dsp_options.database_name = 'silas_labor_newdsp_newstructure.db'; % dsp_options.storage_path = 'Z:\2024\sioe_labor\'; % -% dsp_options.parameters = struct(); -% dsp_options.parameters.mu_dc = [0.005]; +% dsp_options.userParameters = struct(); +% dsp_options.userParameters.mu_dc = [0.005]; % % % === Get Run ID's === % db = DBHandler("pathToDB", [dsp_options.database_path, dsp_options.database_name], "type", "sqlite"); @@ -21,6 +21,7 @@ dsp_options.storage_path = 'W:\labdata\ECOC Silas\ECOC Silas\ecoc_2025\'; dsp_options.server = "192.168.178.192";% "134.245.243.254"; dsp_options.user = "silas"; dsp_options.password = "silas"; +dsp_options.userParameters = struct(); db = DBHandler("dataBase", [dsp_options.dataBase],... "type", dsp_options.database_type,... @@ -40,7 +41,7 @@ fp.where('Runs', 'bitrate','EQUALS', 360e9); [dataTable,~] = db.queryDB(fp, db.getTableFieldNames('Runs')); % === Initialize DataStorage === -wh = DataStorage(dsp_options.parameters); +wh = DataStorage(dsp_options.userParameters); wh.addStorage("ffe_package"); wh.addStorage("mlse_package"); wh.addStorage("vnle_package"); @@ -131,4 +132,4 @@ ylabel(y_var, 'Interpreter', 'none'); legend('show', 'Location', 'best'); grid on; title(sprintf('%s vs. %s', y_var, x_var), 'Interpreter', 'none'); -hold off; \ No newline at end of file +hold off; diff --git a/projects/ECOC_2025_MPI/dsp_run_id.m b/projects/ECOC_2025_MPI/dsp_run_id.m index 5f58f5b..b238727 100644 --- a/projects/ECOC_2025_MPI/dsp_run_id.m +++ b/projects/ECOC_2025_MPI/dsp_run_id.m @@ -4,7 +4,7 @@ arguments run_id options.append_to_db = 0; options.max_occurences = 4; - options.parameters = struct(); + options.userParameters = struct(); options.database_path options.database_name @@ -60,8 +60,8 @@ ffe_buffer_len = 1; smoothing_buffer_length = 4096; smoothing_buffer_update = 224; -% Overwrite default parameters if given in options.parameters -paramStruct = options.parameters; +% Overwrite default parameters if given in options.userParameters +paramStruct = options.userParameters; if ~isempty(paramStruct) paramNames = fieldnames(paramStruct); for i = 1:numel(paramNames) @@ -261,4 +261,4 @@ output.vnle_pf_package = vnle_pf_package; output.vnle_package = vnle_package; output.dbtgt_package = dbtgt_package; -end \ No newline at end of file +end diff --git a/projects/ECOC_2025_MPI/dsp_standalone.m b/projects/ECOC_2025_MPI/dsp_standalone.m index 68c36ea..2b1c7ad 100644 --- a/projects/ECOC_2025_MPI/dsp_standalone.m +++ b/projects/ECOC_2025_MPI/dsp_standalone.m @@ -52,7 +52,7 @@ params.mu_dc = 0.005; futures_list = parallel.FevalFuture.empty(); for id = 1:length(dataTable.run_id) run_id = dataTable.run_id(id); - [out, futures_list(id)] = submit_dsp(run_id, databasePath, database_name, savePath,"parallel",run_par,'max_occurences',num_occ,'paramstruct',params); + [out, futures_list(id)] = submit_dsp(run_id, databasePath, database_name, savePath,"parallel",run_par,'max_occurences',num_occ,'userParameters',params); end % ideal DC tracking @@ -66,7 +66,7 @@ params.mu_dc = 0.005; futures_list = parallel.FevalFuture.empty(); for id = 1:length(dataTable.run_id) run_id = dataTable.run_id(id); - [out, futures_list(id)] = submit_dsp(run_id, databasePath, database_name, savePath,"parallel",run_par,'max_occurences',num_occ,'paramstruct',params); + [out, futures_list(id)] = submit_dsp(run_id, databasePath, database_name, savePath,"parallel",run_par,'max_occurences',num_occ,'userParameters',params); end % DC smoothing @@ -79,7 +79,7 @@ params.mu_dc = 0.00; futures_list = parallel.FevalFuture.empty(); for id = 1:length(dataTable.run_id) run_id = dataTable.run_id(id); - [out, futures_list(id)] = submit_dsp(run_id, databasePath, database_name, savePath,"parallel",run_par,'max_occurences',num_occ,'paramstruct',params); + [out, futures_list(id)] = submit_dsp(run_id, databasePath, database_name, savePath,"parallel",run_par,'max_occurences',num_occ,'userParameters',params); end % only FFE @@ -92,7 +92,7 @@ params.mu_dc = 0.00; futures_list = parallel.FevalFuture.empty(); for id = 1:length(dataTable.run_id) run_id = dataTable.run_id(id); - [out, futures_list(id)] = submit_dsp(run_id, databasePath, database_name, savePath,"parallel",run_par,'max_occurences',num_occ,'paramstruct',params); + [out, futures_list(id)] = submit_dsp(run_id, databasePath, database_name, savePath,"parallel",run_par,'max_occurences',num_occ,'userParameters',params); end diff --git a/projects/ECOC_2025_MPI/dsp_standalone_2.m b/projects/ECOC_2025_MPI/dsp_standalone_2.m index 97007ae..f4baa55 100644 --- a/projects/ECOC_2025_MPI/dsp_standalone_2.m +++ b/projects/ECOC_2025_MPI/dsp_standalone_2.m @@ -45,7 +45,7 @@ params.mu_dc = 0.005; futures_list = parallel.FevalFuture.empty(); for id = 1:length(dataTable.run_id) run_id = dataTable.run_id(id); - [out, futures_list(id)] = submit_dsp(run_id, databasePath, database_name, savePath,"parallel",run_par,'max_occurences',num_occ,'paramstruct',params); + [out, futures_list(id)] = submit_dsp(run_id, databasePath, database_name, savePath,"parallel",run_par,'max_occurences',num_occ,'userParameters',params); end % Extract all ber_mlse values from the vnle_pf_package using cellfun diff --git a/projects/ECOC_2025_MPI/load_signal_standalone.m b/projects/ECOC_2025_MPI/load_signal_standalone.m index c8a801c..b27238a 100644 --- a/projects/ECOC_2025_MPI/load_signal_standalone.m +++ b/projects/ECOC_2025_MPI/load_signal_standalone.m @@ -45,46 +45,12 @@ for i = 1:size(dataTable,1) Scpe_sig_raw.plot("displayname",['Scope Signal (Run ID: ',num2str(dataTable_.run_id)],"fignum",dataTable_.run_id,"clear",0); - Scpe_sig_resampled = Scpe_sig_raw.resample("fs_in",Scpe_sig_raw.fs,"fs_out",fsym); - [~,Scpe_cell,found_sync,test,shifts] = Scpe_sig_resampled.tsynch("reference",Symbols,"fs_ref",dataTable_.symbolrate,"debug_plots",0); - - shifts_mus = shifts./Scpe_sig_resampled.fs .*1e6; - Scpe_sig_resampled = Scpe_sig_resampled.normalize("mode","rms"); - % Scpe_sig_resampled.plot("displayname",['Scope Signal (Run ID: ',num2str(dataTable_.run_id)],"fignum",dataTable_.run_id,"clear",0); - hold on; - % xline(shifts_mus,'HandleVisibility','off'); - - shifts = shifts-shifts(1)+1; - sep_sig = NaN(M,length(Scpe_sig_resampled)); - avg_sig = NaN(M,length(Scpe_sig_resampled)); - for j = 1:size(Scpe_cell,1) - - [sep_sig_,avg_sig_]=showLevelScatter(Scpe_cell{j}.resample("fs_out",Symbols.fs),Symbols,"fignum",400); - s = shifts(j); - sep_sig(:,s+1:s+length(sep_sig_)) = sep_sig_; - avg_sig(:,s+1:s+length(avg_sig_)) = avg_sig_; - - end - disp(num2str(dataTable_.interference_path_length)); + [sep_sig, avg_sig] = showLevelScatter(Scpe_sig_raw, Symbols, ... + "fsym", fsym, ... + "syncFs", 2*fsym, ... + "fignum", 400, ... + "normalize", true, ... + "debug_plots", false); var(sep_sig,0,2,'omitnan') - - xax_in_sec = ((1:length(avg_sig)) / fsym) * 1e6; - figure();hold on; - cols = cbrewer2('Paired',8); - for p = 1:size(avg_sig,1) - sc=scatter(xax_in_sec,sep_sig(p,:),1,'.','MarkerEdgeColor',cols((2*p)-1,:),'MarkerEdgeAlpha',0.2); - end - for p = 1:size(avg_sig,1) - sc=plot(xax_in_sec,avg_sig(p,:),'LineWidth',1,'Color',cols((2*p),:)); - end - - yline(unique(Symbols.signal),'HandleVisibility','off'); - % xline(shifts./ fsym .*1e6,'HandleVisibility','off'); - xlabel('time in $\mu$s'); - ylabel('Normalized Amplitude'); - xlim([0 25]); - ylim([-3 3]); - - drawnow; -end \ No newline at end of file +end diff --git a/projects/ECOC_2025_MPI/submit_dsp.m b/projects/ECOC_2025_MPI/submit_dsp.m index 4c78fd2..50211a8 100644 --- a/projects/ECOC_2025_MPI/submit_dsp.m +++ b/projects/ECOC_2025_MPI/submit_dsp.m @@ -12,7 +12,7 @@ arguments savePath options.parallel (1,1) logical = true options.max_occurences = 1; - options.paramstruct = struct(); + options.userParameters = struct(); end if options.parallel @@ -31,7 +31,7 @@ if options.parallel 'storage_path', savePath, ... 'append_to_db', 1, ... 'max_occurences', options.max_occurences, ... - 'parameters', options.paramstruct ... + 'userParameters', options.userParameters ... ); output = []; @@ -49,7 +49,7 @@ else 'storage_path', savePath, ... 'append_to_db', 1, ... 'max_occurences', options.max_occurences,... - 'parameters', options.paramstruct ... + 'userParameters', options.userParameters ... ); future = []; % No future since it's synchronous diff --git a/projects/MLSE_ML_based/experimental_data.m b/projects/MLSE_ML_based/experimental_data.m index a594a92..1c02517 100644 --- a/projects/MLSE_ML_based/experimental_data.m +++ b/projects/MLSE_ML_based/experimental_data.m @@ -15,7 +15,7 @@ duob_mode = db_mode(strrep(dataTable.db_mode,'"','')); % end % Load and Sync signal data from DB -[Tx_bits, Symbols, Scpe_cell, ~] = loadAndSyncSignalDataFromDb(dataTable, dsp_options); +[Tx_bits, Symbols, Scpe_cell, ~] = loadAndSyncRunSignals(dataTable, dsp_options); % Preprocess signal Scpe_sig = preprocessSignal(Scpe_cell{1}, Symbols, fsym); @@ -160,4 +160,4 @@ set(gca,'YScale','log'); xlabel('Epoch'); ylabel('BER'); title('PAM-4; L=3; SNR=20; AWGN Channel'); -plot(1:200,ml_mlse_equalizer_adap.ber,'DisplayName','Adaptive mu'); \ No newline at end of file +plot(1:200,ml_mlse_equalizer_adap.ber,'DisplayName','Adaptive mu'); diff --git a/projects/minimal_dsp_recipe_example/dsp_architecture_overview.svg b/projects/minimal_dsp_recipe_example/dsp_architecture_overview.svg new file mode 100644 index 0000000..090f84a --- /dev/null +++ b/projects/minimal_dsp_recipe_example/dsp_architecture_overview.svg @@ -0,0 +1,188 @@ + + DSP Job Architecture Overview + Architecture overview for the MATLAB DSP job framework: entry scripts define options and user parameters, submitJobs creates jobs, runBatch executes them, dsp_runid loads inputs and calls a recipe, and recipes use existing EQ kernels. + + + + + + + + + + + + + + DSP Job Architecture Overview + Separation of entry scripts, job execution, input loading, recipe orchestration, existing EQ kernels, and result sinks. + + + ENTRY / USER INTENT + + + JOB CONSTRUCTION AND EXECUTION + + + WORKER INPUT AND RECIPE CALL + + + DSP RECIPE AND EXISTING ALGORITHM KERNELS + + + RESULTS + + + + Entry script + run_minimal_recipe.m or project script + + + dsp_options + mode, recipe, DB settings + append_to_db, debug_plots + + + userParameters + User-defined sweep and overrides + Example: dd_mode = [0, 1] + + + run_ids + Measured or simulated run records + Usually selected by QueryFilter + + + + DataStorage wh + Expands userParameters grid + Defines optional result containers + + + submitJobs + Builds one job per run_id and sweep index + Injects concrete userParameters + Stores matching output fields into wh + + + runBatch + Generic executor + Serial: workerFcn(args{:}) + Parallel: parfeval + fetchNext + + + + dsp_runid + Worker for one job + Normalizes input source + Iterates synchronized occurrences + + + Input source + run_id: loadDspInputFromRunId + file_paths: loadDspInputFromFilePaths + loadAndSyncRunSignals + Returns canonical dspInput struct + + + Canonical dspInput + Scpe_cell, Symbols, Tx_bits + fsym, M, duob_mode, dataTable + One occurrence passed to recipe + + + + Recipe function + dsp_recipe_minimal or project recipe + Orchestrates preprocessing and kernels + + + preprocessSignal + Resample or matched filter + Timing synchronization + Optional filtering / debug plots + + + Existing EQ kernels + ffe(eq, M, rx, symbols, bits, ...) + ml_mlse(eq, M, rx, symbols, bits, ...) + vnle_postfilter_mlse(eq, pf, mlse, M, rx, symbols, bits, ...) + duobinary_target / duobinary_signaling + Return result structs with .metrics and .config + + + + results cell array + + + Warehouse wh.sto fields + + + Optional DB writes + + + + + + + + + + + + + + + + + + + + + + + + + + + + append_to_db = true + + + + Naming rule + userParameters: + only user sweep / + override values + + + Key boundary + runBatch knows jobs, + not DSP. + Recipes know DSP. + + + Recipe output fields such as ffe_package or mlse_package are copied generically by dsp_runid and stored by submitJobs when wh contains matching storage. + diff --git a/projects/minimal_dsp_recipe_example/run_minimal_recipe.m b/projects/minimal_dsp_recipe_example/run_minimal_recipe.m new file mode 100644 index 0000000..c3b91fb --- /dev/null +++ b/projects/minimal_dsp_recipe_example/run_minimal_recipe.m @@ -0,0 +1,130 @@ +% Minimal example for exchangeable DSP recipes. +% +% This script demonstrates the framework wiring only. Set run_ids to one or +% more known runs before executing it. +% +% Signal and parameter flow: +% +% 1) *This script* here entry script defines run selection and DSP options. +% - run_ids selects the measured/simulated records to process. +% - dsp_options contains worker configuration, database connection values, +% the recipe function handle, and the parameter sweep definition. +% +% 1.1) *DataStorage* expands dsp_options.userParameters into concrete jobs. +% Example: +% dsp_options.userParameters.dd_mode = [0, 1] +% creates two jobs for each run_id: +% job 1: userParameters.dd_mode = 0 +% job 2: userParameters.dd_mode = 1 +% +% 2) *submitJobs* builds the job list. +% For each run_id and each Warehouse index it forwards: +% dsp_runid(run_id, "recipe", @dsp_recipe_minimal, ... +% "userParameters", concreteParameterStruct, ...) +% +% 3) *runBatch* executes the job list. +% submitJobs delegates execution to runBatch. runBatch is the generic +% serial/parallel executor; it handles feval/parfeval, progress updates, +% result callbacks, and error callbacks, but it has no DSP, DB, Warehouse, +% or recipe-specific semantics. +% +% 4) *dsp_runid* is the worker. +% It does not own DSP algorithm choices. It only: +% - opens its own DB connection if mode/load/write requires it +% - loads metadata for the run_id +% - loads Tx_bits, Symbols, and received scope data +% - synchronizes/limits Scpe_cell occurrences +% - calls the selected recipe for each occurrence +% +% 5) The *recipe* is the project-specific DSP layer. +% dsp_recipe_minimal receives exactly: +% Scpe_sig_raw: one received scope/synchronized signal occurrence +% Symbols: transmitted symbols/reference +% Tx_bits: transmitted bit reference +% fsym, M, duob_mode: metadata from the run +% userParameters: one concrete parameter struct from the sweep +% +% 6) *Results* flow back as a standard output struct. +% submitJobs auto-creates wh.sto fields from non-empty recipe output +% fields. In this example, output.ffe_package creates wh.sto.ffe_package. + +clearvars; + +% === DSP settings === +% The recipe handle is the only DSP-specific dispatch point in this script. +% Switch it to another function with the same signature to run another DSP +% chain without changing submitJobs, runBatch, DB loading, or Warehouse logic. +dsp_options = struct(); +dsp_options.mode = "run_id"; +dsp_options.recipe = @dsp_recipe_minimal; +dsp_options.append_to_db = false; +dsp_options.max_occurences = 1; +dsp_options.debug_plots = true; + +dsp_options.database_type = "mysql"; +mpi = 1; +if mpi + dsp_options.dataBase = "labor"; + dsp_options.storage_path = "W:\labdata\ECOC Silas\ecoc_2025"; +else + dsp_options.dataBase = "labor_highspeed"; + dsp_options.storage_path = "W:\labdata\sioe_labor\"; +end +dsp_options.server = "192.168.178.192"; +dsp_options.port = 3306; +dsp_options.user = "silas"; +dsp_options.password = "silas"; + +db = DBHandler("dataBase", [dsp_options.dataBase],... + "type", dsp_options.database_type,... + "server", dsp_options.server,... + "user", dsp_options.user, "password", dsp_options.password); + +% Parameters are expanded by DataStorage and interpreted by the recipe. +% submitJobs treats this struct as transport data only; semantic meaning +% belongs to dsp_recipe_minimal or any replacement recipe. +dsp_options.userParameters = struct(); +dsp_options.userParameters.dd_mode = [0,1]; + +fp = QueryFilter(); +if mpi + fp.where('Runs', 'symbolrate', 'EQUALS', 112e9); + fp.where('Runs', 'fiber_length', 'EQUALS', 0); + fp.where('Runs', 'interference_path_length', 'EQUALS', 1000); + fp.where('Runs', 'sir', 'EQUALS', 52); + fp.where('Runs', 'db_mode', 'EQUALS', '"no_db"'); + fp.where('Runs', 'is_mpi', 'EQUALS', 1); + fp.where('Runs', 'pam_level', 'EQUALS', 4); + fp.where('Runs', 'wavelength', 'EQUALS', 1310); + fp.where('Runs', 'v_bias', 'EQUALS', 2.65); +else + % fp.where('Runs', 'pam_level','EQUALS', 4); + % fp.where('Runs', 'bitrate','EQUALS', 360e9); + % fp.where('Runs', 'fiber_length','EQUALS', 2); + % fp.where('Runs', 'is_mpi','EQUALS', 0); + % fp.where('Runs', 'wavelength','EQUALS', 1310); + % fp.where('Runs', 'db_mode','EQUALS', 0); + % fp.where('Runs', 'rop_attenuation','EQUAL', 0); +end + + +[dataTable,~] = db.queryDB(fp, db.getTableFieldNames('Runs')); + +run_ids = dataTable.run_id(1); + +if isempty(run_ids) + error("run_minimal_recipe:MissingRunIds", ... + "Set run_ids to one or more known run IDs before running this example."); +end + +%% === Warehouse setup === +% DataStorage defines the parameter grid. Result containers are created by +% submitJobs from returned recipe output fields. +wh = DataStorage(dsp_options.userParameters); + +%% === Run === +% submitJobs returns the raw per-job results and the filled Warehouse. For a +% single run_id and dd_mode = [0, 1], results is a 2-by-1 cell array. +[results, wh] = submitJobs(run_ids, dsp_options, processingMode.serial, ... + "wh", wh, ... + "waitbar", true); diff --git a/projects/minimal_dsp_recipe_example/run_minimal_runbatch_sweep.m b/projects/minimal_dsp_recipe_example/run_minimal_runbatch_sweep.m new file mode 100644 index 0000000..b5e62d8 --- /dev/null +++ b/projects/minimal_dsp_recipe_example/run_minimal_runbatch_sweep.m @@ -0,0 +1,64 @@ +% Minimal runBatch example without DSP, DB, or Warehouse. +% +% This script shows the generic execution layer directly: +% +% jobs -> runBatch -> workerFcn -> results +% +% runBatch only needs a function handle and a struct array of jobs. Each +% jobs(i).args cell is expanded into the worker function call. + +clearvars; + +% === Define a small calculation sweep === +gainValues = [0.5, 1, 2]; +offsetValues = [-1, 0, 1]; + +jobs = repmat(struct('args', {{}}, 'label', "", 'meta', struct()), ... + 1, numel(gainValues)*numel(offsetValues)); + +jobIndex = 0; +for gi = 1:numel(gainValues) + for oi = 1:numel(offsetValues) + jobIndex = jobIndex + 1; + + gain = gainValues(gi); + offset = offsetValues(oi); + + jobs(jobIndex).args = {gain, offset}; + jobs(jobIndex).label = sprintf("gain %.1f, offset %.1f", gain, offset); + jobs(jobIndex).meta.gainIndex = gi; + jobs(jobIndex).meta.offsetIndex = oi; + jobs(jobIndex).meta.gain = gain; + jobs(jobIndex).meta.offset = offset; + end +end + +% === Execute === +% Switch mode to processingMode.parallel to run the same jobs via parfeval. +results = runBatch(@minimalCalculationWorker, jobs, ... + "mode", processingMode.serial, ... + "waitbar", false, ... + "resultHandler", @printResult); + +% Convert the cell result list into a table for convenient inspection. +resultTable = struct2table([results{:}].'); + +disp(resultTable); + +function result = minimalCalculationWorker(gain, offset) +%MINIMALCALCULATIONWORKER Tiny worker used by this example. +inputSignal = 1:5; +outputSignal = gain .* inputSignal + offset; + +result = struct(); +result.gain = gain; +result.offset = offset; +result.meanOutput = mean(outputSignal); +result.maxOutput = max(outputSignal); +result.outputSignal = outputSignal; +end + +function printResult(result, job, jobIndex) +fprintf('[%s | #%d] mean %.2f, max %.2f\n', ... + job.label, jobIndex, result.meanOutput, result.maxOutput); +end diff --git a/projects/minimal_dsp_recipe_example/workerError.mat b/projects/minimal_dsp_recipe_example/workerError.mat new file mode 100644 index 0000000..7f29705 Binary files /dev/null and b/projects/minimal_dsp_recipe_example/workerError.mat differ diff --git a/workerError.mat b/workerError.mat new file mode 100644 index 0000000..a9bdfa3 Binary files /dev/null and b/workerError.mat differ