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 % do more stuff 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) % 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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%% 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 % 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 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; end count_row = count_row + 1; end count_col = count_col + 1; count_row = 1; end % 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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%% 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 ; 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 ; pm = pm + bm; [alpha(:,n),pm_survivor_fw_idx(:,n)] = max(pm,[],2); % choose lowest path metric as 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 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 ; 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 %%%%% nml_LLP = LLP - max(LLP); expLLP = exp(nml_LLP); %subtract highst value for better numerical stability, LLP's are not always close to zero state_prob = expLLP ./ sum(expLLP); if obj.M == 6 num_bits = 5; % all possible transitions (for now 36, including the "edges" % of the QAM 32 constellation) pam6transitions = combvec(obj.trellis_states,obj.trellis_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,[])'; %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); assert(all(ok1)&all(ok2), 'Some transition amplitude not found in trellis_states') pam6ind = [idx_bit_0, idx2]; tx_bits_pam6_reshaped = reshape(tx_bits,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_bit_0 = pam6bits(:,b)==1; idx_bit_1 = pam6bits(:,b)==0; %--- exact LLR from probabilities P1 = sum(pij(idx_bit_0)); 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 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_bit_0 = (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); % 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)); % 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_1 = (tx_bits(:,k) == 0); %wo sind die 1en idx_bit_0 = (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); % Calculate mutual information for bit position k I0 = mean(log2(1 + exp(llr0(100:end-100)))); % exp(--LLR) = exp(positive) > 1 I1 = mean(log2(1 + exp(-llr1(100:end-100)))); % 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); debug = 0; 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); % 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); % 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); end end end methods (Access=private) end end