SD stuff working for PAM6 and MLSE is also calibrated

Classes/01_transmit/PAMmapper.m

@@ -29,7 +29,7 @@ classdef PAMmapper
 
             obj.levels = obj.get_levels();
 
-            obj.scaling = rms(obj.get_levels());
+            obj.scaling = obj.get_scaling;%rms(obj.get_levels());
 
             obj.eth_style = options.eth_style;
 
@@ -283,6 +283,21 @@ classdef PAMmapper
             end
         end
 
+        function scaling = get_scaling(obj)
+            switch obj.M
+                case 2
+                    scaling = 1;
+                case 4
+                    scaling = sqrt(5);
+                case 6
+                    scaling = sqrt(10);
+                case 8
+                    scaling = sqrt(21);
+                case 16
+                    scaling = 1;
+            end
+        end
+
         function [data_out] = demap_(obj,data_in)
             
             data_in= data_in';
@@ -450,7 +465,7 @@ classdef PAMmapper
         function [Signal_out] = quantize(obj,Signal_in)
 
             constellation = obj.get_levels();
-            constellation = constellation ./ rms(constellation);
+            constellation = constellation ./ obj.scaling;
 
             issignalclass = 0;
             if isa(Signal_in,'Signal')
@@ -481,7 +496,9 @@ classdef PAMmapper
         end
 
         function bitmap = showBitMapping(obj)
-            bitmap = obj.demap([obj.levels ./ obj.scaling]');
+
+            bitmap = obj.demap((obj.levels ./ obj.scaling)');
+
         end
 
     end

Classes/04_DSP/Sequence Detection/MLSE.m

@@ -83,11 +83,6 @@ classdef MLSE < handle
             % impulse respnse to remove from signal
             obj.DIR = flip(obj.DIR); %i.e. -0.2676   -0.0478    1.0000
 
-            % % make the combined impulse-response have net gain = 1
-            % h = obj.DIR(:);
-            % h = h / sum(h);
-            % obj.DIR = h.';
-
             % Normalize the Trellis states to =1 RMS
             obj.trellis_states = obj.trellis_states ./ rms(obj.trellis_states);
 
@@ -122,22 +117,55 @@ classdef MLSE < handle
                 count_row = 1;
             end
 
-            % first: RMS normalization of input data (rms==1)
-            data_in = data_in ./ rms(data_in);
-            data_in = data_in - mean(data_in);
+            % 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
 
-            % 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);
 
-            y_clean = conv( data_ref, flip(obj.DIR), "same" );
-            sigma2 = var( data_in - y_clean );
+            % 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');
+            end
+
+            noise_free_received = (g*noise_free_received + b);
+
             %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
             %%%%% FORWARD PASS (VITERBI -Alpha's) %%%%%
 
@@ -237,12 +265,9 @@ classdef MLSE < handle
             logPstate = LLP - logZ;          % still in log‐domain
             state_prob    = exp(logPstate);      % exact, sums to 1
 
-                % figure
-                % hold on;
-                % for i = 1:obj.M
-                %     scatter(1:length(expLLP),expLLP(i,:),1,'.');
-                % end
-                % scatter(1:length(expLLP),max(expLLP(:,:)),3,'.');
+
+
+            
 
             if obj.M == 6
 
@@ -250,7 +275,8 @@ classdef MLSE < handle
 
                 % all possible transitions (for now 36, including the "edges"
                 % of the QAM 32 constellation)
-                pam6transitions = combvec(obj.trellis_states,obj.trellis_states)'; % pam6transitions =
+                states = PAMmapper(6,0,"eth_style",0).levels;
+                pam6transitions = combvec(states,states)'; % pam6transitions =
                 % [-5    -5;
                 % -3    -5;
                 % -1    -5; ...
@@ -258,22 +284,13 @@ classdef MLSE < handle
                 pam6bits = PAMmapper(6,0,"eth_style",0).demap(reshape(pam6transitions',[],1)./sqrt(10));
 
                 pam6bits = reshape(pam6bits',5,[])';
-                %pam6bits =
-                %  0     0     0     1     0
-                %  0     0     0     1     0
-                %  0     0     0     1     1
-                %  1     0     0     1     1
-                %  1     0     0     1     0
-                %  1     0     0     1     0
-                %  0     0     1     1     0
-                % ....
 
-                [ok1, idx_bit_0] = ismember(pam6transitions(:,1), obj.trellis_states);
-                [ok2, idx2] = ismember(pam6transitions(:,2), obj.trellis_states);
+                [ok1, idx_sym_1] = ismember(pam6transitions(:,1), states);
+                [ok2, idx_sym_2] = ismember(pam6transitions(:,2), states);
                 assert(all(ok1)&all(ok2), 'Some transition amplitude not found in trellis_states')
-                pam6ind = [idx_bit_0, idx2];
+                pam6ind = [idx_sym_1, idx_sym_2];
 
-                tx_bits_pam6_reshaped = reshape(tx_bits,5,[])';
+                tx_bits_pam6_reshaped = reshape(tx_bits,5,[])'; % N x 5
 
                 numPairs = floor(size(LLP,2)/2);
                 LLR_exact     = zeros(numPairs,5);
@@ -294,16 +311,16 @@ classdef MLSE < handle
 
                     % now for each of the 5 bits do exact-LLR or max-log
                     for b = 1:num_bits
-                        idx_bit_0 = pam6bits(:,b)==1;
+                        idx_sym_1 = pam6bits(:,b)==1;
                         idx_bit_1 = pam6bits(:,b)==0;
 
                         %--- exact LLR from probabilities
-                        P1 = sum(pij(idx_bit_0));
+                        P1 = sum(pij(idx_sym_1));
                         P0 = sum(pij(idx_bit_1));
                         LLR_exact(k,b) = log(P1./P0);
 
                         %--- max-log:
-                        LLR_maxlogmap(k,b) = max( Mij(idx_bit_0) ) - max( Mij(idx_bit_1) ); % N x num_bits
+                        LLR_maxlogmap(k,b) = max( Mij(idx_sym_1) ) - max( Mij(idx_bit_1) ); % N x num_bits
                     end
                 end
 
@@ -311,11 +328,11 @@ classdef MLSE < handle
                 for k = 1:num_bits
 
                     idx_bit_1 = (tx_bits_pam6_reshaped(:,k) == 0); %wo sind die 1en
-                    idx_bit_0 = (tx_bits_pam6_reshaped(:,k) == 1); %wo sind die 0en
+                    idx_sym_1 = (tx_bits_pam6_reshaped(:,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_bit_0,k);
+                    llr1 = LLR_exact(idx_sym_1,k);
 
                     % Calculate mutual information for bit position k
                     I0 = mean(log2(1 + exp(llr0))); % exp(--LLR) = exp(positive) > 1
@@ -340,14 +357,14 @@ classdef MLSE < handle
                 for bit_idx = 1:num_bits
 
                     % Find indices where bit is 0 and where it is 1
-                    idx_bit_0 = bit_mapping(:,bit_idx) == 0;
+                    idx_sym_1 = 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_bit_0,:), [], 1);
+                    LLR_maxlogmap(:,bit_idx) = max(LLP(idx_bit_1,:), [], 1) - max(LLP(idx_sym_1,:), [], 1);
 
                     % Sum probabilities over states for which the bit is 1 and 0, respectively.
-                    P0 = sum(state_prob(idx_bit_0, :),1);
+                    P0 = sum(state_prob(idx_sym_1, :),1);
                     P1 = sum(state_prob(idx_bit_1, :),1);
                     LLR_exact(:,bit_idx) = log(P1./P0); % N x num_bits
 
@@ -361,11 +378,11 @@ classdef MLSE < handle
                 for k = 1:num_bits
 
                     idx_bit_1 = (tx_bits(:,k) == 0); %wo sind die 1en
-                    idx_bit_0 = (tx_bits(:,k) == 1); %wo sind die 0en
+                    idx_sym_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_bit_0,k);
+                    llr1 = LLR_exact(idx_sym_1,k);
 
                     % Calculate mutual information for bit position k
                     I0 = mean(log2(1 + exp(llr0))); % exp(--LLR) = exp(positive) > 1
@@ -379,7 +396,7 @@ classdef MLSE < handle
 
             VITERBI_ESTIMATION_SYMBOLS = VITERBI_ESTIMATION_SYMBOLS./rms(VITERBI_ESTIMATION_SYMBOLS);
 
-            debug = 1;
+
             if debug
                 %%% DEBUG PLOT LIKELIHOOD RATIOS %%%
                 figure(115);clf
@@ -403,7 +420,7 @@ classdef MLSE < handle
 
             if debug
                 tx_bits = reshape(tx_bits',[],1);
-
+                fprintf('\n')
                 disp('Start DEBUG MLSE:')
                 % DECIDE based on Viterbi traceback
                 VITERBI_ESTIMATION_SYMBOLS = VITERBI_ESTIMATION_SYMBOLS./rms(VITERBI_ESTIMATION_SYMBOLS);
@@ -443,7 +460,7 @@ classdef MLSE < handle
                 [~,~,ber_fw,~] = calc_ber(rx_bits,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
                 fprintf('FW BER: %.2e \n',ber_fw);
                 disp('Stop DEBUG MLSE:')
-                disp('')
+                fprintf('\n')
             end
 
         end

Functions/Metrics/calc_gmi_bitwise.m

@@ -1,64 +1,135 @@
 function [GMI,NGMI] = calc_gmi_bitwise(test_signal,reference_signal,options)
 %  https://cioffi-group.stanford.edu/doc/book/AppendixG.pdf
+%  Supports PAM2/4/8 (per-symbol) and PAM6 (pairwise: 5 bits per 2 symbols)
 
 arguments(Input)
-    test_signal;
-    reference_signal;
-    options.skip_front = 0;
-    options.skip_end = 0;
-    options.returnErrorLocation = 0;
+    test_signal
+    reference_signal
+    options.skip_front (1,1) double = 0
+    options.skip_end   (1,1) double = 0
+    options.returnErrorLocation (1,1) double = 0 %#ok<NASGU>
 end
 
 options.skip_end   = abs(options.skip_end);
 options.skip_front = abs(options.skip_front);
+assert((options.skip_end+options.skip_front) < numel(test_signal), ...
+    "You can not skip more samples than overall length.");
 
-assert((options.skip_end+options.skip_front)<length(test_signal),"You can not skip more bits than overall length of data! Set skip_front or skip_end to lower value or check data_in");
-
-if isa(reference_signal,'Signal')
-    reference_signal = reference_signal.signal;
-end
-
-if isa(test_signal,'Signal')
-    test_signal = test_signal.signal;
-end
+if isa(reference_signal,'Signal'); reference_signal = reference_signal.signal; end
+if isa(test_signal,'Signal');       test_signal     = test_signal.signal;     end
 
 % TRIM
 [test_signal,reference_signal] = trimseq(test_signal,reference_signal,options.skip_front,options.skip_end);
 
-% CALC AIR
-%%% new implementation of AIR
-% Precompute
+% Common precompute
 N       = length(test_signal);
 M       = numel(unique(reference_signal));
-nBits   = log2(M);
-levels  = PAMmapper(M,0).levels / PAMmapper(M,0).scaling;
-grayBits= PAMmapper(M,0).showBitMapping;
-priors  = ones(1,M)/M;
-
 noise   = test_signal - reference_signal;
-sigma2  = var(noise);
+sigma2  = var(noise);                       % real AWGN variance
+
+% --- Constellation (levels observed)
+constellation = unique(reference_signal);
+M = numel(constellation);
+
+% --- Empirical P_X (uniform is fine too; keep your original behavior)
+N = numel(reference_signal);
+counts = arrayfun(@(c) sum(reference_signal==c), constellation);
+priors = counts / N;
+
+% ----- PAM6: 5 bits mapped to 2 consecutive symbols (36 transitions) -----
+if M == 6
+    % Pair the stream (drop last if odd)
+    numPairs = floor(N/2);
+    if numPairs == 0, GMI = 0; NGMI = 0; return; end
+
+    y1 = test_signal(1:2:2*numPairs);
+    y2 = test_signal(2:2:2*numPairs);
+
+    % Use the same level order as PAMmapper to be consistent with labeling
+    levels = constellation;%PAMmapper(6,0,"eth_style",0).levels;     % 1x6 (actual amplitudes used)
+    % Build all 36 transitions (i,j)
+    [S1,S2] = ndgrid(levels, levels);                  % 6x6
+    trans_pairs = [S1(:), S2(:)];                      % 36x2
+
+    % 5-bit labels for each transition (ETH style).
+    % NOTE: keep the normalization you use with your PAMmapper (often /sqrt(10)).
+    bits72   = PAMmapper(6,0,"eth_style",0).demap( reshape(trans_pairs.',[],1)  );
+    pairBits = reshape(bits72.', 5, []).';             % 36 x 5
+
+    % Transmitted 5-bit labels for each observed pair (from the reference)
+    tx_bits_pair = PAMmapper(6,0,"eth_style",0).demap( reshape(reference_signal(1:2* numPairs).',[],1) );
+    tx_bits_pair = reshape(tx_bits_pair.', 5, []).';   % numPairs x 5
+
+    % Precompute symbol log-priors (uniform => constant; included for shaped inputs)
+    logP = log(priors);
+
+    % Precompute masks for bit=0 / bit=1 in the 6x6 grid
+    mask1 = false(6,6,5);
+    mask0 = false(6,6,5);
+    for b = 1:5
+        tmp = false(6,6);
+        tmp(:) = pairBits(:,b)==1;
+        mask1(:,:,b) = tmp;
+        tmp = false(6,6);
+        tmp(:) = pairBits(:,b)==0;
+        mask0(:,:,b) = tmp;
+    end
+
+    % Exact LLRs via log-sum-exp over the 6x6 grid
+    LLR = zeros(numPairs,5);
+    cst = -0.5*log(2*pi*sigma2);   % cancels in LLR but harmless
+    for k = 1:numPairs
+        li = cst - ((y1(k) - levels).^2)/(2*sigma2) + logP;  % 1x6
+        lj = cst - ((y2(k) - levels).^2)/(2*sigma2) + logP;  % 1x6
+        logW = li(:) + lj(:).';                               % 6x6 (adds logP twice)
+
+        for b = 1:5
+            L1 = logsumexp(logW(mask1(:,:,b)));
+            L0 = logsumexp(logW(mask0(:,:,b)));
+            LLR(k,b) = L1 - L0;                               % natural-log LLR
+        end
+    end
+
+    % Bit-wise MI using consistency relation
+    MI_bits = zeros(1,5);
+    for b = 1:5
+        idx0 = (tx_bits_pair(:,b)==0);
+        idx1 = ~idx0;
+        r0 = LLR(idx0,b);
+        r1 = LLR(idx1,b);
+        I0 = mean( log2(1 + exp( r0)) );     % natural LLR inside exp()
+        I1 = mean( log2(1 + exp(-r1)) );
+        MI_bits(b) = 1 - 0.5*(I0 + I1);
+    end
+
+    % Per-symbol outputs (5 bits per 2 symbols)
+    GMI  = sum(MI_bits)/2;
+    NGMI = GMI / 2.5;
+
+% ----- Generic PAM2/4/8…: per-symbol, Gray bits directly on symbols -----
+else
+    nBits   = log2(M);
+    levels  = PAMmapper(M,0).levels' / PAMmapper(M,0).scaling;  % match PAMmapper ordering
+    grayBits= PAMmapper(M,0).showBitMapping;                    % M x nBits
 
     % Allocate
     LLR     = zeros(N, nBits);
-rxBits  = zeros(N, nBits);
 
     for n = 1:N
         y = test_signal(n);
-    ll_table= -((y - levels).^2)/(2*sigma2) + log(priors);  % M×1
-
-    % Per‐bit LLR
+        ll_table = -((y - levels).^2)/(2*sigma2) + log(priors);   % 1xM (log-domain)
         for b = 1:nBits
             idx0 = grayBits(:,b)==0;
-        idx1 = grayBits(:,b)==1;
+            idx1 = ~idx0;
             L0   = logsumexp(ll_table(idx0));
             L1   = logsumexp(ll_table(idx1));
-        LLR(n,b) = L1 - L0;
+            LLR(n,b) = L1 - L0;                                   % natural-log LLR
         end
     end
 
     tx_bits = PAMmapper(M,0,"eth_style",0).demap(reference_signal);
 
-% Compute GMI (bit‐wise MI averaged over all symbols)
+    % Bit-wise MI
     MI_bits = zeros(1,nBits);
     for b = 1:nBits
         r0 = LLR(tx_bits(:,b)==0,b);
@@ -68,27 +139,20 @@ for b = 1:nBits
         MI_bits(b) = 1 - 0.5*(I0 + I1);
     end
 
-GMI = sum(MI_bits);            % in bits per symbol
-NGMI = GMI / nBits;        % normalized per bit
-
-% Auxiliary nested helper for numerically stable log-sum-exp
-    function s = logsumexp(a)
-        % LOGSUMEXP   Compute log(sum(exp(a))) in a numerically stable way
-        m = max(a);
-        s = m + log(sum(exp(a - m)));
+    GMI  = sum(MI_bits);          % bits / symbol
+    NGMI = GMI / nBits;
 end
 
+% ===== Helpers =====
+function s = logsumexp(a)
+    if isempty(a), s = -inf; return; end
+    m = max(a(:));
+    s = m + log(sum(exp(a(:) - m)));
+end
 
 function [data_,reference_] = trimseq(data,reference,skipstart,skip_end)
-
     data_  = data(skipstart+1:end-skip_end,:);
-
-        delta_bits = length(reference) - length(data);
-
-        skip_end = delta_bits + skip_end;
-
-        reference_ = reference(skipstart+1:end-skip_end,:);
-
+    delta  = numel(reference) - numel(data_);
+    reference_ = reference(skipstart+1:end-(skip_end+delta),:);
 end
-
 end

Functions/Metrics/calc_ngmi.m

@@ -1,125 +1,181 @@
 function [GMI,NGMI] = calc_ngmi(test_signal,reference_signal,options)
-% Silas implementation of (N)GMI calculation according to: J. Cho, L. Schmalen, und P. J. Winzer,
-% „Normalized Generalized Mutual Information as a Forward Error Correction Threshold for Probabilistically Shaped QAM“,
-% in 2017 European Conference on Optical Communication (ECOC), Sep. 2017, doi: 10.1109/ECOC.2017.8345872.
-
-% This implementation assumes the same normal distributed noise for each
-% channel (sigma2 is calculated once for the whole signal)
+% (N)GMI for AWGN, supports PAM2/4/8 (per-symbol) and PAM6 (pair-based, 5 bits / 2 symbols)
 
 arguments(Input)
-    test_signal;
-    reference_signal;
-    options.skip_front = 0;
-    options.skip_end = 0;
-    options.returnErrorLocation = 0;
+    test_signal
+    reference_signal
+    options.skip_front (1,1) double = 0
+    options.skip_end   (1,1) double = 0
+    options.returnErrorLocation (1,1) double = 0 %#ok<NASGU>
 end
 
 options.skip_end   = abs(options.skip_end);
 options.skip_front = abs(options.skip_front);
+assert((options.skip_end+options.skip_front) < numel(test_signal), ...
+    "You can not skip more samples than the overall length.");
 
-assert((options.skip_end+options.skip_front)<length(test_signal),"You can not skip more bits than overall length of data! Set skip_front or skip_end to lower value or check data_in");
+if isa(reference_signal,'Signal'), reference_signal = reference_signal.signal; end
+if isa(test_signal,'Signal'),     test_signal     = test_signal.signal;     end
 
-if isa(reference_signal,'Signal')
-    reference_signal = reference_signal.signal;
-end
-
-if isa(test_signal,'Signal')
-    test_signal = test_signal.signal;
-end
-
-% TRIM
+% --- Trim head/tail consistently
 [test_signal,reference_signal] = trimseq(test_signal,reference_signal,options.skip_front,options.skip_end);
 
-% CALC EVM
-[GMI,NGMI] = calc_ngmi_(test_signal,reference_signal);
+% --- Noise variance (real AWGN): use mean-square error for robustness
+err    = test_signal - reference_signal;
+sigma2 = mean((err(:)).^2);
 
-
-    function [GMI,NGMI] = calc_ngmi_(test_signal,reference_signal)
-
-        assert(length(test_signal) == length(reference_signal),"Sequence length does not match");
-
-        %%% implemented according to [1] J. Cho, L. Schmalen, und P. J. Winzer,
-        % „Normalized Generalized Mutual Information as a Forward Error Correction Threshold for Probabilistically Shaped QAM“,
-        % in 2017 European Conference on Optical Communication (ECOC), Sep. 2017, S. 1–3. doi: 10.1109/ECOC.2017.8345872.
-
-        error_vector = (test_signal-reference_signal);
-        sigma2 = var(error_vector); %noise variance
-
-        %%% Separate Classes
+% --- Constellation (levels observed)
 constellation = unique(reference_signal);
-        received_sd = NaN(numel(constellation),length(reference_signal));
-        lvlcol = cbrewer2('Set1',numel(constellation));
-        for lvl = 1:numel(constellation)
-            %Separate the equalized signal into the
-            %respective levels based on the actually
-            %transmitted level!
-            received_sd(lvl,reference_signal==constellation(lvl)) = test_signal(reference_signal==constellation(lvl));
+M = numel(constellation);
+
+% --- Empirical P_X (uniform is fine too; keep your original behavior)
+N = numel(reference_signal);
+counts = arrayfun(@(c) sum(reference_signal==c), constellation);
+P_X = counts / N;
+
+% --- Dispatch: PAM6 needs pairwise treatment (5 bits across 2 symbols)
+if M==6
+    % ===== PAM6: 5 bits mapped to two consecutive symbols =====
+    % Pair the sequence (drop last symbol if odd)
+    numPairs = floor(N/2);
+    if numPairs == 0
+        GMI = 0; NGMI = 0; return;
+    end
+    y1 = test_signal(1:2:2*numPairs);
+    y2 = test_signal(2:2:2*numPairs);
+    x1 = reference_signal(1:2:2*numPairs);
+    x2 = reference_signal(2:2:2*numPairs);
+
+    % Build the 36 transition table and its 5-bit labels (ETH-style)
+    symbols = sort(constellation); % keep fixed order
+    % all pairs using ndgrid (no dependency on combvec)
+    [S1,S2] = ndgrid(symbols, symbols);             % MxM
+    trans_pairs = [S1(:), S2(:)];                   % (M^2) x 2
+
+    % Map those pairs to 5-bit labels using your PAMmapper (same as in your BCJR code)
+    % NOTE: keep your normalization used by PAMmapper (sqrt(10)) if that's what it expects.
+    bits72 = PAMmapper(6,0,"eth_style",0).demap( reshape(trans_pairs',[],1) );
+    pam6bits = reshape(bits72', 5, []).';           % (M^2) x 5, row-aligned with trans_pairs
+
+    % Precompute mapping (pair -> bit label) and (symbol -> index)
+    [ok1, idx_s1] = ismember(trans_pairs(:,1), symbols);
+    [ok2, idx_s2] = ismember(trans_pairs(:,2), symbols);
+    assert(all(ok1)&all(ok2), 'Transition level not found in symbols.');
+    % For fast masking: linear indices of each (i,j) pair in an MxM grid
+    lin_idx_pairs = sub2ind([M,M], idx_s1, idx_s2); % (M^2) x 1
+
+    % Build pairwise prior P_pair = P_X(i)*P_X(j) on the MxM grid
+    P_pair = P_X(:) * P_X(:).';                     % MxM
+    logP_pair = log(P_pair);                        % for numerical stability
+
+    % Precompute constant term for log-likelihood
+    logc = -0.5*log(2*pi*sigma2);
+
+    % Prepare transmitted 5-bit labels per observed pair (to split llr0/llr1)
+    % Find each (x1(k),x2(k)) in trans_pairs:
+    [tf, row_in_table] = ismember([x1(:), x2(:)], trans_pairs, 'rows');
+    assert(all(tf), 'A reference pair was not found in the transition table.');
+    tx_bits_pair = pam6bits(row_in_table, :);       % numPairs x 5
+
+    % Precompute masks for bit=0 / bit=1 on the MxM grid for each bit position
+    mask1 = false(M,M,5);
+    mask0 = false(M,M,5);
+    for b = 1:5
+        tmp = false(M,M);
+        tmp(lin_idx_pairs) = pam6bits(:,b) == 1;
+        mask1(:,:,b) = tmp;
+        tmp = false(M,M);
+        tmp(lin_idx_pairs) = pam6bits(:,b) == 0;
+        mask0(:,:,b) = tmp;
     end
 
-        N = length(test_signal);           % Number of received samples
+    % Loop pairs: compute exact LLRs via log-sum-exp over the MxM grid
+    LLR_exact = zeros(numPairs,5);
+    for k = 1:numPairs
+        % log-likelihood vectors along each axis
+        li = logc - ((y1(k) - symbols).^2) / (2*sigma2);  % 1xM
+        lj = logc - ((y2(k) - symbols).^2) / (2*sigma2);  % 1xM
+        % joint log-weights over all (i,j)
+        logW = (li(:) + lj(:).') + logP_pair;             % MxM
 
-        M = length(constellation); %P
-        m = log2(M); %bits per symbol
-
-        entries = sum(~isnan(received_sd),2)';
-        P_X = ones(1,M)/M;%entries./N;
-
-        % Parameters
-        symbols = constellation';     % PAM-4 symbols
-
-        gray_bits = PAMmapper(M,0).demap_(constellation);       % Gray coding (bits per symbol)
-
-        % Conditional probability function for AWGN
-        q_Y_given_X = @(y, x) (1 / sqrt(2 * pi * sigma2)) * exp(-(y - x).^2 / (2 * sigma2));
-
-        % Entropy term
-        H_X = -sum(P_X .* log2(P_X)); % Entropy of input distribution
-
-        tx_bits = PAMmapper(M,0,"eth_style",0).demap(reference_signal);
-
-        % GMI computation
-        bit_llr_sum  = 0;
-        for k = 1:N %loop over signal
-            y_k = test_signal(k); % Current received sample
-            % [~, closest_symbol_idx] = min(abs(symbols - y_k)); % Closest symbol index
-
-            for i = 1:m %loop over bit position
-                % Extract i-th bit for each symbol
-                bit_mask = gray_bits(:, i); % Binary column for i-th bit of all symbols
-
-                % matching_symbols = symbols(bit_mask == gray_bits(closest_symbol_idx, i)); %old, based on decision that mihht be wrong
-                matching_symbols = symbols(bit_mask == tx_bits(k,i)); %new, based on tx bits
-
-                % Numerator: Sum over x in x_{b_{k, i}}
-                numerator = sum(q_Y_given_X(y_k, matching_symbols) .* P_X(ismember(symbols, matching_symbols)));
-
-                % Denominator: Sum over all x
-                denominator = sum(q_Y_given_X(y_k, symbols) .* P_X);
-
-                % Logarithmic contribution
-                bit_llr_sum = bit_llr_sum + log2(numerator / denominator);
+        for b = 1:5
+            % log P(bit=1 | y) ∝ logsumexp(logW over mask1)
+            lnum = logsumexp(logW(mask1(:,:,b)));
+            % log P(bit=0 | y) ∝ logsumexp(logW over mask0)
+            lden = logsumexp(logW(mask0(:,:,b)));
+            LLR_exact(k,b) = lnum - lden;                 % natural-log LLR
         end
     end
 
-        % Normalize the noise impact term by N
-        bit_llr_sum = bit_llr_sum / N;
+    % --- Bitwise MI via consistency relation
+    MI = zeros(1,5);
+    for b = 1:5
+        llr_b = LLR_exact(:,b);
+        idx0 = (tx_bits_pair(:,b)==0);
+        idx1 = ~idx0;
+        I0 = mean( log2(1 + exp( llr_b(idx0))) );   % natural LLR inside exp()
+        I1 = mean( log2(1 + exp(-llr_b(idx1))) );
+        MI(b) = 1 - 0.5*(I0 + I1);
+    end
 
-        % GMI
-        GMI = H_X + bit_llr_sum;
+    % GMI per symbol (5 bits per 2 symbols)
+    GMI  = sum(MI)/2;
+    NGMI = GMI / 2.5;
+
+else
+    % ===== Generic PAM2/4/8 etc.: per-symbol bit labeling =====
+    symbols = constellation;                         % 1xM
+    % Gray labels per symbol (your helper)
+    gray_bits = PAMmapper(M,0).demap_(symbols);     % M x log2(M)
+    m = log2(M);
+
+    % Likelihood (real AWGN)
+    q = @(y,x) (1/sqrt(2*pi*sigma2)) * exp(-(y - x).^2/(2*sigma2));
+
+    % Transmitted bits per sample
+    tx_bits = PAMmapper(M,0,"eth_style",0).demap(reference_signal); % N x m
+
+    % Exact bit-LLRs and MI
+    LLR_exact = zeros(N,m);
+    for b = 1:m
+        mask1 = (gray_bits(:,b)==1);
+        mask0 = ~mask1;
+        for k = 1:N
+            yk = test_signal(k);
+            den = sum(q(yk, symbols) .* P_X);
+            num1 = sum(q(yk, symbols(mask1)) .* P_X(mask1));
+            num0 = sum(q(yk, symbols(mask0)) .* P_X(mask0));
+            % Use the ratio of posteriors P(bit=1|y)/P(bit=0|y)
+            LLR_exact(k,b) = log(num1) - log(num0); % natural log
+        end
+    end
+
+    % Bitwise MI
+    MI = zeros(1,m);
+    for b = 1:m
+        llr_b = LLR_exact(:,b);
+        idx0 = (tx_bits(:,b)==0);
+        idx1 = ~idx0;
+        I0 = mean( log2(1 + exp( llr_b(idx0))) );
+        I1 = mean( log2(1 + exp(-llr_b(idx1))) );
+        MI(b) = 1 - 0.5*(I0 + I1);
+    end
+
+    % Per-symbol GMI/NGMI
+    GMI  = sum(MI);
     NGMI = GMI / m;
+end
 
+% ---------- helpers ----------
+function s = logsumexp(a)
+    if isempty(a), s = -inf; return; end
+    amax = max(a(:));
+    s = amax + log(sum(exp(a(:) - amax)));
 end
 
 function [data_,reference_] = trimseq(data,reference,skipstart,skip_end)
-
     data_  = data(skipstart+1:end-skip_end,:);
-
-        delta_bits = length(reference) - length(data);
-
-        skip_end = delta_bits + skip_end;
-
-        reference_ = reference(skipstart+1:end-skip_end,:);
-
+    delta  = numel(reference) - numel(data_);
+    reference_ = reference(skipstart+1:end-(skip_end+delta),:);
 end
-
 end

projects/IMDD_base_system/simulation_bwl.m

@@ -1,6 +1,6 @@
 %%% Run parameters
 % TX
-M               = 4;
+M               = 6;
 fsym            = 112e9;
 
 apply_pulsef    = 1;
@@ -44,9 +44,9 @@ dfe_ = sum(dfe_order)>0;
 
 doub_mode = db_mode.no_db;
 
-rop = [-8];
+rop = [-5];
 bwl = [0.5:0.1:1.5];
-fsym = [128:16:170].*1e9;
+fsym = [72:8:170].*1e9;
 ber_vnle = [];
 ber_mlse = [];
 ber_viterbi = [];
@@ -56,7 +56,7 @@ gmi_vnle = [];
 gmi_mlse = []; 
 gmi_mlse_db = [];
 
-for r = 1:length(fsym)
+parfor r = 1:length(fsym)
 
     Pform       = Pulseformer("fsym",fsym(r),"fdac",4*fsym(r),"pulse","rc","pulselength",16,"alpha",rcalpha);
 
@@ -74,72 +74,56 @@ for r = 1:length(fsym)
         "db_precode",db_precode,"db_encode",db_encode,...
         "mrds_code",0,"mrds_blocklength",512,"duobinary_mode",duob_mode).process();
     
-    
-    % Digi_sig.spectrum("displayname",'Digi Spectrum','fignum',10,'normalizeTo0dB',0);
-  
-    
-    %%%%% AWG
-    % El_sig = M8199A("kover",kover).process(Digi_sig);
-    El_sig = AWG("fdac",fdac,"f_cutoff",fsym(r),"lpf_active",0,"kover",kover,"bit_resolution",12,"upsampling_method","samplehold","precomp_sinc_rolloff",1).process(Digi_sig);
-    % El_sig.spectrum("displayname",'Digi Spectrum','fignum',100,'normalizeTo0dB',0);
-    % El_sig = El_sig.setPower(0,"dBm");
+    El_sig = AWG("fdac",fdac,"f_cutoff",fsym(r),"lpf_active",0,"kover",kover,"bit_resolution",12,"upsampling_method","samplehold","precomp_sinc_rolloff",0).process(Digi_sig);
 
     %%%%% Low-pass el. components %%%%%%
-    % f_nyquist       = fsym/2;
-    % tx_bwl = tx_bw_nyquist.*f_nyquist;
     tx_bwl = 50e9;
     El_sig = Filter('filtdegree',4,"f_cutoff",tx_bwl,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(El_sig);
-    % El_sig.spectrum("displayname",'Digi Spectrum','fignum',100,'normalizeTo0dB',1);
 
     %%%%% Electrical Driver Amplifier %%%%%%
-    El_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",3).process(El_sig);
     El_sig = El_sig.normalize("mode","oneone");
 
     %%%%% MODULATE E/O CONVERSION %%%%%%
     [Opt_sig] = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs,"lambda",laser_wavelength,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",random_key+1).process(El_sig);
     
+    %%%%%% Fiber %%%%%%
     Opt_sig    = Fiber("fsimu",Opt_sig.fs,"fiber_length",link_length/1000,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.07).process(Opt_sig);
 
     %%%%%% ROP %%%%%%
-    Rx_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",rop).process(Opt_sig);
+    Opt_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",rop).process(Opt_sig);
     
     %%%%%% PD Square Law %%%%%%
-    Rx_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20,"nep",1.8e-11).process(Rx_sig);
+    PD_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20,"nep",1.8e-11,"randomkey",random_key).process(Opt_sig);
     
     %%%%%% Low-pass RX (PD, El. Connectors and Scope %%%%%%
-    % rx_bwl = rx_bw_nyquist.*f_nyquist;
     rx_bwl = 50e9;
-    Rx_sig = Filter('filtdegree',4,"f_cutoff",rx_bwl,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(Rx_sig);
+    PD_sig = Filter('filtdegree',4,"f_cutoff",rx_bwl,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(PD_sig);
     
     % %%%%%% Low-pass Scope %%%%%%
     Lp_scpe = Filter('filtdegree',4,"f_cutoff",110e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true);
     
-    % Rx_sig.spectrum("displayname",'Analog Rx Spectrum','fignum',100,'normalizeTo0dB',1);
-    % 
     %%%%%% Scope %%%%%%
     Scpe_sig = Scope("fsimu",fdac*kover,"fadc",fadc,...
         "delay",0,"fixed_delay",0,"filtertype",filtertypes.butterworth,...
         "samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,...
-        "adcresolution",8,"quantbuffer",0.1,'block_dc',1,'lpf_active',0,'H_lpf',Lp_scpe).process(Rx_sig);
+        "adcresolution",8,"quantbuffer",0.1,'block_dc',1,'lpf_active',1,'H_lpf',Lp_scpe).process(PD_sig);
     
-    Scpe_sig_resampled = Scpe_sig.resample("fs_out",2*fsym(r));
+    Scpe_sig_2sps = Scpe_sig.resample("fs_out",2*fsym(r));
     % Scpe_sig_resampled.signal = Scpe_sig_resampled.signal(1:2*length(Symbols));
     
-    [~, Scpe_cell, ~, found_sync] = Scpe_sig_resampled.tsynch("reference", Symbols, "fs_ref", fsym(r), "debug_plots", 1);
-    % Scpe_cell = cell(1);
-    % Scpe_cell{1} = Scpe_sig_resampled;
+    [~, Scpe_cell, ~, found_sync] = Scpe_sig_2sps.tsynch("reference", Symbols, "fs_ref", fsym(r), "debug_plots", 0);
+    Rx_sig = Scpe_cell{1};
 
-    Scpe_cell{1}.spectrum("displayname",'Analog Rx Spectrum','fignum',100,'normalizeTo0dB',1);
+    if 1
+        %Duobinary Targeting
 
         eq_    = EQ("Ne",[vnle_order1,vnle_order2,vnle_order3],"Nb",[dfe_order],"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1);
         mlse_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",M,"trellis_states",PAMmapper(M,0).levels);
 
 
-    if 0
-    %Duobinary Targeting
         db_ref_sequence = Duobinary().encode(Symbols);
         db_ref_constellation = unique(db_ref_sequence.signal);
-    [eq_signal, eq_noise] = eq_.process(Scpe_cell{1},db_ref_sequence);
+        [eq_signal, eq_noise] = eq_.process(Rx_sig,db_ref_sequence);
     
         mlse_.DIR = [1,1];
         [mlse_sig_sd,LLR,gmi_mlse_db(r)] = mlse_.process(eq_signal,Symbols);
@@ -163,7 +147,8 @@ for r = 1:length(fsym)
 
     % FFE or VNLE
     eq_    = EQ("Ne",[vnle_order1,vnle_order2,vnle_order3],"Nb",[dfe_order],"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1);
-    [eq_signal_sd, eq_noise] = eq_.process(Scpe_cell{1}, Symbols);
+    [eq_signal_sd, eq_noise] = eq_.process(Rx_sig, Symbols);
+
     [gmi_gomez(r)]      = calc_air(eq_signal_sd, Symbols, "skip_front", 100, "skip_end", 100);
     [gmi_vnle(r)]       = calc_ngmi(eq_signal_sd,Symbols);
     [gmi_bitwise(r)]    = calc_gmi_bitwise(eq_signal_sd,Symbols);
@@ -171,7 +156,6 @@ for r = 1:length(fsym)
     snr_vnle(r) = calc_snr(Symbols, eq_signal_sd-Symbols);
     eq_signal_sd.plot("displayname",'bla','fignum',118);
 
-    
     % Hard decision on VNLE output
     eq_signal_hd = PAMmapper(M, 0).quantize(eq_signal_sd);
     rx_bits = PAMmapper(M,0,"eth_style",0).demap(eq_signal_hd);
@@ -207,6 +191,7 @@ for r = 1:length(fsym)
         rx_bits = PAMmapper(M,0,"eth_style",0).demap(mlse_sig_hd);
         [~,~,ber_viterbi(r),~] = calc_ber(rx_bits.signal,Tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
         fprintf('Viterbi BER: %.2e \n',ber_viterbi(r));
+
     end
 
     
@@ -224,19 +209,19 @@ ylabel('alpha');
 figure(15);hold on
 plot(fsym.*1e-9,gmi_gomez,'DisplayName','gomez','Marker','x','LineStyle','-','Color',cols(1+d,:));
 plot(fsym.*1e-9,gmi_vnle,'DisplayName','vnle other','Marker','x','LineStyle','-','Color',cols(3+d,:));
-plot(fsym.*1e-9,gmi_bitwise,'DisplayName','bitwise','Marker','x','LineStyle','--','Color',cols(5+d,:));
+% plot(fsym.*1e-9,gmi_bitwise,'DisplayName','bitwise','Marker','x','LineStyle','--','Color',cols(5+d,:));
 plot(fsym.*1e-9,gmi_mlse,'DisplayName','MLSE','Marker','*');
-ylim([0 2.6]);
+plot(fsym.*1e-9,gmi_mlse_db,'DisplayName','DB Output','Marker','*');
+ylim([log2(M)-1 log2(M)]);
 xlabel('Baudrate in GBd');
 ylabel('GMI');
 
-
 figure(13);hold on
 plot(fsym.*1e-9,ber_vnle,'DisplayName','VNLE','Marker','x','LineStyle','-','Color',cols(1+d,:));
 plot(fsym.*1e-9,ber_mlse,'DisplayName','MLSE','Marker','x','LineStyle','-','Color',cols(3+d,:));
 plot(fsym.*1e-9,ber_viterbi,'DisplayName','Viterbi','Marker','x','LineStyle','--','Color',cols(5+d,:));
-% plot(bwl,ber_db_diff_precoded,'DisplayName','MLSE db diff','Marker','.','MarkerSize',15,'LineStyle','-');
-% plot(bwl,ber_db,'DisplayName','MLSE db','Marker','.','MarkerSize',15,'LineStyle','-');
+plot(fsym.*1e-9,ber_db_diff_precoded,'DisplayName','MLSE db diff','Marker','.','MarkerSize',15,'LineStyle','-');
+plot(fsym.*1e-9,ber_db,'DisplayName','MLSE db','Marker','.','MarkerSize',15,'LineStyle','-');
 xlabel('Baudrate in GBd');
 ylabel('BER');
 set(gca, 'yscale', 'log');
@@ -244,7 +229,6 @@ set(gca, 'yscale', 'log');
 legend
 
 
-
 % Auxiliary nested helper for numerically stable log-sum-exp
 function s = logsumexp(a)
     % LOGSUMEXP   Compute log(sum(exp(a))) in a numerically stable way