From f4a22d23a2dd67ecf000ccce4c28437013a277bd Mon Sep 17 00:00:00 2001 From: "silas (home)" Date: Mon, 11 Aug 2025 16:31:09 +0200 Subject: [PATCH] SD stuff working for PAM6 and MLSE is also calibrated --- Classes/01_transmit/PAMmapper.m | 23 +- Classes/04_DSP/Sequence Detection/MLSE.m | 185 ++++++++------- Functions/Metrics/calc_gmi_bitwise.m | 206 ++++++++++------ Functions/Metrics/calc_ngmi.m | 258 +++++++++++++-------- projects/IMDD_base_system/simulation_bwl.m | 170 ++++++-------- 5 files changed, 490 insertions(+), 352 deletions(-) diff --git a/Classes/01_transmit/PAMmapper.m b/Classes/01_transmit/PAMmapper.m index cf4df17..a44d8cd 100644 --- a/Classes/01_transmit/PAMmapper.m +++ b/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 diff --git a/Classes/04_DSP/Sequence Detection/MLSE.m b/Classes/04_DSP/Sequence Detection/MLSE.m index e9d99dc..3883877 100644 --- a/Classes/04_DSP/Sequence Detection/MLSE.m +++ b/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,21 +117,54 @@ 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 + % 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); + + % 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 - y_clean = conv( data_ref, flip(obj.DIR), "same" ); - sigma2 = var( data_in - y_clean ); - inv2s2 = 1/(2*sigma2); + noise_free_received = (g*noise_free_received + b); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%% FORWARD PASS (VITERBI -Alpha's) %%%%% @@ -199,50 +227,47 @@ classdef MLSE < handle end - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %%%%% FORWARD PASS PAM 2,4,8 (Combine Alpha and Beta to yield LLP's) %%%%% + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %%%%% FORWARD PASS PAM 2,4,8 (Combine Alpha and Beta to yield LLP's) %%%%% - %calc the log probabilities (llp's) + %calc the log probabilities (llp's) - for k = 1:length(data_in) + for k = 1:length(data_in) - if k == 1 + if k == 1 - alpha_ = repmat(alpha(:,k)',[nStates,1])'; - beta_ = beta(:,k); + alpha_ = repmat(alpha(:,k)',[nStates,1])'; + beta_ = beta(:,k); - LLP(:,k) = max(alpha_ + beta_,[],2); + LLP(:,k) = max(alpha_ + beta_,[],2); - else + else - alpha_ = repmat(alpha(:,k-1)',[nStates,1])'; - gamma_ = bm_fw(:,:,k)'; - beta_ = beta(:,k); - - LLP(:,k) = max(alpha_ + gamma_,[],1) + beta_'; - end + alpha_ = repmat(alpha(:,k-1)',[nStates,1])'; + gamma_ = bm_fw(:,:,k)'; + beta_ = beta(:,k); + LLP(:,k) = max(alpha_ + gamma_,[],1) + beta_'; end - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %%%%% FORWARD PASS PAM2,4,8 %%%%% + end - 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) + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %%%%% FORWARD PASS PAM2,4,8 %%%%% - % 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 + 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) - % figure - % hold on; - % for i = 1:obj.M - % scatter(1:length(expLLP),expLLP(i,:),1,'.'); - % end - % scatter(1:length(expLLP),max(expLLP(:,:)),3,'.'); + % 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 @@ -250,30 +275,22 @@ 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; ... 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 - % .... + pam6bits = reshape(pam6bits',5,[])'; - [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); @@ -289,33 +306,33 @@ classdef MLSE < handle 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)); + 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_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 - MI = zeros(1, num_bits); + 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 + 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); + 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 @@ -324,7 +341,7 @@ classdef MLSE < handle end GMI = sum(MI); % Total mutual information per symbol - GMI = GMI/2; + GMI = GMI/2; else @@ -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 @@ -356,16 +373,16 @@ classdef MLSE < handle %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%% CALC NGMI %%%%% - MI = zeros(1, num_bits); + 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 + 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); + 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 @@ -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 @@ -389,7 +406,7 @@ classdef MLSE < handle 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; @@ -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 diff --git a/Functions/Metrics/calc_gmi_bitwise.m b/Functions/Metrics/calc_gmi_bitwise.m index 1ce8847..0b9f0b4 100644 --- a/Functions/Metrics/calc_gmi_bitwise.m +++ b/Functions/Metrics/calc_gmi_bitwise.m @@ -1,94 +1,158 @@ 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 end -options.skip_end = abs(options.skip_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) 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); + + for n = 1:N + y = test_signal(n); + ll_table = -((y - levels).^2)/(2*sigma2) + log(priors); % 1xM (log-domain) + for b = 1:nBits + idx0 = grayBits(:,b)==0; + idx1 = ~idx0; + L0 = logsumexp(ll_table(idx0)); + L1 = logsumexp(ll_table(idx1)); + LLR(n,b) = L1 - L0; % natural-log LLR + end + end + + tx_bits = PAMmapper(M,0,"eth_style",0).demap(reference_signal); + + % Bit-wise MI + MI_bits = zeros(1,nBits); for b = 1:nBits - idx0 = grayBits(:,b)==0; - idx1 = grayBits(:,b)==1; - L0 = logsumexp(ll_table(idx0)); - L1 = logsumexp(ll_table(idx1)); - LLR(n,b) = L1 - L0; + r0 = LLR(tx_bits(:,b)==0,b); + r1 = LLR(tx_bits(:,b)==1,b); + I0 = mean( log2(1 + exp( r0)) ); + I1 = mean( log2(1 + exp(-r1)) ); + MI_bits(b) = 1 - 0.5*(I0 + I1); end + + GMI = sum(MI_bits); % bits / symbol + NGMI = GMI / nBits; end -tx_bits = PAMmapper(M,0,"eth_style",0).demap(reference_signal); - -% Compute GMI (bit‐wise MI averaged over all symbols) -MI_bits = zeros(1,nBits); -for b = 1:nBits - r0 = LLR(tx_bits(:,b)==0,b); - r1 = LLR(tx_bits(:,b)==1,b); - I0 = mean( log2(1 + exp(r0)) ); - I1 = mean( log2(1 + exp( -r1)) ); - MI_bits(b) = 1 - 0.5*(I0 + I1); +% ===== Helpers ===== +function s = logsumexp(a) + if isempty(a), s = -inf; return; end + m = max(a(:)); + s = m + log(sum(exp(a(:) - m))); 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))); - 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,:); - - end - -end \ No newline at end of file +function [data_,reference_] = trimseq(data,reference,skipstart,skip_end) + data_ = data(skipstart+1:end-skip_end,:); + delta = numel(reference) - numel(data_); + reference_ = reference(skipstart+1:end-(skip_end+delta),:); +end +end diff --git a/Functions/Metrics/calc_ngmi.m b/Functions/Metrics/calc_ngmi.m index afbbf54..5d7363a 100644 --- a/Functions/Metrics/calc_ngmi.m +++ b/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 end -options.skip_end = abs(options.skip_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) 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 - error_vector = (test_signal-reference_signal); - sigma2 = var(error_vector); %noise variance + % 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 - %%% Separate Classes - 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)); - end + % Precompute constant term for log-likelihood + logc = -0.5*log(2*pi*sigma2); - N = length(test_signal); % Number of received samples - - 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); - end - end - - % Normalize the noise impact term by N - bit_llr_sum = bit_llr_sum / N; - - % GMI - GMI = H_X + bit_llr_sum; - NGMI = GMI / m; + % 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 - 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,:); + % 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 + 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 -end \ No newline at end of file + % --- 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 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 = numel(reference) - numel(data_); + reference_ = reference(skipstart+1:end-(skip_end+delta),:); +end +end diff --git a/projects/IMDD_base_system/simulation_bwl.m b/projects/IMDD_base_system/simulation_bwl.m index f47b0aa..b040f1e 100644 --- a/projects/IMDD_base_system/simulation_bwl.m +++ b/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,104 +74,88 @@ 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); + 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); - - mlse_.DIR = [1,1]; - [mlse_sig_sd,LLR,gmi_mlse_db(r)] = mlse_.process(eq_signal,Symbols); - - mlse_sig_hd = PAMmapper(M,0,"eth_style",0).quantize(mlse_sig_sd); - mlse_sig_hd_precoded = Duobinary().encode(mlse_sig_hd,"M",M); - mlse_sig_hd_precoded = Duobinary().decode(mlse_sig_hd_precoded,"M",M); - - tx_symbols_precoded = Duobinary().encode(Symbols); - tx_symbols_precoded = Duobinary().decode(tx_symbols_precoded); - - tx_bits_precoded = PAMmapper(M,0,"eth_style",0).demap(tx_symbols_precoded); - - rx_bits_mlse = PAMmapper(M,0,"eth_style",0).demap(mlse_sig_hd_precoded); - [~,errors_db_diff_precoded,ber_db_diff_precoded(r),~] = calc_ber(rx_bits_mlse.signal,tx_bits_precoded.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - - %B) Just determine BER - rx_bits_mlse = PAMmapper(M,0,"eth_style",0).demap(mlse_sig_hd); - [bits_mlse,errors_mlse,ber_db(r),~] = calc_ber(rx_bits_mlse.signal,Tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + db_ref_sequence = Duobinary().encode(Symbols); + db_ref_constellation = unique(db_ref_sequence.signal); + [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); + + mlse_sig_hd = PAMmapper(M,0,"eth_style",0).quantize(mlse_sig_sd); + mlse_sig_hd_precoded = Duobinary().encode(mlse_sig_hd,"M",M); + mlse_sig_hd_precoded = Duobinary().decode(mlse_sig_hd_precoded,"M",M); + + tx_symbols_precoded = Duobinary().encode(Symbols); + tx_symbols_precoded = Duobinary().decode(tx_symbols_precoded); + + tx_bits_precoded = PAMmapper(M,0,"eth_style",0).demap(tx_symbols_precoded); + + rx_bits_mlse = PAMmapper(M,0,"eth_style",0).demap(mlse_sig_hd_precoded); + [~,errors_db_diff_precoded,ber_db_diff_precoded(r),~] = calc_ber(rx_bits_mlse.signal,tx_bits_precoded.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + + %B) Just determine BER + rx_bits_mlse = PAMmapper(M,0,"eth_style",0).demap(mlse_sig_hd); + [bits_mlse,errors_mlse,ber_db(r),~] = calc_ber(rx_bits_mlse.signal,Tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); end % 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); - [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); + [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); 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); @@ -181,32 +165,33 @@ for r = 1:length(fsym) if 1 - % Process through postfilter and MLSE - pf_ncoeffs = 1; - pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); - mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); - [mlse_sig_sd,whitened_noise] = pf_.process(eq_signal_sd, eq_noise); - mlse_.DIR = pf_.coefficients; - alpha(r) = pf_.coefficients(2); - [signalclass_hd,LLR,gmi_mlse(r)] = mlse_.process(mlse_sig_sd,Symbols); - mlse_sig_hd = PAMmapper(M, 0, "eth_style", 0).quantize(signalclass_hd); - rx_bits = PAMmapper(M,0,"eth_style",0).demap(mlse_sig_hd); - [~,~,ber_mlse(r),~] = calc_ber(rx_bits.signal,Tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - fprintf('BER: %.2e \n',ber_mlse(r)); - fprintf('GMI MLSE: %.5f \n',gmi_mlse(r)); + % Process through postfilter and MLSE + pf_ncoeffs = 1; + pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); + mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); + [mlse_sig_sd,whitened_noise] = pf_.process(eq_signal_sd, eq_noise); + mlse_.DIR = pf_.coefficients; + alpha(r) = pf_.coefficients(2); + [signalclass_hd,LLR,gmi_mlse(r)] = mlse_.process(mlse_sig_sd,Symbols); + mlse_sig_hd = PAMmapper(M, 0, "eth_style", 0).quantize(signalclass_hd); + rx_bits = PAMmapper(M,0,"eth_style",0).demap(mlse_sig_hd); + [~,~,ber_mlse(r),~] = calc_ber(rx_bits.signal,Tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + fprintf('BER: %.2e \n',ber_mlse(r)); + fprintf('GMI MLSE: %.5f \n',gmi_mlse(r)); + + % Process through postfilter and MLSE + pf_ncoeffs = 1; + pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); + mlse_ = MLSE_viterbi("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); + [mlse_sig_sd,whitened_noise] = pf_.process(eq_signal_sd, eq_noise); + mlse_.DIR = pf_.coefficients; + mlse_sig_sd = mlse_.process(mlse_sig_sd); + mlse_sig_hd = PAMmapper(M, 0, "eth_style", 0).quantize(mlse_sig_sd); + 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)); - % Process through postfilter and MLSE - pf_ncoeffs = 1; - pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); - mlse_ = MLSE_viterbi("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); - [mlse_sig_sd,whitened_noise] = pf_.process(eq_signal_sd, eq_noise); - mlse_.DIR = pf_.coefficients; - mlse_sig_sd = mlse_.process(mlse_sig_sd); - mlse_sig_hd = PAMmapper(M, 0, "eth_style", 0).quantize(mlse_sig_sd); - 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