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

Works at 1 SPS

Browse Source 
NWhitened Signal or FFE Output are okay.
mu is approx. 0.15
delta is 0
order is 4 (higher is not better yet)
This commit is contained in:
Silas Oettinghaus
2025-10-30 14:29:14 +01:00
parent ef4dc53db5
commit ac95aef1e0
1 changed files with 232 additions and 142 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -2,11 +2,11 @@ classdef ML_MLSE < handle
% Implementation of plain and simple FFE.
% 1) Training mode (stable performance when you use NLMS)
% 2) Decision directed mode
%
%LMS: mu in order of 0.0001 for acceptable convergence speed
%NLMS: mu in order of 0.01 for acceptable convergence speed
%RLS: mu is lambda -> 0.99 -> 1 (has a strong dependency on this! use a loop to find out best values)
%
% FFE("epochs_tr",5,"epochs_dd",2,"len_tr",2^13,"mu_dd",mu_dd,"mu_tr",mu_tr,"order",25,"sps",2,"decide",0, "adaption",adaption_method(adaption),"dd_mode",use_dd_mode);
properties
@@ -34,6 +34,20 @@ classdef ML_MLSE < handle
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
end
methods
@@ -51,6 +65,7 @@ classdef ML_MLSE < handle
options.mu_dd = 1e-5;
options.epochs_dd = 5;
options.delta = 0;
options.traceback_depth = 1024;
options.L = 1
@@ -79,12 +94,61 @@ classdef ML_MLSE < handle
warning('Signal length does not fit to reference!');
end
% ==============================================================
% INITIALIZATION (only before final epoch and detection mode)
% ==============================================================
% --- Parameters
obj.S = numel(unique(D.signal)); % alphabet size
obj.Nf = obj.order; % 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(unique(D.signal),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{:}).');
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 all(size(obj.w) ~= [obj.Nf+1,obj.nFeasible])
obj.w = zeros(obj.Nf+1,obj.nFeasible); % filter weights per transition + bias tap
end
% obj.w = randn(obj.Nf+1,obj.nFeasible);
% ==============================================================
% TRAINING
% ==============================================================
% Training Mode
n = obj.len_tr;
training = 1;
obj.equalize(X.signal, D.signal,obj.mu_tr,obj.epochs_dd,n,training);
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;
@@ -109,67 +173,17 @@ classdef ML_MLSE < handle
% ==============================================================
% FFE + Whitening + ML-Based Branch Metric Estimation + Viterbi
% ==============================================================
debug = 0;
% --- Input padding and preallocation
x = [zeros(floor(obj.order/2),1); x; zeros(obj.order,1)];
N_ = N / obj.sps;
y = zeros(N_,1);
y_white = zeros(N_,1);
y = zeros(N,1);
for epoch = 1:epochs
% ==============================================================
% INITIALIZATION (only before final epoch and detection mode)
% ==============================================================
if epoch == epochs
% --- Parameters
S = numel(unique(d)); % alphabet size
Nf = 9; % filter length
Delta = ceil(Nf/2); % delay parameter
nStates = S^obj.L;
nFeasible = nStates*S;
% --- Trellis mapping
obj.DIR = arburg(y-d(1:N_), obj.L);
obj.DIR_flip = flip(obj.DIR);
obj.trellis_states = reshape(unique(d),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));
combs = fliplr(combvec(pre_comb_cell{:}).');
first_sym = combs(:,1);
last_sym = combs(:,end);
nStates = size(combs,1);
% --- Valid transitions
valid = false(nStates);
for from = 1:nStates
for to = 1:nStates
if all(combs(to,2:end) == combs(from,1:end-1))
valid(to,from) = true;
end
end
end
[valid_to, valid_from] = find(valid);
% --- Noise estimation
y_ideal = conv(d(1:N_), obj.DIR(:), "same");
sigma2 = mean(abs(y - y_ideal).^2);
inv2s2 = 1/(2*sigma2);
% --- Allocate vectors and weights
pm = zeros(nStates,1);
w = zeros(Nf,nFeasible); % filter weights per transition
b = zeros(1,nFeasible); % bias terms
v_hat = zeros(1,nFeasible);
v_tilde = zeros(1,nFeasible);
bm_vec = zeros(1,nFeasible);
zi = zeros(max(numel(obj.DIR)-1,0),1);
mu_w = 0.001;
mu_b = 0.001;
end
pm = zeros(obj.nStates,1);
c_hat = zeros(1,obj.nFeasible);
v_tilde = zeros(1,obj.nFeasible);
pred = zeros(N, obj.nStates, 'uint32');
% ==============================================================
% RUNTIME LOOP
@@ -178,97 +192,173 @@ classdef ML_MLSE < handle
for sample = 1:obj.sps:N
symbol = symbol + 1;
% --- FFE output
U = x(obj.order+sample-1:-1:sample);
y(symbol,1) = obj.e.' * U;
k = symbol;
% --- Decision / FFE adaptation
if training
d_hat = d(symbol);
% --- Build Δ-delayed observation window y_k
i1 = k - obj.Nf + 1 + obj.delta;
i2 = k + 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 = (yk.' * obj.w); % [1×nFeasible]
c_hat = c_hat.'; % [nFeasible×1]
% --- Extended path metrics
v_tilde = pm(obj.valid_from_idx) + c_hat; % [nFeasible×1]
% ===== Gradient update (Algorithm 1) =====
% if training
% for current symbol index k -> previous (k-1) and current (k)
if k > obj.L
prev_seq = d(k-obj.L:k-1); % previous state symbols
true_from_state_idx = find(ismember(obj.combs, flip(prev_seq).', 'rows'));% find state indices in trellis
%if ~(true_from_state_idx == true_to_state_idx), warning('Impossible state transition?!'), end
curr_seq = d(k-obj.L+1:k); % next state symbols
true_to_state_idx = find(ismember(obj.combs, flip(curr_seq).', 'rows'));% find state indices in trellis
else
[~,idx] = min(abs(y(symbol) - obj.constellation));
d_hat = obj.constellation(idx);
% not enough history yet
true_from_state_idx = 1;
true_to_state_idx = 1;
end
err = d_hat - y(symbol);
obj.e = obj.e + mu * (err * U);
% --- Whitening + MLSE in last epoch
if epoch == epochs
[y_white(symbol), zi] = filter(obj.DIR,1,y(symbol), zi);
k = symbol;
% --- Build Δ-delayed observation window y_k
i1 = k - Nf + 1 + Delta;
i2 = k + Delta;
buf = y_white(max(1,i1):min(length(y_white),i2));
padL = max(0,1 - i1);
padR = max(0,i2 - length(y_white));
yk = [zeros(padL,1); buf(:); zeros(padR,1)]; % Nf×1
% --- Predict branch metrics for all feasible transitions
v_hat = (yk.' * w) + b; % [1×nFeasible]
v_hat = v_hat.'; % [nFeasible×1]
% --- Extended path metrics
v_tilde = pm(valid_from) + v_hat; % [nFeasible×1]
% ===== Gradient update (Algorithm 1) =====
% if training
% for current symbol index k -> previous (k-1) and current (k)
if k > obj.L
prev_seq = d(k-obj.L:k-1); % previous state symbols
curr_seq = d(k-obj.L+1:k); % next state symbols
% find state indices in trellis
true_from_state = find(ismember(combs, prev_seq.', 'rows'));
true_to_state = find(ismember(combs, curr_seq.', 'rows'));
else
% not enough history yet
true_from_state = 1;
true_to_state = 1;
end
% softmax over -v_tilde, one-hot target t
p = exp(-v_tilde);
p = p./sum(p);
t = zeros(nFeasible,1);
t(valid_from==true_from_state & valid_to==true_to_state) = 1;
delta = t - p; % ∂CE/∂(-v_tilde)
w = w + mu_w * (yk * delta.'); % Nf×nFeasible
b = b + mu_b * delta.'; % 1×nFeasible
% end
% =========================================
% compare–select to next states
pm_next = -inf(nStates,1);
surv_idx = zeros(nStates,1);
for s_to = 1:nStates
mask = (valid_to==s_to);
[pm_next(s_to), arg] = max(-v_tilde(mask)); % max likelihood ↔ min metric
surv_idx(s_to) = valid_from(find(mask,1,'first')-1+arg);
end
pm = pm_next;
% --- Traceback
if mod(symbol,obj.traceback_depth) == 0
[~,viterbi_path(symbol)] = max(pm);
for n = symbol:-1:symbol-obj.traceback_depth+2
viterbi_path(n-1) = surv_idx(viterbi_path(n));
end
end
if 0
disp(['FROM: state',char(num2str(true_from_state_idx)),' : symbol transition', char(num2str(obj.combs(true_from_state_idx,:)))]);
disp(['TO: state',char(num2str(true_to_state_idx)),' : symbol transition', char(num2str(obj.combs(true_to_state_idx,:)))]);
end
% true_from and true_to are (1) -> (-1)
% thus dirac = [1,0,0,0]'
dirac = zeros(obj.nFeasible,1);
dirac(obj.valid_from_idx==true_from_state_idx & obj.valid_to_idx==true_to_state_idx) = 1; % This Dirac delta function δ(s = s∗ k, s′ = s∗ k−1) = 1 if the extended state (s, s′) corresponds to the true realized states (s∗ k, s∗ k−1), and is zero otherwise.
if sum(dirac) == 0
warning('whats happening?!');
end
% softmax over -v_tilde, one-hot target t
% 4 feasible transitions: [0,0][0,1][1,0][1,1] #not in
% order here!
% first round:
% v_tilde is zero
% -> exp(-0) = 1
% -> 1/(1+1+1+1)
% -> 1/4
% -> each transition is equally likely?
% p = exp(-v_tilde);
p = exp(-(v_tilde-max(v_tilde)));
p = p./sum(p); % found in formula (9) and (19)
% 1-0.25 = 0.75
% 0-0.25 = -0.25
% what happens here at the zero? Is this some log
% probability - larger zero is likely; lower zero is
% unlikely? But this is based on known information (training)
dmp = (dirac - p)';
if mod(symbol,128) == 1 && debug
% --- Normalize and compute probabilities
v_norm = v_tilde - max(v_tilde); % numerical stability
probs_lin = exp(-v_norm);
probs_lin = probs_lin ./ sum(probs_lin);
probs_log = -v_norm;
% --- Map back into nStates×nStates grid
probs_mat = nan(obj.nStates, obj.nStates);
probs_mat(obj.valid) = probs_lin; % linear-space probabilities
probs_logmat = nan(obj.nStates, obj.nStates);
probs_logmat(obj.valid) = probs_log; % log-domain scores
% --- Identify the current true transition
[to_idx, from_idx] = find(obj.valid);
cur_idx = find(dirac==1);
cur_to = to_idx(cur_idx);
cur_from = from_idx(cur_idx);
% --- Plot using imagesc (supports hold)
figure(11); clf;
imagesc(probs_logmat);
axis xy; % origin top-left
% colormap(parula);
colorbar;
xlabel('From state');
ylabel('To state');
set(gca,'FontSize',10);
% --- Overlay current transition
hold on;
plot(cur_from, cur_to, 'rs', ...
'MarkerSize', 10, 'LineWidth', 2, 'MarkerFaceColor', 'none');
hold off;
end
if 1 %training
% the correct state gets a high update, we weight this
% with the input signal, from here the signals are not
% "understandable"
% dmp is large for the correct transition to update
% only this one!
dL_Dw = dmp .* (yk); % ∂CE/∂(w) - formula (10)
% from paper: We have observed in simulations that ignoring the derivative of vk−1(s′) during training yields negligible loss after convergence.
% my note: so the update direction is the same for b and w?
%only start with updates when we are inside the signal
if k > obj.L
% actual filter training updates:
obj.w = obj.w - mu * dL_Dw;% Nf×nFeasible
end
end
% compare select
v_tilde_mat = inf(obj.nStates, obj.nStates);
v_tilde_mat(obj.valid) = v_tilde;
[pm, pred(k,:)] = min(v_tilde_mat, [], 2);
pm_sto(:,symbol) = pm;
end
% --- Output reconstructed path
if epoch == epochs && ~training
y_vit = first_sym(viterbi_path);
[~, s_end] = max(pm);
viterbi_path = zeros(N,1,'uint32');
viterbi_path(N) = s_end;
for n = N:-1:2
viterbi_path(n-1) = pred(n, viterbi_path(n));
end
y_vit = obj.first_sym(viterbi_path);
y = obj.first_sym(viterbi_path);
if 1 %debug || training
err = sum(y ~= d(1:length(y)));
ser = err./length(y);
fprintf('Epoch: %d - SER: %.1e \n',epoch, ser);
figure(10);
subplot(2,2,1:2);
heatmap(obj.w);
title('Filter')
subplot(2,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(2,2,4);
% scatter(1:N,pm_sto,1,'.')
plot(1:N,pm_sto)
title('Path Metric Winners')
end
end
end
end
end
end