diff --git a/Classes/00_signals/Signal.m b/Classes/00_signals/Signal.m index 5905556..1377ddc 100644 --- a/Classes/00_signals/Signal.m +++ b/Classes/00_signals/Signal.m @@ -971,8 +971,8 @@ classdef Signal maxA = max(sig(100:end-100))*1.3; minA = min(sig(100:end-100))*1.3; - % maxA = 0.12; - % minA = -0.08; + maxA = 0.12; + minA = -0.08; difference= maxA-minA; @@ -993,7 +993,9 @@ classdef Signal % beautify cm=flip(cbrewer2("RdYlBu",4096)); % cm=flip(cbrewer2("RdBu",4096)); + % cm=flip(cbrewer2("Blues",4096)); cm(1,:) = [1,1,1]; % set zeros to white => clean background + colormap(cm); % colormap('turbo'); % ax.CLim = [0 50]; diff --git a/Classes/04_DSP/Equalizer/ML_MLSE.m b/Classes/04_DSP/Equalizer/ML_MLSE.m index 3f584c3..168b837 100644 --- a/Classes/04_DSP/Equalizer/ML_MLSE.m +++ b/Classes/04_DSP/Equalizer/ML_MLSE.m @@ -1,3 +1,499 @@ +% classdef ML_MLSE < handle +% % ALGORITHM DESCRIBED IN: +% % W. Lanneer and Y. Lefevre, “Machine Learning-Based Pre-Equalizers for +% % Maximum Likelihood Sequence Estimation in High-Speed PONs,” +% % in 2023 31st European Signal Processing Conference +% +% % Further ML Refs: +% % https://machinelearningmastery.com/cross-entropy-for-machine-learning/ +% % https://docs.pytorch.org/docs/stable/generated/torch.nn.CrossEntropyLoss.html +% +% % The central idea is to overcome the (white-) noise assumption within the previously described +% % Viterbi algorithm, more precisely a closed-loop optimization is proposed that finds a suitable +% % filter-set to directly compute the branch metrics c_k (s,s^' ). These can directly be used to +% % carry out the conventional Viterbi algorithm. The system consists of S^L S=F linear FIR filters, +% % combined with one bias coefficient respectively. These filters take the received input samples to +% % compute the branch metrics estimates (c_k ) ̂(s,s^' ) according toThe central idea is to overcome +% % the (white-) noise assumption within the previously described Viterbi algorithm, more precisely +% % a closed-loop optimization is proposed that finds a suitable filter-set to directly compute the +% % branch metrics c_k (s,s^' ). These can directly be used to carry out the conventional Viterbi +% % algorithm. The system consists of S^L S=F linear FIR filters, combined with one bias coefficient +% % respectively. These filters take the received input samples to compute the branch metrics +% % estimates. Finally, the usual Viterbi is carried out... +% +% % Recommended Settings and some findings: +% +% % Requires many training epochs. According to ML people, 100,200 or +% % even up to 1000 epochs are normal for ML-convergence +% +% % The mu parameter _can_ be adaptive - using the cross entropy and when +% % analyzing the isolated training it looks very promisig. However, is +% % later use I found this is not as stable as a fixed learning rate. +% % mu = 0.1 worked good for me +% +% % Longer orders/ filter length are not always better. For me order=11 +% % was good. +% +% % Delay factor (delta) is good when the order is also increased. With +% % order = 11, a delta of =4 shows good results +% +% properties +% sps % usually 2 +% order +% e +% e_tr +% error +% +% len_tr +% mu_tr +% epochs_tr +% +% dd_mode % 1 or 0 to set DD-mode on or off +% mu_dd %weight update in dd mode +% epochs_dd +% +% adaptive_mu +% +% constellation +% +% L %viterbi memory length +% +% alpha +% DIR +% DIR_flip +% trellis_states +% +% traceback_depth +% +% % --- Added internal class variables used later --- +% S +% Nf +% delta +% nStates +% nFeasible +% combs +% first_sym +% last_sym +% valid +% valid_to_idx +% valid_from_idx +% w +% +% % --- New: fast state lookup --- +% true_to_state_idx +% state_dict % containers.Map: key(sequence)->state index +% key_fmt = '%.8g_'; % key format for sequence strings +% nSym % |constellation| +% +% ber = [] +% ce = ones(1,1); +% end +% +% methods +% function obj = ML_MLSE(options) +% arguments(Input) +% +% options.sps = 2; +% options.order = 15; +% +% options.len_tr = 4096; +% options.mu_tr = 0; +% options.epochs_tr = 5; +% +% options.dd_mode = 1; +% options.mu_dd = 1e-5; +% options.epochs_dd = 5; +% +% options.adaptive_mu = 1; +% +% options.delta = 0; +% options.traceback_depth = 1024; +% +% options.L = 1 +% +% end +% +% fn = fieldnames(options); +% for n = 1:numel(fn) +% obj.(fn{n}) = options.(fn{n}); +% end +% +% obj.e = zeros(obj.order,1); +% obj.error = 0; +% end +% +% function [X,X_viterbi] = process(obj, X, D) +% +% % actual processing of the signal (steps 1. - 3.) +% % 1 normalize RMS +% X = X.normalize("mode","rms"); +% +% % Use sorted constellation for deterministic mapping +% obj.constellation = sort(unique(D.signal),'ascend'); +% obj.nSym = numel(obj.constellation); +% +% if length(X)/length(D) ~= obj.sps +% warning('Signal length does not fit to reference!'); +% end +% +% % ============================================================== +% % INITIALIZATION (only before final epoch and detection mode) +% % ============================================================== +% +% % --- Parameters +% obj.S = numel(obj.constellation); % alphabet size +% obj.Nf = obj.order*obj.sps; % filter length +% % obj.delta = 3;%ceil(obj.Nf/2); % delay parameter +% obj.nStates = obj.S^obj.L; +% obj.nFeasible = obj.nStates*obj.S; +% +% % --- Trellis mapping +% obj.trellis_states = reshape(obj.constellation,1,[]); +% pre_comb_mat = repmat(obj.trellis_states, obj.L, 1); +% pre_comb_cell = mat2cell(pre_comb_mat, ones(1,obj.L), size(pre_comb_mat,2)); +% obj.combs = fliplr(combvec(pre_comb_cell{:}).'); % rows: states, columns: [x_k, x_{k-1}, ...] +% obj.first_sym = obj.combs(:,1); +% obj.last_sym = obj.combs(:,end); +% obj.nStates = size(obj.combs,1); +% +% % --- Valid transitions +% obj.valid = false(obj.nStates); +% for from = 1:obj.nStates +% for to = 1:obj.nStates +% if all(obj.combs(to,2:end) == obj.combs(from,1:end-1)) +% obj.valid(to,from) = true; +% end +% end +% end +% [obj.valid_to_idx, obj.valid_from_idx] = find(obj.valid); +% +% % --- Allocate vectors and weights +% % !! IF SHAPE FIT, then we already have smth there an we want +% % to start with the existing fitler-set +% if isempty(obj.w) || any(size(obj.w) ~= [obj.Nf+1,obj.nFeasible]) +% obj.w = zeros(obj.Nf+1,obj.nFeasible); % filter weights per transition + bias tap +% obj.w = randn(obj.Nf+1,obj.nFeasible); +% end +% +% % --- Precompute dictionary for fast state lookup (sequence -> state) +% keys = cell(obj.nStates,1); +% for i = 1:obj.nStates +% keys{i} = obj.seq_key(obj.combs(i,:)); % combs row is already [x_k, x_{k-1}, ...] +% end +% obj.state_dict = containers.Map(keys, 1:obj.nStates); +% +% % ============================================================== +% % TRAINING +% % ============================================================== +% +% % Training Mode +% n = obj.len_tr; +% training = 1; +% obj.equalize(X.signal, D.signal,obj.mu_tr,obj.epochs_tr,n,training); +% obj.e_tr = obj.e; +% +% % ============================================================== +% % DD-Mode / Fixed Mode +% % ============================================================== +% +% % Decision Directed Mode +% n = X.length; +% training = 0; +% [y,y_vit]=obj.equalize(X.signal, D.signal,obj.mu_dd,obj.epochs_dd,n,training); +% +% X_viterbi = X; +% +% X.signal = y; +% X.fs = D.fs; %change sampling frequency of outgoing signal from fdac e.g. 2 sps to symbol spaced = fsym +% lbdesc = [num2str(obj.order),' tap FFE']; +% X = X.logbookentry(lbdesc); % append to logbook +% +% X_viterbi.signal = y_vit; +% X_viterbi.fs = D.fs; %change sampling frequency of outgoing signal from fdac e.g. 2 sps to symbol spaced = fsym +% lbdesc = [num2str(obj.order),'order FFE + PF + Viterbi']; +% X_viterbi = X_viterbi.logbookentry(lbdesc); % append to logbook +% end +% +% function [y,y_ref] = equalize(obj,x,d,mu,epochs,N,training) +% % ============================================================== +% % FFE + Whitening + ML-Based Branch Metric Estimation + Viterbi +% % ============================================================== +% debug = 1; +% showPlots = 1; +% +% % --- Input padding and preallocation +% y = zeros(N,1); +% +% % number of symbol steps in this block +% nSymbols = ceil(N/obj.sps); +% +% for epoch = 1:epochs +% +% % state metrics (log-domain costs): keep as column [nStates×1] +% pm = zeros(obj.nStates,1); % v_{k-1}(s′) +% c_hat = zeros(1,obj.nFeasible); +% v_tilde = zeros(1,obj.nFeasible); +% pred = zeros(nSymbols, obj.nStates, 'uint32'); +% pm_sto = nan(obj.nStates, nSymbols,'like',pm); +% CE_accum = 0; +% +% +% %%% START IDX +% if training +% max_start = length(x) - ( (ceil(N/obj.sps)-1)*obj.sps + 1 ); +% max_start = max(1, max_start); % safety +% start_sample = randi([1, max_start], 1); %rnd training; not really good +% start_sample = 1; +% end_sample = start_sample + (ceil(N/obj.sps)-1)*obj.sps; +% else +% start_sample = 1;%obj.len_tr; +% end_sample = N; +% end +% +% start_symbol = 1 + floor((start_sample - 1)/obj.sps); % ABSOLUTE symbol index +% +% if numel(d) >= obj.L && start_symbol >= obj.L +% init_seq = d(start_symbol-obj.L+1 : start_symbol); % [d_k-L+1 ... d_k] +% true_to_state_idx = obj.state_dict(obj.seq_key(flip(init_seq))); % [d_k ... d_k-L+1] +% else +% % Not enough history – fall back to state 1 +% true_to_state_idx = uint32(1); +% end +% +% symbol = 0; +% for sample = start_sample:obj.sps:end_sample +% symbol = symbol + 1; +% k = symbol; +% sym_idx = start_symbol + (symbol - 1); +% +% % --- Build Δ-delayed observation window y_k +% i1 = sample - obj.Nf + 1 + obj.delta; +% i2 = sample + obj.delta; +% buf = x(max(1,i1):min(length(x),i2)); +% padL = max(0,1 - i1); +% padR = max(0,i2 - length(x)); +% yk = [zeros(padL,1); buf(:); zeros(padR,1)]; % Nf×1 +% yk = [yk;1]; +% +% % --- Predict branch metrics for all feasible transitions: c_hat +% c_hat = (yk.' * obj.w); % [1×nFeasible] +% c_hat = c_hat.'; % [nFeasible×1] +% +% % --- Extended path metrics: v_tilde = pm(from) + c_hat +% % normalize pm to avoid growth (invariant to additive const) +% pm = pm - min(pm); +% v_tilde = pm(obj.valid_from_idx) + c_hat; % [nFeasible×1] +% +% % ===== Gradient update (Algorithm 1) ===== +% +% if 1 %training +% % --- allocate storage once +% if epoch == 1 && symbol == 1 +% obj.true_to_state_idx = ones(ceil(N/obj.sps),1,'uint32'); +% end +% +% % --- previous "to" becomes current "from" +% if symbol > 1 +% true_from_state_idx = obj.true_to_state_idx(symbol-1); +% else +% true_from_state_idx = 1; +% end +% +% % --- compute or reuse "to" state +% if epoch == 1 +% % only compute in first epoch +% if sym_idx >= obj.L +% key_to = obj.seq_key(flip(d(sym_idx-obj.L+1 : sym_idx))); +% if isKey(obj.state_dict, key_to) +% obj.true_to_state_idx(symbol) = obj.state_dict(key_to); +% else +% obj.true_to_state_idx(symbol) = true_from_state_idx; +% end +% else +% obj.true_to_state_idx(symbol) = true_from_state_idx; +% end +% end +% +% % --- reuse cached state from second epoch onward +% true_to_state_idx = obj.true_to_state_idx(symbol); +% +% % --- ensure valid (from,to) +% dirac = zeros(obj.nFeasible,1); +% mask = obj.valid_from_idx==true_from_state_idx & ... +% obj.valid_to_idx ==true_to_state_idx; +% if any(mask) +% dirac(mask) = 1; +% else +% idx = find(obj.valid_from_idx==true_from_state_idx,1,'first'); +% dirac(idx) = 1; +% obj.true_to_state_idx(symbol) = obj.valid_to_idx(idx); +% end +% +% +% +% +% % softmax over -v_tilde (numerically safe shift) +% v_shift = -(v_tilde - min(v_tilde)); % shift to small positive numbers +% v_shift = min(v_shift, 100); % clamp exponent argument (≈ exp(50)=3e21) +% expv = exp(v_shift); +% p = expv ./ (sum(expv) + eps); +% +% % for logging only: +% CE_symbol(symbol) = -log(p(dirac==1) + eps); +% +% if sym_idx > obj.L +% CE_smooth(symbol) = 0.01*CE_symbol(symbol) + 0.99*CE_smooth(symbol-1); +% else +% if epoch > 1 +% CE_smooth(symbol) = obj.ce(end); %use ce from last epoch or =1 for very first round?! +% else +% CE_smooth(symbol) = CE_symbol(symbol); +% end +% end +% +% CE_accum = CE_symbol(symbol) + CE_accum; +% +% +% % gradient term (t - p) +% dmp = (dirac - p)'; % 1×nFeasible +% +% % Per-feature gradient; implicit expansion gives (Nf+1)×nFeasible +% dL_Dw = (yk) .* dmp; +% +% % Start updates only when the ABSOLUTE symbol index has ≥ L history +% if sym_idx >= obj.L +% if obj.adaptive_mu +% mu_eff = CE_smooth(sym_idx); +% mu_eff = max(min(mu_eff, 0.2), 1e-4); +% else +% mu_eff = mu; +% end +% +% obj.w = obj.w - mu_eff .* dL_Dw; % (Nf+1)×nFeasible +% end +% +% % if debug && epoch > 2 +% % figure(100); +% % subplot(4,1,1); +% % heatmap(p'); +% % title('Probs') +% % subplot(4,1,2); +% % heatmap(dmp); +% % title('Update') +% % subplot(4,1,3); +% % heatmap(dL_Dw); +% % title('Update') +% % subplot(4,1,4); +% % heatmap(bj.w); +% % title('Update') +% % +% % end +% +% end +% +% +% +% % --- Compare-Select (matrix form, min of costs) +% v_tilde_mat = inf(obj.nStates, obj.nStates); +% v_tilde_mat(obj.valid) = v_tilde; +% [pm_next, pred(k,:)] = min(v_tilde_mat, [], 2); +% +% % re-center to keep metrics bounded (decision-invariant) +% pm_next = pm_next - min(pm_next); +% +% pm = pm_next; +% pm_sto(:,symbol) = pm; +% end +% +% % --- Traceback (full; you can window with traceback_depth if desired) +% [~, s_end] = min(pm); +% viterbi_path = zeros(symbol,1,'uint32'); +% viterbi_path(symbol) = s_end; +% for n = symbol:-1:2 +% viterbi_path(n-1) = pred(n, viterbi_path(n)); +% end +% +% y_ref = d(start_symbol:end); +% y = obj.first_sym(viterbi_path); +% +% if debug && training +% sym_start = start_symbol; +% sym_end = start_symbol + symbol - 1; +% ref_slice = d(sym_start : sym_end); +% err = sum(y ~= ref_slice(1:numel(y))); +% +% try +% ref_bits = PAMmapper(obj.S,0).demap(ref_slice); +% eq_bits = PAMmapper(obj.S,0).demap(y); +% [~, ~, ber, ~] = calc_ber(ref_bits, eq_bits, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1); +% fprintf('Epoch: %d - BER: %.1e \n',epoch, ber); +% obj.ber(epoch) = ber; +% catch +% ser = err./length(y); +% fprintf('Epoch: %d - SER: %.1e \n',epoch, ser); +% end +% +% obj.ce(epoch) = CE_accum./symbol; +% +% if showPlots +% figure(10);clf +% subplot(3,2,1:2); +% heatmap(obj.w); +% title('Filter') +% +% subplot(3,2,3); +% v_tildemat = NaN(obj.nStates, obj.nStates); +% v_tildemat(obj.valid) = v_tilde; % log-domain scores +% heatmap(v_tildemat); +% title('Path Metrics (v_tilde)') +% +% subplot(3,2,4); +% scatter(1:symbol,pm_sto,1,'.') +% title('Path Metric Winners') +% +% subplot(3,2,5);hold on +% scatter(1:symbol,CE_symbol,1,'.'); +% scatter(1:symbol,CE_smooth,1,'.') +% title('Cross Entropy') +% +% subplot(3,2,6); hold on +% +% % Left y-axis: Cross Entropy (linear) +% yyaxis left +% scatter(1:length(obj.ce), obj.ce, 10, 's', 'filled') +% ylabel('Cross Entropy') +% +% % Right y-axis: BER (logarithmic) +% yyaxis right +% scatter(1:length(obj.ber), obj.ber, 10, 'd', 'filled') +% set(gca, 'YScale', 'log') +% ylabel('BER (log scale)') +% +% xlim([1, epochs]) +% xlabel('Epoch') +% title('Cross Entropy // BER') +% grid on +% +% drawnow +% end +% end +% end +% end +% end +% +% methods (Access=private) +% function k = seq_key(obj, seq) +% % Build a stable key string for a sequence row vector in the *same order as combs rows* ([x_k, x_{k-1}, ...]) +% % Use rounding via sprintf to avoid floating-point issues. +% % seq must be a row vector. +% k = sprintf(obj.key_fmt, seq); +% end +% end +% end + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + classdef ML_MLSE < handle % --------------------------------------------------------------------- % W. Lanneer and Y. Lefevre, @@ -101,7 +597,7 @@ classdef ML_MLSE < handle obj.S = obj.nSym; obj.Nf = obj.order * obj.sps; obj.nStates = obj.S^obj.L; - obj.nFeasible = obj.nStates * obj.S; + obj.nFeasible = obj.nStates * obj.S; %feasible state transitions % --- Trellis mapping obj.trellis_states = reshape(obj.constellation,1,[]); @@ -164,8 +660,8 @@ classdef ML_MLSE < handle % EQUALIZE % ============================================================== function [y,y_ref] = equalize(obj,x,d,mu,epochs,N,training) - debug = 0; - showPlots = 0; + debug = 1; + showPlots = 1; y = zeros(N,1); nSymbols = ceil(N/obj.sps); @@ -241,6 +737,19 @@ classdef ML_MLSE < handle dirac(trans_idx)=1; end + % --- ensure valid (from,to) + if ~any(dirac) + mask = obj.valid_from_idx==true_from_state_idx & ... + obj.valid_to_idx ==true_to_state_idx; + if any(mask) + dirac(mask) = 1; + else + idx = find(obj.valid_from_idx==true_from_state_idx,1,'first'); + dirac(idx) = 1; + obj.true_to_state_idx(symbol) = obj.valid_to_idx(idx); + end + end + % =================================================================== % TRAINING MODE (weight update) % =================================================================== diff --git a/Classes/04_DSP/Timing_Recovery.m b/Classes/04_DSP/Timing_Recovery.m index ba16f10..23cc331 100644 --- a/Classes/04_DSP/Timing_Recovery.m +++ b/Classes/04_DSP/Timing_Recovery.m @@ -9,7 +9,7 @@ classdef Timing_Recovery < handle end methods(Access=public) - function obj = FFE(options) + function obj = Timing_Recovery(options) arguments(Input) options.timing_error_detector = 'Gardner'; @@ -25,12 +25,11 @@ classdef Timing_Recovery < handle obj.(fn{n}) = options.(fn{n}); end - obj.e = zeros(obj.order,1); - obj.error = 0; + end - function data_out = process(obj, data_in) + function [data_out,timing_error] = process(obj, data_in) timing_synchronization = comm.SymbolSynchronizer( ... "TimingErrorDetector", obj.timing_error_detector, ... @@ -39,7 +38,8 @@ classdef Timing_Recovery < handle "NormalizedLoopBandwidth", obj.normalized_loop_bandwidth, ... "DetectorGain", obj.detector_gain); - data_out.signal = timing_synchronization(data_in.signal); + data_out = data_in; + [data_out.signal,timing_error] = timing_synchronization(data_in.signal); end end diff --git a/Functions/EQ_structures/dsp_runid.m b/Functions/EQ_structures/dsp_runid.m index 4d207bd..c10cb5c 100644 --- a/Functions/EQ_structures/dsp_runid.m +++ b/Functions/EQ_structures/dsp_runid.m @@ -95,7 +95,7 @@ try adaption= 1; use_dd_mode = 1; - use_ffe = 0; + use_ffe = 1; use_dfe = 0; use_vnle_mlse = 0; use_dbtgt = 0; @@ -149,6 +149,9 @@ try % Scpe_sig.spectrum("fignum",22233,"normalizeTo0dB",0,"displayname",'Rx'); % Scpe_sig.eye(fsym,M,"fignum",1024); + Scpe_sig.signal = Scpe_sig.signal(1:2*Symbols.length); + Scpe_sig.signal = real(Scpe_sig.signal); + if duob_mode ~= db_mode.db_encoded if use_ffe @@ -246,7 +249,7 @@ try %ML-based MLSE (L=2) mu_ml = 0.01; training_epochs = 100; ml_mlse_equalizer = ML_MLSE("epochs_tr",training_epochs,"epochs_dd",1, ... - "len_tr",length(Scpe_sig),"mu_dd",mu_ml,"mu_tr",mu_ml,"order",11,"sps",2, ... + "len_tr",length(Scpe_sig)/4,"mu_dd",mu_ml,"mu_tr",mu_ml,"order",11,"sps",2, ... "traceback_depth",128,"L",1,"delta",4,"adaptive_mu",0); [ml_mlse_results] = ml_mlse(ml_mlse_equalizer, M, Scpe_sig, Symbols, Tx_bits,"precode_mode",duob_mode); diff --git a/Functions/EQ_structures/ffe.m b/Functions/EQ_structures/ffe.m index df4138e..c7e7c52 100644 --- a/Functions/EQ_structures/ffe.m +++ b/Functions/EQ_structures/ffe.m @@ -157,7 +157,7 @@ function displayAnalysis(eq_noise, eq_signal_sd, rx_signal, eq_, tx_symbols, M, % Initialize figure handles % Corrected line - added tx_symbols as second positional argument - showLevelScatter(rx_signal.resample("fs_out", tx_symbols.fs), tx_symbols, "fignum", 100); + % showLevelScatter(rx_signal.resample("fs_out", tx_symbols.fs), tx_symbols, "fignum", 100); warning off showLevelScatter(eq_signal_sd, tx_symbols, "fignum", 101); @@ -174,19 +174,20 @@ function displayAnalysis(eq_noise, eq_signal_sd, rx_signal, eq_, tx_symbols, M, end try - figure(339); - showEQfilter(eq_.e_tr, eq_signal_sd.fs.*2,"displayname",'training','fignum',339); - showEQfilter(eq_.e, eq_signal_sd.fs.*2,"displayname",'dec. directed','fignum',339); + figure(339);hold on + showEQfilter(eq_.e, eq_signal_sd.fs.*2,"displayname",'training','fignum',339); + % showEQfilter(eq_.e, eq_signal_sd.fs.*2,"displayname",'dec. directed','fignum',339); legend on end - try - figure(240); hold on; plot(pow2db(movmean(eq_.debug_struct.error_tr',100)));ylim([-30,3]);title('error training'); - - figure(241); hold on; plot(pow2db(movmean(eq_.debug_struct.update_tr',100)));title('update step training'); - - figure(242); hold on; plot(pow2db(movmean(eq_.debug_struct.update',1000)));title('update step dd'); - end - % eq_signal_sd.eye(eq_signal_sd.fs,M,"displayname",'Eye','fignum',105); + % try + % figure(240); hold on; plot(pow2db(movmean(eq_.debug_struct.error_tr',100)));ylim([-30,3]);title('error training'); + % + % figure(241); hold on; plot(pow2db(movmean(eq_.debug_struct.update_tr',100)));title('update step training'); + % + % figure(242); hold on; plot(pow2db(movmean(eq_.debug_struct.update',1000)));title('update step dd'); + % end + + eq_signal_sd.eye(eq_signal_sd.fs,M,"displayname",'Eye','fignum',105); end \ No newline at end of file diff --git a/Functions/EQ_structures/ml_mlse.m b/Functions/EQ_structures/ml_mlse.m index f2d468d..04fe317 100644 --- a/Functions/EQ_structures/ml_mlse.m +++ b/Functions/EQ_structures/ml_mlse.m @@ -46,6 +46,12 @@ end ml_mlse_results.config = Equalizerstruct(); eq_small = strip_eq(eq_, 10); +fn = fieldnames(eq_small); +for k = 1:numel(fn) + if issparse(eq_small.(fn{k})) + eq_small.(fn{k}) = full(eq_small.(fn{k})); + end +end json_str = jsonencode(eq_small); ml_mlse_results.config.eq = jsonencode(eq_); @@ -135,6 +141,9 @@ function eq_out = strip_eq(eq_, max_elems) eq_.(props{i}) = []; end end + if issparse(val) + eq_.(props{i}) = find(eq_.(props{i})); + end end eq_out = eq_; end diff --git a/Functions/EQ_structures/vnle_postfilter_mlse.m b/Functions/EQ_structures/vnle_postfilter_mlse.m index c7f10e1..0e6aea5 100644 --- a/Functions/EQ_structures/vnle_postfilter_mlse.m +++ b/Functions/EQ_structures/vnle_postfilter_mlse.m @@ -283,7 +283,7 @@ if ~isempty(postFFE) showEQcoefficients('n1', postFFE.e, "displayname", 'Coefficients', 'fignum', 338); end -showEQcoefficients('n1', eq_.e,'n2', eq_.e2,'n3', eq_.e3, "displayname", 'Coefficients', 'fignum', 339); +% showEQcoefficients('n1', eq_.e1,'n2', eq_.e2,'n3', eq_.e3, "displayname", 'Coefficients', 'fignum', 339); showEQfilter(eq_.e, eq_signal_sd.fs.*2); figure(340); clf; diff --git a/Functions/Theory/Dissertation/mach_zehnder_modulator.m b/Functions/Theory/Dissertation/mach_zehnder_modulator.m new file mode 100644 index 0000000..7b379c2 --- /dev/null +++ b/Functions/Theory/Dissertation/mach_zehnder_modulator.m @@ -0,0 +1,140 @@ +% Minimal MZM transfer-function demo (sinusoidal drive) — aligned with your notation +% +% Implements exactly: +% E_out(t) = E0 * exp(j*w0*t) * exp(-j*w0*L*n_eff/c0) * 1/2 * [ exp(-j*phi1(t)) + rho*exp(-j*phi2(t)) ] +% with phi_{1,2}(t) = pi * v_{1,2}(t)/Vpi +% +% Push-pull: +% v1(t) = +v_drive(t)/2 , v2(t) = -v_drive(t)/2 => phi1 = +pi/2 * v_drive/Vpi, phi2 = -pi/2 * v_drive/Vpi +% +% And the ideal TF (rho=1): +% E_out/E_in = exp(-j*w0*L*n_eff/c0) * cos( (pi/2) * v_drive/Vpi ) +% +% Note: E_in(t) = E0 * exp(j*w0*t) in this script. + +% Parameters +c0 = physconst('lightspeed'); % [m/s] +lambda0 = 1310e-9; % [m] +omega0 = 2*pi*c0/lambda0; + +L = 5e-3; % [m] effective phase section length (set as needed) +n_eff = 2.2; % [-] effective index (set as needed) + +E0 = 1; % field amplitude (arbitrary) +Vpi = 3.2; % [V] half-wave voltage (your V_pi) + +% Drive +f0 = 1e9; % [Hz] +fs = 100e9; % [Hz] +Nper = 1; % number of periods +Vpp = 0.5*Vpi; % [V] peak-to-peak of v_drive(t) + +biasV = 2; % [V] differential bias added to v_drive + +% Analytic +v_ = linspace(-1,2, 2001); +% Field transfer function (amplitude) +Field_mzm_analytic = cos((pi/2)*v_); +% Power transfer function (intensity) +P_mzm_analytic = Field_mzm_analytic.^2; + +% Imbalance factor in YOUR notation: +rho = 1; % rho=1 -> ideal balanced MZM (collapses to ideal TF) + +% Time axis + differential drive voltage v_drive(t) +T = Nper/f0; +t = (0:1/fs:T-1/fs).'; + +v_drive = biasV + (Vpp/2)*sin(2*pi*f0*t); % v_drive(t) (peak = Vpp/2) + +% Push-pull branch voltages (consistent with v_drive = v1 - v2) +v1 = +0.5*v_drive; % arm 1 +v2 = -0.5*v_drive; % arm 2 + +% Phases phi1, phi2 +phi1 = pi * v1 / Vpi; +phi2 = pi * v2 / Vpi; + +% Fields: E_in and E_out (exactly your Eq. (mzm_e_field)) +E_in = E0 .* exp(1i*omega0*t); + +common_phase = exp(-1i * (omega0*L*n_eff/c0)); % exp(-j*omega0*L*n_eff/c0) + +E_out = E0 .* exp(1i*omega0*t) .* common_phase .* 0.5 .* ... + ( exp(-1i*phi1) + rho .* exp(-1i*phi2) ); + +% Transfer function (numerical): E_out/E_in +H_num = E_out ./ E_in; + +% Power (normalized) +Pnorm_num = abs(H_num).^2; % since |E_out/E_in|^2 + + + +% Ideal TF (analytic) for comparison (rho=1, push-pull) +H_ideal = common_phase .* cos( (pi/2) * (v_drive./Vpi) ); + +Pnorm_ideal = abs(H_ideal).^2; +Pnorm_math = cos( (pi/2) * (v_drive./Vpi) ).^2; + + + + + +set(groot, 'defaultLegendInterpreter', 'tex'); +set(groot, 'defaultAxesTickLabelInterpreter', 'tex'); +set(groot, 'defaultTextInterpreter', 'tex'); + +% Normalized voltage axis (multiples of Vpi) +v_norm = v_drive./Vpi; + +colfield = [0,0,0]; %is black +colpow = linspecer(2); +colpow = colpow(1,:); +colvdrive = linspecer(2); +colvdrive = colvdrive(2,:); + +%% SIGNAL IN +figure(1); clf +plot(v_norm,t*1e9, 'LineWidth', 1.0,'Color',colvdrive); grid on; +ylabel('t [ns]'); xlabel('v_{drive}(t)/V_\pi'); +title('Drive voltage (normalized)'); +xlim([min(v_) max(v_)]); + +%% IN/OUT (static transfer) — normalized x-axis + analytic curve +figure(2); clf +plot(v_, Field_mzm_analytic, 'LineWidth', 1.2,'LineStyle','--','Color',colfield); hold on;% analytic power TF +plot(v_, P_mzm_analytic, 'LineWidth', 1.2, 'Color',colpow); hold on;% analytic power TF +% show input time signal +plot(v_norm,-1+t*1e9, 'LineWidth', 1.0,'Color',colvdrive); grid on; +% show output time signal +plot(2+t*1e9, Pnorm_num, 'LineWidth', 1.0,'DisplayName','Intensity', 'Color',colvdrive); hold on; +plot(2+t*1e9, real(H_ideal), '--', 'LineWidth', 1.0,'DisplayName','Field','Color',colfield); hold on; +scatter(v_norm, Pnorm_num, 12, '.', 'LineWidth', 1,'MarkerEdgeColor',colvdrive); +scatter(biasV./Vpi,(cos((pi/2)*biasV./Vpi)^2),10,'Marker','o'); +line([min(v_drive), min(v_drive)]./Vpi,[(cos((pi/2)*min(v_drive)./Vpi)^2), -2],'linewidth',0.5,'color','black','linestyle','--'); +line([max(v_drive) max(v_drive)]./Vpi,[(cos((pi/2)*max(v_drive)./Vpi)^2), -2],'linewidth',0.5,'color','black','linestyle','--'); +xline([min(v_norm) max(v_norm)]) + +grid on; +xlabel('v_{drive}(t)/V_\pi'); ylabel('|E_{out}/E_{in}|^2'); +% legend +xlim([min(v_) max(v_)+1]); +ylim([-1 1]); + +% mat2tikz_improved('C:\Users\Silas\Documents\6971e0b65b380ca6d71c837f\02_IMDD_System\tikz\mzm.tex'); + + +%% +% % FIELD TF (only field here; do not mix power into this figure) +figure(3); clf +% plot(t*1e9, real(H_num), 'LineWidth', 1.0); hold on; +% plot(t*1e9, real(H_ideal), '--', 'LineWidth', 1.0,'DisplayName','Field','Color',colfield); hold on; +plot(t*1e9, Pnorm_num, 'LineWidth', 1.0,'DisplayName','Intensity', 'Color',colpow); hold on; +grid on; +xlabel('t [ns]'); ylabel('Re\{E_{out}/E_{in}\}'); +legend +mat2tikz_improved('C:\Users\Silas\Documents\6971e0b65b380ca6d71c837f\02_IMDD_System\tikz\mzm_out.tex'); + + + diff --git a/Functions/mat2tikz_improved.m b/Functions/mat2tikz_improved.m new file mode 100644 index 0000000..a36b1d1 --- /dev/null +++ b/Functions/mat2tikz_improved.m @@ -0,0 +1,23 @@ +function mat2tikz_improved(filename) +arguments + % Default to the path in your example if no argument is provided + filename (1,1) string = 'C:\Users\Silas\Documents\Dissertation\00_Examples\tikz\textfig.tikz'; +end +cleanfigure; +matlab2tikz(char(filename), ... + 'width','\fwidth', ... + 'height','\fheight', ... + 'showInfo',false, ... + 'extraAxisOptions',{ ... + 'legend style={font=\footnotesize}', ... + 'xlabel style={font=\color{white!15!black},font=\small},',... + 'ylabel style={font=\color{white!15!black},font=\small},',... + 'legend columns=1', ... + 'every axis/.append style={font=\scriptsize}',... + 'legend columns=1',... + 'legend style={at={(0.02,0.98)},font=\footnotesize,draw=black!60,rounded corners=2pt,inner sep=1pt,fill=white,column sep=6pt,anchor= north west}',... + 'legend style={at={(0.02,0.98)},draw=white!0!white,font=\scriptsize,inner sep=0.1pt,fill=white,column sep=1pt,anchor= north west}',... + 'every axis/.append style={font=\scriptsize}',... + }); + +end \ No newline at end of file diff --git a/Libs/Violinplot-Matlab-master/violinplot.m b/Libs/Violinplot-Matlab-master/violinplot_community.m similarity index 100% rename from Libs/Violinplot-Matlab-master/violinplot.m rename to Libs/Violinplot-Matlab-master/violinplot_community.m diff --git a/Libs/wesanderson_colors/WesPalette.m b/Libs/wesanderson_colors/WesPalette.m index 7cd1dcb..13e976f 100644 --- a/Libs/wesanderson_colors/WesPalette.m +++ b/Libs/wesanderson_colors/WesPalette.m @@ -1,10 +1,13 @@ classdef WesPalette % WESPALETTE Wes Anderson color palettes with auto-completion % Usage: - % cmap = WesPalette.Zissou1.rgb() - % cmap = WesPalette.Zissou1.rgb(3) - - % https://github.com/karthik/wesanderson?tab=readme-ov-file + % cmap = WesPalette.Zissou1.rgb() % full palette + % cmap = WesPalette.Zissou1.rgb(3) % 3 colors (discrete default) + % cmap = WesPalette.Zissou1.rgb(12,"discrete") % any n, no interpolation + % cmap = WesPalette.Zissou1.rgb(256,"continuous") % smooth colormap (Lab interpolation) + % + % Requires: + % - colorspace.m (Pascal Getreuer) on MATLAB path for "continuous" mode enumeration BottleRocket1 @@ -34,21 +37,33 @@ classdef WesPalette end methods - function cmap = rgb(obj, n) + function cmap = rgb(obj, n, mode) % Return palette as Nx3 RGB colormap [0–1] + % + % n : number of requested colors (optional) + % mode : "discrete" (default) or "continuous" - hex = obj.hex(); + base_hex = obj.hex(); + base_rgb = WesPalette.hex2rgb(base_hex); - rgb = hex2rgb(hex); + if nargin < 2 || isempty(n) + cmap = base_rgb; + return; + end + if nargin < 3 || isempty(mode) + mode = "discrete"; + end + mode = lower(string(mode)); - if nargin == 2 - if n > size(rgb,1) - error('Requested %d colors, but only %d available.', ... - n, size(rgb,1)) - end - cmap = rgb(1:n,:); + validateattributes(n, {'numeric'}, {'scalar','integer','positive'}, mfilename, 'n'); + if mode ~= "discrete" && mode ~= "continuous" + error('mode must be "discrete" or "continuous".'); + end + + if mode == "discrete" + cmap = WesPalette.sample_discrete(base_rgb, n); else - cmap = rgb; + cmap = WesPalette.interpolate_continuous_lab(base_rgb, n); end end end @@ -56,78 +71,116 @@ classdef WesPalette methods (Access = private) function hex = hex(obj) % Internal HEX storage - switch obj case WesPalette.BottleRocket1 hex = {'#A42820','#5F5647','#9B110E','#3F5151','#4E2A1E','#550307','#0C1707'}; - case WesPalette.BottleRocket2 hex = {'#FAD510','#CB2314','#273046','#354823','#1E1E1E'}; - case {WesPalette.Rushmore1, WesPalette.Rushmore} hex = {'#E1BD6D','#EABE94','#0B775E','#35274A','#F2300F'}; - case WesPalette.Royal1 hex = {'#899DA4','#C93312','#FAEFD1','#DC863B'}; - case WesPalette.Royal2 hex = {'#9A8822','#F5CDB4','#F8AFA8','#FDDDA0','#74A089'}; - case WesPalette.Zissou1 hex = {'#3B9AB2','#78B7C5','#EBCC2A','#E1AF00','#F21A00'}; - case WesPalette.Zissou1Continuous hex = {'#3A9AB2','#6FB2C1','#91BAB6','#A5C2A3','#BDC881', ... '#DCCB4E','#E3B710','#E79805','#EC7A05','#EF5703','#F11B00'}; - case WesPalette.Darjeeling1 hex = {'#FF0000','#00A08A','#F2AD00','#F98400','#5BBCD6'}; - case WesPalette.Darjeeling2 hex = {'#ECCBAE','#046C9A','#D69C4E','#ABDDDE','#000000'}; - case WesPalette.Chevalier1 hex = {'#446455','#FDD262','#D3DDDC','#C7B19C'}; - case WesPalette.FantasticFox1 hex = {'#DD8D29','#E2D200','#46ACC8','#E58601','#B40F20'}; - case WesPalette.Moonrise1 hex = {'#F3DF6C','#CEAB07','#D5D5D3','#24281A'}; - case WesPalette.Moonrise2 hex = {'#798E87','#C27D38','#CCC591','#29211F'}; - case WesPalette.Moonrise3 hex = {'#85D4E3','#F4B5BD','#9C964A','#CDC08C','#FAD77B'}; - case WesPalette.Cavalcanti1 hex = {'#D8B70A','#02401B','#A2A475','#81A88D','#972D15'}; - case WesPalette.GrandBudapest1 hex = {'#F1BB7B','#FD6467','#5B1A18','#D67236'}; - case WesPalette.GrandBudapest2 hex = {'#E6A0C4','#C6CDF7','#D8A499','#7294D4'}; - case WesPalette.IsleofDogs1 hex = {'#9986A5','#79402E','#CCBA72','#0F0D0E','#D9D0D3','#8D8680'}; - case WesPalette.IsleofDogs2 hex = {'#EAD3BF','#AA9486','#B6854D','#39312F','#1C1718'}; - case WesPalette.FrenchDispatch hex = {'#90D4CC','#BD3027','#B0AFA2','#7FC0C6','#9D9C85'}; - case WesPalette.AsteroidCity1 hex = {'#0A9F9D','#CEB175','#E54E21','#6C8645','#C18748'}; - case WesPalette.AsteroidCity2 hex = {'#C52E19','#AC9765','#54D8B1','#B67C3B','#175149','#AF4E24'}; - case WesPalette.AsteroidCity3 hex = {'#FBA72A','#D3D4D8','#CB7A5C','#5785C1'}; end end end + + methods (Static, Access = private) + function rgb = hex2rgb(hex) + % hex: cellstr like {'#RRGGBB', ...} + if isstring(hex), hex = cellstr(hex); end + n = numel(hex); + rgb = zeros(n,3); + for i = 1:n + h = char(hex{i}); + if startsWith(h,'#'), h = h(2:end); end + if numel(h) ~= 6 + error('Invalid HEX color: %s', hex{i}); + end + rgb(i,1) = hex2dec(h(1:2))/255; + rgb(i,2) = hex2dec(h(3:4))/255; + rgb(i,3) = hex2dec(h(5:6))/255; + end + end + + function cmap = sample_discrete(base_rgb, n) + % No interpolation; allow any n by sampling/repeating. + k = size(base_rgb,1); + + if n <= k + idx = round(linspace(1, k, n)); % spread across palette + idx = max(1, min(k, idx)); + cmap = base_rgb(idx,:); + else + reps = floor(n / k); + rmd = mod(n, k); + cmap = [repmat(base_rgb, reps, 1); base_rgb(1:rmd,:)]; + end + end + + function cmap = interpolate_continuous_lab(base_rgb, n) + % Smooth interpolation in Lab using colorspace(). + % Requires colorspace.m by Pascal Getreuer on MATLAB path. + + k = size(base_rgb,1); + if k == 1 + cmap = repmat(base_rgb, n, 1); + return; + end + + % Convert to Lab, interpolate each channel, convert back + lab = colorspace('Lab<-RGB', base_rgb); + + t_base = linspace(0, 1, k); + t_new = linspace(0, 1, n); + + lab_new = zeros(n,3); + for c = 1:3 + lab_new(:,c) = interp1(t_base, lab(:,c), t_new, 'linear'); + end + + rgb_new = colorspace('RGB<-Lab', lab_new); + + % Clamp to displayable gamut + cmap = min(max(rgb_new, 0), 1); + end + end end diff --git a/Libs/wesanderson_colors/minimal_example_wespalette.m b/Libs/wesanderson_colors/minimal_example_wespalette.m index 9fcd99b..f308f27 100644 --- a/Libs/wesanderson_colors/minimal_example_wespalette.m +++ b/Libs/wesanderson_colors/minimal_example_wespalette.m @@ -5,8 +5,8 @@ y2 = 1e0 ./ (1 + exp(-0.4*(x-12))); % NLPN y3 = 1e-6 * 10.^(0.45*x); % RP on gamma y4 = 1e-2 * 10.^(0.18*(x-8)); % RP on beta2 -cmap = WesPalette.AsteroidCity1.rgb(4); -cmap = linspecer(4); +cmap = WesPalette.AsteroidCity1; +% cmap = linspecer(4); figure1=figure(202998);clf;hold on lw = 0.8; ms = 4; plot(x,y1,'LineWidth',lw,'Color',cmap(1,:),'Marker','o','MarkerEdgeColor',cmap(1,:),'MarkerFaceColor',[1,1,1],'MarkerSize',ms); diff --git a/projects/FSO_transmission/first_analysis.m b/projects/FSO_transmission/first_analysis.m index c47a2b1..9a327ab 100644 --- a/projects/FSO_transmission/first_analysis.m +++ b/projects/FSO_transmission/first_analysis.m @@ -1,15 +1,15 @@ -base = "C:\Users\Silas\Nextcloud\Dokumente\02_Ablage_Office\FSO_FP_QCL_60umUTC"; +base = "C:\Users\Silas\Nextcloud4\Dokumente\02_Ablage_Office\FSO_FP_QCL_60umUTC"; mode = 0; %0 oder 1 M = 2; all_files = dir(fullfile(base, "**/*.mat")); if M == 2 - tx_data = load("C:\Users\Silas\Nextcloud\Dokumente\02_Ablage_Office\FSO_FP_QCL_60umUTC\14G_PAM2\tx_info\tx_info_PAM2_14Gbd0.75RRC.mat"); + tx_data = load("C:\Users\Silas\Nextcloud4\Dokumente\02_Ablage_Office\FSO_FP_QCL_60umUTC\14G_PAM2\tx_info\tx_info_PAM2_14Gbd0.75RRC.mat"); filename = fullfile(base, "14G_PAM2\M=2_Rs=1.4e10_Fs=8e10_I=265mA_RoP=46.3mW_L=31m_PS=RRC_rolloff=0.75_Mode=Rise.mat"); elseif M == 4 - tx_data = load("C:\Users\Silas\Nextcloud\Dokumente\02_Ablage_Office\FSO_FP_QCL_60umUTC\6G_PAM4\tx_info\tx_info_PAM4_6Gbd0.6RRC.mat"); + tx_data = load("C:\Users\Silas\Nextcloud4\Dokumente\02_Ablage_Office\FSO_FP_QCL_60umUTC\6G_PAM4\tx_info\tx_info_PAM4_6Gbd0.6RRC.mat"); filename = fullfile(base, "6G_PAM4\M=4_Rs=6e9_Fs=8e10_I=255mA_RoP=42.3mW_L=31m_PS=RRC_rolloff=0.6_Mode=Rise.mat"); end @@ -53,10 +53,14 @@ Bits_ = PM.demap(Symbols); assert(ber == 0); %% For comparison, apply pulsef on Tx Symbols -Pform = Pulseformer("fsym",fsym,"fdac",fs,"pulse","rc","pulselength",16,"alpha",rolloff); -Digi_sig_compare = Pform.process(Symbols); -MF = Pulseformer("fsym",fsym,"fdac",2*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff); -Rx_sig_compare = MF.process(Digi_sig_compare); +Pform = Pulseformer("fsym",fsym,"fdac",fs,"pulse","rrc","pulselength",16,"alpha",rolloff); +Digi_sig_tx_compare = Pform.process(Symbols); + +Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff); +Digi_sync = Pform.process(Symbols); + +MF = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff); +Rx_sig_compare = MF.process(Digi_sig_tx_compare); %% @@ -68,53 +72,59 @@ scoperead_volts = (traceData.RawData - traceData.YReference) * traceData.YIncrem demystified = isequal(traceData.YData,scoperead_volts); assert(demystified); +timesig_compare = [0:1:datas.tr.lastData(1).trace.ch3.Points-1] ./ fs; +timesig = datas.tr.lastData(1).trace.ch3.XData; +assert(isequal(traceData.YData,scoperead_volts)); + Scope_sig = Electricalsignal(traceData.YData,"fs",fs); -Scope_sig.plot("displayname",'raw','fignum',100); -Scope_sig.spectrum("displayname",'raw','fignum',101) - +%% % 1) matched filter -% pulse is symmetric, hence we can use pulsef firectly as matched filter. +% pulse is symmetric, hence we can use pulsef directly as matched filter. % It feels off (bit I think correct) that the fsym is now the output freq.!! % -> output 2 sps to omit timing recovery!? -Pform = Pulseformer("fsym",fsym,"fdac",2*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff,"matched",1); -Rx_matched = Pform.process(Scope_sig); -Rx_matched.spectrum("displayname",'Signal after matched filter','fignum',1); +apply_matched_filter = 1; +k = 1; +if apply_matched_filter + + Pform = Pulseformer("fsym",fsym,"fdac",k*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff,"matched",1); + Rx_matched = Pform.process(Scope_sig); +else + + Rx_matched = Filter('filtdegree',4,"f_cutoff",fsym*0.5,"fs",Scope_sig.fs,"filterType",filtertypes.gaussian,"active",true).process(Scope_sig); + Rx_matched = Rx_matched.resample("fs_out",k*fsym); +end +Rx_matched.spectrum(); %% -sys = comm.SymbolSynchronizer('TimingErrorDetector', 'Gardner (non-data-aided)', ... - 'SamplesPerSymbol', 2, ... - 'DampingFactor', 0.7, ... - 'NormalizedLoopBandwidth', 0.01); -Rx_symbolsync = Rx_matched; -[Rx_symbolsync.signal, timing_error] = sys(Rx_matched.signal); +coefficients = arburg(Rx_matched.signal,25); -plot(timing_error); % If this is a ramp, you have drift! +figure() +[h,w] = freqz(1,coefficients,Rx_matched.length,"whole",Rx_matched.fs); +h = h/max(abs(h)); +hold on +w_ = (w - Rx_matched.fs/2); +plot(w_.*1e-9,20*log10(fftshift(abs(h))),'DisplayName',['Burg Coeffs: ', num2str(round(coefficients,2)), ' '],'LineWidth',2); -%% timing sync -> at this point we still have no symbol timing recovery, we -% % try to do this with 2sps EQ! - -[~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_symbolsync.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1); - - -%% not working.. -Rx_synced = Rx_synced_cell{1}; -len_tr = 4096*2; -mu_ffe1 = 0.0001; -mu_ffe2 = 0.0008; -mu_ffe3 = 0.001; -mu_dc = 0.005; -mu_ffe = [mu_ffe1 mu_ffe3 mu_ffe3]; -mu_dfe = 0.0004; -duob_mode = db_mode.no_db; - -Rx_synced.plot("displayname",'RX: Matched+Sync+2sps','fignum',2); - -Digi_sig_compare.normalize("mode","rms").spectrum("displayname",'Tx: RC-shaped','fignum',1,'normalizeTo0dB',0); -Rx_sig_compare.normalize("mode","rms").spectrum("displayname",'Tx: RC-shaped + matched filtered ','fignum',1,'normalizeTo0dB',0); -Rx_synced.normalize("mode","rms").spectrum("displayname",'RX: matched filtered + synced','fignum',1,'normalizeTo0dB',0); +%% Timing Rec +apply_timing_rec = 1; +if apply_timing_rec + [Rx_symbolsync, timing_error] = Timing_Recovery("timing_error_detector",'Gardner (non-data-aided)','sps',k,'damping_factor',0.1,'normalized_loop_bandwidth',0.1,'detector_gain',2.7).process(Rx_matched); + figure();plot(timing_error); + Rx_symbolsync.fs = fsym; + % Tsynch + [~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_symbolsync.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 0); + Rx_synced = Rx_synced_cell{1}; + sps = 1; +else + % Tsynch + [~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_matched.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 0); + Rx_synced = Rx_synced_cell{1}; + Rx_synced = Rx_synced.resample("fs_out",2*fsym); + sps = 2; +end if M == 2 ber_in_paper = 10^(-2.6); %fig 3a) 4 Gb/s MWIR FSO Transmission using Directly Modulated QCL and an Uncooled UTC-PD at Room-Temperature @@ -122,50 +132,121 @@ elseif M == 4 ber_in_paper = 10^(-2.5); end -%% -------------------- FFE -------------------- -% requires some more digging what is going on :-) -eq_ffe = EQ("Ne",[50, 5, 5],"Nb",[2,0,0], ... - "training_length",len_tr,"training_loops",5,"dd_loops",5, ... - "K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... + +%% + +Rx_synced = Rx_synced_cell{1}; + + +% -------------------- FFE -------------------- +% requires some more digging what is going on :-) +eq_ffe = EQ("Ne",[50, 1, 1],"Nb",[2,0,0], ... + "training_length",512,"training_loops",5,"dd_loops",5, ... + "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... "FFEmu",0,"plotfinal",0,"ideal_dfe",1); -ffe_results = ffe(eq_ffe,M,Rx_synced,Symbols,Bits, ... - "precode_mode",duob_mode,'showAnalysis',0,"postFFE",[], ... - "eth_style_symbol_mapping",mapping_style); +vars = logspace(-4,-3,36); -% ffe_results.metrics.print -fprintf('My EQ: %.1e \n',ffe_results.metrics.BER); +parfor i = 1:numel(vars) + + len_tr = 4096; + mu_ffe1 = 0.01;% mus(i);%0.0001; + mu_ffe2 = 0.0008; + mu_ffe3 = 0.001; + mu_dc = 0.005; + mu_ffe = [mu_ffe1 mu_ffe2 mu_ffe3]; + mu_dfe = vars(i); + duob_mode = db_mode.no_db; + + % requires some more digging what is going on :-) + eq_ffe_1 = EQ("Ne",[150, 1, 0],"Nb",[50,0,0], ... + "training_length",len_tr,"training_loops",5,"dd_loops",5, ... + "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... + "FFEmu",0,"plotfinal",0,"ideal_dfe",0); + + eq_ffe_2 = FFE("epochs_tr",1,"epochs_dd",vars(i),"len_tr",4096,"mu_dd",vars(i),"mu_tr",vars(i),"order",999,"sps",1,"decide",0, "adaption",adaption_method.nlms,"dd_mode",0); + % eq_ffe_2 = FFE_DFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"ffe_mu_dd",1e-5,"dfe_mu_dd",mus(i),"ffe_mu_tr",0,"dfe_mu_tr",0,"ffe_order",50,"dfe_order",10,"sps",1,"decide",1); + + + ffe_results = ffe(eq_ffe_1,M,Rx_synced,Symbols,Bits, ... + "precode_mode",duob_mode,'showAnalysis',1,"postFFE",[], ... + "eth_style_symbol_mapping",mapping_style); + + ffe_results.metrics.BER + bers(i) = ffe_results.metrics.BER; +end + +figure(); +plot(vars,bers); +yline(ber_in_paper) +beautifyBERplot(); +% fprintf('Paper: %.1e \n \n',ber_in_paper); +ffe_results.metrics.print("description",'FFE'); +fprintf('FFE: %.1e \n',ffe_results.metrics.BER); + %% -------------------- VNLE + MLSE -------------------- +len_tr = 4096; +mu_ffe1 = 0.0001;% mus(i);%0.0001; +mu_ffe2 = 0.0008; +mu_ffe3 = 0.001; +mu_dc = 0.005; +mu_ffe = [mu_ffe1 mu_ffe2 mu_ffe3]; +mu_dfe = 0.0004; +duob_mode = db_mode.no_db; -pf_ncoeffs = 1; -eq_v = EQ("Ne",[100, 5, 5],"Nb",[0, 0, 0], ... +pf_ncoeffs = 4; + +vars = 1:7; +bers = zeros(size(vars)); +parfor i = 1:numel(vars) + + eqv = EQ("Ne",[200, 1, 0],"Nb",[2, 0, 0], ... "training_length",len_tr,"training_loops",5,"dd_loops",5, ... - "K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... + "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... "FFEmu",0,"plotfinal",0,"ideal_dfe",1); -pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); -mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0,"eth_style",mapping_style).levels); - -[vnle_results, mlse_results] = vnle_postfilter_mlse(eq_v, pf_, mlse_, M, Rx_synced, Symbols, Bits, ... - "precode_mode", duob_mode, 'showAnalysis', 1, "postFFE", [], "eth_style_symbol_mapping", mapping_style); - -mlse_results.metrics.print("description",'MLSE') -fprintf('My EQ: %.1e \n',mlse_results.metrics.BER); -fprintf('Paper: %.1e \n \n',ber_in_paper); + pf_ncoeffs = vars(i); + pf = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); + mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0,"eth_style",mapping_style).levels); + [vnle_results, mlse_results] = vnle_postfilter_mlse(eqv, pf, mlse_, M, Rx_synced, Symbols, Bits, ... + "precode_mode", duob_mode, 'showAnalysis', 1, "postFFE", [], "eth_style_symbol_mapping", mapping_style); + fprintf('Paper: %.1e \n \n',ber_in_paper); + vnle_results.metrics.print("description",'VNLE'); + mlse_results.metrics.print("description",'MLSE'); + bers(i) = mlse_results.metrics.BER; +end +%% +figure();hold on +plot(vars,bers_ffe,'DisplayName','FFE [200,0,0] + PF + MLSE'); +plot(vars,bers_vnle,'DisplayName','VNLE [200,1,0] + PF + MLSE'); +plot(vars,bers_vnledfe,'DisplayName','VNLE [200,1,0] + DFE [2] + PF + MLSE'); +plot(vars,bers_vnledfe_ideal,'DisplayName','VNLE [200,1,0] + ideal DFE [2] + PF + MLSE'); +yline(ber_in_paper); +yline([2e-2, 4.85e-3, 3.8e-3, 2,2e-4],'LineWidth',2,'Color',[0.8,0.8,0.8],'LineStyle',':','HandleVisibility','off'); +ylim([1e-5,0.1]); +beautifyBERplot(); %% -------------------- DB target -------------------- mlse_db_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",M,'trellis_states',PAMmapper(M,0).levels); eq_ = EQ("Ne",[50, 5, 5],"Nb",[0,0,0],"training_length",len_tr,"training_loops",5,"dd_loops",5, ... - "K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1); + "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1); dbt_results = duobinary_target(eq_,mlse_db_, M, Rx_synced, Symbols, Bits, ... "precode_mode", duob_mode, 'showAnalysis', 0, "postFFE", [],"eth_style_symbol_mapping",mapping_style); -dbt_results.metrics.print("description",'Duobinary'); -mlse_results.metrics.print -fprintf('My EQ: %.1e \n',dbt_results.metrics.BER); fprintf('Paper: %.1e \n \n',ber_in_paper); +dbt_results.metrics.print("description",'Duobinary'); +%% + +%ML-based MLSE (L=2) +mu_ml = 0.01; training_epochs = 100; +ml_mlse_equalizer = ML_MLSE("epochs_tr",training_epochs,"epochs_dd",1, ... + "len_tr",len_tr,"mu_dd",mu_ml,"mu_tr",mu_ml,"order",11,"sps",1, ... + "traceback_depth",128,"L",3,"delta",4,"adaptive_mu",0); + +[ml_mlse_results] = ml_mlse(ml_mlse_equalizer, M, Rx_synced, Symbols, Bits,"precode_mode",duob_mode); +ml_mlse_results.metrics.print("description",'ML '); \ No newline at end of file diff --git a/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/Copy_of_FIGURE_WAVELENGTH.m b/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/Copy_of_FIGURE_WAVELENGTH.m index 21f881f..986fcdb 100644 --- a/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/Copy_of_FIGURE_WAVELENGTH.m +++ b/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/Copy_of_FIGURE_WAVELENGTH.m @@ -6,7 +6,7 @@ db = DBHandler("dataBase", "labor_highspeed", "type", database_type); pam_levels = [4, 6, 8]; % three tiles bitrate_set = 360e9; -fiberL = 10; +fiberL = 2; fields = [ db.getTableFieldNames('power_state_info'); @@ -96,7 +96,7 @@ end %% ============================================================ % PLOT — 1×3 (PAM-4, PAM-6, PAM-8) % ============================================================ -fig = figure(9110); clf; +fig = figure(9112); clf; tiledlayout(1,3,'TileSpacing','compact','Padding','compact'); lw = 1.8; @@ -223,19 +223,19 @@ ylabel(''); end -pos = 1e3.*[2.7770 1.2017 1.4000 0.3200]; -set(fig, 'Position', pos); +% pos = 1e3.*[2.7770 1.2017 1.4000 0.3200]; +% set(fig, 'Position', pos); %% === EXPORT === -outfile = 'C:\Users\Silas\Documents\latex\JLT_400G copy\media\matlab2tikz\wavelength_analysis.tikz'; -matlab2tikz(outfile, ... - 'width','\fwidth', ... - 'height','\fheight', ... - 'showInfo',false, ... - 'extraAxisOptions',{ ... - 'legend style={font=\footnotesize}', ... - 'legend columns=1' ... - }); +% outfile = 'C:\Users\Silas\Documents\latex\JLT_400G copy\media\matlab2tikz\wavelength_analysis.tikz'; +% matlab2tikz(outfile, ... +% 'width','\fwidth', ... +% 'height','\fheight', ... +% 'showInfo',false, ... +% 'extraAxisOptions',{ ... +% 'legend style={font=\footnotesize}', ... +% 'legend columns=1' ... +% }); diff --git a/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/FIGURE_EYES.m b/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/FIGURE_EYES.m index c22e90f..0d953af 100644 --- a/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/FIGURE_EYES.m +++ b/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/FIGURE_EYES.m @@ -2,7 +2,7 @@ dsp_options.storage_path = 'Z:\2024\sioe_labor\'; dsp_options.max_occurences = 1; database = DBHandler("dataBase", 'labor_highspeed', "type", 'mysql' ); -rate = [300e9]; + cols = cbrewer2('BuPu',25); cols = [cols(end-10:2:end,:)]; cols = cbrewer2('Set1',6); @@ -12,10 +12,10 @@ fig=figure(fignum);clf; dbmode = 0; - % 1 - PAM 4 with preemphasis fp = QueryFilter(); -M = 6; +M = 8; +rate = [360e9]; fp.where('Runs', 'pam_level','EQUALS', M); fp.where('Runs', 'bitrate','EQUALS', rate);%360,390 fp.where('Runs', 'fiber_length','EQUALS', 2); @@ -23,7 +23,7 @@ fp.where('Runs', 'wavelength','EQUALS', 1310); fp.where('Runs', 'db_mode','EQUALS', dbmode); fp.where('Runs', 'rop_attenuation','EQUAL', 0); -[dataTable,~] = db.queryDB(fp, database.getTableFieldNames('Runs')); +[dataTable,~] = database.queryDB(fp, database.getTableFieldNames('Runs')); dataTable = queryRunid(dataTable.run_id, database); fsym = dataTable.symbolrate; @@ -33,14 +33,36 @@ duob_mode = db_mode(strrep(dataTable.db_mode,'"','')); % Load and Sync signal data from DB [Tx_bits, Symbols, Scpe_cell, ~] = loadAndSyncSignalDataFromDb(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'); +%%%%%% SNR CHEAT - Avges the measured signal occurences found after correlation in "tsynch" %%%%%% +average_signals = 1; +if average_signals + Scpe_sig_avg = Scpe_sig_syncd; + scope_mean = zeros(size(Scpe_cell{1}.signal)); + for n=1:numel(Scpe_cell) + scope_mean = scope_mean + Scpe_cell{n}.signal; + end + scope_mean = scope_mean ./ n; + Scpe_sig_avg.signal = scope_mean; + + Scpe_sig_avg.spectrum("displayname","Scope PSD","fignum",20,"normalizeTo0dB",1); + Scpe_sig_avg.plot("displayname","Scope raw signal","fignum",27,"clear",1); + Scpe_sig_avg = Scpe_sig_avg.*1.25; + Scpe_sig_avg.eye(fsym,M,"fignum",rate.*1e-9*M,"displayname",' Eye of AVG Signal'); +end + % Preprocess signal -Scpe_sig = preprocessSignal(Scpe_cell{1}, Symbols, fsym); +Scpe_sig = preprocessSignal(Scpe_sig_avg, Symbols, fsym); Scpe_sig.eye(fsym,M,"fignum",M*10); + + %% === EXPORT TO TIKZ === -% outfile = ['C:\Users\Silas\Documents\latex\JLT_400G copy\media\matlab2tikz\eye_pam_',num2str(M),'.tikz']; -% outfile = ['C:\Users\Silas\Documents\latex\JLT_400G copy\media\matlab2tikz\vnle_optimization.tikz']; + +% outfile = ['C:\Users\Silas\Documents\latex\JLT_400G_submission\media\matlab2tikz\eye_pam_',num2str(M),'-2.tikz']; +% % outfile = ['C:\Users\Silas\Documents\latex\JLT_400G copy\media\matlab2tikz\vnle_optimization.tikz']; % matlab2tikz(outfile, ... % 'width','\fwidth', ... % 'height','\fheight', ... diff --git a/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/FIGURE_introduction.m b/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/FIGURE_introduction.m index 2213c8a..147c7d9 100644 --- a/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/FIGURE_introduction.m +++ b/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/FIGURE_introduction.m @@ -1,83 +1,139 @@ +tablename = 'C:\Users\Silas\Documents\latex\JLT_400G_submission\HighSpeedExperiments_oneandonly_csv.csv'; +% Returns a Table +data = readtable(tablename,"Delimiter",';','DecimalSeparator',','); + + %% ============================================================ % PLOT % ============================================================ +%% 1. DATA EXTRACTION & SETUP +%% 1. DATA EXTRACTION & SETUP +raw_M = data.M; +raw_baud = data.BaudRate; +raw_net = data.NetRate; +raw_codes = string(data.ZoteroCode); +raw_names = string(data.Name); +raw_band = string(data.Band); + +% Filter Valid Data +target_M = [2, 4, 6, 8]; +validIdx = ismember(raw_M, target_M) & ~isnan(raw_baud) & ~isnan(raw_net); + +Mvals = raw_M(validIdx); +baud = raw_baud(validIdx); +netrate = raw_net(validIdx); +codes = raw_codes(validIdx); +names = raw_names(validIdx); +bands = raw_band(validIdx); + +pam_list = target_M; +colors = flip(cbrewer2('SET1',4)); + +%% 2. PLOT (For Visual Check only) figure; hold on; -ms = 32; % scatter size -lw = 0.8; % line width - -for k = 1:4 % PAM-2/4/6/8 +ms = 20; +lw = 0.5; +for k = 1:length(pam_list) M = pam_list(k); idxPam = (Mvals == M); - - % Extract for this PAM + x = baud(idxPam); y = netrate(idxPam); + b = bands(idxPam); n = names(idxPam); - - % Get color for this PAM format col = colors(k,:); - - % ----- LEGEND FLAG (only add one entry per PAM) ----- + firstLegend = true; - - % ---- PLOT ALL POINTS (marker based on publication) ---- + for i = 1:sum(idxPam) - % marker selection by publication - pubIdx = find(pub_list == n(i), 1); - marker = markerlist{mod(pubIdx-1, nMarkers) + 1}; - + % Marker Logic + ms = 20; + if strcmpi(b(i), 'O') + marker = 'o'; + elseif strcmpi(b(i), 'C') + marker = 'd'; + else + marker = 's'; + end + if strcmpi(n(i), 'THIS WORK') + marker = 'pentagram'; + ms = 100; + end + + % Plot Scatter if firstLegend - h = scatter(x(i), y(i), ms, ... - 'Marker', marker, ... - 'MarkerEdgeColor', col, ... - 'MarkerFaceColor', col, ... + scatter(x(i), y(i), ms, 'Marker', marker, ... + 'MarkerEdgeColor', col, 'MarkerFaceColor', col, ... 'DisplayName', sprintf('PAM-%d', M)); firstLegend = false; else - h = scatter(x(i), y(i), ms, ... - 'Marker', marker, ... - 'MarkerEdgeColor', col, ... - 'MarkerFaceColor', col, ... + scatter(x(i), y(i), ms, 'Marker', marker, ... + 'MarkerEdgeColor', col, 'MarkerFaceColor', col, ... 'HandleVisibility','off'); end - - % ====== CUSTOM DATATIP CONTENT ====== - dt = h.DataTipTemplate; - dt.DataTipRows(1).Label = 'Baud rate'; - dt.DataTipRows(2).Label = 'Net rate'; - - % Add publication name - dt.DataTipRows(end+1) = dataTipTextRow('Publication', n(i)); - - end - - % ---- Fit (PAM-specific) ---- - valid = ~isnan(x) & ~isnan(y); - if sum(valid) >= 3 - p = polyfit(x(valid), y(valid), 2); - xfit = linspace(min(x(valid)), max(x(valid)), 200); - yfit = polyval(p, xfit); - - plot(xfit, yfit, ':', ... - 'LineWidth', lw, ... - 'Color', col, ... - 'HandleVisibility', 'off'); % do NOT add to legend + + % Fit lines + if length(x) >= 3 + [p, S, mu] = polyfit(x, y, 2); + xfit = linspace(min(x), max(x), 200); + yfit = polyval(p, xfit, S, mu); + plot(xfit, yfit, '-', 'LineWidth', lw, 'Color', col, 'HandleVisibility', 'off'); end end -grid on; +grid on; box on; xlabel('Baud rate [GBd]'); ylabel('Net rate [Gb/s]'); +% title('Check Command Window for TikZ Code'); +% legend('Location','northwest'); -legend('Location','northwest'); -set(gca,'FontSize',11); +%% 3. GENERATE TIKZ ANNOTATION CODE +% This prints the manual \draw commands to the console +%% GENERATE TIKZ ANNOTATION CODE +% This prints the manual \draw commands to the console + +%% GENERATE TIKZ ANNOTATION CODE (Colored Borders + Tiny Font) +%% GENERATE TIKZ ANNOTATION CODE (No Arrow, Close Text) +fprintf('\n\n%% ===========================================================\n'); +fprintf('%% COPY THE FOLLOWING LINES INTO YOUR .TEX FILE \n'); +fprintf('%% (Paste them just before \\end{axis})\n'); +fprintf('%% ===========================================================\n\n'); + +for i = 1:length(baud) + bx = baud(i); + by = netrate(i); + key = codes(i); + M_val = Mvals(i); + + % --- PLACEMENT LOGIC --- + if M_val == 8 + % PAM-8: Place Top-Left + % 'south east' anchor means the text's bottom-right corner touches the coordinate + % shift moves it slightly up and left to clear the marker + anchorStr = 'south east'; + shiftStr = 'shift={(-3pt, 3pt)}'; + else + % Others: Place Bottom-Right + % 'north west' anchor means the text's top-left corner touches the coordinate + % shift moves it slightly down and right + anchorStr = 'north west'; + shiftStr = 'shift={(3pt, -3pt)}'; + end + + % --- PRINT COMMAND --- + % Uses \node directly at the coordinate (axis cs:...) + fprintf('\\node[anchor=%s, %s, font=\\tiny, fill=white, inner sep=1pt] at (axis cs:%.2f, %.2f) {\\cite{%s}};\n', ... + anchorStr, shiftStr, bx, by, key); +end +fprintf('\n') %% === EXPORT === -outfile = 'C:\Users\Silas\Documents\latex\JLT_400G copy\media\matlab2tikz\highspeedresults.tikz'; +outfile = 'C:\Users\Silas\Documents\latex\JLT_400G_submission\media\matlab2tikz\highspeedresults_test.tikz'; + matlab2tikz(outfile, ... 'width','\fwidth', ... 'height','\fheight', ... diff --git a/projects/HighSpeedExperiment_2024/Auswertung_JLT/run_dsp_from_db.m b/projects/HighSpeedExperiment_2024/Auswertung_JLT/run_dsp_from_db.m index 8f96ba2..ab31451 100644 --- a/projects/HighSpeedExperiment_2024/Auswertung_JLT/run_dsp_from_db.m +++ b/projects/HighSpeedExperiment_2024/Auswertung_JLT/run_dsp_from_db.m @@ -42,14 +42,14 @@ fp = QueryFilter(); % fp.where('Runs', 'run_id','EQUALS', 2776); M = 6; fp.where('Runs', 'pam_level','EQUALS', M); -% fp.where('Runs', 'bitrate','EQUALS', 390e9);%360,390 +fp.where('Runs', 'bitrate','EQUALS', 360e9);%360,390 % fp.where('Runs', 'symbolrate','EQUALS', 195e9); -fp.where('Runs', 'fiber_length','EQUALS', 10); +fp.where('Runs', 'fiber_length','EQUALS', 2); fp.where('Runs', 'is_mpi','EQUALS', 0); % fp.where('Runs', 'interference_path_length','EQUALS', 1000); % fp.where('Runs', 'loop_id','GREATER_THAN', 11); % fp.where('Runs', 'sir','EQUALS',18); -% fp.where('Runs', 'wavelength','EQUALS', 1310); +fp.where('Runs', 'wavelength','EQUALS', 1310); fp.where('Runs', 'db_mode','EQUALS', 0); fp.where('Runs', 'rop_attenuation','EQUAL', 0); % fp.where('Runs', 'power_pd_in','LESS_THAN', 7); @@ -70,7 +70,7 @@ wh.addStorage("mlmlse_package"); %% === RUN IT === -[results,wh] = submitJobs(dataTable.run_id(:), dsp_options, "parallel", 'wh', wh, 'waitbar', true); +[results,wh] = submitJobs(dataTable.run_id(:), dsp_options, "serial", 'wh', wh, 'waitbar', true); diff --git a/projects/HighSpeedExperiment_2024/a_minimal_example.m b/projects/HighSpeedExperiment_2024/a_minimal_example.m index 966d5b9..2a94e13 100644 --- a/projects/HighSpeedExperiment_2024/a_minimal_example.m +++ b/projects/HighSpeedExperiment_2024/a_minimal_example.m @@ -4,70 +4,77 @@ dsp_options.storage_path = 'Z:\2024\sioe_labor\'; dsp_options.max_occurences = 1; db = DBHandler("dataBase", 'labor_highspeed', "type", 'mysql' ); -fp = QueryFilter(); - -fp.where('Runs','fiber_length','EQUALS', 2); -fp.where('Runs','wavelength','EQUALS', 1310); -fp.where('Runs','bitrate','EQUALS', 300e9); -fp.where('Runs','pam_level','EQUALS', 4); -fp.where('Runs','rop_attenuation','EQUALS', 0); -fp.where('Runs','is_mpi','EQUALS', 0); -fp.where('Runs', 'db_mode','EQUALS', 0); -% fields = db.getTableFieldNames('Runs'); -% [dataTable,~] = db.queryDB(fp, fields); -fields = db.getTableFieldNames('power_state_info'); -fields = [fields; db.getTableFieldNames('Runs')]; -fields = [fields; db.getTableFieldNames('dashboard_ungrouped_alltime')]; %dashboard_ungrouped_after_nov_2025 dashboard_ungrouped_aug_nov_2025 -fields = unique(fields); -[dataTable,~] = db.queryDB(fp, fields); - - -fsym = dataTable(1,:).symbolrate; -M = double(dataTable(1,:).pam_level); -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); - -% Preprocess signal -Scpe_sig = preprocessSignal(Scpe_cell{1}, Symbols, fsym); - -% Show spectrum -Scpe_sig.spectrum("fignum",1,"displayname",'Rx') - -% meta = struct(); -% meta.varnames = dataTable.Properties.VariableNames; +% fp = QueryFilter(); % -% for k = 1:numel(meta.varnames) -% v = meta.varnames{k}; -% col = dataTable.(v); +% fp.where('Runs','fiber_length','EQUALS', 2); +% fp.where('Runs','wavelength','EQUALS', 1310); +% fp.where('Runs','bitrate','EQUALS', 300e9); +% fp.where('Runs','pam_level','EQUALS', 4); +% fp.where('Runs','rop_attenuation','EQUALS', 0); +% fp.where('Runs','is_mpi','EQUALS', 0); +% fp.where('Runs', 'db_mode','EQUALS', 0); +% % fields = db.getTableFieldNames('Runs'); +% % [dataTable,~] = db.queryDB(fp, fields); % -% if isnumeric(col) || islogical(col) -% meta.(v) = col; -% elseif isstring(col) -% meta.(v) = cellstr(col); -% elseif iscellstr(col) -% meta.(v) = col; -% else -% error("Unsupported table column type: %s", class(col)) -% end -% end -% exp_data.metadata = meta; -exp_data = struct(); -exp_data.metadata = dataTable; -exp_data.tx_bits = Tx_bits.signal; -exp_data.tx_signal = Symbols.signal; -exp_data.rx_signal_2sps = Scpe_sig.signal; +run_ids = [ 993 1205 1413 1623 1833 2043 2253 2628 2836 2958 3000 3042 3323 5098 5225 5267 5309]; -fname = dataTable(1,:).rx_raw_path; -[~, filename, ext] = fileparts(fname); -filename = strrep(filename,"_raw_signal",""); -filename = filename + ext; -savepath = fullfile('F:\2024\sioe_labor\export_skuehl\',filename); -save(savepath,'exp_data','-v7.3'); +for id = run_ids + + fp = QueryFilter(); + fp.where('Runs', 'run_id','EQUALS', id); + fields = db.getTableFieldNames('power_state_info'); + fields = [fields; db.getTableFieldNames('Runs')]; + fields = [fields; db.getTableFieldNames('dashboard_ungrouped_alltime')]; %dashboard_ungrouped_after_nov_2025 dashboard_ungrouped_aug_nov_2025 + fields = unique(fields); + [dataTable,~] = db.queryDB(fp, fields); + + fsym = dataTable(1,:).symbolrate; + M = double(dataTable(1,:).pam_level); + 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); + + % Preprocess signal + Scpe_sig = preprocessSignal(Scpe_cell{1}, Symbols, fsym); + + % Show spectrum + Scpe_sig.spectrum("fignum",1,"displayname",'Rx') + + meta = struct(); + meta.varnames = dataTable.Properties.VariableNames; + + for k = 1:numel(meta.varnames) + v = meta.varnames{k}; + col = dataTable.(v); + + if isnumeric(col) || islogical(col) + meta.(v) = col; + elseif isstring(col) + meta.(v) = cellstr(col); + elseif iscellstr(col) + meta.(v) = col; + else + error("Unsupported table column type: %s", class(col)) + end + end + exp_data.metadata = meta; + + exp_data = struct(); + exp_data.metadata = dataTable; + exp_data.tx_bits = Tx_bits.signal; + exp_data.tx_signal = Symbols.signal; + exp_data.rx_signal_2sps = Scpe_sig.signal; + + fname = dataTable(1,:).rx_raw_path; + [~, filename, ext] = fileparts(fname); + filename = strrep(filename,"_raw_signal",""); + filename = filename + ext; + savepath = fullfile('F:\2024\sioe_labor\export_skuehl\',filename); + save(savepath,'exp_data','-v7.3'); +end diff --git a/projects/HighSpeedExperiment_2024/baudrate_evaluation.m b/projects/HighSpeedExperiment_2024/baudrate_evaluation.m index 0a48e57..100fa1d 100644 --- a/projects/HighSpeedExperiment_2024/baudrate_evaluation.m +++ b/projects/HighSpeedExperiment_2024/baudrate_evaluation.m @@ -1,5 +1,5 @@ -filename = "C:\Users\sioe\Documents\High_Speed_Measurement_2024\baudrate_sweep_b2b\PAMX_b2b_baudrate20241024_210648_wh.mat"; +filename = "F:\2024\sioe\High Speed Messungen Oktober\baudrate_sweep_b2b\PAMX_b2b_baudrate20241024_210648_wh_final.mat"; a = load(filename); wh = a.obj; diff --git a/projects/IMDD_base_system/minimal_example.m b/projects/IMDD_base_system/minimal_example.m index 89f24ff..32e1ef8 100644 --- a/projects/IMDD_base_system/minimal_example.m +++ b/projects/IMDD_base_system/minimal_example.m @@ -8,7 +8,7 @@ fdac = 256e9; fadc = 256e9; random_key = 1; -rcalpha = 0.05; +rcalpha = 0.6; kover = 16; duob_mode = db_mode.no_db; @@ -25,7 +25,7 @@ tx_bw_nyquist = 0.8; link_length = 1; % RX -rop = -8; +rop = -9; rx_bw_nyquist = 0.8; vnle_order1 = 50; @@ -51,7 +51,7 @@ mu_ffe = [mu_ffe1 mu_ffe3 mu_ffe3]; mu_dfe = 0.0004; -Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",16,"alpha",rcalpha); +Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"alpha",rcalpha); [Digi_sig,Symbols,Tx_bits] = PAMsource(... "fsym",fsym,"M",M,"order",18,"useprbs",0,... @@ -65,7 +65,7 @@ Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",1 %%%%% AWG El_sig = M8199A("kover",kover).process(Digi_sig); % El_sig = AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",0,"kover",kover,"bit_resolution",12,"upsampling_method","samplehold","precomp_sinc_rolloff",1).process(Digi_sig); -El_sig.spectrum("displayname",'Digi Spectrum','fignum',100,'normalizeTo0dB',0); +El_sig.spectrum("displayname",'Digi Spectrum','fignum',1,'normalizeTo0dB',1); % El_sig = El_sig.setPower(0,"dBm"); %%%%% Electrical Driver Amplifier %%%%%% @@ -101,32 +101,53 @@ Scpe_sig = Scope("fsimu",fdac*kover,"fadc",fadc,... "delay",0,"fixed_delay",0,"filtertype",filtertypes.butterworth,... "samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,... "adcresolution",8,"quantbuffer",0.1,'block_dc',1,'lpf_active',1,'H_lpf',Lp_scpe).process(Rx_sig); +%% -%%%%%% Sample to 2x fsym %%%%%% -Scpe_sig = Scpe_sig.resample("fs_out",2*fsym); -Scpe_sig.signal = Scpe_sig.signal(1:2*length(Symbols)); +% 1) matched filter +% pulse is symmetric, hence we can use pulsef firectly as matched filter. +% It feels off (bit I think correct) that the fsym is now the output freq.!! +% -> output 2 sps to omit timing recovery!? +Pform = Pulseformer("fsym",fsym,"fdac",2*fsym,"pulse","rrc","pulselength",16,"alpha",rcalpha,"matched",1); +Scpe_sig = Pform.process(Scpe_sig); +Scpe_sig.spectrum("displayname",'Signal after matched filter','fignum',1,'normalizeTo0dB',1); +% +% %% +% %%%%%% Sample to 2x fsym %%%%%% +% Scpe_sig = Scpe_sig.resample("fs_out",2*fsym); +% Scpe_sig.signal = Scpe_sig.signal(1:2*length(Symbols)); + +%% %%%%%% Sync Rx signal with reference %%%%%% -[Scpe_sig,~] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym,"debug_plots",0); +[Scpe_sig,~] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym,"debug_plots",1); Scpe_sig.spectrum("displayname",'Opt Spectrum','fignum',11,'normalizeTo0dB',1); -Scpe_sig = Filter('filtdegree',4,"f_cutoff",Symbols.fs.*0.5,"fs",Scpe_sig.fs,"filterType",filtertypes.gaussian,"active",true).process(Scpe_sig); +% Scpe_sig = Filter('filtdegree',4,"f_cutoff",Symbols.fs.*0.5,"fs",Scpe_sig.fs,"filterType",filtertypes.gaussian,"active",true).process(Scpe_sig); Scpe_sig = Scpe_sig - mean(Scpe_sig.signal); +Scpe_sig.signal = Scpe_sig.signal(1:2*length(Symbols)); + +%% % -------------------- FFE -------------------- ffe_order = [50, 0, 0]; -eq_ffe = EQ("Ne",ffe_order,"Nb",[0,0,0], ... +eq_ = EQ("Ne",ffe_order,"Nb",[2,0,0], ... "training_length",len_tr,"training_loops",5,"dd_loops",5, ... "K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... "FFEmu",0,"plotfinal",0,"ideal_dfe",0); -output.ffe_results = ffe(eq_ffe,M,Scpe_sig,Symbols,Tx_bits, ... - "precode_mode",duob_mode,'showAnalysis',0,"postFFE",[], ... +% eq_ = FFE("epochs_tr",4,"epochs_dd",5,"len_tr",4096,"mu_dd",0.01,"mu_tr",0.01,"order",50,"sps",2,"decide",0, "adaption",adaption_method.nlms,"dd_mode",1); +eq_ = FFE_DFE("epochs_tr",5,"epochs_dd",5,"len_tr",512,"ffe_mu_dd",1e-4,"dfe_mu_dd",5e-4,"ffe_mu_tr",0,"dfe_mu_tr",0,"ffe_order",99,"dfe_order",99,"sps",2,"decide",0); + + +output.ffe_results = ffe(eq_,M,Scpe_sig,Symbols,Tx_bits, ... + "precode_mode",duob_mode,'showAnalysis',1,"postFFE",[], ... "eth_style_symbol_mapping",0); output.ffe_results.metrics.print +%% + % -------------------- VNLE + MLSE -------------------- pf_ncoeffs = 1; ffe_order3 = [50, 5, 5]; @@ -139,8 +160,11 @@ pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); [output.vnle_results, output.mlse_results] = vnle_postfilter_mlse(eq_v, pf_, mlse_, M, Scpe_sig, Symbols, Tx_bits, ... - "precode_mode", duob_mode, 'showAnalysis', 0, "postFFE", [], "eth_style_symbol_mapping", 0); + "precode_mode", duob_mode, 'showAnalysis', 1, "postFFE", [], "eth_style_symbol_mapping", 0); +output.mlse_results.metrics.print + +%% % -------------------- DB target -------------------- mlse_db_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",M,'trellis_states',PAMmapper(M,0).levels); diff --git a/projects/Lab_analysis/data_analysis_and_plot_scripts/lambda_vs_spectrum.m b/projects/Lab_analysis/data_analysis_and_plot_scripts/lambda_vs_spectrum.m new file mode 100644 index 0000000..c14d43d --- /dev/null +++ b/projects/Lab_analysis/data_analysis_and_plot_scripts/lambda_vs_spectrum.m @@ -0,0 +1,88 @@ +wh_aeon = load("C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\Lab_analysis\aeon_soa_measurement_lambda_plaser_pump.mat"); +wh_aeon = wh_aeon.wh; +wh_thor = load("C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\Lab_analysis\thorlabs_pdfa_measurement_lambda_plaser_pump.mat"); +wh_thor = wh_thor.wh; + +% Silas' custom "warehouse" datatype +wh_aeon.showInfo; +wh_thor.showInfo; + +%% PLOT OSA SPECTRA + +figure();hold on +cols = linspecer(4); +cols = cbrewer2('Paired',8); +subplot(1,2,1);hold on; +wh = wh_thor; +ccnt = 1; + +pumps = wh.parameter.pump.values; +plasers = wh.parameter.laserpower.values; +scnt = 1; +for p = 1:numel(plasers) + + for pmp = 1:numel(pumps) + subplot(numel(plasers),numel(pumps),scnt);hold on; + title(sprintf('Laser: %d dBm; Pump: %d %',plasers(p),pumps(pmp)),"Interpreter","tex") + for lambda = wh.parameter.lambda.values + + spectrum_osa = wh_thor.getStoValue('spectrum_osa',plasers(p),lambda,pumps(pmp)); + wavelength_osa = wh_thor.getStoValue('wavelength_osa',plasers(p),lambda,pumps(pmp)); + plot(wavelength_osa,spectrum_osa,'DisplayName',sprintf('P_{in}: %d dB ',plasers(p)),'Color',cols(ccnt,:)); + + + spectrum_osa = wh_aeon.getStoValue('spectrum_osa',plasers(p),lambda,pumps(pmp)); + wavelength_osa = wh_aeon.getStoValue('wavelength_osa',plasers(p),lambda,pumps(pmp)); + plot(wavelength_osa,spectrum_osa,'DisplayName',sprintf('P_{in}: %d dB ',plasers(p)),'Color',cols(ccnt+1,:)); + + + ylim([-60, 20]); + beautifyBERplot("logscale",false,"setmarkers",0,"setcolors",0); + end + scnt = scnt +1; + end + ccnt = ccnt+2; +end + + +xlabel('Wavelength [nm]'); +ylabel('OSNR [dB]') + + +%% + +figure();hold on +cols = linspecer(4); +subplot(1,2,1);hold on; +wh = wh_thor; +ccnt = 1; + +pumps = wh.parameter.pump.values; +plasers = wh.parameter.laserpower.values; +scnt = 1; +for p = 1:numel(plasers) + + for pmp = 1:numel(pumps) + subplot(numel(plasers),numel(pumps),scnt);hold on; + + + ase_noise = wh_thor.getStoValue('pase_osa',plasers(p),wh.parameter.lambda.values,pumps(pmp)); + signal = wh_thor.getStoValue('psig_osa',plasers(p),wh.parameter.lambda.values,pumps(pmp)); + plot(wh.parameter.lambda.values,ase_noise,'DisplayName',sprintf('P_{in}: %d dB ',plasers(p)),'Color',cols(ccnt,:)); + + + ase_noise = wh_aeon.getStoValue('pase_osa',plasers(p),wh.parameter.lambda.values,pumps(pmp)); + signal = wh_aeon.getStoValue('psig_osa',plasers(p),wh.parameter.lambda.values,pumps(pmp)); + plot(wh.parameter.lambda.values,ase_noise,'DisplayName',sprintf('P_{in}: %d dB ',plasers(p)),'Color',cols(ccnt,:),'LineStyle','-'); + + ylim([-60, 20]); + beautifyBERplot("logscale",false,"setmarkers",0,"setcolors",0); + + scnt = scnt +1; + end + ccnt = ccnt+1; +end + + +xlabel('Wavelength [nm]'); +ylabel('OSNR [dB]') diff --git a/projects/ML_based_MLSE/analyze_filter_length.m b/projects/ML_based_MLSE/analyze_filter_length.m index 616b715..1b31d6b 100644 --- a/projects/ML_based_MLSE/analyze_filter_length.m +++ b/projects/ML_based_MLSE/analyze_filter_length.m @@ -5,9 +5,9 @@ M = 4; randkey = 1; % --- Parameter sweep -order_range = 2:3:11; % FFE order -delta_range = 0:2:4; % delta -SNR_dB = 20; +order_range = 5:5:50; % FFE order +delta_range = 0:5:20; % delta +SNR_dB = 30; % --- Prepare bit sequence order_bits = 19; @@ -21,10 +21,14 @@ Symbols = PAMmapper(M,0).map(Bits); Symbols.fs = 200e9; % --- Channel (minimal ISI + AWGN) -h = [0.3 0.9 0.3]; h = h/norm(h); +h = abs([0.3 0.9 0.3]); h = h/norm(h); + +% h = [1 -1.67085330039878 1.17918163282514 -0.805210559745616 0.571564213123367 -0.296337147529674 0.00649773445209780 0.0854177610195952 -0.0576009020965258 0.0520994427061551 -0.0624586034913656 0.0553280962699552 -0.00705582559925755 -0.0336399056707792 0.0706903719452810 -0.0334124287931977 0.0131699455037966 0.0587431373842994 -0.0515902976066452 0.00647904355473619 0.0137506750904990 -0.0547974515885928 0.00994735499340592 -0.0135513582534086 -0.00463322575007739 0.0277311946101940]; +% h = h/norm(h); symbols_filt = Symbols.filter(h,1); symbols_noi = symbols_filt; symbols_noi.signal = awgn(symbols_filt.signal,SNR_dB,'measured'); +symbols_noi.spectrum(); % --- Generate all parameter pairs [O,D] = ndgrid(order_range, delta_range); @@ -38,7 +42,7 @@ ce_vec = nan(size(pairs,1),1); ce_training = nan(size(pairs,1),training_len); % --- Parallel loop over parameter pairs -parfor k = 1:size(pairs,1) +for k = 1:size(pairs,1) order_k = pairs(k,1); delta_k = pairs(k,2); @@ -89,7 +93,7 @@ end beautifyBERplot ylabel('BER'); xlabel('Filter Order [N]'); title('BER vs. Filter order'); -ylim([1e-4, 0.1]); +% ylim([1e-4, 0.1]); yline(3.8e-3,'HandleVisibility','off'); yline(2.2e-4,'HandleVisibility','off'); diff --git a/projects/WDM/WDM_auswertung.m b/projects/WDM/WDM_auswertung.m index 98b42b9..1435d1f 100644 --- a/projects/WDM/WDM_auswertung.m +++ b/projects/WDM/WDM_auswertung.m @@ -1,5 +1,5 @@ -base = "C:\Users\Silas\Nextcloud\Cluster"; +base = "C:\Users\Silas\Nextcloud4\Cluster"; all_files = dir(fullfile(base, "**/*.mat")); schemes = ["co","pair","alt","seg"]; @@ -13,8 +13,8 @@ T = table('Size',[0 10], ... rx = "^WDM_(?\d{8})_(?