diff --git a/Classes/00_signals/Signal.m b/Classes/00_signals/Signal.m index 5b7a710..f56c09d 100644 --- a/Classes/00_signals/Signal.m +++ b/Classes/00_signals/Signal.m @@ -927,7 +927,7 @@ classdef Signal options.displayname = ""; end - mode = 1; + mode = 2; histpoints = 2048; %% verticale resolution histpoints = floor(histpoints/2)*2+1; %% to have the eye digram centered around one point make the vertical resolution uneven diff --git a/Classes/01_transmit/PAMmapper.m b/Classes/01_transmit/PAMmapper.m index a44d8cd..6c1276d 100644 --- a/Classes/01_transmit/PAMmapper.m +++ b/Classes/01_transmit/PAMmapper.m @@ -462,10 +462,19 @@ classdef PAMmapper end - function [Signal_out] = quantize(obj,Signal_in) + function [Signal_out] = quantize(obj,Signal_in,options) + arguments + obj + Signal_in Signal + options.custom_const = []; + end - constellation = obj.get_levels(); - constellation = constellation ./ obj.scaling; + if isempty(options.custom_const) + constellation = obj.get_levels(); + constellation = constellation ./ obj.scaling; + else + constellation = options.custom_const; + end issignalclass = 0; if isa(Signal_in,'Signal') diff --git a/Classes/03_electrical/CTLE.m b/Classes/03_electrical/CTLE.m new file mode 100644 index 0000000..84647fb --- /dev/null +++ b/Classes/03_electrical/CTLE.m @@ -0,0 +1,82 @@ +classdef CTLE < handle + + properties(Access=public) + Aac_dB + Adc_dB + f_p1 + f_p2 + plot + end + + methods(Access=public) + function obj = CTLE(options) + arguments + options.Aac_dB = 0; + options.Adc_dB = -6; + options.f_p1 = 1.5e9; + options.f_p2 = 5e9; + options.plot = 0; + end + + fn = fieldnames(options); + for n = 1:numel(fn) + obj.(fn{n}) = options.(fn{n}); + end + + end + + function data_out = process(obj, data_in) + + x = data_in.signal(:); + Fs = data_in.fs; + N = length(x); + + % CTLE Transfer Function Parameters + A_ac = 10^(obj.Aac_dB/20); + A_dc = 10^(obj.Adc_dB/20); + w1 = 2*pi*obj.f_p1; + w2 = 2*pi*obj.f_p2; + + % Frequency Domain Conversion + X = fft(x); + + % Generating Frequency Axis In The Range Of [-Fs/2,Fs/2] + f_fft = (0:N-1).' * (Fs/N); + f_signed = f_fft; + idxNeg = f_signed > Fs/2; + f_signed(idxNeg) = f_signed(idxNeg) - Fs; + f_pos = abs(f_signed); + w_pos = 2*pi*f_pos; + s_pos = 1j*w_pos; + + % Calculate CTLE Transfer Function + H_pos = A_ac*w2 .* (s_pos + (A_dc/A_ac)*w1) ./ ((s_pos + w1).*(s_pos + w2)); + + % Enforce H(-f) = conj(H(f)) For Negative Bins + H = H_pos; + H(idxNeg) = conj(H_pos(idxNeg)); + + % Apply CTLE + Y = H .* X; + + % Calculate Time Domain Signal + y = ifft(Y, 'symmetric'); + data_out = data_in; + data_out.signal = y; + data_out.fs = Fs; + + % Plot + if obj.plot + % plot only positive frequencies up to Fs/2 + k = 1:floor(N/2)+1; + fplot = f_fft(k); + Hplot = H(k); + + figure; + semilogx(fplot, 20*log10(abs(Hplot)+1e-15)); + grid on; xlabel('frequency (Hz)'); ylabel('magnitude (dB)'); + title('CTLE magnitude on FFT grid'); + end + end + end +end \ No newline at end of file diff --git a/Classes/03_electrical/Electrical_Hybrid.m b/Classes/03_electrical/Electrical_Hybrid.m new file mode 100644 index 0000000..c1fa1b1 --- /dev/null +++ b/Classes/03_electrical/Electrical_Hybrid.m @@ -0,0 +1,130 @@ +classdef Electrical_Hybrid < handle + + properties(Access=public) + file_path + plot = 0 + + % If true: perform digital residual-echo cancellation using known TX + cancel_echo = 1 + + % If true: in addition return v_hyb and v_echo_est + return_intermediates = 1 + end + + methods(Access=public) + + function obj = Electrical_Hybrid(options) + arguments + options.file_path = '' + options.plot = 0 + options.cancel_echo = 1 + options.return_intermediates = 1 + end + + fn = fieldnames(options); + for n = 1:numel(fn) + obj.(fn{n}) = options.(fn{n}); + end + end + + function [v_fe_rec, v_hyb, v_echo_est] = process(obj, v_ne_tx, v_fe_tx) + % v_ne_tx : Near-End TX signal object (Port 1) + % v_fe_tx : Far-End TX signal object (Port 2) + % + % Output: + % v_hyb : physical hybrid differential output (Port4 - Port3) + % v_echo_est : estimated residual echo due to local TX only + % v_fe_rec : v_hyb - v_echo_est (if cancel_echo enabled), else v_hyb + + % --- Load S-Parameters (.s4p) --- + net = sparameters(obj.file_path); + f_s = net.Frequencies(:); + S4 = net.Parameters; + + % Extract needed S-parameters for differential output: + % V3 = S31*V1 + S32*V2 + % V4 = S41*V1 + S42*V2 + % Vhyb = V4 - V3 = (S41-S31)*V1 + (S42-S32)*V2 + S31_s = squeeze(S4(3,1,:)); + S41_s = squeeze(S4(4,1,:)); + S32_s = squeeze(S4(3,2,:)); + S42_s = squeeze(S4(4,2,:)); + + % --- Time-domain signals --- + x1 = v_ne_tx.signal(:); + x2 = v_fe_tx.signal(:); + + if length(x2) ~= length(x1) + error('Near-end and far-end signals must have the same length.'); + end + + N = length(x1); + Fs = v_ne_tx.fs; + + % --- Use zero padding to avoid circular convolution artifacts --- + Nfft = 2^nextpow2(2*N); % robust choice + + % --- FFT --- + V1 = fft(x1, Nfft); + V2 = fft(x2, Nfft); + + % --- Frequency axis for interpolation (signed then abs) --- + f_fft = (0:Nfft-1).' * (Fs/Nfft); + f_signed = f_fft; + idxNeg = f_signed > Fs/2; + f_signed(idxNeg) = f_signed(idxNeg) - Fs; % (-Fs/2, Fs/2] + f_pos = abs(f_signed); + + % --- Interpolate S-parameters onto f_pos --- + S31 = interp1(f_s, S31_s, f_pos, 'linear', 'extrap'); + S41 = interp1(f_s, S41_s, f_pos, 'linear', 'extrap'); + S32 = interp1(f_s, S32_s, f_pos, 'linear', 'extrap'); + S42 = interp1(f_s, S42_s, f_pos, 'linear', 'extrap'); + + % Hermitian symmetry for real time-domain response: + % For negative frequencies enforce conj symmetry. + S31(idxNeg) = conj(S31(idxNeg)); + S41(idxNeg) = conj(S41(idxNeg)); + S32(idxNeg) = conj(S32(idxNeg)); + S42(idxNeg) = conj(S42(idxNeg)); + + % --- Physical hybrid differential output --- + % Vhyb = (S41-S31)*V1 + (S42-S32)*V2 + He = (S41 - S31); % residual echo transfer from local TX + Hr = (S42 - S32); % transfer from far-end TX to output + + V_hyb = He .* V1 + Hr .* V2; + + % --- Residual echo estimate (digital canceller model) --- + V_echo = He .* V1; + + % --- Back to time-domain (take first N samples after padding) --- + v_hyb_full = ifft(V_hyb, 'symmetric'); + v_echo_full = ifft(V_echo, 'symmetric'); + + v_hyb = v_hyb_full(1:N); + v_echo_est = v_echo_full(1:N); + + % --- Optional cancellation --- + if obj.cancel_echo + v_fe_rec = v_hyb - v_echo_est; + else + v_fe_rec = v_hyb; + end + + if ~obj.return_intermediates + v_hyb = []; + v_echo_est = []; + end + + % --- Bring output in the correct form --- + v_fe_rec = Informationsignal(v_fe_rec,"fs",v_fe_tx.fs); + v_hyb = Informationsignal(v_hyb,"fs",v_fe_tx.fs); + v_echo_est = Informationsignal(v_echo_est,"fs",v_fe_tx.fs); + + if obj.plot + rfplot(net) + end + end + end +end \ No newline at end of file diff --git a/Classes/03_electrical/Electrical_Trace.m b/Classes/03_electrical/Electrical_Trace.m new file mode 100644 index 0000000..8231929 --- /dev/null +++ b/Classes/03_electrical/Electrical_Trace.m @@ -0,0 +1,31 @@ +classdef Electrical_Trace < handle + + properties(Access=public) + file_path + end + + methods(Access=public) + function obj = Electrical_Trace(options) + arguments(Input) + + options.file_path = 'C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\COM Test\Mellitzz\TA_6002_6003_FX_B6_C6_B7_C7_Terminated.s4p' + + end + + fn = fieldnames(options); + for n = 1:numel(fn) + obj.(fn{n}) = options.(fn{n}); + end + + + + end + + function [data_out,timing_error] = process(obj, data_in) + + S = + + end + end +end + diff --git a/Classes/03_electrical/Electrical_Trace_BiDi.m b/Classes/03_electrical/Electrical_Trace_BiDi.m new file mode 100644 index 0000000..468c73f --- /dev/null +++ b/Classes/03_electrical/Electrical_Trace_BiDi.m @@ -0,0 +1,114 @@ +classdef Electrical_Trace_BiDi < handle + + properties(Access=public) + file_path + fsym + rolloff + K_over + plot + test + S_test + end + + methods(Access=public) + function obj = Electrical_Trace_BiDi(options) + arguments(Input) + + options.file_path = 'C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\COM Test\Mellitzz\TA_6002_6003_FX_B6_C6_B7_C7_Terminated.s4p' + options.fsym = 0; + options.rolloff = 0; + options.K_over = 1; + options.plot = 0; + options.test = 0; + options.S_test = [0,1;1,0]; + end + + fn = fieldnames(options); + for n = 1:numel(fn) + obj.(fn{n}) = options.(fn{n}); + end + + + + end + + function [b_1,b_2] = process(obj, a_1, a_2) + + % Rx Signal Calculation Using S-Paramters For a 2-Port Network + % [B_1(f); B_2(f)] = [S_11(f), S_12(f); S_21(f), S_22(f)] * [A_1(f); A_2(f)] + + % Initialize Rx Time Domain Signals + b_1 = a_2; + b_2 = a_1; + + % Loading and Extract S-Paramters + net = sparameters(obj.file_path); + f_s = net.Frequencies(:); + + if ~obj.test + S2 = net.Parameters(1:2, 1:2, :); + else + S2 = repmat(obj.S_test,1,1,size(f_s,1)); + end + + S11_s = squeeze(S2(1,1,:)); + S12_s = squeeze(S2(1,2,:)); + S21_s = squeeze(S2(2,1,:)); + S22_s = squeeze(S2(2,2,:)); + + % Setup Time Domain Signals + x1 = a_1.signal(:); + x2 = a_2.signal(:); + N = length(x1); + assert(length(x2)==N, 'a_1 and a_2 must have same length'); + + % Extract Time Domain Signal Frequencies + if obj.fsym == 0 && obj.rolloff == 0 && obj.K_over == 1 + Fs = a_1.fs; + else + Fs = (1+obj.rolloff)*obj.fsym; + end + + % Calculate Frequency Domain Signals + A1 = fft(x1); + A2 = fft(x2); + + % Calculate Frequency Axis [-Fs/2,...,Fs/2] + f_fft = (0:N-1).' * (Fs/N); + f_signed = f_fft; + idxNeg = f_signed > Fs/2; + f_signed(idxNeg) = f_signed(idxNeg) - Fs; % Now in (-Fs/2, Fs/2] + + % Absolute Value Used for Interpolation + f_pos = abs(f_signed); + + % Interpolate S-Parameters onto f_pos + S11_p = interp1(f_s, S11_s, f_pos, 'linear', 'extrap'); + S12_p = interp1(f_s, S12_s, f_pos, 'linear', 'extrap'); + S21_p = interp1(f_s, S21_s, f_pos, 'linear', 'extrap'); + S22_p = interp1(f_s, S22_s, f_pos, 'linear', 'extrap'); + + % Apply Hermitian Symmetry: S(-f) = conj(S(f)) + S11 = S11_p; S12 = S12_p; S21 = S21_p; S22 = S22_p; + S11(idxNeg) = conj(S11_p(idxNeg)); + S12(idxNeg) = conj(S12_p(idxNeg)); + S21(idxNeg) = conj(S21_p(idxNeg)); + S22(idxNeg) = conj(S22_p(idxNeg)); + + % Compute Rx Frequency Domain Signals + B1 = S11 .* A1 + S12 .* A2; + B2 = S21 .* A1 + S22 .* A2; + + % Calculate Rx Time Domain Signals + b_1.signal = ifft(B1, 'symmetric'); + b_2.signal = ifft(B2, 'symmetric'); + + if obj.plot + rfplot(net) + end + + + end + end +end + diff --git a/Classes/04_DSP/Equalizer/EQ.m b/Classes/04_DSP/Equalizer/EQ.m index 1be73b2..bd19004 100644 --- a/Classes/04_DSP/Equalizer/EQ.m +++ b/Classes/04_DSP/Equalizer/EQ.m @@ -480,7 +480,7 @@ classdef EQ < handle output_vec_weighted(m) = xbar; fb_sym = xbar; - output_vec(m) = xbar; + % output_vec(m) = xbar; end if obj.Nb(1) > 0 diff --git a/Classes/04_DSP/Timing Recovery/MaxVar_Timing_Recovery.m b/Classes/04_DSP/Timing Recovery/MaxVar_Timing_Recovery.m index d6517d5..f7eb502 100644 --- a/Classes/04_DSP/Timing Recovery/MaxVar_Timing_Recovery.m +++ b/Classes/04_DSP/Timing Recovery/MaxVar_Timing_Recovery.m @@ -103,6 +103,7 @@ classdef MaxVar_Timing_Recovery < handle % y = [y 0]; % end data_out.signal = y.'; + data_out.fs = obj.fsym; end diff --git a/Datatypes/db_decoder.m b/Datatypes/db_decoder.m new file mode 100644 index 0000000..78f09d5 --- /dev/null +++ b/Datatypes/db_decoder.m @@ -0,0 +1,8 @@ +classdef db_decoder < int32 + + enumeration + sequencedetection (0) % use MLSE for decoding + memoryless (1) % use modulo + end + +end \ No newline at end of file diff --git a/Documentations/Demystifying Duobinary.pdf b/Documentations/Demystifying Duobinary.pdf new file mode 100644 index 0000000..8423506 Binary files /dev/null and b/Documentations/Demystifying Duobinary.pdf differ diff --git a/Functions/EQ_structures/duobinary_target.m b/Functions/EQ_structures/duobinary_target.m index 36150dd..c24ffca 100644 --- a/Functions/EQ_structures/duobinary_target.m +++ b/Functions/EQ_structures/duobinary_target.m @@ -1,153 +1,187 @@ function [db_results] = duobinary_target(eq_, mlse_,M, rx_signal, tx_symbols, tx_bits, options) -arguments - eq_ - mlse_ - M - rx_signal - tx_symbols - tx_bits - options.precode_mode db_mode - options.showAnalysis = 0; - options.eth_style_symbol_mapping = 0; - options.postFFE = []; -end + arguments + eq_ + mlse_ + M + rx_signal + tx_symbols + tx_bits + options.precode_mode db_mode + options.showAnalysis = 0; + options.eth_style_symbol_mapping = 0; + options.postFFE = []; + options.decoding_mode db_decoder = db_decoder.sequencedetection; + end + + + %Duobinary Targeting + db_ref_sequence = Duobinary().encode(tx_symbols); + db_ref_constellation = unique(db_ref_sequence.signal); + [eq_signal, eq_noise] = eq_.process(rx_signal,db_ref_sequence); + + if ~isempty(options.postFFE) + [eq_signal,eq_noise] = options.postFFE.process(eq_signal,db_ref_sequence); + end + % + + switch options.decoding_mode + case db_decoder.sequencedetection %MLSE + mlse_.DIR = [1,1]; + if isa(mlse_,'MLSE_viterbi') + mlse_sig_sd = mlse_.process(eq_signal); + else + [mlse_sig_sd,LLR,GMI_MLSE] = mlse_.process(eq_signal,tx_symbols); + end + pam_sig_hd = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).quantize(mlse_sig_sd); + case db_decoder.memoryless %DB Target FFE + % Hard decision on FFE output + eq_signal_hd = PAMmapper(M, 0).quantize(eq_signal,'custom_const',db_ref_constellation.'); + eq_signal_hd = Duobinary().decode(eq_signal_hd); + pam_sig_hd = eq_signal_hd; + % tx_symbols = Duobinary().encode(tx_symbols); + % tx_symbols = Duobinary().decode(tx_symbols.*db_const_scale_factor); + end + + % precoding to mitigate error propagation, most prominently used in + % combination with duobinary signaling to avoid catastrophic error + % behavior (see J.W.M. Bergmans, Digital Baseband Transmission and Recording -> partial response signaling) + switch options.precode_mode + + case db_mode.no_db + % TX Data is not precoded: + + % A) Emulate diff precoding + switch options.decoding_mode + case db_decoder.sequencedetection %MLSE + mlse_sig_hd_precoded = Duobinary().encode(pam_sig_hd); + mlse_sig_hd_precoded = Duobinary().decode(mlse_sig_hd_precoded); + case db_decoder.memoryless %DB Target FFE + % mlse_sig_hd_precoded = Duobinary().encode(pam_sig_hd); + % mlse_sig_hd_precoded = Duobinary().decode(mlse_sig_hd_precoded); + mlse_sig_hd_precoded = pam_sig_hd; + end -%Duobinary Targeting -db_ref_sequence = Duobinary().encode(tx_symbols); -db_ref_constellation = unique(db_ref_sequence.signal); -[eq_signal, eq_noise] = eq_.process(rx_signal,db_ref_sequence); + tx_symbols_precoded = Duobinary().encode(tx_symbols); + tx_symbols_precoded = Duobinary().decode(tx_symbols_precoded); -if ~isempty(options.postFFE) - [eq_signal,eq_noise] = options.postFFE.process(eq_signal,db_ref_sequence); -end + tx_bits_precoded = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(tx_symbols_precoded); + rx_bits_mlse = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(mlse_sig_hd_precoded); + + [~,errors_db_diff_precoded,ber_db_diff_precoded,~] = calc_ber(rx_bits_mlse.signal,tx_bits_precoded.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + + %B) Just determine BER + rx_bits_mlse = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(pam_sig_hd); + tx_bits = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(tx_symbols); + [bits_mlse,errors_mlse,ber_db,~] = calc_ber(rx_bits_mlse.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + + case db_mode.db_precoded + + % Daten SIND TATSÄCHLICH precoded auf TX Seite: + + % A) Decode at Rx if no DB targeting was applied (we are in VNLE or MLSE EQ structure here! + switch options.decoding_mode + case db_decoder.sequencedetection %MLSE + mlse_sig_hd_decoded = Duobinary().encode(pam_sig_hd,"M",M); + mlse_sig_hd_decoded = Duobinary().decode(mlse_sig_hd_decoded,"M",M); + case db_decoder.memoryless %DB Target FFE + mlse_sig_hd_decoded = pam_sig_hd; + end -mlse_.DIR = [1,1]; -% + tx_symbols_precoded = Duobinary().encode(tx_symbols,"M",M); + tx_symbols_precoded = Duobinary().decode(tx_symbols_precoded,"M",M); -if isa(mlse_,'MLSE_viterbi') - mlse_sig_sd = mlse_.process(eq_signal); -else - [mlse_sig_sd,LLR,GMI_MLSE] = mlse_.process(eq_signal,tx_symbols); -end + tx_bits_precoded = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(tx_symbols_precoded); + rx_bits_mlse_decoded = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(mlse_sig_hd_decoded); -mlse_sig_hd = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).quantize(mlse_sig_sd); - -% precoding to mitigate error propagation, most prominently used in -% combination with duobinary signaling to avoid catastrophic error -% behavior (see J.W.M. Bergmans, Digital Baseband Transmission and Recording -> partial response signaling) -switch options.precode_mode - - case db_mode.no_db - % TX Data is not precoded: - - % A) Emulate diff precoding - mlse_sig_hd_precoded = Duobinary().encode(mlse_sig_hd,"M",M); - mlse_sig_hd_precoded = Duobinary().decode(mlse_sig_hd_precoded,"M",M); - - tx_symbols_precoded = Duobinary().encode(tx_symbols); - tx_symbols_precoded = Duobinary().decode(tx_symbols_precoded); - - tx_bits_precoded = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(tx_symbols_precoded); - rx_bits_mlse = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(mlse_sig_hd_precoded); - - [~,errors_db_diff_precoded,ber_db_diff_precoded,~] = calc_ber(rx_bits_mlse.signal,tx_bits_precoded.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - - %B) Just determine BER - rx_bits_mlse = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(mlse_sig_hd); - [bits_mlse,errors_mlse,ber_db,~] = calc_ber(rx_bits_mlse.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - - case db_mode.db_precoded - - % Daten SIND TATSÄCHLICH precoded auf TX Seite: - - % A) Decode at Rx if no DB targeting was applied (we are in VNLE or MLSE EQ structure here! - mlse_sig_hd_decoded = Duobinary().encode(mlse_sig_hd,"M",M); - mlse_sig_hd_decoded = Duobinary().decode(mlse_sig_hd_decoded,"M",M); - rx_bits_mlse_decoded = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(mlse_sig_hd_decoded); - [~,errors_db_diff_precoded,ber_db_diff_precoded,a] = calc_ber(rx_bits_mlse_decoded.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - burst_db_precoded = count_error_bursts(a, 40); - % B) Omit the Coding by comparing with demapped TX symbol sequence - - tx_bits = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(tx_symbols); - rx_bits_mlse = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(mlse_sig_hd); - [bits_db,errors_db,ber_db,a] = calc_ber(rx_bits_mlse.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - burst_db = count_error_bursts(a, 40); - - cols = linspecer(8); - figure();hold on; - stem(1:40,burst_db,'LineWidth',1,'Color',cols(4,:),'Marker','_','DisplayName','w/o diff. precoder'); - stem(1:40,burst_db_precoded,'LineWidth',1,'Color',cols(3,:),'Marker','.','LineStyle','-','DisplayName','w diff. precoder'); - xlabel('Bit Error Burst Length') - ylabel('Occurence') - set(gca, 'yscale', 'log'); -end - -% M = numel(unique(tx_symbols.signal)); -rx_bits = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(mlse_sig_hd); - -[bits_db,errors_db,ber_db,errorIndice_db] = calc_ber(rx_bits.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - - - - - -alpha = arburg(eq_noise.signal,1);%pf_.coefficients(2); -alpha = alpha(2); -if isa(mlse_,'MLSE_viterbi') - gmi_mlse = NaN; - air_mlse = NaN; -else - gmi_mlse = GMI_MLSE; - air_mlse = tx_symbols.fs .* floor(log2(double(M))*10)/10 .* gmi_mlse ./ log2(double(M)); -end - -db_results = struct(); -db_results.metrics = Metricstruct; -db_results.metrics.result_id = NaN; -db_results.metrics.run_id = NaN; -db_results.metrics.eqParam_id = NaN; -db_results.metrics.date_of_processing = datetime('now'); -db_results.metrics.BER = ber_db; -db_results.metrics.numBits = bits_db; -db_results.metrics.numBitErr = errors_db; -db_results.metrics.BER_precoded = ber_db_diff_precoded; -db_results.metrics.numBitErr_precoded = errors_db_diff_precoded; -db_results.metrics.GMI = gmi_mlse; -db_results.metrics.AIR = air_mlse; -db_results.metrics.MLSE_dir = mlse_.DIR; -db_results.metrics.Alpha = alpha; - -% Create DB results structure -db_results.config = Equalizerstruct(); -eq_.e = []; -eq_.e2 = []; -eq_.e3 = []; -db_results.config.eq = jsonencode(eq_); -% mlse_.DIR = []; -db_results.config.mlse = jsonencode(mlse_); -db_results.config.equalizer_structure = int32(equalizer_structure.vnle_db_mlse); -db_results.config.comment = 'function: Duobinary tgt. (VNLE -> MLSE)'; - -if options.showAnalysis - - eq_signal.eye(eq_signal.fs,M,"fignum",249); - - - eq_noise = eq_noise - mean(eq_noise.signal); - - rx_signal.spectrum("normalizeTo0dB",1,"fignum",250,"displayname","Rx Spectrum"); - - Duobinary().encode(tx_symbols).spectrum("normalizeTo0dB",1,"fignum",10,"displayname","DB encoded reference"); - - showEQNoisePSD(eq_noise,"fignum",250,"displayname",'Duobinary Target Noise after Equalization'); - - fprintf('DB tgt BER: %.2e \n',ber_db); - - figure(341); clf; - showLevelHistogram(eq_signal, db_ref_sequence, "fignum", 341); - -end + [~,errors_db_diff_precoded,ber_db_diff_precoded,a] = calc_ber(rx_bits_mlse_decoded.signal,tx_bits_precoded.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + burst_db_precoded = count_error_bursts(a, 40); + % B) Omit the Coding by comparing with demapped TX symbol sequence + + tx_bits = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(tx_symbols); + rx_bits_mlse = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(pam_sig_hd); + [bits_db,errors_db,ber_db,a] = calc_ber(rx_bits_mlse.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + burst_db = count_error_bursts(a, 40); + + cols = linspecer(8); + figure();hold on; + stem(1:40,burst_db,'LineWidth',1,'Color',cols(4,:),'Marker','_','DisplayName','w/o diff. precoder'); + stem(1:40,burst_db_precoded,'LineWidth',1,'Color',cols(3,:),'Marker','.','LineStyle','-','DisplayName','w diff. precoder'); + xlabel('Bit Error Burst Length') + ylabel('Occurence') + set(gca, 'yscale', 'log'); + end + + % M = numel(unique(tx_symbols.signal)); + rx_bits = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(pam_sig_hd); + + [bits_db,errors_db,ber_db,errorIndice_db] = calc_ber(rx_bits.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + + alpha = arburg(eq_noise.signal,1);%pf_.coefficients(2); + alpha = alpha(2); + + switch options.decoding_mode + case db_decoder.sequencedetection %MLSE + if isa(mlse_,'MLSE_viterbi') + gmi_mlse = NaN; + air_mlse = NaN; + else + gmi_mlse = GMI_MLSE; + air_mlse = tx_symbols.fs .* floor(log2(double(M))*10)/10 .* gmi_mlse ./ log2(double(M)); + end + case db_decoder.memoryless %DB Target FFE + % [gmi] = calc_air(eq_signal_sd, tx_symbols, "skip_front", 10000, "skip_end", 10000); + [gmi] = calc_ngmi(eq_signal,tx_symbols); + gmi_mlse = max(gmi,0); + air_mlse = tx_symbols.fs .* floor(log2(double(M))*10)/10 .* gmi ./ log2(double(M)); + end + + db_results = struct(); + db_results.metrics = Metricstruct; + db_results.metrics.result_id = NaN; + db_results.metrics.run_id = NaN; + db_results.metrics.eqParam_id = NaN; + db_results.metrics.date_of_processing = datetime('now'); + db_results.metrics.BER = ber_db; + db_results.metrics.numBits = bits_db; + db_results.metrics.numBitErr = errors_db; + db_results.metrics.BER_precoded = ber_db_diff_precoded; + db_results.metrics.numBitErr_precoded = errors_db_diff_precoded; + db_results.metrics.GMI = gmi_mlse; + db_results.metrics.AIR = air_mlse; + db_results.metrics.MLSE_dir = mlse_.DIR; + db_results.metrics.Alpha = alpha; + + % Create DB results structure + db_results.config = Equalizerstruct(); + eq_.e = []; + eq_.e2 = []; + eq_.e3 = []; + db_results.config.eq = jsonencode(eq_); + % mlse_.DIR = []; + db_results.config.mlse = jsonencode(mlse_); + db_results.config.equalizer_structure = int32(equalizer_structure.vnle_db_mlse); + db_results.config.comment = 'function: Duobinary tgt. (VNLE -> MLSE)'; + + if options.showAnalysis + + eq_signal.eye(eq_signal.fs,M,"fignum",249); + + + eq_noise = eq_noise - mean(eq_noise.signal); + + rx_signal.spectrum("normalizeTo0dB",1,"fignum",250,"displayname","Rx Spectrum"); + + Duobinary().encode(tx_symbols).spectrum("normalizeTo0dB",1,"fignum",10,"displayname","DB encoded reference"); + + showEQNoisePSD(eq_noise,"fignum",250,"displayname",'Duobinary Target Noise after Equalization'); + + fprintf('DB tgt BER: %.2e \n',ber_db); + + figure(341); clf; + showLevelHistogram(eq_signal, db_ref_sequence, "fignum", 341); + + end end \ No newline at end of file diff --git a/Functions/EQ_structures/ffe_db.m b/Functions/EQ_structures/ffe_db.m new file mode 100644 index 0000000..eaf88d6 --- /dev/null +++ b/Functions/EQ_structures/ffe_db.m @@ -0,0 +1,211 @@ +function [ffe_results] = ffe_db(eq_, M, rx_signal, tx_symbols, tx_bits, options) +% FFE Processes signals through FFE equalizer +% +% Inputs: +% eq_ - Equalizer object +% M - Modulation order +% rx_signal - Received signal +% tx_symbols - Transmitted symbols +% tx_bits - Transmitted bits +% options - Optional parameters +% +% Outputs: +% ffe_results - Results from FFE processing + +arguments + eq_ + M + rx_signal + tx_symbols + tx_bits + options.db_target = 0; + options.precode_mode db_mode + options.showAnalysis = 0; + options.eth_style_symbol_mapping = 0; + options.postFFE = []; + options.database = []; +end + +%% Process signals through equalizer +% FFE or VNLE +if options.db_target + tx_symbols_ref = Duobinary().encode(tx_symbols); + db_ref_constellation = unique(tx_symbols_ref.signal); + [eq_signal_sd, eq_noise] = eq_.process(rx_signal, tx_symbols_ref); +else + [eq_signal_sd, eq_noise] = eq_.process(rx_signal, tx_symbols); +end + +% Apply post-FFE if provided +if ~isempty(options.postFFE) + tic + [eq_signal_sd, eq_noise] = options.postFFE.process(eq_signal_sd, tx_symbols); + toc +end + +try + ch_coefficients = arburg(eq_noise.signal,1); + channel_alpha = ch_coefficients(2); +end + +% Hard decision on FFE output +if options.db_target + eq_signal_hd = PAMmapper(M, 0).quantize(eq_signal_sd,'custom_const',db_ref_constellation.'); +else + eq_signal_hd = PAMmapper(M, 0).quantize(eq_signal_sd); +end + +if options.db_target + eq_signal_hd = Duobinary().decode(eq_signal_hd); + tx_symbols = Duobinary().encode(tx_symbols); + tx_symbols = Duobinary().decode(tx_symbols); +end + +%% Calculate BER based on precoding mode +[bits, errors, ber, error_pos, errors_precoded, ber_precoded] = calculateBER(eq_signal_hd, tx_symbols, tx_bits, options.precode_mode, M, options.eth_style_symbol_mapping); + +%% Calculate performance metrics +[snr, snr_lvl] = calc_snr(tx_symbols.signal, eq_noise.signal); +% [gmi] = calc_air(eq_signal_sd, tx_symbols, "skip_front", 10000, "skip_end", 10000); +[gmi] = calc_ngmi(eq_signal_sd,tx_symbols); +gmi = max(gmi,0); + +air = tx_symbols.fs .* floor(log2(double(M))*10)/10 .* gmi ./ log2(double(M)); +[evm_total, evm_lvl] = calc_evm(eq_signal_sd, tx_symbols); +[std_total, std_lvl] = calc_std(eq_signal_sd, tx_symbols); +[std_rxraw_total, std_rxraw_lvl] = calc_std(rx_signal.resample("fs_out", tx_symbols.fs), tx_symbols); + +%% Display analysis if requested +if options.showAnalysis + displayAnalysis(eq_noise, eq_signal_sd, rx_signal, eq_, tx_symbols, M, options.postFFE); +end + + +%% Prepare output structure +% Determine postFFE order +if ~isempty(options.postFFE) + npostFFE = options.postFFE.order; +else + npostFFE = 0; +end + +% Create FFE results structure +ffe_results = struct(); +try + eq_.e = []; + eq_.e2 = []; + eq_.e3 = []; + eq_.b = []; + eq_.b2 = []; + eq_.b3 = []; +end + +ffe_results.config = Equalizerstruct(); +ffe_results.config.eq = jsonencode(eq_); +ffe_results.config.equalizer_structure = int32(equalizer_structure.ffe); +ffe_results.config.comment = 'function: ffe'; + +ffe_results.metrics = Metricstruct; +ffe_results.metrics.result_id = NaN; +ffe_results.metrics.run_id = NaN; +ffe_results.metrics.eqParam_id = NaN; +ffe_results.metrics.date_of_processing = datetime('now'); +ffe_results.metrics.BER = ber; +ffe_results.metrics.numBits = bits; +ffe_results.metrics.numBitErr = errors; +ffe_results.metrics.BER_precoded = ber_precoded; +ffe_results.metrics.numBitErr_precoded = errors_precoded; +ffe_results.metrics.SNR = snr; +ffe_results.metrics.SNR_level = snr_lvl; +ffe_results.metrics.STD = std_total; +ffe_results.metrics.STD_level = std_lvl; +ffe_results.metrics.STDrx = std_rxraw_total; +ffe_results.metrics.STDrx_level = std_rxraw_lvl; +ffe_results.metrics.GMI = gmi; +ffe_results.metrics.AIR = air; +ffe_results.metrics.EVM = evm_total; +ffe_results.metrics.EVM_level = evm_lvl; +ffe_results.metrics.Alpha = channel_alpha; + + +end + +%% Helper Functions +function [bits, errors, ber, error_pos, errors_precoded, ber_precoded] = calculateBER(eq_signal_hd, tx_symbols, tx_bits, precode_mode, M, eth_style) +% Calculate BER based on precoding mode +mapper = PAMmapper(M, 0, "eth_style", eth_style); + +switch precode_mode + case db_mode.no_db + % TX Data is not precoded + % A) Emulate diff precoding + eq_signal_hd_precoded = Duobinary().encode(eq_signal_hd, "M", M); + eq_signal_hd_precoded = Duobinary().decode(eq_signal_hd_precoded, "M", M); + + tx_symbols_precoded = Duobinary().encode(tx_symbols); + tx_symbols_precoded = Duobinary().decode(tx_symbols_precoded); + + tx_bits_precoded = mapper.demap(tx_symbols_precoded); + + rx_bits = mapper.demap(eq_signal_hd_precoded); + [~, errors_precoded, ber_precoded, ~] = calc_ber(rx_bits.signal, tx_bits_precoded.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1); + + % B) Just determine BER + rx_bits = mapper.demap(eq_signal_hd); + tx_bits = mapper.demap(tx_symbols); + [bits, errors, ber, error_pos] = calc_ber(rx_bits.signal, tx_bits.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1); + + case db_mode.db_precoded + % Data is precoded on TX side + % A) Decode at Rx if no DB targeting was applied + eq_signal_hd_decoded = Duobinary().encode(eq_signal_hd, "M", M); + eq_signal_hd_decoded = Duobinary().decode(eq_signal_hd_decoded, "M", M); + rx_bits_decoded = mapper.demap(eq_signal_hd_decoded); + [~, errors_precoded, ber_precoded, ~] = calc_ber(rx_bits_decoded.signal, tx_bits.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1); + + % B) Omit the Coding by comparing with demapped TX symbol sequence + tx_bits_demapped = mapper.demap(tx_symbols); + rx_bits = mapper.demap(eq_signal_hd); + [bits, errors, ber, error_pos] = calc_ber(rx_bits.signal, tx_bits_demapped.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1); + end +end + +function displayAnalysis(eq_noise, eq_signal_sd, rx_signal, eq_, tx_symbols, M, postFFE) + % Display analysis plots and metrics + + % Initialize figure handles + % Corrected line - added tx_symbols as second positional argument + % showLevelScatter(rx_signal.resample("fs_out", tx_symbols.fs), tx_symbols, "fignum", 100); + + warning off + showLevelScatter(eq_signal_sd, tx_symbols, "fignum", 101); + figure(gcf);hold on; plot(((1:length(eq_noise.signal)) / eq_noise.fs) * 1e6,movmean(eq_noise.signal,2000,1), 'LineWidth',3,'Color','black') + warning on + + showLevelHistogram(eq_signal_sd, tx_symbols, "fignum", 102); + + showEQNoisePSD(eq_noise, "fignum", 103, "displayname", 'Residual Noise after FFE'); + + % Figure 2: Post-FFE coefficients (if available) + if ~isempty(postFFE) + showEQcoefficients('n1', postFFE.e, "displayname", 'Coefficients', 'fignum', 104); + end + + try + figure(339);hold on + showEQfilter(eq_.e, eq_signal_sd.fs.*2,"displayname",'training','fignum',339); + % showEQfilter(eq_.e, eq_signal_sd.fs.*2,"displayname",'dec. directed','fignum',339); + legend on + end + + % try + % figure(240); hold on; plot(pow2db(movmean(eq_.debug_struct.error_tr',100)));ylim([-30,3]);title('error training'); + % + % figure(241); hold on; plot(pow2db(movmean(eq_.debug_struct.update_tr',100)));title('update step training'); + % + % figure(242); hold on; plot(pow2db(movmean(eq_.debug_struct.update',1000)));title('update step dd'); + % end + + eq_signal_sd.eye(eq_signal_sd.fs,M,"displayname",'Eye','fignum',105); + +end \ No newline at end of file diff --git a/projects/ASTRO Electrical Link/Electrical_Link_GPT.m b/projects/ASTRO Electrical Link/Electrical_Link_GPT.m new file mode 100644 index 0000000..a3da40e --- /dev/null +++ b/projects/ASTRO Electrical Link/Electrical_Link_GPT.m @@ -0,0 +1,195 @@ +%========================================================================== +% High-Speed Link Modell (Dispersion + Reflexionen + Crosstalk) +% Eigenständiges MATLAB-Skript mit Testsequenz (PRBS), Kanalmodell im f-Bereich +% und Auswertung (Eye + Spektren). +% +% Idee: +% - Hauptkanal: H11(f) = exp(-gamma(f)*l) +% - Crosstalk: H21(f) = k_xt * j*(f/fc) / (1 + j*(f/fc)) * exp(-gamma_xt(f)*l) +% (Highpass-artig: typisches FEXT-ähnliches Verhalten; simpel aber nützlich) +% - Reflexionen: optionaler 2-Tap Echo-Term (vereinfachtes Stub/Mismatch-Modell) +% +% Keine Toolboxes zwingend nötig (außer eye diagram -> optional). +%========================================================================== + +clear; close all; clc; + +%% ----------------------- Simulationseinstellungen ----------------------- +Rb = 10e9; % Bitrate [b/s] +UI = 1/Rb; % Unit Interval [s] +os = 32; % Oversampling +Fs = Rb*os; % Abtastrate [Hz] +Ts = 1/Fs; + +Nbits = 4096; % Anzahl Bits pro Lane +A = 1.0; % NRZ-Amplitude (0/1 -> +/-A) + +% PRBS-Seed +rng(1); + +%% ----------------------- Testsequenzen (Aggressor + Victim) ------------- +% Victim-Lane Daten +b_v = randi([0 1], Nbits, 1); +x_v = A*(2*b_v - 1); % NRZ: -A/+A + +% Aggressor-Lane Daten (unabhängig) +b_a = randi([0 1], Nbits, 1); +x_a = A*(2*b_a - 1); + +% Upsampling (NRZ-Rechteck) +x_v_up = upsample(x_v, os); +x_a_up = upsample(x_a, os); + +% Simple Zero-Order Hold (rechteckig halten) +x_v_up = filter(ones(os,1), 1, x_v_up); +x_a_up = filter(ones(os,1), 1, x_a_up); + +% Am Anfang Filter-Transient entfernen / zentrieren +x_v_up = x_v_up(os:end); +x_a_up = x_a_up(os:end); + +% Signal-Länge +N = length(x_v_up); + +% Optional: leichte TX-RiseTime-Abbildung (1st order LP) +% (hilft, unrealistische unendliche Bandbreite zu vermeiden) +f_tx_3dB = 0.35/(35e-12); % grob: 35ps 10-90% -> f3dB ~ 0.35/tr +[x_v_up, x_a_up] = apply_1pole_lp(x_v_up, x_a_up, Fs, f_tx_3dB); + +%% ----------------------- Frequenzachse (FFT) ---------------------------- +% FFT-Länge als Power-of-two (schneller, weniger Circular-Artefakte) +Nfft = 2^nextpow2(N*2); +f = (0:Nfft-1).'*Fs/Nfft; % 0 ... Fs*(1-1/Nfft) +% Für Transferfunktionen brauchen wir auch negative Frequenzen implizit: +% MATLAB-FFT nutzt diese Darstellung automatisch, solange H hermitesch ist. +% Wir definieren H für alle FFT-Bins konsistent mit reellen Signalen. + +%% ----------------------- RLGC-Modell (Dispersion & Loss) ----------------- +% Leitungslänge +l = 0.40; % [m] + +% Per-unit-length Parameter (Beispielwerte, grob PCB-ähnlich) +L = 3.2e-7; % [H/m] +C = 1.6e-10; % [F/m] +G0 = 0; % [S/m] (DC) +tanD = 0.012; % Verlustfaktor (vereinfachtes Dielektrikum-Modell) +% G(f) ~ omega*C*tanD +omega = 2*pi*f; +G = G0 + omega.*C*tanD; + +% Skin-Effekt / Rauheit: R(f) ~ Rdc + k*sqrt(f) +Rdc = 0.05; % [Ohm/m] +kR = 0.20; % [Ohm/m/sqrt(Hz)] * sqrt(Hz) -> Ohm/m +R = Rdc + kR*sqrt(max(f,1)); % max(..,1) vermeidet sqrt(0) + +% gamma(f) und Z0(f) +gamma = sqrt( (R + 1j*omega*L) .* (G + 1j*omega*C) ); % [1/m] +Z0 = sqrt( (R + 1j*omega*L) ./ (G + 1j*omega*C) ); % [Ohm] + +% Hauptpfad-Transferfunktion (perfekt terminiert) +H11 = exp(-gamma*l); + +%% ----------------------- Crosstalk-Modell (Aggressor -> Victim) ---------- +% Sehr vereinfachtes FEXT-ähnliches Modell: wächst mit f, sättigt +k_xt = 0.04; % Crosstalk-Stärke (skalieren) +fc = 8e9; % Eckfrequenz der Koppelcharakteristik + +% optional etwas andere Verluste auf dem Koppelpfad +gamma_xt = gamma .* (1.05 + 0.00j); + +Hxt_shape = (1j*(f/fc)) ./ (1 + 1j*(f/fc)); % ~ Highpass +H21 = k_xt .* Hxt_shape .* exp(-gamma_xt*l); + +%% ----------------------- Reflexions/Echo-Modell (optional) --------------- +% Einfache 2-Tap-Approximation: Hecho = 1 + a*exp(-j*2*pi*f*tau) +% (z.B. Stub/Via-Resonanz als Echo) +useEcho = true; +a_echo = 0.18; % Echo-Amplitude (|a|<1) +tau = 120e-12; % Echo-Delay [s] + +if useEcho + Hecho = 1 + a_echo*exp(-1j*2*pi*f*tau); + H11 = H11 .* Hecho; +end + +%% ----------------------- Kanal-Anwendung per FFT ------------------------- +% Zero-padding auf Nfft +Xv = fft(x_v_up, Nfft); +Xa = fft(x_a_up, Nfft); + +% Victim-Ausgang = Hauptpfad * Victim + Crosstalkpfad * Aggressor +Yv = H11 .* Xv + H21 .* Xa; + +y_v = real(ifft(Yv, Nfft)); +y_v = y_v(1:N); % zurück auf Originallänge + +%% ----------------------- Auswertung: Eye Diagramm ------------------------ +% Eye über viele UIs (ohne Toolbox: eigene Darstellung) +plot_eye(y_v, os, 2); % 2 UI Breite + +%% ----------------------- Auswertung: Frequenzgänge ----------------------- +figure('Name','Transferfunktionen'); +subplot(2,1,1); +plot(f/1e9, 20*log10(abs(H11)+1e-15)); grid on; +xlabel('f [GHz]'); ylabel('|H_{11}| [dB]'); +title('Hauptkanal (inkl. Dispersion/Loss + optional Echo)'); + +subplot(2,1,2); +plot(f/1e9, 20*log10(abs(H21)+1e-15)); grid on; +xlabel('f [GHz]'); ylabel('|H_{21}| [dB]'); +title('Crosstalk-Pfad (Aggressor -> Victim)'); + +%% ----------------------- Optional: Bitfehlersicht / Sampling ------------- +% Einfaches Sampling in der UI-Mitte (ohne Taktwiedergewinnung) +% (nur grobe Demo, kein echter CDR) +sampleOffset = round(os/2); +idx = sampleOffset:os:(sampleOffset + (Nbits-2)*os); +samples = y_v(idx); + +% Hard Decision +bh = samples > 0; +ber = mean(bh(:) ~= b_v(1:length(bh))); + +fprintf('--- Ergebnis ---\n'); +fprintf('Bitrate: %.2f Gb/s, Oversampling: %d, Fs: %.2f GHz\n', Rb/1e9, os, Fs/1e9); +fprintf('Echo: %d (a=%.3f, tau=%.1f ps)\n', useEcho, a_echo, tau*1e12); +fprintf('Crosstalk k_xt=%.3f, fc=%.2f GHz\n', k_xt, fc/1e9); +fprintf('Grobe BER (ohne CDR/EQ): %.3e\n', ber); + +%========================================================================== +% Lokale Hilfsfunktionen +%========================================================================== + +function [x1o, x2o] = apply_1pole_lp(x1, x2, Fs, f3dB) + % 1. Ordnung Tiefpass als TX-Bandbegrenzung (optional) + if isempty(f3dB) || f3dB<=0 + x1o = x1; x2o = x2; return; + end + w = 2*pi*f3dB; + a = exp(-w/Fs); + b = 1-a; + x1o = filter(b, [1 -a], x1); + x2o = filter(b, [1 -a], x2); +end + +function plot_eye(y, os, nUI) + % Eye-Diagramm ohne Toolbox: Segmente von nUI*os sampeln und überlagern + Lseg = nUI*os; + nSeg = floor(length(y)/os) - nUI - 2; + if nSeg < 10 + warning('Zu wenig Daten für Eye.'); + return; + end + + figure('Name','Eye Diagramm'); + hold on; grid on; + t = (0:Lseg-1)/os; % in UI + for k = 1:nSeg + i0 = (k-1)*os + 1; + seg = y(i0:i0+Lseg-1); + plot(t, seg); + end + xlabel('Zeit [UI]'); + ylabel('Amplitude [V (rel.)]'); + title(sprintf('Eye Diagramm (Breite = %d UI, Oversampling = %d)', nUI, os)); +end \ No newline at end of file diff --git a/projects/ASTRO Electrical Link/SBD_Transmission.m b/projects/ASTRO Electrical Link/SBD_Transmission.m new file mode 100644 index 0000000..f87e3c7 --- /dev/null +++ b/projects/ASTRO Electrical Link/SBD_Transmission.m @@ -0,0 +1,116 @@ +%% +clear all; +close all; + +% General Parameters +M = 2; +fsym = 20e9; +K_over = 8; +fs = K_over*fsym; +pulselength_mf = 16; +rolloff_mf = 0.1; +N = 18; +Pform = Pulseformer("fsym",fsym,"fdac",fs,"pulse","rc","pulselength",16,"alpha",0.1); +duob_mode = db_mode.no_db; + +% FFE Parameters +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; +ffe_order = [300, 0, 0]; + +% Toggles +hybrid = 1; +ctle = 1; +timing_recovery = 0; +ffe = 0; +plots = 1; + +% Generate Data +[a_1,a_1_Symbols,a_1_Tx_bits] = PAMsource(... + "fsym",fsym,"M",M,"order",N,"useprbs",0,... + "fs_out",fs,... + "applyclipping",0,"clipfactor",1.5,... + "applypulseform",1,"pulseformer",Pform,... + "randkey",1,... + 'duobinary_mode',duob_mode,... + "mrds_code",0,"mrds_blocklength",512).process(); + +[a_2,a_2_Symbols,a_2_Tx_bits] = PAMsource(... + "fsym",fsym,"M",M,"order",18,"useprbs",0,... + "fs_out",fs,... + "applyclipping",0,"clipfactor",1.5,... + "applypulseform",1,"pulseformer",Pform,... + "randkey",1,... + 'duobinary_mode',duob_mode,... + "mrds_code",0,"mrds_blocklength",512).process(); + +% S-Parameter Calculation +file_path = "C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\COM Test\Mellitz\host_pkg_top_50mm_max_skew_cable_module_pin_pad_fext1.s4p"; +plot_parameters = 1; +test = 0; +S_test = [0,1;1,0]; +[b_1, b_2] = Electrical_Trace_BiDi('file_path',file_path,'plot',plot_parameters,'test',test,'S_test',S_test).process(a_1,a_2); + +% Hybrid +if hybrid + [~,b_1,~] = Electrical_Hybrid("file_path",file_path).process(a_1,b_1); + [~,b_2,~] = Electrical_Hybrid("file_path",file_path).process(a_2,b_2); +end + +% CTLE +if ctle + b_2.normalize("mode","rms").plot("displayname",'Before CTLE','fignum',9); + b_2.normalize("mode","rms").spectrum("displayname",'Before CTLE','normalizeTo0dB',1,'fignum',10); + b_2.normalize("mode","rms").eye(fsym,2,'displayname','Before CTLE','fignum',11); + + b_2 = CTLE('Aac_dB',0,'Adc_dB',-5,'f_p1',7.5e9,'f_p2',10e9,'plot',1).process(b_2); + + b_2.normalize("mode","rms").plot("displayname",'After CTLE','fignum',9); + b_2.normalize("mode","rms").spectrum("displayname",'After CTLE','normalizeTo0dB',1,'fignum',10); + b_2.normalize("mode","rms").eye(fsym,2,'displayname','After CTLE','fignum',12); +end + +% Timing Recovery +if timing_recovery + b_1 = MaxVar_Timing_Recovery('mode',0,'fsym',fsym,'fadc',K_over*fsym,'num_tau',K_over*128,'sps',K_over,'comp_signal',0,'comp_mode',0).process(b_1); + b_2 = MaxVar_Timing_Recovery('mode',0,'fsym',fsym,'fadc',K_over*fsym,'num_tau',K_over*128,'sps',K_over,'comp_signal',0,'comp_mode',0).process(b_2); +end + +% FFE +if ffe + eq_ffe = EQ("Ne",ffe_order,"Nb",[0,0,0], ... + "training_length",len_tr,"training_loops",5,"dd_loops",5, ... + "K",1,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... + "FFEmu",0,"plotfinal",0,"ideal_dfe",1); + + output.ffe_results = ffe(eq_ffe,M,b_2,a_1_Symbols,a_1_Tx_bits, ... + "precode_mode",duob_mode,'showAnalysis',1,"postFFE",[], ... + "eth_style_symbol_mapping",0); + + output.ffe_results.metrics.print +end + +% Plots +if plots + % Time Domain Signals + a_1.normalize("mode","rms").plot("displayname",'1st Port Tx','fignum',3); + b_2.normalize("mode","rms").plot("displayname",'2nd Port Rx','fignum',3); + + a_2.normalize("mode","rms").plot("displayname",'2nd Port Tx','fignum',4); + b_1.normalize("mode","rms").plot("displayname",'1st Port Rx','fignum',4); + + % Spectra + a_1.normalize("mode","rms").spectrum("displayname",'1st Port Tx','normalizeTo0dB',1,'fignum',5); + b_2.normalize("mode","rms").spectrum("displayname",'2nd Port Rx','normalizeTo0dB',1,'fignum',5); + + a_2.normalize("mode","rms").spectrum("displayname",'2nd Port Tx','normalizeTo0dB',1,'fignum',6); + b_1.normalize("mode","rms").spectrum("displayname",'1st Port Rx','normalizeTo0dB',1,'fignum',6); + + figure;plot(xcorr(a_1.signal,b_2.signal)); + figure;plot(xcorr(a_2.signal,b_1.signal)); +end \ No newline at end of file diff --git a/projects/FSO_transmission/Evaluation Scripts/db_test_i_dont_get_it.m b/projects/FSO_transmission/Evaluation Scripts/db_test_i_dont_get_it.m new file mode 100644 index 0000000..5f5450d --- /dev/null +++ b/projects/FSO_transmission/Evaluation Scripts/db_test_i_dont_get_it.m @@ -0,0 +1,26 @@ +% Generate Random NRZ Signal +fsym = 50e9; +K_over = 4; +fs = K_over*fsym; +pulselength_mf = 16; +rolloff_mf = 0.1; +N = 2^18; + +test_signal = (2*randi([0 1], 1, N) - 1).'; +test_signal = Electricalsignal(test_signal,"fs",fsym); + +% test_signal_en = Duobinary().encode(test_signal); +% test_signal_de = Duobinary().decode(test_signal_en); +% +% test_signal_de.normalize("mode","rms").spectrum("displayname",'After Enc-Dec','fignum',2); +% test_signal_de.normalize("mode","rms").plot("displayname",'After Enc-Dec','fignum',3); + +test_signal_pre = Duobinary().precode(test_signal); + +test_signal.normalize("mode","rms").spectrum("displayname",'Before Pre','fignum',2); +test_signal.normalize("mode","rms").plot("displayname",'Before Pre','fignum',3); + +test_signal_pre.normalize("mode","rms").spectrum("displayname",'After Pre','fignum',2); +test_signal_pre.normalize("mode","rms").plot("displayname",'After Pre','fignum',3); + + diff --git a/projects/FSO_transmission/Evaluation Scripts/evalscript_ffe_with_and_without_db_target.m b/projects/FSO_transmission/Evaluation Scripts/evalscript_ffe_with_and_without_db_target.m new file mode 100644 index 0000000..e902578 --- /dev/null +++ b/projects/FSO_transmission/Evaluation Scripts/evalscript_ffe_with_and_without_db_target.m @@ -0,0 +1,74 @@ +%% +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 = 4; +trlength = 4096*4; +x = 225:10:285; +BER_PAM_2 = []; +post_only = 0; +filter_length_vec = 140; +num_signal = 11; +our_signal = 1; +weighted_DFE = 0; + +for eq_method = 1 + + for db_target = 0:1 + 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, filter_length_vec, num_signal, our_signal, weighted_DFE, db_target); + BER_PAM_2 = [BER_PAM_2, BER_run]; + end + save('C:\Users\magf\Desktop\Desktop\Projekte\FSO\Data\BER_PAM_4_FFE_Only_DB_Target_' + string(db_target) + '.mat', 'x', 'BER_PAM_2') + BER_PAM_2 = []; + end + +end + +%% +x = 225:10:285; +for k = 0:1 + BER = load('C:\Users\magf\Desktop\Desktop\Projekte\FSO\Data\BER_PAM_4_FFE_Only_DB_Target_' + string(k) + '.mat'); + BER = BER.BER_PAM_2; + + figure(202120) + plot(x, BER, '-o','LineWidth',1.75); + hold on +end +if M == 2 + old_BER = [-1.5, -1.75, -2, -2.3, -2.6, -2.25, -1.95]; +elseif M == 4 + old_BER = [-1.6, -1.85, -2.2, -2.45, -2.3, -2, -1.4]; +end +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 without DB Target - 200 Taps', 'FFE with DB Target - 200 Taps', 'BER Paper', ... + 'Interpreter','latex', ... + 'Location','southwest', 'FontSize', 14) + +% FEC Labels direkt im Plot +text(221,2.2e-2,'o-FEC','Color','k','FontSize', 14, 'Interpreter','latex') +text(221,3.3e-3,'HD-FEC','Color','b','FontSize', 14, 'Interpreter','latex') +text(221,5.4e-3,'KP4+Hamming','Color','g','FontSize', 14,'Interpreter','latex') +text(221,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 1]) +set(gca,'YScale','log') +% beautifyBERplot +hold off \ No newline at end of file diff --git a/projects/FSO_transmission/Evaluation Scripts/evalscript_ffe_with_and_without_db_target_baud_rate_sweep.m b/projects/FSO_transmission/Evaluation Scripts/evalscript_ffe_with_and_without_db_target_baud_rate_sweep.m new file mode 100644 index 0000000..4fa1e73 --- /dev/null +++ b/projects/FSO_transmission/Evaluation Scripts/evalscript_ffe_with_and_without_db_target_baud_rate_sweep.m @@ -0,0 +1,69 @@ +%% +current = 225:10:285; +current = string(current); + +power = [8.4, 8.4, 42.3, 8.4, 8.4]; +power = string(power); + +num_pf_coeff = 4; +taps_ffe = [200, 0, 0]; +taps_dfe = [0, 0, 0]; +M = 4; +trlength = 4096*4; +x = 4:1:8; +BER_PAM_2 = []; +post_only = 0; +filter_length_vec = 140; +num_signal = 11; +our_signal = 1; +weighted_DFE = 0; + +for eq_method = 1 + + for db_target = 0:1 + for baud_rate = 4:1:8 + BER_run = first_analysis_baud_rate_sweep(string(baud_rate), power(baud_rate-3), num_pf_coeff, taps_ffe, taps_dfe, M, trlength, eq_method, filter_length_vec, num_signal, db_target); + BER_PAM_2 = [BER_PAM_2, BER_run]; + end + save('C:\Users\magf\Desktop\Desktop\Projekte\FSO\Data\BER_PAM_4_FFE_Only_DB_Target_Baud_Rate_Sweep' + string(db_target) + '.mat', 'x', 'BER_PAM_2') + BER_PAM_2 = []; + end + +end + +%% +x = 4:1:8; +for k = 0:1 + BER = load('C:\Users\magf\Desktop\Desktop\Projekte\FSO\Data\BER_PAM_4_FFE_Only_DB_Target_Baud_Rate_Sweep' + string(k) + '.mat'); + BER = BER.BER_PAM_2; + + figure(202120) + 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 without DB Target - 200 Taps', 'FFE with DB Target - 200 Taps', ... + 'Interpreter','latex', ... + 'Location','southwest', 'FontSize', 14) + +% FEC Labels direkt im Plot +text(221,2.2e-2,'o-FEC','Color','k','FontSize', 14, 'Interpreter','latex') +text(221,3.3e-3,'HD-FEC','Color','b','FontSize', 14, 'Interpreter','latex') +text(221,5.4e-3,'KP4+Hamming','Color','g','FontSize', 14,'Interpreter','latex') +text(221,2.4e-4,'KP4','Color','r','FontSize', 14,'Interpreter','latex') + +xlabel('Symbol Rate [GBd]', 'Interpreter','latex') +ylabel('BER', 'Interpreter','latex') +% title('BER for PAM-2', 'Interpreter','latex') + +grid minor +ylim([1e-2 1]) +set(gca,'YScale','log') +% beautifyBERplot +hold off \ No newline at end of file diff --git a/projects/FSO_transmission/Evaluation Scripts/evalscript_pam_2_20G.m b/projects/FSO_transmission/Evaluation Scripts/evalscript_pam_2_20G.m new file mode 100644 index 0000000..26a63e0 --- /dev/null +++ b/projects/FSO_transmission/Evaluation Scripts/evalscript_pam_2_20G.m @@ -0,0 +1,149 @@ +%% +M = 4; +current = 255; +power = [8.4, 8.4, 42.3, 8.4, 8.4, 8.4, 8.4]; +rolloff = 0.6; +baudrate = 4:1:10; +baudrate_e = ["4e9", "5e9", "6e9", "7e9", "8e9", "9e9", "1e10"]; + +num_pf_coeff = 4; +taps_ffe = [300, 0, 0]; +taps_dfe = [0, 0, 0]; +trlength = 4096*4; +filter_length = 210; + +x = 4:1:10; +num_currents = length(x); +num_signals = 18; +num_eq_methods = 4; +BER_PAM_2 = zeros(num_currents,num_signals); +save_name = "BER_PAM_4_Waitbar_Test"; + +% ---------------- Progressbar Setup ---------------- +totalRuns = num_eq_methods * num_currents * num_signals; +runCount = 0; + +h = waitbar(0, 'Starte Simulation...', 'Name', 'FSO Simulation Progress'); +tStart = tic; + +% Optional: nur alle N Schritte UI updaten (reduziert Overhead) +updateEvery = 1; % z.B. 10 setzen, wenn es extrem viele Iterationen werden +% --------------------------------------------------- + +try + for eq_method = 1:1:4 + num_current_iteration = 1; + for i = 1:num_currents + for num_signal = 1:1:num_signals + + % --- Dein eigentlicher Run --- + try + BER_run = first_analysis_fso(current, power(i), baudrate(i), baudrate_e(i), rolloff, ... + num_pf_coeff, taps_ffe, taps_dfe, M, trlength, eq_method, filter_length, num_signal); + catch + BER_run = NaN; + end + BER_PAM_2(num_current_iteration, num_signal) = BER_run; + % ----------------------------- + + % --- Progress updaten --- + runCount = runCount + 1; + + if mod(runCount, updateEvery) == 0 || runCount == 1 || runCount == totalRuns + frac = runCount / totalRuns; + + elapsed = toc(tStart); + if frac > 0 + remaining = elapsed * (1/frac - 1); + else + remaining = NaN; + end + + msg = sprintf(['EQ %d/4 | i %d/%d | Signal %d/%d\n' ... + 'Gesamt %d/%d (%.1f %%) | ETA ~ %.1f min'], ... + eq_method, i, num_currents, num_signal, num_signals, ... + runCount, totalRuns, 100*frac, remaining/60); + + waitbar(frac, h, msg); + + % Falls jemand das waitbar-Fenster schließt: sauber abbrechen + if ~ishandle(h) + error('Progressbar wurde geschlossen. Abbruch durch Nutzer.'); + end + end + % ------------------------ + + end + num_current_iteration = num_current_iteration + 1; + end + + save('C:\Users\magf\Desktop\Desktop\Projekte\FSO\Data\' + save_name + "_" + string(eq_method) + '.mat', 'x', 'BER_PAM_2') + BER_PAM_2 = zeros(num_currents, num_signals); + end + +catch ME + % Falls irgendwas schiefgeht: waitbar schließen und Fehler weiterwerfen + if exist('h','var') && ishandle(h) + close(h); + end + rethrow(ME); +end + +% Clean exit +if exist('h','var') && ishandle(h) + close(h); +end + +%% +x = 4:1:10; +save_name = "BER_PAM_4_Waitbar_Test"; +figure(202120); clf; hold on; grid on; +for k = 1:4 + BER_min = zeros(1,length(x)); + BER_avg = zeros(1,length(x)); + BER_max = zeros(1,length(x)); + + tmp = load('C:\Users\magf\Desktop\Desktop\Projekte\FSO\Data\' + save_name + "_" + string(k) + '.mat'); + BER = tmp.BER_PAM_2; + + for i = 1:length(x) + BERs = BER(i,:); + BER_min(i) = min(BERs); + BER_avg(i) = mean(BERs); + BER_max(i) = max(BERs); + end + + err_low = BER_avg - BER_min; + err_high = BER_max - BER_avg; + errorbar(x, BER_avg, err_low, err_high, '-o', 'LineWidth', 1.75); +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(221,2.2e-2,'o-FEC','Color','k','FontSize', 14, 'Interpreter','latex') +text(221,3.3e-3,'HD-FEC','Color','b','FontSize', 14, 'Interpreter','latex') +text(221,5.4e-3,'KP4+Hamming','Color','g','FontSize', 14,'Interpreter','latex') +text(231,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 index 20b7b72..70a615c 100644 --- a/projects/FSO_transmission/Evaluation Scripts/evalscript_pam_4.m +++ b/projects/FSO_transmission/Evaluation Scripts/evalscript_pam_4.m @@ -5,7 +5,7 @@ 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]; +% taps_dfe = [0, 0, 0]; M = 4; trlength = 4096*4; x = 225:10:285; @@ -13,20 +13,25 @@ BER_PAM_4 = []; filter_length = 210; num_signal = 11; our_signal = 1; +weighted_DFE = [1, 0.1, 0.1]; -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, filter_length, num_signal, our_signal); - BER_PAM_4 = [BER_PAM_4, BER_run]; - save_and_append('meineDB.sqlite', 'BER_PAM_4_Save_and_Append', eq_method, num_signal, BER_run, x) +for m = 0:5:40 + for eq_method = 2 + for j = 4 + BER_run = first_analysis_ber(current(j), power(j), num_pf_coeff, taps_ffe, [m, 0, 0], M, trlength, eq_method, filter_length, num_signal, our_signal, weighted_DFE); + BER_PAM_4 = [BER_PAM_4, BER_run]; + % save_and_append('meineDB.sqlite', 'BER_PAM_4_Save_and_Append', eq_method, num_signal, BER_run, x) + end end - BER_PAM_4 = []; end +save('C:\Users\magf\Desktop\Desktop\Projekte\FSO\Data\BER_PAM_4_Weighted_DFE_' + string(k) + '.mat'); +BER_PAM_4 = []; + %% x = 0:5:40; for k = 2 - BER = load('C:\Users\magf\Desktop\Desktop\Projekte\FSO\Data\BER_PAM_4_DFE_Tap_Sweep_' + string(k) + '.mat'); + BER = load('C:\Users\magf\Desktop\Desktop\Projekte\FSO\Data\BER_PAM_4_Weighted_DFE_' + string(k) + '.mat'); BER_values = BER.BER_PAM_4; figure(202120) @@ -53,7 +58,7 @@ text(23,2.6e-3,'HD-FEC','Color','b','FontSize', 14, 'Interpreter','latex') text(23.5,6.5e-3,'KP4+Hamming','Color','g','FontSize', 14,'Interpreter','latex') text(23,1.4e-4,'KP4','Color','r','FontSize', 14,'Interpreter','latex') -xlabel('Laser Bias Current [mA]', 'Interpreter','latex', 'FontSize', 14) +xlabel('Number of First Order DFE Taps', 'Interpreter','latex', 'FontSize', 14) ylabel('BER', 'Interpreter','latex', 'FontSize', 14) % title('BER for PAM-4', 'Interpreter','latex') diff --git a/projects/FSO_transmission/Evaluation Scripts/evalscript_tests.m b/projects/FSO_transmission/Evaluation Scripts/evalscript_tests.m new file mode 100644 index 0000000..32857dd --- /dev/null +++ b/projects/FSO_transmission/Evaluation Scripts/evalscript_tests.m @@ -0,0 +1,18 @@ +%% +M = 2; +current = 265; +power = 8.4; +rolloff = 0.6; +baudrate = 20; +baudrate_e = "2e10"; + +num_pf_coeff = 4; +taps_ffe = [200, 0, 0]; +taps_dfe = [0, 0, 0]; +trlength = 4096*2; +filter_length = 140; +num_signal = 1; +eq_method = 1; + +BER_run = first_analysis_fso(current, power, baudrate, baudrate_e, rolloff, ... + num_pf_coeff, taps_ffe, taps_dfe, M, trlength, eq_method, filter_length, num_signal); diff --git a/projects/FSO_transmission/Evaluation Scripts/first_analysis_20G.m b/projects/FSO_transmission/Evaluation Scripts/first_analysis_20G.m new file mode 100644 index 0000000..4c7cd36 --- /dev/null +++ b/projects/FSO_transmission/Evaluation Scripts/first_analysis_20G.m @@ -0,0 +1,222 @@ +function BER_value = first_analysis_fso(current, power, rolloff, num_pf_coeff, taps_ffe, taps_dfe, M, trlength, eq_method, filter_length, num_signal, our_signal, weighted_DFE_mode, db_target) + %% + 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, "20G_PAM2\tx_info\tx_info_PAM2_20Gbd" + rolloff + "RRC.mat"); + filename = fullfile(base, "20G_PAM2\M=2_Rs=2e10_Fs=8e10_I=" + current + "mA_RoP=" + power + "mW_L=31m_PS=RRC_rolloff=" + rolloff + "_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=" + current + "mA_RoP=" + power + "mW_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) + + 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{num_signal}; + + 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{num_signal}; + + % 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 = Rx_Time_Rec.normalize('mode','rms'); + Rx_tr_mf = Rx_tr_mf.normalize('mode','rms'); + + % 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.. + 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 = 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",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); + + ffe_results = ffe_db(eq_ffe,M,Rx_synced,Symbols,Bits, ... + "precode_mode",duob_mode,'showAnalysis',0,"postFFE",[], ... + "eth_style_symbol_mapping",mapping_style,'db_target',db_target); + + % 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",0,'weighted_DFE',0,'weighted_DFE_d_min',0.5, ... + 'weighted_DFE_mode','I2','weighted_DFE_I_mode',weighted_DFE_mode,'PDFE_coefficient',0.01); + 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,"decoding_mode","memoryless"); + + 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",filter_length,"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_baud_rate_sweep.m b/projects/FSO_transmission/Evaluation Scripts/first_analysis_baud_rate_sweep.m index 26d6ca2..23831ee 100644 --- a/projects/FSO_transmission/Evaluation Scripts/first_analysis_baud_rate_sweep.m +++ b/projects/FSO_transmission/Evaluation Scripts/first_analysis_baud_rate_sweep.m @@ -1,4 +1,4 @@ -function [BER, Channel_Alpha] = first_analysis_baud_rate_sweep(baud_rate, power, num_pf_coeff, taps_ffe, taps_dfe, M, trlength, method, filter_length, num_signal) +function [BER, Channel_Alpha] = first_analysis_baud_rate_sweep(baud_rate, power, num_pf_coeff, taps_ffe, taps_dfe, M, trlength, method, filter_length, num_signal, db_target) %% close all base = "C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\FSO Equalizer\Sweep Data\"; @@ -138,14 +138,14 @@ function [BER, Channel_Alpha] = first_analysis_baud_rate_sweep(baud_rate, power, if method == 1 %% -------------------- FFE -------------------- % requires some more digging what is going on :-) - eq_ffe = EQ("Ne",[500, 0, 0],"Nb",[0, 0, 0], ... + eq_ffe = 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); ffe_results = ffe(eq_ffe,M,Rx_synced,Symbols,Bits, ... "precode_mode",duob_mode,'showAnalysis',0,"postFFE",[], ... - "eth_style_symbol_mapping",mapping_style); + "eth_style_symbol_mapping",mapping_style,'db_target',db_target); % ffe_results.metrics.print fprintf('My EQ: %.1e \n',ffe_results.metrics.BER); diff --git a/projects/FSO_transmission/Evaluation Scripts/first_analysis_ber.m b/projects/FSO_transmission/Evaluation Scripts/first_analysis_ber.m index 3c4abb3..dc8fc0c 100644 --- a/projects/FSO_transmission/Evaluation Scripts/first_analysis_ber.m +++ b/projects/FSO_transmission/Evaluation Scripts/first_analysis_ber.m @@ -1,4 +1,4 @@ -function BER_value = first_analysis_ber(current, power, num_pf_coeff, taps_ffe, taps_dfe, M, trlength, eq_method, filter_length, num_signal, our_signal, weighted_DFE) +function BER_value = first_analysis_ber(current, power, num_pf_coeff, taps_ffe, taps_dfe, M, trlength, eq_method, filter_length, num_signal, our_signal, weighted_DFE, db_target) %% close all @@ -154,14 +154,14 @@ function BER_value = first_analysis_ber(current, power, num_pf_coeff, taps_ffe, %% -------------------- FFE -------------------- % requires some more digging what is going on :-) - eq_ffe = EQ("Ne",[500, 0, 0],"Nb",[0, 0, 0], ... + eq_ffe = 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); ffe_results = ffe(eq_ffe,M,Rx_synced,Symbols,Bits, ... "precode_mode",duob_mode,'showAnalysis',0,"postFFE",[], ... - "eth_style_symbol_mapping",mapping_style); + "eth_style_symbol_mapping",mapping_style,'db_target',db_target); % ffe_results.metrics.print fprintf('My EQ: %.1e \n',ffe_results.metrics.BER); diff --git a/projects/FSO_transmission/Evaluation Scripts/first_analysis_fso.m b/projects/FSO_transmission/Evaluation Scripts/first_analysis_fso.m new file mode 100644 index 0000000..a8eae50 --- /dev/null +++ b/projects/FSO_transmission/Evaluation Scripts/first_analysis_fso.m @@ -0,0 +1,180 @@ +function BER_value = first_analysis_fso(current, power, baudrate, baudrate_e, rolloff, ... + num_pf_coeff, taps_ffe, taps_dfe, M, trlength, eq_method, filter_length, num_signal) + + %% Close Remaining Figures + close all + + %% Load Data + base = "C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\FSO Equalizer\Sweep Data\"; + mode = 0; %0 oder 1 + + tx_data_path = fullfile(base, string(baudrate) + "G_PAM" + string(M) + "\tx_info\tx_info_PAM" + string(M) + "_" + string(baudrate) + "Gbd" + string(rolloff) + "RRC.mat"); + filename = fullfile(base, string(baudrate) + "G_PAM" + string(M) + "\M=" + string(M) + "_Rs=" + baudrate_e + "_Fs=8e10_I=" + string(current) + "mA_RoP=" + string(power) + "mW_L=31m_PS=RRC_rolloff=" + string(rolloff) + "_Mode=Rise.mat"); + + 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')); + + %% Generate 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); + + %% Matched Filter + 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); + + % Apply Matched Filter + 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 Synchronization + [~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_matched.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1); + Rx_matched_1 = Rx_synced_cell{num_signal}; + + %% Timing Recovery + 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 = Rx_Time_Rec.normalize('mode','rms'); + + %% DSP Parameters + Rx_synced = Rx_Time_Rec; + 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; + + 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 + + %% EQ Methods + 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",0,'weighted_DFE',0,'weighted_DFE_d_min',0.5, ... + 'weighted_DFE_mode','I2','weighted_DFE_I_mode',0,'PDFE_coefficient',0.01); + 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,"decoding_mode","sequencedetection"); + + 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",filter_length,"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_tr_mf.m b/projects/FSO_transmission/Evaluation Scripts/first_analysis_tr_mf.m index aadbf1f..5add8a4 100644 --- a/projects/FSO_transmission/Evaluation Scripts/first_analysis_tr_mf.m +++ b/projects/FSO_transmission/Evaluation Scripts/first_analysis_tr_mf.m @@ -189,7 +189,7 @@ for our_signal = 1 end % Rx_synced = Rx_Time_Rec; % Rx_synced = Rx_synced_cell{1}; - len_tr = 4096*4; + len_tr = 4096*2; mu_ffe1 = 0.0001; mu_ffe2 = 0.0008; mu_ffe3 = 0.001; @@ -212,7 +212,7 @@ for our_signal = 1 end %% -------------------- FFE -------------------- - % requires some more digging what is going on :-) + % % 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, ... @@ -222,7 +222,7 @@ for our_signal = 1 % % ffe_results = ffe(eq_ffe,M,Rx_synced,Symbols,Bits, ... % "precode_mode",duob_mode,'showAnalysis',0,"postFFE",[], ... - % "eth_style_symbol_mapping",mapping_style); + % "eth_style_symbol_mapping",mapping_style,'db_target',1); % % if our_signal % fprintf('Our signal: %.1e \n',ffe_results.metrics.BER); @@ -233,29 +233,43 @@ for our_signal = 1 % end %% -------------------- VNLE + MLSE -------------------- - pf_ncoeffs = 4; - eq_v = EQ("Ne",[300, 0, 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",0, ... - 'weighted_DFE',1,'weighted_DFE_d_min',0.5,'weighted_DFE_mode','I2','weighted_DFE_I_mode',[1,0.1,0.1], ... - 'PDFE_coefficient',0.01); - % 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); + % 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",0, ... + % 'weighted_DFE',2,'weighted_DFE_d_min',0.9,'weighted_DFE_mode','R1','weighted_DFE_I_mode',[1,0.1,0.1], ... + % 'PDFE_coefficient',0.01); + % % 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 - [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); + %% -------------------- DB target -------------------- + mlse_db_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",M,'trellis_states',PAMmapper(M,0).levels); - 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 + eq_ = 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); + + 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,'decoding_mode',0); + + 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; %% -------------------- ML-based MLSE (L=2) -------------------- % ml_mlse_equalizer = ML_MLSE("epochs_tr",100,"epochs_dd",1, ... diff --git a/projects/FSO_transmission/Evaluation Scripts/plotscript_pam_2_conf.m b/projects/FSO_transmission/Evaluation Scripts/plotscript_pam_2_conf.m new file mode 100644 index 0000000..4e65995 --- /dev/null +++ b/projects/FSO_transmission/Evaluation Scripts/plotscript_pam_2_conf.m @@ -0,0 +1,83 @@ +figure(202120); clf; hold on +set(gcf, 'WindowState', 'maximized'); + +x = 225:10:285; +shift = 0.5; +x = [x;x-shift;x+shift;x;x]; + +hErr = gobjects(1,4); % Errorbars -> für Legende +hLine = gobjects(1,4); % Linien durch Mittelwerte (nicht in Legende) + +for k = 1:4 + BER = load("C:\Users\magf\Desktop\Desktop\Projekte\FSO\Data\BER_PAM_2_Optimal" + string(k) + ".mat"); + BER_all = BER.BER_num_signal_sweep_matrix; + BER_traces = BER_all((k-1)*7+1:k*7,:); + + BER_min = zeros(1,7); BER_avg = zeros(1,7); BER_max = zeros(1,7); + for l = 1:7 + BER_min(l) = min(BER_traces(l,:)); + BER_avg(l) = mean(BER_traces(l,:)); + BER_max(l) = max(BER_traces(l,:)); + end + + % Log-sicherer Floor + yFloor = 1e-12; + BER_min(BER_min <= 0 | isnan(BER_min) | isinf(BER_min)) = yFloor; + BER_max(BER_max <= 0 | isnan(BER_max) | isinf(BER_max)) = yFloor; + BER_max = max(BER_max, BER_min); + + % Errorbar-Parameter (asymmetrisch) + y = BER_avg; + yneg = BER_avg - BER_min; + ypos = BER_max - BER_avg; + + % Errorbars (Punkt + Intervall) + hErr(k) = errorbar(x(k,:), y, yneg, ypos, 'o', ... + 'LineWidth', 2.75, 'MarkerSize', 8, 'CapSize', 10); + + % Gestrichelte Linie durch die Mittelwerte (gleiche Farbe) + hLine(k) = plot(x(k,:), y, '--o', ... + 'LineWidth', 2.75, ... + 'MarkerSize', 8, ... + 'Color', hErr(k).Color); +end + +% Paper-Kurve +old_BER = 10.^([-1.5, -1.75, -2, -2.3, -2.6, -2.25, -1.95]); +hPaper = plot(x(k+1,:), old_BER, '-o', 'LineWidth', 2.75, 'MarkerSize', 8); + +% FEC-Linien +yline(2e-2, ':k', 'LineWidth',1.2); +yline(3.8e-3,':b', 'LineWidth',1.2); +yline(4.85e-3,':g','LineWidth',1.2); +yline(2.2e-4,':r', 'LineWidth',1.2); + +% FEC Labels (FontSize = 22) +text(240, 2.3e-2,'o-FEC', 'Color','k','FontSize',22,'Interpreter','latex') +text(222, 3e-3,'HD-FEC', 'Color','b','FontSize',22,'Interpreter','latex') +text(221, 5.7e-3,'KP4+Hamming', 'Color','g','FontSize',22,'Interpreter','latex') +text(232, 2.5e-4,'KP4', 'Color','r','FontSize',22,'Interpreter','latex') + +% Achsen +set(gca,'YScale','log'); grid minor +ylim([1e-4 5e-1]) + +% Schriftgrößen Achsen +ax = gca; +ax.FontSize = 22; +ax.LabelFontSizeMultiplier = 1; + +xlabel('Laser Bias Current [mA]', 'Interpreter','latex', 'FontSize', 22) +ylabel('BER', 'Interpreter','latex', 'FontSize', 22) + +% Legende: Errorbar-Handles verwenden (Marker + vertikale Linie in der Legende) +lgd = legend([hErr(1) hErr(2) hErr(3) hErr(4) hPaper], ... + {'FFE','FFE+PF+MLSE','DB','ML-MLSE','BER Paper'}, ... + 'Interpreter','latex','Location','southwest'); +lgd.FontSize = 22; + +% Figure-Style +set(gcf, 'Color', 'w'); +set(gcf, 'Renderer', 'painters'); + +hold off \ No newline at end of file diff --git a/projects/FSO_transmission/Evaluation Scripts/plotscript_pam_4_baud_rate_sweep_conf.m b/projects/FSO_transmission/Evaluation Scripts/plotscript_pam_4_baud_rate_sweep_conf.m new file mode 100644 index 0000000..1172913 --- /dev/null +++ b/projects/FSO_transmission/Evaluation Scripts/plotscript_pam_4_baud_rate_sweep_conf.m @@ -0,0 +1,80 @@ +figure(202120); clf; hold on +set(gcf, 'WindowState', 'maximized'); % automatisch Vollbild + +x = 4:1:8; +shift = 0.05; +x = [x;x-shift;x;x+shift]; + +hErr = gobjects(1,4); % Errorbars -> für Legende +hLine = gobjects(1,4); % Linien durch Mittelwerte (nicht in Legende) + +for k = 1:4 + BER = load("C:\Users\magf\Desktop\Desktop\Projekte\FSO\Data\BER_PAM_4_Baud_Rate_Sweep_Optimal_" + string(k) + ".mat"); + BER_all = BER.BER_num_signal_sweep_matrix; + BER_traces = BER_all((k-1)*5+1:k*5,:); + + BER_min = zeros(1,5); BER_avg = zeros(1,5); BER_max = zeros(1,5); + for l = 1:5 + BER_min(l) = min(BER_traces(l,:)); + BER_avg(l) = mean(BER_traces(l,:)); + BER_max(l) = max(BER_traces(l,:)); + end + + % Floor für Log-Plot (muss > 0 sein) + yFloor = 1e-7; + BER_min(BER_min <= 0 | isnan(BER_min) | isinf(BER_min)) = yFloor; + BER_max(BER_max <= 0 | isnan(BER_max) | isinf(BER_max)) = yFloor; + BER_max = max(BER_max, BER_min); + + % Errorbar-Parameter (asymmetrisch) + y = BER_avg; + yneg = BER_avg - BER_min; + ypos = BER_max - BER_avg; + + % Errorbars (Punkt + Intervall) + hErr(k) = errorbar(x(k,:), y, yneg, ypos, 'o', ... + 'LineWidth', 2.5, 'MarkerSize', 8, 'CapSize', 10); + + % Gestrichelte Linie durch die Mittelwerte (gleiche Farbe) + hLine(k) = plot(x(k,:), y, '--o', ... + 'LineWidth', 2.5, ... + 'MarkerSize', 8, ... + 'Color', hErr(k).Color); +end + +% FEC-Linien +yline(2e-2, ':k', 'LineWidth',1.2); +yline(3.8e-3,':b', 'LineWidth',1.2); +yline(4.85e-3,':g','LineWidth',1.2); +yline(2.2e-4,':r', 'LineWidth',1.2); + +% FEC Labels (bleiben bei 14) +text(3.1,2.4e-2,'o-FEC','Color','k','FontSize',22,'Interpreter','latex') +text(3.1,2.7e-3,'HD-FEC','Color','b','FontSize',22,'Interpreter','latex') +text(3.1,6.3e-3,'KP4+Hamming','Color','g','FontSize',22,'Interpreter','latex') +text(3.1,2.7e-4,'KP4','Color','r','FontSize',22,'Interpreter','latex') + +% Achsen +set(gca,'YScale','log'); grid minor +ylim([5e-7 5e-1]) +xlim([3 9]) + +% ---- Schriftgrößen: Achsen + Legende = 22 ---- +ax = gca; +ax.FontSize = 22; +ax.LabelFontSizeMultiplier = 1; + +xlabel('Baud Rate [GBd]', 'Interpreter','latex', 'FontSize', 22) +ylabel('BER', 'Interpreter','latex', 'FontSize', 22) + +% ---- Legende: Errorbar-Handles verwenden (Marker + vertikale Linie) ---- +lgd = legend([hErr(1) hErr(2) hErr(3) hErr(4)], ... + {'FFE','FFE+PF+MLSE','DB','ML-MLSE'}, ... + 'Interpreter','latex','Location','southwest'); +lgd.FontSize = 22; + +% ---- Figure "Presentation"-Style ---- +set(gcf, 'Color', 'w'); +set(gcf, 'Renderer', 'painters'); + +hold off \ No newline at end of file diff --git a/projects/FSO_transmission/Evaluation Scripts/plotscript_pam_4_conf.m b/projects/FSO_transmission/Evaluation Scripts/plotscript_pam_4_conf.m new file mode 100644 index 0000000..9a60152 --- /dev/null +++ b/projects/FSO_transmission/Evaluation Scripts/plotscript_pam_4_conf.m @@ -0,0 +1,84 @@ +figure(202120); clf; hold on +set(gcf, 'WindowState', 'maximized'); + +x = 225:10:285; +shift = 0.5; +x = [x;x-shift;x+shift;x;x]; + +hErr = gobjects(1,4); % Errorbars -> für Legende +hLine = gobjects(1,4); % Linien durch Mittelwerte (nicht in Legende) + +for k = 1:4 + BER = load("C:\Users\magf\Desktop\Desktop\Projekte\FSO\Data\BER_PAM_4_Optimal" + string(k) + ".mat"); + BER_all = BER.BER_num_signal_sweep_matrix; + BER_traces = BER_all((k-1)*7+1:k*7,:); + + BER_min = zeros(1,7); BER_avg = zeros(1,7); BER_max = zeros(1,7); + for l = 1:7 + BER_min(l) = min(BER_traces(l,:)); + BER_avg(l) = mean(BER_traces(l,:)); + BER_max(l) = max(BER_traces(l,:)); + end + + % Log-sicherer Floor (falls min=0) + yFloor = 1e-12; + BER_min(BER_min <= 0 | isnan(BER_min) | isinf(BER_min)) = yFloor; + BER_max(BER_max <= 0 | isnan(BER_max) | isinf(BER_max)) = yFloor; + BER_max = max(BER_max, BER_min); + + % Errorbar-Parameter (asymmetrisch) + y = BER_avg; + yneg = BER_avg - BER_min; + ypos = BER_max - BER_avg; + + % Errorbars (Punkt + Intervall) + hErr(k) = errorbar(x(k,:), y, yneg, ypos, 'o', ... + 'LineWidth', 2.75, 'MarkerSize', 8, 'CapSize', 10); + + % Gestrichelte Linie durch die Mittelwerte (gleiche Farbe) + hLine(k) = plot(x(k,:), y, '--o', ... + 'LineWidth', 2.75, ... + 'MarkerSize', 8, ... + 'Color', hErr(k).Color); +end + +% Paper-Kurve +old_BER = 10.^([-1.6, -1.85, -2.2, -2.45, -2.3, -2, -1.4]); +hPaper = plot(x(k+1,:), old_BER, '-o', 'LineWidth', 2.75, 'MarkerSize', 8); + +% FEC-Linien +yline(2e-2, ':k', 'LineWidth',1.2); +yline(3.8e-3,':b', 'LineWidth',1.2); +yline(4.85e-3,':g','LineWidth',1.2); +yline(2.2e-4,':r', 'LineWidth',1.2); + +% FEC Labels (wie bei dir, FontSize = 16) +text(221, 2.3e-2,'o-FEC', 'Color','k','FontSize',22,'Interpreter','latex') +text(221, 3e-3,'HD-FEC', 'Color','b','FontSize',22,'Interpreter','latex') +text(221, 6e-3,'KP4+Hamming', 'Color','g','FontSize',22,'Interpreter','latex') +text(221, 2.5e-4,'KP4', 'Color','r','FontSize',22,'Interpreter','latex') + +% Achsen +set(gca,'YScale','log'); grid minor +ylim([1e-5 5e-1]) + +% ---- Schriftgrößen: Achsen + Legende = 22 ---- +ax = gca; +ax.FontSize = 22; +ax.LabelFontSizeMultiplier = 1; + +xlabel('Laser Bias Current [mA]', 'Interpreter','latex', 'FontSize', 22) +ylabel('BER', 'Interpreter','latex', 'FontSize', 22) + +% Legende: Errorbar-Handles verwenden (Marker + vertikale Linie) +lgd = legend([hErr(1) hErr(2) hErr(3) hErr(4) hPaper], ... + {'FFE','FFE+PF+MLSE','DB','ML-MLSE','BER Paper'}, ... + 'Interpreter','latex','Location','southwest'); +lgd.FontSize = 22; + +% ---- Figure "Presentation"-Style ---- +set(gcf, 'Color', 'w'); +set(gcf, 'Units', 'pixels', 'Position', [100 100 1400 900]); +set(gcf, 'Renderer', 'painters'); + +hold off \ No newline at end of file diff --git a/projects/IMDD_base_system/imdd_it.m b/projects/IMDD_base_system/imdd_it.m index 887532b..cd451f0 100644 --- a/projects/IMDD_base_system/imdd_it.m +++ b/projects/IMDD_base_system/imdd_it.m @@ -1,4 +1,4 @@ - +close all; if 1 @@ -8,23 +8,26 @@ if 1 % uloops.laser_wavelength = [1293,1297.5,1302,1306.5,1310,1313.4,1318,1322.7,1327.4]; uloops.laser_wavelength = [1293]; uloops.M = [4]; - uloops.link_length = [0:2:10]; % 1,2,3,5,6,8,10 + uloops.link_length = 1; + % uloops.link_length = [0:2:10]; % 1,2,3,5,6,8,10 uloops.alpha = [0]; - uloops.duob_mode = db_mode.no_db; - + uloops.duob_mode = db_mode.db_precoded; + uloops.decoding_mode = "memoryless"; + wh = DataStorage(uloops); wh.addStorage("ber"); - + wh = submit_handle(@imdd_model,wh,"parallel",0); - + end %% figure hold on for alpha = uloops.alpha - a=wh.getStoValue('ber',1, [300].*1e9 , 1293, 4, uloops.link_length,alpha); - ffe = cellfun(@(x) x.ffe_results.metrics.BER, a); + a=wh.getStoValue('ber',1, [300].*1e9 , 1293, 4, uloops.link_length,alpha,uloops.duob_mode,uloops.decoding_mode); + % ffe = cellfun(@(x) x.ffe_results.metrics.BER, a); + ffe = cellfun(@(x) x.dbt_results.metrics.BER, a); plot(uloops.link_length,ffe,'DisplayName',sprintf('Alpha: %d',alpha),'LineStyle','-','HandleVisibility','on'); end @@ -38,20 +41,20 @@ ylabel('BER'); -% +% % wh_ana = wh_master; -% +% % cols = cbrewer2('Paired',8); -% +% % figure() -% +% % for precomp = [0,1] % wavelength=uloops.laser_wavelength; % for m = [6] -% +% % baudrate = wh_ana.parameter.bitrate.values; -% -% +% +% % %VNLE % precode = 1; % a = wh_ana.getStoValue('ber',precomp, precode, baudrate , wavelength, m, uloops.link_length); @@ -60,7 +63,7 @@ ylabel('BER'); % ber_mlse_pc = cellfun(@(x) x.vnle_pf_package{1,1}.ber_mlse, a); % %DB % ber_dbtgt_pc = cellfun(@(x) x.dbtgt_package{1,1}.ber, a); -% +% % precode = 0; % a = wh_ana.getStoValue('ber',precomp, precode, baudrate , wavelength, m, uloops.link_length); % ber_vnle = cellfun(@(x) x.vnle_pf_package{1,1}.ber_vnle, a); @@ -68,16 +71,16 @@ ylabel('BER'); % ber_mlse = cellfun(@(x) x.vnle_pf_package{1,1}.ber_mlse, a); % %DB % ber_dbtgt = cellfun(@(x) x.dbtgt_package{1,1}.ber, a); -% +% % if precomp % legndname1 = ['Pre-Emphasis']; % else % legndname1 = ['No Pre-Emphasis']; % end -% -% +% +% % baudrate = floor( uloops.bitrate.*1e-9 ./log2(m) ) .* 2.5 .* 1e9; -% +% % subplot(1,3,1) % hold on % title(sprintf('%d km | %d nm | PAM %d',uloops.link_length,wavelength,m)); @@ -93,9 +96,9 @@ ylabel('BER'); % beautifyBERplot() % xlabel('Bit Rate in Gbps'); % ylabel('BER'); -% -% -% +% +% +% % subplot(1,3,2) % hold on % % title(sprintf('%d km | %d nm | PAM %d',uloops.link_length,wavelength,m)); @@ -111,8 +114,8 @@ ylabel('BER'); % beautifyBERplot() % xlabel('Bit Rate in Gbps'); % ylabel('BER'); -% -% +% +% % subplot(1,3,3) % hold on % % title(sprintf('%d km | %d nm | PAM %d',uloops.link_length,wavelength,m)); @@ -183,44 +186,44 @@ ylabel('BER'); % end -% +% % m = 6; % ir = [2,2.5,3]; % cols = linspecer(6); % baudrate_gather = []; % for i = 1:3 % m = uloops.M(i); -% -% %%% GET VNLE VALS +% +% %%% GET VNLE VALS % precode = 0; % precomp = 1; % a = wh_master.getStoValue('ber',precomp, precode, uloops.bitrate , uloops.laser_wavelength, m, uloops.link_length); % ber_vnle(i,:) = cellfun(@(x) x.vnle_pf_package{1,1}.ber_vnle, a); -% -% %%% GET DB VALS +% +% %%% GET DB VALS % precode = 1; % precomp = 0; % a = wh_master.getStoValue('ber',precomp, precode, uloops.bitrate , uloops.laser_wavelength, m, uloops.link_length); % ber_db(i,:) = cellfun(@(x) x.dbtgt_package{1,1}.ber, a); -% -% %%% GET MLSE VALS +% +% %%% GET MLSE VALS % precode = 0; % precomp = 0; % a = wh_master.getStoValue('ber',precomp, precode, uloops.bitrate , uloops.laser_wavelength, m, uloops.link_length); % ber_mlse(i,:) = cellfun(@(x) x.vnle_pf_package{1,1}.ber_mlse, a); -% +% % inf_rate_pam(i,:) = cellfun(@(x) x.vnle_pf_package{1,1}.air, a); % inf_rate_pam(i,:) = inf_rate_pam(i,:)./log2(m); -% +% % bitrate = floor( uloops.bitrate.*1e-9 ./log2(m) ) .* ir(i) .* 1e9; % baudrate = floor( uloops.bitrate.*1e-9 ./log2(m) ) .* 1e9; % baudrate_gather = union(baudrate_gather,baudrate); % baudrate_ticks = 100:20:240; % bitrate_ticks = 300:30:480; -% +% % tp = TransmissionPerformance; % netRatesVNLE = tp.calculateNetRate(bitrate, 'NGMI', inf_rate_pam(i,:), 'BER', ber_vnle(i,:)); -% +% % %%% NGMI % figure(12) % hold on @@ -231,8 +234,8 @@ ylabel('BER'); % beautifyBERplot() % xticks(baudrate_ticks); % xlim([min(baudrate_ticks) max(baudrate_ticks)]); -% -% +% +% % %%% AIR % figure(14) % hold on @@ -243,7 +246,7 @@ ylabel('BER'); % beautifyBERplot() % xticks(baudrate_ticks); % xlim([min(baudrate_ticks) max(baudrate_ticks)]); -% +% % %%% RATES % figure(16) % hold on @@ -262,7 +265,7 @@ ylabel('BER'); % ylabel('Net Bitrate in Gbps') % beautifyBERplot() % ylim([250 410]) -% +% % %%% CODE OVERHEAD IN % % figure(18) % hold on @@ -275,7 +278,7 @@ ylabel('BER'); % xlabel('Baud rate in GBd'); % ylabel('FEC Overhead in %') % beautifyBERplot() -% +% % %%% CLASSIC BER % figure(22) % subplot(1,4,i) @@ -298,19 +301,19 @@ ylabel('BER'); % subplot(1,4,4) % hold on % if m == 4 -% +% % plot(bitrate*1e-9,ber_db(i,:),'DisplayName',sprintf('Diff. Code + DB tgt.'),'Color',cols(i,:),'LineStyle','--','HandleVisibility','on','Marker','diamond'); -% +% % elseif m == 6 -% +% % plot(bitrate*1e-9,ber_vnle(i,:),'DisplayName',sprintf('VNLE + PF + MLSE'),'Color',cols(i,:),'LineStyle','-','HandleVisibility','on','Marker','o'); % % plot(bitrate*1e-9,ber_mlse(i,:),'DisplayName',sprintf('MLSE',m),'Color',cols(i,:),'LineStyle',':','HandleVisibility',hv,'Marker','square'); -% +% % elseif m ==8 -% +% % plot(bitrate*1e-9,ber_vnle(i,:),'DisplayName',sprintf('Tx precomp + VNLE'),'Color',cols(i,:),'LineStyle','-','HandleVisibility','on','Marker','o'); % % plot(bitrate*1e-9,ber_mlse(i,:),'DisplayName',sprintf('MLSE',m),'Color',cols(i,:),'LineStyle',':','HandleVisibility',hv,'Marker','square'); -% +% % end % yline(4.85e-3,'HandleVisibility','off'); % yline(2e-2,'HandleVisibility','off'); @@ -321,12 +324,12 @@ ylabel('BER'); % % ylabel('BER') % beautifyBERplot() % set(gca, 'YScale', 'log'); -% -% +% +% % end -% -% -% +% +% +% % figure() % title(sprintf('%d km | 1310 nm | PAM %d | VNLE',uloops.link_length,uloops.M)); % hold on diff --git a/projects/IMDD_base_system/imdd_model.m b/projects/IMDD_base_system/imdd_model.m index 0735636..0386187 100644 --- a/projects/IMDD_base_system/imdd_model.m +++ b/projects/IMDD_base_system/imdd_model.m @@ -62,7 +62,7 @@ mu_dfe = 0.0004; dfe_ = sum(dfe_order)>0; -duob_mode = db_mode.no_db; +% duob_mode = db_mode.no_db; %%% change specific parameter if given in varargin % Parse optional input arguments @@ -77,7 +77,7 @@ if ~isempty(varargin) else eval([fields{i}, ' = ', 'var_s.(fields{',num2str(i),'})' , ';']); end - + end else error('Optional variables should be passed as a struct.'); @@ -132,7 +132,7 @@ El_sig = El_sig .* scaling; %%%%% MODULATE E/O CONVERSION %%%%%% [Opt_sig] = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs,"lambda",laser_wavelength,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",random_key+1,"alpha",alpha).process(El_sig); -Opt_sig.spectrum("displayname",'Opt Spectrum','fignum',10,'normalizeTo0dB',1); +% Opt_sig.spectrum("displayname",'Opt Spectrum','fignum',10,'normalizeTo0dB',1); Opt_sig = Fiber("fsimu",Opt_sig.fs,"fiber_length",link_length,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.07).process(Opt_sig); @@ -165,7 +165,7 @@ Scpe_sig = Scpe_sig.resample("fs_out",2*fsym); Scpe_sig.signal = Scpe_sig.signal(1:2*length(Symbols)); %%%%%% Sync Rx signal with reference %%%%%% -[Scpe_sig,~] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym,"debug_plots",1); +[Scpe_sig,~] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym,"debug_plots",0); Scpe_sig = Filter('filtdegree',4,"f_cutoff",Symbols.fs.*0.5,"fs",Scpe_sig.fs,"filterType",filtertypes.gaussian,"active",true).process(Scpe_sig); @@ -180,57 +180,60 @@ if 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",0); - - output.ffe_results = ffe(eq_ffe,M,Scpe_sig,Symbols,Tx_bits, ... - "precode_mode",duob_mode,'showAnalysis',0,"postFFE",[], ... - "eth_style_symbol_mapping",0); - - output.ffe_results.metrics.print -end -if 0 - % -------------------- DFE -------------------- - eq_dfe = EQ("Ne",ffe_order,"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",0); - - output.dfe_results = ffe(eq_dfe,M,Scpe_sig,Symbols,Tx_bits, ... - "precode_mode",duob_mode,'showAnalysis',0,"postFFE",[], ... - "eth_style_symbol_mapping",0); - - output.dfe_results.metrics.print("description",'DFE'); - -end -if 0 - % -------------------- VNLE + MLSE -------------------- - pf_ncoeffs = 1; - ffe_order3 = [50, 5, 5]; - eq_v = EQ("Ne",ffe_order3,"Nb",dfe_order, ... - "training_length",len_tr,"training_loops",5,"dd_loops",5, ... - "K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... - "FFEmu",0,"plotfinal",0,"ideal_dfe",1); - pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); - - mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); - - [output.vnle_results, output.mlse_results] = vnle_postfilter_mlse(eq_v, pf_, mlse_, M, Scpe_sig, Symbols, Tx_bits, ... - "precode_mode", duob_mode, 'showAnalysis', 0, "postFFE", [], "eth_style_symbol_mapping", 0); -end -if 1 - + + % % -------------------- FFE -------------------- + % ffe_order = [50, 0, 0]; + % eq_ffe = EQ("Ne",ffe_order,"Nb",[0,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",0); + % + % output.ffe_results = ffe(eq_ffe,M,Scpe_sig,Symbols,Tx_bits, ... + % "precode_mode",duob_mode,'showAnalysis',0,"postFFE",[], ... + % "eth_style_symbol_mapping",0); + % + % output.ffe_results.metrics.print + + % % -------------------- DFE -------------------- + % eq_dfe = EQ("Ne",ffe_order,"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",0); + % + % output.dfe_results = ffe(eq_dfe,M,Scpe_sig,Symbols,Tx_bits, ... + % "precode_mode",duob_mode,'showAnalysis',0,"postFFE",[], ... + % "eth_style_symbol_mapping",0); + % + % output.dfe_results.metrics.print("description",'DFE'); + + + % % -------------------- VNLE + MLSE -------------------- + % pf_ncoeffs = 1; + % ffe_order3 = [200, 0, 0]; + % eq_v = EQ("Ne",ffe_order3,"Nb",dfe_order, ... + % "training_length",len_tr,"training_loops",5,"dd_loops",5, ... + % "K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... + % "FFEmu",0,"plotfinal",0,"ideal_dfe",1); + % pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); + % + % mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); + % + % [output.vnle_results, output.mlse_results] = vnle_postfilter_mlse(eq_v, pf_, mlse_, M, Scpe_sig, Symbols, Tx_bits, ... + % "precode_mode", duob_mode, 'showAnalysis', 0, "postFFE", [], "eth_style_symbol_mapping", 0); + + % -------------------- DB target -------------------- mlse_db_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",M,'trellis_states',PAMmapper(M,0).levels); - ffe_order = [50, 5, 5]; + ffe_order = [50, 0, 0]; eq_ = EQ("Ne",ffe_order,"Nb",dfe_order,"training_length",len_tr,"training_loops",5,"dd_loops",5, ... "K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1); output.dbt_results = duobinary_target(eq_,mlse_db_, M, Scpe_sig, Symbols, Tx_bits, ... - "precode_mode", duob_mode, 'showAnalysis', 0, "postFFE", []); - + "precode_mode", duob_mode, 'showAnalysis', 0, "postFFE", [], "decoding_mode", decoding_mode); + output.dbt_results.metrics.print("description",'Duobinary'); - - + + disp('- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ') fprintf('\n') -end - + end \ No newline at end of file diff --git a/projects/db_minimal_example.m b/projects/db_minimal_example.m new file mode 100644 index 0000000..efc7816 --- /dev/null +++ b/projects/db_minimal_example.m @@ -0,0 +1,71 @@ +% einstellungen +M = 2; +apply_precode_at_tx = 0; + +% daten erzeugen +bitpattern = []; +s = RandStream('twister','Seed',1); +for i = 1:log2(M) + N = 2^(17-1); %length of prbs + bitpattern(:,i) = randi(s,[0 1], N, 1); +end + +if M == 6 + bitpattern = reshape(bitpattern',[],1); + bitpattern = bitpattern(1:end-mod(length(bitpattern),5)); +end + +tx_bits = Informationsignal(bitpattern); +tx_symbols = PAMmapper(M,0).map(tx_bits); + +if apply_precode_at_tx + + % Precode + tx_symbols_precoded = Duobinary().precode(tx_symbols); + % Entschiedene Symbole codieren: d_DB(n) = d(n) + d(n-1) (im Fall von PAM4 7 level [0 1 2 3 4 5 6]) + symbols_db = Duobinary().encode(tx_symbols_precoded); + % Entschiedene codierte Symbole decodieren: d_dec(n) = d_DB(n) mod4 + rx_symbols = Duobinary().decode(symbols_db); + bits_rx = PAMmapper(M,0).demap(rx_symbols); + + figure(1) + clf + hold on + stairs(tx_symbols.signal(1:100),'DisplayName','Tx Symbols','LineStyle','-','LineWidth',2); + stairs(rx_symbols.signal(1:100),'DisplayName','Rx Symbols','LineStyle',':','LineWidth',1); + legend + ylim([-2 2]) + + + [~,~,ber,~] = calc_ber(tx_bits.signal,bits_rx.signal,"skip_front",0,"skip_end",0,"returnErrorLocation",1); + disp(['BER: ',sprintf('%.1E',ber)]); + assert(ber == 0) +else + + + % Entschiedene Symbole codieren: d_DB(n) = d(n) + d(n-1) (im Fall von PAM4 7 level [0 1 2 3 4 5 6]) + symbols_db = Duobinary().encode(tx_symbols); + + % Entschiedene codierte Symbole decodieren: d_dec(n) = d_DB(n) mod4 + rx_symbols = Duobinary().decode(symbols_db); + bits_rx = PAMmapper(M,0).demap(rx_symbols); + + % ref symbole precoden: + tx_symbols_ref = Duobinary().encode(tx_symbols); + tx_symbols_ref = Duobinary().decode(tx_symbols_ref); + tx_bits_ref = PAMmapper(M,0).demap(tx_symbols_ref); + + % step plot to check visually if symbols overlap + figure(1) + clf + hold on + stairs(tx_symbols_ref.signal(1:100),'DisplayName','Tx Symbols','LineStyle','-','LineWidth',2); + stairs(rx_symbols.signal(1:100),'DisplayName','Rx Symbols','LineStyle',':','LineWidth',1); + legend + ylim([-2 2]) + + + [~,~,ber,~] = calc_ber(tx_bits_ref.signal,bits_rx.signal,"skip_front",1,"skip_end",1,"returnErrorLocation",1); + disp(['BER: ',sprintf('%.1E',ber)]); + assert(ber == 0) +end \ No newline at end of file diff --git a/test/test_db_minimal_example.m b/test/test_db_minimal_example.m new file mode 100644 index 0000000..55ee52a --- /dev/null +++ b/test/test_db_minimal_example.m @@ -0,0 +1,65 @@ + +M = 4; +apply_precode_at_tx = 1; + +bitpattern = []; +s = RandStream('twister','Seed',1); +for i = 1:log2(M) + N = 2^(17-1); %length of prbs + bitpattern(:,i) = randi(s,[0 1], N, 1); +end + +if M == 6 + bitpattern = reshape(bitpattern',[],1); + bitpattern = bitpattern(1:end-mod(length(bitpattern),5)); +end + +bits = Informationsignal(bitpattern); + +symbols = PAMmapper(M,0).map(bits); + +if apply_precode_at_tx + symbols_tx = Duobinary().precode(symbols); +else + symbols_tx = symbols; +end +disp(['Tx Sequenz: -- RMS:',sprintf('%.1f',rms(symbols_tx.signal)),' - - Levels -',num2str(numel(unique(symbols_tx.signal)))]); +unique(symbols_tx.signal) +disp('- - - - - - - - - -'); + +symbols_tx.signal = awgn(symbols_tx.signal,20,"measured",1); +% show2Dconstellation(symbols_tx,symbols_tx,"displayname",'VNLE Out','fignum',2241); + + +if apply_precode_at_tx + % Entschiedene Symbole codieren: d_DB(n) = d(n) + d(n-1) (im Fall von PAM4 7 level [0 1 2 3 4 5 6]) + symbols_db = Duobinary().encode(symbols_tx); + + disp(['DB encoded -- RMS:',sprintf('%.1f',rms(symbols_db.signal)),' - - Levels -',num2str(numel(unique(symbols_db.signal)))]); + unique(symbols_db.signal) + disp('- - - - - - - - - -'); + + % Entschiedene codierte Symbole decodieren: d_dec(n) = d_DB(n) mod4 + symbols_rx = Duobinary().decode(symbols_db); +else + symbols_db = Duobinary().encode(symbols_tx); + symbols_rx = Duobinary().decode(symbols_db); +end + +% Vergleichen von b(n) und d_dec(n) +bits_rx = PAMmapper(M,0).demap(symbols_rx); +disp(['Wieder normal -- RMS:',sprintf('%.1f',rms(symbols_rx.signal)),' - - Levels -',num2str(numel(unique(symbols_rx.signal)))]); +unique(symbols_rx.signal) +disp('- - - - - - - - - -'); + + +[~,~,ber,~] = calc_ber(bits.signal,bits_rx.signal,"skip_front",10,"skip_end",10,"returnErrorLocation",1); + +disp(['BER: ',sprintf('%.1E',ber),' - - PAM-',num2str(M)]); + +figure() +subplot(1,2,1) +histogram(symbols_tx.signal,100,'Normalization','count') + +subplot(1,2,2) +histogram(symbols_db.signal,100,'Normalization','count') \ No newline at end of file