From d099efea030a595d685db5ddbeae8e0a3b232007 Mon Sep 17 00:00:00 2001 From: Silas Oettinghaus Date: Mon, 17 Feb 2025 21:31:12 +0100 Subject: [PATCH] Add new functions --- .../04_DSP/Sequence Detection/MLSE_viterbi.m | 177 +++++++++++++++ Classes/04_DSP/TransmissionPerformance.m | 214 ++++++++++++++++++ Classes/DataBaseHandler/copyStylingFrom.m | 76 +++++++ Functions/EQ_visuals/analyzeEQperformance.m | 150 ++++++++++++ Functions/EQ_visuals/showEQNoisePSD.m | 44 ++++ Functions/EQ_visuals/showEQNoiseSNR.m | 70 ++++++ Functions/EQ_visuals/showEQTimeSignal.m | 54 +++++ Functions/EQ_visuals/showEQcoefficients.m | 58 +++++ Functions/EQ_visuals/showErrorBurstCount.m | 47 ++++ .../EQ_visuals/showLevelConfusionMatrix.m | 27 +++ Functions/EQ_visuals/showLevelHistogram.m | 50 ++++ Functions/Metrics/air_garcia_implementation.m | 156 +++++++++++++ Functions/Metrics/calc_air.m | 48 ++++ Functions/Metrics/calc_ber.m | 52 +++++ Functions/Metrics/calc_evm.m | 66 ++++++ Functions/Metrics/calc_ngmi.m | 123 ++++++++++ .../imdd_simulation_minimal.m | 176 -------------- projects/A_minimal_example/run_loops.m | 23 -- .../A_minimal_example/submit_simulations.m | 72 ------ .../auswertung/analyze_mu_dc.m | 89 ++++++++ projects/IMDD_base_system/imdd_it.m | 2 +- 21 files changed, 1502 insertions(+), 272 deletions(-) create mode 100644 Classes/04_DSP/Sequence Detection/MLSE_viterbi.m create mode 100644 Classes/04_DSP/TransmissionPerformance.m create mode 100644 Classes/DataBaseHandler/copyStylingFrom.m create mode 100644 Functions/EQ_visuals/analyzeEQperformance.m create mode 100644 Functions/EQ_visuals/showEQNoisePSD.m create mode 100644 Functions/EQ_visuals/showEQNoiseSNR.m create mode 100644 Functions/EQ_visuals/showEQTimeSignal.m create mode 100644 Functions/EQ_visuals/showEQcoefficients.m create mode 100644 Functions/EQ_visuals/showErrorBurstCount.m create mode 100644 Functions/EQ_visuals/showLevelConfusionMatrix.m create mode 100644 Functions/EQ_visuals/showLevelHistogram.m create mode 100644 Functions/Metrics/air_garcia_implementation.m create mode 100644 Functions/Metrics/calc_air.m create mode 100644 Functions/Metrics/calc_ber.m create mode 100644 Functions/Metrics/calc_evm.m create mode 100644 Functions/Metrics/calc_ngmi.m delete mode 100644 projects/A_minimal_example/imdd_simulation_minimal.m delete mode 100644 projects/A_minimal_example/run_loops.m delete mode 100644 projects/A_minimal_example/submit_simulations.m create mode 100644 projects/HighSpeedExperiment_2024/auswertung/analyze_mu_dc.m diff --git a/Classes/04_DSP/Sequence Detection/MLSE_viterbi.m b/Classes/04_DSP/Sequence Detection/MLSE_viterbi.m new file mode 100644 index 0000000..9a67728 --- /dev/null +++ b/Classes/04_DSP/Sequence Detection/MLSE_viterbi.m @@ -0,0 +1,177 @@ +classdef MLSE_viterbi < 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_viterbi(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 = process(obj,signalclass) + + data_in = signalclass.signal; + + data_out = obj.process_(data_in); + + signalclass.signal = data_out; + + end + + function data_out = process_(obj,data_in) + + + % 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 + + + %%%%%%%%%%%%%%%%%%%%%%%%%%%% WORKING + % 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 + + % OLD + % initilaize the output vector + data_out = NaN(size(data_in)); + + sum_path_metrics = zeros(length(states),length(states)); + + % first trellis path + % euclidian distance as path metric + path_metrics = (abs(repmat(data_in(1),size(noise_free_received)) - noise_free_received)).^2; + sum_path_metrics = sum_path_metrics + path_metrics; % calculation of all possible sum path metrics + [sum_path_metrics_res(:,1),path_idx(:,1)] = min(sum_path_metrics,[],2); % find the best path to each state, store sum path metric and predecessor for each state + + % remaining trellis paths + for n = 2:length(data_in) + sum_path_metrics = repmat(sum_path_metrics_res(:,n-1).',length(states),1); + % path_metrics = (repmat(data_in(n),size(noise_free_received)) - noise_free_received).^2;%.*prob_mat; + path_metrics = (abs(repmat(data_in(n),size(noise_free_received)) - noise_free_received)).^2; + sum_path_metrics = sum_path_metrics + path_metrics; + [sum_path_metrics_res(:,n),path_idx(:,n)] = min(sum_path_metrics,[],2); + end + + %% trace back + ideal_path = NaN(1,length(data_in)+1); + % find ideal trellis path by going through the trellis + % backwards + [~,ideal_path(length(data_in)+1)] = min(sum_path_metrics_res(:,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 + ideal_path(h) = path_idx(ideal_path(h+1),h); + end + + idx_out = ideal_path(2:length(data_in)+1); + %%%%%%%%%%%%%%%%%%%%%%%%%%%% WORKING + + + + + if obj.duobinary_output + %use duobinary encoder, output is already scaled inside this + %one + data_out = Duobinary().encode(first_sym(idx_out)); + + else + % + data_out(1:length(data_in)) = first_sym(idx_out); + % scale to rms = 1 using the standard sqrt() expressions + if obj.M == 4 + data_out = data_out./sqrt(5); + elseif obj.M == 6 + data_out = data_out./sqrt(10); + elseif obj.M == 8 + data_out = data_out./sqrt(21); + end + + end + + + end + end + + methods (Access=private) + % Cant be seen from outside! So put all your functions here that can/ + % shall not be called from outside + + + end +end diff --git a/Classes/04_DSP/TransmissionPerformance.m b/Classes/04_DSP/TransmissionPerformance.m new file mode 100644 index 0000000..49ca0c7 --- /dev/null +++ b/Classes/04_DSP/TransmissionPerformance.m @@ -0,0 +1,214 @@ +classdef TransmissionPerformance + % TransmissionPerformance + % + % This class calculates the best possible net data rate (in bits/s) + % from a given gross rate (bits/s) and a measured channel quality + % parameter (either bit error ratio (BER) or NGMI). The class uses + % lookup tables for three different FEC modes: + % + % 1. SD+HD concatenated (uses overall code rate and an NGMI threshold) + % 2. SD+HD concatenated with punctured LDPC (uses overall code rate and an NGMI threshold) + % 3. HD-only (uses a code rate and a BER threshold) + % + % The basic algorithm is as follows: + % 1. Given a measured quality (NGMI or BER) for each measurement, find in + % the table the first (i.e. “best”) entry that is satisfied by the + % measured value. For NGMI the condition is measured NGMI >= threshold, + % while for BER the condition is measured BER <= threshold. + % 2. Use the corresponding overall code rate to compute: + % + % net rate = gross rate * code rate. + % + % USAGE EXAMPLE (array inputs): + % + % tp = TransmissionPerformance; + % % Suppose we have 3 measurements: + % % Gross rates: 10, 15, and 20 Gbps + % % Measured NGMI: 0.92, 0.88, and 0.85 + % % Measured BER: 1e-4, 7e-3, and 2e-2 + % netrates = tp.calculatenetrate([10e9 15e9 20e9], ... + % 'NGMI', [0.92 0.88 0.85], ... + % 'BER', [1e-4 7e-3 2e-2]); + % + % The returned structure netrates will contain the net rates (and the + % used code rates and thresholds) for each mode in vector form. + + properties (Constant) + %% Lookup Table for SD+HD concatenated (using NGMI) + % Overall code rate and NGMI threshold for each row. + OVERALLCODERATES_SDHD = [0.7519, 0.7602, 0.7684, 0.7766, 0.7850, ... + 0.7932, 0.8014, 0.8098, 0.8180, 0.8262, ... + 0.8345, 0.8428, 0.8510, 0.8593, 0.8676, ... + 0.8733, 0.8790, 0.8848, 0.8905, 0.8962, ... + 0.9019, 0.9077, 0.9134, 0.9191, 0.9248, ... + 0.9306, 0.9363, 0.9420, 0.9477, 0.9535]; + NGMITHRESHOLDS_SDHD = [0.8116, 0.8167, 0.8241, 0.8317, 0.8401, ... + 0.8459, 0.8512, 0.8574, 0.8685, 0.8746, ... + 0.8829, 0.8892, 0.8958, 0.9022, 0.9090,... + 0.9150, 0.9210, 0.9270, 0.9330, 0.9390, ... + 0.9450, 0.9510, 0.9570, 0.9630, 0.9690, ... + 0.9750, 0.9810, 0.9870, 0.9930, 0.9990]; + + ISPUNCT_LDPI = [zeros(1,15),ones(1,15)]; + + %% Lookup Table for HD-only mode (using BER) + % For HD-only, the table gives the code rate and the maximum acceptable BER. + CODE_RATES_HD = [0.7500, 0.8000, 0.8123, 0.8333, 0.8571, ... + 0.8750, 0.8889, 0.9000, 0.9091, 0.9167, 0.9412]; + BERTHRESHOLDS_HD = [2.12e-2, 1.76e-2, 1.62e-2, 1.44e-2, 1.25e-2, ... + 1.03e-2, 9.29e-3, 8.33e-3, 7.54e-3, 7.04e-3, 4.70e-3]; + + %% LUT for KP4-FEC and Inner Code https://grouper.ieee.org/groups/802/3/dj/public/23_03/patra_3dj_01b_2303.pdf + CODE_RATE_KP4_AND_INNER = [0.885799]; + BERTHRESHOLDS_KP4_AND_INNER = 4.85e-3; + + + + end + + methods + function netrates = calculateNetRate(obj, grossRate, varargin) + % calculatenetrate Calculate the net data rate(s) from measured data. + % + % netrates = tp.calculatenetrate(grossRate, 'BER', berValue, 'NGMI', ngmiValue) + % + % INPUTS: + % grossRate - (scalar or vector) gross data rate(s) [bits/s] + % + % Optional name/value pairs: + % 'BER' - (scalar or vector) measured bit error ratio(s) + % 'NGMI' - (scalar or vector) measured NGMI value(s) + % + % At least one of 'BER' or 'NGMI' must be provided. + % + % OUTPUT: + % netrates - a structure with the following fields (if available): + % .SDHD - for SD+HD concatenated (using NGMI) + % .SDHD_LDPC - for SD+HD with punctured LDPC (using NGMI) + % .HD - for HD-only (using BER) + % + % Each sub-structure contains fields: + % .CodeRate - the chosen overall code rate from the table (vector) + % .Threshold - the threshold value used from the table (vector) + % . - net rate computed as grossRate * CodeRate (vector) + + % Parse inputs. + p = inputParser; + addRequired(p, 'grossRate', @(x) isnumeric(x)); + addParameter(p, 'BER', [], @(x) isnumeric(x)); + addParameter(p, 'NGMI', [], @(x) isnumeric(x)); + parse(p, grossRate, varargin{:}); + + ber = p.Results.BER; + ngmi = p.Results.NGMI; + + if isempty(ber) && isempty(ngmi) + error('At least one of ''BER'' or ''NGMI'' must be provided.'); + end + + % Determine the number of measurements. + numMeasurements = max([numel(grossRate), numel(ngmi), numel(ber)]); + + % Broadcast scalars if needed. + if isscalar(grossRate) && numMeasurements > 1 + grossRate = repmat(grossRate, 1, numMeasurements); + elseif numel(grossRate) ~= numMeasurements + error('grossRate must be scalar or have %d elements.', numMeasurements); + end + + if ~isempty(ngmi) + if isscalar(ngmi) && numMeasurements > 1 + ngmi = repmat(ngmi, 1, numMeasurements); + elseif numel(ngmi) ~= numMeasurements + error('NGMI must be scalar or have %d elements.', numMeasurements); + end + end + + if ~isempty(ber) + if isscalar(ber) && numMeasurements > 1 + ber = repmat(ber, 1, numMeasurements); + elseif numel(ber) ~= numMeasurements + error('BER must be scalar or have %d elements.', numMeasurements); + end + end + + % Initialize the output structure. + netrates = struct; + if ~isempty(ngmi) + netrates.SDHD.GrossRate = NaN(1, numMeasurements); + netrates.SDHD.NetRate = NaN(1, numMeasurements); + netrates.SDHD.punctLDPI = NaN(1, numMeasurements); + netrates.SDHD.CodeRate = NaN(1, numMeasurements); + netrates.SDHD.Threshold = NaN(1, numMeasurements); + end + if ~isempty(ber) + netrates.HD.GrossRate = NaN(1, numMeasurements); + netrates.HD.NetRate = NaN(1, numMeasurements); + netrates.HD.CodeRate = NaN(1, numMeasurements); + netrates.HD.Threshold = NaN(1, numMeasurements); + + netrates.KP4_hamming.GrossRate = NaN(1, numMeasurements); + netrates.KP4_hamming.NetRate = NaN(1, numMeasurements); + netrates.KP4_hamming.CodeRate = NaN(1, numMeasurements); + netrates.KP4_hamming.Threshold = NaN(1, numMeasurements); + end + + % Process each measurement individually. + for i = 1:numMeasurements + + % --- NGMI-based modes --- + if ~isempty(ngmi) + % SD+HD concatenated mode. + idx = find(obj.NGMITHRESHOLDS_SDHD <= ngmi(i), 1, 'last'); + if ~isempty(idx) + codeRate = obj.OVERALLCODERATES_SDHD(idx); + netrates.SDHD.NetRate(i) = grossRate(i) * codeRate; + netrates.SDHD.GrossRate(i) = grossRate(i) ; + netrates.SDHD.CodeRate(i) = codeRate; + netrates.SDHD.Threshold(i) = obj.NGMITHRESHOLDS_SDHD(idx); + netrates.SDHD.punctLDPI(i) = obj.ISPUNCT_LDPI(idx); + end + + end + + % --- BER-based mode (HD-only) --- + if ~isempty(ber) + % Loop through the BER thresholds from the best (highest code rate) + % to the worst until the measured BER is acceptable. + idxBER = []; + for j = length(obj.BERTHRESHOLDS_HD):-1:1 + if ber(i) <= obj.BERTHRESHOLDS_HD(j) + idxBER = j; + break; + end + end + if ~isempty(idxBER) + codeRate = obj.CODE_RATES_HD(idxBER); + netrates.HD.NetRate(i) = grossRate(i) * codeRate; + netrates.HD.GrossRate(i) = grossRate(i) ; + netrates.HD.CodeRate(i) = codeRate; + netrates.HD.Threshold(i) = obj.BERTHRESHOLDS_HD(idxBER); + end + idxBER = []; + for j = length(obj.BERTHRESHOLDS_KP4_AND_INNER):-1:1 + if ber(i) <= obj.BERTHRESHOLDS_KP4_AND_INNER(j) + idxBER = j; + break; + end + end + if ~isempty(idxBER) + codeRate = obj.CODE_RATE_KP4_AND_INNER(idxBER); + netrates.KP4_hamming.NetRate(i) = grossRate(i) * codeRate; + netrates.KP4_hamming.GrossRate(i) = grossRate(i) ; + netrates.KP4_hamming.CodeRate(i) = codeRate; + netrates.KP4_hamming.Threshold(i) = obj.BERTHRESHOLDS_KP4_AND_INNER(idxBER); + end + + + end + + + end + end + end +end diff --git a/Classes/DataBaseHandler/copyStylingFrom.m b/Classes/DataBaseHandler/copyStylingFrom.m new file mode 100644 index 0000000..aaf110f --- /dev/null +++ b/Classes/DataBaseHandler/copyStylingFrom.m @@ -0,0 +1,76 @@ +function copyStylingFrom(figNumSource, figNumTgt) + % Get handles to the source and target figures + sourceFig = figure(figNumSource); + targetFig = figure(figNumTgt); + + % Get axes of source and target figures + sourceAxes = findall(sourceFig, 'type', 'axes'); + targetAxes = findall(targetFig, 'type', 'axes'); + + % Ensure the number of axes match + if length(sourceAxes) ~= length(targetAxes) + error('Number of axes in source and target figures must be the same.'); + end + + % Loop through each pair of axes and copy styling properties + for i = 1:length(sourceAxes) + copyAxesProperties(sourceAxes(i), targetAxes(i)); + end + + % Apply general figure properties if desired + targetFig.Color = sourceFig.Color; % Background color +end + +function copyAxesProperties(sourceAx, targetAx) + % List of properties to copy from source to target axes + propsToCopy = {'XColor', 'YColor', 'ZColor', 'FontSize', 'FontName', ... + 'GridColor', 'GridLineStyle', 'MinorGridColor', 'Box', ... + 'XGrid', 'YGrid', 'ZGrid', 'XMinorGrid', 'YMinorGrid', 'ZMinorGrid', ... + 'LineWidth', 'TitleFontSizeMultiplier', 'LabelFontSizeMultiplier'}; + + % Copy properties from source to target + for i = 1:length(propsToCopy) + try + targetAx.(propsToCopy{i}) = sourceAx.(propsToCopy{i}); + catch + % Skip property if it doesn't exist or can't be copied + end + end + + % Copy axis labels and titles + targetAx.Title.String = sourceAx.Title.String; + targetAx.XLabel.String = sourceAx.XLabel.String; + targetAx.YLabel.String = sourceAx.YLabel.String; + targetAx.ZLabel.String = sourceAx.ZLabel.String; + + % Copy children elements like lines, patches, etc. + sourceChildren = allchild(sourceAx); + targetChildren = allchild(targetAx); + + % Ensure the number of children elements match + if length(sourceChildren) ~= length(targetChildren) + warning('Number of elements in source and target axes differ. Styling may not be applied completely.'); + end + + % Copy properties of children (like lines, patches, etc.), except colors and legends + for i = 1:min(length(sourceChildren), length(targetChildren)) + copyObjectProperties(sourceChildren(i), targetChildren(i)); + end +end + +function copyObjectProperties(sourceObj, targetObj) + % List of common properties to copy for plot elements (lines, patches, etc.) + propsToCopy = {'LineStyle', 'LineWidth', 'Marker', 'MarkerSize', ... + 'MarkerEdgeColor', 'MarkerFaceColor', 'DisplayName'}; + + % Copy properties from source to target, excluding colors + for i = 1:length(propsToCopy) + try + if ~contains(propsToCopy{i}, 'Color') % Skip color properties + targetObj.(propsToCopy{i}) = sourceObj.(propsToCopy{i}); + end + catch + % Skip property if it doesn't exist or can't be copied + end + end +end diff --git a/Functions/EQ_visuals/analyzeEQperformance.m b/Functions/EQ_visuals/analyzeEQperformance.m new file mode 100644 index 0000000..0325fe0 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Functions/EQ_visuals/showEQNoisePSD.m b/Functions/EQ_visuals/showEQNoisePSD.m new file mode 100644 index 0000000..51b12d7 --- /dev/null +++ b/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 diff --git a/Functions/EQ_visuals/showEQNoiseSNR.m b/Functions/EQ_visuals/showEQNoiseSNR.m new file mode 100644 index 0000000..8b517fd --- /dev/null +++ b/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 diff --git a/Functions/EQ_visuals/showEQTimeSignal.m b/Functions/EQ_visuals/showEQTimeSignal.m new file mode 100644 index 0000000..52d01b3 --- /dev/null +++ b/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 + diff --git a/Functions/EQ_visuals/showEQcoefficients.m b/Functions/EQ_visuals/showEQcoefficients.m new file mode 100644 index 0000000..3794212 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Functions/EQ_visuals/showErrorBurstCount.m b/Functions/EQ_visuals/showErrorBurstCount.m new file mode 100644 index 0000000..cb9e016 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Functions/EQ_visuals/showLevelConfusionMatrix.m b/Functions/EQ_visuals/showLevelConfusionMatrix.m new file mode 100644 index 0000000..61ea976 --- /dev/null +++ b/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 + diff --git a/Functions/EQ_visuals/showLevelHistogram.m b/Functions/EQ_visuals/showLevelHistogram.m new file mode 100644 index 0000000..40f31f4 --- /dev/null +++ b/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 + diff --git a/Functions/Metrics/air_garcia_implementation.m b/Functions/Metrics/air_garcia_implementation.m new file mode 100644 index 0000000..c71a26a --- /dev/null +++ b/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 + diff --git a/Functions/Metrics/calc_air.m b/Functions/Metrics/calc_air.m new file mode 100644 index 0000000..c42ab2f --- /dev/null +++ b/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) nice for looping parameter sets - - curFolder = pwd; - funcFolder=fileparts(mfilename('fullpath')); - if ~isempty(funcFolder) - cd(funcFolder); - end - - % TX - M = 4; - fsym = 112e9; - bitrate = 224e9; - apply_pulsef = 1; - fdac = 256e9; - fadc = 256e9; - random_key = 1; - db_precode = 0; - db_encode = 0; - rcalpha = 0.05; - kover = 16; - vbias_rel = 0.5; - u_pi = 2.9; - vbias = -vbias_rel*u_pi; - - laser_wavelength = 1310; - laser_linewidth = 0; - tx_bw_nyquist = 0.95; - - % Channel - link_length = 1; %km - - % RX - rop = -7; - rx_bw_nyquist = 0.9; - - % EQ - vnle_order=[50,7,7]; - dfe_order = [2 0 0]; - - len_tr = 4096*2; - mu_ffe = [0.0004 0.0004 0.0004]; - mu_dfe = 0.0004; - mu_dc = 0.00; - - % Replace optional input arguments if there are any - if ~isempty(varargin) - var_s = varargin{1}; - if isstruct(var_s) - fields = fieldnames(var_s); - for i = 1:numel(fields) - eval([fields{i}, ' = ', num2str( var_s.(fields{i}) ), ';']); - fprintf("%s <-- %.2f \n", fields{i}, var_s.(fields{i})); - end - else - error('Optional variables should be passed as a struct.'); - end - end - - fsym_ = floor( bitrate*1e-9./log2(M) ).*1e9; - if fsym_ ~= fsym - fsym = fsym_; - fprintf('Adapted symbolrate to %d GBd, to match provided bitrate of %d GBit/s using PAM %d \n',fsym.*1e-9,bitrate.*1e-9, M); - end - f_nyquist = fsym/2; - - %%%% TX Signal - Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"rrcalpha",rcalpha); - - [Digi_sig,Tx_symbols,Tx_bits] = PAMsource(... - "fsym",fsym,"M",M,"order",18,"useprbs",1,... - "fs_out",fdac,... - "applyclipping",0,"clipfactor",1.5,... - "applypulseform",apply_pulsef,"pulseformer",Pform,... - "randkey",random_key,... - "db_precode",db_precode,"db_encode",db_encode,... - "mrds_code",0,"mrds_blocklength",512).process(); - - Digi_sig.spectrum("displayname",'Digital Spectrum at TX','fignum',10,'normalizeTo0dB',1); - - %%%%% AWG - % El_sig = M8199A("kover",kover).process(Digi_sig); - El_sig = AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",1,"kover",kover,"bit_resolution",12).process(Digi_sig); - % El_sig.spectrum("displayname",'Digi Spectrum','fignum',100,'normalizeTo0dB',0); - El_sig = El_sig.setPower(0,"dBm"); - - El_sig.eye(fsym,M,"fignum",11,'displayname','Rx Eye'); - - %%%%% Low-pass el. components %%%%%% - El_sig = Filter('filtdegree',4,"f_cutoff",75e9,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(El_sig); - - %%%%% Electrical Driver Amplifier %%%%%% - El_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",3).process(El_sig); - - %%%%% 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); - - %%%%%% PD Square Law %%%%%% - Rx_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20,"nep",1.8e-11).process(Rx_sig); - - Rx_sig.eye(fsym,M,"fignum",30,'displayname','Rx Eye'); - - %%%%%% Low-pass RX (PD, El. Connectors and Scope %%%%%% - Rx_sig = Filter('filtdegree',4,"f_cutoff",75e9,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(Rx_sig); - - %%%%%% Low-pass inside Scope %%%%%% - Lp_scpe = Filter('filtdegree',4,"f_cutoff",110e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true); - - %%%%%% 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',1,'H_lpf',Lp_scpe).process(Rx_sig); - - Scpe_sig.spectrum("displayname",'Digital (256 GSa/s) Rx Spectrum','fignum',10,'normalizeTo0dB',1); - - Scpe_sig.normalize("mode","rms").plot("displayname",'Digital (256 GSa/s) Rx Spectrum','fignum',21,'clear',1); - - %%%%%% Sample to 2x fsym %%%%%% - Scpe_sig = Scpe_sig.resample("fs_in",fadc,"fs_out",2*fsym); - - %%%%%% Sync Rx signal with reference %%%%%% - [Scpe_sig,S] = Scpe_sig.tsynch("reference",Tx_symbols,"fs_ref",fsym); - - - %%% EQUALIZING - - %FFE or VNLE - [Eq_signal,Eq_noise] = EQ("Ne",vnle_order,"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).process(Scpe_sig,Tx_symbols); - - Rx_bits = PAMmapper(M,0).demap(Eq_signal); - - [~,numErrors,ber,~] = calc_ber(Rx_bits.signal,Tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - - Eq_signal.normalize("mode","rms").plot("displayname",'Digital (256 GSa/s) Rx Spectrum','fignum',21,'clear',0); - - %%% Visualize EQ Stuff - - figure(40); - clf - title(sprintf('PAM %d after EQ ; BER: %1.2e',M, ber )); - constellation = unique(Tx_symbols.signal); - received = NaN(numel(constellation),length(Tx_symbols)); - for lvl = 1:numel(constellation) - %Separate the equalized signal into the - %respective levels based on the actually - %transmitted level! - received(lvl,Tx_symbols.signal==constellation(lvl)) = Eq_signal.signal(Tx_symbols.signal==constellation(lvl)); - intermediate = received(lvl,:); - cnt(lvl) = numel(intermediate(~isnan(intermediate))); - hold on - histogram(received(lvl,:),1000,"EdgeAlpha",0,'DisplayName',['Lvl ',num2str(lvl),' | ',num2str(cnt(lvl)),' entries']); - end - legend - - fprintf('BER: %.2e \n',ber); - - autoArrangeFigures; - - if ~isempty(curFolder) - cd(curFolder); - end - -end \ No newline at end of file diff --git a/projects/A_minimal_example/run_loops.m b/projects/A_minimal_example/run_loops.m deleted file mode 100644 index 23245d1..0000000 --- a/projects/A_minimal_example/run_loops.m +++ /dev/null @@ -1,23 +0,0 @@ -uloops = struct; -uloops.bitrate = [300,330,360,390,420,450,480].*1e9; -uloops.M = [4,6,8]; - -wh = DataStorage(uloops); -wh.addStorage("ber"); - -fprintf("Let's simulate %d configuratiuons! \n",wh.getLastLinIndice); - -wh = submit_simulations(wh,"parallel",1); - -for m_ = uloops.M - ber_row = wh.getStoValue('ber', uloops.bitrate,m_); - - figure(2024) - hold on - plot(uloops.bitrate.*1e-9,ber_row,'DisplayName',sprintf('PAM %d',m_)); -end - -yline(3.8e-3,'LineWidth',2,'DisplayName','3.8e-3'); -yline(2e-2,'LineWidth',2,'LineStyle','--','DisplayName','2e-2'); -beautifyBERplot() -legend \ No newline at end of file diff --git a/projects/A_minimal_example/submit_simulations.m b/projects/A_minimal_example/submit_simulations.m deleted file mode 100644 index 32509c5..0000000 --- a/projects/A_minimal_example/submit_simulations.m +++ /dev/null @@ -1,72 +0,0 @@ -function wh = submit_simulations(wh,options) - -arguments - wh - options.parallel = 1; -end - -%%% 2) SUBMIT SIMULATION -% Initialize job results -if options.parallel - if isempty(gcp('nocreate')) - parpool; - end - results = parallel.FevalFuture.empty(); -else - results = []; -end - - -lin_idx = 1; -for lin_idx = 1:wh.getLastLinIndice - optionalVars = struct(); - if ~isempty(wh.getDimension) - % Build the optionalVars struct - [parametervalues,parameternames]=wh.getPhysIndicesByLinIndex(lin_idx); - for pidx = 1:numel(parameternames) - optionalVars.(parameternames{pidx}) = parametervalues{pidx}; - end - end - - %%% SIMULATION HERE - if options.parallel - - numOutputs = 1; - results(lin_idx) = parfeval(@imdd_simulation_minimal, numOutputs, optionalVars); - - else - - finalresults{lin_idx} = imdd_simulation_minimal(optionalVars); - wh.addValueToStorageByLinIdx(finalresults{lin_idx}, 'ber', lin_idx); - - end -end - - - -if options.parallel - - %%% 4) Setup waitbar - h = waitbar(0, 'Processing Simulations...'); - % Helper function to compute progress - updateWaitbar = @(~) waitbar(mean(arrayfun(@(f) strcmp(f.State, 'finished'), results)), h); - - fprintf('Fetching results... \n'); - - % Update the waitbar after each simulation - updateWaitbarFutures = afterEach(results, updateWaitbar, 0); - - % Close the waitbar after all simulations complete - afterAll(updateWaitbarFutures, @(~) delete(h), 0); - - %%% 7) Fetch final results after all computations - fetchOutputs(results); - - for ridx = 1:length(results) - wh.addValueToStorageByLinIdx(results(ridx).OutputArguments{1}, 'ber', ridx); - end - -end - - -end \ No newline at end of file diff --git a/projects/HighSpeedExperiment_2024/auswertung/analyze_mu_dc.m b/projects/HighSpeedExperiment_2024/auswertung/analyze_mu_dc.m new file mode 100644 index 0000000..ce35867 --- /dev/null +++ b/projects/HighSpeedExperiment_2024/auswertung/analyze_mu_dc.m @@ -0,0 +1,89 @@ + + +if 0 + uloops = struct; + uloops.precomp = [0,1]; + uloops.db_precode = [0,1]; + uloops.bitrate = [300,330,360,390,420,450,480].*1e9; %[300,330,360,390,420,450,480] + uloops.laser_wavelength = [1293,1302,1310,1318,1327.4]; + uloops.M = [4,6,8]; + uloops.link_length = [2]; % 1,2,3,5,6,8,10 + wh = DataStorage(uloops); + wh.addStorage("ber"); + wh = submit_simulations(wh,"parallel",1,"simulation_mode",0); + save('wh_2km',"wh"); +end + + + +for wavelength = wh.parameter.laser_wavelength.values + + cols = linspecer(6);%cbrewer2('Set2',10); + figcnt = 0; + figWidth = 21; % Full-width for IEEE double-column papers (~7 inches) + figHeight = 6; % Adjust height as needed (~3.5 inches) + figure('Units','centimeters','Position', [1 1 figWidth figHeight],'PaperUnits','centimeters','PaperPosition', [0 0 figWidth figHeight]) + tiledlayout(1,3, 'Padding', 'compact', 'TileSpacing', 'compact'); + for m = uloops.M + + + + figcnt = figcnt+1; + % subplot(1,3,figcnt) + nexttile + hold on + title(sprintf('%d km | %d nm | PAM %d',uloops.link_length,wavelength,m)); + + precomp = 1; + db_precode = 0; + a = wh.getStoValue('ber',precomp, db_precode, uloops.bitrate , wavelength, m, uloops.link_length); + ber_vnle = cellfun(@(x) x.vnle_pf_package{1,1}.ber_vnle, a); + ber_vnle_cell = cellfun(@(s) cellfun(@(p) p.ber_vnle, s.vnle_pf_package, 'UniformOutput', true), a, 'UniformOutput', false); + ber_vnle_best = cellfun(@(c) min(c), ber_vnle_cell); + plot(uloops.bitrate.*1e-9,ber_vnle_best,'DisplayName',sprintf('Tx precomp + VNLE'),'Color',cols(3,:),'LineStyle','-','HandleVisibility','on'); + xticks(uloops.bitrate.*1e-9); + xlim([min(uloops.bitrate.*1e-9) max(uloops.bitrate.*1e-9)]); + + precomp = 0; %0 + db_precode = 1; + a = wh.getStoValue('ber',precomp, db_precode, uloops.bitrate , wavelength, m, uloops.link_length); + % ber_db = cellfun(@(x) x.dbtgt_package{1,1}.ber, a); + ber_db_cell = cellfun(@(s) cellfun(@(p) p.ber, s.dbtgt_package, 'UniformOutput', true), a, 'UniformOutput', false); + ber_db_best = cellfun(@(c) min(c), ber_db_cell); + plot(uloops.bitrate.*1e-9,ber_db_best,'DisplayName',sprintf('DB tgt. + MLSE',uloops.link_length,uloops.M),'Color',cols(1,:),'LineStyle','-','HandleVisibility','on'); + xticks(uloops.bitrate.*1e-9); + xlim([min(uloops.bitrate.*1e-9) max(uloops.bitrate.*1e-9)]); + + precomp = 0; %0 + db_precode = 0; %1 + a = wh.getStoValue('ber',precomp, db_precode, uloops.bitrate , wavelength, m, uloops.link_length); + % ber_mlse = cellfun(@(x) x.vnle_pf_package{1,1}.ber_mlse, a); + ber_mlse_cell = cellfun(@(s) cellfun(@(p) p.ber_mlse, s.vnle_pf_package, 'UniformOutput', true), a, 'UniformOutput', false); + ber_mlse_best = cellfun(@(c) min(c), ber_mlse_cell); + plot(uloops.bitrate.*1e-9,ber_mlse_best,'DisplayName',sprintf('VNLE + 1 tap post-filter + MLSE',uloops.link_length,uloops.M),'Color',cols(4,:),'LineStyle','-','HandleVisibility','on'); + xticks(uloops.bitrate.*1e-9); + xlim([min(uloops.bitrate.*1e-9) max(uloops.bitrate.*1e-9)]); + + set(gca, 'YScale', 'log'); + ylim([8e-5 0.3]); + yline([4.8e-3, 2e-2],'HandleVisibility','off','LineWidth',1,'LineStyle','--','Color',[0.1 0.1 0.1]); + % legend + beautifyBERplot() + xlabel('Bit Rate in Gbps'); + ylabel('BER'); + + if m ==4 + text(310,6.8e-3,"4.8e-3","FontSize",10,"Interpreter","latex") + text(310,3e-2,"2e-2","FontSize",10,"Interpreter","latex") + end + + % text(0.5,1,sprintf('%d km %d nm PAM %d',uloops.link_length,wavelength,m),... + % 'Units', 'normalized',"FontSize",10,"Interpreter","latex","BackgroundColor",[1 1 1],"EdgeColor",[0 0 0],'HorizontalAlignment','center','VerticalAlignment','top') + end + + lgd = legend; + % Place the legend underneath the tiled layout + lgd.NumColumns = 3; + lgd.Layout.Tile = 'south'; + +end \ No newline at end of file diff --git a/projects/IMDD_base_system/imdd_it.m b/projects/IMDD_base_system/imdd_it.m index b69a8c0..d0f4084 100644 --- a/projects/IMDD_base_system/imdd_it.m +++ b/projects/IMDD_base_system/imdd_it.m @@ -14,7 +14,7 @@ if 1 wh = DataStorage(uloops); wh.addStorage("ber"); - wh = submit_simulations(wh,"parallel",1,"simulation_mode",0); + wh = submit_simulations(wh,"parallel",0,"simulation_mode",1); end wh_ana = wh;