From 005e821131ff8007b9e0907cf588abdd4bf6bcd2 Mon Sep 17 00:00:00 2001 From: magf Date: Mon, 2 Feb 2026 10:56:12 +0100 Subject: [PATCH] Added: - Class folder 'Timing Recovery' with different timing recoveries - Minimal example for the timing recovery on the FSO data - New evaluation scripts in the FSO project folder --- .../Timing Recovery/Godard_Timing_Recovery.m | 112 +++++++ .../Timing Recovery/MaxVar_Timing_Recovery.m | 111 +++++++ Classes/04_DSP/Timing Recovery/Time_Shifter.m | 39 +++ .../04_DSP/Timing Recovery/Timing_Recovery.m | 49 +++ .../Timing Recovery/Timing_Recovery_GPT.m | 77 +++++ .../Timing Recovery/Timing_Recovery_Move_It.m | 77 +++++ .../Evaluation Scripts/evalscript_FFE_DFE.m | 36 +++ .../Evaluation Scripts/evalscript_ffe.m | 57 ++++ ...ript_optimize_timing_recovery_parameters.m | 68 ++++ .../Evaluation Scripts/evalscript_pam_2.m | 66 ++++ .../Evaluation Scripts/evalscript_pam_4.m | 61 ++++ .../evalscript_pam_4_baud_rate_sweep.m | 77 +++++ .../evalscript_timing_shift.m | 9 + .../Evaluation Scripts/first_analysis_2.m | 190 ++++++++++++ .../first_analysis_baud_rate_sweep.m | 202 ++++++++++++ .../Evaluation Scripts/first_analysis_ber.m | 207 +++++++++++++ .../Evaluation Scripts/first_analysis_mf.m | 293 ++++++++++++++++++ .../first_analysis_optimize.m | 177 +++++++++++ .../first_analysis_time_rec.m | 160 ++++++++++ .../first_analysis_time_rec_pam2.m | 171 ++++++++++ .../first_analysis_time_rec_pam4.m | 189 +++++++++++ .../first_analysis_time_shift.m | 221 +++++++++++++ .../Evaluation Scripts/first_analysis_tr_mf.m | 271 ++++++++++++++++ .../first_analysis_tr_mf_pam2.m | 195 ++++++++++++ .../first_analysis_tr_mf_pam_2.m | 272 ++++++++++++++++ .../ml_mlse_complexity_plot.m | 11 + .../plotscript_baud_rate_sweep_method.m | 36 +++ .../FSO_timing_recovery_minimal_example.m | 251 +++++++++++++++ 28 files changed, 3685 insertions(+) create mode 100644 Classes/04_DSP/Timing Recovery/Godard_Timing_Recovery.m create mode 100644 Classes/04_DSP/Timing Recovery/MaxVar_Timing_Recovery.m create mode 100644 Classes/04_DSP/Timing Recovery/Time_Shifter.m create mode 100644 Classes/04_DSP/Timing Recovery/Timing_Recovery.m create mode 100644 Classes/04_DSP/Timing Recovery/Timing_Recovery_GPT.m create mode 100644 Classes/04_DSP/Timing Recovery/Timing_Recovery_Move_It.m create mode 100644 projects/FSO_transmission/Evaluation Scripts/evalscript_FFE_DFE.m create mode 100644 projects/FSO_transmission/Evaluation Scripts/evalscript_ffe.m create mode 100644 projects/FSO_transmission/Evaluation Scripts/evalscript_optimize_timing_recovery_parameters.m create mode 100644 projects/FSO_transmission/Evaluation Scripts/evalscript_pam_2.m create mode 100644 projects/FSO_transmission/Evaluation Scripts/evalscript_pam_4.m create mode 100644 projects/FSO_transmission/Evaluation Scripts/evalscript_pam_4_baud_rate_sweep.m create mode 100644 projects/FSO_transmission/Evaluation Scripts/evalscript_timing_shift.m create mode 100644 projects/FSO_transmission/Evaluation Scripts/first_analysis_2.m create mode 100644 projects/FSO_transmission/Evaluation Scripts/first_analysis_baud_rate_sweep.m create mode 100644 projects/FSO_transmission/Evaluation Scripts/first_analysis_ber.m create mode 100644 projects/FSO_transmission/Evaluation Scripts/first_analysis_mf.m create mode 100644 projects/FSO_transmission/Evaluation Scripts/first_analysis_optimize.m create mode 100644 projects/FSO_transmission/Evaluation Scripts/first_analysis_time_rec.m create mode 100644 projects/FSO_transmission/Evaluation Scripts/first_analysis_time_rec_pam2.m create mode 100644 projects/FSO_transmission/Evaluation Scripts/first_analysis_time_rec_pam4.m create mode 100644 projects/FSO_transmission/Evaluation Scripts/first_analysis_time_shift.m create mode 100644 projects/FSO_transmission/Evaluation Scripts/first_analysis_tr_mf.m create mode 100644 projects/FSO_transmission/Evaluation Scripts/first_analysis_tr_mf_pam2.m create mode 100644 projects/FSO_transmission/Evaluation Scripts/first_analysis_tr_mf_pam_2.m create mode 100644 projects/FSO_transmission/Evaluation Scripts/ml_mlse_complexity_plot.m create mode 100644 projects/FSO_transmission/Evaluation Scripts/plotscript_baud_rate_sweep_method.m create mode 100644 projects/FSO_transmission/FSO_timing_recovery_minimal_example.m diff --git a/Classes/04_DSP/Timing Recovery/Godard_Timing_Recovery.m b/Classes/04_DSP/Timing Recovery/Godard_Timing_Recovery.m new file mode 100644 index 0000000..12eeeaf --- /dev/null +++ b/Classes/04_DSP/Timing Recovery/Godard_Timing_Recovery.m @@ -0,0 +1,112 @@ +classdef Godard_Timing_Recovery < handle + + properties(Access=public) + mode + num_blocks + fft_length + sps + rolloff + mu + Ki + end + + methods(Access=public) + function obj = Godard_Timing_Recovery(options) + arguments(Input) + + options.mode = 0; + options.num_blocks = 1; + options.fft_length = 1024; + options.sps = 2; + options.rolloff = 0; + options.mu = 0; + options.Ki = 1e-3; + + end + + fn = fieldnames(options); + for n = 1:numel(fn) + obj.(fn{n}) = options.(fn{n}); + end + + %obj-Initialization here% + + end + + function [data_out, tau_hat] = process(obj, data_in) + + data_out = data_in; + output_vector = zeros(length(data_in),1); + block_length = length(data_in)/obj.num_blocks; + for i = 1:obj.num_blocks + i_start = ((i-1)*block_length)+1; + i_end = i*block_length; + x = data_in.signal(i_start:i_end); + beta = obj.rolloff; + eta = obj.sps; + R = fft(x,obj.fft_length); + R_full = fft(x); + + if obj.mode == 0 % Classic Godard + % Calculate time shift + k0 = (1:obj.fft_length/2).'; + idx1 = k0; + idx2 = k0+(obj.fft_length/2); + tau_hat = sum(imag(R(idx1) .* conj(R(idx2)))); + elseif obj.mode == 1 || obj.mode == 2 % Modified Godard 1 + % Calculate Shift for the received signal + shiftBins = (1 - 1/eta) * obj.fft_length; + if abs(shiftBins - round(shiftBins)) > 1e-12 + disp('Warning: shiftBins=(1-1/eta)*fft_length is non-integer. Choose compatible values for fft_length and eta.'); + end + shiftBins = round(shiftBins); + + % Calculate upper and lower bounds + kStart = ((1-beta)/(2*eta)) * obj.fft_length; + kEnd = ((1+beta)/(2*eta)) * obj.fft_length - 1; + if abs(kStart - round(kStart)) > 1e-12 || abs(kEnd - round(kEnd)) > 1e-12 + disp('Warning: kStart/kEnd are non-integer. Choose compatible values for fft_length, eta, and beta.'); + end + kStart = round(kStart); + kEnd = round(kEnd); + + k0 = (kStart:kEnd).'; + idx1 = k0 + 1; + idx2 = mod(k0 + shiftBins, obj.fft_length) + 1; + + % Calculate time shift + if obj.mode == 1 + tau_hat = sum(imag(R(idx1) .* conj(R(idx2)))); + elseif obj.mode == 2 + tau_hat = sum(angle(R(idx1))-angle(R(idx2))); + end + end + + + % % Normalization + % denom = floor(log10(abs(tau_hat))); + % tau_hat = tau_hat / 10^denom; + + % Calculate mu using a first-order loop filter + mu_block = obj.mu + obj.Ki * tau_hat; + + % % Shifting the signal in time domain using interpolation + % x_original = linspace(0,length(x)-1,length(x)).'; + % x_new = x_original + mu_block; + % output_vector(i_start:i_end) = interp1(x_original,x,x_new,'linear','extrap'); + + % Shifting the signal in frequency domain + k = (0:block_length-1).'; + phaseRamp = exp(-1j * 2*pi * (k/block_length) * mu_block); + R_shifted = R_full .* phaseRamp; + + output_vector(i_start:i_end) = real(ifft(R_shifted)); + + end + data_out.signal = output_vector; + end + + + end +end + diff --git a/Classes/04_DSP/Timing Recovery/MaxVar_Timing_Recovery.m b/Classes/04_DSP/Timing Recovery/MaxVar_Timing_Recovery.m new file mode 100644 index 0000000..d6517d5 --- /dev/null +++ b/Classes/04_DSP/Timing Recovery/MaxVar_Timing_Recovery.m @@ -0,0 +1,111 @@ +classdef MaxVar_Timing_Recovery < handle + + properties(Access=public) + mode + sps + fsym + fadc + num_tau + comp_signal + comp_mode + end + + methods(Access=public) + function obj = MaxVar_Timing_Recovery(options) + arguments(Input) + + options.mode = 0; + options.sps = 2; + options.fsym = 32e9; + options.fadc = 80e9; + options.num_tau = 64; + options.comp_signal = 0; + options.comp_mode = 0; + end + + fn = fieldnames(options); + for n = 1:numel(fn) + obj.(fn{n}) = options.(fn{n}); + end + + %obj-Initialization here% + + end + + function [data_out] = process(obj, data_in) + + data_out = data_in; + x = data_in.signal; + + if obj.comp_mode + our_signal = data_in; + their_signal = obj.comp_signal; + end + + if obj.mode == 0 + vars = zeros(obj.sps,1); + for phi = 1:obj.sps + + % Test variance for different start samples + y_phi = x(phi:obj.sps:end); + vars(phi) = var(y_phi); + + % Comparison to reference signal + if obj.comp_mode + our_signal.signal = x(phi:obj.sps:end); + our_signal.fs = 6e9; + our_signal.normalize("mode","rms").plot("displayname",['Our signal, ' num2str(vars(phi))],'fignum',1231+phi); + their_signal.normalize("mode","rms").plot("displayname",'Their signal','fignum',1231+phi); + end + + end + + % Choose signal configuration with the maximum variance + [~,phi_opt] = max(vars); + + y = x(phi_opt:obj.sps:end).'; + + elseif obj.mode == 1 + % Create a grid with different (sub)sample starting points + T = 1/obj.fsym; + t = (0:length(x)-1)/obj.fadc; + tauGrid = linspace(0, T, obj.num_tau); + + % Interpolate the signal starting from every defined point + % and calculate the MMSE/variance + vars = zeros(size(tauGrid)); + for i = 1:numel(tauGrid) + tau = tauGrid(i); + % tk = linspace(tau, t(end), length(x)/obj.sps); + tk = tau : T : t(end)+tau; + xk = interp1(t, x, tk, 'linear', 'extrap'); + if obj.comp_mode % MMSE + e = xk.' - obj.comp_signal.signal; + vars(i) = mean(abs(e).^2); + else % Variance + vars(i) = var(xk, 1); + end + end + + if obj.comp_mode % MMSE + [~,idx] = min(vars); + else % Variance + [~,idx] = max(vars); + end + tau_opt = tauGrid(idx); + + % Choosing the signal at optimum + % tk = linspace(tau_opt, t(end), length(x)/obj.sps); + tk = tau_opt : T : t(end)+tau_opt; + y = interp1(t, x, tk, 'linear', 'extrap'); + end + % if length(y) ~= length(x)/obj.sps + % y = [y 0]; + % end + data_out.signal = y.'; + end + + + end +end + diff --git a/Classes/04_DSP/Timing Recovery/Time_Shifter.m b/Classes/04_DSP/Timing Recovery/Time_Shifter.m new file mode 100644 index 0000000..f08e3fe --- /dev/null +++ b/Classes/04_DSP/Timing Recovery/Time_Shifter.m @@ -0,0 +1,39 @@ +classdef Time_Shifter < handle + + properties(Access=public) + value + end + + methods(Access=public) + function obj = Time_Shifter(options) + arguments(Input) + + options.value = 0; + + end + + fn = fieldnames(options); + for n = 1:numel(fn) + obj.(fn{n}) = options.(fn{n}); + end + + %obj-Initialization here% + + end + + function [data_out, tau_hat] = process(obj, data_in) + + data_out = data_in; + x = data_in.signal; + R = fft(x); + k = (0:length(R)-1).'; + phaseRamp = exp(-1j * 2*pi * (k/length(R)) * obj.value); + R_shifted = R .* phaseRamp; + data_out.signal = real(ifft(R_shifted)); + + end + + + end +end + diff --git a/Classes/04_DSP/Timing Recovery/Timing_Recovery.m b/Classes/04_DSP/Timing Recovery/Timing_Recovery.m new file mode 100644 index 0000000..3999ad7 --- /dev/null +++ b/Classes/04_DSP/Timing Recovery/Timing_Recovery.m @@ -0,0 +1,49 @@ +classdef Timing_Recovery < handle + + properties(Access=public) + modulation + timing_error_detector + sps + damping_factor + normalized_loop_bandwidth + detector_gain + end + + methods(Access=public) + function obj = Timing_Recovery(options) + arguments(Input) + + options.modulation = 'PAM/PSK/QAM'; + options.timing_error_detector = 'Gardner (non-data-aided)'; + options.sps = 2; + options.damping_factor = 1.0; + options.normalized_loop_bandwidth = 0.01; + options.detector_gain = 2.7; + + end + + fn = fieldnames(options); + for n = 1:numel(fn) + obj.(fn{n}) = options.(fn{n}); + end + + %obj-Initialization here% + + end + + function [data_out, timing_error] = process(obj, data_in) + + timing_synchronization = comm.SymbolSynchronizer( ... + "Modulation", obj.modulation,... + "TimingErrorDetector", obj.timing_error_detector, ... + "SamplesPerSymbol", obj.sps, ... + "DampingFactor", obj.damping_factor, ... + "NormalizedLoopBandwidth", obj.normalized_loop_bandwidth, ... + "DetectorGain", obj.detector_gain); + data_out = data_in; + [data_out.signal, timing_error] = timing_synchronization(data_in.signal); + + end + end +end + diff --git a/Classes/04_DSP/Timing Recovery/Timing_Recovery_GPT.m b/Classes/04_DSP/Timing Recovery/Timing_Recovery_GPT.m new file mode 100644 index 0000000..33d0509 --- /dev/null +++ b/Classes/04_DSP/Timing Recovery/Timing_Recovery_GPT.m @@ -0,0 +1,77 @@ +classdef Timing_Recovery_GPT < handle + + properties(Access=public) + sps + muGrid + end + + methods(Access=public) + function obj = Timing_Recovery_GPT(options) + arguments(Input) + + options.sps = 2; + options.muGrid = 0; + + end + + fn = fieldnames(options); + for n = 1:numel(fn) + obj.(fn{n}) = options.(fn{n}); + end + + %obj-Initialization here% + + end + + function [data_out, mu_best, score] = process(obj, data_in) + %MAXVARTIMINGSYNC Choose sampling phase mu that maximizes variance of downsampled symbols. + % + % data_in : matched-filtered samples (complex or real), length N + % sps : samples per symbol (here typically 2) + % muGrid : candidate fractional offsets in [0,1) + % + % data_out : symbol-rate samples (length floor(N/sps)) + % mu_best: chosen fractional offset + % score : variance score for each mu in muGrid + + data_out = data_in; + + x = data_in.signal(:); + N = length(x); + Ns = floor(N/obj.sps); + + if nargin < 3 || isempty(obj.muGrid) + obj.muGrid = linspace(0, 0.99, 101); % 0..0.99 in ~0.01 steps + end + + % Symbol indices (1-based sample positions) + n0 = 1; % start sample index + k = (0:Ns-1).'; + tBase = n0 + k*obj.sps; % integer times (1, 1+sps, ...) + + score = zeros(numel(obj.muGrid),1); + + for m = 1:numel(obj.muGrid) + mu = obj.muGrid(m); + t = tBase + mu; + + % Linear fractional sampling + y = interp1(1:N, x, t, 'linear', 'extrap'); + + % For PAM, maximize variance of real part (or abs if you prefer) + yr = real(y); + score(m) = var(yr, 1); % use population variance (normalization doesn't matter for argmax) + end + + % Pick best mu + [~, idx] = max(score); + mu_best = obj.muGrid(idx); + + % Resample with best mu + t = tBase + mu_best; + data_out.signal = interp1(1:N, x, t, 'linear', 'extrap'); + + end + end +end + diff --git a/Classes/04_DSP/Timing Recovery/Timing_Recovery_Move_It.m b/Classes/04_DSP/Timing Recovery/Timing_Recovery_Move_It.m new file mode 100644 index 0000000..bbd0578 --- /dev/null +++ b/Classes/04_DSP/Timing Recovery/Timing_Recovery_Move_It.m @@ -0,0 +1,77 @@ +classdef Timing_Recovery_Move_It < handle + + properties(Access=public) + f_sim + gamma + end + + methods(Access=public) + function obj = Timing_Recovery_Move_It(options) + arguments(Input) + + options.f_sim = 14e9; + options.gamma = 0.1; + + end + + fn = fieldnames(options); + for n = 1:numel(fn) + obj.(fn{n}) = options.(fn{n}); + end + + %obj-Initialization here% + + end + + function data_out = process(obj, data_in) + + e = NaN(size(data_in.signal)); + T = 1/obj.f_sim; + t = zeros(size(data_in.signal)); + t(2) = T/2; + mu = zeros(size(data_in.signal)); + +% n = 2; + for k = 3:2:length(data_in.signal) + data_in.signal(k-1) = data_in.signal(k-1)*(1-mu(k-2)) + data_in.signal(k)*mu(k-2); + data_in.signal(k) = data_in.signal(k)*(1-mu(k-1)) + data_in.signal(k+1)*mu(k-1); + + % TED + e(k) = (data_in.signal(k-2)-data_in.signal(k))*data_in.signal(k-1); + e(k+1) = e(k); + + % TED Mueller Mueller +% e(k) = ref_in(n-1)*data_in.signal(k) - ref_in(n)*data_in.signal(k-2); +% e(k+1) = e(k); +% n = n+1; +% + % interpolator control +% t(k) = t(k-1) + T/2 + obj.gamma/2*e(k); +% t(k+1) = t(k) + T/2 + obj.gamma/2*e(k+1); + t(k) = t(k-1) + T/2 + obj.gamma*e(k)*T/2; + t(k+1) = t(k) + T/2 + obj.gamma*e(k+1)*T/2; +% t(k) = k*T/2 + obj.gamma*e(k)*T/2; +% t(k+1) = k*T/2 + obj.gamma*e(k+1)*T/2; + + % interpolator + mu(k) = t(k)/(T/2) - round(t(k)/(T/2)); + mu(k+1) = t(k+1)/(T/2) - round(t(k+1)/(T/2)); + + thres = 0.7; + if mu(k)-mu(k-1) > thres + mu(k) = mu(k) - 1; + elseif mu(k) - mu(k-1) < -thres + mu(k) = mu(k) + 1; + end + + if mu(k+1)-mu(k) > thres + mu(k+1) = mu(k+1) - 1; + elseif mu(k+1) - mu(k) < -thres + mu(k+1) = mu(k+1) + 1; + end + end + data_out = data_in; + end + end +end + diff --git a/projects/FSO_transmission/Evaluation Scripts/evalscript_FFE_DFE.m b/projects/FSO_transmission/Evaluation Scripts/evalscript_FFE_DFE.m new file mode 100644 index 0000000..29069ce --- /dev/null +++ b/projects/FSO_transmission/Evaluation Scripts/evalscript_FFE_DFE.m @@ -0,0 +1,36 @@ +taps_ffe = [200, 3, 2]; +taps_dfe = [5, 2, 1]; +trlen = 4096*2; +trloops = 5; +ddloops = 5; +K_FFE = 1; +DCmu = 0.005; +mu_ffe_values = [0.001, 0.0005, 0.001]; +mu_dfe_values = 0.0007; +DDmu = [mu_ffe_values, mu_dfe_values]; +DFEmu = 0.004; +FFEmu = 0; +plot_final = 0; +idealdfe = 1; +num_pf_coeffs = 1; +damp_factor = 1; +norm_loop_bw = 0.01; +det_gain = 2.7; + +step_size = 0.0001; + +BER = []; +taps_ffe_and_dfe = zeros(1,6); +for i = 1:6 + for j = 1:200 + taps_ffe_and_dfe(i) = j; + BER_run = first_analysis_2(taps_ffe_and_dfe(1:3), taps_ffe_and_dfe(4:6), trlen, trloops, ddloops, ... + K_FFE, DCmu, DDmu, DFEmu, FFEmu, plot_final, idealdfe, num_pf_coeffs, ... + damp_factor, norm_loop_bw, det_gain); + BER = [BER, BER_run]; + close all + end + [~,index_minimum_ber] = min(BER); + taps_ffe_and_dfe(i) = index_minimum_ber; +end +disp(taps_ffe_and_dfe) \ No newline at end of file diff --git a/projects/FSO_transmission/Evaluation Scripts/evalscript_ffe.m b/projects/FSO_transmission/Evaluation Scripts/evalscript_ffe.m new file mode 100644 index 0000000..8b55700 --- /dev/null +++ b/projects/FSO_transmission/Evaluation Scripts/evalscript_ffe.m @@ -0,0 +1,57 @@ +%% +ffe_start = 0; +ffe_step = 5; +ffe_end = 30; +ffe_first_order = ffe_start:ffe_step:ffe_end; + +current = 225:10:285; +current = string(current); +power = [28.02, 33.2, 37.5, 42.3, 46.3, 49.3, 53.4]; +power = string(power); +num_pf_coeff = 4; +% taps_ffe = [300, 0, 0]; +taps_dfe = [0, 0, 0]; +M = 4; +trlength = 4096*4; +x = 225:10:285; +BER_PAM_4 = []; + +for i = 1:length(ffe_first_order) + for eq_method = 2 + for j = 4 + BER_run = first_analysis_ber(current(j), power(j), num_pf_coeff, [300, 15, ffe_first_order(i)], taps_dfe, M, trlength, eq_method); + BER_PAM_4 = [BER_PAM_4, BER_run]; + end + end +end +save('C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\Silas DSP\imdd_simulation\projects\FSO_transmission\BER_PAM_4_FFE_Second_Order_Tap_Sweep.mat', 'ffe_first_order', 'BER_PAM_4') +BER_PAM_4 = []; + +%% +BER = load('C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\Silas DSP\imdd_simulation\projects\FSO_transmission\BER_PAM_4_FFE_Second_Order_Tap_Sweep.mat'); +BER = BER.BER_PAM_4; + +figure(202120) +plot(ffe_first_order, BER, '-o','LineWidth',1.75); +hold on + +h1 = yline(2e-2, ':k', 'LineWidth',1.5); +h2 = yline(3.8e-3,':b', 'LineWidth',1.5); +h3 = yline(4.85e-3,':g', 'LineWidth',1.5); +h4 = yline(2.2e-4,':r', 'LineWidth',1.5); + +% FEC Labels direkt im Plot +text(286,2.2e-2,'o-FEC','Color','k','FontSize', 14, 'Interpreter','latex') +text(285,3.3e-3,'HD-FEC','Color','b','FontSize', 14, 'Interpreter','latex') +text(282.5,5.4e-3,'KP4+Hamming','Color','g','FontSize', 14,'Interpreter','latex') +text(287,2.4e-4,'KP4','Color','r','FontSize', 14,'Interpreter','latex') + +xlabel('Number of DFE First Order Taps', 'Interpreter','latex') +ylabel('BER', 'Interpreter','latex') +% title('BER for PAM-4', 'Interpreter','latex') + +grid minor +% ylim([1e-4 5e-1]) +set(gca,'YScale','log') +% beautifyBERplot +hold off \ No newline at end of file diff --git a/projects/FSO_transmission/Evaluation Scripts/evalscript_optimize_timing_recovery_parameters.m b/projects/FSO_transmission/Evaluation Scripts/evalscript_optimize_timing_recovery_parameters.m new file mode 100644 index 0000000..f867c5a --- /dev/null +++ b/projects/FSO_transmission/Evaluation Scripts/evalscript_optimize_timing_recovery_parameters.m @@ -0,0 +1,68 @@ +%% Fixe Simulationsparameter +baud_rate = 6; +power = 42.3; +num_pf_coeff = 4; +M = 4; +trlength = 4096*2; + +%% === Suchräume für Tap-Werte (selbst bestimmbar) ================= +% WICHTIG: wähle Grenzen passend zur Skalierung deiner Taps. +% Beispiel: FFE Taps vielleicht grob im Bereich [-500 .. 500], +% DFE Taps eher kleiner, z.B. [-50 .. 50] – anpassen! +ffeRange = [0, 500]; +dfeRange = [0, 200]; + +vars = [ + optimizableVariable('ffe1', ffeRange) + optimizableVariable('ffe2', ffeRange) + optimizableVariable('ffe3', ffeRange) + + optimizableVariable('dfe1', dfeRange) + optimizableVariable('dfe2', dfeRange) + optimizableVariable('dfe3', dfeRange) +]; + +%% === Objective =================================================== +objectiveFcn = @(x) objectiveWrapperTaps(x, ... + baud_rate, power, num_pf_coeff, M, trlength); + +%% === bayesopt ==================================================== +results = bayesopt(objectiveFcn, vars, ... + 'AcquisitionFunctionName','expected-improvement-plus', ... + 'MaxObjectiveEvaluations', 60, ... + 'IsObjectiveDeterministic', false, ... + 'ExplorationRatio', 0.5, ... + 'Verbose', 1, ... + 'PlotFcn', {@plotObjectiveModel,@plotMinObjective}); + +%% === Beste Lösung ausgeben ====================================== +bestX = results.XAtMinObjective; +bestBER = results.MinObjective; + +bestFFE = [bestX.ffe1, bestX.ffe2, bestX.ffe3]; +bestDFE = [bestX.dfe1, bestX.dfe2, bestX.dfe3]; + +fprintf('\n=== Bestes Ergebnis ===\n'); +fprintf('BER : %.3e\n', bestBER); +fprintf('taps_ffe : [%g %g %g]\n', bestFFE); +fprintf('taps_dfe : [%g %g %g]\n\n', bestDFE); + +%% ================================================================= +function ber = objectiveWrapperTaps(x, baud_rate, power, num_pf_coeff, M, trlength) + taps_ffe = double([x.ffe1, x.ffe2, x.ffe3]); + taps_dfe = double([x.dfe1, x.dfe2, x.dfe3]); + + % Optional: falls du z.B. Monotonicität / Struktur erzwingen willst, hier. + % taps_dfe(1) = 0; % Beispiel: ersten DFE-Tap fixieren + + try + ber = first_analysis_optimize(baud_rate, power, num_pf_coeff, taps_ffe, taps_dfe, M, trlength); + + if isempty(ber) || ~isscalar(ber) || ~isfinite(ber) || ber < 0 + ber = 1; % Penalty + end + catch ME + warning("Objective failed: %s", ME.message); + ber = 1; % Penalty + end +end diff --git a/projects/FSO_transmission/Evaluation Scripts/evalscript_pam_2.m b/projects/FSO_transmission/Evaluation Scripts/evalscript_pam_2.m new file mode 100644 index 0000000..b0310bd --- /dev/null +++ b/projects/FSO_transmission/Evaluation Scripts/evalscript_pam_2.m @@ -0,0 +1,66 @@ +%% +current = 225:10:285; +current = string(current); +power = [28.02, 33.2, 37.5, 42.3, 46.3, 49.3, 53.4]; +power = string(power); +num_pf_coeff = 4; +taps_ffe = [200, 0, 0]; +taps_dfe = [0, 0, 0]; +M = 2; +trlength = 4096*2; +x = 225:10:285; +BER_PAM_2 = []; +post_only = 0; + +for eq_method = 2:4 + + if ~post_only + for j = 1:length(current) + BER_run = first_analysis_ber(current(j), power(j), num_pf_coeff, taps_ffe, taps_dfe, M, trlength, eq_method); + BER_PAM_2 = [BER_PAM_2, BER_run]; + end + + save('C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\Silas DSP\imdd_simulation\projects\FSO_transmission\BER_PAM_2_' + string(eq_method) + '.mat', 'x', 'BER_PAM_2') + end + +end + +%% +x = 225:10:285; +for k = 1:4 + BER = load('C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\Silas DSP\imdd_simulation\projects\FSO_transmission\BER_PAM_2_' + string(k) + '.mat'); + BER = BER.BER_PAM_2; + + figure(202120) + plot(x, BER, '-o','LineWidth',1.75); + hold on +end +old_BER = [-1.5, -1.75, -2, -2.3, -2.6, -2.25, -1.95]; +old_BER = 10.^(old_BER); +plot(x, old_BER, '-o','LineWidth',1.75) + +h1 = yline(2e-2, ':k', 'LineWidth',1.5); +h2 = yline(3.8e-3,':b', 'LineWidth',1.5); +h3 = yline(4.85e-3,':g', 'LineWidth',1.5); +h4 = yline(2.2e-4,':r', 'LineWidth',1.5); + +% Legende NUR für Kurven +legend('FFE', 'FFE+PF+MLSE', 'DB', 'ML-MLSE', 'BER Paper', ... + 'Interpreter','latex', ... + 'Location','southwest', 'FontSize', 14) + +% FEC Labels direkt im Plot +text(280,2.2e-2,'o-FEC','Color','k','FontSize', 14, 'Interpreter','latex') +text(280,3.3e-3,'HD-FEC','Color','b','FontSize', 14, 'Interpreter','latex') +text(280,5.4e-3,'KP4+Hamming','Color','g','FontSize', 14,'Interpreter','latex') +text(280,2.4e-4,'KP4','Color','r','FontSize', 14,'Interpreter','latex') + +xlabel('Laser Bias Current [mA]', 'Interpreter','latex') +ylabel('BER', 'Interpreter','latex') +% title('BER for PAM-2', 'Interpreter','latex') + +grid minor +ylim([1e-4 5e-1]) +set(gca,'YScale','log') +% beautifyBERplot +hold off \ No newline at end of file diff --git a/projects/FSO_transmission/Evaluation Scripts/evalscript_pam_4.m b/projects/FSO_transmission/Evaluation Scripts/evalscript_pam_4.m new file mode 100644 index 0000000..cf3b301 --- /dev/null +++ b/projects/FSO_transmission/Evaluation Scripts/evalscript_pam_4.m @@ -0,0 +1,61 @@ +%% +current = 225:10:285; +current = string(current); +power = [28.02, 33.2, 37.5, 42.3, 46.3, 49.3, 53.4]; +power = string(power); +num_pf_coeff = 4; +taps_ffe = [300, 0, 0]; +taps_dfe = [5, 0, 0]; +M = 4; +trlength = 4096*4; +x = 225:10:285; +BER_PAM_4 = []; + +for eq_method = 2 + for j = 4 + BER_run = first_analysis_ber(current(j), power(j), num_pf_coeff, taps_ffe, taps_dfe, M, trlength, eq_method); + BER_PAM_4 = [BER_PAM_4, BER_run]; + end + save('C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\Silas DSP\imdd_simulation\projects\FSO_transmission\BER_PAM_4_DFE_' + string(eq_method) + '.mat', 'x', 'BER_PAM_4') + BER_PAM_4 = []; +end + +%% +x = 225:10:285; +for k = 1:4 + BER = load('C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\Silas DSP\imdd_simulation\projects\FSO_transmission\BER_PAM_4_' + string(k) + '.mat'); + BER = BER.BER_PAM_4; + + figure(202120) + plot(x, BER, '-o','LineWidth',1.75); + hold on +end +old_BER = [-1.6, -1.85, -2.2, -2.45, -2.3, -2, -1.4]; +old_BER = 10.^(old_BER); +plot(x, old_BER, '-o','LineWidth',1.75) + +h1 = yline(2e-2, ':k', 'LineWidth',1.5); +h2 = yline(3.8e-3,':b', 'LineWidth',1.5); +h3 = yline(4.85e-3,':g', 'LineWidth',1.5); +h4 = yline(2.2e-4,':r', 'LineWidth',1.5); + +% Legende NUR für Kurven +legend('FFE', 'FFE+PF+MLSE', 'DB', 'ML-MLSE', 'BER Paper', ... + 'Interpreter','latex', ... + 'Location','southwest', 'FontSize', 14) + +% FEC Labels direkt im Plot +text(286,2.2e-2,'o-FEC','Color','k','FontSize', 14, 'Interpreter','latex') +text(285,3.3e-3,'HD-FEC','Color','b','FontSize', 14, 'Interpreter','latex') +text(282.5,5.4e-3,'KP4+Hamming','Color','g','FontSize', 14,'Interpreter','latex') +text(287,2.4e-4,'KP4','Color','r','FontSize', 14,'Interpreter','latex') + +xlabel('Laser Bias Current [mA]', 'Interpreter','latex') +ylabel('BER', 'Interpreter','latex') +% title('BER for PAM-4', 'Interpreter','latex') + +grid minor +ylim([1e-4 5e-1]) +set(gca,'YScale','log') +% beautifyBERplot +hold off \ No newline at end of file diff --git a/projects/FSO_transmission/Evaluation Scripts/evalscript_pam_4_baud_rate_sweep.m b/projects/FSO_transmission/Evaluation Scripts/evalscript_pam_4_baud_rate_sweep.m new file mode 100644 index 0000000..0575787 --- /dev/null +++ b/projects/FSO_transmission/Evaluation Scripts/evalscript_pam_4_baud_rate_sweep.m @@ -0,0 +1,77 @@ +baud_rate = 4:1:8; +baud_rate = string(baud_rate); +power = [8.4, 8.4, 42.3, 8.4, 8.4]; +power = string(power); +num_pf_coeff = 4; +taps_ffe = [300, 0, 0]; +taps_dfe = [0, 0, 0]; +M = 4; +trlength = 4096*4; +x = 4:1:8; +BER_PAM_4 = []; +Alpha_PAM_4 = []; + +for method = 1:4 + for i = 1:5 + [BER_run, ~] = first_analysis_baud_rate_sweep(baud_rate(i), power(i), num_pf_coeff, taps_ffe, taps_dfe, M, trlength, method); + BER_PAM_4 = [BER_PAM_4, BER_run]; + % Alpha_PAM_4 = [Alpha_PAM_4, Alpha_run]; + end + save('C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\Silas DSP\imdd_simulation\projects\FSO_transmission\New Baud Rate Sweep Data\BER_PAM_4_' + string(method) + '.mat', 'x', 'BER_PAM_4') + BER_PAM_4 = []; +end + +%% BER +for k = 1:4 + BER = load('C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\Silas DSP\imdd_simulation\projects\FSO_transmission\New Baud Rate Sweep Data\BER_PAM_4_' + string(k) + '.mat'); + BER = BER.BER_PAM_4; + % x = BER.x; + % BER_Alpha = load(filename); + % BER = BER_Alpha.BER_PAM_4; + + % close all + figure(239) + plot(x, BER, '-o','LineWidth',1.75) + hold on +end + +% old_BER = [-4.4, -4, -3.3, -2.7, -2.5, -2.475]; +% old_BER = 10.^(old_BER); +% plot(x, old_BER, '--o','LineWidth',1) + +h1 = yline(2e-2, ':k', 'LineWidth',1.5); +h2 = yline(3.8e-3,':b', 'LineWidth',1.5); +h3 = yline(4.85e-3,':g', 'LineWidth',1.5); +h4 = yline(2.2e-4,':r', 'LineWidth',1.5); + +% Legende NUR für Kurven +legend('FFE','FFE+PF+MLSE','DB','ML-MLSE',... + 'Interpreter','latex', ... + 'Location','southwest', 'FontSize', 14) + +% FEC Labels direkt im Plot +text(8.2,2.3e-2,'o-FEC','Color','k','FontSize', 14, 'Interpreter','latex') +text(8.2,3e-3,'HD-FEC','Color','b','FontSize', 14, 'Interpreter','latex') +text(8.2,5.6e-3,'KP4+Hamming','Color','g','FontSize', 14,'Interpreter','latex') +text(8.2,2.5e-4,'KP4','Color','r','FontSize', 14,'Interpreter','latex') + +xlabel('Symbol Rate [GBd]', 'Interpreter','latex') +ylabel('BER', 'Interpreter','latex') +% title('BER for PAM-4', 'Interpreter','latex') + +grid minor +% ylim([1e-4 5e-2]) +xlim([3 9]) +set(gca,'YScale','log') +hold off + +%% Channel Alpha +% Alpha = BER_Alpha.Alpha_PAM_4; +% figure +% plot(x, Alpha, '--o','LineWidth',1) +% hold on +% xlabel('Symbol Rate [GBd]', 'Interpreter','latex') +% ylabel('Channel Alpha', 'Interpreter','latex') +% title('Channel Alpha for PAM-4 - 300-Tap-FFE - 1 Postfilter Coefficients - 8192 Training Symbols', 'Interpreter','latex') +% grid minor +% hold off \ No newline at end of file diff --git a/projects/FSO_transmission/Evaluation Scripts/evalscript_timing_shift.m b/projects/FSO_transmission/Evaluation Scripts/evalscript_timing_shift.m new file mode 100644 index 0000000..435b540 --- /dev/null +++ b/projects/FSO_transmission/Evaluation Scripts/evalscript_timing_shift.m @@ -0,0 +1,9 @@ +time_shift = -2:0.1:2; +BER = zeros(1,length(time_shift)); +for i = 1:length(time_shift) + BER_value = first_analysis_time_shift(time_shift(i)); + BER(i) = BER_value; + close all +end +figure; +plot(time_shift,BER) \ No newline at end of file diff --git a/projects/FSO_transmission/Evaluation Scripts/first_analysis_2.m b/projects/FSO_transmission/Evaluation Scripts/first_analysis_2.m new file mode 100644 index 0000000..c3a44f0 --- /dev/null +++ b/projects/FSO_transmission/Evaluation Scripts/first_analysis_2.m @@ -0,0 +1,190 @@ +function BER_value = first_analysis_2(taps_ffe, taps_dfe, trlen, trloops, ddloops, K_FFE, DCmu, DDmu, DFEmu, FFEmu, plot_final, idealdfe, num_pf_coeffs, damp_factor, norm_loop_bw, det_gain) + %% + base = "C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\FSO Equalizer\FSO_FP_QCL_60umUTC\"; + mode = 0; %0 oder 1 + M = 2; + + all_files = dir(fullfile(base, "**/*.mat")); + + if M == 2 + tx_data_path = fullfile(base, "14G_PAM2\tx_info\tx_info_PAM2_14Gbd0.75RRC.mat"); + filename = fullfile(base, "14G_PAM2\M=2_Rs=1.4e10_Fs=8e10_I=265mA_RoP=46.3mW_L=31m_PS=RRC_rolloff=0.75_Mode=Rise.mat"); + elseif M == 4 + tx_data_path = fullfile(base, "6G_PAM4\tx_info\tx_info_PAM4_6Gbd0.6RRC.mat"); + filename = fullfile(base, "6G_PAM4\M=4_Rs=6e9_Fs=8e10_I=255mA_RoP=42.3mW_L=31m_PS=RRC_rolloff=0.6_Mode=Rise.mat"); + end + + if mode == 1 + [f, p] = uigetfile(fullfile(base, "**/*.mat")); + if f~=0 + filename = fullfile(p,f); + end + end + + tx_data = load(tx_data_path); + datas = load(filename); + + %% + str = filename; + M_ = str2double(regexp(str, 'M=([^_]+)', 'tokens', 'once')); + assert(M==M_); + fsym = str2double(regexp(str, 'Rs=([^_]+)', 'tokens', 'once')); + fs = str2double(regexp(str, 'Fs=([^_]+)', 'tokens', 'once')); + I = sscanf(char(regexp(str, 'I=([^_]+)', 'tokens', 'once')), '%f'); + rop = sscanf(char(regexp(str, 'RoP=([^_]+)', 'tokens', 'once')), '%f'); + L = sscanf(char(regexp(str, 'L=([^_]+)', 'tokens', 'once')), '%f'); + pulseshape = string( regexp(str, 'PS=([^_]+)', 'tokens', 'once')); + rolloff = str2double(regexp(str, 'rolloff=([^_]+)', 'tokens', 'once')); + mode = string( regexp(str, 'Mode=([^\.]+)', 'tokens', 'once')); + + %% + % Tx data + + Bits = Informationsignal(tx_data.tx_data,"fs",fsym); + Symbols = Informationsignal(real(tx_data.tx_PAM_sym),"fs",fsym); + + mapping_style = M==4; % Pam2 is like move-it; PAM-4 is different, same mapping like ETH peopled used in Zurich... hence the "eth_style" argument here and there + PM = PAMmapper(M,0,"eth_style",mapping_style); % one should rename "eth style" as this is simply a different mapping scheme + + Symbols_ = PM.map(Bits) .* PM.scaling; + assert(isequal(Symbols.signal,Symbols_.signal)); + + Bits_ = PM.demap(Symbols); + [bits,errors,ber,errorIndice] = calc_ber(Bits_.signal,Bits.signal); + assert(ber == 0); + + %% For comparison, apply pulsef on Tx Symbols + Pform = Pulseformer("fsym",fsym,"fdac",fs,"pulse","rc","pulselength",16,"alpha",rolloff); + Digi_sig_compare = Pform.process(Symbols); + MF = Pulseformer("fsym",fsym,"fdac",2*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff); + Rx_sig_compare = MF.process(Digi_sig_compare); + + %% + + % Rx Data + traceData = datas.tr.lastData(2).trace.ch3; + + %FYI: Voltage=(RawData−YReference)×YIncrement+YOrigin + scoperead_volts = (traceData.RawData - traceData.YReference) * traceData.YIncrement + traceData.YOrigin; + demystified = isequal(traceData.YData,scoperead_volts); + assert(demystified); + + Scope_sig = Electricalsignal(traceData.YData,"fs",fs); + + Scope_sig.plot("displayname",'raw','fignum',100); + Scope_sig.spectrum("displayname",'raw','fignum',101) + + % 1) matched filter + % pulse is symmetric, hence we can use pulsef firectly as matched filter. + % It feels off (bit I think correct) that the fsym is now the output freq.!! + % -> output 2 sps to omit timing recovery!? + Matched_Filter = Pulseformer("fsym",fsym,"fdac",2*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff,"matched",1); + Rx_matched = Matched_Filter.process(Scope_sig); + Rx_matched.spectrum("displayname",'Signal after matched filter','fignum',1); + + % timing sync -> at this point we still have no symbol timing recovery, we + % try to do this with 2sps EQ! + [~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_matched.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1); + Rx_matched_1 = Rx_synced_cell{1}; + + % Rx_Time_Rec = Rx_matched; + [Rx_Time_Rec, Timing_Error] = Timing_Recovery("modulation", 'PAM/PSK/QAM', "timing_error_detector",'Gardner (non-data-aided)','sps',2, ... + 'damping_factor',damp_factor,'normalized_loop_bandwidth',norm_loop_bw,'detector_gain',det_gain).process(Rx_matched_1); + figure;plot(Timing_Error); + % Rx_Time_Rec.signal = resample(Rx_Time_Rec.signal, 28e9, 14e9); + + %% not working.. + Rx_synced = Rx_Time_Rec; + Rx_synced.fs = 14e9; + % Rx_synced = Rx_synced_cell{1}; + % len_tr = trlen; + % mu_ffe1 = mu_ffe_values(1); + % mu_ffe2 = mu_ffe_values(2); + % mu_ffe3 = mu_ffe_values(3); + % mu_dc = 0.005; + % mu_ffe = [mu_ffe1 mu_ffe3 mu_ffe3]; + % mu_dfe = mu_dfe_values; + duob_mode = db_mode.no_db; + + Rx_synced.plot("displayname",'RX: Matched+Sync+2sps','fignum',2); + + Digi_sig_compare.normalize("mode","rms").spectrum("displayname",'Tx: RC-shaped','fignum',1,'normalizeTo0dB',0); + Rx_sig_compare.normalize("mode","rms").spectrum("displayname",'Tx: RC-shaped + matched filtered ','fignum',1,'normalizeTo0dB',0); + Rx_synced.normalize("mode","rms").spectrum("displayname",'RX: matched filtered + synced','fignum',1,'normalizeTo0dB',0); + + if M == 2 + ber_in_paper = 10^(-2.6); %fig 3a) 4 Gb/s MWIR FSO Transmission using Directly Modulated QCL and an Uncooled UTC-PD at Room-Temperature + elseif M == 4 + ber_in_paper = 10^(-2.5); + end + + %% -------------------- FFE -------------------- + % % requires some more digging what is going on :-) + % eq_ffe = EQ("Ne",taps_ffe,"Nb",taps_dfe, ... + % "training_length",trlen,"training_loops",trloops,"dd_loops",ddloops, ... + % "K",K_FFE,"DCmu",DCmu,"DDmu",DDmu,"DFEmu",DFEmu, ... + % "FFEmu",FFEmu,"plotfinal",plot_final,"ideal_dfe",idealdfe); + % + % ffe_results = ffe(eq_ffe,M,Rx_synced,Symbols,Bits, ... + % "precode_mode",duob_mode,'showAnalysis',0,"postFFE",[], ... + % "eth_style_symbol_mapping",mapping_style); + % + % % ffe_results.metrics.print + % fprintf('My EQ: %.1e \n',ffe_results.metrics.BER); + % fprintf('Paper: %.1e \n \n',ber_in_paper); + % BER_value = ffe_results.metrics.BER; + + %% -------------------- VNLE + MLSE -------------------- + + pf_ncoeffs = num_pf_coeffs; + eq_v = EQ("Ne",taps_ffe,"Nb",taps_dfe, ... + "training_length",trlen,"training_loops",trloops,"dd_loops",ddloops, ... + "K",K_FFE,"DCmu",DCmu,"DDmu",DDmu,"DFEmu",DFEmu, ... + "FFEmu",FFEmu,"plotfinal",plot_final,"ideal_dfe",idealdfe); + pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); + mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0,"eth_style",mapping_style).levels); + + [vnle_results, mlse_results] = vnle_postfilter_mlse(eq_v, pf_, mlse_, M, Rx_synced, Symbols, Bits, ... + "precode_mode", duob_mode, 'showAnalysis', 1, "postFFE", [], "eth_style_symbol_mapping", mapping_style); + + mlse_results.metrics.print + fprintf('My EQ: %.1e \n',mlse_results.metrics.BER); + fprintf('Paper: %.1e \n \n',ber_in_paper); + BER_value = mlse_results.metrics.BER; + + %% -------------------- DB target -------------------- + % mlse_db_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",M,'trellis_states',PAMmapper(M,0).levels); + % + % eq_ = EQ("Ne",taps_ffe,"Nb",taps_dfe, ... + % "training_length",trlen,"training_loops",trloops,"dd_loops",ddloops, ... + % "K",K_FFE,"DCmu",DCmu,"DDmu",DDmu,"DFEmu",DFEmu, ... + % "FFEmu",FFEmu,"plotfinal",plot_final,"ideal_dfe",idealdfe); + % + % dbt_results = duobinary_target(eq_,mlse_db_, M, Rx_synced, Symbols, Bits, ... + % "precode_mode", duob_mode, 'showAnalysis', 0, "postFFE", [],"eth_style_symbol_mapping",mapping_style); + % + % dbt_results.metrics.print("description",'Duobinary'); + % mlse_results.metrics.print + % fprintf('My EQ: %.1e \n',dbt_results.metrics.BER); + % fprintf('Paper: %.1e \n \n',ber_in_paper); + + %% -------------------- ML-based MLSE (L=2) -------------------- + % ml_mlse_equalizer = ML_MLSE("epochs_tr",100,"epochs_dd",1, ... + % "len_tr",length(Rx_synced)/2,"mu_dd",0.03,"mu_tr",0.03,"order",30,"sps",1, ... + % "traceback_depth",256,"L",2,"delta",10,"adaptive_mu",0); + % + % [ml_mlse_results] = ml_mlse(ml_mlse_equalizer, M, Rx_synced, Symbols, Bits,"precode_mode",duob_mode,"eth_style_symbol_mapping",mapping_style); + % fprintf('ML-based MLSE:'); + % fprintf('My EQ: %.1e \n',ml_mlse_results.metrics.BER); + % fprintf('Paper: %.1e \n \n',ber_in_paper); + % + % % -------------------- Post-FFE -------------------- + % eq_post = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",len_tr,"mu_dd",1e-4,"mu_tr",0,"order",1001,"sps",1,"decide",0); + % post_ffe_results = ffe(eq_ffe,M,Rx_synced,Symbols,Bits, ... + % "precode_mode",duob_mode,'showAnalysis',0,"postFFE",[], ... + % "eth_style_symbol_mapping",mapping_style); + % fprintf('Post-FFE:'); + % fprintf('My EQ: %.1e \n',post_ffe_results.metrics.BER); + % fprintf('Paper: %.1e \n \n',ber_in_paper); + +end \ No newline at end of file diff --git a/projects/FSO_transmission/Evaluation Scripts/first_analysis_baud_rate_sweep.m b/projects/FSO_transmission/Evaluation Scripts/first_analysis_baud_rate_sweep.m new file mode 100644 index 0000000..0fe859c --- /dev/null +++ b/projects/FSO_transmission/Evaluation Scripts/first_analysis_baud_rate_sweep.m @@ -0,0 +1,202 @@ +function [BER, Channel_Alpha] = first_analysis_baud_rate_sweep(baud_rate, power, num_pf_coeff, taps_ffe, taps_dfe, M, trlength, method) + %% + close all + base = "C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\FSO Equalizer\Sweep Data\"; + mode = 0; %0 oder 1 + % M = 2; + + all_files = dir(fullfile(base, "**/*.mat")); + + if M == 2 + tx_data_path = fullfile(base, "14G_PAM2\tx_info\tx_info_PAM2_14Gbd0.75RRC.mat"); + filename = fullfile(base, "14G_PAM2\M=2_Rs=1.4e10_Fs=8e10_I=" + current + "mA_RoP=" + power + "mW_L=31m_PS=RRC_rolloff=0.75_Mode=Rise.mat"); + elseif M == 4 + tx_data_path = fullfile(base, baud_rate + "G_PAM4\tx_info\tx_info_PAM4_" + baud_rate + "Gbd0.6RRC.mat"); + filename = fullfile(base, baud_rate + "G_PAM4\M=4_Rs="+ baud_rate + "e9_Fs=8e10_I=255mA_RoP=" + power + "mW_L=31m_PS=RRC_rolloff=0.6_Mode=Rise.mat"); + end + + if mode == 1 + [f, p] = uigetfile(fullfile(base, "**/*.mat")); + if f~=0 + filename = fullfile(p,f); + end + end + + tx_data = load(tx_data_path); + datas = load(filename); + + %% + str = filename; + M_ = str2double(regexp(str, 'M=([^_]+)', 'tokens', 'once')); + assert(M==M_); + fsym = str2double(regexp(str, 'Rs=([^_]+)', 'tokens', 'once')); + fs = str2double(regexp(str, 'Fs=([^_]+)', 'tokens', 'once')); + I = sscanf(char(regexp(str, 'I=([^_]+)', 'tokens', 'once')), '%f'); + rop = sscanf(char(regexp(str, 'RoP=([^_]+)', 'tokens', 'once')), '%f'); + L = sscanf(char(regexp(str, 'L=([^_]+)', 'tokens', 'once')), '%f'); + pulseshape = string( regexp(str, 'PS=([^_]+)', 'tokens', 'once')); + rolloff = str2double(regexp(str, 'rolloff=([^_]+)', 'tokens', 'once')); + mode = string( regexp(str, 'Mode=([^\.]+)', 'tokens', 'once')); + + %% + % Tx data + + Bits = Informationsignal(tx_data.tx_data,"fs",fsym); + Symbols = Informationsignal(real(tx_data.tx_PAM_sym),"fs",fsym); + + mapping_style = M==4; % Pam2 is like move-it; PAM-4 is different, same mapping like ETH peopled used in Zurich... hence the "eth_style" argument here and there + PM = PAMmapper(M,0,"eth_style",mapping_style); % one should rename "eth style" as this is simply a different mapping scheme + + Symbols_ = PM.map(Bits) .* PM.scaling; + assert(isequal(Symbols.signal,Symbols_.signal)); + + Bits_ = PM.demap(Symbols); + [bits,errors,ber,errorIndice] = calc_ber(Bits_.signal,Bits.signal); + assert(ber == 0); + + %% For comparison, apply pulsef on Tx Symbols + Pform = Pulseformer("fsym",fsym,"fdac",fs,"pulse","rrc","pulselength",16,"alpha",rolloff); + Digi_sig_compare = Pform.process(Symbols); + MF = Pulseformer("fsym",fsym,"fdac",2*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff); + Rx_sig_compare = MF.process(Digi_sig_compare); + + %% + + % Rx Data + traceData = datas.tr.lastData(2).trace.ch3; + + %FYI: Voltage=(RawData−YReference)×YIncrement+YOrigin + scoperead_volts = (traceData.RawData - traceData.YReference) * traceData.YIncrement + traceData.YOrigin; + demystified = isequal(traceData.YData,scoperead_volts); + assert(demystified); + + Scope_sig = Electricalsignal(traceData.YData,"fs",fs); + + Scope_sig.plot("displayname",'raw','fignum',100); + Scope_sig.spectrum("displayname",'raw','fignum',101) + + Kov = 14; + + Scope_sig = Scope_sig.resample('fs_in',fs,'fs_out',Kov*fsym); + + % 1) matched filter + % pulse is symmetric, hence we can use pulsef firectly as matched filter. + % It feels off (bit I think correct) that the fsym is now the output freq.!! + % -> output 2 sps to omit timing recovery!? + Matched_Filter = Pulseformer("fsym",fsym,"fdac",Kov*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff,"matched",1); + Rx_matched = Matched_Filter.process(Scope_sig); + Rx_matched.spectrum("displayname",'Signal after matched filter','fignum',1); + + % timing sync -> at this point we still have no symbol timing recovery, we + % try to do this with 2sps EQ! + [~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_matched.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1); + Rx_matched_1 = Rx_synced_cell{1}; + + % Rx_Time_Rec = Rx_matched; + % Rx_Time_Rec = Timing_Recovery_Move_It('f_sim', 28e9, 'gamma', 0.1).process(Rx_matched_1); + + Time_Rec = 1; + if Time_Rec + Rx_Time_Rec = MaxVar_Timing_Recovery('mode',0,'fsym',fsym,'fadc',Kov*fsym,'num_tau',Kov*128,'sps',Kov,'comp_signal',0,'comp_mode',0).process(Rx_matched_1); + sps = 1; + else + sps = 2; + end + + Rx_Time_Rec.fs = fsym; + + % Rx_Time_Rec.signal = resample(Rx_Time_Rec.signal, 12e9, 6e9); + + % Rx_matched_1.plot("fignum",231231) + % Rx_Time_Rec.plot("fignum",231231) + + %% not working.. + Rx_synced = Rx_Time_Rec; + % Rx_synced = Rx_synced_cell{1}; + len_tr = trlength; + mu_ffe1 = 0.0001; + mu_ffe2 = 0.0008; + mu_ffe3 = 0.001; + mu_dc = 0.004; + mu_ffe = [mu_ffe1 mu_ffe3 mu_ffe3]; + mu_dfe = 0.0004; + duob_mode = db_mode.no_db; + + Rx_synced.plot("displayname",'RX: Matched+Sync+2sps','fignum',103); + Rx_synced.spectrum("displayname",'RX: Matched+Sync+2sps','fignum',104); + + Digi_sig_compare.normalize("mode","rms").spectrum("displayname",'Tx: RC-shaped','fignum',1,'normalizeTo0dB',1); + Rx_sig_compare.normalize("mode","rms").spectrum("displayname",'Tx: RC-shaped + matched filtered ','fignum',1,'normalizeTo0dB',1); + Rx_synced.normalize("mode","rms").spectrum("displayname",'RX: matched filtered + synced','fignum',1,'normalizeTo0dB',1); + + if M == 2 + ber_in_paper = 10^(-2.6); %fig 3a) 4 Gb/s MWIR FSO Transmission using Directly Modulated QCL and an Uncooled UTC-PD at Room-Temperature + elseif M == 4 + ber_in_paper = 10^(-2.5); + end + + if method == 1 + %% -------------------- FFE -------------------- + % requires some more digging what is going on :-) + eq_ffe = EQ("Ne",[500, 0, 0],"Nb",[0, 0, 0], ... + "training_length",len_tr,"training_loops",5,"dd_loops",5, ... + "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... + "FFEmu",0,"plotfinal",0,"ideal_dfe",1); + + ffe_results = ffe(eq_ffe,M,Rx_synced,Symbols,Bits, ... + "precode_mode",duob_mode,'showAnalysis',0,"postFFE",[], ... + "eth_style_symbol_mapping",mapping_style); + + % ffe_results.metrics.print + fprintf('My EQ: %.1e \n',ffe_results.metrics.BER); + fprintf('Paper: %.1e \n \n',ber_in_paper); + BER = ffe_results.metrics.BER; + + elseif method == 2 + %% -------------------- VNLE + MLSE -------------------- + pf_ncoeffs = num_pf_coeff; + eq_v = EQ("Ne",taps_ffe,"Nb",taps_dfe, ... + "training_length",len_tr,"training_loops",5,"dd_loops",5, ... + "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... + "FFEmu",0,"plotfinal",0,"ideal_dfe",1); + pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); + mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0,"eth_style",mapping_style).levels); + + [vnle_results, mlse_results] = vnle_postfilter_mlse(eq_v, pf_, mlse_, M, Rx_synced, Symbols, Bits, ... + "precode_mode", duob_mode, 'showAnalysis', 1, "postFFE", [], "eth_style_symbol_mapping", mapping_style); + + mlse_results.metrics.print + fprintf('My EQ: %.1e \n',mlse_results.metrics.BER); + fprintf('Paper: %.1e \n \n',ber_in_paper); + BER = mlse_results.metrics.BER; + % Channel_Alpha = mlse_results.metrics.Alpha; + + elseif method == 3 + %% -------------------- ML-based MLSE (L=2) -------------------- + ml_mlse_equalizer = ML_MLSE("epochs_tr",150,"epochs_dd",1, ... + "len_tr",length(Rx_synced),"mu_dd",0.03,"mu_tr",0.03,"order",80,"sps",sps, ... + "traceback_depth",256,"L",1,"delta",4,"adaptive_mu",0); + + [ml_mlse_results] = ml_mlse(ml_mlse_equalizer, M, Rx_synced, Symbols, Bits,"precode_mode",duob_mode,"eth_style_symbol_mapping",mapping_style); + fprintf('ML-based MLSE:\n'); + fprintf('My EQ: %.1e \n',ml_mlse_results.metrics.BER); + fprintf('Paper: %.1e \n \n',ber_in_paper); + BER = ml_mlse_results.metrics.BER; + + elseif method == 4 + %% -------------------- DB target -------------------- + mlse_db_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",M,'trellis_states',PAMmapper(M,0).levels); + + eq_ = EQ("Ne",taps_ffe,"Nb",taps_dfe,"training_length",len_tr,"training_loops",5,"dd_loops",5, ... + "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1); + + dbt_results = duobinary_target(eq_,mlse_db_, M, Rx_synced, Symbols, Bits, ... + "precode_mode", duob_mode, 'showAnalysis', 0, "postFFE", [],"eth_style_symbol_mapping",mapping_style); + + dbt_results.metrics.print("description",'Duobinary'); + fprintf('My EQ: %.1e \n',dbt_results.metrics.BER); + fprintf('Paper: %.1e \n \n',ber_in_paper); + BER = dbt_results.metrics.BER; + end + Channel_Alpha = 0; +end \ No newline at end of file diff --git a/projects/FSO_transmission/Evaluation Scripts/first_analysis_ber.m b/projects/FSO_transmission/Evaluation Scripts/first_analysis_ber.m new file mode 100644 index 0000000..bc02b58 --- /dev/null +++ b/projects/FSO_transmission/Evaluation Scripts/first_analysis_ber.m @@ -0,0 +1,207 @@ +function BER_value = first_analysis_ber(current, power, num_pf_coeff, taps_ffe, taps_dfe, M, trlength, eq_method) + %% + close all + + %% + base = "C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\FSO Equalizer\FSO_FP_QCL_60umUTC\"; + mode = 0; %0 oder 1 + % M = 2; + + all_files = dir(fullfile(base, "**/*.mat")); + + if M == 2 + tx_data_path = fullfile(base, "14G_PAM2\tx_info\tx_info_PAM2_14Gbd0.75RRC.mat"); + filename = fullfile(base, "14G_PAM2\M=2_Rs=1.4e10_Fs=8e10_I=" + current + "mA_RoP=" + power + "mW_L=31m_PS=RRC_rolloff=0.75_Mode=Rise.mat"); + elseif M == 4 + tx_data_path = fullfile(base, "6G_PAM4\tx_info\tx_info_PAM4_6Gbd0.6RRC.mat"); + filename = fullfile(base, "6G_PAM4\M=4_Rs=6e9_Fs=8e10_I=" + current + "mA_RoP=" + power + "mW_L=31m_PS=RRC_rolloff=0.6_Mode=Rise.mat"); + end + + if mode == 1 + [f, p] = uigetfile(fullfile(base, "**/*.mat")); + if f~=0 + filename = fullfile(p,f); + end + end + + tx_data = load(tx_data_path); + datas = load(filename); + + %% + str = filename; + M_ = str2double(regexp(str, 'M=([^_]+)', 'tokens', 'once')); + assert(M==M_); + fsym = str2double(regexp(str, 'Rs=([^_]+)', 'tokens', 'once')); + fs = str2double(regexp(str, 'Fs=([^_]+)', 'tokens', 'once')); + I = sscanf(char(regexp(str, 'I=([^_]+)', 'tokens', 'once')), '%f'); + rop = sscanf(char(regexp(str, 'RoP=([^_]+)', 'tokens', 'once')), '%f'); + L = sscanf(char(regexp(str, 'L=([^_]+)', 'tokens', 'once')), '%f'); + pulseshape = string( regexp(str, 'PS=([^_]+)', 'tokens', 'once')); + rolloff = str2double(regexp(str, 'rolloff=([^_]+)', 'tokens', 'once')); + mode = string( regexp(str, 'Mode=([^\.]+)', 'tokens', 'once')); + + %% + % Tx data + + Bits = Informationsignal(tx_data.tx_data,"fs",fsym); + Symbols = Informationsignal(real(tx_data.tx_PAM_sym),"fs",fsym); + + mapping_style = M==4; % Pam2 is like move-it; PAM-4 is different, same mapping like ETH peopled used in Zurich... hence the "eth_style" argument here and there + PM = PAMmapper(M,0,"eth_style",mapping_style); % one should rename "eth style" as this is simply a different mapping scheme + + Symbols_ = PM.map(Bits) .* PM.scaling; + assert(isequal(Symbols.signal,Symbols_.signal)); + + Bits_ = PM.demap(Symbols); + [bits,errors,ber,errorIndice] = calc_ber(Bits_.signal,Bits.signal); + assert(ber == 0); + + %% For comparison, apply pulsef on Tx Symbols + Pform = Pulseformer("fsym",fsym,"fdac",fs,"pulse","rrc","pulselength",16,"alpha",rolloff); + Digi_sig_compare = Pform.process(Symbols); + MF = Pulseformer("fsym",fsym,"fdac",2*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff); + Rx_sig_compare = MF.process(Digi_sig_compare); + + %% + + % Rx Data + traceData = datas.tr.lastData(2).trace.ch3; + + %FYI: Voltage=(RawData−YReference)×YIncrement+YOrigin + scoperead_volts = (traceData.RawData - traceData.YReference) * traceData.YIncrement + traceData.YOrigin; + demystified = isequal(traceData.YData,scoperead_volts); + assert(demystified); + + Scope_sig = Electricalsignal(traceData.YData,"fs",fs); + + Scope_sig.plot("displayname",'raw','fignum',100); + Scope_sig.spectrum("displayname",'raw','fignum',101) + + Kov = 14; + + Scope_sig = Scope_sig.resample('fs_in', fs, 'fs_out', Kov*fsym); + + % 1) matched filter + % pulse is symmetric, hence we can use pulsef firectly as matched filter. + % It feels off (bit I think correct) that the fsym is now the output freq.!! + % -> output 2 sps to omit timing recovery!? + Matched_Filter = Pulseformer("fsym",fsym,"fdac",Kov*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff,"matched",1); + Rx_matched = Matched_Filter.process(Scope_sig); + Rx_matched.spectrum("displayname",'Signal after matched filter','fignum',1); + + % timing sync -> at this point we still have no symbol timing recovery, we + % try to do this with 2sps EQ! + [~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_matched.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1); + Rx_matched_1 = Rx_synced_cell{1}; + + % Rx_Time_Rec = Rx_matched; + % Rx_Time_Rec = Timing_Recovery_Move_It('f_sim', 28e9, 'gamma', 0.1).process(Rx_matched_1); + + Time_Rec = 1; + if Time_Rec + Rx_Time_Rec = MaxVar_Timing_Recovery('mode',0,'fsym',fsym,'fadc',Kov*fsym,'num_tau',Kov*128,'sps',Kov,'comp_signal',0,'comp_mode',0).process(Rx_matched_1); + sps = 1; + else + sps = Kov; + end + Rx_Time_Rec.fs = fsym; + + % Rx_Time_Rec.signal = resample(Rx_Time_Rec.signal, 12e9, 6e9); + + % Rx_matched_1.plot("fignum",231231) + % Rx_Time_Rec.plot("fignum",231231) + + %% not working.. + Rx_synced = Rx_Time_Rec; + % Rx_synced = Rx_synced_cell{1}; + len_tr = trlength; + mu_ffe1 = 0.0001; + mu_ffe2 = 0.0008; + mu_ffe3 = 0.001; + mu_dc = 0.004; + mu_ffe = [mu_ffe1 mu_ffe3 mu_ffe3]; + mu_dfe = 0.0004; + duob_mode = db_mode.no_db; + + Rx_synced.plot("displayname",'RX: Matched+Sync+2sps','fignum',103); + Rx_synced.spectrum("displayname",'RX: Matched+Sync+2sps','fignum',104); + + Digi_sig_compare.normalize("mode","rms").spectrum("displayname",'Tx: RC-shaped','fignum',1,'normalizeTo0dB',1); + Rx_sig_compare.normalize("mode","rms").spectrum("displayname",'Tx: RC-shaped + matched filtered ','fignum',1,'normalizeTo0dB',1); + Rx_synced.normalize("mode","rms").spectrum("displayname",'RX: matched filtered + synced','fignum',1,'normalizeTo0dB',1); + + if M == 2 + ber_in_paper = 10^(-2.6); %fig 3a) 4 Gb/s MWIR FSO Transmission using Directly Modulated QCL and an Uncooled UTC-PD at Room-Temperature + elseif M == 4 + ber_in_paper = 10^(-2.5); + end + + if eq_method == 1 + + %% -------------------- FFE -------------------- + % requires some more digging what is going on :-) + eq_ffe = EQ("Ne",[500, 0, 0],"Nb",[0,0,0], ... + "training_length",len_tr,"training_loops",5,"dd_loops",5, ... + "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... + "FFEmu",0,"plotfinal",0,"ideal_dfe",1); + + ffe_results = ffe(eq_ffe,M,Rx_synced,Symbols,Bits, ... + "precode_mode",duob_mode,'showAnalysis',0,"postFFE",[], ... + "eth_style_symbol_mapping",mapping_style); + + % ffe_results.metrics.print + fprintf('My EQ: %.1e \n',ffe_results.metrics.BER); + fprintf('Paper: %.1e \n \n',ber_in_paper); + BER_value = ffe_results.metrics.BER; + + elseif eq_method == 2 + + %% -------------------- VNLE + MLSE -------------------- + + pf_ncoeffs = num_pf_coeff; + eq_v = EQ("Ne",taps_ffe,"Nb",taps_dfe, ... + "training_length",len_tr,"training_loops",5,"dd_loops",5, ... + "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... + "FFEmu",0,"plotfinal",0,"ideal_dfe",1); + pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); + mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0,"eth_style",mapping_style).levels); + + [vnle_results, mlse_results] = vnle_postfilter_mlse(eq_v, pf_, mlse_, M, Rx_synced, Symbols, Bits, ... + "precode_mode", duob_mode, 'showAnalysis', 1, "postFFE", [], "eth_style_symbol_mapping", mapping_style); + + mlse_results.metrics.print + fprintf('My EQ: %.1e \n',mlse_results.metrics.BER); + fprintf('Paper: %.1e \n \n',ber_in_paper); + BER_value = mlse_results.metrics.BER; + + elseif eq_method == 3 + + %% -------------------- DB target -------------------- + mlse_db_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",M,'trellis_states',PAMmapper(M,0).levels); + + eq_ = EQ("Ne",taps_ffe,"Nb",taps_dfe,"training_length",len_tr,"training_loops",5,"dd_loops",5, ... + "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1); + + dbt_results = duobinary_target(eq_,mlse_db_, M, Rx_synced, Symbols, Bits, ... + "precode_mode", duob_mode, 'showAnalysis', 0, "postFFE", [],"eth_style_symbol_mapping",mapping_style); + + dbt_results.metrics.print("description",'Duobinary'); + fprintf('My EQ: %.1e \n',dbt_results.metrics.BER); + fprintf('Paper: %.1e \n \n',ber_in_paper); + BER_value = dbt_results.metrics.BER; + + elseif eq_method == 4 + + %% -------------------- ML-based MLSE (L=2) -------------------- + ml_mlse_equalizer = ML_MLSE("epochs_tr",150,"epochs_dd",1, ... + "len_tr",length(Rx_synced),"mu_dd",0.03,"mu_tr",0.03,"order",80,"sps",1, ... + "traceback_depth",256,"L",1,"delta",4,"adaptive_mu",0); + + [ml_mlse_results] = ml_mlse(ml_mlse_equalizer, M, Rx_synced, Symbols, Bits,"precode_mode",duob_mode,"eth_style_symbol_mapping",mapping_style); + fprintf('ML-based MLSE:\n'); + fprintf('My EQ: %.1e \n',ml_mlse_results.metrics.BER); + fprintf('Paper: %.1e \n \n',ber_in_paper); + BER_value = ml_mlse_results.metrics.BER; + + end +end \ No newline at end of file diff --git a/projects/FSO_transmission/Evaluation Scripts/first_analysis_mf.m b/projects/FSO_transmission/Evaluation Scripts/first_analysis_mf.m new file mode 100644 index 0000000..b2a6027 --- /dev/null +++ b/projects/FSO_transmission/Evaluation Scripts/first_analysis_mf.m @@ -0,0 +1,293 @@ + +<<<<<<< HEAD +base = "C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\FSO Equalizer\FSO_FP_QCL_60umUTC\"; +======= +base = "C:\Users\Silas\Nextcloud4\Dokumente\02_Ablage_Office\FSO_FP_QCL_60umUTC"; +>>>>>>> 719e5508e776c18e3dcf72a954c71bdedda179e0 +mode = 0; %0 oder 1 +M = 2; + +all_files = dir(fullfile(base, "**/*.mat")); + +if M == 2 +<<<<<<< HEAD + tx_data_path = fullfile(base, "14G_PAM2\tx_info\tx_info_PAM2_14Gbd0.75RRC.mat"); + filename = fullfile(base, "14G_PAM2\M=2_Rs=1.4e10_Fs=8e10_I=265mA_RoP=46.3mW_L=31m_PS=RRC_rolloff=0.75_Mode=Rise.mat"); +elseif M == 4 + tx_data_path = fullfile(base, "6G_PAM4\tx_info\tx_info_PAM4_6Gbd0.6RRC.mat"); +======= + tx_data = load("C:\Users\Silas\Nextcloud4\Dokumente\02_Ablage_Office\FSO_FP_QCL_60umUTC\14G_PAM2\tx_info\tx_info_PAM2_14Gbd0.75RRC.mat"); + filename = fullfile(base, "14G_PAM2\M=2_Rs=1.4e10_Fs=8e10_I=265mA_RoP=46.3mW_L=31m_PS=RRC_rolloff=0.75_Mode=Rise.mat"); +elseif M == 4 + tx_data = load("C:\Users\Silas\Nextcloud4\Dokumente\02_Ablage_Office\FSO_FP_QCL_60umUTC\6G_PAM4\tx_info\tx_info_PAM4_6Gbd0.6RRC.mat"); +>>>>>>> 719e5508e776c18e3dcf72a954c71bdedda179e0 + filename = fullfile(base, "6G_PAM4\M=4_Rs=6e9_Fs=8e10_I=255mA_RoP=42.3mW_L=31m_PS=RRC_rolloff=0.6_Mode=Rise.mat"); +end + +if mode == 1 + [f, p] = uigetfile(fullfile(base, "**/*.mat")); + if f~=0 + filename = fullfile(p,f); + end +end + +tx_data = load(tx_data_path); +datas = load(filename); + +%% +str = filename; +M_ = str2double(regexp(str, 'M=([^_]+)', 'tokens', 'once')); +assert(M==M_); +fsym = str2double(regexp(str, 'Rs=([^_]+)', 'tokens', 'once')); +fs = str2double(regexp(str, 'Fs=([^_]+)', 'tokens', 'once')); +I = sscanf(char(regexp(str, 'I=([^_]+)', 'tokens', 'once')), '%f'); +rop = sscanf(char(regexp(str, 'RoP=([^_]+)', 'tokens', 'once')), '%f'); +L = sscanf(char(regexp(str, 'L=([^_]+)', 'tokens', 'once')), '%f'); +pulseshape = string( regexp(str, 'PS=([^_]+)', 'tokens', 'once')); +rolloff = str2double(regexp(str, 'rolloff=([^_]+)', 'tokens', 'once')); +mode = string( regexp(str, 'Mode=([^\.]+)', 'tokens', 'once')); + +%% +% Tx data + +Bits = Informationsignal(tx_data.tx_data,"fs",fsym); +Symbols = Informationsignal(real(tx_data.tx_PAM_sym),"fs",fsym); + +mapping_style = M==4; % Pam2 is like move-it; PAM-4 is different, same mapping like ETH peopled used in Zurich... hence the "eth_style" argument here and there +PM = PAMmapper(M,0,"eth_style",mapping_style); % one should rename "eth style" as this is simply a different mapping scheme + +Symbols_ = PM.map(Bits) .* PM.scaling; +assert(isequal(Symbols.signal,Symbols_.signal)); + +Bits_ = PM.demap(Symbols); +[bits,errors,ber,errorIndice] = calc_ber(Bits_.signal,Bits.signal); +assert(ber == 0); + +%% For comparison, apply pulsef on Tx Symbols +Pform = Pulseformer("fsym",fsym,"fdac",fs,"pulse","rrc","pulselength",16,"alpha",rolloff); +Digi_sig_tx_compare = Pform.process(Symbols); + +Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff); +Digi_sync = Pform.process(Symbols); + +MF = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff); +Rx_sig_compare = MF.process(Digi_sig_tx_compare); + +%% + +% Rx Data +traceData = datas.tr.lastData(2).trace.ch3; + +%FYI: Voltage=(RawData−YReference)×YIncrement+YOrigin +scoperead_volts = (traceData.RawData - traceData.YReference) * traceData.YIncrement + traceData.YOrigin; +demystified = isequal(traceData.YData,scoperead_volts); +assert(demystified); + +timesig_compare = [0:1:datas.tr.lastData(1).trace.ch3.Points-1] ./ fs; +timesig = datas.tr.lastData(1).trace.ch3.XData; +assert(isequal(traceData.YData,scoperead_volts)); + +Scope_sig = Electricalsignal(traceData.YData,"fs",fs); + +%% +% 1) matched filter +% pulse is symmetric, hence we can use pulsef directly as matched filter. +% It feels off (bit I think correct) that the fsym is now the output freq.!! +% -> output 2 sps to omit timing recovery!? +apply_matched_filter = 0; +k = 4; +if apply_matched_filter + + Pform = Pulseformer("fsym",fsym,"fdac",k*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff,"matched",1); + Rx_matched = Pform.process(Scope_sig); +else + + Rx_matched = Filter('filtdegree',4,"f_cutoff",fsym*0.5,"fs",Scope_sig.fs,"filterType",filtertypes.gaussian,"active",true).process(Scope_sig); + Rx_matched = Rx_matched.resample("fs_out",k*fsym); +end +Rx_matched.spectrum(); + +<<<<<<< HEAD +%% + +sys = comm.SymbolSynchronizer('TimingErrorDetector', 'Gardner (non-data-aided)', ... + 'SamplesPerSymbol', 2, ... + 'DampingFactor', 0.7, ... + 'NormalizedLoopBandwidth', 0.01); +Rx_symbolsync = Rx_matched; +[Rx_symbolsync.signal, timing_error] = sys(Rx_matched.signal); + +plot(timing_error); % If this is a ramp, you have drift! + +%% timing sync -> at this point we still have no symbol timing recovery, we +% % try to do this with 2sps EQ! + +[~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_symbolsync.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1); +% Rx_matched_1 = Rx_synced_cell{1}; +% +% % Rx_Time_Rec = Rx_matched; +% [Rx_Time_Rec, Timing_Error] = Timing_Recovery("modulation", 'PAM/PSK/QAM', "timing_error_detector",'Gardner (non-data-aided)','sps',2,'damping_factor',1,'normalized_loop_bandwidth',0.01,'detector_gain',2.7).process(Rx_matched_1); +% figure;plot(Timing_Error); +% % Rx_Time_Rec.signal = resample(Rx_Time_Rec.signal, 28e9, 14e9); + +%% not working.. +% Rx_synced = Rx_Time_Rec; +Rx_synced = Rx_synced_cell{1}; +len_tr = 4096*2; +mu_ffe1 = 0.0001; +mu_ffe2 = 0.0008; +mu_ffe3 = 0.001; +mu_dc = 0.005; +mu_ffe = [mu_ffe1 mu_ffe3 mu_ffe3]; +mu_dfe = 0.0004; +duob_mode = db_mode.no_db; + +Rx_synced.plot("displayname",'RX: Matched+Sync+2sps','fignum',2); + +Digi_sig_compare.normalize("mode","rms").spectrum("displayname",'Tx: RC-shaped','fignum',1,'normalizeTo0dB',0); +Rx_sig_compare.normalize("mode","rms").spectrum("displayname",'Tx: RC-shaped + matched filtered ','fignum',1,'normalizeTo0dB',0); +Rx_synced.normalize("mode","rms").spectrum("displayname",'RX: matched filtered + synced','fignum',1,'normalizeTo0dB',0); +======= +%% Timing Rec +apply_timing_rec = 1; +if apply_timing_rec + [Rx_symbolsync, timing_error] = Timing_Recovery("timing_error_detector",'Gardner (non-data-aided)','sps',k,'damping_factor',0.1,'normalized_loop_bandwidth',0.1,'detector_gain',2.7).process(Rx_matched); + figure();plot(timing_error); + Rx_symbolsync.fs = fsym; + % Tsynch + [~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_symbolsync.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 0); + Rx_synced = Rx_synced_cell{1}; + sps = 1; +else + % Tsynch + [~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_matched.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 0); + Rx_synced = Rx_synced_cell{1}; + Rx_synced = Rx_synced.resample("fs_out",2*fsym); + sps = 2; +end +>>>>>>> 719e5508e776c18e3dcf72a954c71bdedda179e0 + +if M == 2 + ber_in_paper = 10^(-2.6); %fig 3a) 4 Gb/s MWIR FSO Transmission using Directly Modulated QCL and an Uncooled UTC-PD at Room-Temperature +elseif M == 4 + ber_in_paper = 10^(-2.5); +end + + +%% + +Rx_synced = Rx_synced_cell{1}; + + +% -------------------- FFE -------------------- +% requires some more digging what is going on :-) +eq_ffe = EQ("Ne",[50, 1, 1],"Nb",[2,0,0], ... + "training_length",512,"training_loops",5,"dd_loops",5, ... + "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... + "FFEmu",0,"plotfinal",0,"ideal_dfe",1); + +vars = logspace(-4,-3,36); + +parfor i = 1:numel(vars) + + len_tr = 4096; + mu_ffe1 = 0.01;% mus(i);%0.0001; + mu_ffe2 = 0.0008; + mu_ffe3 = 0.001; + mu_dc = 0.005; + mu_ffe = [mu_ffe1 mu_ffe2 mu_ffe3]; + mu_dfe = vars(i); + duob_mode = db_mode.no_db; + + % requires some more digging what is going on :-) + eq_ffe_1 = EQ("Ne",[150, 1, 0],"Nb",[50,0,0], ... + "training_length",len_tr,"training_loops",5,"dd_loops",5, ... + "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... + "FFEmu",0,"plotfinal",0,"ideal_dfe",0); + + eq_ffe_2 = FFE("epochs_tr",1,"epochs_dd",vars(i),"len_tr",4096,"mu_dd",vars(i),"mu_tr",vars(i),"order",999,"sps",1,"decide",0, "adaption",adaption_method.nlms,"dd_mode",0); + % eq_ffe_2 = FFE_DFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"ffe_mu_dd",1e-5,"dfe_mu_dd",mus(i),"ffe_mu_tr",0,"dfe_mu_tr",0,"ffe_order",50,"dfe_order",10,"sps",1,"decide",1); + + + ffe_results = ffe(eq_ffe_1,M,Rx_synced,Symbols,Bits, ... + "precode_mode",duob_mode,'showAnalysis',1,"postFFE",[], ... + "eth_style_symbol_mapping",mapping_style); + + ffe_results.metrics.BER + bers(i) = ffe_results.metrics.BER; +end + +figure(); +plot(vars,bers); +yline(ber_in_paper) +beautifyBERplot(); +% +fprintf('Paper: %.1e \n \n',ber_in_paper); +ffe_results.metrics.print("description",'FFE'); +fprintf('FFE: %.1e \n',ffe_results.metrics.BER); + + +%% -------------------- VNLE + MLSE -------------------- +len_tr = 4096; +mu_ffe1 = 0.0001;% mus(i);%0.0001; +mu_ffe2 = 0.0008; +mu_ffe3 = 0.001; +mu_dc = 0.005; +mu_ffe = [mu_ffe1 mu_ffe2 mu_ffe3]; +mu_dfe = 0.0004; +duob_mode = db_mode.no_db; + +pf_ncoeffs = 4; + +vars = 1:7; +bers = zeros(size(vars)); +parfor i = 1:numel(vars) + + eqv = EQ("Ne",[200, 1, 0],"Nb",[2, 0, 0], ... + "training_length",len_tr,"training_loops",5,"dd_loops",5, ... + "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... + "FFEmu",0,"plotfinal",0,"ideal_dfe",1); + pf_ncoeffs = vars(i); + pf = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); + mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0,"eth_style",mapping_style).levels); + [vnle_results, mlse_results] = vnle_postfilter_mlse(eqv, pf, mlse_, M, Rx_synced, Symbols, Bits, ... + "precode_mode", duob_mode, 'showAnalysis', 1, "postFFE", [], "eth_style_symbol_mapping", mapping_style); + fprintf('Paper: %.1e \n \n',ber_in_paper); + vnle_results.metrics.print("description",'VNLE'); + mlse_results.metrics.print("description",'MLSE'); + bers(i) = mlse_results.metrics.BER; +end +%% +figure();hold on +plot(vars,bers_ffe,'DisplayName','FFE [200,0,0] + PF + MLSE'); +plot(vars,bers_vnle,'DisplayName','VNLE [200,1,0] + PF + MLSE'); +plot(vars,bers_vnledfe,'DisplayName','VNLE [200,1,0] + DFE [2] + PF + MLSE'); +plot(vars,bers_vnledfe_ideal,'DisplayName','VNLE [200,1,0] + ideal DFE [2] + PF + MLSE'); +yline(ber_in_paper); +yline([2e-2, 4.85e-3, 3.8e-3, 2,2e-4],'LineWidth',2,'Color',[0.8,0.8,0.8],'LineStyle',':','HandleVisibility','off'); +ylim([1e-5,0.1]); +beautifyBERplot(); + +%% -------------------- DB target -------------------- +mlse_db_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",M,'trellis_states',PAMmapper(M,0).levels); + +eq_ = EQ("Ne",[50, 5, 5],"Nb",[0,0,0],"training_length",len_tr,"training_loops",5,"dd_loops",5, ... + "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1); + +dbt_results = duobinary_target(eq_,mlse_db_, M, Rx_synced, Symbols, Bits, ... + "precode_mode", duob_mode, 'showAnalysis', 0, "postFFE", [],"eth_style_symbol_mapping",mapping_style); + +fprintf('Paper: %.1e \n \n',ber_in_paper); +dbt_results.metrics.print("description",'Duobinary'); + + +%% + +%ML-based MLSE (L=2) +mu_ml = 0.01; training_epochs = 100; +ml_mlse_equalizer = ML_MLSE("epochs_tr",training_epochs,"epochs_dd",1, ... + "len_tr",len_tr,"mu_dd",mu_ml,"mu_tr",mu_ml,"order",11,"sps",1, ... + "traceback_depth",128,"L",3,"delta",4,"adaptive_mu",0); + +[ml_mlse_results] = ml_mlse(ml_mlse_equalizer, M, Rx_synced, Symbols, Bits,"precode_mode",duob_mode); +ml_mlse_results.metrics.print("description",'ML '); \ No newline at end of file diff --git a/projects/FSO_transmission/Evaluation Scripts/first_analysis_optimize.m b/projects/FSO_transmission/Evaluation Scripts/first_analysis_optimize.m new file mode 100644 index 0000000..5d5b4f9 --- /dev/null +++ b/projects/FSO_transmission/Evaluation Scripts/first_analysis_optimize.m @@ -0,0 +1,177 @@ +function BER = first_analysis_optimize(baud_rate, power, num_pf_coeff, taps_ffe, taps_dfe, M, trlength) + %% + base = "C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\FSO Equalizer\Sweep Data\"; + mode = 0; %0 oder 1 + % M = 2; + + all_files = dir(fullfile(base, "**/*.mat")); + + if M == 2 + tx_data_path = fullfile(base, "14G_PAM2\tx_info\tx_info_PAM2_14Gbd0.75RRC.mat"); + filename = fullfile(base, "14G_PAM2\M=2_Rs=1.4e10_Fs=8e10_I=" + current + "mA_RoP=" + power + "mW_L=31m_PS=RRC_rolloff=0.75_Mode=Rise.mat"); + elseif M == 4 + tx_data_path = fullfile(base, baud_rate + "G_PAM4\tx_info\tx_info_PAM4_" + baud_rate + "Gbd0.6RRC.mat"); + filename = fullfile(base, baud_rate + "G_PAM4\M=4_Rs="+ baud_rate + "e9_Fs=8e10_I=255mA_RoP=" + power + "mW_L=31m_PS=RRC_rolloff=0.6_Mode=Rise.mat"); + end + + if mode == 1 + [f, p] = uigetfile(fullfile(base, "**/*.mat")); + if f~=0 + filename = fullfile(p,f); + end + end + + tx_data = load(tx_data_path); + datas = load(filename); + + %% + str = filename; + M_ = str2double(regexp(str, 'M=([^_]+)', 'tokens', 'once')); + assert(M==M_); + fsym = str2double(regexp(str, 'Rs=([^_]+)', 'tokens', 'once')); + fs = str2double(regexp(str, 'Fs=([^_]+)', 'tokens', 'once')); + I = sscanf(char(regexp(str, 'I=([^_]+)', 'tokens', 'once')), '%f'); + rop = sscanf(char(regexp(str, 'RoP=([^_]+)', 'tokens', 'once')), '%f'); + L = sscanf(char(regexp(str, 'L=([^_]+)', 'tokens', 'once')), '%f'); + pulseshape = string( regexp(str, 'PS=([^_]+)', 'tokens', 'once')); + rolloff = str2double(regexp(str, 'rolloff=([^_]+)', 'tokens', 'once')); + mode = string( regexp(str, 'Mode=([^\.]+)', 'tokens', 'once')); + + %% + % Tx data + + Bits = Informationsignal(tx_data.tx_data,"fs",fsym); + Symbols = Informationsignal(real(tx_data.tx_PAM_sym),"fs",fsym); + + mapping_style = M==4; % Pam2 is like move-it; PAM-4 is different, same mapping like ETH peopled used in Zurich... hence the "eth_style" argument here and there + PM = PAMmapper(M,0,"eth_style",mapping_style); % one should rename "eth style" as this is simply a different mapping scheme + + Symbols_ = PM.map(Bits) .* PM.scaling; + assert(isequal(Symbols.signal,Symbols_.signal)); + + Bits_ = PM.demap(Symbols); + [bits,errors,ber,errorIndice] = calc_ber(Bits_.signal,Bits.signal); + assert(ber == 0); + + %% For comparison, apply pulsef on Tx Symbols + Pform = Pulseformer("fsym",fsym,"fdac",fs,"pulse","rrc","pulselength",16,"alpha",rolloff); + Digi_sig_compare = Pform.process(Symbols); + MF = Pulseformer("fsym",fsym,"fdac",2*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff); + Rx_sig_compare = MF.process(Digi_sig_compare); + + %% + + % Rx Data + traceData = datas.tr.lastData(2).trace.ch3; + + %FYI: Voltage=(RawData−YReference)×YIncrement+YOrigin + scoperead_volts = (traceData.RawData - traceData.YReference) * traceData.YIncrement + traceData.YOrigin; + demystified = isequal(traceData.YData,scoperead_volts); + assert(demystified); + + Scope_sig = Electricalsignal(traceData.YData,"fs",fs); + + Scope_sig.plot("displayname",'raw','fignum',100); + Scope_sig.spectrum("displayname",'raw','fignum',101) + + % 1) matched filter + % pulse is symmetric, hence we can use pulsef firectly as matched filter. + % It feels off (bit I think correct) that the fsym is now the output freq.!! + % -> output 2 sps to omit timing recovery!? + Matched_Filter = Pulseformer("fsym",fsym,"fdac",2*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff,"matched",1); + Rx_matched = Matched_Filter.process(Scope_sig); + Rx_matched.spectrum("displayname",'Signal after matched filter','fignum',1); + + % timing sync -> at this point we still have no symbol timing recovery, we + % try to do this with 2sps EQ! + [~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_matched.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1); + Rx_matched_1 = Rx_synced_cell{1}; + + % Rx_Time_Rec = Rx_matched; + % Rx_Time_Rec = Timing_Recovery_Move_It('f_sim', 28e9, 'gamma', 0.1).process(Rx_matched_1); + + Time_Rec = 1; + if Time_Rec + [Rx_Time_Rec, Timing_Error] = Timing_Recovery("modulation", 'PAM/PSK/QAM', "timing_error_detector",'Gardner (non-data-aided)','sps',2,'damping_factor',0.0123793,'normalized_loop_bandwidth',4.50865e-06,'detector_gain',3.33311).process(Rx_matched_1); + sps = 1; + else + sps = 2; + end + + baud_rate_num = str2double(baud_rate); + Rx_Time_Rec.fs = baud_rate_num*10^9; + + % Rx_Time_Rec.signal = resample(Rx_Time_Rec.signal, 12e9, 6e9); + + % Rx_matched_1.plot("fignum",231231) + % Rx_Time_Rec.plot("fignum",231231) + + %% not working.. + Rx_synced = Rx_Time_Rec; + % Rx_synced = Rx_synced_cell{1}; + len_tr = trlength; + mu_ffe1 = 0.0001; + mu_ffe2 = 0.0008; + mu_ffe3 = 0.001; + mu_dc = 0.004; + mu_ffe = [mu_ffe1 mu_ffe3 mu_ffe3]; + mu_dfe = 0.0004; + duob_mode = db_mode.no_db; + + Rx_synced.plot("displayname",'RX: Matched+Sync+2sps','fignum',103); + Rx_synced.spectrum("displayname",'RX: Matched+Sync+2sps','fignum',104); + + Digi_sig_compare.normalize("mode","rms").spectrum("displayname",'Tx: RC-shaped','fignum',1,'normalizeTo0dB',1); + Rx_sig_compare.normalize("mode","rms").spectrum("displayname",'Tx: RC-shaped + matched filtered ','fignum',1,'normalizeTo0dB',1); + Rx_synced.normalize("mode","rms").spectrum("displayname",'RX: matched filtered + synced','fignum',1,'normalizeTo0dB',1); + + if M == 2 + ber_in_paper = 10^(-2.6); %fig 3a) 4 Gb/s MWIR FSO Transmission using Directly Modulated QCL and an Uncooled UTC-PD at Room-Temperature + elseif M == 4 + ber_in_paper = 10^(-2.5); + end + + %% -------------------- FFE -------------------- + % % requires some more digging what is going on :-) + % eq_ffe = EQ("Ne",[400, 0, 0],"Nb",[0,0,0], ... + % "training_length",len_tr,"training_loops",5,"dd_loops",5, ... + % "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... + % "FFEmu",0,"plotfinal",0,"ideal_dfe",1); + % + % ffe_results = ffe(eq_ffe,M,Rx_synced,Symbols,Bits, ... + % "precode_mode",duob_mode,'showAnalysis',0,"postFFE",[], ... + % "eth_style_symbol_mapping",mapping_style); + % + % % ffe_results.metrics.print + % fprintf('My EQ: %.1e \n',ffe_results.metrics.BER); + % fprintf('Paper: %.1e \n \n',ber_in_paper); + % BER_value = ffe_results.metrics.BER; + + %% -------------------- VNLE + MLSE -------------------- + + pf_ncoeffs = num_pf_coeff; + eq_v = EQ("Ne",taps_ffe,"Nb",taps_dfe, ... + "training_length",len_tr,"training_loops",5,"dd_loops",5, ... + "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... + "FFEmu",0,"plotfinal",0,"ideal_dfe",1); + pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); + mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0,"eth_style",mapping_style).levels); + + [vnle_results, mlse_results] = vnle_postfilter_mlse(eq_v, pf_, mlse_, M, Rx_synced, Symbols, Bits, ... + "precode_mode", duob_mode, 'showAnalysis', 1, "postFFE", [], "eth_style_symbol_mapping", mapping_style); + + mlse_results.metrics.print + fprintf('My EQ: %.1e \n',mlse_results.metrics.BER); + fprintf('Paper: %.1e \n \n',ber_in_paper); + BER = mlse_results.metrics.BER; + + %% -------------------- ML-based MLSE (L=2) -------------------- + % ml_mlse_equalizer = ML_MLSE("epochs_tr",150,"epochs_dd",1, ... + % "len_tr",length(Rx_synced),"mu_dd",0.03,"mu_tr",0.03,"order",80,"sps",1, ... + % "traceback_depth",256,"L",1,"delta",4,"adaptive_mu",0); + % + % [ml_mlse_results] = ml_mlse(ml_mlse_equalizer, M, Rx_synced, Symbols, Bits,"precode_mode",duob_mode,"eth_style_symbol_mapping",mapping_style); + % fprintf('ML-based MLSE:\n'); + % fprintf('My EQ: %.1e \n',ml_mlse_results.metrics.BER); + % fprintf('Paper: %.1e \n \n',ber_in_paper); +end \ No newline at end of file diff --git a/projects/FSO_transmission/Evaluation Scripts/first_analysis_time_rec.m b/projects/FSO_transmission/Evaluation Scripts/first_analysis_time_rec.m new file mode 100644 index 0000000..22622ec --- /dev/null +++ b/projects/FSO_transmission/Evaluation Scripts/first_analysis_time_rec.m @@ -0,0 +1,160 @@ +%% +base = "C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\FSO Equalizer\FSO_FP_QCL_60umUTC\"; +mode = 0; %0 oder 1 +M = 4; + +all_files = dir(fullfile(base, "**/*.mat")); + +if M == 2 + tx_data_path = fullfile(base, "14G_PAM2\tx_info\tx_info_PAM2_14Gbd0.75RRC.mat"); + filename = fullfile(base, "14G_PAM2\M=2_Rs=1.4e10_Fs=8e10_I=265mA_RoP=46.3mW_L=31m_PS=RRC_rolloff=0.75_Mode=Rise.mat"); +elseif M == 4 + tx_data_path = fullfile(base, "6G_PAM4\tx_info\tx_info_PAM4_6Gbd0.6RRC.mat"); + filename = fullfile(base, "6G_PAM4\M=4_Rs=6e9_Fs=8e10_I=255mA_RoP=42.3mW_L=31m_PS=RRC_rolloff=0.6_Mode=Rise.mat"); +end + +if mode == 1 + [f, p] = uigetfile(fullfile(base, "**/*.mat")); + if f~=0 + filename = fullfile(p,f); + end +end + +tx_data = load(tx_data_path); +datas = load(filename); + +%% +str = filename; +M_ = str2double(regexp(str, 'M=([^_]+)', 'tokens', 'once')); +assert(M==M_); +fsym = str2double(regexp(str, 'Rs=([^_]+)', 'tokens', 'once')); +fs = str2double(regexp(str, 'Fs=([^_]+)', 'tokens', 'once')); +I = sscanf(char(regexp(str, 'I=([^_]+)', 'tokens', 'once')), '%f'); +rop = sscanf(char(regexp(str, 'RoP=([^_]+)', 'tokens', 'once')), '%f'); +L = sscanf(char(regexp(str, 'L=([^_]+)', 'tokens', 'once')), '%f'); +pulseshape = string( regexp(str, 'PS=([^_]+)', 'tokens', 'once')); +rolloff = str2double(regexp(str, 'rolloff=([^_]+)', 'tokens', 'once')); +mode = string( regexp(str, 'Mode=([^\.]+)', 'tokens', 'once')); + +%% +% Tx data + +Bits = Informationsignal(tx_data.tx_data,"fs",fsym); +Symbols = Informationsignal(real(tx_data.tx_PAM_sym),"fs",fsym); + +mapping_style = M==4; % Pam2 is like move-it; PAM-4 is different, same mapping like ETH peopled used in Zurich... hence the "eth_style" argument here and there +PM = PAMmapper(M,0,"eth_style",mapping_style); % one should rename "eth style" as this is simply a different mapping scheme + +Symbols_ = PM.map(Bits) .* PM.scaling; +assert(isequal(Symbols.signal,Symbols_.signal)); + +Bits_ = PM.demap(Symbols); +[bits,errors,ber,errorIndice] = calc_ber(Bits_.signal,Bits.signal); +assert(ber == 0); + +%% For comparison, apply pulsef on Tx Symbols +Pform = Pulseformer("fsym",fsym,"fdac",fs,"pulse","rrc","pulselength",16,"alpha",rolloff); +Digi_sig_compare = Pform.process(Symbols); +MF = Pulseformer("fsym",fsym,"fdac",2*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff); +Rx_sig_compare = MF.process(Digi_sig_compare); + +%% + +% Rx Data +traceData = datas.tr.lastData(2).trace.ch3; + +%FYI: Voltage=(RawData−YReference)×YIncrement+YOrigin +scoperead_volts = (traceData.RawData - traceData.YReference) * traceData.YIncrement + traceData.YOrigin; +demystified = isequal(traceData.YData,scoperead_volts); +assert(demystified); + +Scope_sig = Electricalsignal(traceData.YData,"fs",fs); + +Scope_sig.plot("displayname",'raw','fignum',100); +Scope_sig.spectrum("displayname",'raw','fignum',101) + +% 1) matched filter +% pulse is symmetric, hence we can use pulsef firectly as matched filter. +% It feels off (bit I think correct) that the fsym is now the output freq.!! +% -> output 2 sps to omit timing recovery!? +Matched_Filter = Pulseformer("fsym",fsym,"fdac",2*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff,"matched",1); +Rx_matched = Matched_Filter.process(Scope_sig); +Rx_matched.spectrum("displayname",'Signal after matched filter','fignum',1); + +% timing sync -> at this point we still have no symbol timing recovery, we +% try to do this with 2sps EQ! +[~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_matched.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1); +Rx_matched_1 = Rx_synced_cell{1}; + +% Rx_Time_Rec = Rx_matched; +% Rx_Time_Rec = Timing_Recovery_Move_It('f_sim', 28e9, 'gamma', 0.1).process(Rx_matched_1); +[Rx_Time_Rec, Timing_Error] = Timing_Recovery("modulation", 'PAM/PSK/QAM', "timing_error_detector",'Gardner (non-data-aided)','sps',2,'damping_factor',1,'normalized_loop_bandwidth',0.01,'detector_gain',2.7).process(Rx_matched_1); + +%% not working.. +Rx_synced = Rx_Time_Rec; +% Rx_synced = Rx_synced_cell{1}; +len_tr = 4096*2; +mu_ffe1 = 0.0001; +mu_ffe2 = 0.0008; +mu_ffe3 = 0.001; +mu_dc = 0.004; +mu_ffe = [mu_ffe1 mu_ffe3 mu_ffe3]; +mu_dfe = 0.0004; +duob_mode = db_mode.no_db; +sps = 1; + +Rx_synced.plot("displayname",'RX: Matched+Sync+2sps','fignum',103); +Rx_synced.spectrum("displayname",'RX: Matched+Sync+2sps','fignum',104); + +Digi_sig_compare.normalize("mode","rms").spectrum("displayname",'Tx: RC-shaped','fignum',1,'normalizeTo0dB',1); +Rx_sig_compare.normalize("mode","rms").spectrum("displayname",'Tx: RC-shaped + matched filtered ','fignum',1,'normalizeTo0dB',1); +Rx_synced.normalize("mode","rms").spectrum("displayname",'RX: matched filtered + synced','fignum',1,'normalizeTo0dB',1); + +if M == 2 + ber_in_paper = 10^(-2.6); %fig 3a) 4 Gb/s MWIR FSO Transmission using Directly Modulated QCL and an Uncooled UTC-PD at Room-Temperature +elseif M == 4 + ber_in_paper = 10^(-2.5); +end + +%% -------------------- FFE -------------------- +% % requires some more digging what is going on :-) +% eq_ffe = EQ("Ne",[50, 5, 5],"Nb",[2,0,0], ... +% "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); +% +% ffe_results = ffe(eq_ffe,M,Rx_synced,Symbols,Bits, ... +% "precode_mode",duob_mode,'showAnalysis',0,"postFFE",[], ... +% "eth_style_symbol_mapping",mapping_style); +% +% % ffe_results.metrics.print +% fprintf('My EQ: %.1e \n',ffe_results.metrics.BER); +% fprintf('Paper: %.1e \n \n',ber_in_paper); +% BER_value = ffe_results.metrics.BER; + +%% -------------------- VNLE + MLSE -------------------- + +pf_ncoeffs = 1; +eq_v = EQ("Ne",[200, 3, 2],"Nb",[5, 2, 1], ... + "training_length",len_tr,"training_loops",5,"dd_loops",5, ... + "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... + "FFEmu",0,"plotfinal",0,"ideal_dfe",1); +pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); +mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0,"eth_style",mapping_style).levels); + +[vnle_results, mlse_results] = vnle_postfilter_mlse(eq_v, pf_, mlse_, M, Rx_synced, Symbols, Bits, ... + "precode_mode", duob_mode, 'showAnalysis', 1, "postFFE", [], "eth_style_symbol_mapping", mapping_style); + +mlse_results.metrics.print +fprintf('My EQ: %.1e \n',mlse_results.metrics.BER); +fprintf('Paper: %.1e \n \n',ber_in_paper); + +%% -------------------- ML-based MLSE (L=2) -------------------- +% ml_mlse_equalizer = ML_MLSE("epochs_tr",150,"epochs_dd",1, ... +% "len_tr",length(Rx_synced),"mu_dd",0.03,"mu_tr",0.03,"order",80,"sps",1, ... +% "traceback_depth",256,"L",1,"delta",4,"adaptive_mu",0); +% +% [ml_mlse_results] = ml_mlse(ml_mlse_equalizer, M, Rx_synced, Symbols, Bits,"precode_mode",duob_mode,"eth_style_symbol_mapping",mapping_style); +% fprintf('ML-based MLSE:\n'); +% fprintf('My EQ: %.1e \n',ml_mlse_results.metrics.BER); +% fprintf('Paper: %.1e \n \n',ber_in_paper); \ No newline at end of file diff --git a/projects/FSO_transmission/Evaluation Scripts/first_analysis_time_rec_pam2.m b/projects/FSO_transmission/Evaluation Scripts/first_analysis_time_rec_pam2.m new file mode 100644 index 0000000..e681f2e --- /dev/null +++ b/projects/FSO_transmission/Evaluation Scripts/first_analysis_time_rec_pam2.m @@ -0,0 +1,171 @@ +%% +base = "C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\FSO Equalizer\FSO_FP_QCL_60umUTC\"; +mode = 0; %0 oder 1 +M = 2; + +all_files = dir(fullfile(base, "**/*.mat")); + +if M == 2 + tx_data_path = fullfile(base, "14G_PAM2\tx_info\tx_info_PAM2_14Gbd0.75RRC.mat"); + filename = fullfile(base, "14G_PAM2\M=2_Rs=1.4e10_Fs=8e10_I=265mA_RoP=46.3mW_L=31m_PS=RRC_rolloff=0.75_Mode=Rise.mat"); +elseif M == 4 + tx_data_path = fullfile(base, "6G_PAM4\tx_info\tx_info_PAM4_6Gbd0.6RRC.mat"); + filename = fullfile(base, "6G_PAM4\M=4_Rs=6e9_Fs=8e10_I=255mA_RoP=42.3mW_L=31m_PS=RRC_rolloff=0.6_Mode=Rise.mat"); +end + +if mode == 1 + [f, p] = uigetfile(fullfile(base, "**/*.mat")); + if f~=0 + filename = fullfile(p,f); + end +end + +tx_data = load(tx_data_path); +datas = load(filename); + +%% +str = filename; +M_ = str2double(regexp(str, 'M=([^_]+)', 'tokens', 'once')); +assert(M==M_); +fsym = str2double(regexp(str, 'Rs=([^_]+)', 'tokens', 'once')); +fs = str2double(regexp(str, 'Fs=([^_]+)', 'tokens', 'once')); +I = sscanf(char(regexp(str, 'I=([^_]+)', 'tokens', 'once')), '%f'); +rop = sscanf(char(regexp(str, 'RoP=([^_]+)', 'tokens', 'once')), '%f'); +L = sscanf(char(regexp(str, 'L=([^_]+)', 'tokens', 'once')), '%f'); +pulseshape = string( regexp(str, 'PS=([^_]+)', 'tokens', 'once')); +rolloff = str2double(regexp(str, 'rolloff=([^_]+)', 'tokens', 'once')); +mode = string( regexp(str, 'Mode=([^\.]+)', 'tokens', 'once')); + +%% +% Tx data + +Bits = Informationsignal(tx_data.tx_data,"fs",fsym); +Symbols = Informationsignal(real(tx_data.tx_PAM_sym),"fs",fsym); + +mapping_style = M==4; % Pam2 is like move-it; PAM-4 is different, same mapping like ETH peopled used in Zurich... hence the "eth_style" argument here and there +PM = PAMmapper(M,0,"eth_style",mapping_style); % one should rename "eth style" as this is simply a different mapping scheme + +Symbols_ = PM.map(Bits) .* PM.scaling; +assert(isequal(Symbols.signal,Symbols_.signal)); + +Bits_ = PM.demap(Symbols); +[bits,errors,ber,errorIndice] = calc_ber(Bits_.signal,Bits.signal); +assert(ber == 0); + +%% For comparison, apply pulsef on Tx Symbols +Pform = Pulseformer("fsym",fsym,"fdac",fs,"pulse","rrc","pulselength",16,"alpha",rolloff); +Digi_sig_compare = Pform.process(Symbols); +MF = Pulseformer("fsym",fsym,"fdac",2*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff); +Rx_sig_compare = MF.process(Digi_sig_compare); + +%% + +% Rx Data +traceData = datas.tr.lastData(2).trace.ch3; + +%FYI: Voltage=(RawData−YReference)×YIncrement+YOrigin +scoperead_volts = (traceData.RawData - traceData.YReference) * traceData.YIncrement + traceData.YOrigin; +demystified = isequal(traceData.YData,scoperead_volts); +assert(demystified); + +Scope_sig = Electricalsignal(traceData.YData,"fs",fs); + +Scope_sig.plot("displayname",'raw','fignum',100); +Scope_sig.spectrum("displayname",'raw','fignum',101) + +% 1) matched filter +% pulse is symmetric, hence we can use pulsef firectly as matched filter. +% It feels off (bit I think correct) that the fsym is now the output freq.!! +% -> output 2 sps to omit timing recovery!? +Matched_Filter = Pulseformer("fsym",fsym,"fdac",2*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff,"matched",1); +Rx_matched = Matched_Filter.process(Scope_sig); +Rx_matched.spectrum("displayname",'Signal after matched filter','fignum',1); + +% timing sync -> at this point we still have no symbol timing recovery, we +% try to do this with 2sps EQ! +[~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_matched.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1); +Rx_matched_1 = Rx_synced_cell{1}; + +Time_Rec = 1; +if Time_Rec + [Rx_Time_Rec, Timing_Error] = Timing_Recovery("modulation", 'PAM/PSK/QAM', "timing_error_detector",'Gardner (non-data-aided)','sps',2,'damping_factor',1,'normalized_loop_bandwidth',0.01,'detector_gain',2.7).process(Rx_matched_1); + sps = 1; +else + sps = 2; +end + +if Time_Rec == 1 && M == 2 + Rx_Time_Rec.fs = 14e9; +elseif Time_Rec == 1 && M == 4 + Rx_Time_Rec.fs = 6e9; +end + +%% not working.. +Rx_synced = Rx_Time_Rec; +% Rx_synced = Rx_synced_cell{1}; +len_tr = 4096*4; +mu_ffe1 = 0.0001; +mu_ffe2 = 0.0008; +mu_ffe3 = 0.001; +mu_dc = 0.004; +mu_ffe = [mu_ffe1 mu_ffe3 mu_ffe3]; +mu_dfe = 0.0004; +duob_mode = db_mode.no_db; +sps = 1; + +Rx_synced.plot("displayname",'RX: Matched+Sync+2sps','fignum',103); +Rx_synced.spectrum("displayname",'RX: Matched+Sync+2sps','fignum',104); + +Digi_sig_compare.normalize("mode","rms").spectrum("displayname",'Tx: RC-shaped','fignum',1,'normalizeTo0dB',1); +Rx_sig_compare.normalize("mode","rms").spectrum("displayname",'Tx: RC-shaped + matched filtered ','fignum',1,'normalizeTo0dB',1); +Rx_synced.normalize("mode","rms").spectrum("displayname",'RX: matched filtered + synced','fignum',1,'normalizeTo0dB',1); + +if M == 2 + ber_in_paper = 10^(-2.6); %fig 3a) 4 Gb/s MWIR FSO Transmission using Directly Modulated QCL and an Uncooled UTC-PD at Room-Temperature +elseif M == 4 + ber_in_paper = 10^(-2.5); +end + +%% -------------------- FFE -------------------- +% % requires some more digging what is going on :-) +% eq_ffe = EQ("Ne",[200, 0, 0],"Nb",[0,0,0], ... +% "training_length",len_tr,"training_loops",5,"dd_loops",5, ... +% "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... +% "FFEmu",0,"plotfinal",0,"ideal_dfe",1); +% +% ffe_results = ffe(eq_ffe,M,Rx_synced,Symbols,Bits, ... +% "precode_mode",duob_mode,'showAnalysis',0,"postFFE",[], ... +% "eth_style_symbol_mapping",mapping_style); +% +% % ffe_results.metrics.print +% fprintf('My EQ: %.1e \n',ffe_results.metrics.BER); +% fprintf('Paper: %.1e \n \n',ber_in_paper); +% BER_value = ffe_results.metrics.BER; + +%% -------------------- VNLE + MLSE -------------------- + +pf_ncoeffs = 4; +eq_v = EQ("Ne",[250, 0, 0],"Nb",[7, 0, 0], ... + "training_length",len_tr,"training_loops",5,"dd_loops",5, ... + "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... + "FFEmu",0,"plotfinal",0,"ideal_dfe",0, ... + 'weighted_DFE',1,'weighted_DFE_d_min',0.5,'weighted_DFE_mode','I2','weighted_DFE_I_mode',[5,0.5,0.6]); +pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); +mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0,"eth_style",mapping_style).levels); + +[vnle_results, mlse_results] = vnle_postfilter_mlse(eq_v, pf_, mlse_, M, Rx_synced, Symbols, Bits, ... + "precode_mode", duob_mode, 'showAnalysis', 1, "postFFE", [], "eth_style_symbol_mapping", mapping_style); + +mlse_results.metrics.print +fprintf('My EQ: %.1e \n',mlse_results.metrics.BER); +fprintf('Paper: %.1e \n \n',ber_in_paper); + +%% -------------------- ML-based MLSE (L=2) -------------------- +% ml_mlse_equalizer = ML_MLSE("epochs_tr",400,"epochs_dd",1, ... +% "len_tr",len_tr,"mu_dd",0.03,"mu_tr",0.03,"order",100,"sps",sps, ... +% "traceback_depth",256,"L",2,"delta",4,"adaptive_mu",0); +% +% [ml_mlse_results] = ml_mlse(ml_mlse_equalizer, M, Rx_synced, Symbols, Bits,"precode_mode",duob_mode,"eth_style_symbol_mapping",mapping_style); +% fprintf('ML-based MLSE:\n'); +% fprintf('My EQ: %.1e \n',ml_mlse_results.metrics.BER); +% fprintf('Paper: %.1e \n \n',ber_in_paper); \ No newline at end of file diff --git a/projects/FSO_transmission/Evaluation Scripts/first_analysis_time_rec_pam4.m b/projects/FSO_transmission/Evaluation Scripts/first_analysis_time_rec_pam4.m new file mode 100644 index 0000000..2881ccf --- /dev/null +++ b/projects/FSO_transmission/Evaluation Scripts/first_analysis_time_rec_pam4.m @@ -0,0 +1,189 @@ +%% +base = "C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\FSO Equalizer\Sweep Data\"; +mode = 0; %0 oder 1 +M = 4; + +baud_rate = '8'; +baud_rate_num = str2double(baud_rate); + +all_files = dir(fullfile(base, "**/*.mat")); + +if M == 2 + tx_data_path = fullfile(base, baud_rate + "G_PAM2\tx_info\tx_info_PAM2_" + baud_rate + "Gbd0.75RRC.mat"); + filename = fullfile(base, baud_rate + "G_PAM2\M=2_Rs=1.4e10_Fs=8e10_I=265mA_RoP=46.3mW_L=31m_PS=RRC_rolloff=0.75_Mode=Rise.mat"); +elseif M == 4 + tx_data_path = fullfile(base, baud_rate + "G_PAM4\tx_info\tx_info_PAM4_" + baud_rate + "Gbd0.6RRC.mat"); + filename = fullfile(base, baud_rate + "G_PAM4\M=4_Rs=" + baud_rate + "e9_Fs=8e10_I=255mA_RoP=8.4mW_L=31m_PS=RRC_rolloff=0.6_Mode=Rise.mat"); +end + +if mode == 1 + [f, p] = uigetfile(fullfile(base, "**/*.mat")); + if f~=0 + filename = fullfile(p,f); + end +end + +tx_data = load(tx_data_path); +datas = load(filename); + +%% +str = filename; +M_ = str2double(regexp(str, 'M=([^_]+)', 'tokens', 'once')); +assert(M==M_); +fsym = str2double(regexp(str, 'Rs=([^_]+)', 'tokens', 'once')); +fs = str2double(regexp(str, 'Fs=([^_]+)', 'tokens', 'once')); +I = sscanf(char(regexp(str, 'I=([^_]+)', 'tokens', 'once')), '%f'); +rop = sscanf(char(regexp(str, 'RoP=([^_]+)', 'tokens', 'once')), '%f'); +L = sscanf(char(regexp(str, 'L=([^_]+)', 'tokens', 'once')), '%f'); +pulseshape = string( regexp(str, 'PS=([^_]+)', 'tokens', 'once')); +rolloff = str2double(regexp(str, 'rolloff=([^_]+)', 'tokens', 'once')); +mode = string( regexp(str, 'Mode=([^\.]+)', 'tokens', 'once')); + +%% +% Tx data + +Bits = Informationsignal(tx_data.tx_data,"fs",fsym); +Symbols = Informationsignal(real(tx_data.tx_PAM_sym),"fs",fsym); + +mapping_style = M==4; % Pam2 is like move-it; PAM-4 is different, same mapping like ETH peopled used in Zurich... hence the "eth_style" argument here and there +PM = PAMmapper(M,0,"eth_style",mapping_style); % one should rename "eth style" as this is simply a different mapping scheme + +Symbols_ = PM.map(Bits) .* PM.scaling; +assert(isequal(Symbols.signal,Symbols_.signal)); + +Bits_ = PM.demap(Symbols); +[bits,errors,ber,errorIndice] = calc_ber(Bits_.signal,Bits.signal); +assert(ber == 0); + +%% For comparison, apply pulsef on Tx Symbols +Pform = Pulseformer("fsym",fsym,"fdac",fs,"pulse","rrc","pulselength",16,"alpha",rolloff); +Digi_sig_compare = Pform.process(Symbols); +MF = Pulseformer("fsym",fsym,"fdac",2*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff); +Rx_sig_compare = MF.process(Digi_sig_compare); + +%% + +% Rx Data +traceData = datas.tr.lastData(2).trace.ch3; + +%FYI: Voltage=(RawData−YReference)×YIncrement+YOrigin +scoperead_volts = (traceData.RawData - traceData.YReference) * traceData.YIncrement + traceData.YOrigin; +demystified = isequal(traceData.YData,scoperead_volts); +assert(demystified); + +Scope_sig = Electricalsignal(traceData.YData,"fs",fs); + +Scope_sig.plot("displayname",'raw','fignum',100); +Scope_sig.spectrum("displayname",'raw','fignum',101) + +% 1) matched filter +% pulse is symmetric, hence we can use pulsef firectly as matched filter. +% It feels off (bit I think correct) that the fsym is now the output freq.!! +% -> output 2 sps to omit timing recovery!? +Matched_Filter = Pulseformer("fsym",fsym,"fdac",2*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff,"matched",1); +Rx_matched = Matched_Filter.process(Scope_sig); +Rx_matched.spectrum("displayname",'Signal after matched filter','fignum',1); + +% timing sync -> at this point we still have no symbol timing recovery, we +% try to do this with 2sps EQ! +[~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_matched.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1); +Rx_matched_1 = Rx_synced_cell{1}; + +% Rx_Time_Rec = Rx_matched; +% Rx_Time_Rec = Timing_Recovery_Move_It('f_sim', 28e9, 'gamma', 0.1).process(Rx_matched_1); + +Time_Rec = 1; +if Time_Rec + [Rx_Time_Rec, Timing_Error] = Timing_Recovery("modulation", 'PAM/PSK/QAM', "timing_error_detector",'Gardner (non-data-aided)','sps',2,'damping_factor',1,'normalized_loop_bandwidth',0.01,'detector_gain',2.7).process(Rx_matched_1); + sps = 1; +else + sps = 2; +end + +Rx_Time_Rec.fs = baud_rate_num * 10^9; + +% Rx_Time_Rec.signal = resample(Rx_Time_Rec.signal, 12e9, 6e9); + +% Rx_matched_1.plot("fignum",231231) +% Rx_Time_Rec.plot("fignum",231231) + +%% not working.. +Rx_synced = Rx_Time_Rec; +% Rx_synced = Rx_synced_cell{1}; +len_tr = 4096*2; +mu_ffe1 = 0.0001; +mu_ffe2 = 0.0008; +mu_ffe3 = 0.001; +mu_dc = 0.004; +mu_ffe = [mu_ffe1 mu_ffe3 mu_ffe3]; +mu_dfe = 0.0004; +duob_mode = db_mode.no_db; + +Rx_synced.plot("displayname",'RX: Matched+Sync+2sps','fignum',103); +Rx_synced.spectrum("displayname",'RX: Matched+Sync+2sps','fignum',104); + +Digi_sig_compare.normalize("mode","rms").spectrum("displayname",'Tx: RC-shaped','fignum',1,'normalizeTo0dB',1); +Rx_sig_compare.normalize("mode","rms").spectrum("displayname",'Tx: RC-shaped + matched filtered ','fignum',1,'normalizeTo0dB',1); +Rx_synced.normalize("mode","rms").spectrum("displayname",'RX: matched filtered + synced','fignum',1,'normalizeTo0dB',1); + +if M == 2 + ber_in_paper = 10^(-2.6); %fig 3a) 4 Gb/s MWIR FSO Transmission using Directly Modulated QCL and an Uncooled UTC-PD at Room-Temperature +elseif M == 4 + ber_in_paper = 10^(-2.5); +end + +%% -------------------- FFE -------------------- +% % requires some more digging what is going on :-) +% eq_ffe = EQ("Ne",[400, 0, 0],"Nb",[0,0,0], ... +% "training_length",len_tr,"training_loops",5,"dd_loops",5, ... +% "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... +% "FFEmu",0,"plotfinal",0,"ideal_dfe",1); +% +% ffe_results = ffe(eq_ffe,M,Rx_synced,Symbols,Bits, ... +% "precode_mode",duob_mode,'showAnalysis',0,"postFFE",[], ... +% "eth_style_symbol_mapping",mapping_style); +% +% % ffe_results.metrics.print +% fprintf('My EQ: %.1e \n',ffe_results.metrics.BER); +% fprintf('Paper: %.1e \n \n',ber_in_paper); +% BER_value = ffe_results.metrics.BER; + +%% -------------------- VNLE + MLSE -------------------- + +pf_ncoeffs = 4; +eq_v = EQ("Ne",[300, 0, 0],"Nb",[0, 0, 0], ... + "training_length",len_tr,"training_loops",5,"dd_loops",5, ... + "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... + "FFEmu",0,"plotfinal",0,"ideal_dfe",1); +pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); +mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0,"eth_style",mapping_style).levels); + +[vnle_results, mlse_results] = vnle_postfilter_mlse(eq_v, pf_, mlse_, M, Rx_synced, Symbols, Bits, ... + "precode_mode", duob_mode, 'showAnalysis', 1, "postFFE", [], "eth_style_symbol_mapping", mapping_style); + +mlse_results.metrics.print +fprintf('My EQ: %.1e \n',mlse_results.metrics.BER); +fprintf('Paper: %.1e \n \n',ber_in_paper); + +%% -------------------- ML-based MLSE (L=2) -------------------- +% ml_mlse_equalizer = ML_MLSE("epochs_tr",200,"epochs_dd",1, ... +% "len_tr",length(Rx_synced)/2,"mu_dd",0,"mu_tr",0.03,"order",100,"sps",sps, ... +% "traceback_depth",256,"L",3,"delta",4,"adaptive_mu",0); +% +% [ml_mlse_results] = ml_mlse(ml_mlse_equalizer, M, Rx_synced, Symbols, Bits,"precode_mode",duob_mode,"eth_style_symbol_mapping",mapping_style); +% fprintf('ML-based MLSE:\n'); +% fprintf('My EQ: %.1e \n',ml_mlse_results.metrics.BER); +% fprintf('Paper: %.1e \n \n',ber_in_paper); + +%% -------------------- DB target -------------------- +% mlse_db_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",M,'trellis_states',PAMmapper(M,0).levels); +% +% eq_ = EQ("Ne",[50, 5, 5],"Nb",[0,0,0],"training_length",len_tr,"training_loops",5,"dd_loops",5, ... +% "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1); +% +% dbt_results = duobinary_target(eq_,mlse_db_, M, Rx_synced, Symbols, Bits, ... +% "precode_mode", duob_mode, 'showAnalysis', 0, "postFFE", [],"eth_style_symbol_mapping",mapping_style); +% +% dbt_results.metrics.print("description",'Duobinary'); +% fprintf('My EQ: %.1e \n',dbt_results.metrics.BER); +% fprintf('Paper: %.1e \n \n',ber_in_paper); \ No newline at end of file diff --git a/projects/FSO_transmission/Evaluation Scripts/first_analysis_time_shift.m b/projects/FSO_transmission/Evaluation Scripts/first_analysis_time_shift.m new file mode 100644 index 0000000..162d6a6 --- /dev/null +++ b/projects/FSO_transmission/Evaluation Scripts/first_analysis_time_shift.m @@ -0,0 +1,221 @@ +function BER = first_analysis_time_shift(time_shift) + %% + base = "C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\FSO Equalizer\FSO_FP_QCL_60umUTC\"; + mode = 0; %0 oder 1 + M = 4; + + all_files = dir(fullfile(base, "**/*.mat")); + + % data_tr_mf is the already recovered and filtered data + if M == 2 + tx_data_path = fullfile(base, "14G_PAM2\tx_info\tx_info_PAM2_14Gbd0.75RRC.mat"); + filename = fullfile(base, "14G_PAM2\M=2_Rs=1.4e10_Fs=8e10_I=265mA_RoP=46.3mW_L=31m_PS=RRC_rolloff=0.75_Mode=Rise.mat"); + data_tr_mf = load("C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\FSO Equalizer\FSO_FP_QCL_60umUTC\Already Recovered and Filtered\AfterSync_M=2_Rs=1.4e10_Fs=8e10_I=265mA_RoP=46.3mW_L=31m_PS=RRC_rolloff=0.75_Mode=Rise.mat"); + elseif M == 4 + tx_data_path = fullfile(base, "6G_PAM4\tx_info\tx_info_PAM4_6Gbd0.6RRC.mat"); + filename = fullfile(base, "6G_PAM4\M=4_Rs=6e9_Fs=8e10_I=255mA_RoP=42.3mW_L=31m_PS=RRC_rolloff=0.6_Mode=Rise.mat"); + data_tr_mf = load("C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\FSO Equalizer\FSO_FP_QCL_60umUTC\Already Recovered and Filtered\AfterSync_M=4_Rs=6e9_Fs=8e10_I=255mA_RoP=42.3mW_L=31m_PS=RRC_rolloff=0.6_Mode=Rise.mat"); + end + + if mode == 1 + [f, p] = uigetfile(fullfile(base, "**/*.mat")); + if f~=0 + filename = fullfile(p,f); + end + end + + tx_data = load(tx_data_path); + datas = load(filename); + + %% + str = filename; + M_ = str2double(regexp(str, 'M=([^_]+)', 'tokens', 'once')); + assert(M==M_); + fsym = str2double(regexp(str, 'Rs=([^_]+)', 'tokens', 'once')); + fs = str2double(regexp(str, 'Fs=([^_]+)', 'tokens', 'once')); + I = sscanf(char(regexp(str, 'I=([^_]+)', 'tokens', 'once')), '%f'); + rop = sscanf(char(regexp(str, 'RoP=([^_]+)', 'tokens', 'once')), '%f'); + L = sscanf(char(regexp(str, 'L=([^_]+)', 'tokens', 'once')), '%f'); + pulseshape = string( regexp(str, 'PS=([^_]+)', 'tokens', 'once')); + rolloff = str2double(regexp(str, 'rolloff=([^_]+)', 'tokens', 'once')); + mode = string( regexp(str, 'Mode=([^\.]+)', 'tokens', 'once')); + + %% + % Tx data + + Bits = Informationsignal(tx_data.tx_data,"fs",fsym); + Symbols = Informationsignal(real(tx_data.tx_PAM_sym),"fs",fsym); + + mapping_style = M==4; % Pam2 is like move-it; PAM-4 is different, same mapping like ETH peopled used in Zurich... hence the "eth_style" argument here and there + PM = PAMmapper(M,0,"eth_style",mapping_style); % one should rename "eth style" as this is simply a different mapping scheme + + Symbols_ = PM.map(Bits) .* PM.scaling; + assert(isequal(Symbols.signal,Symbols_.signal)); + + Bits_ = PM.demap(Symbols); + [bits,errors,ber,errorIndice] = calc_ber(Bits_.signal,Bits.signal); + assert(ber == 0); + + %% For comparison, apply pulsef on Tx Symbols + Pform = Pulseformer("fsym",fsym,"fdac",fs,"pulse","rrc","pulselength",16,"alpha",rolloff); + Digi_sig_compare = Pform.process(Symbols); + MF = Pulseformer("fsym",fsym,"fdac",2*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff); + Rx_sig_compare = MF.process(Digi_sig_compare); + + %% + + % Rx Data + traceData = datas.tr.lastData(2).trace.ch3; + + %FYI: Voltage=(RawData−YReference)×YIncrement+YOrigin + scoperead_volts = (traceData.RawData - traceData.YReference) * traceData.YIncrement + traceData.YOrigin; + demystified = isequal(traceData.YData,scoperead_volts); + assert(demystified); + + Scope_sig = Electricalsignal(traceData.YData,"fs",fs); + + Scope_sig.plot("displayname",'raw','fignum',100); + Scope_sig.spectrum("displayname",'raw','fignum',101) + + % 1) matched filter + % pulse is symmetric, hence we can use pulsef firectly as matched filter. + % It feels off (bit I think correct) that the fsym is now the output freq.!! + % -> output 2 sps to omit timing recovery!? + Matched_Filter = Pulseformer("fsym",fsym,"fdac",2*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff,"matched",1); + Rx_matched = Matched_Filter.process(Scope_sig); + Rx_matched.spectrum("displayname",'Signal after matched filter','fignum',1); + + data_tr_mf = Electricalsignal(data_tr_mf.Results, "fs", fsym); + + % timing sync -> at this point we still have no symbol timing recovery, we + % try to do this with 2sps EQ! + [~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_matched.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1); + Rx_matched_1 = Rx_synced_cell{1}; + Rx_matched_original = Rx_synced_cell{1}; + Rx_matched_original.signal = resample(Rx_matched_original.signal,1,2); + + [~,Rx_synced_cell_tr_mf,inverted_tr_mf,sequenceFound_tr_mf,sequenceStarts_tr_mf] = data_tr_mf.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1); + Rx_tr_mf = Rx_synced_cell_tr_mf{1}; + + Time_Rec = 1; + if Time_Rec + % [Rx_Time_Rec, Timing_Error] = Timing_Recovery("modulation", 'PAM/PSK/QAM', "timing_error_detector",'Gardner (non-data-aided)','sps',2,'damping_factor',1,'normalized_loop_bandwidth',1e-4,'detector_gain',2.7).process(Rx_matched_1); + Rx_Time_Rec = Time_Shifter('value',time_shift).process(Rx_matched_1); + Rx_Time_Rec.signal = resample(Rx_Time_Rec.signal,1,2); + % figure(2222222) + % plot(Timing_Error) + sps = 1; + else + sps = 2; + end + + if Time_Rec == 1 && M == 2 + Rx_Time_Rec.fs = 14e9; + elseif Time_Rec == 1 && M == 4 + Rx_Time_Rec.fs = 6e9; + end + + Rx_Time_Rec.spectrum('normalizeTo0dB',1,"displayname",'Our signal','fignum',101111); + % Rx_matched_original.spectrum('normalizeTo0dB',1,"displayname",'Our signal','fignum',101111); + Rx_tr_mf.spectrum('normalizeTo0dB',1,"displayname",'Their signal','fignum',101111); + + Rx_Time_Rec_plot = Rx_Time_Rec; + % Rx_matched_original_plot = Rx_matched_original; + Rx_tr_mf_plot = Rx_tr_mf; + Rx_Time_Rec_plot.normalize("mode","rms").plot("displayname",'Our signal','fignum',101311); + % Rx_matched_original_plot.normalize("mode","rms").plot("displayname",'Original signal','fignum',101311); + Rx_tr_mf_plot.normalize("mode","rms").plot("displayname",'Their signal','fignum',101311); + + %% not working.. + % Use our or their signal + for our_signal = 1 + if our_signal + Rx_synced = Rx_Time_Rec; + else + Rx_synced = Rx_tr_mf; + end + % Rx_synced = Rx_Time_Rec; + % Rx_synced = Rx_synced_cell{1}; + len_tr = 4096*2; + mu_ffe1 = 0.0001; + mu_ffe2 = 0.0008; + mu_ffe3 = 0.001; + mu_dc = 0.004; + mu_ffe = [mu_ffe1 mu_ffe3 mu_ffe3]; + mu_dfe = 0.0004; + duob_mode = db_mode.no_db; + + Rx_synced.plot("displayname",'RX: Matched+Sync+2sps','fignum',103); + Rx_synced.spectrum("displayname",'RX: Matched+Sync+2sps','fignum',104); + + Digi_sig_compare.normalize("mode","rms").spectrum("displayname",'Tx: RC-shaped','fignum',1,'normalizeTo0dB',1); + Rx_sig_compare.normalize("mode","rms").spectrum("displayname",'Tx: RC-shaped + matched filtered ','fignum',1,'normalizeTo0dB',1); + Rx_synced.normalize("mode","rms").spectrum("displayname",'RX: matched filtered + synced','fignum',1,'normalizeTo0dB',1); + + if M == 2 + ber_in_paper = 10^(-2.6); %fig 3a) 4 Gb/s MWIR FSO Transmission using Directly Modulated QCL and an Uncooled UTC-PD at Room-Temperature + elseif M == 4 + ber_in_paper = 10^(-2.5); + end + + %% -------------------- FFE -------------------- + % % requires some more digging what is going on :-) + % eq_ffe = EQ("Ne",[500, 0, 0],"Nb",[0,0,0], ... + % "training_length",len_tr,"training_loops",5,"dd_loops",5, ... + % "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... + % "FFEmu",0,"plotfinal",0,"ideal_dfe",1); + % + % % eq_ffe = FFE_DFE('ffe_order',99,'dfe_order',99,'len_tr',len_tr,'epochs_tr',5,'epochs_dd',5,'sps',sps,'decide',0); + % + % ffe_results = ffe(eq_ffe,M,Rx_synced,Symbols,Bits, ... + % "precode_mode",duob_mode,'showAnalysis',0,"postFFE",[], ... + % "eth_style_symbol_mapping",mapping_style); + % + % if our_signal + % fprintf('Our signal: %.1e \n',ffe_results.metrics.BER); + % fprintf('Paper: %.1e \n \n',ber_in_paper); + % else + % fprintf('Their signal: %.1e \n',ffe_results.metrics.BER); + % fprintf('Paper: %.1e \n \n',ber_in_paper); + % end + + %% -------------------- VNLE + MLSE -------------------- + pf_ncoeffs = 4; + eq_v = EQ("Ne",[100, 0, 0],"Nb",[0, 0, 0], ... + "training_length",len_tr,"training_loops",5,"dd_loops",5, ... + "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... + "FFEmu",0,"plotfinal",0,"ideal_dfe",1, ... + 'weighted_DFE',0,'weighted_DFE_d_min',0.5,'weighted_DFE_mode','R2','weighted_DFE_I_mode',[5,0.5,0.6]); + % eq_v = FFE_DFE('ffe_order',200,'dfe_order',0,'len_tr',len_tr,'epochs_tr',5,'epochs_dd',5, ... + % 'ffe_mu_dd',mu_ffe,'ffe_mu_tr',0,'dfe_mu_dd',mu_dfe,'dfe_mu_tr',0.005,'sps',sps,'decide',0); + pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); + mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0,"eth_style",mapping_style).levels); + + [vnle_results, mlse_results] = vnle_postfilter_mlse(eq_v, pf_, mlse_, M, Rx_synced, Symbols, Bits, ... + "precode_mode", duob_mode, 'showAnalysis', 1, "postFFE", [], "eth_style_symbol_mapping", mapping_style); + + mlse_results.metrics.print + if our_signal + fprintf('Our Signal: %.1e \n',mlse_results.metrics.BER); + fprintf('Paper: %.1e \n \n',ber_in_paper); + else + fprintf('Their Signal: %.1e \n',mlse_results.metrics.BER); + fprintf('Paper: %.1e \n \n',ber_in_paper); + end + BER = mlse_results.metrics.BER; + + %% -------------------- ML-based MLSE (L=2) -------------------- + % ml_mlse_equalizer = ML_MLSE("epochs_tr",100,"epochs_dd",1, ... + % "len_tr",length(Rx_synced),"mu_dd",0.03,"mu_tr",0.03,"order",11,"sps",sps, ... + % "traceback_depth",256,"L",4,"delta",4,"adaptive_mu",0); + % + % [ml_mlse_results] = ml_mlse(ml_mlse_equalizer, M, Rx_synced, Symbols, Bits,"precode_mode",duob_mode,"eth_style_symbol_mapping",mapping_style); + % if our_signal + % fprintf('Our Signal: %.1e \n',ml_mlse_results.metrics.BER); + % fprintf('Paper: %.1e \n \n',ber_in_paper); + % else + % fprintf('Their EQ: %.1e \n',ml_mlse_results.metrics.BER); + % fprintf('Paper: %.1e \n \n',ber_in_paper); + % end + end +end \ No newline at end of file diff --git a/projects/FSO_transmission/Evaluation Scripts/first_analysis_tr_mf.m b/projects/FSO_transmission/Evaluation Scripts/first_analysis_tr_mf.m new file mode 100644 index 0000000..fe8f371 --- /dev/null +++ b/projects/FSO_transmission/Evaluation Scripts/first_analysis_tr_mf.m @@ -0,0 +1,271 @@ +%% +clear all; +close all; + +%% +base = "C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\FSO Equalizer\FSO_FP_QCL_60umUTC\"; +mode = 0; %0 oder 1 +M = 4; + +all_files = dir(fullfile(base, "**/*.mat")); + +% data_tr_mf is the already recovered and filtered data +if M == 2 + tx_data_path = fullfile(base, "14G_PAM2\tx_info\tx_info_PAM2_14Gbd0.75RRC.mat"); + filename = fullfile(base, "14G_PAM2\M=2_Rs=1.4e10_Fs=8e10_I=265mA_RoP=46.3mW_L=31m_PS=RRC_rolloff=0.75_Mode=Rise.mat"); + data_tr_mf = load("C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\FSO Equalizer\FSO_FP_QCL_60umUTC\Already Recovered and Filtered\AfterSync_M=2_Rs=1.4e10_Fs=8e10_I=265mA_RoP=46.3mW_L=31m_PS=RRC_rolloff=0.75_Mode=Rise.mat"); +elseif M == 4 + tx_data_path = fullfile(base, "6G_PAM4\tx_info\tx_info_PAM4_6Gbd0.6RRC.mat"); + filename = fullfile(base, "6G_PAM4\M=4_Rs=6e9_Fs=8e10_I=255mA_RoP=42.3mW_L=31m_PS=RRC_rolloff=0.6_Mode=Rise.mat"); + data_tr_mf = load("C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\FSO Equalizer\FSO_FP_QCL_60umUTC\Already Recovered and Filtered\AfterSync_M=4_Rs=6e9_Fs=8e10_I=255mA_RoP=42.3mW_L=31m_PS=RRC_rolloff=0.6_Mode=Rise.mat"); +end + +if mode == 1 + [f, p] = uigetfile(fullfile(base, "**/*.mat")); + if f~=0 + filename = fullfile(p,f); + end +end + +tx_data = load(tx_data_path); +datas = load(filename); + +%% +str = filename; +M_ = str2double(regexp(str, 'M=([^_]+)', 'tokens', 'once')); +assert(M==M_); +fsym = str2double(regexp(str, 'Rs=([^_]+)', 'tokens', 'once')); +fs = str2double(regexp(str, 'Fs=([^_]+)', 'tokens', 'once')); +I = sscanf(char(regexp(str, 'I=([^_]+)', 'tokens', 'once')), '%f'); +rop = sscanf(char(regexp(str, 'RoP=([^_]+)', 'tokens', 'once')), '%f'); +L = sscanf(char(regexp(str, 'L=([^_]+)', 'tokens', 'once')), '%f'); +pulseshape = string( regexp(str, 'PS=([^_]+)', 'tokens', 'once')); +rolloff = str2double(regexp(str, 'rolloff=([^_]+)', 'tokens', 'once')); +mode = string( regexp(str, 'Mode=([^\.]+)', 'tokens', 'once')); + +%% +% Tx data + +Bits = Informationsignal(tx_data.tx_data,"fs",fsym); +Symbols = Informationsignal(real(tx_data.tx_PAM_sym),"fs",fsym); + +mapping_style = M==4; % Pam2 is like move-it; PAM-4 is different, same mapping like ETH peopled used in Zurich... hence the "eth_style" argument here and there +PM = PAMmapper(M,0,"eth_style",mapping_style); % one should rename "eth style" as this is simply a different mapping scheme + +Symbols_ = PM.map(Bits) .* PM.scaling; +assert(isequal(Symbols.signal,Symbols_.signal)); + +Bits_ = PM.demap(Symbols); +[bits,errors,ber,errorIndice] = calc_ber(Bits_.signal,Bits.signal); +assert(ber == 0); + +%% For comparison, apply pulsef on Tx Symbols +Pform = Pulseformer("fsym",fsym,"fdac",fs,"pulse","rrc","pulselength",16,"alpha",rolloff); +Digi_sig_compare = Pform.process(Symbols); + +MF = Pulseformer("fsym",fsym,"fdac",2*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff); +Rx_sig_compare = MF.process(Digi_sig_compare); + +%% + +% Rx Data +traceData = datas.tr.lastData(2).trace.ch3; + +%FYI: Voltage=(RawData−YReference)×YIncrement+YOrigin +scoperead_volts = (traceData.RawData - traceData.YReference) * traceData.YIncrement + traceData.YOrigin; +demystified = isequal(traceData.YData,scoperead_volts); +assert(demystified); + +Scope_sig = Electricalsignal(traceData.YData,"fs",fs); + +Scope_sig.plot("displayname",'raw','fignum',100); +Scope_sig.spectrum("displayname",'raw','fignum',101) + +%Calculate Transfer Function +Tx_spectrum = Digi_sig_compare; +Rx_spectrum = Scope_sig; + +[~,Rx_synced_spectrum,inverted_spectrum,sequenceFound_spectrum,sequenceStarts_spectrum] = Rx_spectrum.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1); +Rx_spectrum = Rx_synced_spectrum{1}; + +Tx_spectrum = Tx_spectrum.resample("fs_out",fsym); +Tx_spectrum.signal = fft(Tx_spectrum.signal); + +Rx_spectrum = Rx_spectrum.resample("fs_out",fsym); +Rx_spectrum.signal = fft(Rx_spectrum.signal); + +H_transfer = Rx_spectrum.signal./Tx_spectrum.signal; +H_inv = 1./H_transfer; + +%Number of Samples/Symbol after Matched Filter +Kov = 14; + +Scope_sig = Scope_sig.resample('fs_in',fs,'fs_out',Kov*fsym); + +% 1) matched filter +% pulse is symmetric, hence we can use pulsef firectly as matched filter. +% It feels off (bit I think correct) that the fsym is now the output freq.!! +% -> output 2 sps to omit timing recovery!? +Matched_Filter = Pulseformer("fsym",fsym,"fdac",Kov*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff,"matched",1); +Rx_matched = Matched_Filter.process(Scope_sig); +Rx_matched.spectrum("displayname",'Signal after matched filter','fignum',1); +% Rx_matched = Scope_sig; +% Rx_matched = Rx_matched.resample('fs_out',Kov*fsym); + +data_tr_mf = Electricalsignal(data_tr_mf.Results, "fs", fsym); +[~,Rx_synced_cell_tr_mf,inverted_tr_mf,sequenceFound_tr_mf,sequenceStarts_tr_mf] = data_tr_mf.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1); +Rx_tr_mf = Rx_synced_cell_tr_mf{1}; + +% timing sync -> at this point we still have no symbol timing recovery, we +% try to do this with 2sps EQ! + +% Rx_matched = Scope_sig; +% Rx_matched = Rx_matched.resample("fs_out",2*fsym); + +[~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_matched.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1); +Rx_matched_1 = Rx_synced_cell{1}; + +Rx_matched_original = Rx_synced_cell{1}; +Rx_matched_original.signal = resample(Rx_matched_original.signal,1,Kov); +Rx_matched_original.fs = fsym; + +% Rx_matched_1 = Rx_matched_1.resample("fs_in",Kov*fsym,'fs_out',2*fsym); + +Time_Rec = 1; +if Time_Rec + % [Rx_Time_Rec, Timing_Error] = Timing_Recovery("modulation", 'PAM/PSK/QAM', "timing_error_detector",'Gardner (non-data-aided)','sps',Kov,'damping_factor',1,'normalized_loop_bandwidth',1e-4,'detector_gain',2.7).process(Rx_matched_1); + + [Rx_Time_Rec, Timing_Error_MG] = Godard_Timing_Recovery('mode',3,'num_blocks',1,'fft_length',length(Rx_matched_1),'sps',Kov,'rolloff',0.6,'mu',-0.2,'Ki',1e-4).process(Rx_matched_1); + Rx_Time_Rec = Rx_Time_Rec.resample('fs_in',Kov*fsym,'fs_out',fsym); + + % Rx_Time_Rec = MaxVar_Timing_Recovery('mode',0,'fsym',fsym,'fadc',Kov*fsym,'num_tau',Kov*128,'sps',Kov,'comp_signal',Rx_tr_mf,'comp_mode',0).process(Rx_matched_1); + + sps = 1; +else + Rx_Time_Rec = Rx_matched_1; + sps = 1; +end + +Rx_Time_Rec.fs = fsym; + +% h = gaussdesign(0.4,64,1); +% Rx_Time_Rec.signal = filtfilt(h,1,Rx_Time_Rec.signal); + +% Rx_com = Rx_Time_Rec; +% Rx_com = Rx_com.resample('fs_in',fsym,'fs_out',2*fsym); +% Rx_com.spectrum("displayname",'After TR','normalizeTo0dB',1,'fignum',101114); +% Rx_matched_1.spectrum("displayname",'Before TR','normalizeTo0dB',1,'fignum',101114); +% Rx_tr_mf_2sps = Rx_tr_mf; +% Rx_tr_mf_2sps = Rx_tr_mf_2sps.resample('fs_in',fsym,'fs_out',2*fsym); +% Rx_tr_mf_2sps.spectrum("displayname",'Their signal after TR','normalizeTo0dB',1,'fignum',101114); + +% Rx_Time_Rec = Rx_Time_Rec.resample('fs_out',6e9); +% Rx_tr_mf = Rx_tr_mf.resample('fs_out',fsym); + +% Amax = 10; +% H_inv = min(abs(H_inv), Amax) .* exp(1j*angle(H_inv)); +% Rx_Time_Rec.signal = fft(Rx_Time_Rec.signal); +% Rx_Time_Rec.signal = ifft(Rx_Time_Rec.signal .* H_inv); + +Rx_Time_Rec = Rx_Time_Rec.normalize('mode','rms'); +Rx_tr_mf = Rx_tr_mf.normalize('mode','rms'); + +Rx_Time_Rec.spectrum("displayname",'Our signal','normalizeTo0dB',1,'fignum',101111); +% Rx_matched_original.spectrum('normalizeTo0dB',1,"displayname",'Our signal','fignum',101111); +Rx_tr_mf.spectrum("displayname",'Their signal','normalizeTo0dB',1,'fignum',101111); + +Rx_Time_Rec.normalize("mode","rms").plot("displayname",'Our signal','fignum',101311); +% Rx_matched_original.normalize("mode","rms").plot("displayname",'Original signal','fignum',101311); +Rx_tr_mf.normalize("mode","rms").plot("displayname",'Their signal','fignum',101311); + +%% not working.. +% Use our or their signal +for our_signal = 1 + if our_signal + Rx_synced = Rx_Time_Rec; + else + Rx_synced = Rx_tr_mf; + end + % Rx_synced = Rx_Time_Rec; + % Rx_synced = Rx_synced_cell{1}; + len_tr = 4096*2; + mu_ffe1 = 0.0001; + mu_ffe2 = 0.0008; + mu_ffe3 = 0.001; + mu_dc = 0.004; + mu_ffe = [mu_ffe1 mu_ffe3 mu_ffe3]; + mu_dfe = 0.0004; + duob_mode = db_mode.no_db; + + Rx_synced.plot("displayname",'RX: Matched+Sync+2sps','fignum',103); + Rx_synced.spectrum("displayname",'RX: Matched+Sync+2sps','fignum',104); + + Digi_sig_compare.normalize("mode","rms").spectrum("displayname",'Tx: RC-shaped','fignum',1,'normalizeTo0dB',1); + Rx_sig_compare.normalize("mode","rms").spectrum("displayname",'Tx: RC-shaped + matched filtered ','fignum',1,'normalizeTo0dB',1); + Rx_synced.normalize("mode","rms").spectrum("displayname",'RX: matched filtered + synced','fignum',1,'normalizeTo0dB',1); + + if M == 2 + ber_in_paper = 10^(-2.6); %fig 3a) 4 Gb/s MWIR FSO Transmission using Directly Modulated QCL and an Uncooled UTC-PD at Room-Temperature + elseif M == 4 + ber_in_paper = 10^(-2.5); + end + + %% -------------------- FFE -------------------- + % requires some more digging what is going on :-) + % eq_ffe = EQ("Ne",[150, 0, 0],"Nb",[0,0,0], ... + % "training_length",len_tr,"training_loops",5,"dd_loops",5, ... + % "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... + % "FFEmu",0,"plotfinal",0,"ideal_dfe",1); + % + % % eq_ffe = FFE_DFE('ffe_order',99,'dfe_order',99,'len_tr',len_tr,'epochs_tr',5,'epochs_dd',5,'sps',sps,'decide',0); + % + % ffe_results = ffe(eq_ffe,M,Rx_synced,Symbols,Bits, ... + % "precode_mode",duob_mode,'showAnalysis',0,"postFFE",[], ... + % "eth_style_symbol_mapping",mapping_style); + % + % if our_signal + % fprintf('Our signal: %.1e \n',ffe_results.metrics.BER); + % fprintf('Paper: %.1e \n \n',ber_in_paper); + % else + % fprintf('Their signal: %.1e \n',ffe_results.metrics.BER); + % fprintf('Paper: %.1e \n \n',ber_in_paper); + % end + + %% -------------------- VNLE + MLSE -------------------- + pf_ncoeffs = 4; + eq_v = EQ("Ne",[300, 0, 0],"Nb",[0, 0, 0], ... + "training_length",len_tr,"training_loops",5,"dd_loops",5, ... + "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... + "FFEmu",0,"plotfinal",0,"ideal_dfe",1, ... + 'weighted_DFE',0,'weighted_DFE_d_min',0.5,'weighted_DFE_mode','R2','weighted_DFE_I_mode',[5,0.5,0.6]); + % eq_v = FFE_DFE('ffe_order',300,'dfe_order',5,'len_tr',len_tr,'epochs_tr',5,'epochs_dd',5, ... + % 'ffe_mu_dd',mu_ffe,'ffe_mu_tr',0,'dfe_mu_dd',mu_dfe,'dfe_mu_tr',0.005,'sps',sps,'decide',0); + pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); + mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0,"eth_style",mapping_style).levels); + + [vnle_results, mlse_results] = vnle_postfilter_mlse(eq_v, pf_, mlse_, M, Rx_synced, Symbols, Bits, ... + "precode_mode", duob_mode, 'showAnalysis', 1, "postFFE", [], "eth_style_symbol_mapping", mapping_style); + + mlse_results.metrics.print + if our_signal + fprintf('Our Signal: %.1e \n',mlse_results.metrics.BER); + fprintf('Paper: %.1e \n \n',ber_in_paper); + else + fprintf('Their Signal: %.1e \n',mlse_results.metrics.BER); + fprintf('Paper: %.1e \n \n',ber_in_paper); + end + + %% -------------------- ML-based MLSE (L=2) -------------------- + % ml_mlse_equalizer = ML_MLSE("epochs_tr",100,"epochs_dd",1, ... + % "len_tr",length(Rx_synced),"mu_dd",0.03,"mu_tr",0.03,"order",11,"sps",sps, ... + % "traceback_depth",256,"L",4,"delta",4,"adaptive_mu",0); + % + % [ml_mlse_results] = ml_mlse(ml_mlse_equalizer, M, Rx_synced, Symbols, Bits,"precode_mode",duob_mode,"eth_style_symbol_mapping",mapping_style); + % if our_signal + % fprintf('Our Signal: %.1e \n',ml_mlse_results.metrics.BER); + % fprintf('Paper: %.1e \n \n',ber_in_paper); + % else + % fprintf('Their EQ: %.1e \n',ml_mlse_results.metrics.BER); + % fprintf('Paper: %.1e \n \n',ber_in_paper); + % end +end \ No newline at end of file diff --git a/projects/FSO_transmission/Evaluation Scripts/first_analysis_tr_mf_pam2.m b/projects/FSO_transmission/Evaluation Scripts/first_analysis_tr_mf_pam2.m new file mode 100644 index 0000000..e6f4c45 --- /dev/null +++ b/projects/FSO_transmission/Evaluation Scripts/first_analysis_tr_mf_pam2.m @@ -0,0 +1,195 @@ +%% +base = "C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\FSO Equalizer\FSO_FP_QCL_60umUTC\"; +mode = 0; %0 oder 1 +M = 2; + +all_files = dir(fullfile(base, "**/*.mat")); + +if M == 2 + tx_data_path = fullfile(base, "14G_PAM2\tx_info\tx_info_PAM2_14Gbd0.75RRC.mat"); + filename = fullfile(base, "14G_PAM2\M=2_Rs=1.4e10_Fs=8e10_I=265mA_RoP=46.3mW_L=31m_PS=RRC_rolloff=0.75_Mode=Rise.mat"); + % Already recovered and filtered data + data_tr_mf = load("C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\FSO Equalizer\FSO_FP_QCL_60umUTC\Already Recovered and Filtered\AfterSync_M=2_Rs=1.4e10_Fs=8e10_I=265mA_RoP=46.3mW_L=31m_PS=RRC_rolloff=0.75_Mode=Rise.mat"); +elseif M == 4 + tx_data_path = fullfile(base, "6G_PAM4\tx_info\tx_info_PAM4_6Gbd0.6RRC.mat"); + filename = fullfile(base, "6G_PAM4\M=4_Rs=6e9_Fs=8e10_I=255mA_RoP=42.3mW_L=31m_PS=RRC_rolloff=0.6_Mode=Rise.mat"); + data_tr_mf = load("C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\FSO Equalizer\FSO_FP_QCL_60umUTC\Already Recovered and Filtered\AfterSync_M=4_Rs=6e9_Fs=8e10_I=255mA_RoP=42.3mW_L=31m_PS=RRC_rolloff=0.6_Mode=Rise.mat"); +end + +if mode == 1 + [f, p] = uigetfile(fullfile(base, "**/*.mat")); + if f~=0 + filename = fullfile(p,f); + end +end + +tx_data = load(tx_data_path); +datas = load(filename); + +%% +str = filename; +M_ = str2double(regexp(str, 'M=([^_]+)', 'tokens', 'once')); +assert(M==M_); +fsym = str2double(regexp(str, 'Rs=([^_]+)', 'tokens', 'once')); +fs = str2double(regexp(str, 'Fs=([^_]+)', 'tokens', 'once')); +I = sscanf(char(regexp(str, 'I=([^_]+)', 'tokens', 'once')), '%f'); +rop = sscanf(char(regexp(str, 'RoP=([^_]+)', 'tokens', 'once')), '%f'); +L = sscanf(char(regexp(str, 'L=([^_]+)', 'tokens', 'once')), '%f'); +pulseshape = string( regexp(str, 'PS=([^_]+)', 'tokens', 'once')); +rolloff = str2double(regexp(str, 'rolloff=([^_]+)', 'tokens', 'once')); +mode = string( regexp(str, 'Mode=([^\.]+)', 'tokens', 'once')); + +%% +% Tx data + +Bits = Informationsignal(tx_data.tx_data,"fs",fsym); +Symbols = Informationsignal(real(tx_data.tx_PAM_sym),"fs",fsym); + +mapping_style = M==4; % Pam2 is like move-it; PAM-4 is different, same mapping like ETH peopled used in Zurich... hence the "eth_style" argument here and there +PM = PAMmapper(M,0,"eth_style",mapping_style); % one should rename "eth style" as this is simply a different mapping scheme + +Symbols_ = PM.map(Bits) .* PM.scaling; +assert(isequal(Symbols.signal,Symbols_.signal)); + +Bits_ = PM.demap(Symbols); +[bits,errors,ber,errorIndice] = calc_ber(Bits_.signal,Bits.signal); +assert(ber == 0); + +%% For comparison, apply pulsef on Tx Symbols +Pform = Pulseformer("fsym",fsym,"fdac",fs,"pulse","rrc","pulselength",16,"alpha",rolloff); +Digi_sig_compare = Pform.process(Symbols); +MF = Pulseformer("fsym",fsym,"fdac",2*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff); +Rx_sig_compare = MF.process(Digi_sig_compare); + +%% + +% Rx Data +traceData = datas.tr.lastData(2).trace.ch3; + +%FYI: Voltage=(RawData−YReference)×YIncrement+YOrigin +scoperead_volts = (traceData.RawData - traceData.YReference) * traceData.YIncrement + traceData.YOrigin; +demystified = isequal(traceData.YData,scoperead_volts); +assert(demystified); + +Scope_sig = Electricalsignal(traceData.YData,"fs",fs); + +Scope_sig.plot("displayname",'raw','fignum',100); +Scope_sig.spectrum("displayname",'raw','fignum',101) + +% 1) matched filter +% pulse is symmetric, hence we can use pulsef firectly as matched filter. +% It feels off (bit I think correct) that the fsym is now the output freq.!! +% -> output 2 sps to omit timing recovery!? +Matched_Filter = Pulseformer("fsym",fsym,"fdac",2*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff,"matched",1); +Rx_matched = Matched_Filter.process(Scope_sig); +Rx_matched.spectrum("displayname",'Signal after matched filter','fignum',1); + +data_tr_mf = Electricalsignal(data_tr_mf.Results, "fs", fsym); + +% timing sync -> at this point we still have no symbol timing recovery, we +% try to do this with 2sps EQ! +[~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_matched.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1); +Rx_matched_1 = Rx_synced_cell{1}; + +[~,Rx_synced_cell_tr_mf,inverted_tr_mf,sequenceFound_tr_mf,sequenceStarts_tr_mf] = data_tr_mf.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1); +Rx_tr_mf = Rx_synced_cell_tr_mf{1}; + + +Time_Rec = 1; +if Time_Rec + [Rx_Time_Rec, Timing_Error] = Timing_Recovery("modulation", 'PAM/PSK/QAM', "timing_error_detector",'Gardner (non-data-aided)','sps',2,'damping_factor',1,'normalized_loop_bandwidth',0.01,'detector_gain',2.7).process(Rx_matched_1); + sps = 1; +else + sps = 2; +end + +if Time_Rec == 1 && M == 2 + Rx_Time_Rec.fs = 14e9; +elseif Time_Rec == 1 && M == 4 + Rx_Time_Rec.fs = 6e9; +end + +Rx_Time_Rec.spectrum("displayname",'Signal after matched filter, synchronization, and timing recovery','fignum',101111); +Rx_tr_mf.spectrum("displayname",'Signal after matched filter, synchronization, and timing recovery','fignum',101111); + +Rx_Time_Rec_plot = Rx_Time_Rec; +Rx_tr_mf_plot = Rx_tr_mf; +Rx_Time_Rec_plot.normalize("mode","rms").plot("displayname",'Signal after matched filter, synchronization, and timing recovery','fignum',101311); +Rx_tr_mf_plot.normalize("mode","rms").plot("displayname",'Signal after matched filter, synchronization, and timing recovery','fignum',101311); + +%% not working.. +% Use our or their signal +our_signal = 0; +if our_signal + Rx_synced = Rx_Time_Rec; +else + Rx_synced = Rx_tr_mf; +end +% Rx_synced = Rx_Time_Rec; +% Rx_synced = Rx_synced_cell{1}; +len_tr = 4096*2; +mu_ffe1 = 0.0001; +mu_ffe2 = 0.0008; +mu_ffe3 = 0.001; +mu_dc = 0.004; +mu_ffe = [mu_ffe1 mu_ffe3 mu_ffe3]; +mu_dfe = 0.0004; +duob_mode = db_mode.no_db; +sps = 1; + +Rx_synced.plot("displayname",'RX: Matched+Sync+2sps','fignum',103); +Rx_synced.spectrum("displayname",'RX: Matched+Sync+2sps','fignum',104); + +Digi_sig_compare.normalize("mode","rms").spectrum("displayname",'Tx: RC-shaped','fignum',1,'normalizeTo0dB',1); +Rx_sig_compare.normalize("mode","rms").spectrum("displayname",'Tx: RC-shaped + matched filtered ','fignum',1,'normalizeTo0dB',1); +Rx_synced.normalize("mode","rms").spectrum("displayname",'RX: matched filtered + synced','fignum',1,'normalizeTo0dB',1); + +if M == 2 + ber_in_paper = 10^(-2.6); %fig 3a) 4 Gb/s MWIR FSO Transmission using Directly Modulated QCL and an Uncooled UTC-PD at Room-Temperature +elseif M == 4 + ber_in_paper = 10^(-2.5); +end + +%% -------------------- FFE -------------------- +% % requires some more digging what is going on :-) +% eq_ffe = EQ("Ne",[200, 0, 0],"Nb",[0,0,0], ... +% "training_length",len_tr,"training_loops",5,"dd_loops",5, ... +% "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... +% "FFEmu",0,"plotfinal",0,"ideal_dfe",1); +% +% ffe_results = ffe(eq_ffe,M,Rx_synced,Symbols,Bits, ... +% "precode_mode",duob_mode,'showAnalysis',0,"postFFE",[], ... +% "eth_style_symbol_mapping",mapping_style); +% +% % ffe_results.metrics.print +% fprintf('My EQ: %.1e \n',ffe_results.metrics.BER); +% fprintf('Paper: %.1e \n \n',ber_in_paper); +% BER_value = ffe_results.metrics.BER; + +%% -------------------- VNLE + MLSE -------------------- + +pf_ncoeffs = 4; +eq_v = EQ("Ne",[200, 0, 0],"Nb",[0, 0, 0], ... + "training_length",len_tr,"training_loops",5,"dd_loops",5, ... + "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... + "FFEmu",0,"plotfinal",0,"ideal_dfe",1, ... + 'weighted_DFE',0,'weighted_DFE_d_min',0.5,'weighted_DFE_mode','I2','weighted_DFE_I_mode',[5,0.5,0.6]); +pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); +mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0,"eth_style",mapping_style).levels); + +[vnle_results, mlse_results] = vnle_postfilter_mlse(eq_v, pf_, mlse_, M, Rx_synced, Symbols, Bits, ... + "precode_mode", duob_mode, 'showAnalysis', 1, "postFFE", [], "eth_style_symbol_mapping", mapping_style); + +mlse_results.metrics.print +fprintf('My EQ: %.1e \n',mlse_results.metrics.BER); +fprintf('Paper: %.1e \n \n',ber_in_paper); + +%% -------------------- ML-based MLSE (L=2) -------------------- +% ml_mlse_equalizer = ML_MLSE("epochs_tr",400,"epochs_dd",1, ... +% "len_tr",len_tr,"mu_dd",0.03,"mu_tr",0.03,"order",100,"sps",sps, ... +% "traceback_depth",256,"L",2,"delta",4,"adaptive_mu",0); +% +% [ml_mlse_results] = ml_mlse(ml_mlse_equalizer, M, Rx_synced, Symbols, Bits,"precode_mode",duob_mode,"eth_style_symbol_mapping",mapping_style); +% fprintf('ML-based MLSE:\n'); +% fprintf('My EQ: %.1e \n',ml_mlse_results.metrics.BER); +% fprintf('Paper: %.1e \n \n',ber_in_paper); \ No newline at end of file diff --git a/projects/FSO_transmission/Evaluation Scripts/first_analysis_tr_mf_pam_2.m b/projects/FSO_transmission/Evaluation Scripts/first_analysis_tr_mf_pam_2.m new file mode 100644 index 0000000..d27d6e1 --- /dev/null +++ b/projects/FSO_transmission/Evaluation Scripts/first_analysis_tr_mf_pam_2.m @@ -0,0 +1,272 @@ +%% +clear all; +close all; + +%% +base = "C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\FSO Equalizer\FSO_FP_QCL_60umUTC\"; +mode = 0; %0 oder 1 +M = 2; + +all_files = dir(fullfile(base, "**/*.mat")); + +% data_tr_mf is the already recovered and filtered data +if M == 2 + tx_data_path = fullfile(base, "14G_PAM2\tx_info\tx_info_PAM2_14Gbd0.75RRC.mat"); + filename = fullfile(base, "14G_PAM2\M=2_Rs=1.4e10_Fs=8e10_I=265mA_RoP=46.3mW_L=31m_PS=RRC_rolloff=0.75_Mode=Rise.mat"); + data_tr_mf = load("C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\FSO Equalizer\FSO_FP_QCL_60umUTC\Already Recovered and Filtered\AfterSync_M=2_Rs=1.4e10_Fs=8e10_I=265mA_RoP=46.3mW_L=31m_PS=RRC_rolloff=0.75_Mode=Rise.mat"); +elseif M == 4 + tx_data_path = fullfile(base, "6G_PAM4\tx_info\tx_info_PAM4_6Gbd0.6RRC.mat"); + filename = fullfile(base, "6G_PAM4\M=4_Rs=6e9_Fs=8e10_I=255mA_RoP=42.3mW_L=31m_PS=RRC_rolloff=0.6_Mode=Rise.mat"); + data_tr_mf = load("C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\FSO Equalizer\FSO_FP_QCL_60umUTC\Already Recovered and Filtered\AfterSync_M=4_Rs=6e9_Fs=8e10_I=255mA_RoP=42.3mW_L=31m_PS=RRC_rolloff=0.6_Mode=Rise.mat"); +end + +if mode == 1 + [f, p] = uigetfile(fullfile(base, "**/*.mat")); + if f~=0 + filename = fullfile(p,f); + end +end + +tx_data = load(tx_data_path); +datas = load(filename); + +%% +str = filename; +M_ = str2double(regexp(str, 'M=([^_]+)', 'tokens', 'once')); +assert(M==M_); +fsym = str2double(regexp(str, 'Rs=([^_]+)', 'tokens', 'once')); +fs = str2double(regexp(str, 'Fs=([^_]+)', 'tokens', 'once')); +I = sscanf(char(regexp(str, 'I=([^_]+)', 'tokens', 'once')), '%f'); +rop = sscanf(char(regexp(str, 'RoP=([^_]+)', 'tokens', 'once')), '%f'); +L = sscanf(char(regexp(str, 'L=([^_]+)', 'tokens', 'once')), '%f'); +pulseshape = string( regexp(str, 'PS=([^_]+)', 'tokens', 'once')); +rolloff = str2double(regexp(str, 'rolloff=([^_]+)', 'tokens', 'once')); +mode = string( regexp(str, 'Mode=([^\.]+)', 'tokens', 'once')); + +%% +% Tx data + +Bits = Informationsignal(tx_data.tx_data,"fs",fsym); +Symbols = Informationsignal(real(tx_data.tx_PAM_sym),"fs",fsym); + +mapping_style = M==4; % Pam2 is like move-it; PAM-4 is different, same mapping like ETH peopled used in Zurich... hence the "eth_style" argument here and there +PM = PAMmapper(M,0,"eth_style",mapping_style); % one should rename "eth style" as this is simply a different mapping scheme + +Symbols_ = PM.map(Bits) .* PM.scaling; +assert(isequal(Symbols.signal,Symbols_.signal)); + +Bits_ = PM.demap(Symbols); +[bits,errors,ber,errorIndice] = calc_ber(Bits_.signal,Bits.signal); +assert(ber == 0); + +%% For comparison, apply pulsef on Tx Symbols +Pform = Pulseformer("fsym",fsym,"fdac",fs,"pulse","rrc","pulselength",16,"alpha",rolloff); +Digi_sig_compare = Pform.process(Symbols); + +MF = Pulseformer("fsym",fsym,"fdac",2*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff); +Rx_sig_compare = MF.process(Digi_sig_compare); + +%% + +% Rx Data +traceData = datas.tr.lastData(2).trace.ch3; + +%FYI: Voltage=(RawData−YReference)×YIncrement+YOrigin +scoperead_volts = (traceData.RawData - traceData.YReference) * traceData.YIncrement + traceData.YOrigin; +demystified = isequal(traceData.YData,scoperead_volts); +assert(demystified); + +Scope_sig = Electricalsignal(traceData.YData,"fs",fs); + +Scope_sig.plot("displayname",'raw','fignum',100); +Scope_sig.spectrum("displayname",'raw','fignum',101) + +%Calculate Transfer Function +Tx_spectrum = Digi_sig_compare; +Rx_spectrum = Scope_sig; + +[~,Rx_synced_spectrum,inverted_spectrum,sequenceFound_spectrum,sequenceStarts_spectrum] = Rx_spectrum.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1); +Rx_spectrum = Rx_synced_spectrum{1}; + +Tx_spectrum = Tx_spectrum.resample("fs_out",fsym); +Tx_spectrum.signal = fft(Tx_spectrum.signal); + +Rx_spectrum = Rx_spectrum.resample("fs_out",fsym); +Rx_spectrum.signal = fft(Rx_spectrum.signal); + +H_transfer = Rx_spectrum.signal./Tx_spectrum.signal; +H_inv = 1./H_transfer; + +%Number of Samples/Symbol after Matched Filter +Kov = 14; + +Scope_sig = Scope_sig.resample('fs_in',fs,'fs_out',Kov*fsym); + +% 1) matched filter +% pulse is symmetric, hence we can use pulsef firectly as matched filter. +% It feels off (bit I think correct) that the fsym is now the output freq.!! +% -> output 2 sps to omit timing recovery!? +Matched_Filter = Pulseformer("fsym",fsym,"fdac",Kov*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff,"matched",1); +Rx_matched = Matched_Filter.process(Scope_sig); +Rx_matched.spectrum("displayname",'Signal after matched filter','fignum',1); +% Rx_matched = Scope_sig; +% Rx_matched = Rx_matched.resample('fs_out',Kov*fsym); + +data_tr_mf = Electricalsignal(data_tr_mf.Results, "fs", fsym); +[~,Rx_synced_cell_tr_mf,inverted_tr_mf,sequenceFound_tr_mf,sequenceStarts_tr_mf] = data_tr_mf.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1); +Rx_tr_mf = Rx_synced_cell_tr_mf{1}; + +% timing sync -> at this point we still have no symbol timing recovery, we +% try to do this with 2sps EQ! + +% Rx_matched = Scope_sig; +% Rx_matched = Rx_matched.resample("fs_out",2*fsym); + +[~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_matched.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1); +Rx_matched_1 = Rx_synced_cell{1}; + +Rx_matched_original = Rx_synced_cell{1}; +Rx_matched_original.signal = resample(Rx_matched_original.signal,1,Kov); +Rx_matched_original.fs = fsym; + +% Rx_matched_1 = Rx_matched_1.resample("fs_in",Kov*fsym,'fs_out',2*fsym); + +Time_Rec = 1; +if Time_Rec + % [Rx_Time_Rec, Timing_Error] = Timing_Recovery("modulation", 'PAM/PSK/QAM', "timing_error_detector",'Gardner (non-data-aided)','sps',Kov,'damping_factor',1,'normalized_loop_bandwidth',1e-4,'detector_gain',2.7).process(Rx_matched_1); + % Rx_matched_1.spectrum('normalizeTo0dB',1,"displayname",'Signal before Timing Recovery','fignum',1023); + % [Rx_Time_Rec, Timing_Error_MG] = Godard_Timing_Recovery('mode',3,'num_blocks',1,'fft_length',length(Rx_matched_1),'sps',Kov,'rolloff',0.6,'mu',-0.2,'Ki',1e-4).process(Rx_matched_1); + Rx_Time_Rec = MaxVar_Timing_Recovery('mode',0,'fsym',fsym,'fadc',Kov*fsym,'num_tau',Kov*16,'sps',Kov,'comp_signal',Rx_tr_mf,'comp_mode',0).process(Rx_matched_1); + % Rx_Time_Rec.spectrum('normalizeTo0dB',1,"displayname",'Signal after Timing Recovery','fignum',1023); + % idx = 1:2:length(Rx_Time_Rec); + % Rx_Time_Rec.signal = Rx_Time_Rec.signal(idx); + % Rx_Time_Rec.signal = resample(Rx_Time_Rec.signal,1,Kov); + % Rx_Time_Rec = Rx_matched_1; + % Rx_Time_Rec.signal = resample(Rx_Time_Rec.signal,1,Kov); + % figure(2222222) + % plot(Timing_Error) + sps = 1; +else + Rx_Time_Rec = Rx_matched_1; + sps = 1; +end +Rx_Time_Rec.fs = fsym; + +Rx_com = Rx_Time_Rec; +Rx_com = Rx_com.resample('fs_in',fsym,'fs_out',2*fsym); +Rx_com.spectrum("displayname",'After TR','normalizeTo0dB',1,'fignum',101114); +Rx_matched_1.spectrum("displayname",'Before TR','normalizeTo0dB',1,'fignum',101114); +Rx_tr_mf_2sps = Rx_tr_mf; +Rx_tr_mf_2sps = Rx_tr_mf_2sps.resample('fs_in',fsym,'fs_out',2*fsym); +Rx_tr_mf_2sps.spectrum("displayname",'Their signal after TR','normalizeTo0dB',1,'fignum',101114); + +% Rx_Time_Rec = Rx_Time_Rec.resample('fs_out',6e9); +% Rx_tr_mf = Rx_tr_mf.resample('fs_out',fsym); + +% Amax = 10; +% H_inv = min(abs(H_inv), Amax) .* exp(1j*angle(H_inv)); +% Rx_Time_Rec.signal = fft(Rx_Time_Rec.signal); +% Rx_Time_Rec.signal = ifft(Rx_Time_Rec.signal .* H_inv); + +Rx_Time_Rec = Rx_Time_Rec.normalize('mode','rms'); +Rx_tr_mf = Rx_tr_mf.normalize('mode','rms'); + +Rx_Time_Rec.spectrum('normalizeTo0dB',1,"displayname",'Our signal','fignum',101111); +% Rx_matched_original.spectrum('normalizeTo0dB',1,"displayname",'Our signal','fignum',101111); +Rx_tr_mf.spectrum('normalizeTo0dB',1,"displayname",'Their signal','fignum',101111); + +Rx_Time_Rec.normalize("mode","rms").plot("displayname",'Our signal','fignum',101311); +% Rx_matched_original.normalize("mode","rms").plot("displayname",'Original signal','fignum',101311); +Rx_tr_mf.normalize("mode","rms").plot("displayname",'Their signal','fignum',101311); + +%% not working.. +% Use our or their signal +for our_signal = 1 + if our_signal + Rx_synced = Rx_Time_Rec; + else + Rx_synced = Rx_tr_mf; + end + % Rx_synced = Rx_Time_Rec; + % Rx_synced = Rx_synced_cell{1}; + len_tr = 4096*2; + mu_ffe1 = 0.0001; + mu_ffe2 = 0.0008; + mu_ffe3 = 0.001; + mu_dc = 0.004; + mu_ffe = [mu_ffe1 mu_ffe3 mu_ffe3]; + mu_dfe = 0.0004; + duob_mode = db_mode.no_db; + + Rx_synced.plot("displayname",'RX: Matched+Sync+2sps','fignum',103); + Rx_synced.spectrum("displayname",'RX: Matched+Sync+2sps','fignum',104); + + Digi_sig_compare.normalize("mode","rms").spectrum("displayname",'Tx: RC-shaped','fignum',1,'normalizeTo0dB',1); + Rx_sig_compare.normalize("mode","rms").spectrum("displayname",'Tx: RC-shaped + matched filtered ','fignum',1,'normalizeTo0dB',1); + Rx_synced.normalize("mode","rms").spectrum("displayname",'RX: matched filtered + synced','fignum',1,'normalizeTo0dB',1); + + if M == 2 + ber_in_paper = 10^(-2.6); %fig 3a) 4 Gb/s MWIR FSO Transmission using Directly Modulated QCL and an Uncooled UTC-PD at Room-Temperature + elseif M == 4 + ber_in_paper = 10^(-2.5); + end + + %% -------------------- FFE -------------------- + % % requires some more digging what is going on :-) + % eq_ffe = EQ("Ne",[500, 0, 0],"Nb",[0,0,0], ... + % "training_length",len_tr,"training_loops",5,"dd_loops",5, ... + % "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... + % "FFEmu",0,"plotfinal",0,"ideal_dfe",1); + % + % % eq_ffe = FFE_DFE('ffe_order',99,'dfe_order',99,'len_tr',len_tr,'epochs_tr',5,'epochs_dd',5,'sps',sps,'decide',0); + % + % ffe_results = ffe(eq_ffe,M,Rx_synced,Symbols,Bits, ... + % "precode_mode",duob_mode,'showAnalysis',0,"postFFE",[], ... + % "eth_style_symbol_mapping",mapping_style); + % + % if our_signal + % fprintf('Our signal: %.1e \n',ffe_results.metrics.BER); + % fprintf('Paper: %.1e \n \n',ber_in_paper); + % else + % fprintf('Their signal: %.1e \n',ffe_results.metrics.BER); + % fprintf('Paper: %.1e \n \n',ber_in_paper); + % end + + %% -------------------- VNLE + MLSE -------------------- + pf_ncoeffs = 4; + eq_v = EQ("Ne",[250, 0, 0],"Nb",[0, 0, 0], ... + "training_length",len_tr,"training_loops",5,"dd_loops",5, ... + "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... + "FFEmu",0,"plotfinal",0,"ideal_dfe",1, ... + 'weighted_DFE',0,'weighted_DFE_d_min',0.5,'weighted_DFE_mode','R2','weighted_DFE_I_mode',[5,0.5,0.6]); + % eq_v = FFE_DFE('ffe_order',200,'dfe_order',0,'len_tr',len_tr,'epochs_tr',5,'epochs_dd',5, ... + % 'ffe_mu_dd',mu_ffe,'ffe_mu_tr',0,'dfe_mu_dd',mu_dfe,'dfe_mu_tr',0.005,'sps',sps,'decide',0); + pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); + mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0,"eth_style",mapping_style).levels); + + [vnle_results, mlse_results] = vnle_postfilter_mlse(eq_v, pf_, mlse_, M, Rx_synced, Symbols, Bits, ... + "precode_mode", duob_mode, 'showAnalysis', 1, "postFFE", [], "eth_style_symbol_mapping", mapping_style); + + mlse_results.metrics.print + if our_signal + fprintf('Our Signal: %.1e \n',mlse_results.metrics.BER); + fprintf('Paper: %.1e \n \n',ber_in_paper); + else + fprintf('Their Signal: %.1e \n',mlse_results.metrics.BER); + fprintf('Paper: %.1e \n \n',ber_in_paper); + end + + %% -------------------- ML-based MLSE (L=2) -------------------- + % ml_mlse_equalizer = ML_MLSE("epochs_tr",100,"epochs_dd",1, ... + % "len_tr",length(Rx_synced),"mu_dd",0.03,"mu_tr",0.03,"order",11,"sps",sps, ... + % "traceback_depth",256,"L",4,"delta",4,"adaptive_mu",0); + % + % [ml_mlse_results] = ml_mlse(ml_mlse_equalizer, M, Rx_synced, Symbols, Bits,"precode_mode",duob_mode,"eth_style_symbol_mapping",mapping_style); + % if our_signal + % fprintf('Our Signal: %.1e \n',ml_mlse_results.metrics.BER); + % fprintf('Paper: %.1e \n \n',ber_in_paper); + % else + % fprintf('Their EQ: %.1e \n',ml_mlse_results.metrics.BER); + % fprintf('Paper: %.1e \n \n',ber_in_paper); + % end +end \ No newline at end of file diff --git a/projects/FSO_transmission/Evaluation Scripts/ml_mlse_complexity_plot.m b/projects/FSO_transmission/Evaluation Scripts/ml_mlse_complexity_plot.m new file mode 100644 index 0000000..d30035b --- /dev/null +++ b/projects/FSO_transmission/Evaluation Scripts/ml_mlse_complexity_plot.m @@ -0,0 +1,11 @@ +N = 1:4; +S = 4; +L = 1:4; +Complexity_N = (N+1).*S.^(1+1); +Complexity_L = (1+1).*S.^(L+1); + +figure; +plot(N,Complexity_N) +hold on +plot(L,Complexity_L) +hold off diff --git a/projects/FSO_transmission/Evaluation Scripts/plotscript_baud_rate_sweep_method.m b/projects/FSO_transmission/Evaluation Scripts/plotscript_baud_rate_sweep_method.m new file mode 100644 index 0000000..ae3a1b5 --- /dev/null +++ b/projects/FSO_transmission/Evaluation Scripts/plotscript_baud_rate_sweep_method.m @@ -0,0 +1,36 @@ +x = 4:8; +for method = 1:4 + filename = fullfile('C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\Silas DSP\imdd_simulation\projects\FSO_transmission\Data', ['BER_PAM_4_baud_rate_method' num2str(method) '.mat']); + data = load(filename); + BER = data.BER_PAM_4; + figure(22113344) + plot(x,BER,'-o','LineWidth',1.75); + hold on +end + +h1 = yline(2e-2, ':k', 'LineWidth',1.5); +h2 = yline(3.8e-3,':b', 'LineWidth',1.5); +h3 = yline(4.85e-3,':g', 'LineWidth',1.5); +h4 = yline(2.2e-4,':r', 'LineWidth',1.5); + +% Legende NUR für Kurven +legend('FFE - 300 Taps','FFE+PF+MLSE - 300 FFE Taps - 4 PF Coefficients','ML-MLSE - 1000 epochs - 100th order - Memory Length = 2','DB - 300 FFE Taps',... + 'Interpreter','latex', ... + 'Location','southeast', 'FontSize', 14) + +% FEC Labels direkt im Plot +text(6.85,2.3e-2,'o-FEC','Color','k','FontSize', 14, 'Interpreter','latex') +text(7,3e-3,'HD-FEC','Color','b','FontSize', 14, 'Interpreter','latex') +text(7,5.6e-3,'KP4+Hamming','Color','g','FontSize', 14,'Interpreter','latex') +text(7,2.5e-4,'KP4','Color','r','FontSize', 14,'Interpreter','latex') + +xlabel('Symbol Rate [GBd]', 'Interpreter','latex') +ylabel('BER', 'Interpreter','latex') +% title('BER for PAM-4', 'Interpreter','latex') + +grid minor +% ylim([1e-4 5e-2]) +xlim([3 9]) +set(gca,'YScale','log') +hold off +% beautifyBERplot \ No newline at end of file diff --git a/projects/FSO_transmission/FSO_timing_recovery_minimal_example.m b/projects/FSO_transmission/FSO_timing_recovery_minimal_example.m new file mode 100644 index 0000000..36884ff --- /dev/null +++ b/projects/FSO_transmission/FSO_timing_recovery_minimal_example.m @@ -0,0 +1,251 @@ +%% +clear all; +close all; + +%% Choose a fitting base leading to the FSO Data +base = "C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\FSO Equalizer\FSO_FP_QCL_60umUTC\"; +mode = 0; %0 oder 1 +M = 4; + +all_files = dir(fullfile(base, "**/*.mat")); + +% data_tr_mf is the already recovered and filtered data +if M == 2 + tx_data_path = fullfile(base, "14G_PAM2\tx_info\tx_info_PAM2_14Gbd0.75RRC.mat"); + filename = fullfile(base, "14G_PAM2\M=2_Rs=1.4e10_Fs=8e10_I=265mA_RoP=46.3mW_L=31m_PS=RRC_rolloff=0.75_Mode=Rise.mat"); + data_tr_mf = load("C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\FSO Equalizer\FSO_FP_QCL_60umUTC\Already Recovered and Filtered\AfterSync_M=2_Rs=1.4e10_Fs=8e10_I=265mA_RoP=46.3mW_L=31m_PS=RRC_rolloff=0.75_Mode=Rise.mat"); +elseif M == 4 + tx_data_path = fullfile(base, "6G_PAM4\tx_info\tx_info_PAM4_6Gbd0.6RRC.mat"); + filename = fullfile(base, "6G_PAM4\M=4_Rs=6e9_Fs=8e10_I=255mA_RoP=42.3mW_L=31m_PS=RRC_rolloff=0.6_Mode=Rise.mat"); + data_tr_mf = load("C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\FSO Equalizer\FSO_FP_QCL_60umUTC\Already Recovered and Filtered\AfterSync_M=4_Rs=6e9_Fs=8e10_I=255mA_RoP=42.3mW_L=31m_PS=RRC_rolloff=0.6_Mode=Rise.mat"); +end + +if mode == 1 + [f, p] = uigetfile(fullfile(base, "**/*.mat")); + if f~=0 + filename = fullfile(p,f); + end +end + +tx_data = load(tx_data_path); +datas = load(filename); + +%% +str = filename; +M_ = str2double(regexp(str, 'M=([^_]+)', 'tokens', 'once')); +assert(M==M_); +fsym = str2double(regexp(str, 'Rs=([^_]+)', 'tokens', 'once')); +fs = str2double(regexp(str, 'Fs=([^_]+)', 'tokens', 'once')); +I = sscanf(char(regexp(str, 'I=([^_]+)', 'tokens', 'once')), '%f'); +rop = sscanf(char(regexp(str, 'RoP=([^_]+)', 'tokens', 'once')), '%f'); +L = sscanf(char(regexp(str, 'L=([^_]+)', 'tokens', 'once')), '%f'); +pulseshape = string( regexp(str, 'PS=([^_]+)', 'tokens', 'once')); +rolloff = str2double(regexp(str, 'rolloff=([^_]+)', 'tokens', 'once')); +mode = string( regexp(str, 'Mode=([^\.]+)', 'tokens', 'once')); + +%% +% Tx data + +Bits = Informationsignal(tx_data.tx_data,"fs",fsym); +Symbols = Informationsignal(real(tx_data.tx_PAM_sym),"fs",fsym); + +mapping_style = M==4; % Pam2 is like move-it; PAM-4 is different, same mapping like ETH peopled used in Zurich... hence the "eth_style" argument here and there +PM = PAMmapper(M,0,"eth_style",mapping_style); % one should rename "eth style" as this is simply a different mapping scheme + +Symbols_ = PM.map(Bits) .* PM.scaling; +assert(isequal(Symbols.signal,Symbols_.signal)); + +Bits_ = PM.demap(Symbols); +[bits,errors,ber,errorIndice] = calc_ber(Bits_.signal,Bits.signal); +assert(ber == 0); + +%% For comparison, apply pulsef on Tx Symbols +Pform = Pulseformer("fsym",fsym,"fdac",fs,"pulse","rrc","pulselength",16,"alpha",rolloff); +Digi_sig_compare = Pform.process(Symbols); + +MF = Pulseformer("fsym",fsym,"fdac",2*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff); +Rx_sig_compare = MF.process(Digi_sig_compare); + +%% + +% Rx Data +traceData = datas.tr.lastData(2).trace.ch3; + +%FYI: Voltage=(RawData−YReference)×YIncrement+YOrigin +scoperead_volts = (traceData.RawData - traceData.YReference) * traceData.YIncrement + traceData.YOrigin; +demystified = isequal(traceData.YData,scoperead_volts); +assert(demystified); + +Scope_sig = Electricalsignal(traceData.YData,"fs",fs); + +Scope_sig.plot("displayname",'raw','fignum',100); +Scope_sig.spectrum("displayname",'raw','fignum',101) + +%Calculate Transfer Function +Tx_spectrum = Digi_sig_compare; +Rx_spectrum = Scope_sig; + +[~,Rx_synced_spectrum,inverted_spectrum,sequenceFound_spectrum,sequenceStarts_spectrum] = Rx_spectrum.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1); +Rx_spectrum = Rx_synced_spectrum{1}; + +Tx_spectrum = Tx_spectrum.resample("fs_out",fsym); +Tx_spectrum.signal = fft(Tx_spectrum.signal); + +Rx_spectrum = Rx_spectrum.resample("fs_out",fsym); +Rx_spectrum.signal = fft(Rx_spectrum.signal); + +H_transfer = Rx_spectrum.signal./Tx_spectrum.signal; +H_inv = 1./H_transfer; + +%Number of Samples/Symbol after Matched Filter +Kov = 14; + +Scope_sig = Scope_sig.resample('fs_in',fs,'fs_out',Kov*fsym); + +% 1) matched filter +% pulse is symmetric, hence we can use pulsef firectly as matched filter. +% It feels off (bit I think correct) that the fsym is now the output freq.!! +% -> output 2 sps to omit timing recovery!? +Matched_Filter = Pulseformer("fsym",fsym,"fdac",Kov*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff,"matched",1); +Rx_matched = Matched_Filter.process(Scope_sig); +Rx_matched.spectrum("displayname",'Signal after matched filter','fignum',1); + +% Loading and synchronizing their matched filtered and timing recovered data +data_tr_mf = Electricalsignal(data_tr_mf.Results, "fs", fsym); +[~,Rx_synced_cell_tr_mf,inverted_tr_mf,sequenceFound_tr_mf,sequenceStarts_tr_mf] = data_tr_mf.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1); +Rx_tr_mf = Rx_synced_cell_tr_mf{1}; + +% Timing sync -> at this point we still have no symbol timing recovery, we +% try to do this with 2sps EQ! +[~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_matched.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1); +Rx_matched_1 = Rx_synced_cell{1}; + +% Timing recovery +Time_Rec = 1; +Timing_Mode = 3; +if Time_Rec + if Timing_Mode == 1 % Zero-Crossing/Gardner/Early-Late/Müller-Mueller Timing Recovery + [Rx_Time_Rec, Timing_Error] = Timing_Recovery("modulation", 'PAM/PSK/QAM', "timing_error_detector",'Gardner (non-data-aided)','sps',Kov,'damping_factor',1,'normalized_loop_bandwidth',1e-4,'detector_gain',2.7).process(Rx_matched_1); + elseif Timing_Mode == 2 % Godard Timing Recovery (resampling is required as the Godard timing recovery does not change the number of samples per symbol) + [Rx_Time_Rec, Timing_Error_MG] = Godard_Timing_Recovery('mode',3,'num_blocks',1,'fft_length',length(Rx_matched_1),'sps',Kov,'rolloff',0.6,'mu',-0.2,'Ki',1e-4).process(Rx_matched_1); + Rx_Time_Rec = Rx_Time_Rec.resample('fs_in',Kov*fsym,'fs_out',fsym); + elseif Timing_Mode == 3 % Maximum Variance Timing Recovery + Rx_Time_Rec = MaxVar_Timing_Recovery('mode',0,'fsym',fsym,'fadc',Kov*fsym,'num_tau',Kov*128,'sps',Kov,'comp_signal',Rx_tr_mf,'comp_mode',0).process(Rx_matched_1); + end +else + Rx_Time_Rec = Rx_matched_1.resample('fs_in',Kov*fsym,'fs_out',fsym); +end +sps = 1; +Rx_Time_Rec.fs = fsym; + +% % Compensate using the inverse transfer function +% Amax = 10; +% H_inv = min(abs(H_inv), Amax) .* exp(1j*angle(H_inv)); +% Rx_Time_Rec.signal = fft(Rx_Time_Rec.signal); +% Rx_Time_Rec.signal = ifft(Rx_Time_Rec.signal .* H_inv); + +% Normalization +Rx_Time_Rec = Rx_Time_Rec.normalize('mode','rms'); +Rx_tr_mf = Rx_tr_mf.normalize('mode','rms'); + +% Compare spectra of our and their time signal +Rx_Time_Rec.spectrum("displayname",'Our signal','normalizeTo0dB',1,'fignum',101111); +% Rx_matched_original.spectrum('normalizeTo0dB',1,"displayname",'Our signal','fignum',101111); +Rx_tr_mf.spectrum("displayname",'Their signal','normalizeTo0dB',1,'fignum',101111); + +% Compare our and their time signal +Rx_Time_Rec.normalize("mode","rms").plot("displayname",'Our signal','fignum',101311); +% Rx_matched_original.normalize("mode","rms").plot("displayname",'Original signal','fignum',101311); +Rx_tr_mf.normalize("mode","rms").plot("displayname",'Their signal','fignum',101311); + +%% +% Use our or their signal +for our_signal = 1 + if our_signal + Rx_synced = Rx_Time_Rec; + else + Rx_synced = Rx_tr_mf; + end + % Rx_synced = Rx_Time_Rec; + % Rx_synced = Rx_synced_cell{1}; + len_tr = 4096*2; + mu_ffe1 = 0.0001; + mu_ffe2 = 0.0008; + mu_ffe3 = 0.001; + mu_dc = 0.004; + mu_ffe = [mu_ffe1 mu_ffe3 mu_ffe3]; + mu_dfe = 0.0004; + duob_mode = db_mode.no_db; + + Rx_synced.plot("displayname",'RX: Matched+Sync+2sps','fignum',103); + Rx_synced.spectrum("displayname",'RX: Matched+Sync+2sps','fignum',104); + + Digi_sig_compare.normalize("mode","rms").spectrum("displayname",'Tx: RC-shaped','fignum',1,'normalizeTo0dB',1); + Rx_sig_compare.normalize("mode","rms").spectrum("displayname",'Tx: RC-shaped + matched filtered ','fignum',1,'normalizeTo0dB',1); + Rx_synced.normalize("mode","rms").spectrum("displayname",'RX: matched filtered + synced','fignum',1,'normalizeTo0dB',1); + + if M == 2 + ber_in_paper = 10^(-2.6); %fig 3a) 4 Gb/s MWIR FSO Transmission using Directly Modulated QCL and an Uncooled UTC-PD at Room-Temperature + elseif M == 4 + ber_in_paper = 10^(-2.5); + end + + %% -------------------- FFE -------------------- + % requires some more digging what is going on :-) + % eq_ffe = EQ("Ne",[150, 0, 0],"Nb",[0,0,0], ... + % "training_length",len_tr,"training_loops",5,"dd_loops",5, ... + % "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... + % "FFEmu",0,"plotfinal",0,"ideal_dfe",1); + % + % % eq_ffe = FFE_DFE('ffe_order',99,'dfe_order',99,'len_tr',len_tr,'epochs_tr',5,'epochs_dd',5,'sps',sps,'decide',0); + % + % ffe_results = ffe(eq_ffe,M,Rx_synced,Symbols,Bits, ... + % "precode_mode",duob_mode,'showAnalysis',0,"postFFE",[], ... + % "eth_style_symbol_mapping",mapping_style); + % + % if our_signal + % fprintf('Our signal: %.1e \n',ffe_results.metrics.BER); + % fprintf('Paper: %.1e \n \n',ber_in_paper); + % else + % fprintf('Their signal: %.1e \n',ffe_results.metrics.BER); + % fprintf('Paper: %.1e \n \n',ber_in_paper); + % end + + %% -------------------- VNLE + MLSE -------------------- + pf_ncoeffs = 4; + eq_v = EQ("Ne",[300, 0, 0],"Nb",[0, 0, 0], ... + "training_length",len_tr,"training_loops",5,"dd_loops",5, ... + "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... + "FFEmu",0,"plotfinal",0,"ideal_dfe",1, ... + 'weighted_DFE',0,'weighted_DFE_d_min',0.5,'weighted_DFE_mode','R2','weighted_DFE_I_mode',[5,0.5,0.6]); + % eq_v = FFE_DFE('ffe_order',300,'dfe_order',5,'len_tr',len_tr,'epochs_tr',5,'epochs_dd',5, ... + % 'ffe_mu_dd',mu_ffe,'ffe_mu_tr',0,'dfe_mu_dd',mu_dfe,'dfe_mu_tr',0.005,'sps',sps,'decide',0); + pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); + mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0,"eth_style",mapping_style).levels); + + [vnle_results, mlse_results] = vnle_postfilter_mlse(eq_v, pf_, mlse_, M, Rx_synced, Symbols, Bits, ... + "precode_mode", duob_mode, 'showAnalysis', 1, "postFFE", [], "eth_style_symbol_mapping", mapping_style); + + mlse_results.metrics.print + if our_signal + fprintf('Our Signal: %.1e \n',mlse_results.metrics.BER); + fprintf('Paper: %.1e \n \n',ber_in_paper); + else + fprintf('Their Signal: %.1e \n',mlse_results.metrics.BER); + fprintf('Paper: %.1e \n \n',ber_in_paper); + end + + %% -------------------- ML-based MLSE (L=2) -------------------- + % ml_mlse_equalizer = ML_MLSE("epochs_tr",100,"epochs_dd",1, ... + % "len_tr",length(Rx_synced),"mu_dd",0.03,"mu_tr",0.03,"order",11,"sps",sps, ... + % "traceback_depth",256,"L",4,"delta",4,"adaptive_mu",0); + % + % [ml_mlse_results] = ml_mlse(ml_mlse_equalizer, M, Rx_synced, Symbols, Bits,"precode_mode",duob_mode,"eth_style_symbol_mapping",mapping_style); + % if our_signal + % fprintf('Our Signal: %.1e \n',ml_mlse_results.metrics.BER); + % fprintf('Paper: %.1e \n \n',ber_in_paper); + % else + % fprintf('Their EQ: %.1e \n',ml_mlse_results.metrics.BER); + % fprintf('Paper: %.1e \n \n',ber_in_paper); + % end +end \ No newline at end of file