diff --git a/Classes/00_signals/Electricalsignal.m b/Classes/00_signals/Electricalsignal.m index 7da0e69..e9b4bf7 100644 --- a/Classes/00_signals/Electricalsignal.m +++ b/Classes/00_signals/Electricalsignal.m @@ -108,9 +108,6 @@ classdef Electricalsignal < Signal obj.signal = obj.signal*scling; end - - - end end diff --git a/Classes/00_signals/Opticalsignal.m b/Classes/00_signals/Opticalsignal.m index b39511e..214e9e2 100644 --- a/Classes/00_signals/Opticalsignal.m +++ b/Classes/00_signals/Opticalsignal.m @@ -70,8 +70,7 @@ classdef Opticalsignal < Signal s = mean( abs(obj.signal-mean(obj.signal)).^2 ); cspr = 10*log10(c / s); - - + end diff --git a/Classes/00_signals/Signal.m b/Classes/00_signals/Signal.m index a21e71c..5905556 100644 --- a/Classes/00_signals/Signal.m +++ b/Classes/00_signals/Signal.m @@ -98,7 +98,6 @@ classdef Signal end - %% CONVERT TO Opticalsignal function o_sig = Opticalsignal(obj, options) @@ -143,7 +142,7 @@ classdef Signal arguments obj options.fignum = randi(1000) - options.displayname = []; + options.displayname = ''; options.timeframe = 0; options.clear = 0; options.color = []; @@ -273,7 +272,6 @@ classdef Signal SignalCopy = []; ModifierName = class(CallingModifier); - cell = {SignalType , TimeStamp , Length , SignalPower(1) , Nase, SignalCopy, ModifierName, ModifierCopy, Description}; obj.logbook = [obj.logbook; cell]; @@ -363,8 +361,7 @@ classdef Signal options.fft_length = 2^(nextpow2(length(obj.signal))-9); end - - + if options.normalizeToNyquist == 0 [p_lin,f_Hz] = pwelch(obj.signal, hanning(options.fft_length), ... options.fft_length/2, options.fft_length, ... @@ -488,7 +485,6 @@ classdef Signal end - function move_it_spectrum(obj,options) arguments obj @@ -630,7 +626,6 @@ classdef Signal case power_notation.W %pow = pow % Watt - end end @@ -668,7 +663,7 @@ classdef Signal %% PAPR of signal function papr_db = papr_db(obj) - %PAPR The peak-to-average power ratio (PAPR) is the peak amplitude squared (giving the peak power) + % PAPR The peak-to-average power ratio (PAPR) is the peak amplitude squared (giving the peak power) % divided by the RMS value squared (giving the average power).[1] It is the square of the crest factor. % papr = max(abs(timesignal))^2 / rms(timesignal)^2; ODER papr = peak2rms(sig)^2; @@ -722,7 +717,7 @@ classdef Signal end - %% + %% function [obj,S,inverted,sequenceFound,sequenceStarts] = tsynch(obj,options) % time sync and cut arguments @@ -732,14 +727,11 @@ classdef Signal options.debug_plots = 0; end - - S = {}; inverted = -1; sequenceFound = 0; sequenceStarts = []; - %normalize the signal a = obj.normalize("mode","oneone").signal; @@ -765,11 +757,11 @@ classdef Signal pkpos = sort(pkpos); - if mean(w) > 10 || mean(p) > 10 - return - else - sequenceFound = 1; - end + % if mean(w) > 15 || mean(p) > 15 + % return + % else + % sequenceFound = 1; + % end if options.debug_plots figure(121212);clf @@ -777,8 +769,6 @@ classdef Signal findpeaks(abs(co./max(co)),'MinPeakDistance',length(b)/2,'MinPeakHeight',0.2,'NPeaks',maxpeaknum,'SortStr','descend') end - - shifts = lags(pkpos); sequenceStarts = shifts; shifts = shifts(shifts>=0); @@ -831,7 +821,6 @@ classdef Signal end - %% function obj = filter(obj,a,b) @@ -921,6 +910,7 @@ classdef Signal end + %% function eye(obj,fsym,M,options) @@ -981,8 +971,8 @@ classdef Signal maxA = max(sig(100:end-100))*1.3; minA = min(sig(100:end-100))*1.3; - maxA = 0.12; - minA = -0.08; + % maxA = 0.12; + % minA = -0.08; difference= maxA-minA; diff --git a/Classes/01_transmit/Pulseformer.m b/Classes/01_transmit/Pulseformer.m index 4efe336..d137924 100644 --- a/Classes/01_transmit/Pulseformer.m +++ b/Classes/01_transmit/Pulseformer.m @@ -43,7 +43,7 @@ classdef Pulseformer function signalclass_out = process(obj,signalclass_in) % actual processing of the signal (steps 1. - 3.) - signalclass_in.signal = obj.process_(signalclass_in.signal); + signalclass_in.signal = obj.process_(signalclass_in); % append to logbook lbdesc = 'Applied Pulseshaping'; @@ -65,109 +65,171 @@ classdef Pulseformer % Cant be seen from outside! So put all your functions here that can/ % shall not be called from outside - function data_out = process_(obj,data_in) - %METHOD1 Summary of this method goes here - % Detailed explanation goes here - arguments(Input) - obj - data_in - end - - arguments(Output) - data_out - end - - if ~rem(obj.fdac,obj.fsym) - %ist ein Vielfaches - sps = obj.fdac / obj.fsym; - p = sps; - q = 1; - else - %ist kein Vielfaches - p = obj.fsym / gcd(obj.fdac, obj.fsym); %upsampling p->->-> - q = obj.fdac/ gcd(obj.fdac, obj.fsym); %downsampling <-q - sps= q; %sps während dem pulse shaping - end - - if obj.pulse == pulseform.rc - filtertype = 'normal'; - elseif pulseform.rrc - filtertype = 'sqrt'; - end - - %Bau das Filter (hier rc) - racos_len = obj.pulselength*2; - h = rcosdesign(obj.alpha,racos_len,sps,filtertype); - % h = h./ max(h); - - if obj.matched - h = conj(fliplr(h)); - end - - manual_cyclic_convolution = 0; - upfirdn_convolution = 1; - - if manual_cyclic_convolution - - % Apply filter the long way (from move_it) - data_in = data_in'; - blen = length(data_in)*sps; - - % oversample symbol sequence - symbolov=zeros(size(data_in,1),blen); - symbolov(:,1:sps:blen-sps+1)=data_in; - H=fft(h,blen); - - % Convolution of Bit sequence with impulse response - data_out=ifft( fft(symbolov.') .* repmat( H,size(data_in,1),1 ).' ).'; - data_out = circshift(data_out,[0 -(obj.pulselength*sps)]); +function data_out = process_(obj, data_in_signal) + % Extract the incoming sampling rate + % Safety check: If fs is missing (e.g. raw symbols), assume it is fsym + if isprop(data_in_signal, 'fs') && ~isempty(data_in_signal.fs) + f_in = data_in_signal.fs; + else + f_in = obj.fsym; + end - if rem(obj.fdac,obj.fsym) - data_out = data_out(1:q:end); - end + f_out = obj.fdac; - end + % 1. Calculate Resampling Factors (P and Q) + % We need rational approximation: f_out/f_in = p/q + [p, q] = rat(f_out / f_in); - if upfirdn_convolution + % 2. Calculate SPS for the Filter Design + % The filter operates at the INTERMEDIATE rate (f_in * p). + % We need to know how many samples represent one symbol AT THAT RATE. + fs_intermediate = f_in * p; + sps_filter = fs_intermediate / obj.fsym; - %Apply Filter using Matlab build in fctn. + % 3. Filter Design + if obj.pulse == pulseform.rc + filtertype = 'normal'; + elseif obj.pulse == pulseform.rrc % assuming enum logic holds + filtertype = 'sqrt'; + end - data_out_ = upfirdn(data_in,h,p,q); + % Standard RRC Design + % Note: rcosdesign sps must be integer? Usually yes, but for polyphase it can handle it. + % If sps_filter is not integer, rcosdesign might complain. + % For your setup (powers of 2), it will likely be integer. + racos_len = obj.pulselength; % span in symbols + h = rcosdesign(obj.alpha, racos_len, sps_filter, filtertype); - %cut signal, which is longer due to fir filter - st = round(p/q*racos_len/2); %we need to cut y_out - en = round(st + (length(data_in)*p/q)); - data_out = data_out_(st:en); + % Matched Filter Flip (Complex Conjugate Time Reversal) + if obj.matched + h = conj(fliplr(h)); + end - end + % 4. Processing + % Apply upfirdn using the calculated P and Q + data_out_ = upfirdn(data_in_signal.signal, h, p, q); - if upfirdn_convolution && manual_cyclic_convolution - figure() - subplot(2,1,1) - title("Convolution vs. Upfirdn and Cut") - hold on - % plot(data_out_(1:200),'DisplayName','Matlab upfirdn'); - plot(data_out(1:200),'DisplayName','By Hand cyclic convolution') - subplot(2,1,2) - hold on - plot(data_out(1:2000),'DisplayName','OUT'); - plot(data_in(1:2000),'DisplayName','IN'); - end + % 5. Trim Tail (Group Delay Correction) + % The delay of linear phase filter is (N-1)/2 samples @ intermediate rate + delay_samples_intermediate = (length(h) - 1) / 2; + + % Convert delay to output samples + delay_samples_out = delay_samples_intermediate / q; + + % We usually want to trim the "start" transient + st = floor(delay_samples_out) + 1; + + % Calculate expected output length + len_out = ceil(length(data_in_signal.signal) * p / q); + + % Cut + data_out = data_out_(st : st + len_out - 1); +end - %scaling?! see pulsef module line 696 -% scale = max(max([abs(real(data_out)) abs(imag(data_out))])); %find max value from real and imag part -% data_out = data_out./scale; - - % data_out = data_out'; - - %Check output integrity - if abs(round(p/q * length(data_in)) - length(data_out)) > 4 - warning('Check signal length after pulse shaping'); - %disp('Check signal length after pulse shaping'); - end - - - end +% +% function data_out = process_(obj,data_in) +% %METHOD1 Summary of this method goes here +% % Detailed explanation goes here +% arguments(Input) +% obj +% data_in +% end +% +% arguments(Output) +% data_out +% end +% +% if ~rem(obj.fdac,obj.fsym) +% %ist ein Vielfaches +% sps = obj.fdac / obj.fsym; +% p = sps; +% q = 1; +% else +% %ist kein Vielfaches +% p = obj.fsym / gcd(obj.fdac, obj.fsym); %upsampling p->->-> +% q = obj.fdac/ gcd(obj.fdac, obj.fsym); %downsampling <-q +% sps= q; %sps während dem pulse shaping +% end +% +% if obj.pulse == pulseform.rc +% filtertype = 'normal'; +% elseif pulseform.rrc +% filtertype = 'sqrt'; +% end +% +% %Bau das Filter (hier rc) +% racos_len = obj.pulselength*2; +% h = rcosdesign(obj.alpha,racos_len,sps,filtertype); +% % h = h./ max(h); +% +% if obj.matched +% h = conj(fliplr(h)); +% end +% +% manual_cyclic_convolution = 0; +% upfirdn_convolution = 1; +% +% if manual_cyclic_convolution +% +% % Apply filter the long way (from move_it) +% data_in = data_in'; +% blen = length(data_in)*sps; +% +% % oversample symbol sequence +% symbolov=zeros(size(data_in,1),blen); +% symbolov(:,1:sps:blen-sps+1)=data_in; +% H=fft(h,blen); +% +% % Convolution of Bit sequence with impulse response +% data_out=ifft( fft(symbolov.') .* repmat( H,size(data_in,1),1 ).' ).'; +% data_out = circshift(data_out,[0 -(obj.pulselength*sps)]); +% +% if rem(obj.fdac,obj.fsym) +% data_out = data_out(1:q:end); +% end +% +% end +% +% if upfirdn_convolution +% +% %Apply Filter using Matlab build in fctn. +% +% data_out_ = upfirdn(data_in,h,p,q); +% +% %cut signal, which is longer due to fir filter +% st = round(p/q*racos_len/2); %we need to cut y_out +% en = round(st + (length(data_in)*p/q)); +% data_out = data_out_(st:en); +% +% end +% +% if upfirdn_convolution && manual_cyclic_convolution +% figure() +% subplot(2,1,1) +% title("Convolution vs. Upfirdn and Cut") +% hold on +% % plot(data_out_(1:200),'DisplayName','Matlab upfirdn'); +% plot(data_out(1:200),'DisplayName','By Hand cyclic convolution') +% subplot(2,1,2) +% hold on +% plot(data_out(1:2000),'DisplayName','OUT'); +% plot(data_in(1:2000),'DisplayName','IN'); +% end +% +% %scaling?! see pulsef module line 696 +% % scale = max(max([abs(real(data_out)) abs(imag(data_out))])); %find max value from real and imag part +% % data_out = data_out./scale; +% +% % data_out = data_out'; +% +% %Check output integrity +% if abs(round(p/q * length(data_in)) - length(data_out)) > 4 +% warning('Check signal length after pulse shaping'); +% %disp('Check signal length after pulse shaping'); +% end +% +% +% end end diff --git a/Classes/02_optical/EML.m b/Classes/02_optical/EML.m index 32a902c..4a97686 100644 --- a/Classes/02_optical/EML.m +++ b/Classes/02_optical/EML.m @@ -15,6 +15,7 @@ classdef EML u_pi randomkey randomstream + alpha %on instance creation field @@ -38,6 +39,7 @@ classdef EML options.lambda; options.power; options.linewidth = 0; + options.alpha = 0; options.ampl_imbal = 0; options.pha_imbal = 0; options.bias; @@ -101,9 +103,21 @@ classdef EML end %modulate the laserfield with the electrical signal - opt_out = obj.externalmodulation(laserfield,elec_in); + laserfield = obj.externalmodulation(laserfield,elec_in); + % add chirp + opt_out = obj.chirp(laserfield); + end + + function chirped_field = chirp(obj,laserfield) + + % Chirp + p = abs(laserfield.^2); + derv_p = [0; diff(p)]; + delta_phi = derv_p./(4*pi*p).*obj.alpha; + delta_phi = cumsum(delta_phi); + chirped_field = laserfield.*exp(1i*2*pi*delta_phi); end diff --git a/Classes/04_DSP/Timing_Recovery.m b/Classes/04_DSP/Timing_Recovery.m new file mode 100644 index 0000000..ba16f10 --- /dev/null +++ b/Classes/04_DSP/Timing_Recovery.m @@ -0,0 +1,47 @@ +classdef Timing_Recovery < handle + + properties(Access=public) + timing_error_detector + sps + damping_factor + normalized_loop_bandwidth + detector_gain + end + + methods(Access=public) + function obj = FFE(options) + arguments(Input) + + options.timing_error_detector = 'Gardner'; + options.sps = 2; + options.damping_factor = 1.0; + options.normalized_loop_bandwidth = 0.005; + options.detector_gain = 1; + + end + + fn = fieldnames(options); + for n = 1:numel(fn) + obj.(fn{n}) = options.(fn{n}); + end + + obj.e = zeros(obj.order,1); + obj.error = 0; + + end + + function data_out = process(obj, data_in) + + timing_synchronization = comm.SymbolSynchronizer( ... + "TimingErrorDetector", obj.timing_error_detector, ... + "SamplesPerSymbol", obj.sps, ... + "DampingFactor", obj.damping_factor, ... + "NormalizedLoopBandwidth", obj.normalized_loop_bandwidth, ... + "DetectorGain", obj.detector_gain); + + data_out.signal = timing_synchronization(data_in.signal); + + end + end +end + diff --git a/Classes/DataBaseHandler/Metricstruct.m b/Classes/DataBaseHandler/Metricstruct.m index ef00675..fe9a818 100644 --- a/Classes/DataBaseHandler/Metricstruct.m +++ b/Classes/DataBaseHandler/Metricstruct.m @@ -8,9 +8,9 @@ classdef Metricstruct date_of_processing (1,1) datetime = datetime('now') numBits (1,1) double {mustBeInteger, mustBeNonnegative} = 0 - BER (1,1) double {mustBeNumeric, mustBeNonnegative, mustBeLessThanOrEqual(BER,1)} = 0 + BER (1,1) double {mustBeNumeric, mustBeNonnegative} = 0 numBitErr (1,1) double {mustBeInteger, mustBeNonnegative} = 0 - BER_precoded (1,1) double {mustBeNumeric, mustBeNonnegative, mustBeLessThanOrEqual(BER_precoded,1)} = 0 + BER_precoded (1,1) double {mustBeNumeric, mustBeNonnegative} = 0 numBitErr_precoded (1,1) double {mustBeInteger, mustBeNonnegative} = 0 SNR (1,1) double {mustBeNumeric} = NaN diff --git a/Classes/Moveit_wrapper.m b/Classes/Moveit_wrapper.m index 621ab67..8c7ea8a 100644 --- a/Classes/Moveit_wrapper.m +++ b/Classes/Moveit_wrapper.m @@ -27,7 +27,7 @@ classdef Moveit_wrapper < handle % Ensure the function exists if ~exist(obj.moveit_function_name, 'file') - error('Function "%s" does not exist.', obj.moveit_function_name); + error('Function "%s" does not exist. The move-it module must be on path for Matlab, otherwise the wrapper can not call it...', obj.moveit_function_name); end % Step 1: Get default parameters by calling moveit module diff --git a/Functions/EQ_structures/ffe.m b/Functions/EQ_structures/ffe.m index 6cfefc5..df4138e 100644 --- a/Functions/EQ_structures/ffe.m +++ b/Functions/EQ_structures/ffe.m @@ -131,11 +131,11 @@ switch precode_mode 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", 30000, "skip_end", 150, "returnErrorLocation", 1); + [~, 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); - [bits, errors, ber, error_pos] = calc_ber(rx_bits.signal, tx_bits.signal, "skip_front", 30000, "skip_end", 150, "returnErrorLocation", 1); + [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 @@ -143,12 +143,12 @@ switch precode_mode 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", 30000, "skip_end", 150, "returnErrorLocation", 1); + [~, 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", 30000, "skip_end", 150, "returnErrorLocation", 1); + [bits, errors, ber, error_pos] = calc_ber(rx_bits.signal, tx_bits_demapped.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1); end end diff --git a/Functions/Metrics/air_garcia_implementation.m b/Functions/Metrics/air_garcia_implementation.m index c71a26a..c560bf5 100644 --- a/Functions/Metrics/air_garcia_implementation.m +++ b/Functions/Metrics/air_garcia_implementation.m @@ -11,7 +11,6 @@ if nargin == 4 M_training = []; end - % if input is complex, separate into real and imaginary parts if any(imag(x(:))~=0) || any(imag(r(:))~=0) x = [real(x); imag(x)]; diff --git a/projects/FSO_transmission/first_analysis.m b/projects/FSO_transmission/first_analysis.m new file mode 100644 index 0000000..c47a2b1 --- /dev/null +++ b/projects/FSO_transmission/first_analysis.m @@ -0,0 +1,171 @@ + +base = "C:\Users\Silas\Nextcloud\Dokumente\02_Ablage_Office\FSO_FP_QCL_60umUTC"; +mode = 0; %0 oder 1 +M = 2; + +all_files = dir(fullfile(base, "**/*.mat")); + +if M == 2 + tx_data = load("C:\Users\Silas\Nextcloud\Dokumente\02_Ablage_Office\FSO_FP_QCL_60umUTC\14G_PAM2\tx_info\tx_info_PAM2_14Gbd0.75RRC.mat"); + filename = fullfile(base, "14G_PAM2\M=2_Rs=1.4e10_Fs=8e10_I=265mA_RoP=46.3mW_L=31m_PS=RRC_rolloff=0.75_Mode=Rise.mat"); +elseif M == 4 + tx_data = load("C:\Users\Silas\Nextcloud\Dokumente\02_Ablage_Office\FSO_FP_QCL_60umUTC\6G_PAM4\tx_info\tx_info_PAM4_6Gbd0.6RRC.mat"); + filename = fullfile(base, "6G_PAM4\M=4_Rs=6e9_Fs=8e10_I=255mA_RoP=42.3mW_L=31m_PS=RRC_rolloff=0.6_Mode=Rise.mat"); +end + +if mode == 1 + [f, p] = uigetfile(fullfile(base, "**/*.mat")); + if f~=0 + filename = fullfile(p,f); + end +end + + +datas = load(filename); + +%% +str = filename; +M_ = str2double(regexp(str, 'M=([^_]+)', 'tokens', 'once')); +assert(M==M_); +fsym = str2double(regexp(str, 'Rs=([^_]+)', 'tokens', 'once')); +fs = str2double(regexp(str, 'Fs=([^_]+)', 'tokens', 'once')); +I = sscanf(char(regexp(str, 'I=([^_]+)', 'tokens', 'once')), '%f'); +rop = sscanf(char(regexp(str, 'RoP=([^_]+)', 'tokens', 'once')), '%f'); +L = sscanf(char(regexp(str, 'L=([^_]+)', 'tokens', 'once')), '%f'); +pulseshape = string( regexp(str, 'PS=([^_]+)', 'tokens', 'once')); +rolloff = str2double(regexp(str, 'rolloff=([^_]+)', 'tokens', 'once')); +mode = string( regexp(str, 'Mode=([^\.]+)', 'tokens', 'once')); + +%% +% Tx data + +Bits = Informationsignal(tx_data.tx_data,"fs",fsym); +Symbols = Informationsignal(real(tx_data.tx_PAM_sym),"fs",fsym); + +mapping_style = M==4; % Pam2 is like move-it; PAM-4 is different, same mapping like ETH peopled used in Zurich... hence the "eth_style" argument here and there +PM = PAMmapper(M,0,"eth_style",mapping_style); % one should rename "eth style" as this is simply a different mapping scheme + +Symbols_ = PM.map(Bits) .* PM.scaling; +assert(isequal(Symbols.signal,Symbols_.signal)); + +Bits_ = PM.demap(Symbols); +[bits,errors,ber,errorIndice] = calc_ber(Bits_.signal,Bits.signal); +assert(ber == 0); + +%% For comparison, apply pulsef on Tx Symbols +Pform = Pulseformer("fsym",fsym,"fdac",fs,"pulse","rc","pulselength",16,"alpha",rolloff); +Digi_sig_compare = Pform.process(Symbols); +MF = Pulseformer("fsym",fsym,"fdac",2*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff); +Rx_sig_compare = MF.process(Digi_sig_compare); + +%% + +% Rx Data +traceData = datas.tr.lastData(2).trace.ch3; + +%FYI: Voltage=(RawData−YReference)×YIncrement+YOrigin +scoperead_volts = (traceData.RawData - traceData.YReference) * traceData.YIncrement + traceData.YOrigin; +demystified = isequal(traceData.YData,scoperead_volts); +assert(demystified); + +Scope_sig = Electricalsignal(traceData.YData,"fs",fs); + +Scope_sig.plot("displayname",'raw','fignum',100); +Scope_sig.spectrum("displayname",'raw','fignum',101) + +% 1) matched filter +% pulse is symmetric, hence we can use pulsef firectly as matched filter. +% It feels off (bit I think correct) that the fsym is now the output freq.!! +% -> output 2 sps to omit timing recovery!? +Pform = Pulseformer("fsym",fsym,"fdac",2*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff,"matched",1); +Rx_matched = Pform.process(Scope_sig); +Rx_matched.spectrum("displayname",'Signal after matched filter','fignum',1); + + +%% + +sys = comm.SymbolSynchronizer('TimingErrorDetector', 'Gardner (non-data-aided)', ... + 'SamplesPerSymbol', 2, ... + 'DampingFactor', 0.7, ... + 'NormalizedLoopBandwidth', 0.01); +Rx_symbolsync = Rx_matched; +[Rx_symbolsync.signal, timing_error] = sys(Rx_matched.signal); + +plot(timing_error); % If this is a ramp, you have drift! + +%% timing sync -> at this point we still have no symbol timing recovery, we +% % try to do this with 2sps EQ! + +[~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_symbolsync.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1); + + +%% not working.. +Rx_synced = Rx_synced_cell{1}; +len_tr = 4096*2; +mu_ffe1 = 0.0001; +mu_ffe2 = 0.0008; +mu_ffe3 = 0.001; +mu_dc = 0.005; +mu_ffe = [mu_ffe1 mu_ffe3 mu_ffe3]; +mu_dfe = 0.0004; +duob_mode = db_mode.no_db; + +Rx_synced.plot("displayname",'RX: Matched+Sync+2sps','fignum',2); + +Digi_sig_compare.normalize("mode","rms").spectrum("displayname",'Tx: RC-shaped','fignum',1,'normalizeTo0dB',0); +Rx_sig_compare.normalize("mode","rms").spectrum("displayname",'Tx: RC-shaped + matched filtered ','fignum',1,'normalizeTo0dB',0); +Rx_synced.normalize("mode","rms").spectrum("displayname",'RX: matched filtered + synced','fignum',1,'normalizeTo0dB',0); + +if M == 2 + ber_in_paper = 10^(-2.6); %fig 3a) 4 Gb/s MWIR FSO Transmission using Directly Modulated QCL and an Uncooled UTC-PD at Room-Temperature +elseif M == 4 + ber_in_paper = 10^(-2.5); +end + +%% -------------------- FFE -------------------- +% requires some more digging what is going on :-) +eq_ffe = EQ("Ne",[50, 5, 5],"Nb",[2,0,0], ... + "training_length",len_tr,"training_loops",5,"dd_loops",5, ... + "K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... + "FFEmu",0,"plotfinal",0,"ideal_dfe",1); + +ffe_results = ffe(eq_ffe,M,Rx_synced,Symbols,Bits, ... + "precode_mode",duob_mode,'showAnalysis',0,"postFFE",[], ... + "eth_style_symbol_mapping",mapping_style); + +% ffe_results.metrics.print +fprintf('My EQ: %.1e \n',ffe_results.metrics.BER); +fprintf('Paper: %.1e \n \n',ber_in_paper); + +%% -------------------- VNLE + MLSE -------------------- + +pf_ncoeffs = 1; +eq_v = EQ("Ne",[100, 5, 5],"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",1); +pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); +mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0,"eth_style",mapping_style).levels); + +[vnle_results, mlse_results] = vnle_postfilter_mlse(eq_v, pf_, mlse_, M, Rx_synced, Symbols, Bits, ... + "precode_mode", duob_mode, 'showAnalysis', 1, "postFFE", [], "eth_style_symbol_mapping", mapping_style); + +mlse_results.metrics.print("description",'MLSE') +fprintf('My EQ: %.1e \n',mlse_results.metrics.BER); +fprintf('Paper: %.1e \n \n',ber_in_paper); + +%% -------------------- DB target -------------------- +mlse_db_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",M,'trellis_states',PAMmapper(M,0).levels); + +eq_ = EQ("Ne",[50, 5, 5],"Nb",[0,0,0],"training_length",len_tr,"training_loops",5,"dd_loops",5, ... + "K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1); + +dbt_results = duobinary_target(eq_,mlse_db_, M, Rx_synced, Symbols, Bits, ... + "precode_mode", duob_mode, 'showAnalysis', 0, "postFFE", [],"eth_style_symbol_mapping",mapping_style); + +dbt_results.metrics.print("description",'Duobinary'); +mlse_results.metrics.print +fprintf('My EQ: %.1e \n',dbt_results.metrics.BER); +fprintf('Paper: %.1e \n \n',ber_in_paper); + + diff --git a/projects/IMDD_base_system/imdd_it.m b/projects/IMDD_base_system/imdd_it.m index bd2461d..a495d74 100644 --- a/projects/IMDD_base_system/imdd_it.m +++ b/projects/IMDD_base_system/imdd_it.m @@ -1,117 +1,136 @@ -% basePath = 'C:\Users\Silas\Documents\MATLAB\Datensätze\sioe_labor\'; -% db = DBHandler("pathToDB",[basePath,'silas_labor.db']); + if 1 + uloops = struct; uloops.precomp = [1]; - uloops.db_precode = [0]; - uloops.bitrate = [224].*1e9; %[300,330,360,390,420,450,480] [224,336,360,390,420,448] for MPI + uloops.bitrate = [300].*1e9; %[300,330,360,390,420,450,480] [224,336,360,390,420,448] for MPI % uloops.laser_wavelength = [1293,1297.5,1302,1306.5,1310,1313.4,1318,1322.7,1327.4]; - uloops.laser_wavelength = [1310]; + uloops.laser_wavelength = [1293]; uloops.M = [4]; - uloops.link_length = [1]; % 1,2,3,5,6,8,10 - uloops.interference_attenuation = [0,3,6,9,12,15,18,21,24,27,30,45]; + uloops.link_length = [0:2:10]; % 1,2,3,5,6,8,10 + uloops.alpha = [0]; + wh = DataStorage(uloops); wh.addStorage("ber"); - % wh = submit_simulations(wh,"parallel",0,"simulation_mode",0); wh = submit_handle(@imdd_model,wh,"parallel",0); end -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); - ber_vnle_pc = cellfun(@(x) x.vnle_pf_package{1,1}.ber_vnle, a); - %MLSE - 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); - %MLSE - 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)); - title(sprintf('PAM %d',m)); - plot(baudrate.*1e-9,ber_vnle,'DisplayName',['Pre-Emphasis: ', num2str(precomp), '| Diff.-Code: 0'],'Color',cols(1+precomp,:),'LineStyle','-','HandleVisibility','on'); - plot(baudrate.*1e-9,ber_vnle_pc,'DisplayName',['Pre-Emphasis: ', num2str(precomp), '| Diff.-Code: 1'],'Color',cols(1+precomp,:),'LineStyle','-.','HandleVisibility','on'); - xticks(baudrate.*1e-9); - set(gca, 'YScale', 'log'); - ylim([5e-5 0.4]); - xlim([min(baudrate(2).*1e-9), max(baudrate.*1e-9) ]); - yline([3.8e-3, 2e-2],'HandleVisibility','off'); - legend - 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)); - title(sprintf('PAM %d',m)); - plot(baudrate.*1e-9,ber_dbtgt,'DisplayName',['Pre-Emphasis: ', num2str(precomp), '| Diff.-Code: 0'],'Color',cols(5+precomp,:),'LineStyle','-','HandleVisibility','on'); - plot(baudrate.*1e-9,ber_dbtgt_pc,'DisplayName',['Pre-Emphasis: ', num2str(precomp), '| Diff.-Code: 1'],'Color',cols(5+precomp,:),'LineStyle','-.','HandleVisibility','on'); - xticks(baudrate.*1e-9); - set(gca, 'YScale', 'log'); - ylim([5e-5 0.4]); - xlim([min(baudrate(2).*1e-9), max(baudrate.*1e-9) ]); - yline([3.8e-3, 2e-2],'HandleVisibility','off'); - legend - 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)); - title(sprintf('PAM %d',m)); - plot(baudrate.*1e-9,ber_mlse,'DisplayName',['Pre-Emphasis: ', num2str(precomp), '| Diff.-Code: 0'],'Color',cols(3+precomp,:),'LineStyle','-','HandleVisibility','on'); - plot(baudrate.*1e-9,ber_mlse_pc,'DisplayName',['Pre-Emphasis: ', num2str(precomp), '| Diff.-Code: 1'],'Color',cols(3+precomp,:),'LineStyle','-.','HandleVisibility','on'); - xticks(baudrate.*1e-9); - set(gca, 'YScale', 'log'); - ylim([5e-5 0.4]); - xlim([min(baudrate(2).*1e-9), max(baudrate.*1e-9) ]); - yline([3.8e-3, 2e-2],'HandleVisibility','off'); - legend - beautifyBERplot() - xlabel('Bit Rate in Gbps'); - ylabel('BER'); - 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); + plot(uloops.link_length,ffe,'DisplayName',sprintf('Alpha: %d',alpha),'LineStyle','-','HandleVisibility','on'); end -% -% + +set(gca, 'YScale', 'log'); +ylim([5e-5 0.4]); +yline([3.8e-3, 2e-2],'HandleVisibility','off'); +legend +beautifyBERplot() +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); +% ber_vnle_pc = cellfun(@(x) x.vnle_pf_package{1,1}.ber_vnle, a); +% %MLSE +% 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); +% %MLSE +% 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)); +% title(sprintf('PAM %d',m)); +% plot(baudrate.*1e-9,ber_vnle,'DisplayName',['Pre-Emphasis: ', num2str(precomp), '| Diff.-Code: 0'],'Color',cols(1+precomp,:),'LineStyle','-','HandleVisibility','on'); +% plot(baudrate.*1e-9,ber_vnle_pc,'DisplayName',['Pre-Emphasis: ', num2str(precomp), '| Diff.-Code: 1'],'Color',cols(1+precomp,:),'LineStyle','-.','HandleVisibility','on'); +% xticks(baudrate.*1e-9); +% set(gca, 'YScale', 'log'); +% ylim([5e-5 0.4]); +% xlim([min(baudrate(2).*1e-9), max(baudrate.*1e-9) ]); +% yline([3.8e-3, 2e-2],'HandleVisibility','off'); +% legend +% 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)); +% title(sprintf('PAM %d',m)); +% plot(baudrate.*1e-9,ber_dbtgt,'DisplayName',['Pre-Emphasis: ', num2str(precomp), '| Diff.-Code: 0'],'Color',cols(5+precomp,:),'LineStyle','-','HandleVisibility','on'); +% plot(baudrate.*1e-9,ber_dbtgt_pc,'DisplayName',['Pre-Emphasis: ', num2str(precomp), '| Diff.-Code: 1'],'Color',cols(5+precomp,:),'LineStyle','-.','HandleVisibility','on'); +% xticks(baudrate.*1e-9); +% set(gca, 'YScale', 'log'); +% ylim([5e-5 0.4]); +% xlim([min(baudrate(2).*1e-9), max(baudrate.*1e-9) ]); +% yline([3.8e-3, 2e-2],'HandleVisibility','off'); +% legend +% 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)); +% title(sprintf('PAM %d',m)); +% plot(baudrate.*1e-9,ber_mlse,'DisplayName',['Pre-Emphasis: ', num2str(precomp), '| Diff.-Code: 0'],'Color',cols(3+precomp,:),'LineStyle','-','HandleVisibility','on'); +% plot(baudrate.*1e-9,ber_mlse_pc,'DisplayName',['Pre-Emphasis: ', num2str(precomp), '| Diff.-Code: 1'],'Color',cols(3+precomp,:),'LineStyle','-.','HandleVisibility','on'); +% xticks(baudrate.*1e-9); +% set(gca, 'YScale', 'log'); +% ylim([5e-5 0.4]); +% xlim([min(baudrate(2).*1e-9), max(baudrate.*1e-9) ]); +% yline([3.8e-3, 2e-2],'HandleVisibility','off'); +% legend +% beautifyBERplot() +% xlabel('Bit Rate in Gbps'); +% ylabel('BER'); +% end +% end +% % +% % % cols = linspecer(7);%cbrewer2('Set2',10); % % @@ -163,163 +182,163 @@ end % end - -m = 6; -ir = [2,2.5,3]; -cols = linspecer(6); -baudrate_gather = []; -for i = 1:3 - m = uloops.M(i); - - %%% 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 - 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 - 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 - title(sprintf('Performance at 1310 nm')); - plot(baudrate.*1e-9,inf_rate_pam(i,:),'DisplayName',sprintf('NGMI; PAM %d',m),'Color',cols(i,:),'LineStyle','-'); - xlabel('Baud rate in GBd'); - ylabel('NGMI') - beautifyBERplot() - xticks(baudrate_ticks); - xlim([min(baudrate_ticks) max(baudrate_ticks)]); - - - %%% AIR - figure(14) - hold on - title(sprintf('Performance at 1310 nm')); - plot(baudrate.*1e-9,inf_rate_pam(i,:).*bitrate.*1e-9,'DisplayName',sprintf('AIR; PAM %d',m),'Color',cols(i,:),'LineStyle','-'); - xlabel('Baud rate in GBd'); - ylabel('AIR'); - beautifyBERplot() - xticks(baudrate_ticks); - xlim([min(baudrate_ticks) max(baudrate_ticks)]); - - %%% RATES - figure(16) - hold on - title(sprintf('Performance at 1310 nm')); - if i == 1 - hv = 'on'; - else - hv = 'off'; - end - plot(baudrate*1e-9,netRatesVNLE.SDHD.NetRate.*1e-9,'DisplayName',sprintf('SD+HD FEC',m),'Color',cols(i,:),'LineStyle','-','HandleVisibility',hv,'Marker','o'); - plot(baudrate.*1e-9,netRatesVNLE.HD.NetRate.*1e-9,'DisplayName',sprintf('HD FEC',m),'Color',cols(i,:),'LineStyle',':','HandleVisibility',hv,'Marker','diamond'); - plot(baudrate.*1e-9,netRatesVNLE.KP4_hamming.NetRate*1e-9,'DisplayName',sprintf('KP4+Hamming',m),'Color',cols(i,:),'LineStyle','-.','HandleVisibility',hv,'Marker','square'); - xticks(baudrate_ticks); - xlim([min(baudrate_ticks) max(baudrate_ticks)]); - xlabel('Baud rate in GBd'); - ylabel('Net Bitrate in Gbps') - beautifyBERplot() - ylim([250 410]) - - %%% CODE OVERHEAD IN % - figure(18) - hold on - title(sprintf('Performance at 1310 nm')); - plot(baudrate*1e-9,100.*(1-netRatesVNLE.SDHD.CodeRate)./netRatesVNLE.SDHD.CodeRate,'DisplayName',sprintf('SD+HD FEC',m),'Color',cols(i,:),'LineStyle','-','HandleVisibility',hv,'Marker','o'); - plot(baudrate.*1e-9,100.*(1-netRatesVNLE.HD.CodeRate)./netRatesVNLE.HD.CodeRate,'DisplayName',sprintf('HD FEC',m),'Color',cols(i,:),'LineStyle',':','HandleVisibility',hv,'Marker','diamond'); - plot(baudrate.*1e-9,100.*(1-netRatesVNLE.KP4_hamming.CodeRate)./netRatesVNLE.KP4_hamming.CodeRate,'DisplayName',sprintf('KP4+Hamming',m),'Color',cols(i,:),'LineStyle','-.','HandleVisibility',hv,'Marker','square'); - xticks(baudrate_ticks); - xlim([min(baudrate_ticks) max(baudrate_ticks)]); - xlabel('Baud rate in GBd'); - ylabel('FEC Overhead in %') - beautifyBERplot() - - %%% CLASSIC BER - figure(22) - subplot(1,4,i) - hold on - plot(baudrate*1e-9,ber_vnle(i,:),'DisplayName',sprintf('Tx precomp + VNLE',m),'Color',cols(i,:),'LineStyle','-','HandleVisibility','on','Marker','o'); - plot(baudrate*1e-9,ber_mlse(i,:),'DisplayName',sprintf('VNLE + PF + MLSE',m),'Color',cols(i,:),'LineStyle',':','HandleVisibility','on','Marker','square'); - plot(baudrate*1e-9,ber_db(i,:),'DisplayName',sprintf('Diff. Code + DB tgt.',m),'Color',cols(i,:),'LineStyle','--','HandleVisibility','on','Marker','diamond'); - yline(4.85e-3,'HandleVisibility','off'); - yline(2e-2,'HandleVisibility','off'); - xticks(baudrate*1e-9); - xlim([min(baudrate*1e-9) max(baudrate*1e-9)]); - ylim([1e-4 0.3]); - xlabel('Baudrate in GBd'); - if i == 1 - ylabel('BER') - end - beautifyBERplot() - set(gca, 'YScale', 'log'); -legend - 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'); - xticks(bitrate_ticks); - xlim([min(bitrate_ticks) max(bitrate_ticks)]); - ylim([1e-4 0.3]); - xlabel('Gross Bitrate in Gbps'); - % 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 -line([min(uloops.bitrate.*1e-9) max(uloops.bitrate.*1e-9)],[min(uloops.bitrate.*1e-9) max(uloops.bitrate.*1e-9)],'Color',[.7,.7,.7],'Marker','none','Handlevisibility','off'); -plot(uloops.bitrate.*1e-9,cellfun(@min, inf_rate_vnle).*uloops.bitrate./log2(uloops.M).*1e-9,'DisplayName',sprintf('AIR'),'Color',cols(1,:),'LineStyle',':'); -plot(uloops.bitrate.*1e-9,netRatesVNLE.SDHD.NetRate.*1e-9,'DisplayName',sprintf('SD+HD'),'Color',cols(2,:),'LineStyle',':'); -plot(uloops.bitrate.*1e-9,netRatesVNLE.HD.NetRate.*1e-9,'DisplayName',sprintf('HD'),'Color',cols(3,:),'LineStyle',':'); -plot(uloops.bitrate.*1e-9,netRatesVNLE.KP4_hamming.NetRate.*1e-9,'DisplayName',sprintf('KP4+Hamming'),'Color',cols(4,:),'LineStyle',':'); -beautifyBERplot() -xlim([min(uloops.bitrate.*1e-9) max(uloops.bitrate.*1e-9)]) -xlabel('Gross Bitrate in Gbps'); -ylabel('Net Bitrate in Gbps'); -legend +% +% m = 6; +% ir = [2,2.5,3]; +% cols = linspecer(6); +% baudrate_gather = []; +% for i = 1:3 +% m = uloops.M(i); +% +% %%% 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 +% 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 +% 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 +% title(sprintf('Performance at 1310 nm')); +% plot(baudrate.*1e-9,inf_rate_pam(i,:),'DisplayName',sprintf('NGMI; PAM %d',m),'Color',cols(i,:),'LineStyle','-'); +% xlabel('Baud rate in GBd'); +% ylabel('NGMI') +% beautifyBERplot() +% xticks(baudrate_ticks); +% xlim([min(baudrate_ticks) max(baudrate_ticks)]); +% +% +% %%% AIR +% figure(14) +% hold on +% title(sprintf('Performance at 1310 nm')); +% plot(baudrate.*1e-9,inf_rate_pam(i,:).*bitrate.*1e-9,'DisplayName',sprintf('AIR; PAM %d',m),'Color',cols(i,:),'LineStyle','-'); +% xlabel('Baud rate in GBd'); +% ylabel('AIR'); +% beautifyBERplot() +% xticks(baudrate_ticks); +% xlim([min(baudrate_ticks) max(baudrate_ticks)]); +% +% %%% RATES +% figure(16) +% hold on +% title(sprintf('Performance at 1310 nm')); +% if i == 1 +% hv = 'on'; +% else +% hv = 'off'; +% end +% plot(baudrate*1e-9,netRatesVNLE.SDHD.NetRate.*1e-9,'DisplayName',sprintf('SD+HD FEC',m),'Color',cols(i,:),'LineStyle','-','HandleVisibility',hv,'Marker','o'); +% plot(baudrate.*1e-9,netRatesVNLE.HD.NetRate.*1e-9,'DisplayName',sprintf('HD FEC',m),'Color',cols(i,:),'LineStyle',':','HandleVisibility',hv,'Marker','diamond'); +% plot(baudrate.*1e-9,netRatesVNLE.KP4_hamming.NetRate*1e-9,'DisplayName',sprintf('KP4+Hamming',m),'Color',cols(i,:),'LineStyle','-.','HandleVisibility',hv,'Marker','square'); +% xticks(baudrate_ticks); +% xlim([min(baudrate_ticks) max(baudrate_ticks)]); +% xlabel('Baud rate in GBd'); +% ylabel('Net Bitrate in Gbps') +% beautifyBERplot() +% ylim([250 410]) +% +% %%% CODE OVERHEAD IN % +% figure(18) +% hold on +% title(sprintf('Performance at 1310 nm')); +% plot(baudrate*1e-9,100.*(1-netRatesVNLE.SDHD.CodeRate)./netRatesVNLE.SDHD.CodeRate,'DisplayName',sprintf('SD+HD FEC',m),'Color',cols(i,:),'LineStyle','-','HandleVisibility',hv,'Marker','o'); +% plot(baudrate.*1e-9,100.*(1-netRatesVNLE.HD.CodeRate)./netRatesVNLE.HD.CodeRate,'DisplayName',sprintf('HD FEC',m),'Color',cols(i,:),'LineStyle',':','HandleVisibility',hv,'Marker','diamond'); +% plot(baudrate.*1e-9,100.*(1-netRatesVNLE.KP4_hamming.CodeRate)./netRatesVNLE.KP4_hamming.CodeRate,'DisplayName',sprintf('KP4+Hamming',m),'Color',cols(i,:),'LineStyle','-.','HandleVisibility',hv,'Marker','square'); +% xticks(baudrate_ticks); +% xlim([min(baudrate_ticks) max(baudrate_ticks)]); +% xlabel('Baud rate in GBd'); +% ylabel('FEC Overhead in %') +% beautifyBERplot() +% +% %%% CLASSIC BER +% figure(22) +% subplot(1,4,i) +% hold on +% plot(baudrate*1e-9,ber_vnle(i,:),'DisplayName',sprintf('Tx precomp + VNLE',m),'Color',cols(i,:),'LineStyle','-','HandleVisibility','on','Marker','o'); +% plot(baudrate*1e-9,ber_mlse(i,:),'DisplayName',sprintf('VNLE + PF + MLSE',m),'Color',cols(i,:),'LineStyle',':','HandleVisibility','on','Marker','square'); +% plot(baudrate*1e-9,ber_db(i,:),'DisplayName',sprintf('Diff. Code + DB tgt.',m),'Color',cols(i,:),'LineStyle','--','HandleVisibility','on','Marker','diamond'); +% yline(4.85e-3,'HandleVisibility','off'); +% yline(2e-2,'HandleVisibility','off'); +% xticks(baudrate*1e-9); +% xlim([min(baudrate*1e-9) max(baudrate*1e-9)]); +% ylim([1e-4 0.3]); +% xlabel('Baudrate in GBd'); +% if i == 1 +% ylabel('BER') +% end +% beautifyBERplot() +% set(gca, 'YScale', 'log'); +% legend +% 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'); +% xticks(bitrate_ticks); +% xlim([min(bitrate_ticks) max(bitrate_ticks)]); +% ylim([1e-4 0.3]); +% xlabel('Gross Bitrate in Gbps'); +% % 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 +% line([min(uloops.bitrate.*1e-9) max(uloops.bitrate.*1e-9)],[min(uloops.bitrate.*1e-9) max(uloops.bitrate.*1e-9)],'Color',[.7,.7,.7],'Marker','none','Handlevisibility','off'); +% plot(uloops.bitrate.*1e-9,cellfun(@min, inf_rate_vnle).*uloops.bitrate./log2(uloops.M).*1e-9,'DisplayName',sprintf('AIR'),'Color',cols(1,:),'LineStyle',':'); +% plot(uloops.bitrate.*1e-9,netRatesVNLE.SDHD.NetRate.*1e-9,'DisplayName',sprintf('SD+HD'),'Color',cols(2,:),'LineStyle',':'); +% plot(uloops.bitrate.*1e-9,netRatesVNLE.HD.NetRate.*1e-9,'DisplayName',sprintf('HD'),'Color',cols(3,:),'LineStyle',':'); +% plot(uloops.bitrate.*1e-9,netRatesVNLE.KP4_hamming.NetRate.*1e-9,'DisplayName',sprintf('KP4+Hamming'),'Color',cols(4,:),'LineStyle',':'); +% beautifyBERplot() +% xlim([min(uloops.bitrate.*1e-9) max(uloops.bitrate.*1e-9)]) +% xlabel('Gross Bitrate in Gbps'); +% ylabel('Net Bitrate in Gbps'); +% legend % plot(uloops.bitrate.*1e-9,ber_vnle,'DisplayName',sprintf('NGMI MLSE; %d km',len),'Color',cols(1,:),'LineStyle','-'); diff --git a/projects/IMDD_base_system/imdd_model.m b/projects/IMDD_base_system/imdd_model.m index d400350..5a89483 100644 --- a/projects/IMDD_base_system/imdd_model.m +++ b/projects/IMDD_base_system/imdd_model.m @@ -19,19 +19,13 @@ fdac = 256e9; fadc = 256e9; random_key = 1; -interference_attenuation = 0; -is_mpi = 1; - -precomp = 0; -db_precode = 0; - -db_encode = 0; - rcalpha = 0.05; kover = 16; + vbias_rel = 0.5; -u_pi = 2.9; +u_pi = 3; vbias = -vbias_rel*u_pi; + laser_wavelength = 1293; laser_linewidth = 0; tx_bw_nyquist = 0.8; @@ -40,7 +34,7 @@ tx_bw_nyquist = 0.8; link_length = 1; % RX -rop = -5; +rop = -8; rx_bw_nyquist = 0.8; vnle_order1 = 50; @@ -68,7 +62,7 @@ mu_dfe = 0.0004; dfe_ = sum(dfe_order)>0; -doub_mode = db_mode.no_db; +duob_mode = db_mode.no_db; %%% change specific parameter if given in varargin % Parse optional input arguments @@ -90,43 +84,6 @@ if ~isempty(varargin) end end -if doub_mode ~= db_mode.db_encoded - if precomp == 0 && db_precode == 1 - doub_mode = db_mode.db_precoded; - - db_precode = 1; % preceded data (in my measurement set, this corresponds to low precomp too!) - discard_precode = 0; % - emulate_precode = 0; - legendentry = 'low precomp; precoded'; - disp('low precomp; precoded') - elseif precomp == 1 && db_precode == 1 - doub_mode = db_mode.db_emulate; - - db_precode = 0; % preceded data (in my measurement set, this corresponds to low precomp too!) - discard_precode = 0; % - emulate_precode = 1; - legendentry = 'high precomp; precoded'; - disp('high precomp; precoded') - elseif precomp == 0 && db_precode == 0 - doub_mode = db_mode.db_discard; - - db_precode = 1; % preceded data (in my measurement set, this corresponds to low precomp too!) - discard_precode = 1; % - emulate_precode = 0; - legendentry = 'no precomp; not precoded'; - disp('no precomp; not precoded') - elseif precomp == 1 && db_precode == 0 - doub_mode = db_mode.no_db; - - db_precode = 0; % preceded data (in my measurement set, this corresponds to low precomp too!) - discard_precode = 0; % - emulate_precode = 0; - legendentry = 'high precomp; not precoded'; - disp('high precomp; not precoded') - end -else - -end fsym_ = floor( bitrate*1e-9./log2(M) ).*1e9; @@ -135,239 +92,141 @@ if fsym_ ~= fsym % fprintf('Adapted symbolrate to %d GBd, to match provided bitrate of %d GBit/s using PAM %d \n',fsym.*1e-9,bitrate.*1e-9, M); end + + + + f_nyquist = fsym/2; -%%% run the simulation or measurement or ... -if simulation_mode +Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",16,"alpha",rcalpha); - Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",16,"alpha",rcalpha); - rcalpha = 1; - Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"alpha",rcalpha); +[Digi_sig,Symbols,Tx_bits] = PAMsource(... + "fsym",fsym,"M",M,"order",18,"useprbs",0,... + "fs_out",fdac,... + "applyclipping",0,"clipfactor",1.5,... + "applypulseform",apply_pulsef,"pulseformer",Pform,... + "randkey",random_key,... + 'duobinary_mode',duob_mode,... + "mrds_code",0,"mrds_blocklength",512).process(); - db_precode = 0; - db_encode = 0; - apply_pulsef = 1; - [Digi_sig,Symbols,Tx_bits] = PAMsource(... - "fsym",fsym,"M",M,"order",18,"useprbs",0,... - "fs_out",fdac,... - "applyclipping",0,"clipfactor",1.5,... - "applypulseform",apply_pulsef,"pulseformer",Pform,... - "randkey",random_key,... - "db_precode",db_precode,"db_encode",db_encode,... - "mrds_code",0,"mrds_blocklength",512).process(); +Digi_sig.spectrum("displayname",'Digi Spectrum','fignum',10,'normalizeTo0dB',1); - Digi_sig.spectrum("displayname",'Digi Spectrum','fignum',10,'normalizeTo0dB',1); +%%%%% AWG +% El_sig = M8199A("kover",kover).process(Digi_sig); +El_sig = AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",0,"kover",kover,"bit_resolution",12,"upsampling_method","samplehold","precomp_sinc_rolloff",1).process(Digi_sig); +% El_sig.spectrum("displayname",'Digi Spectrum','fignum',100,'normalizeTo0dB',0); +% El_sig = El_sig.setPower(0,"dBm"); +%%%%% Low-pass el. components %%%%%% +tx_bwl = tx_bw_nyquist.*f_nyquist; +% tx_bwl = 80e9; +El_sig = Filter('filtdegree',4,"f_cutoff",tx_bwl,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(El_sig); +% El_sig.spectrum("displayname",'Digi Spectrum','fignum',100,'normalizeTo0dB',1); +%%%%% Electrical Driver Amplifier %%%%%% +% El_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",3).process(El_sig); +El_sig = El_sig.normalize("mode","oneone"); +scaling = 0.6*(u_pi/2-abs(vbias-u_pi/2)); +El_sig = El_sig .* scaling; - %%%%% AWG - % El_sig = M8199A("kover",kover).process(Digi_sig); - El_sig = AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",0,"kover",kover,"bit_resolution",12,"upsampling_method","samplehold","precomp_sinc_rolloff",1).process(Digi_sig); - % El_sig.spectrum("displayname",'Digi Spectrum','fignum',100,'normalizeTo0dB',0); - % El_sig = El_sig.setPower(0,"dBm"); +%%%%% 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); - %%%%% Low-pass el. components %%%%%% - tx_bwl = tx_bw_nyquist.*f_nyquist; - % tx_bwl = 80e9; - El_sig = Filter('filtdegree',4,"f_cutoff",tx_bwl,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(El_sig); - % El_sig.spectrum("displayname",'Digi Spectrum','fignum',100,'normalizeTo0dB',1); +Opt_sig.spectrum("displayname",'Opt Spectrum','fignum',10,'normalizeTo0dB',1); - %%%%% Electrical Driver Amplifier %%%%%% - El_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",3).process(El_sig); - El_sig = El_sig.normalize("mode","oneone"); +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); - %%%%% MODULATE E/O CONVERSION %%%%%% - [Opt_sig] = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs,"lambda",laser_wavelength,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",random_key+1).process(El_sig); +%%%%%% ROP %%%%%% +Rx_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",rop).process(Opt_sig); - Opt_sig = Fiber("fsimu",Opt_sig.fs,"fiber_length",link_length/1000,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.07).process(Opt_sig); +%%%%%% PD Square Law %%%%%% +Rx_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20,"nep",1.8e-11).process(Rx_sig); - %%%%%% ROP %%%%%% - Rx_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",rop).process(Opt_sig); +%%%%%% Low-pass RX (PD, El. Connectors and Scope %%%%%% +rx_bwl = rx_bw_nyquist.*f_nyquist; +% rx_bwl = 80e9; +Rx_sig = Filter('filtdegree',4,"f_cutoff",rx_bwl,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(Rx_sig); - %%%%%% PD Square Law %%%%%% - Rx_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20,"nep",1.8e-11).process(Rx_sig); +% %%%%%% Low-pass Scope %%%%%% +Lp_scpe = Filter('filtdegree',4,"f_cutoff",110e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true); - %%%%%% Low-pass RX (PD, El. Connectors and Scope %%%%%% - rx_bwl = rx_bw_nyquist.*f_nyquist; - % rx_bwl = 80e9; - Rx_sig = Filter('filtdegree',4,"f_cutoff",rx_bwl,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(Rx_sig); +% Rx_sig.spectrum("displayname",'Analog Rx Spectrum','fignum',100,'normalizeTo0dB',1); - % %%%%%% Low-pass Scope %%%%%% - Lp_scpe = Filter('filtdegree',4,"f_cutoff",110e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true); - - % Rx_sig.spectrum("displayname",'Analog Rx Spectrum','fignum',100,'normalizeTo0dB',1); - - %%%%%% Scope %%%%%% - Scpe_sig = Scope("fsimu",fdac*kover,"fadc",fadc,... - "delay",0,"fixed_delay",0,"filtertype",filtertypes.butterworth,... - "samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,... - "adcresolution",8,"quantbuffer",0.1,'block_dc',1,'lpf_active',1,'H_lpf',Lp_scpe).process(Rx_sig); - - Scpe_cell{1} = Scpe_sig; - -else - profile on - basePath = 'C:\Users\Silas\Documents\MATLAB\Datensätze\sioe_labor\'; - database = DBHandler("pathToDB",[basePath,'silas_labor.db']); - profile off - basePath = 'C:\Users\Silas\Documents\MATLAB\Datensätze\sioe_labor\'; - useGui = 0; - % db = DBHandler("pathToDB",[basePath,'silas_labor.db']); - filterParams = database.tables; - % filterParams.Runs.run_id = 2958; % no db - % filterParams.Runs.run_id = 2937; % no db - filterParams.Configurations = struct( ... - 'bitrate', bitrate, ... - 'db_mode', db_precode+db_encode, ... - 'fiber_length', link_length, ... - 'interference_attenuation', [], ... - 'interference_path_length', [], ... - 'is_mpi', is_mpi, ... - 'pam_level', M, ... - 'precomp_amp', [], ... - 'rop_attenuation', 0, ... - 'symbolrate', [], ... - 'v_awg', [], ... - 'v_bias', [], ... - 'wavelength', laser_wavelength ... - ); - - selectedFields = {'Runs.run_id','Runs.tx_bits_path', 'Runs.tx_symbols_path', 'Runs.rx_sync_path','Runs.rx_raw_path',... - 'Configurations.db_mode','Configurations.pam_level','Configurations.bitrate','Configurations.symbolrate','Configurations.fiber_length','Configurations.wavelength','Configurations.precomp_amp','Measurements.power_rop','Configurations.v_bias',... - 'Configurations.interference_attenuation'}; - - [dataTable,sql_query] = database.queryDB(filterParams, selectedFields); - [~, uniqueIdx] = unique(dataTable.run_id); % Get unique run_id indices - dataTable = dataTable(uniqueIdx,:); % Extract unique configurations for each run_id - fprintf('Found %d entries for requested Configuration. IDs are: %s \n \n',size(dataTable,1),jsonencode(dataTable.run_id(1:min(size(dataTable,1),100)))); - - Tx_bits = load([basePath, char(dataTable.tx_bits_path(end))]); - Tx_bits = Tx_bits.Bits; - - Symbols = load([basePath, char(dataTable.tx_symbols_path(end))]); - Symbols = Symbols.Symbols; - - Scpe_load = load([basePath, char(dataTable.rx_sync_path(end))]); - Scpe_cell = Scpe_load.S; - - - % Raw_signal = load([basePath, char(dataTable.rx_raw_path(1))]); - % Raw_signal.Scpe_sig_raw.plot("displayname",'0db atten','fignum',10101) - % Raw_signal = Raw_signal.Scpe_sig_raw; - % - % Raw_signal = Filter('filtdegree',4,"f_cutoff",Symbols.fs.*0.55,"fs",Raw_signal.fs,"filterType",filtertypes.gaussian,"active",true).process(Raw_signal); - % - % Scpe_cell{1}.eye(fsym,M,"displayname",'eye','fignum',227); - % - % Raw_signal.spectrum("normalizeTo0dB",0,"fignum",11,"fft_length",2^12); - % Raw_signal.move_it_spectrum("fignum",334); - % Raw_signal.move_it_spectrum("fignum",334); - - fsym = Symbols.fs; - -end - -if db_precode - Symbols_precoded = Symbols; -end +%%%%%% Scope %%%%%% +Scpe_sig = Scope("fsimu",fdac*kover,"fadc",fadc,... + "delay",0,"fixed_delay",0,"filtertype",filtertypes.butterworth,... + "samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,... + "adcresolution",8,"quantbuffer",0.1,'block_dc',1,'lpf_active',1,'H_lpf',Lp_scpe).process(Rx_sig); output = struct(); -vnle_pf_package = {}; -vnle_dfe_package = {}; -dbtgt_package = {}; + +%%%%%% Sample to 2x fsym %%%%%% +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 = Filter('filtdegree',4,"f_cutoff",Symbols.fs.*0.5,"fs",Scpe_sig.fs,"filterType",filtertypes.gaussian,"active",true).process(Scpe_sig); + +Scpe_sig = Scpe_sig - mean(Scpe_sig.signal); + +%%% EQUALIZING -proc_occ = min(1,length(Scpe_cell)); -for occ = 1%:proc_occ +% -------------------- 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); - Scpe_sig = Scpe_cell{occ}; +output.ffe_results = ffe(eq_ffe,M,Scpe_sig,Symbols,Tx_bits, ... + "precode_mode",duob_mode,'showAnalysis',0,"postFFE",[], ... + "eth_style_symbol_mapping",0); - %%%%%% Sample to 2x fsym %%%%%% - Scpe_sig = Scpe_sig.resample("fs_out",2*fsym); +output.ffe_results.metrics.print - %%%%%% Sync Rx signal with reference %%%%%% - [Scpe_sig,~] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym); +% -------------------- 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); - Scpe_sig = Filter('filtdegree',4,"f_cutoff",Symbols.fs.*0.5,"fs",Scpe_sig.fs,"filterType",filtertypes.gaussian,"active",true).process(Scpe_sig); +output.dfe_results = ffe(eq_dfe,M,Scpe_sig,Symbols,Tx_bits, ... + "precode_mode",duob_mode,'showAnalysis',0,"postFFE",[], ... + "eth_style_symbol_mapping",0); - Scpe_sig = Scpe_sig - mean(Scpe_sig.signal); - % - % Pform = Pulseformer("fsym",Scpe_sig.fs,"fdac",2*fsym,"pulse","rrc","pulselength",16,"alpha",rcalpha,"matched",0); - % - % Scpe_sig_matched = Pform.process(Scpe_sig); - % - % Scpe_sig.spectrum("normalizeTo0dB",0,"fignum",336,"displayname","scope "); - % Scpe_sig_matched.spectrum("normalizeTo0dB",0,"fignum",336,"displayname","matched"); - - %%% EQUALIZING - - - % eq_mlse = FFE_DCremoval("epochs_tr",5,"epochs_dd",5,"len_tr",len_tr,"mu_dd",mu_ffe(1),"mu_tr",0,"order",ffe_order(1),"sps",2,"decide",0,"dc_buffer_len",1,"mu_dc",0.05); - % eq_mlse = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",len_tr,"mu_dd",mu_ffe(1),"mu_tr",0,"order",ffe_order(1),"sps",2,"decide",0); - % eq_mlse = FFE_DCremoval("epochs_tr",5,"epochs_dd",5,"len_tr",len_tr,"mu_dd",mu_ffe(1),"mu_tr",0,"order",ffe_order(1),"sps",2,"decide",0,"dc_buffer_len",512,"mu_dc",0.05); - - mu_ffe = [mu_ffe1 mu_ffe2 mu_ffe3]; - vnle_order=[vnle_order1,vnle_order2,vnle_order3]; - - % %%%%% VNLE + DFE %%%% - if 0 - - eq_vnle_dfe = EQ("Ne",vnle_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); - eq_2 = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",2001,"sps",1,"decide",0); - - [result] = vnle(eq_vnle_dfe,M,Scpe_sig,Symbols,Tx_bits,"precode_mode",doub_mode,"showAnalysis",1,"postFFE",[]); - vnle_dfe_package{occ} = result; - - end - %%%%% VNLE + PF + MLSE %%%% - if 1 - - % len_tr = length(Symbols)-1000; - eq_vnle_ = EQ("Ne",vnle_order,"Nb",dfe_order,"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1); - % eq_vnle_ = VNLE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",[0.0004 0.0005 0.0006],"mu_tr",0,"order",vnle_order,"sps",2,"decide",0); - pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); - mlse_ = MLSE_viterbi("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); - - [result] = vnle_postfilter_mlse(eq_vnle_,pf_,mlse_,M,Scpe_sig,Symbols,Tx_bits,"precode_mode",doub_mode,'showAnalysis',1); - vnle_pf_package{occ} = result; - - end +output.dfe_results.metrics.print("description",'DFE'); - %%%%% Duobinary Targeting %%%% - if 1 +% -------------------- 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_db = MLSE_viterbi("DIR",[1,1],"duobinary_output",0,"M",M,"trellis_states",PAMmapper(M,0).levels); - eq_db = EQ("Ne",vnle_order,"Nb",dfe_order,"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1); +mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); - [result] = duobinary_target(eq_db, mlse_db, M, Scpe_sig, Symbols, Tx_bits, "precode_mode", doub_mode,'showAnalysis',0); - dbtgt_package{occ} = result; - - - end - - %%%%%% %db signaling => db encoded %%%%% - if 0 - mlse_db_enc = MLSE_viterbi("DIR",[1,1],"duobinary_output",0,"M",M,"trellis_states",PAMmapper(M,0).levels); - eq_db_enc = EQ("Ne",vnle_order,"Nb",dfe_order,"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1); - [result] = duobinary_signaling(eq_db_enc, mlse_db_enc,M, Scpe_sig ,Symbols, Tx_bits); - dbenc_package{occ} = result; - end +[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); - % autoArrangeFigures; - disp('- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ') - fprintf('\n') - +% -------------------- DB target -------------------- +mlse_db_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",M,'trellis_states',PAMmapper(M,0).levels); +ffe_order = [50, 5, 5]; +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", []); -end +output.dbt_results.metrics.print("description",'Duobinary'); -output.vnle_dfe_package = vnle_dfe_package; -output.vnle_pf_package = vnle_pf_package; -output.dbtgt_package = dbtgt_package; -if ~isempty(curFolder) - cd(curFolder); -end +disp('- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ') +fprintf('\n') end \ No newline at end of file diff --git a/projects/IMDD_base_system/minimal_example.m b/projects/IMDD_base_system/minimal_example.m new file mode 100644 index 0000000..89f24ff --- /dev/null +++ b/projects/IMDD_base_system/minimal_example.m @@ -0,0 +1,156 @@ +% minimal example IM/DD + +M = 4; +fsym = 180e9; + +apply_pulsef = 1; +fdac = 256e9; +fadc = 256e9; +random_key = 1; + +rcalpha = 0.05; +kover = 16; + +duob_mode = db_mode.no_db; + +vbias_rel = 0.5; +u_pi = 3; +vbias = -vbias_rel*u_pi; + +laser_wavelength = 1293; +laser_linewidth = 0; +tx_bw_nyquist = 0.8; + +% Channel +link_length = 1; + +% RX +rop = -8; +rx_bw_nyquist = 0.8; + +vnle_order1 = 50; +vnle_order2 = 7; +vnle_order3 = 7; + +vnle_order=[vnle_order1,vnle_order2,vnle_order3]; +dfe_order = [0 0 0]; + +pf_ncoeffs = 1; + +alpha = 0; + +len_tr = 4096*2; + +mu_ffe1 = 0.0001; +mu_ffe2 = 0.0008; +mu_ffe3 = 0.001; +mu_dc = 0.005; +% mu_dc = 0; + +mu_ffe = [mu_ffe1 mu_ffe3 mu_ffe3]; +mu_dfe = 0.0004; + + +Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",16,"alpha",rcalpha); + +[Digi_sig,Symbols,Tx_bits] = PAMsource(... + "fsym",fsym,"M",M,"order",18,"useprbs",0,... + "fs_out",fdac,... + "applyclipping",0,"clipfactor",1.5,... + "applypulseform",apply_pulsef,"pulseformer",Pform,... + "randkey",random_key,... + 'duobinary_mode',duob_mode,... + "mrds_code",0,"mrds_blocklength",512).process(); + +%%%%% AWG +El_sig = M8199A("kover",kover).process(Digi_sig); +% El_sig = AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",0,"kover",kover,"bit_resolution",12,"upsampling_method","samplehold","precomp_sinc_rolloff",1).process(Digi_sig); +El_sig.spectrum("displayname",'Digi Spectrum','fignum',100,'normalizeTo0dB',0); +% El_sig = El_sig.setPower(0,"dBm"); + +%%%%% Electrical Driver Amplifier %%%%%% +% El_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",3).process(El_sig); +El_sig = El_sig.normalize("mode","oneone"); +scaling = 0.6*(u_pi/2-abs(vbias-u_pi/2)); +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.eye(fsym,M,"displayname",'eye adter modulator','fignum',2026); + +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); + +%%%%%% ROP %%%%%% +Rx_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",rop).process(Opt_sig); + +%%%%%% PD Square Law %%%%%% +Rx_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20,"nep",1.8e-11).process(Rx_sig); + +%%%%%% Low-pass RX (PD, El. Connectors and Scope %%%%%% +rx_bwl = 80e9; +Rx_sig = Filter('filtdegree',4,"f_cutoff",rx_bwl,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(Rx_sig); + +% %%%%%% Low-pass Scope %%%%%% +Lp_scpe = Filter('filtdegree',4,"f_cutoff",110e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true); + +%%%%%% Scope %%%%%% +Scpe_sig = Scope("fsimu",fdac*kover,"fadc",fadc,... + "delay",0,"fixed_delay",0,"filtertype",filtertypes.butterworth,... + "samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,... + "adcresolution",8,"quantbuffer",0.1,'block_dc',1,'lpf_active',1,'H_lpf',Lp_scpe).process(Rx_sig); + +%%%%%% Sample to 2x fsym %%%%%% +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",0); +Scpe_sig.spectrum("displayname",'Opt Spectrum','fignum',11,'normalizeTo0dB',1); + +Scpe_sig = Filter('filtdegree',4,"f_cutoff",Symbols.fs.*0.5,"fs",Scpe_sig.fs,"filterType",filtertypes.gaussian,"active",true).process(Scpe_sig); + +Scpe_sig = Scpe_sig - mean(Scpe_sig.signal); + +% -------------------- 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 + +% -------------------- 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); + + +% -------------------- DB target -------------------- +mlse_db_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",M,'trellis_states',PAMmapper(M,0).levels); +ffe_order = [50, 5, 5]; +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", []); + +output.dbt_results.metrics.print("description",'Duobinary'); + + + diff --git a/projects/IMDD_base_system/simulation_bwl_2.m b/projects/IMDD_base_system/simulation_bwl_2.m index dec8036..e336ad2 100644 --- a/projects/IMDD_base_system/simulation_bwl_2.m +++ b/projects/IMDD_base_system/simulation_bwl_2.m @@ -59,7 +59,7 @@ Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",1 db_precode = 0; db_encode = 0; -duob_mode = db_mode.db_precoded; +duob_mode = db_mode.no_db; apply_pulsef = 1; [Digi_sig,Symbols,Tx_bits] = PAMsource(... "fsym",fsym,"M",M,"order",18,"useprbs",0,... @@ -72,13 +72,12 @@ apply_pulsef = 1; Digi_sig.spectrum("displayname",'Digi Spectrum','fignum',10,'normalizeTo0dB',1); -%% proof of concept -Symbols_db = Duobinary().encode(Symbols); -mim_decoded = Duobinary().decode(Symbols_db,"M",M); -rx_bits_mim_decoded = PAMmapper(M,0,"eth_style",0).demap(mim_decoded); -rx_bits_mim_decoded_.signal = circshift(rx_bits_mim_decoded.signal,0); -[~,~,ber_mim_decode,~] = calc_ber(rx_bits_mim_decoded_.signal,Tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); -fprintf('BER mim: %.2e \n',ber_mim_decode); +% %% proof of concept memoryless inverse mapping (direct db targeting and decoding) +% Symbols_db = Duobinary().encode(Symbols); +% mim_decoded = Duobinary().decode(Symbols_db,"M",M); +% rx_bits_mim_decoded = PAMmapper(M,0,"eth_style",0).demap(mim_decoded); +% [~,~,ber_mim_decode,~] = calc_ber(rx_bits_mim_decoded_.signal,Tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); +% fprintf('BER mim: %.2e \n',ber_mim_decode); %% diff --git a/projects/Messung_Zürich/submit_handle.m b/projects/IMDD_base_system/submit_handle.m similarity index 100% rename from projects/Messung_Zürich/submit_handle.m rename to projects/IMDD_base_system/submit_handle.m diff --git a/projects/WDM/WDM_auswertung.m b/projects/WDM/WDM_auswertung.m index 1f3b949..98b42b9 100644 --- a/projects/WDM/WDM_auswertung.m +++ b/projects/WDM/WDM_auswertung.m @@ -1,248 +1,63 @@ +base = "C:\Users\Silas\Nextcloud\Cluster"; +all_files = dir(fullfile(base, "**/*.mat")); + +schemes = ["co","pair","alt","seg"]; + +% Preallocate as table (minimal + convenient) +T = table('Size',[0 10], ... + 'VariableTypes', ["string","string","string","datetime","double","double","double","double","double","double"], ... + 'VariableNames', ["folder","file","scheme","date","node","jobid","L_km","Nch","df_GHz","alpha"]); + +% filename parser +rx = "^WDM_(?\d{8})_(?