Add new functions

Functions/EQ_visuals/analyzeEQperformance.m

@@ -0,0 +1,150 @@
+function analyzeEQperformance(ref_bits,ref_symbols,rx_signal,eq_signal,eq_decisions,fsym,M,options)
+arguments
+    ref_bits 
+    ref_symbols 
+    rx_signal 
+    eq_signal 
+    eq_decisions 
+    fsym 
+    M 
+    options.postfilterclass
+    options.eqclass
+    options.mlseclass
+    options.db_precoded
+    options.displayname
+end
+
+% toolkit to visualize stuff related to DSP of IM/DD
+%%% Preps
+if isempty(eq_decisions)
+    eq_decisions = PAMmapper(M,0).quantize(eq_signal);
+end
+
+%%% Demap eqlzd signal to determine BER
+if options.db_precoded
+    eq_hd = PAMmapper(M,0).quantize(eq_signal);
+    eq_hd = Duobinary().encode(eq_hd);
+    eq_hd = Duobinary().decode(eq_hd,"M",M);
+    rx_bits = PAMmapper(M,0).demap(eq_hd);
+else
+    rx_bits = PAMmapper(M,0).demap(eq_signal);
+end
+
+[~,numerr,ber_sd,errpos] = calc_ber(rx_bits.signal,ref_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
+% fprintf('SD BER: %.2e \n',ber_sd);
+
+%%% Demap provided decisions to determine BER
+rx_bits = PAMmapper(M,0).demap(eq_decisions);
+[~,numerr,ber_mlse,errpos] = calc_ber(rx_bits.signal,ref_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
+% fprintf('MLSE BER: %.2e \n',ber_mlse);
+
+
+
+%%% Noise prior to DSP (is thius accurate with resampling?)
+rx_resampled = rx_signal.normalize("mode","rms").resample("fs_out",fsym);
+rx_noise = rx_resampled-ref_symbols;
+
+%%% Noise after to soft-decision DSP 
+eq_noise = eq_signal-ref_symbols;
+
+
+col = cbrewer2('paired',12);
+lines = numel(findall(figure(200), 'Type', 'Line'))+1;
+darkcoloridx = max(mod(2*lines,12),2);
+lightcoloridx = max(mod(2*lines-1,12),1);
+
+%%% Separate Classes
+constellation = unique(ref_symbols.signal);
+received_sd = NaN(numel(constellation),length(ref_symbols));
+received_hd = NaN(numel(constellation),length(ref_symbols));
+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,ref_symbols.signal==constellation(lvl)) = eq_signal.signal(ref_symbols.signal==constellation(lvl));
+    received_hd(lvl,ref_symbols.signal==constellation(lvl)) = eq_decisions.signal(ref_symbols.signal==constellation(lvl));
+end
+
+%%% bursts
+% Find differences between consecutive elements
+diff_indices = diff(errpos);
+
+% Identify the start of new sequences (when the difference is not 1)
+sequence_starts = [1, find(diff_indices ~= 1) + 1]; % Include the first index
+sequence_ends = [sequence_starts(2:end) - 1, length(errpos)]; % Calculate end indices
+
+% Initialize burst count and print bursts longer than 10
+burst_len = 1:10;
+burst_count = zeros(length(burst_len),1);
+for t = 1:numel(burst_len)
+    for i = 1:length(sequence_starts)
+        % Extract current sequence
+        current_burst = errpos(sequence_starts(i):sequence_ends(i));
+        % Check if the sequence length matches criterion
+        if length(current_burst) == burst_len(t)
+            burst_count(t) = burst_count(t) + 1;
+        end
+    end
+end
+
+burst_symbols = burst_count .* burst_len';
+burst_rate =  burst_symbols ;%./ length(rx_signal);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%%% Rx Spectrum %200
+rx_signal.spectrum("displayname",sprintf('Rx %d GBd PAM%d',fsym.*1e-9,M),'fignum',200,'normalizeTo0dB',1,'color',col(darkcoloridx,:));
+xline([-fsym/2,fsym/2].*1e-9,'Color',col(mod(lines,12)+2,:),'HandleVisibility','off');
+ylim([-20,3]);
+%%% EQ Spectrum
+% 
+
+%%% EQ Time Series %210
+showEQTimeSignal(eq_signal,ref_symbols)
+
+
+%%% EQ Noise Spectrum + inverted Postfilter %220 
+showEQNoisePSD(eq_noise,options.postfilterclass.burg_coeff,"fignum",220,"color",col(darkcoloridx,:));
+
+%%% EQ SNR Spectrum %230
+fft_length = 2^13;
+[s_lin,w] = pwelch(eq_signal.signal,hanning(fft_length),fft_length/2,fft_length,eq_signal.fs,"centered","psd","mean");
+[n_lin,w] = pwelch(eq_noise.signal,hanning(fft_length),fft_length/2,fft_length,eq_noise.fs,"centered","psd","mean");
+w = w.*1e-9;
+snr_dbm = 10*log10(s_lin./n_lin);
+
+figure(230)
+hold on
+plot(w,snr_dbm,'DisplayName','SNR after EQ','LineWidth',1,'Color',col(darkcoloridx,:));
+xlabel("Frequency in GHz");
+edgetick = 2^(nextpow2(eq_signal.fs*1e-9));
+xticks(-edgetick:16:edgetick);
+xlim([-128 128]);
+ylim([-20,35]);
+grid minor
+yticks(-200:10:100);
+grid on; grid minor;
+legend('Interpreter','none');
+title('Noise of soft decision signal (not MLSE)');
+
+%%% FFE histogram %240
+showLevelHistogram(eq_signal,ref_symbols)
+
+%%% Confusion Matrix
+showLevelConfusionMatrix(eq_decisions,ref_symbols,"fignum",250)
+
+
+%%% Burst Count
+figure(260)
+hold on
+plot(burst_len,burst_rate,'Marker','x','Color',col(darkcoloridx,:),"displayname",sprintf('Rx %d GBd PAM%d; %s',fsym.*1e-9,M,options.displayname));
+xlabel('length of error burst');
+ylabel('occurences')
+grid on
+set(gca, 'YScale', 'log');
+
+
+autoArrangeFigures
+
+end

Functions/EQ_visuals/showEQNoisePSD.m

@@ -0,0 +1,44 @@
+function showEQNoisePSD(eq_noise, options)
+arguments
+    eq_noise 
+    options.postfilter_taps = NaN
+    options.fignum (1,1) double = NaN % Default to NaN if not provided
+    options.displayname (1,:) char = '' % Default to an empty string if not provided
+    options.color = [0.2157    0.4941    0.7216];
+end
+
+    % Determine the figure number to use or create a new figure
+    if isnan(options.fignum)
+        fig = figure; % Create a new figure and get its handle
+    else
+        fig = figure(options.fignum); % Use the specified figure number
+    end
+    hold on
+    ax = gca;
+    % N = numel(ax.Children);
+    N = sum(arrayfun(@(x) strcmp(x.LineStyle, '-'), ax.Children));
+    cmap = linspecer(8);
+    options.color = cmap(mod(N, size(cmap, 1)) + 1, :);
+
+    % Ensure the figure is ready before calling spectrum
+    eq_noise.spectrum("displayname", options.displayname, "fignum", fig.Number, "normalizeTo0dB", 1,"color",options.color);
+    title('Noise of soft decision signal (not MLSE)')
+
+    if ~isnan(options.postfilter_taps)
+        % Hold on to the figure for further plotting
+        hold on;
+    
+        % Compute the frequency response of the postfilter
+        [h, w] = freqz(1, options.postfilter_taps, length(eq_noise), "whole", eq_noise.fs);
+        h = h / max(abs(h)); % Normalize the filter response
+    
+        % Adjust frequency axis to center at 0
+        w_ = (w - eq_noise.fs / 2);
+    
+        % Plot the inverted postfilter response
+        plot(w_ * 1e-9, 20 * log10(fftshift(abs(h))), 'DisplayName', ['Burg Coeffs: ', num2str(round(options.postfilter_taps, 2)), ' '], 'LineWidth', 1,'Color',options.color,'LineStyle','--');
+    
+        % Ensure a legend is displayed
+        legend('show');
+    end
+end

Functions/EQ_visuals/showEQNoiseSNR.m

@@ -0,0 +1,70 @@
+function showEQNoiseSNR(tx_signal, rx_signal, options)
+arguments
+    tx_signal
+    rx_signal 
+    options.fs_tx
+    options.fs_rx
+    options.fignum (1,1) double = NaN % Default to NaN if not provided
+    options.displayname (1,:) char = '' % Default to an empty string if not provided
+    options.color = [0.2157    0.4941    0.7216];
+end
+
+    % Determine the figure number to use or create a new figure
+    if isnan(options.fignum)
+        fig = figure; % Create a new figure and get its handle
+    else
+        fig = figure(options.fignum); % Use the specified figure number
+    end
+
+    if isa(tx_signal,'Signal')
+        options.fs_tx = tx_signal.fs;
+        tx_signal = tx_signal.signal;
+    end
+    if isa(rx_signal,'Signal')
+        options.fs_rx = rx_signal.fs;
+        rx_signal = rx_signal.signal;
+    end
+
+
+    hold on
+    ax = gca;
+
+    % N = numel(ax.Children);
+    N = sum(arrayfun(@(x) strcmp(x.LineStyle, '-'), ax.Children));
+    cmap = linspecer(8);
+    options.color = cmap(mod(N, size(cmap, 1)) + 1, :);
+
+    % Ensure the figure is ready before calling spectrum
+    title('SNR of received Signal')
+
+    fft_length = 2^(nextpow2(length(tx_signal))-7);
+    
+    [s_lin,w] = pwelch(tx_signal,hanning(fft_length),fft_length/2,fft_length,options.fs_tx,"centered","psd","mean");
+    [n_lin,w] = pwelch(rx_signal,hanning(fft_length),fft_length/2,fft_length,options.fs_rx,"centered","psd","mean");
+
+    w = w.*1e-9;
+    snr_dbm = 10*log10(s_lin./n_lin);
+
+    % figure(231)
+    hold on
+    plot(w,snr_dbm,'DisplayName','SNR','LineWidth',0.5,'Color',options.color);
+    xlabel("Frequency in GHz");
+
+    edgetick = 2^(nextpow2(options.fs_tx*1e-9));
+    ticks = -edgetick:16:edgetick;
+    xticks(ticks);
+
+    [~,b]=min(abs((-edgetick:16:edgetick)-max(w)));
+    xlim([-ticks(b+1) ticks(b+1)]);
+    
+    max_snr = ceil(max(snr_dbm)/10)*10;
+    min_snr = floor(min(snr_dbm)/10)*10;
+    ylim([min_snr,max_snr]);
+    yticks(-200:10:100);
+
+    grid on; grid minor;
+    legend('Interpreter','none');
+    title('Noise of soft decision signal (not MLSE)');
+
+
+end

Functions/EQ_visuals/showEQTimeSignal.m

@@ -0,0 +1,54 @@
+function showEQTimeSignal(eq_signal,ref_symbols,options)
+
+    arguments
+        eq_signal
+        ref_symbols
+        options.fignum (1,1) double = NaN % Default to NaN if not provided
+        options.displayname (1,:) char = '' % Default to an empty string if not provided
+        options.color = [0.2157    0.4941    0.7216];
+    end
+    % Determine the figure number to use or create a new figure
+    if isnan(options.fignum)
+        fig = figure; % Create a new figure and get its handle
+    else
+        fig = figure(options.fignum); % Use the specified figure number
+    end
+
+
+    M = numel(unique(ref_symbols.signal));
+    lvlcol = cbrewer2('Set1',M);
+    col = cbrewer2('paired',2);
+
+    eq_decisions = PAMmapper(M,0).quantize(eq_signal);
+    rx_bits = PAMmapper(M,0).demap(eq_signal);
+    ref_bits = PAMmapper(M,0).demap(ref_symbols);
+
+    %%% Separate Classes
+    constellation = unique(ref_symbols.signal);
+    received_sd = NaN(numel(constellation),length(ref_symbols));
+    received_hd = NaN(numel(constellation),length(ref_symbols));
+    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,ref_symbols.signal==constellation(lvl)) = eq_signal.signal(ref_symbols.signal==constellation(lvl));
+        received_hd(lvl,ref_symbols.signal==constellation(lvl)) = eq_decisions.signal(ref_symbols.signal==constellation(lvl));
+    end
+
+    [~,numerr,ber_sd,errpos] = calc_ber(rx_bits.signal,ref_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
+
+    %%% EQ Time Series %210
+    eq_signal.plot("fignum",fig.Number,"displayname",'Equalized Signal','color',col(1,:),'clear',1);
+    hold on
+
+    try
+        yline(PAMmapper(M,0).get_demodulation_thresholds,'HandleVisibility','off','LineStyle','--');
+    
+        for c = 1:M
+            scatter(errpos./eq_signal.fs,received_sd(c,errpos),1,'x','MarkerEdgeColor',lvlcol(c,:),'LineWidth',2,'DisplayName',sprintf('Tx Lvl: %d',c));
+        end                 
+    end
+
+end
+

Functions/EQ_visuals/showEQcoefficients.m

@@ -0,0 +1,58 @@
+function showEQcoefficients(n1, n2, n3, options)
+    % Show filter coefficients as stem plot
+    % n1, n2, and n3 in different subplots
+    % Scale all y-axis to -1 and 1
+
+    arguments
+        n1
+        n2
+        n3
+        options.fignum (1,1) double = NaN % Default to NaN if not provided
+        options.displayname (1,:) char = '' % Default to an empty string if not provided
+        options.color = [0.2157, 0.4941, 0.7216];
+        options.clf = 0; % Clear figure before plotting new
+    end
+
+    % Determine the figure number to use or create a new figure
+    if isnan(options.fignum)
+        fig = figure; % Create a new figure and get its handle
+    else
+        fig = figure(options.fignum); % Use the specified figure number
+    end
+
+    if options.clf
+        clf(fig); % Clear the figure if requested
+    end
+
+    hold on
+    ax = gca;
+    N = numel(ax.Children);
+
+    % Set up a colormap for consistent coloring
+    cmap = linspecer(8);
+    options.color = cmap(mod(N, size(cmap, 1)) + 1, :);
+
+    % Create subplots for n1, n2, n3
+    for i = 1:3
+        subplot(3, 1, i);
+        switch i
+            case 1
+                stem(n1, 'Color', options.color, 'LineWidth', 1,'Marker','.','MarkerSize',10);
+                title(sprintf('1st order Filter Coefficients: %d',numel(n1)));
+            case 2
+                stem(n2, 'Color', options.color, 'LineWidth', 1,'Marker','.','MarkerSize',10);
+                title(sprintf('2nd order Filter Coefficients: %d',numel(n2)));
+            case 3
+                stem(n3, 'Color', options.color, 'LineWidth', 1,'Marker','.','MarkerSize',10);
+                title(sprintf('3rd order Filter Coefficients: %d',numel(n3)));
+        end
+        ylim([-1, 1]); % Scale y-axis to -1 and 1
+        grid on;
+        grid minor
+        xlabel('Coefficient Index');
+        ylabel('Amplitude');
+    end
+
+    % Ensure the layout is tight for better visibility
+    sgtitle('Filter Coefficients'); % Overall title
+end

Functions/EQ_visuals/showErrorBurstCount.m

@@ -0,0 +1,47 @@
+function showErrorBurstCount(eq_signal,ref_symbols,options)
+arguments
+    eq_signal
+    ref_symbols
+    options.fignum (1,1) double = NaN % Default to NaN if not provided
+    options.displayname (1,:) char = '' % Default to an empty string if not provided
+end
+
+    % Determine the figure number to use or create a new figure
+    if isnan(options.fignum)
+        fig = figure; % Create a new figure and get its handle
+    else
+        fig = figure(options.fignum); % Use the specified figure number
+    end
+
+    if numel(unique(eq_signal.signal)) > 20
+        M = numel(unique(ref_symbols.signal));
+        eq_signal = PAMmapper(M,0).quantize(eq_signal);
+    end
+    
+    diff_indices = eq_signal.signal == ref_symbols.signal;
+
+
+    % Identify the start of new sequences (when the difference is not 1)
+    sequence_starts = [1, find(diff_indices ~= 1) + 1]; % Include the first index
+    sequence_ends = [sequence_starts(2:end) - 1, length(errpos)]; % Calculate end indices
+
+    % Initialize burst count and print bursts longer than 10
+    burst_len = 1:10;
+    burst_count = zeros(length(burst_len),1);
+    for t = 1:numel(burst_len)
+        for i = 1:length(sequence_starts)
+            % Extract current sequence
+            current_burst = errpos(sequence_starts(i):sequence_ends(i));
+            % Check if the sequence length matches criterion
+            if length(current_burst) == burst_len(t)
+                burst_count(t) = burst_count(t) + 1;
+            end
+        end
+    end
+
+    burst_symbols = burst_count .* burst_len';
+    burst_rate =  burst_symbols ;%./ length(rx_signal);
+
+
+
+end

Functions/EQ_visuals/showLevelConfusionMatrix.m

@@ -0,0 +1,27 @@
+function showLevelConfusionMatrix(decided_symbols,ref_symbols,options)
+arguments
+    decided_symbols
+    ref_symbols
+    options.M = NaN
+    options.fignum (1,1) double = NaN % Default to NaN if not provided
+    options.displayname (1,:) char = '' % Default to an empty string if not provided
+end
+
+% Determine the figure number to use or create a new figure
+if isnan(options.fignum)
+    fig = figure; % Create a new figure and get its handle
+else
+    fig = figure(options.fignum); % Use the specified figure number
+end
+
+if length(unique(decided_symbols.signal))>20
+    assert(~isnan(options.M),'Provide either decided symbol sequence or modulation order M (PAM-4 -> M=4)');
+    decided_symbols = PAMmapper(options.M,0).quantize(decided_symbols);
+end
+
+%%% Confusion Matrix
+cm = confusionchart(ref_symbols.signal,decided_symbols.signal,'RowSummary','row-normalized','ColumnSummary','column-normalized','Title','Confusion Matrix','XLabel','Decisions','YLabel','Transmitted');
+
+
+end
+

Functions/EQ_visuals/showLevelHistogram.m

@@ -0,0 +1,50 @@
+function showLevelHistogram(eq_signal,ref_symbols,options)
+arguments
+    eq_signal
+    ref_symbols
+    options.fignum (1,1) double = NaN % Default to NaN if not provided
+    options.displayname (1,:) char = '' % Default to an empty string if not provided
+end
+
+    if isa(eq_signal,'Signal')
+        eq_signal = eq_signal.signal;
+    end
+    if isa(ref_symbols,'Signal')
+        ref_symbols = ref_symbols.signal;
+    end
+
+    % Determine the figure number to use or create a new figure
+    if isnan(options.fignum)
+        fig = figure; % Create a new figure and get its handle
+    else
+        fig = figure(options.fignum); % Use the specified figure number
+    end
+
+    %%% Separate Classes
+    constellation = unique(ref_symbols);
+    received_sd = NaN(numel(constellation),length(ref_symbols));
+    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,ref_symbols==constellation(lvl)) = eq_signal(ref_symbols==constellation(lvl));
+    end
+    
+    
+    
+    %%% FFE histogram 
+    clf
+    for lvl = 1:numel(constellation)
+        intermediate = received_sd(lvl,:);
+        cnt(lvl) = round(numel(intermediate(~isnan(intermediate)))./length(eq_signal),3).*100;
+        hold on
+        histogram(received_sd(lvl,:),1000,"EdgeAlpha",0,'DisplayName',['Lvl ',num2str(lvl),' | ',num2str(cnt(lvl)),' %'],'FaceColor',lvlcol(lvl,:),'Normalization','pdf');
+    end
+    legend
+    grid on
+
+
+
+end
+

Functions/Metrics/air_garcia_implementation.m

@@ -0,0 +1,156 @@
+function air = air_garcia_implementation(x,r,idx_tx,Px,M_training)
+
+MAX_MEMORY = 200e6; % maximum allowed size for a matrix
+
+if nargin == 3
+    Px = [];
+    M_training = [];
+end
+
+if nargin == 4
+    M_training = [];
+end
+
+
+% if input is complex, separate into real and imaginary parts
+if any(imag(x(:))~=0) || any(imag(r(:))~=0)
+    x = [real(x); imag(x)];
+    r = [real(r); imag(r)];
+end
+
+D = size(x, 1); % D = 2 if complex x
+N = size(x, 2); % number of constellation points
+M = size(r, 2); % number of samples
+
+% set default training set size
+if isempty(M_training)
+    M_training = ceil(0.3*M);
+end
+
+M_testing = M - M_training;
+
+% Training: estimate parameters of the conditionally Gaussian model
+% sort according to transmit index
+[idx_tx_training, idx_sort] = sort(idx_tx(1:M_training));
+r_training = r(:, idx_sort);
+i_bounds = zeros(1, N+1);
+
+% compute conditional means and covariance matrices
+C_n = zeros(D, D, N);
+det_n = zeros(1, N);
+
+for n=1:N
+    % find how many times x(:, n) was transmitted and update i_bounds
+    N_current_x = find(idx_tx_training((i_bounds(n)+1):end)==n, 1, 'last');
+    if isempty(N_current_x), N_current_x=0; end
+    i_bounds(n+1) = i_bounds(n) + N_current_x;
+
+    if N_current_x > 0
+        % Compute mu_n=E[Y|X=x_n] according to Eq. (14) and store it in
+        % x(:, n) to save space
+        x(:, n) = sum(r_training(:, (i_bounds(n)+1):i_bounds(n+1)), 2)/(i_bounds(n+1)-i_bounds(n));
+
+        % compute C_n=cov[Y|X=x_n] according to Eq. (15)
+        r_meanfree = r_training(:, (i_bounds(n)+1):i_bounds(n+1)) - x(:, n);
+        C_n(:, :, n) = (r_meanfree*r_meanfree')/(i_bounds(n+1)-i_bounds(n));
+        % store also the determinant of C(:, :, n)
+        det_n(n) = det(C_n(:, :, n));
+
+        % if the determinant is 0, or if the matrix is badly conditioned,
+        % regularize by adding a small identity matrix. Note that we do
+        % need the check for 0 determinant, in case a cloud has exactly 0
+        % variance according to the training set
+        if det_n(n)==0 || cond(C_n(:, :, n))>1e16
+            C_n(:, :, n) = C_n(:, :, n) + 5 * eps * eye(D);
+            det_n(n) = (5*eps)^D;
+        end
+
+    end
+end
+
+% uniform input pmf Px if not provided
+if isempty(Px)
+    Px = repmat(1/N, [1, N]);
+end
+
+
+% extract testing set and sort it according to transmit index
+[idx_tx_testing, idx_sort] = sort(idx_tx((M_training+1):M));
+r_testing = r(:, M_training+idx_sort);
+
+% computation of h(Y|X)
+h_Y_X = 0;
+i_bounds_testing = zeros(1, N+1);
+% loop over constellation points to compute h(Y|X)
+for n = 1:N
+    % find how many times x(:, n) was transmitted and update
+    % i_bounds_testing
+    N_current_x = find(idx_tx_testing((i_bounds_testing(n)+1):end)==n, 1, 'last');
+    if isempty(N_current_x), N_current_x=0; end
+    i_bounds_testing(n+1) = i_bounds_testing(n) + N_current_x;
+    % add the corresponding contribution to the mutual information (two
+    % first lines of Eq. (17)). This, together with
+    % D/2*log2(2*pi) after the end of the loop, gives h(Y|X)
+    h_Y_X = h_Y_X + N_current_x * log2(det_n(n))/2+...
+        sum(sum(conj(r_testing(:, (i_bounds_testing(n)+1):i_bounds_testing(n+1))-x(:, n)).*(C_n(:, :, n)\(r_testing(:, (i_bounds_testing(n)+1):i_bounds_testing(n+1))-x(:, n)))))/2/log(2);
+
+    
+end
+h_Y_X = D/2*log2(2*pi) + h_Y_X/M_testing;
+
+% When computing log(py), we might run out of memory. If necessary, we
+% doe the computation in blocks
+logpy = zeros(1, M_testing);
+BLOCK_SIZE = floor(MAX_MEMORY/N);
+N_blocks = ceil(M_testing/BLOCK_SIZE);
+
+% loop over blocks of symbols. This loop can be replaced by parfor to allow
+% parallel computation
+for i_block = 1:N_blocks
+
+    logpy_cur = zeros(1, M_testing);
+
+    % beginning of block
+    i_start = (i_block-1) * BLOCK_SIZE + 1;
+    % end of block
+    i_end = min(M_testing, i_block*BLOCK_SIZE);
+    % block size
+    current_block_size = i_end-i_start+1;
+
+    % compute exponents of third line of (17)
+    exponents = zeros(N, current_block_size);
+    for n = 1:N
+        exponents(n, :) = -log(det_n(n))/2-real(sum(conj(r_testing(:, i_start:i_end)-x(:, n)).*(C_n(:, :, n)\(r_testing(:, i_start:i_end)-x(:, n))), 1))/2;
+        %sum über 2 einträge von r
+    end
+
+    % compute third line of Eq. (17). Use a custom function
+    % that computes log(sum(exp(x))) avoiding overflow errors
+    logpy_cur(i_start:i_end) = math_logsumexp(log(Px(:))+exponents, 1);
+    logpy = logpy + logpy_cur;
+end
+
+% output entropy h(Y)
+h_Y = D/2*log2(2*pi) - mean(logpy)/log(2);%log basis change
+
+% compute mutual information
+air = h_Y - h_Y_X;
+
+
+end
+
+
+function [y] = math_logsumexp(x, dim)
+%[y] = math_logsumexp(x, dim)
+%  Computes log(sum(exp(x), dim)), avoiding overflow errors when one of the
+%  x is large.
+
+if nargin<2 || isempty(dim)
+    m = max(x);
+    y = m + log(sum(exp(x-m)));
+else
+    m = max(x, [], dim);
+    y = m + log(sum(exp(x-m), dim));
+end
+end
+

Functions/Metrics/calc_air.m

@@ -0,0 +1,48 @@
+function [ach_inf_rate] = calc_air(test_signal,reference_signal,options)
+% Calculation of AIR acc. to J. Kozesnik, „Numerically Computing Achievable Rates of  Memoryless Channels“, Francisco Javier Garcıa-Gomez,  doi: 10.1007/978-94-009-9857-5.
+% Implementation is not accessible, I mailed TUM to get the code...
+
+arguments(Input)
+    test_signal;
+    reference_signal;
+    options.skip_front = 0;
+    options.skip_end = 0;
+    options.returnErrorLocation = 0;
+end
+
+options.skip_end = abs(options.skip_end);
+options.skip_front = abs(options.skip_front);
+
+assert((options.skip_end+options.skip_front)<length(test_signal),"You can not skip more bits than overall length of data! Set skip_front or skip_end to lower value or check data_in");
+
+if isa(reference_signal,'Signal')
+    reference_signal = reference_signal.signal;
+end
+
+if isa(test_signal,'Signal')
+    test_signal = test_signal.signal;
+end
+
+% TRIM
+[test_signal,reference_signal]=trimseq(test_signal,reference_signal,options.skip_front,options.skip_end);
+
+% CALC EVM
+%%% new implementation of AIR
+constellation = unique(reference_signal);
+reference_idx = arrayfun(@(x) find(constellation == x, 1), reference_signal);
+ach_inf_rate = air_garcia_implementation(constellation',test_signal',reference_idx');
+
+
+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

Functions/Metrics/calc_ber.m

@@ -0,0 +1,52 @@
+function [bits,errors,ber,errorIndice] = calc_ber(data_in,data_ref,options)
+
+arguments(Input)
+    data_in;
+    data_ref;
+    options.skip_front = 0;
+    options.skip_end = 0;
+    options.returnErrorLocation = 0;
+end
+options.skip_end = abs(options.skip_end);
+options.skip_front = abs(options.skip_front);
+
+assert((options.skip_end+options.skip_front)<length(data_in),"You can not skip more bits than overall length of data! Set skip_front or skip_end to lower value or check data_in");
+
+errorIndice= [];
+
+% trim sequence to given start; stop. trim reference if signal is shorter
+[data_in,data_ref]=trimseq(data_in,data_ref,options.skip_front,options.skip_end);
+
+if length(data_ref) == length(data_in)
+
+    bits = numel(data_in(:,options.skip_front+1:end));
+
+    if options.returnErrorLocation == 0
+        errors = sum( data_in ~= data_ref,"all" );
+    else
+        errorIndice = sum(data_in ~= data_ref,1);
+        errors = sum(errorIndice ,"all" );
+        [~,errorIndice] = find(errorIndice~=0);
+        errorIndice = errorIndice+options.skip_front;
+    end
+
+    % Determine BER
+    ber = sum(errors)/sum(bits);
+
+else
+    error('Sequence length does not match');
+end
+
+    function [data_,reference_]=trimseq(data,reference,skipstart,skip_end)
+
+        data_  =  logical(data(skipstart+1:end-skip_end,:))';
+
+        delta_bits = length(reference) - length(data);
+
+        skip_end = delta_bits + skip_end;
+
+        reference_ = logical(reference(skipstart+1:end-skip_end,:))';
+
+    end
+
+end

Functions/Metrics/calc_evm.m

@@ -0,0 +1,66 @@
+function [evm_total,evm_lvl] = calc_evm(test_signal,reference_signal,options)
+
+arguments(Input)
+    test_signal;
+    reference_signal;
+    options.skip_front = 0;
+    options.skip_end = 0;
+    options.returnErrorLocation = 0;
+end
+
+options.skip_end = abs(options.skip_end);
+options.skip_front = abs(options.skip_front);
+
+assert((options.skip_end+options.skip_front)<length(test_signal),"You can not skip more bits than overall length of data! Set skip_front or skip_end to lower value or check data_in");
+
+if isa(reference_signal,'Signal')
+    reference_signal = reference_signal.signal;
+end
+
+if isa(test_signal,'Signal')
+    test_signal = test_signal.signal;
+end
+
+% TRIM
+[test_signal,reference_signal]=trimseq(test_signal,reference_signal,options.skip_front,options.skip_end);
+
+% CALC EVM
+[evm_total,evm_lvl] = calc_evm_(test_signal,reference_signal);
+
+
+function [evm_total,evm_lvl] = calc_evm_(test_signal,reference_signal)
+
+    assert(length(test_signal) == length(reference_signal),"Sequence length does not match");
+
+    error_vector = (test_signal-reference_signal);
+
+    %%% Overall EVM
+    evm_total = rms(error_vector);
+
+    try
+        %%% Per Level EVM
+        k = unique(reference_signal);
+        for lvl = 1:length(k)
+            lvl_errors = error_vector(reference_signal==k(lvl));
+            evm_lvl(lvl) =  rms(lvl_errors);
+        end
+    catch
+        evm_lvl = NaN;
+        warning('No EVM per level calculated')
+    end
+
+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

Functions/Metrics/calc_ngmi.m

@@ -0,0 +1,123 @@
+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)
+
+arguments(Input)
+    test_signal;
+    reference_signal;
+    options.skip_front = 0;
+    options.skip_end = 0;
+    options.returnErrorLocation = 0;
+end
+
+options.skip_end = abs(options.skip_end);
+options.skip_front = abs(options.skip_front);
+
+assert((options.skip_end+options.skip_front)<length(test_signal),"You can not skip more bits than overall length of data! Set skip_front or skip_end to lower value or check data_in");
+
+if isa(reference_signal,'Signal')
+    reference_signal = reference_signal.signal;
+end
+
+if isa(test_signal,'Signal')
+    test_signal = test_signal.signal;
+end
+
+% TRIM
+[test_signal,reference_signal]=trimseq(test_signal,reference_signal,options.skip_front,options.skip_end);
+
+% CALC EVM
+[GMI,NGMI] = calc_ngmi_(test_signal,reference_signal);
+
+
+    function [GMI,NGMI] = calc_ngmi_(test_signal,reference_signal)
+
+        assert(length(test_signal) == length(reference_signal),"Sequence length does not match");
+
+        %%% implemented according to [1] 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, S. 1–3. doi: 10.1109/ECOC.2017.8345872.
+
+        error_vector = (test_signal-reference_signal);
+        sigma2 = var(error_vector); %noise variance
+
+        %%% 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
+
+        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 = 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
+
+        % GMI computation
+        noise_impact_term = 0;
+        for k = 1:N
+            y_k = test_signal(k); % Current received sample
+            [~, closest_symbol_idx] = min(abs(symbols - y_k)); % Closest symbol index
+            closest_symbol = symbols(closest_symbol_idx);     % Closest symbol
+
+            for i = 1:m
+                % 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));
+
+                % 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
+                noise_impact_term = noise_impact_term + log2(numerator / denominator);
+            end
+        end
+
+        % Normalize the noise impact term by N
+        noise_impact_term = noise_impact_term / N;
+
+        % GMI
+        GMI = H_X + noise_impact_term;
+        NGMI = GMI / 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