classdef MLSE < handle %MLSE calculates the most probable sequence for an input signal with given/ known channel impulse response of any length properties(Access=public) M %PAM-M DIR trellis_states duobinary_output end methods (Access=public) function obj = MLSE(options) %NAME Construct an instance of this class % Detailed explanation goes here arguments options.M double = 4; options.DIR double = [1]; options.trellis_states double = [-3 -1 1 3]; options.duobinary_output logical = false; end % fn = fieldnames(options); for n = 1:numel(fn) try obj.(fn{n}) = options.(fn{n}); end end end function [signalclass_hd,LLR,GMI] = process(obj,signalclass,ref_symbolclass) data_in = signalclass.signal; shape_in = size(data_in); data_ref = ref_symbolclass.signal; [data_out_hd,LLR,GMI] = obj.process_(data_in,data_ref); try data_out_hd = reshape(data_out_hd,shape_in(1),shape_in(2)); catch warning('output reshaping failed after MLSE'); end signalclass_hd = signalclass; signalclass_hd.signal = data_out_hd; % signalclass_sd = signalclass; % signalclass_sd.signal = data_out_sd; end function [VITERBI_ESTIMATION_SYMBOLS,LLR_maxlogmap,GMI] = process_(obj,data_in,data_ref) debug = 1; trellis_state_mode = 0; % 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 = 0; % 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 DIR_nonzero = find(obj.DIR ~= 0); if DIR_nonzero(1) > 1 obj.DIR(1:DIR_nonzero(1)-1) = []; end if isscalar(obj.DIR) obj.DIR = [0 obj.DIR]; end % impulse respnse to remove from signal obj.DIR = flip(obj.DIR); %i.e. -0.2676 -0.0478 1.0000 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 pre_comb_mat = repmat(obj.trellis_states,length(obj.DIR)-1,1); pre_comb_cell = mat2cell(pre_comb_mat,ones(1,size(pre_comb_mat,1)),size(pre_comb_mat,2)); combs = fliplr(combvec(pre_comb_cell{:}).'); first_sym = combs(:,1); % das ist das älteste/ trailing Symbol aus der sequenz last_sym = combs(:,end); %hiermit wird entschieden/ das ist das cursor symbol am ende der sequenz 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 % 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 end end h = flip(obj.DIR(:)).'; data_in = data_in(:); y_ideal = conv(data_ref(:), h, "same"); 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 % 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) %%%%% % Initialize the output vector pm = zeros(nStates,nStates); bm_fw = zeros(nStates,nStates,length(data_in)); % first start is evaluated without ISI/ wihout the full Impulse response % so simply use the constellation here bm = -(data_in(1) - last_sym).^2 * inv2s2; pm = pm + bm; [alpha(:,1),pm_survivor_fw_idx(:,1)] = max(pm,[],2); pm = repmat(alpha(:,1).',nStates,1); bm_fw(:,:,1) = pm; % Forward Recursion (FSM Computation) 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 (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; end % we can now get the best path as min viterbi_path = NaN(1,length(data_in)); % find ideal trellis path by going through the trellis backwards [~,viterbi_path(length(data_in))] = max(alpha(:,length(data_in))); for n = length(data_in):-1:2 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); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%% BACKWARD PASS (Beta's) %%%%% % Initialize the output vector pm = zeros(nStates,nStates); beta = zeros(nStates,length(data_in)); pm_survivor_bw_idx = zeros(nStates,length(data_in)); bm_bw = zeros(nStates,nStates,length(data_in)); % starting with the state that has the lowest sum path % metric, follow the stored information about the % predecessor 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) bm_bw(:,:,h) = bm; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%% FORWARD PASS PAM 2,4,8 (Combine Alpha and Beta to yield LLP's) %%%%% %calc the log probabilities (llp's) for k = 1:length(data_in) if k == 1 alpha_ = repmat(alpha(:,k)',[nStates,1])'; beta_ = beta(:,k); LLP(:,k) = max(alpha_ + beta_,[],2); else alpha_ = repmat(alpha(:,k-1)',[nStates,1])'; gamma_ = bm_fw(:,:,k)'; beta_ = beta(:,k); LLP(:,k) = max(alpha_ + gamma_,[],1) + beta_'; end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%% 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) % compute symbol‐posteriors from LLP in the log‐domain: amax = max(LLP,[],1); logZ = amax + log(sum(exp(LLP - amax), 1)); logPstate = LLP - logZ; % still in log‐domain state_prob = exp(logPstate); % exact, sums to 1 if obj.M == 6 num_bits = 5; % all possible transitions (for now 36, including the "edges" % of the QAM 32 constellation) states = PAMmapper(6,0,"eth_style",0).levels; pam6transitions = combvec(states,states)'; % pam6transitions = % [-5 -5; % -3 -5; % -1 -5; ... pam6bits = PAMmapper(6,0,"eth_style",0).demap(reshape(pam6transitions',[],1)./sqrt(10)); pam6bits = reshape(pam6bits',5,[])'; [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_sym_1, idx_sym_2]; tx_bits_pam6_reshaped = reshape(tx_bits,5,[])'; % N x 5 numPairs = floor(size(LLP,2)/2); LLR_exact = zeros(numPairs,5); LLR_maxlogmap = zeros(numPairs,5); for k = 1:numPairs symbol1 = 2*k-1; symbol2 = 2*k; LLP1 = LLP (:,symbol1); LLP2 = LLP (:,symbol2); prob1 = state_prob(:,symbol1); prob2 = state_prob(:,symbol2); % 36 joint‐metrics M = log P(i)*P(j) = L1(i)+L2(j) Mij = LLP1(pam6ind(:,1)) + LLP2(pam6ind(:,2)); pij = prob1(pam6ind(:,1)) .* prob2(pam6ind(:,2)); % now for each of the 5 bits do exact-LLR or max-log for b = 1:num_bits idx_sym_1 = pam6bits(:,b)==1; idx_bit_1 = pam6bits(:,b)==0; %--- exact LLR from probabilities 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_sym_1) ) - max( Mij(idx_bit_1) ); % N x num_bits end end MI = zeros(1, num_bits); for k = 1:num_bits idx_bit_1 = (tx_bits_pam6_reshaped(:,k) == 0); %wo sind die 1en 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_sym_1,k); % Calculate mutual information for bit position k I0 = mean(log2(1 + exp(llr0))); % exp(--LLR) = exp(positive) > 1 I1 = mean(log2(1 + exp(-llr1))); % exp(-+LLR) = exp(negative) < 1 MI(k) = 1 - 0.5 * (I0 + I1); end GMI = sum(MI); % Total mutual information per symbol GMI = GMI/2; else % Number of symbols and bits per symbol num_bits = log2(length(obj.trellis_states)); % 2 bits per symbol % bit_mapping = PAMmapper(length(obj.trellis_states),0,"eth_style",0).demap(first_sym./rms(first_sym)); bit_mapping = PAMmapper(length(obj.trellis_states),0,"eth_style",0).showBitMapping; % Initialize LLR storage LLR_maxlogmap = zeros(length(data_in),num_bits); LLR_exact = zeros(length(data_in),num_bits); % Compute bit-wise LLRs 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_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); % Sum probabilities over states for which the bit is 1 and 0, respectively. 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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%% CALC NGMI %%%%% MI = zeros(1, num_bits); LLR_exact = LLR_exact; for k = 1:num_bits 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_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 I1 = mean(log2(1 + exp(-llr1))); % exp(-+LLR) = exp(negative) < 1 MI(k) = 1 - 0.5 * (I0 + I1); end GMI = sum(MI); % Total mutual information for 2 symbols end VITERBI_ESTIMATION_SYMBOLS = VITERBI_ESTIMATION_SYMBOLS./rms(VITERBI_ESTIMATION_SYMBOLS); if debug %%% DEBUG PLOT LIKELIHOOD RATIOS %%% figure(115);clf subplot(2,1,1) for bit = 1:num_bits hold on; histogram(LLR_exact(:,bit),1000,"DisplayName",sprintf('Actual LLR of Bit Pos %d',bit),'LineStyle','none','FaceAlpha',0.4); end legend subplot(2,1,2) for bit = 1:num_bits hold on; histogram(LLR_maxlogmap(:,bit),1000,"DisplayName",sprintf('Max Log LLR of Bit Pos %d',bit),'LineStyle','none','FaceAlpha',0.4); end legend end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%% CHECK BER's %%%%% 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); rx_bits = PAMmapper(obj.M,0,"eth_style",0).demap(VITERBI_ESTIMATION_SYMBOLS'); rx_bits = reshape(rx_bits',[],1); [~,numErr,ber_viterbi,~] = calc_ber(rx_bits,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); 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); [~,~,ber_llr,~] = calc_ber(bit_stream,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); fprintf('LLR MaxLogMAP BER : %.2e \n',ber_llr); % Convert LLR values to a hard-decision bit stream bit_stream = LLR_exact > 0; %ratio separates lower or higher than =0 -> simply decode for the negative value bit_stream = reshape(bit_stream',[],1); [~,~,ber_llr,~] = calc_ber(bit_stream,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); fprintf('LLR BER : %.2e \n',ber_llr); % DECIDE based on lowest LLP index and check BER [~,llp_based_state_seq]=max(LLP); LLP_EST(1:length(data_in)) = first_sym(llp_based_state_seq); LLP_EST = LLP_EST./rms(LLP_EST); rx_bits = PAMmapper(obj.M,0,"eth_style",0).demap(LLP_EST'); rx_bits = reshape(rx_bits',[],1); [~,~,ber_llp,~] = calc_ber(rx_bits,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); fprintf('LLP BER : %.2e \n',ber_llp); % directly decide based on the FW path metrics [~,fw_direct_state_seq]=max(alpha); FW_EST(1:length(data_in)) = first_sym(fw_direct_state_seq); FW_EST = FW_EST./rms(FW_EST); 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); disp('Stop DEBUG MLSE:') fprintf('\n') end end function s = logsumexp(a,dim) % returns log(sum(exp(a),dim)) safely amax = max(a,[],dim); s = amax + log(sum(exp(a - amax), dim)); end end end