diff --git a/Classes/00_signals/Signal.m b/Classes/00_signals/Signal.m index 5905556..8fef781 100644 --- a/Classes/00_signals/Signal.m +++ b/Classes/00_signals/Signal.m @@ -172,9 +172,11 @@ classdef Signal hold on; if isempty(options.color) - plot(t* 1e6, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Marker', '.', 'LineStyle','none', 'MarkerSize', 0.1); + % plot(t* 1e6, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Marker', '.', 'LineStyle','none', 'MarkerSize', 0.1); + plot(t* 1e6, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1); else - plot(t* 1e6, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Marker', '.', 'LineStyle','none', 'MarkerSize', 0.1,'Color',options.color); + % plot(t* 1e6, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Marker', '.', 'LineStyle','none', 'MarkerSize', 0.1,'Color',options.color); + plot(t* 1e6, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Color',options.color); end % 2 c) % - xlabel if not already here: time in readable format (1 ms and not 1e-3 s) diff --git a/Classes/04_DSP/Equalizer/EQ.m b/Classes/04_DSP/Equalizer/EQ.m index a03d66b..f9d58d2 100644 --- a/Classes/04_DSP/Equalizer/EQ.m +++ b/Classes/04_DSP/Equalizer/EQ.m @@ -10,6 +10,10 @@ classdef EQ < handle training_length %Number of training symbols training_loops %Number of loops through sequence for training mode ideal_dfe %Error free DFE decisions + weighted_DFE %Weighted DFE on/off + weighted_DFE_mode %Weighted DFE mode + weighted_DFE_d_min %d_min threshold parameter for weighted DFE mode 1 + weighted_DFE_I_mode %[a_s, b_s, I_max]-parameters for the weighted DFE DB_aim %Aim at duobinary output sequence @@ -68,6 +72,10 @@ classdef EQ < handle options.training_length = 1024 %Number of training symbols options.training_loops = 1 %Number of loops through sequence for training mode options.ideal_dfe = 0 %Error free DFE decisions + options.weighted_DFE = 0; + options.weighted_DFE_mode = 'R1'; + options.weighted_DFE_d_min = 0.5; + options.weighted_DFE_I_mode = [5,0.5,0.6]; options.DB_aim %Aim at duobinary output sequence @@ -438,6 +446,43 @@ classdef EQ < handle dd_out(k) = constellation_in_(dd_idx); end + % Implementation of a weighted DFE in + % order to prevent error propagation. + % For further details, study [1], chapter 3.2.2 - + % Modifications of DFE + + if obj.weighted_DFE + % define new constellations + const = unique(ref_in); + + % determine reliability factor gamma_k + if output_vec(m) > min(const) && output_vec(m) < max(const) + gamma_k = 1 - abs(output_vec(m) - dd_out(k)); + else + gamma_k = 1; + end + + % select mode + if strcmp(obj.weighted_DFE_mode,'R1') + if gamma_k >= obj.weighted_DFE_d_min + f_gamma_k = 1; + else + f_gamma_k = 0; + end + elseif strcmp(obj.weighted_DFE_mode,'R2') + f_gamma_k = gamma_k; + elseif strcmp(obj.weighted_DFE_mode,'I1') + nom = 1-exp(-obj.weighted_DFE_I_mode(1)*((gamma_k/obj.weighted_DFE_I_mode(2))-1)); + denom = 1+exp(-obj.weighted_DFE_I_mode(1)*((gamma_k/obj.weighted_DFE_I_mode(2))-1)); + f_gamma_k = (1/2)*((nom/denom) - 1); + elseif strcmp(obj.weighted_DFE_mode,'I2') + nom = 1-exp(-obj.weighted_DFE_I_mode(1)*((gamma_k/obj.weighted_DFE_I_mode(2))-1)); + denom = 1+exp(-obj.weighted_DFE_I_mode(1)*((gamma_k/obj.weighted_DFE_I_mode(2))-1)); + f_gamma_k = (obj.weighted_DFE_I_mode(3)/2)*((nom/denom) - 1); + end + % calculate weighted output + output_vec(m) = f_gamma_k.*dd_out(k)+(1-f_gamma_k).*output_vec(m); + end if obj.Nb(1) > 0 dd_DFE(2:end) = dd_DFE(1:end-1); @@ -739,3 +784,6 @@ classdef EQ < handle end end +% References +% [1] T. J. Wettlin, “Experimental Evaluation of Advanced Digital Signal Processing for Intra-Datacenter Systems using Direct-Detection,” 2023. [Online]. Available: https://nbn-resolving.org/urn:nbn:de:gbv:8:3-2023-00703-8 + diff --git a/Classes/04_DSP/Equalizer/ML_MLSE.m b/Classes/04_DSP/Equalizer/ML_MLSE.m index 3f584c3..9965525 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, @@ -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..3999ad7 100644 --- a/Classes/04_DSP/Timing_Recovery.m +++ b/Classes/04_DSP/Timing_Recovery.m @@ -1,6 +1,7 @@ classdef Timing_Recovery < handle properties(Access=public) + modulation timing_error_detector sps damping_factor @@ -9,14 +10,15 @@ classdef Timing_Recovery < handle end methods(Access=public) - function obj = FFE(options) + function obj = Timing_Recovery(options) arguments(Input) - options.timing_error_detector = 'Gardner'; + options.modulation = 'PAM/PSK/QAM'; + options.timing_error_detector = 'Gardner (non-data-aided)'; options.sps = 2; options.damping_factor = 1.0; - options.normalized_loop_bandwidth = 0.005; - options.detector_gain = 1; + options.normalized_loop_bandwidth = 0.01; + options.detector_gain = 2.7; end @@ -25,21 +27,21 @@ classdef Timing_Recovery < handle obj.(fn{n}) = options.(fn{n}); end - obj.e = zeros(obj.order,1); - obj.error = 0; + %obj-Initialization here% end - function data_out = process(obj, data_in) + function [data_out, timing_error] = process(obj, data_in) timing_synchronization = comm.SymbolSynchronizer( ... + "Modulation", obj.modulation,... "TimingErrorDetector", obj.timing_error_detector, ... "SamplesPerSymbol", obj.sps, ... "DampingFactor", obj.damping_factor, ... "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/ml_mlse.m b/Functions/EQ_structures/ml_mlse.m index f2d468d..80e28cf 100644 --- a/Functions/EQ_structures/ml_mlse.m +++ b/Functions/EQ_structures/ml_mlse.m @@ -46,8 +46,20 @@ 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); +fn = fieldnames(eq_); +for k = 1:numel(fn) + if issparse(eq_.(fn{k})) + eq_.(fn{k}) = full(eq_.(fn{k})); + end +end ml_mlse_results.config.eq = jsonencode(eq_); ml_mlse_results.config.equalizer_structure = int32(equalizer_structure.ml_mlse); ml_mlse_results.config.comment = 'function: ML-based MLSE'; @@ -81,6 +93,7 @@ end function [bits, errors, ber, error_pos, errors_precoded, ber_precoded] = calculateBER(eq_signal_hd, tx_symbols, tx_bits, precode_mode, M, eth_style) % Calculate BER based on precoding mode mapper = PAMmapper(M, 0, "eth_style", eth_style); +skip_front = 150; switch precode_mode case db_mode.no_db @@ -95,11 +108,11 @@ switch precode_mode tx_bits_precoded = mapper.demap(tx_symbols_precoded); rx_bits = mapper.demap(eq_signal_hd_precoded); - [~, errors_precoded, ber_precoded, ~] = calc_ber(rx_bits.signal, tx_bits_precoded.signal, "skip_front", 30000, "skip_end", 150, "returnErrorLocation", 1); + [~, errors_precoded, ber_precoded, ~] = calc_ber(rx_bits.signal, tx_bits_precoded.signal, "skip_front", skip_front, "skip_end", 150, "returnErrorLocation", 1); % B) Just determine BER rx_bits = mapper.demap(eq_signal_hd); - [bits, errors, ber, error_pos] = calc_ber(rx_bits.signal, tx_bits.signal, "skip_front", 30000, "skip_end", 150, "returnErrorLocation", 1); + [bits, errors, ber, error_pos] = calc_ber(rx_bits.signal, tx_bits.signal, "skip_front", skip_front, "skip_end", 150, "returnErrorLocation", 1); case db_mode.db_precoded % Data is precoded on TX side @@ -107,12 +120,12 @@ switch precode_mode eq_signal_hd_decoded = Duobinary().encode(eq_signal_hd, "M", M); eq_signal_hd_decoded = Duobinary().decode(eq_signal_hd_decoded, "M", M); rx_bits_decoded = mapper.demap(eq_signal_hd_decoded); - [~, errors_precoded, ber_precoded, ~] = calc_ber(rx_bits_decoded.signal, tx_bits.signal, "skip_front", 30000, "skip_end", 150, "returnErrorLocation", 1); + [~, errors_precoded, ber_precoded, ~] = calc_ber(rx_bits_decoded.signal, tx_bits.signal, "skip_front", skip_front, "skip_end", 150, "returnErrorLocation", 1); % B) Omit the Coding by comparing with demapped TX symbol sequence tx_bits_demapped = mapper.demap(tx_symbols); rx_bits = mapper.demap(eq_signal_hd); - [bits, errors, ber, error_pos] = calc_ber(rx_bits.signal, tx_bits_demapped.signal, "skip_front", 30000, "skip_end", 150, "returnErrorLocation", 1); + [bits, errors, ber, error_pos] = calc_ber(rx_bits.signal, tx_bits_demapped.signal, "skip_front", skip_front, "skip_end", 150, "returnErrorLocation", 1); end end diff --git a/projects/FSO_transmission/first_analysis.m b/projects/FSO_transmission/first_analysis.m index c47a2b1..42ea7b7 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\magf\Desktop\Desktop\MATLAB-Zeugs\FSO Equalizer\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_path = fullfile(base, "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_path = fullfile(base, "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 @@ -20,7 +20,7 @@ if mode == 1 end end - +tx_data = load(tx_data_path); datas = load(filename); %% @@ -97,9 +97,15 @@ plot(timing_error); % If this is a ramp, you have drift! % % try to do this with 2sps EQ! [~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_symbolsync.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1); - +% Rx_matched_1 = Rx_synced_cell{1}; +% +% % Rx_Time_Rec = Rx_matched; +% [Rx_Time_Rec, Timing_Error] = Timing_Recovery("modulation", 'PAM/PSK/QAM', "timing_error_detector",'Gardner (non-data-aided)','sps',2,'damping_factor',1,'normalized_loop_bandwidth',0.01,'detector_gain',2.7).process(Rx_matched_1); +% figure;plot(Timing_Error); +% % Rx_Time_Rec.signal = resample(Rx_Time_Rec.signal, 28e9, 14e9); %% not working.. +% Rx_synced = Rx_Time_Rec; Rx_synced = Rx_synced_cell{1}; len_tr = 4096*2; mu_ffe1 = 0.0001; diff --git a/projects/ML_based_MLSE/theoretic_channel_evaluation.m b/projects/ML_based_MLSE/theoretic_channel_evaluation.m index 3eb6f09..3e1efcc 100644 --- a/projects/ML_based_MLSE/theoretic_channel_evaluation.m +++ b/projects/ML_based_MLSE/theoretic_channel_evaluation.m @@ -48,36 +48,36 @@ for i = 1:numel(SNR_dB) symbols_noi = symbols_filt; symbols_noi.signal = awgn(symbols_filt.signal, SNR_dB(i), 'measured'); % AWGN with given SNR - % Sequence Est L=5 - mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels,'scale_mode',0,'trellis_exclusion',0,'trellis_state_mode',2,'debug',0,'DIR',h); - mlse_.DIR = h; - [y_mlse] = mlse_.process(symbols_noi,Symbols); - mlse_bits = PAMmapper(M, 0, "eth_style", 0).demap(y_mlse); - [~, ~, ber_mlse_l5(i), ~] = calc_ber(mlse_bits.signal, Bits.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1); - fprintf('MLSE L5: %.2e \n',ber_mlse_l5(i)); - - % 2nd Approach - mu_lms = 0.0005; - pf_ncoeffs = 1; - eq_ = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",2^13,"mu_dd",mu_lms,"mu_tr",mu_lms,"order",16,"sps",1,"dd_mode",1,"adaption_technique","lms"); - pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); - - % FFE - [y_ffe, ffe_noise] = eq_.process(symbols_noi, Symbols); - - Eq_bits = PAMmapper(M, 0, "eth_style", 0).demap(y_ffe); - [~, ~, ber_ffe(i), ~] = calc_ber(Eq_bits.signal, Bits.signal, "skip_front", 0, "skip_end", 0, "returnErrorLocation", 1); - fprintf('FFE: %.2e \n',ber_ffe(i)); - - % Postfilter - [y_white,~] = pf_.process(y_ffe, ffe_noise); - - % Sequence Est - mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels,'scale_mode',0,'trellis_exclusion',0,'trellis_state_mode',2,'debug',0,'DIR',pf_.coefficients); - [y_mlse] = mlse_.process(y_white,Symbols); - mlse_bits = PAMmapper(M, 0, "eth_style", 0).demap(y_mlse); - [~, errors, ber_nwf_mlse_l2(i), errpos] = calc_ber(mlse_bits.signal, Bits.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1); - fprintf('MLSE: %.2e \n',ber_nwf_mlse_l2(i)); + % % Sequence Est L=5 + % mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels,'scale_mode',0,'trellis_exclusion',0,'trellis_state_mode',2,'debug',0,'DIR',h); + % mlse_.DIR = h; + % [y_mlse] = mlse_.process(symbols_noi,Symbols); + % mlse_bits = PAMmapper(M, 0, "eth_style", 0).demap(y_mlse); + % [~, ~, ber_mlse_l5(i), ~] = calc_ber(mlse_bits.signal, Bits.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1); + % fprintf('MLSE L5: %.2e \n',ber_mlse_l5(i)); + % + % % 2nd Approach + % mu_lms = 0.0005; + % pf_ncoeffs = 1; + % eq_ = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",2^13,"mu_dd",mu_lms,"mu_tr",mu_lms,"order",16,"sps",1,"dd_mode",1,"adaption_technique","lms"); + % pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); + % + % % FFE + % [y_ffe, ffe_noise] = eq_.process(symbols_noi, Symbols); + % + % Eq_bits = PAMmapper(M, 0, "eth_style", 0).demap(y_ffe); + % [~, ~, ber_ffe(i), ~] = calc_ber(Eq_bits.signal, Bits.signal, "skip_front", 0, "skip_end", 0, "returnErrorLocation", 1); + % fprintf('FFE: %.2e \n',ber_ffe(i)); + % + % % Postfilter + % [y_white,~] = pf_.process(y_ffe, ffe_noise); + % + % % Sequence Est + % mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels,'scale_mode',0,'trellis_exclusion',0,'trellis_state_mode',2,'debug',0,'DIR',pf_.coefficients); + % [y_mlse] = mlse_.process(y_white,Symbols); + % mlse_bits = PAMmapper(M, 0, "eth_style", 0).demap(y_mlse); + % [~, errors, ber_nwf_mlse_l2(i), errpos] = calc_ber(mlse_bits.signal, Bits.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1); + % fprintf('MLSE: %.2e \n',ber_nwf_mlse_l2(i)); % ML-base MLSE L=2 adaptive_mu = 0; diff --git a/test/duobinary_minimal_example.m b/test/duobinary_minimal_example.m index c39183e..d61c044 100644 --- a/test/duobinary_minimal_example.m +++ b/test/duobinary_minimal_example.m @@ -29,6 +29,9 @@ para.skip =0; para.bruijn = 0; para.reset_prms = 0; para.method = 1; +para.pcs = 0; +para.shape_para = 0; +para.rng_num = 0; data_in = []; global loop;