M = 6;

bitpattern = [];
s = RandStream('twister','Seed',1);
for i = 1:log2(M)
    N = 2^(17-1); %length of prbs
    bitpattern(:,i) = randi(s,[0 1], N, 1);
end

if M == 6
    bitpattern = reshape(bitpattern',[],1);
    bitpattern = bitpattern(1:end-mod(length(bitpattern),5));
end

bits = Informationsignal(bitpattern);

symbols = PAMmapper(M,0).map(bits);

bits_rx = PAMmapper(M,0).demap(symbols);
[~,~,ber_direct,~] = calc_ber(bits.signal,bits_rx.signal,"skip_front",0,"skip_end",0,"returnErrorLocation",1);
assert(ber_direct==0,'Mapping is wrong');

nBursts = 0;
% No Precoding %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%SEND DIRECTLY
symbols_tx = symbols;

symbols_rx = introduce_symbol_errors(symbols_tx, 1, 10, nBursts, 42);

%RECEIVE BRANCH (do nothing special)
bits_rx = PAMmapper(M,0).demap(symbols_rx);

[~,~,ber,errpos] = calc_ber(bits.signal,bits_rx.signal,"skip_front",0,"skip_end",0,"returnErrorLocation",1);
disp(['BER normal: - ',sprintf('%.1E',ber),' - - PAM-',num2str(M)]);
bursts_normal = count_error_bursts(errpos, 20);


% Precode Emulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%SEND DIRECTLY
symbols_tx = symbols;

symbols_rx = introduce_symbol_errors(symbols_tx, 1, 10, nBursts, 42);

%REFERENCE BRACH
symbols_db = Duobinary().encode(symbols_tx);
symbols_tx_emu = Duobinary().decode(symbols_db);
bits_tx_emu = PAMmapper(M,0).demap(symbols_tx_emu);

% symbols_rx = introduce_symbol_errors(symbols_tx, 1, 10, 200, 42);

%RECEIVE BRANCH
symbols_db = Duobinary().encode(symbols_rx);
symbols_rx_emu = Duobinary().decode(symbols_db);
bits_rx = PAMmapper(M,0).demap(symbols_rx_emu);

[~,~,ber_precode_emulation,errpos_precode_emulation] = calc_ber(bits_tx_emu.signal,bits_rx.signal,"skip_front",0,"skip_end",0,"returnErrorLocation",1);
disp(['BER precode emulation: ',sprintf('%.1E',ber_precode_emulation),' - - PAM-',num2str(M)]);
bursts_precode_emulation = count_error_bursts(errpos_precode_emulation, 20);




% Precode at Tx   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%SEND PRECODED DATA
symbols_tx_prec = Duobinary().precode(symbols);

symbols_rx_prec = introduce_symbol_errors(symbols_tx_prec, 1, 10, nBursts, 42);

%RECEIVE BRANCH
symbols_db = Duobinary().encode(symbols_rx_prec);
symbols_rx_prec = Duobinary().decode(symbols_db);
bits_rx = PAMmapper(M,0).demap(symbols_rx_prec);

[~,~,ber_precoded,errpos_precoded] = calc_ber(bits.signal,bits_rx.signal,"skip_front",0,"skip_end",0,"returnErrorLocation",1);
disp(['BER precoded: ',sprintf('%.1E',ber_precoded),' - - PAM-',num2str(M)]);
burst_precoded = count_error_bursts(errpos_precoded, 20);


% Precode at Tx but omit at Rx  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%SEND PRECODED DATA
symbols_tx_prec = Duobinary().precode(symbols);
bits_tx_prec = PAMmapper(M,0).demap(symbols_tx_prec);

symbols_rx_omit = introduce_symbol_errors(symbols_tx_prec, 1, 10, nBursts, 42);

%RECEIVE BRANCH
bits_rx = PAMmapper(M,0).demap(symbols_rx_omit);

[~,~,ber_omit,errpos_omit] = calc_ber(bits_tx_prec.signal,bits_rx.signal,"skip_front",0,"skip_end",0,"returnErrorLocation",1);
disp(['BER (omit precode): ',sprintf('%.1E',ber_omit),' - - PAM-',num2str(M)]);
burst_omit = count_error_bursts(errpos_omit, 20);

if 0
    cols = linspecer(8);
    figure();hold on;
    stem(1:20,bursts_normal,'LineWidth',2,'Color',cols(4,:),'Marker','_','DisplayName','w/o diff. precoder');
    stem(1:20,bursts_precode_emulation,'LineWidth',2,'Color',cols(3,:),'Marker','.','LineStyle','-','DisplayName','emulated precoder');
    stem(1:20,burst_precoded,'LineWidth',1,'Color',cols(6,:),'Marker','_','DisplayName','w/ diff. precoder');
    stem(1:20,burst_omit,'LineWidth',1,'Color',cols(5,:),'Marker','.','LineStyle',':','DisplayName','omit precoder');
    xlabel('Bit Error Burst Length')
    ylabel('Occurence')
    set(gca, 'yscale', 'log');
end


%% State Analysis
signal_to_analyze = symbols_tx_emu;
x = signal_to_analyze.signal(:);
levels = sort(unique(x)).';              % or provide known 1x6 level values

[~,ix] = min(abs(x - levels),[],2);
x = levels(ix);                          % snapped/quantized

%% TRANSITION COUNTS & PROBABILITIES
K = numel(levels);
% map to state indices 1..K
[tf, idx] = ismember(x, levels);        
idx = idx(:);
from = idx(1:end-1);
to   = idx(2:end);
from = idx(1:2:end);
to   = idx(2:2:end);

% counts C(from,to)
C = accumarray([from,to], 1, [K K], @sum, 0);
% row-stochastic transition matrix P(to|from)
rowSums = sum(C,2);
P = C ./ max(rowSums,1);

%% 1) HEATMAP (which transitions are more probable?)
figure('Name','Transition Probabilities (to | from)');
h = heatmap(levels, levels, P, 'Colormap', parula, 'ColorbarVisible','on');
colormap(gca,[[1,1,1];flip(cbrewer2('Spectral',100))]);clim([0,ceil(max(P(:))*10)/10]);
h.XLabel = 'From state (level)'; 
h.YLabel = 'To state (level)';
h.Title  = 'P(to | from)';

%% 2) WEIGHTED TRANSITION GRAPH
% Use dtmc if you have Econometrics Toolbox:
mc = dtmc(P, 'StateNames', string(levels.*PAMmapper(M,0).get_scaling));
figure('Name','Markov Graph (dtmc)');
gp = graphplot(mc, 'ColorEdges',true, 'LabelEdges',true);


function symbols = introduce_symbol_errors(symbols, j, maxBurstLen, nBursts, seed)
%INTRODUCE_SYMBOL_ERRORS injects bursty level errors into symbols.signal.
% symbols.signal : column/row vector of quantized levels (exactly one of 6 values)
% j              : max level step per sample (default 1)
% maxBurstLen    : maximum burst length (default 8)
% nBursts        : number of bursts to insert (default ~1% of length)
% seed           : RNG seed (optional)

    if nargin < 2 || isempty(j),            j = 1; end
    if nargin < 3 || isempty(maxBurstLen),  maxBurstLen = 8; end
    x = symbols.signal(:);
    N = numel(x);
    if nargin < 4 || isempty(nBursts),      nBursts = max(1, round(0.01*N)); end
    if nargin >= 5 && ~isempty(seed),       rng(seed); end

    % known levels and index mapping
    lvls = sort(unique(x)).';
    K = numel(lvls);

    [~, idx] = ismember(x, lvls);   % idx in 1..6

    used = false(N,1);              % avoid overlapping bursts
    burst_ranges = zeros(nBursts,2);

    for b = 1:nBursts
        % pick start not inside an existing burst
        s = randi(N);
        while used(s), s = randi(N); end
        L = randi(maxBurstLen);
        e = min(N, s+L-1);

        % mark used range
        used(s:e) = true;
        burst_ranges(b,:) = [s e];

        % choose one direction for the whole burst: -1 (down) or +1 (up)
        dir = randi([0 1])*2 - 1;

        % apply level errors within the burst
        for t = s:e
            k = idx(t);  % current level index (1..6)

            % force inward movement at edges; prevents "flipping" to opposite edge
            if k == 1 && dir == -1, dir = +1; end
            if k == K && dir == +1, dir = -1; end

            step = randi([1 j]);              % 1..j steps
            kNew = k + dir*step;

            % clamp to [1,K], no wrap-around
            if kNew < 1, kNew = 1; elseif kNew > K, kNew = K; end

            % if clamped to the same edge repeatedly, flip direction to keep changing
            if kNew == k
                dir = -dir;
                kNew = max(1, min(K, k + dir*step));
            end

            idx(t) = kNew;
        end
    end

    x_err = lvls(idx);
    symbols.signal = reshape(x_err, size(symbols.signal));  % preserve original shape

end