silas/imdd_silas
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

- Fixed 'ML_MLSE'

Browse Source 
- Added weighted DFE to 'EQ' (work in progress).
- Added the folders 'Evaluation Scripts' and 'Data' in 'projects\FSO_transmission', which contain different evaluations and measured data regarding the FSO transmission data.
This commit is contained in:
magf
2026-01-20 13:36:21 +01:00
parent 3b44315eea
commit 616945bf27
8 changed files with 635 additions and 52 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
Classes/00_signals/Signal.m
View File

@@ -172,9 +172,11 @@ classdef Signal
hold on; hold on;
if isempty(options.color) 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 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 end
% 2 c) % 2 c)
% - xlabel if not already here: time in readable format (1 ms and not 1e-3 s) % - xlabel if not already here: time in readable format (1 ms and not 1e-3 s)

48
Classes/04_DSP/Equalizer/EQ.m
View File

@@ -10,6 +10,10 @@ classdef EQ < handle
training_length %Number of training symbols training_length %Number of training symbols
training_loops %Number of loops through sequence for training mode training_loops %Number of loops through sequence for training mode
ideal_dfe %Error free DFE decisions 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 DB_aim %Aim at duobinary output sequence
@@ -68,6 +72,10 @@ classdef EQ < handle
options.training_length = 1024 %Number of training symbols options.training_length = 1024 %Number of training symbols
options.training_loops = 1 %Number of loops through sequence for training mode options.training_loops = 1 %Number of loops through sequence for training mode
options.ideal_dfe = 0 %Error free DFE decisions 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 options.DB_aim %Aim at duobinary output sequence
@@ -438,6 +446,43 @@ classdef EQ < handle
dd_out(k) = constellation_in_(dd_idx); dd_out(k) = constellation_in_(dd_idx);
end 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 if obj.Nb(1) > 0
dd_DFE(2:end) = dd_DFE(1:end-1); dd_DFE(2:end) = dd_DFE(1:end-1);
@@ -739,3 +784,6 @@ classdef EQ < handle
end end
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

513
Classes/04_DSP/Equalizer/ML_MLSE.m
View File

@@ -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 classdef ML_MLSE < handle
% --------------------------------------------------------------------- % ---------------------------------------------------------------------
% W. Lanneer and Y. Lefevre, % W. Lanneer and Y. Lefevre,
@@ -164,8 +660,8 @@ classdef ML_MLSE < handle
% EQUALIZE % EQUALIZE
% ============================================================== % ==============================================================
function [y,y_ref] = equalize(obj,x,d,mu,epochs,N,training) function [y,y_ref] = equalize(obj,x,d,mu,epochs,N,training)
debug = 0; debug = 1;
showPlots = 0; showPlots = 1;
y = zeros(N,1); y = zeros(N,1);
nSymbols = ceil(N/obj.sps); nSymbols = ceil(N/obj.sps);
@@ -241,6 +737,19 @@ classdef ML_MLSE < handle
dirac(trans_idx)=1; dirac(trans_idx)=1;
end 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) % TRAINING MODE (weight update)
% =================================================================== % ===================================================================

20
Classes/04_DSP/Timing_Recovery.m
View File

@@ -1,6 +1,7 @@
classdef Timing_Recovery < handle classdef Timing_Recovery < handle
properties(Access=public) properties(Access=public)
modulation
timing_error_detector timing_error_detector
sps sps
damping_factor damping_factor
@@ -9,14 +10,15 @@ classdef Timing_Recovery < handle
end end
methods(Access=public) methods(Access=public)
function obj = FFE(options) function obj = Timing_Recovery(options)
arguments(Input) arguments(Input)
options.timing_error_detector = 'Gardner'; options.modulation = 'PAM/PSK/QAM';
options.timing_error_detector = 'Gardner (non-data-aided)';
options.sps = 2; options.sps = 2;
options.damping_factor = 1.0; options.damping_factor = 1.0;
options.normalized_loop_bandwidth = 0.005; options.normalized_loop_bandwidth = 0.01;
options.detector_gain = 1; options.detector_gain = 2.7;
end end
@@ -25,21 +27,21 @@ classdef Timing_Recovery < handle
obj.(fn{n}) = options.(fn{n}); obj.(fn{n}) = options.(fn{n});
end end
obj.e = zeros(obj.order,1); %obj-Initialization here%
obj.error = 0;
end end
function data_out = process(obj, data_in) function [data_out, timing_error] = process(obj, data_in)
timing_synchronization = comm.SymbolSynchronizer( ... timing_synchronization = comm.SymbolSynchronizer( ...
"Modulation", obj.modulation,...
"TimingErrorDetector", obj.timing_error_detector, ... "TimingErrorDetector", obj.timing_error_detector, ...
"SamplesPerSymbol", obj.sps, ... "SamplesPerSymbol", obj.sps, ...
"DampingFactor", obj.damping_factor, ... "DampingFactor", obj.damping_factor, ...
"NormalizedLoopBandwidth", obj.normalized_loop_bandwidth, ... "NormalizedLoopBandwidth", obj.normalized_loop_bandwidth, ...
"DetectorGain", obj.detector_gain); "DetectorGain", obj.detector_gain);
data_out = data_in;
data_out.signal = timing_synchronization(data_in.signal); [data_out.signal, timing_error] = timing_synchronization(data_in.signal);
end end
end end

21
Functions/EQ_structures/ml_mlse.m
View File

@@ -46,8 +46,20 @@ end
ml_mlse_results.config = Equalizerstruct(); ml_mlse_results.config = Equalizerstruct();
eq_small = strip_eq(eq_, 10); 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); 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.eq = jsonencode(eq_);
ml_mlse_results.config.equalizer_structure = int32(equalizer_structure.ml_mlse); ml_mlse_results.config.equalizer_structure = int32(equalizer_structure.ml_mlse);
ml_mlse_results.config.comment = 'function: ML-based 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) 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 % Calculate BER based on precoding mode
mapper = PAMmapper(M, 0, "eth_style", eth_style); mapper = PAMmapper(M, 0, "eth_style", eth_style);
skip_front = 150;
switch precode_mode switch precode_mode
case db_mode.no_db case db_mode.no_db
@@ -95,11 +108,11 @@ switch precode_mode
tx_bits_precoded = mapper.demap(tx_symbols_precoded); tx_bits_precoded = mapper.demap(tx_symbols_precoded);
rx_bits = mapper.demap(eq_signal_hd_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 % B) Just determine BER
rx_bits = mapper.demap(eq_signal_hd); 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 case db_mode.db_precoded
% Data is precoded on TX side % 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().encode(eq_signal_hd, "M", M);
eq_signal_hd_decoded = Duobinary().decode(eq_signal_hd_decoded, "M", M); eq_signal_hd_decoded = Duobinary().decode(eq_signal_hd_decoded, "M", M);
rx_bits_decoded = mapper.demap(eq_signal_hd_decoded); 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 % B) Omit the Coding by comparing with demapped TX symbol sequence
tx_bits_demapped = mapper.demap(tx_symbols); tx_bits_demapped = mapper.demap(tx_symbols);
rx_bits = mapper.demap(eq_signal_hd); 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
end end

16
projects/FSO_transmission/first_analysis.m
View File

@@ -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 mode = 0; %0 oder 1
M = 2; M = 2;
all_files = dir(fullfile(base, "**/*.mat")); all_files = dir(fullfile(base, "**/*.mat"));
if M == 2 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"); 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 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"); 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 end
@@ -20,7 +20,7 @@ if mode == 1
end end
end end
tx_data = load(tx_data_path);
datas = load(filename); 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! % % 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_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.. %% not working..
% Rx_synced = Rx_Time_Rec;
Rx_synced = Rx_synced_cell{1}; Rx_synced = Rx_synced_cell{1};
len_tr = 4096*2; len_tr = 4096*2;
mu_ffe1 = 0.0001; mu_ffe1 = 0.0001;

60
projects/ML_based_MLSE/theoretic_channel_evaluation.m
View File

@@ -48,36 +48,36 @@ for i = 1:numel(SNR_dB)
symbols_noi = symbols_filt; symbols_noi = symbols_filt;
symbols_noi.signal = awgn(symbols_filt.signal, SNR_dB(i), 'measured'); % AWGN with given SNR symbols_noi.signal = awgn(symbols_filt.signal, SNR_dB(i), 'measured'); % AWGN with given SNR
% Sequence Est L=5 % % 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_ = 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; % mlse_.DIR = h;
[y_mlse] = mlse_.process(symbols_noi,Symbols); % [y_mlse] = mlse_.process(symbols_noi,Symbols);
mlse_bits = PAMmapper(M, 0, "eth_style", 0).demap(y_mlse); % 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); % [~, ~, 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)); % fprintf('MLSE L5: %.2e \n',ber_mlse_l5(i));
%
% 2nd Approach % % 2nd Approach
mu_lms = 0.0005; % mu_lms = 0.0005;
pf_ncoeffs = 1; % 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"); % 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); % pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1);
%
% FFE % % FFE
[y_ffe, ffe_noise] = eq_.process(symbols_noi, Symbols); % [y_ffe, ffe_noise] = eq_.process(symbols_noi, Symbols);
%
Eq_bits = PAMmapper(M, 0, "eth_style", 0).demap(y_ffe); % 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); % [~, ~, 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)); % fprintf('FFE: %.2e \n',ber_ffe(i));
%
% Postfilter % % Postfilter
[y_white,~] = pf_.process(y_ffe, ffe_noise); % [y_white,~] = pf_.process(y_ffe, ffe_noise);
%
% Sequence Est % % 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); % 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); % [y_mlse] = mlse_.process(y_white,Symbols);
mlse_bits = PAMmapper(M, 0, "eth_style", 0).demap(y_mlse); % 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); % [~, 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)); % fprintf('MLSE: %.2e \n',ber_nwf_mlse_l2(i));
% ML-base MLSE L=2 % ML-base MLSE L=2
adaptive_mu = 0; adaptive_mu = 0;

3
test/duobinary_minimal_example.m
View File

@@ -29,6 +29,9 @@ para.skip =0;
para.bruijn = 0; para.bruijn = 0;
para.reset_prms = 0; para.reset_prms = 0;
para.method = 1; para.method = 1;
para.pcs = 0;
para.shape_para = 0;
para.rng_num = 0;
data_in = []; data_in = [];
global loop; global loop;