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,signalclass_sd] = process(obj,signalclass,ref_symbolclass) data_in = signalclass.signal; data_ref = ref_symbolclass.signal; [data_out_hd,data_out_sd] = obj.process_(data_in,data_ref); signalclass_hd = signalclass; signalclass_hd.signal = data_out_hd; signalclass_sd = signalclass; signalclass_sd.signal = data_out_sd; end function [VITERBI_ESTIMATION_SYMBOLS,soft_decisions] = 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; tx_bits = PAMmapper(numel(constellation),0).demap(data_ref); % impulse respnse i.e. [0.5, 1.0000] obj.DIR = flip(obj.DIR); % RMS normalization of input data data_in = data_in ./ rms(data_in); % 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{:}).'); % Save first and last symbol of each state first_sym = combs(:,1); last_sym = combs(:,end); states = sum(combs,2); % 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(length(states),length(states)); count_row = 1; count_col = 1; for l1 = 1:length(states) for l2 = 1:length(states) 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 % 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 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 %%%%% % Initialize the output vector pm = zeros(length(states),length(states)); bm_fw = zeros(length(states),length(states),length(data_in)); % Forward Recursion (FSM Computation) for n = 1:length(data_in) bm = abs(data_in(n) - noise_free_received).^2; pm = pm + bm; [pm_survivor_fw(:,n),pm_survivor_fw_idx(:,n)] = min(pm,[],2); % choose lowest path metric as new state pm = repmat(pm_survivor_fw(:,n).',length(states),1); % update pm (chosen state to 2nd dimension -> FROM state) bm_fw(:,:,n) = bm; end % we can now get the best path as min best_fw_path = NaN(1,length(data_in)+1); % find ideal trellis path by going through the trellis backwards [~,best_fw_path(length(data_in)+1)] = min(pm_survivor_fw(:,length(data_in))); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%% BACKWARD PASS %%%%% % Initialize the output vector pm = zeros(length(states),length(states)); pm_survivor_bw = zeros(length(states),length(data_in)+1); pm_survivor_bw_idx = zeros(length(states),length(data_in)+1); bm_bw = zeros(length(states),length(states),length(data_in)+1); % 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 best_fw_path(h) = pm_survivor_fw_idx(best_fw_path(h+1),h); bm = abs(data_in(h) - noise_free_received).^2; pm = pm + bm.'; [pm_survivor_bw(:,h),pm_survivor_bw_idx(:,h)] = min(pm,[],2); % choose lowest path metric as new state pm = repmat(pm_survivor_bw(:,h).',length(states),1); % update pm (chosen state to 2nd dimension -> FROM state) bm_bw(:,:,h) = bm; end VITERBI_ESTIMATION_IDX(1:length(data_in)) = first_sym(best_fw_path(2:end)); VITERBI_ESTIMATION_SYMBOLS(1:length(data_in)) = constellation(best_fw_path(2:end)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%% FORWARD PASS %%%%% %calc the log probabilities (llp's) for k = 2:length(data_in) fsm = repmat(pm_survivor_fw(:,k-1)',[length(states),1]); %pm_survivor_fw size: length(states)xlength(sequence) bm = bm_fw(:,:,k); %bm_fw size: length(states)xlength(states)xlength(sequence) bsm = repmat(pm_survivor_bw(:,k)',[length(states),1]); %pm_survivor_bw size: length(states)xlength(sequence)+1 llp(:,k) = min(fsm+bm+bsm,[],2); end % directly decide based on lowest LLP index [~,llp_based_state_seq]=min(llp); LLP_EST(1:length(data_in)) = constellation(llp_based_state_seq); rx_bits = PAMmapper(numel(constellation),0).demap(LLP_EST'); [~,~,ber_llp,~] = calc_ber(rx_bits,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); fprintf('LLP BER : %.2e \n',ber_llp); % [~,~,ber_llp,~] = calc_ber(circshift(rx_bits,1),tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); % fprintf('LLP BER +1: %.2e \n',ber_llp); % [~,~,ber_llp,~] = calc_ber(circshift(rx_bits,-1),tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); % fprintf('LLP BER -1: %.2e \n',ber_llp); %%%% DECIDE based on Viterbi traceback rx_bits = PAMmapper(numel(constellation),0).demap(VITERBI_ESTIMATION_SYMBOLS'); [~,~,ber_viterbi,~] = calc_ber(rx_bits,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); fprintf('Viterbi BER: %.2e \n',ber_viterbi); % [~,~,ber_viterbi,~] = calc_ber(circshift(rx_bits,1),tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); % fprintf('Viterbi BER: %.2e \n',ber_viterbi); % directly decide based on the FW path metrics [~,fw_direct_state_seq]=min(pm_survivor_fw); FW_EST(1:length(data_in)) = constellation(fw_direct_state_seq); rx_bits = PAMmapper(numel(constellation),0).demap(FW_EST'); [~,~,ber_fw,~] = calc_ber(rx_bits,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); fprintf('FW BER: %.2e \n',ber_fw); % directly decide based on the BW path metrics [~,bw_direct_state_seq]=min(pm_survivor_bw); BW_EST(1:length(data_in)) = constellation(bw_direct_state_seq(2:end)); rx_bits = PAMmapper(numel(constellation),0).demap(BW_EST'); [~,~,ber_bw,~] = calc_ber(rx_bits,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); fprintf('BW BER: %.2e \n',ber_bw); % [~,~,ber_viterbi,~] = calc_ber(circshift(rx_bits,1),tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); % fprintf('BW BER: %.2e \n',ber_viterbi); % [~,~,ber_viterbi,~] = calc_ber(circshift(rx_bits,-1),tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); % fprintf('BW BER: %.2e \n',ber_viterbi); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% tx_symbolpos = zeros(numel(constellation),length(data_ref)); est_symbolpos = zeros(numel(constellation),length(data_ref)); for lvl = 1:numel(constellation) tx_symbolpos(lvl,data_ref==constellation(lvl)) = 1; est_symbolpos(lvl,VITERBI_ESTIMATION_IDX==first_sym(lvl)) = 1; est_symbolpos_llm(lvl,llp_based_state_seq==lvl) = 1; end est_symbolpos= logical(est_symbolpos); tx_symbolpos = logical(tx_symbolpos); est_symbolpos_llm = logical(est_symbolpos_llm); figure(299);clf;hold on; llp_filt = NaN(length(states),length(data_in)); llp_ = llp - min(llp,[],1); for c = 2%:numel(constellation) for c2 = 1:3%numel(constellation) idx = est_symbolpos_llm(c,:); llp_filt(c2,idx) = (llp(c2,idx)); plotidx = 2:200000; scatter(plotidx,llp_filt(c2,plotidx),1,'o','LineWidth',1,'DisplayName',sprintf('LLP of Symbol %d | %d was transmitted',c2,c)); end end % Assume llp is a 4 x N matrix (4 PAM-4 symbols, N time steps) [numSymbols, numTimeSteps] = size(llp); P_symbols = zeros(numSymbols, numTimeSteps); % To hold the soft output probabilities for k = 1:numTimeSteps % Compute the exponentials. (Use -llp_shifted if llp's are costs.) exponents = -llp_(:, k); % Normalize to form a probability vector. P_symbols(:, k) = exponents / sum(exponents); end % Optionally, if you prefer a single soft output value per symbol (an expectation), % define your PAM-4 constellation levels, e.g.: soft_output = sum(P_symbols .* constellation, 1); % 1 x N vector of soft outputs %%%% BITWISE LLR's ??? %%%%% figure(100); clf hold on title('LLR between inner and outer bits') bitmapping = PAMmapper(numel(constellation),0).demap(constellation); for bp = 1:size(bitmapping,2) pos_bitone = find(bitmapping(:,bp)==1); pos_bitzero = find(bitmapping(:,bp)~=1); llr_bits(bp,:) = min(llp(pos_bitone,:),[],1) - min(llp(pos_bitzero,:),[],1); subplot(size(bitmapping,2),1,bp) scatter(1:length(data_ref),llr_bits(bp,:),1,'.'); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%% % From: Log-Likelihood Probabilities LLP --> To: Log-Likelihood Ratios LLR for s = 1:length(states)-1 llr(s,:) = llp_(s,:) - llp_(s+1,:); % subtract llp's of successive const. points to get llr's end % Build Sets of LLR's that contain only those values at % timepoint k, where the symbols: % a) were actually transmitted: set "S" % b) were decoded by viterbi: set "SC" S = NaN(length(states)-1,length(data_in)); SC = NaN(length(states)-1,length(data_in)); llp_inf = llp_; llp_inf(llp_==0) = Inf; [~,scnd_idx] = min(llp_inf,[],1) ; for s = 1:length(states)-1 %%% SET "S" % find all time indices k, where const point s was actually transmitted indice = tx_symbolpos(s,:)==1; S(s,indice) = llr(s,indice); %calc mean of all llr's where symbol s was transmitted K(s,1) = mean(S(s,indice),'omitnan'); % find all time indices k, where const. point s+1 was actually transmitted indice_plusone = tx_symbolpos(s+1,:)==1; S(s,indice_plusone) = llr(s,indice_plusone); K(s,2) = mean(S(s,indice_plusone),'omitnan'); %%% SET "SC" % find all time indices k, where symbol s was decoded % idx: find positions in time where a symbol s was decoded idx = est_symbolpos_llm(s,:); % idx 2: find positions where the strongest competitor is % s+1, i.e. the second best llp is at s+1 idx2 = scnd_idx == s+1; SC(s,idx&idx2) = llr(s,idx&idx2); KC(s,1) = mean(SC(s,idx&idx2),'omitnan'); % find all time indices k, where symbol s+1 was decoded % idx: find positions in time where a symbol s+1 was decoded idx = est_symbolpos_llm(s+1,:); idx2 = scnd_idx == s; % idx 2: find positions where the strongest competitor is % s, i.e. the second best llp is at s SC(s,idx&idx2) = llr(s,idx&idx2); KC(s,2) = mean(SC(s,idx&idx2),'omitnan'); % scale the sets, using the average (?) of the llrcn(s,:) = SC(s,:) * ( constellation(s+1)-constellation(s) ) ./ ( K(s,2)-K(s,1)) + decisionLevels(s); end figure(111) title("Bla") clf for s = 1:length(states)-1 figure(1111) hold on title(sprintf('llrcn = llp %d - llp %d',s, s+1,s, s+1)); scatter(1:length(llrcn),llrcn(s,:),1,'.'); % STUFENARTIGE LLR'S UND LLP'S %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure(111) subplot(1,length(states),s) hold on title(sprintf('llr = llp %d - llp %d \n SC=llr(tx sym = %d or %d)',s, s+1,s, s+1)); scatter(1:length(llr),llr(s,:),1,'.'); scatter(1:length(SC),SC(s,:),1,'.'); yline([K(s,1),K(s,2)]); low = min(min(llr)); hi = max(max(llr)); ylim([low,hi]); subplot(1,length(states),length(states)) hold on scatter(1:length(llr),llr(s,:),1,'.'); ylim([low,hi]); % HISTOGRAMME %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure(112) subplot(length(states),1,s) hold on title(sprintf('histogram of llr; between const point %d - %d',s, s+1)); histogram(llrcn(s,:),10000,'EdgeAlpha',0) % histogram(SC(s,:),10000,'EdgeAlpha',0) xlim([-20, 20]); subplot(length(states),1,length(states)) hold on histogram(llrcn(s,:),10000,'EdgeAlpha',0) % histogram(SC(s,:),10000,'EdgeAlpha',0) % xlim([-20, 20]); end softdecisions = mean(llrcn,1,'omitnan'); distance = abs(VITERBI_ESTIMATION_SYMBOLS-llrcn); [win_cost,win_idx] = min(distance,[],1); for i = 1:length(data_in) soft_decisions(i) = llrcn(win_idx(i),i); end soft_decisions = max(llrcn,[],1); showLevelHistogram(soft_decisions,data_ref) figure() % scatter(1:length(data_in),VITERBI_ESTIMATION_SYMBOLS,1,'.'); scatter(1:length(data_in),soft_decisions,1,'.'); MLM_ESTIMATION(1:length(data_in)) = PAMmapper(numel(constellation),0).quantize(soft_decisions)'; rx_bits = PAMmapper(numel(constellation),0).demap(MLM_ESTIMATION'); [~,~,ber_mlm,~] = calc_ber(rx_bits,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); fprintf('MLM BER: %.2e \n',ber_mlm); end end methods (Access=private) end end