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

BCJR implementation

Browse Source 
WDM code added (Pol Cont., Opt MUX/DEMUX, Opt Atten, DP_Fiber) -> the codebase is not optimized to always work with dp signals!
This commit is contained in:
Silas Oettinghaus
2025-09-17 13:58:58 +02:00
parent f4a22d23a2
commit 4099f6820f
37 changed files with 3643 additions and 593 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
Classes/04_DSP/Coding/Duobinary.m
View File

@@ -53,6 +53,7 @@ classdef Duobinary
% bk(k+1) =mod(data(k)-bk(k),M);
% end
% THIS WAS USED!
for k = 2:numel(data)
bk(k) = mod(data(k)-bk(k-1),M);
end

261
Classes/04_DSP/Sequence Detection/MLSE.m
View File

@@ -29,9 +29,6 @@ classdef MLSE < handle
obj.(fn{n}) = options.(fn{n});
end
end
% do more stuff
end
function [signalclass_hd,LLR,GMI] = process(obj,signalclass,ref_symbolclass)
@@ -57,6 +54,26 @@ classdef MLSE < handle
function [VITERBI_ESTIMATION_SYMBOLS,LLR_maxlogmap,GMI] = process_(obj,data_in,data_ref)
debug = 0;
trellis_state_mode = 2; % General: States should match the target states of the prev. EQ (EQ's job was to reduce the error between signal and the target)
% 0 = use provided states (MUST provide the correct states);
% 1 = normalize to = 1 rms;
% 2 = use target symbols;
% 3 = use statistical levels
% 3 analyzes avg of rx signal levels - can help with nonlinear impairments
trellis_exclusion = 0; % PAM-6 only (only if data is NOT precoded!)
scale_mode = 2; % scale_mode:
% 0 = no scaling,
% 1 = RMS→scale MODEL,
% 2 = MMSE/time-corr→scale MODEL,
% 3 = RMS→scale DATA,
% 4 = MMSE/time-corr→scale DATA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% PREPARATIONS %%%%%%%%
% remove unnecessary zeros at start of impulse response to keep
% number of trellis states minimal
@@ -69,22 +86,42 @@ classdef MLSE < handle
obj.DIR = [0 obj.DIR];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% PREPARATIONS %%%%%%%%
%%%% Separate the equalized signal into the respective levels based on the actually transmitted level
constellation = unique(data_ref);
decisionLevels = (constellation(1:end-1) + constellation(2:end)) / 2;
N = length(data_in);
tx_bits = PAMmapper(obj.M,0,"eth_style",0).demap(data_ref);
% impulse respnse to remove from signal
obj.DIR = flip(obj.DIR); %i.e. -0.2676 -0.0478 1.0000
% Normalize the Trellis states to =1 RMS
obj.trellis_states = obj.trellis_states ./ rms(obj.trellis_states);
tx_bits = PAMmapper(obj.M,0,"eth_style",0).demap(data_ref);
% Trellis States
if trellis_state_mode == 1 % Normalize the Trellis states to =1 RMS
obj.trellis_states = obj.trellis_states ./ rms(obj.trellis_states);
elseif trellis_state_mode == 2 %simply use the states from the ref signal (should be a robust option)
obj.trellis_states = reshape(unique(data_ref),size(obj.trellis_states));
elseif trellis_state_mode == 3 %use_statistical_levels
%%%% Separate the equalized signal into the respective levels based on the actually transmitted level
constellation = unique(data_ref);
% find actual levels from rx signal
symbols_for_lvl = NaN(numel(constellation),length(data_ref));
for l = 1:numel(constellation)
level_amplitude = constellation(l);
symbols_for_lvl(l,data_ref==level_amplitude) = data_in(data_ref==level_amplitude);
end
%replace the trellis states
avg_levels = mean(symbols_for_lvl,2,'omitnan');
obj.trellis_states = sort(avg_levels)';
%also replace the whole ref signal (PAM-M) levels
[~, idx] = ismember(data_ref, unique(data_ref));
data_ref = avg_levels(idx);
end
% seems to be the only way to use combvec for a flexible amount
% of vectors. 'combs' contains all trellis states
@@ -97,74 +134,84 @@ classdef MLSE < handle
states = sum(combs,2);
nStates = length(last_sym);
% Calculate all possible input symbols for the desired impulse
% response. Row number is the index of the previous state,
% column number is the index of the next state
% noise free received == branch metrics
noise_free_received = zeros(nStates,nStates);
count_row = 1;
count_col = 1;
for l1 = 1:nStates
for l2 = 1:nStates
if sum(combs(l2,2:end) == combs(l1,1:end-1)) == size(combs,2)-1
noise_free_received(count_row,count_col) = sum(combs(l2,:).*obj.DIR(end:-1:2)) + last_sym(l1)*obj.DIR(1);
else
noise_free_received(count_row,count_col) = inf;
% % Calculate all possible input symbols for the desired impulse
% % response. Row number is the index of the previous state,
% % column number is the index of the next state
% % noise free received == branch metrics
% assumes: last_sym = combs(:,end); % already defined earlier
levels = sort(unique(obj.trellis_states(:)).');
edges = [levels(1) levels(end)]; % edge levels (0 and 5 in PAM6)
noise_free_received = inf(nStates,nStates); % rows: to, cols: from
edge_edge_mask = false(nStates,nStates); % rows: to, cols: from
for from = 1:nStates
for to = 1:nStates
% valid transition if shift-register overlap holds
if all(combs(to,2:end) == combs(from,1:end-1))
% noiseless sample for the 'to' state reached from 'from'
noise_free_received(to,from) = ...
dot(combs(to,:), obj.DIR(end:-1:2)) + last_sym(from)*obj.DIR(1);
% mark edge→edge candidate (to be excluded only on even→odd steps)
edge_edge_mask(to,from) = ...
(last_sym(from)==edges(1) || last_sym(from)==edges(2)) && ...
(last_sym(to) ==edges(1) || last_sym(to) ==edges(2));
end
count_row = count_row + 1;
end
count_col = count_col + 1;
count_row = 1;
end
% Was used earlier from Tom Wettlin etc. Sufficient for Viterbi
% but the BCJR is sensitive to the scaling, so I use another
% apporach
% % first: RMS normalization of input data (rms==1)
% data_in = data_in ./ rms(data_in);
%
% % then, match amplitude levels of input signal to those of the calculated ideal symbols
% % i.e. match the rms values of data_in to noise_free_received (rms=1.xx)
% if isreal(data_in)
% if obj.M == round(obj.M)
% data_in = data_in * rms(noise_free_received(noise_free_received ~= inf),'all','omitnan');
% end
% end
% quick and dirty alignment with xcorr
y = data_in;
x = data_ref;
h = flip(obj.DIR(:)).';
y_ideal = conv(x, h, "same");
[c,lags] = xcorr(y, y_ideal, 64); % small max lag is enough
[~,ix] = max(abs(c));
lag = lags(ix);
y_ideal = circshift(y_ideal, lag);
data_in = data_in(:);
y_ideal = conv(data_ref(:), h, "same");
% center, skip transients, rm mean
skip = max(100, numel(obj.DIR)+50);
y_prep = y(1+skip:end-skip) - mean(y(1+skip:end-skip));
y_ideal_prep = y_ideal(1+skip:end-skip) - mean(y_ideal(1+skip:end-skip));
% one-tap LS/MMSE gain and optional bias
g = (y_ideal_prep' * y_prep) / (y_ideal_prep' * y_ideal_prep); % complex allowed
b = mean(data_in) - g*mean(y_ideal); % intercept (if you keep means)
dp = dot(y_ideal_prep, y_prep - g*y_ideal_prep); %shall be close to zero
sigma2 = mean(abs(y_prep - g*y_ideal_prep).^2);
inv2s2 = 1/(2*sigma2);
debug = 0;
if debug
figure(100); hold on
plot(1:length(y),y,'DisplayName','Y: Received Signal','LineStyle','none','Marker','.','MarkerSize',1);
plot(1:length(y_ideal),y_ideal,'DisplayName','Y ideal: x * DIR','LineStyle','none','Marker','.','MarkerSize',1);
yline(noise_free_received(:),'DisplayName','All Transition States','HandleVisibility','off','Color','red');
yline(g*noise_free_received(:)+ b,'DisplayName','Scaled Transition States','HandleVisibility','off');
switch scale_mode
case 0
g = 1; b = 0;
case 1 % RMS: scale model to data
g = rms(data_in)/rms(y_ideal); b = mean(data_in) - g*mean(y_ideal);
case 2 % MMSE/time-corr: scale states to data
[c,lags] = xcorr(data_in(:), y_ideal, 64);
[~,ix] = max(abs(c));
lag = lags(ix);
y_ideal = circshift(y_ideal, lag);
mu_y = mean(data_in(:));
mu_i = mean(y_ideal);
y_c = data_in(:)-mu_y;
yi_c = y_ideal-mu_i;
g = (yi_c'*y_c)/(yi_c'*yi_c);
b = mu_y - g*mu_i;
case 3 % RMS flipped: scale data to model
gd = rms(y_ideal)/rms(data_in); bd = mean(y_ideal) - gd*mean(data_in);
data_in = gd*data_in + bd;
g = 1; b = 0;
case 4 % MMSE/time-corr flipped: scale data to states
[c,lags] = xcorr(data_in(:), y_ideal(:), 64);
[~,ix] = max(abs(c));
lag = lags(ix);
y_ideal = circshift(y_ideal(:), lag);
mu_y = mean(data_in(:));
mu_i = mean(y_ideal);
y_c = data_in(:) - mu_y; % data_in centered
yi_c = y_ideal - mu_i; % ideal centered
g = (y_c' * yi_c) / (y_c' * y_c);
b = mu_i - g * mu_y;
data_in = g * data_in(:) + b;
g = 1; b = 0;
end
noise_free_received = (g*noise_free_received + b);
% apply (g,b) to model and compute common sigma
noise_free_received = g*noise_free_received + b;
last_sym = g*last_sym + b;
sigma2 = mean(abs(data_in - (g*y_ideal + b)).^2);
inv2s2 = 1/(2*sigma2);
if debug
figure(100); clf; hold on
showLevelScatter(data_in, data_ref, "fignum", 100);
yline(noise_free_received(:), 'DisplayName','Transition States','Color','red','HandleVisibility','off');
yline(obj.trellis_states(:), 'DisplayName','Transition States','Color','green','LineWidth',2,'HandleVisibility','off')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% FORWARD PASS (VITERBI -Alpha's) %%%%%
@@ -185,8 +232,14 @@ classdef MLSE < handle
for n = 2:length(data_in)
bm = -(data_in(n) - noise_free_received).^2 * inv2s2;
% exclude edge→edge transitions only for even->odd steps && PAM-6
if mod(n,2) == 0 && obj.M == 6 && trellis_exclusion
bm(edge_edge_mask) = -Inf;
end
pm = pm + bm;
[alpha(:,n),pm_survivor_fw_idx(:,n)] = max(pm,[],2); % choose lowest path metric as new state
[alpha(:,n),pm_survivor_fw_idx(:,n)] = max(pm,[],2); % choose lowest path metric as new state (get min distance for all state transitions towards a new state)
pm = repmat(alpha(:,n).',nStates,1); % update pm (chosen state to 2nd dimension -> FROM state)
bm_fw(:,:,n) = bm;
@@ -202,6 +255,18 @@ classdef MLSE < handle
viterbi_path(n-1) = pm_survivor_fw_idx(viterbi_path(n),n);
end
if debug
alpha_ = alpha - min(alpha) + eps;
figure();hold on;
n = 10;
scatter(1:n,obj.trellis_states(repmat([1:numel(obj.trellis_states)]',1,n)),abs(alpha_(:,end-n+1:end)),'Marker','o','LineWidth',1);
scatter(1:n,obj.trellis_states(viterbi_path(end-n+1:end)),500,'Marker','x','LineWidth',1,'MarkerEdgeColor','green');
% scatter(1:n,data_ref(end-n+1:end),500,'Marker','x','LineWidth',1,'MarkerEdgeColor','red');
yticks(obj.trellis_states);
ylim([min(obj.trellis_states)-1 max(obj.trellis_states)+1]);
end
VITERBI_ESTIMATION_SYMBOLS(1:length(data_in)) = first_sym(viterbi_path);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -219,6 +284,12 @@ classdef MLSE < handle
for h = length(data_in)-1:-1:1
bm = -(data_in(h+1) - noise_free_received).^2 * inv2s2;
% exclude edge→edge transitions only for even->odd steps && PAM-6
if mod(h+1, 2) == 0 && obj.M == 6 && trellis_exclusion
bm(edge_edge_mask) = -Inf;
end
pm = pm + bm.';
[beta(:,h),pm_survivor_bw_idx(:,h)] = max(pm,[],2); % choose lowest path metric as new state
pm = repmat(beta(:,h).',nStates,1); % update pm (chosen state to 2nd dimension -> FROM state)
@@ -255,6 +326,7 @@ classdef MLSE < handle
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% FORWARD PASS PAM2,4,8 %%%%%
% These are interchangeable... second is chatgpt:
nml_LLP = LLP - max(LLP); %subtract highest value for better numerical stability, LLP's are not always close to zero
expLLP = exp(nml_LLP);
state_prob = expLLP ./ sum(expLLP); % sums to one (or numerically close to one)
@@ -266,9 +338,6 @@ classdef MLSE < handle
state_prob = exp(logPstate); % exact, sums to 1
if obj.M == 6
num_bits = 5;
@@ -357,17 +426,19 @@ classdef MLSE < handle
for bit_idx = 1:num_bits
% Find indices where bit is 0 and where it is 1
idx_sym_1 = bit_mapping(:,bit_idx) == 0;
idx_bit_0 = bit_mapping(:,bit_idx) == 0;
idx_bit_1 = bit_mapping(:,bit_idx) == 1;
% Sum over log-probabilities (Max-Log approximation: using max instead of sum)
LLR_maxlogmap(:,bit_idx) = max(LLP(idx_bit_1,:), [], 1) - max(LLP(idx_sym_1,:), [], 1);
% Sum over log-probabilities
% Max-Log approximation: using max instead of sum)
LLR_maxlogmap(:,bit_idx) = max(LLP(idx_bit_1,:), [], 1) - max(LLP(idx_bit_0,:), [], 1);
% Sum probabilities over states for which the bit is 1 and 0, respectively.
P0 = sum(state_prob(idx_sym_1, :),1);
P0 = sum(state_prob(idx_bit_0, :),1);
P1 = sum(state_prob(idx_bit_1, :),1);
LLR_exact(:,bit_idx) = log(P1./P0); % N x num_bits
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -377,12 +448,12 @@ classdef MLSE < handle
LLR_exact = LLR_exact;
for k = 1:num_bits
idx_bit_1 = (tx_bits(:,k) == 0); %wo sind die 1en
idx_sym_1 = (tx_bits(:,k) == 1); %wo sind die 0en
idx_bit_0 = (tx_bits(:,k) == 0); %wo sind die 1en
idx_bit_1 = (tx_bits(:,k) == 1); %wo sind die 0en
%LLR's for all actually transmitted ones or zeros
llr0 = LLR_exact(idx_bit_1,k);
llr1 = LLR_exact(idx_sym_1,k);
llr0 = LLR_exact(idx_bit_0,k);
llr1 = LLR_exact(idx_bit_1,k);
% Calculate mutual information for bit position k
I0 = mean(log2(1 + exp(llr0))); % exp(--LLR) = exp(positive) > 1
@@ -396,7 +467,6 @@ classdef MLSE < handle
VITERBI_ESTIMATION_SYMBOLS = VITERBI_ESTIMATION_SYMBOLS./rms(VITERBI_ESTIMATION_SYMBOLS);
if debug
%%% DEBUG PLOT LIKELIHOOD RATIOS %%%
figure(115);clf
@@ -415,6 +485,8 @@ classdef MLSE < handle
legend
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% CHECK BER's %%%%%
@@ -430,6 +502,13 @@ classdef MLSE < handle
fprintf('Viterbi BER: %.2e \n',ber_viterbi);
fprintf('Viterbi Errors = %d\n', numErr);
pairs = reshape(VITERBI_ESTIMATION_SYMBOLS,2,[]).';
levels = sort(unique(VITERBI_ESTIMATION_SYMBOLS));
isedge = ismember(pairs, [levels(1) levels(end)]);
isforbidden = sum(isedge,2)==2;
fprintf('Found %d forbidden transitions (even→odd edges).\n', nnz(isforbidden));
% Convert LLR values to a hard-decision bit stream
bit_stream = LLR_maxlogmap > 0; %ratio separates lower or higher than =0 -> simply decode for the negative values
bit_stream = reshape(bit_stream',[],1);
@@ -458,7 +537,7 @@ classdef MLSE < handle
rx_bits = PAMmapper(obj.M,0,"eth_style",0).demap(FW_EST');
rx_bits = reshape(rx_bits',[],1);
[~,~,ber_fw,~] = calc_ber(rx_bits,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
fprintf('FW BER: %.2e \n',ber_fw);
fprintf('FW BER : %.2e \n',ber_fw);
disp('Stop DEBUG MLSE:')
fprintf('\n')
end
@@ -471,7 +550,5 @@ classdef MLSE < handle
s = amax + log(sum(exp(a - amax), dim));
end
end
end

52
Classes/04_DSP/TransmissionPerformance.m
View File

@@ -42,6 +42,7 @@ classdef TransmissionPerformance
0.8733, 0.8790, 0.8848, 0.8905, 0.8962, ...
0.9019, 0.9077, 0.9134, 0.9191, 0.9248, ...
0.9306, 0.9363, 0.9420, 0.9477, 0.9535];
NGMITHRESHOLDS_SDHD = [0.8116, 0.8167, 0.8241, 0.8317, 0.8401, ...
0.8459, 0.8512, 0.8574, 0.8685, 0.8746, ...
0.8829, 0.8892, 0.8958, 0.9022, 0.9090,...
@@ -60,7 +61,7 @@ classdef TransmissionPerformance
%% LUT for KP4-FEC and Inner Code https://grouper.ieee.org/groups/802/3/dj/public/23_03/patra_3dj_01b_2303.pdf
CODE_RATE_KP4_AND_INNER = [0.885799];
CODE_RATE_KP4_AND_INNER = [0.885799]; %https://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=9979198 -> 199.8 / 224 = 0.891 als code rate
BERTHRESHOLDS_KP4_AND_INNER = 4.85e-3;
% https://www.ieee802.org/3/bs/public/14_11/parthasarathy_3bs_01a_1114.pdf
@@ -153,11 +154,21 @@ classdef TransmissionPerformance
netrates.SDHD.Threshold = NaN(1, numMeasurements);
end
if ~isempty(ber)
netrates.STAIR.GrossRate = NaN(1, numMeasurements);
netrates.STAIR.NetRate = NaN(1, numMeasurements);
netrates.STAIR.CodeRate = NaN(1, numMeasurements);
netrates.STAIR.Threshold = NaN(1, numMeasurements);
netrates.HD.GrossRate = NaN(1, numMeasurements);
netrates.HD.NetRate = NaN(1, numMeasurements);
netrates.HD.CodeRate = NaN(1, numMeasurements);
netrates.HD.Threshold = NaN(1, numMeasurements);
netrates.KP4.GrossRate = NaN(1, numMeasurements);
netrates.KP4.NetRate = NaN(1, numMeasurements);
netrates.KP4.CodeRate = NaN(1, numMeasurements);
netrates.KP4.Threshold = NaN(1, numMeasurements);
netrates.KP4_hamming.GrossRate = NaN(1, numMeasurements);
netrates.KP4_hamming.NetRate = NaN(1, numMeasurements);
netrates.KP4_hamming.CodeRate = NaN(1, numMeasurements);
@@ -200,11 +211,47 @@ classdef TransmissionPerformance
end
if ~isempty(idxBER)
codeRate = obj.CODE_RATES_HD(idxBER);
netrates.STAIR.NetRate(i) = grossRate(i) * codeRate;
netrates.STAIR.GrossRate(i) = grossRate(i) ;
netrates.STAIR.CodeRate(i) = codeRate;
netrates.STAIR.Threshold(i) = obj.BERTHRESHOLDS_HD(idxBER);
end
% HD FEC 3,8e-3
idxBER = [];
for j = length(obj.BERTHRESHOLDS_HDFEC):-1:1
if ber(i) <= obj.BERTHRESHOLDS_HDFEC(j)
idxBER = j;
break;
end
end
if ~isempty(idxBER)
codeRate = obj.CODE_RATE_HDFEC(idxBER);
netrates.HD.NetRate(i) = grossRate(i) * codeRate;
netrates.HD.GrossRate(i) = grossRate(i) ;
netrates.HD.CodeRate(i) = codeRate;
netrates.HD.Threshold(i) = obj.BERTHRESHOLDS_HD(idxBER);
netrates.HD.Threshold(i) = obj.CODE_RATE_HDFEC(idxBER);
end
% KP4
idxBER = [];
for j = length(obj.BERTHRESHOLDS_KP4):-1:1
if ber(i) <= obj.BERTHRESHOLDS_KP4(j)
idxBER = j;
break;
end
end
if ~isempty(idxBER)
codeRate = obj.CODE_RATE_KP4(idxBER);
netrates.KP4.NetRate(i) = grossRate(i) * codeRate;
netrates.KP4.GrossRate(i) = grossRate(i) ;
netrates.KP4.CodeRate(i) = codeRate;
netrates.KP4.Threshold(i) = obj.CODE_RATE_KP4(idxBER);
end
% KP4 + inner Hamming
idxBER = [];
for j = length(obj.BERTHRESHOLDS_KP4_AND_INNER):-1:1
if ber(i) <= obj.BERTHRESHOLDS_KP4_AND_INNER(j)
@@ -220,6 +267,7 @@ classdef TransmissionPerformance
netrates.KP4_hamming.Threshold(i) = obj.BERTHRESHOLDS_KP4_AND_INNER(idxBER);
end
% O FEC
idxBER = [];
for j = length(obj.BERTHRESHOLDS_O_FEC):-1:1
if ber(i) <= obj.BERTHRESHOLDS_O_FEC(j)