From ed179534074f0dbf2e63081c101ef1ab4e59cf5f Mon Sep 17 00:00:00 2001 From: Silas Oettinghaus Date: Fri, 19 Apr 2024 10:31:36 +0200 Subject: [PATCH] Update April --- .DS_Store | Bin 6148 -> 0 bytes Classes/00_signals/Electricalsignal.m | 24 +- Classes/00_signals/Informationsignal.m | 17 +- Classes/00_signals/Opticalsignal.m | 46 +- Classes/00_signals/Signal.m | 107 ++- Classes/01_transmit/AWG.m | 33 +- Classes/01_transmit/M8196A.m | 35 + Classes/01_transmit/M8199B.m | 34 + Classes/01_transmit/PAMmapper.m | 51 +- Classes/01_transmit/Pulseformer.m | 20 +- Classes/02_etc/Filter.m | 195 +++-- Classes/02_optical/EML.m | 9 +- Classes/02_optical/Fiber.m | 2 +- Classes/03_receive/Photodiode.m | 28 +- Classes/03_receive/Scope.m | 27 +- Classes/04_DSP/EQ.m | 9 +- Classes/04_DSP/EQ_silas.m | 95 ++- Classes/04_DSP/EQ_silas_plain.m | 764 +++++++++-------- .../EQ_silas_sliding_window_dc_removal.m | 408 +++++++++ Functions/calc_ber.m | 57 +- MPI_offline_analysis/mpi_analysis_pam4.m | 120 --- .../IMDD_base_system/filter_debug.m | 140 +-- .../IMDD_base_system/imddmodel.m | 451 +++++----- .../IMDD_base_system/run_loop.m | 416 ++++----- .../IMDD_base_system/setup_simulation.m | 798 +++++++++--------- projects/MPI_April/mpi_simulation_cspr.m | 198 +++++ projects/MPI_April/plot_analysis_window.m | 70 ++ projects/MPI_April/simulate_tx_sig.m | 55 ++ .../MPI_offline_analysis/mpi_analysis_pam4.m | 444 ++++++++++ .../mpi_analysis_pam4_old_notworking.m | 188 +++++ 30 files changed, 3261 insertions(+), 1580 deletions(-) delete mode 100644 .DS_Store create mode 100644 Classes/01_transmit/M8196A.m create mode 100644 Classes/01_transmit/M8199B.m create mode 100644 Classes/04_DSP/EQ_silas_sliding_window_dc_removal.m delete mode 100644 MPI_offline_analysis/mpi_analysis_pam4.m rename filter_debug.m => projects/IMDD_base_system/filter_debug.m (94%) rename imddmodel.m => projects/IMDD_base_system/imddmodel.m (96%) rename run_loop.m => projects/IMDD_base_system/run_loop.m (96%) rename setup_simulation.m => projects/IMDD_base_system/setup_simulation.m (97%) create mode 100644 projects/MPI_April/mpi_simulation_cspr.m create mode 100644 projects/MPI_April/plot_analysis_window.m create mode 100644 projects/MPI_April/simulate_tx_sig.m create mode 100644 projects/MPI_offline_analysis/mpi_analysis_pam4.m create mode 100644 projects/MPI_offline_analysis/mpi_analysis_pam4_old_notworking.m diff --git a/.DS_Store b/.DS_Store deleted file mode 100644 index 6f25ac248c3bbfe0a1cbc90dc313136dc8fd974d..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 6148 zcmeHKUuzmc5T8|J_K<)NEi`WneHH4(Bo_K|HMSohSK5cRSkDtZa^8w})(4T^dqTfY z-`h_R=(oss$m}1BL4`a>C^7@HzuDQDUH7-gt^)w9H4R;WDgba)g#`zj4Wh=Cu1Udq z%7{$QkwF9!5=evDS~NNSM+Vf{Iq)8)cxpXZIlpuq#4^CqU%@zz(yU(p%L>KPLHUp$ z@(Taz59Q3yy4f^sbtj)y?}Lj(!~iis4D5mdd%f_=E+(Tnhyh|?uNYAGgFsd2SF5_W&a;>sH0i)>@xkoM%vLB&ua5Hz9S+PhNG&ly3_N6@ zXu4JP{=eH^|36Hk9x*@+>=gs7)b`sg+>(7;*EXl#T7mivRf&qr4gO0(LqEkBtDfQ) bs9Ml3l!543%nhOkh5iU=8mJ)#o|J)m!Y6V& diff --git a/Classes/00_signals/Electricalsignal.m b/Classes/00_signals/Electricalsignal.m index 2c067f2..2d28cfe 100644 --- a/Classes/00_signals/Electricalsignal.m +++ b/Classes/00_signals/Electricalsignal.m @@ -11,25 +11,39 @@ classdef Electricalsignal < Signal %ELECTRICALSIGNAL Construct an instance of this class % Detailed explanation goes here arguments - signal + signal options.fs options.logbook end obj = obj@Signal(signal); - + fn = fieldnames(options); for l = 1:numel(fn) - obj.(fn{l}) = options.(fn{l}); + obj.(fn{l}) = options.(fn{l}); end end - function pow = power(obj) - pow = mean(abs(obj.signal),"all") ; + pow = mean( abs(obj.signal).^2 ) ; + pow = pow2db(pow)+30; %dbm + + end + + function cspr = cspr(obj) + + carrier_power_dbm = pow2db( abs(mean(obj.signal)).^2 )+30; % dB -> +30 -> dBm + signal_power_dbm = obj.power; + + carr_lin = abs(mean(obj.signal)).^2; + sign_lin = mean(abs(obj.signal).^2); + + cspr_lin = pow2db(carr_lin/sign_lin); + + cspr =carrier_power_dbm-signal_power_dbm; end diff --git a/Classes/00_signals/Informationsignal.m b/Classes/00_signals/Informationsignal.m index 1e32d05..afa7462 100644 --- a/Classes/00_signals/Informationsignal.m +++ b/Classes/00_signals/Informationsignal.m @@ -3,15 +3,26 @@ classdef Informationsignal < Signal % Detailed explanation goes here properties - - + fs end methods - function obj = Informationsignal(signal) + function obj = Informationsignal(signal,options) %OPTICALSIGNAL Construct an instance of this class % Detailed explanation goes here + arguments + signal + options.fs + options.logbook + end + obj = obj@Signal(signal); + + fn = fieldnames(options); + + for l = 1:numel(fn) + obj.(fn{l}) = options.(fn{l}); + end end diff --git a/Classes/00_signals/Opticalsignal.m b/Classes/00_signals/Opticalsignal.m index d63c452..7afc4b5 100644 --- a/Classes/00_signals/Opticalsignal.m +++ b/Classes/00_signals/Opticalsignal.m @@ -32,11 +32,53 @@ classdef Opticalsignal < Signal end + %% + function [obj,delay_n] = delay(obj,options) + + arguments + obj Opticalsignal + options.delay_meter double = 0 + options.delay_samples = 0 + end + + if options.delay_meter ~= 0 + delay_t = options.delay_meter /(physconst("LightSpeed")/1.4677); + delay_n = delay_t .* obj.fs; + else + delay_n = options.delay_samples; + end + + obj.signal=delayseq(obj.signal,delay_n); + + end + function pow = power(obj) - pow = mean(abs(obj.signal.^2)) ; - % pow = pow2db(pow)+30; + pow = mean( abs(obj.signal).^2 ) ; + pow = pow2db(pow)+30; %dbm + + end + + function cspr = cspr(obj) + + carrier_power_dbm = pow2db( abs(mean(obj.signal)).^2 )+30; % dB -> +30 -> dBm + signal_power_dbm = obj.power; + + carr_lin = abs(mean(obj.signal)).^2; + sign_lin = mean(abs(obj.signal).^2); + + cspr_lin = pow2db(carr_lin/sign_lin); + + cspr =carrier_power_dbm-signal_power_dbm; + + end + + function papr = papr(obj) + + average_power = obj.power; + peak_power = pow2db(max(abs(obj.signal.^2)))+30; + papr = peak_power - average_power; end end diff --git a/Classes/00_signals/Signal.m b/Classes/00_signals/Signal.m index 39be34e..05a309d 100644 --- a/Classes/00_signals/Signal.m +++ b/Classes/00_signals/Signal.m @@ -12,6 +12,7 @@ classdef Signal %SIGNAL Construct an instance of this class % Detailed explanation goes here obj.signal = signal; + obj.signal = obj.signal; SignalType = []; TimeStamp = []; @@ -26,13 +27,18 @@ classdef Signal end %% CONVERT TO INFORMATIONSIGNAL - function [i_sig, varargout] = Informationsignal(obj) + function [i_sig, varargout] = Informationsignal(obj,options) + arguments + obj + options.fs + options.logbook + end if isa(obj,'Electricalsignal') - %convert to optical + %convert to information varargout{1} = obj.fs; - i_sig = Informationsignal(obj.signal); + i_sig = Informationsignal(obj.signal,"fs",options.fs,"logbook",options.logbook); elseif isa(obj,'Opticalsignal') @@ -109,12 +115,32 @@ classdef Signal end %% Add signals from one signal to another, the first object will sustain - function Sum = plus(X,Y) - Sum = X; %first input object will sustain - Sum.signal = X.signal + Y.signal; - end + function Sum = plus(X,y) - %% Display length + if isa(y,'Signal') + Sum = X; + Sum.signal = X.signal + y.signal; + elseif isnumeric(y) + Sum = X; + Sum.signal = X.signal + y; + end + + end + + function Product = times(X,y) + + if isa(y,'Signal') + Product = X; + Product.signal = X.signal .* y.signal; + elseif isnumeric(y) + Product = X; + Product.signal = X.signal .* y; + end + + end + + + %% Display length function return_length = length(obj) %METHOD1 Summary of this method goes here % Detailed explanation goes here @@ -142,7 +168,7 @@ classdef Signal end - %% Resample Signal + %% Resample Signal function obj = resample(obj,options) arguments @@ -157,6 +183,8 @@ classdef Signal obj = obj.logbookentry(desc); + obj.fs = options.fs_out; + end %% @@ -262,6 +290,7 @@ classdef Signal end + %% function obj = normalize(obj,options) @@ -270,7 +299,7 @@ classdef Signal options.mode normalization_mode = normalization_mode.rms end - switch options.mode + switch options.mode case normalization_mode.rms obj.signal = obj.signal/sqrt(mean(abs(obj.signal).^2,"all")); case normalization_mode.oneone @@ -283,20 +312,58 @@ classdef Signal function [obj,delay_n] = delay(obj,options) arguments - obj Opticalsignal - options.delay_meter double = 0 + obj Signal + options.delay_samples = 0 end - delay_t = options.delay_meter /(physconst("LightSpeed")/1.4677); - - delay_n = round(delay_t .* obj.fs); - - % finally circshift the signal -% obj.signal=circshift(obj.signal,delay_n); - - obj.signal=[zeros(delay_n,1); obj.signal(1:end-delay_n) ]; + obj.signal=delayseq(obj.signal,options.delay_samples); end + + %% + function [obj,D,cuts] = tsynch(obj,options) + % time sync and cut + arguments + obj Signal + options.reference Signal + options.fs_ref = 0; + end + + %normalize the signal + a = obj.normalize("mode","oneone").signal; + + %resample the reference + q = obj.fs/options.fs_ref; + b = options.reference.resample("fs_in",options.fs_ref,"fs_out",obj.fs).normalize("mode","oneone").signal; + + %estimate delay between signals + [co,lags] = xcorr(a,b,100); + [~,pos] = max(co); + D = lags(pos); + + if D == 100 + warning('Check Sync: Delay is found to be 100 which is the max. windowlength is xcorr function!') + end + + % delay by lagging samples + obj = obj.delay("delay_samples",-D); + + cuts = obj.length-(options.reference.length*q); + + if cuts > 10 + warning('Check Sync: Signal difference larger than 10.') + end + + if cuts < 0 + % warning('Check Sync: Reference Signal shorter than signal to sync.') + else + % then cut out the length of the reference + obj.signal = obj.signal(1:end-cuts); + end + + + end + end end diff --git a/Classes/01_transmit/AWG.m b/Classes/01_transmit/AWG.m index f7a8cd4..47ecd23 100644 --- a/Classes/01_transmit/AWG.m +++ b/Classes/01_transmit/AWG.m @@ -4,7 +4,6 @@ classdef AWG properties(Access=public) - preset kover %oversampling factor e.g. 16 repetitions %repeat the signal to generate a longer sequence? fdac %needed @@ -32,7 +31,6 @@ classdef AWG % Detailed explanation goes here arguments - options.preset = 'none'; options.kover = 16; options.repetitions = 1; options.normalize = 1; @@ -45,6 +43,7 @@ classdef AWG options.lpf_active = 0; options.lpf_type = 0; options.f_cutoff = 32e9; + options.H_lpf Filter end fn = fieldnames(options); @@ -54,40 +53,34 @@ classdef AWG end end - if options.preset == "M8196A" - % M8196A (92GBd) https://www.keysight.com/us/en/product/M8196A/92-gsa-s-arbitrary-waveform-generators.html - obj.dac_max = 0.5; - obj.dac_min = -.5; - obj.f_cutoff = 50e9; - elseif options.preset == "M8199B" - %https://www.keysight.com/us/en/assets/3120-1465/data-sheets/M8199A-128-256-GSa-s-Arbitrary-Waveform-Generator.pdf - % obj.fdac = 256e9; - obj.dac_max = 0.5; - obj.dac_min = -.5; - obj.f_cutoff = 80e9; - end - end function signalclass_out = process(obj,signalclass_in) - len_in = length(signalclass_in.signal); - % actual processing of the signal (steps 1. - 3.) + % 1-3. actual processing of the signal (normalize->quantize->sample hold) signalclass_in.signal = obj.process_(signalclass_in.signal); % 4. Apply LPF on the signal if obj.lpf_active - lpf = Filter('filtdegree',5,"f_cutoff",obj.f_cutoff,"fsamp",obj.kover*obj.fdac,"filterType",obj.lpf_type); - signalclass_in = lpf.process(signalclass_in); + if isa(obj.H_lpf,'Filter') + %4.A) user already specified a complete filter class when + % initializing the AWG module + signalclass_in = obj.H_lpf.process(signalclass_in); + else + %4.B) just use a standard filter + lpf = Filter('filtdegree',5,"f_cutoff",obj.f_cutoff,"fsamp",obj.kover*obj.fdac,"filterType",obj.lpf_type); + signalclass_in = lpf.process(signalclass_in); + end end % cast the inform. signal to electrical signal signalclass_in = Electricalsignal(signalclass_in,"fs",obj.fdac*obj.kover,"logbook",signalclass_in.logbook); % append to logbook - lbdesc = ['AWG preset ', obj.preset, 'k_over:',num2str(obj.kover),'. f_dac:',num2str(obj.fdac*1e-9),'GHz. Resolution:',num2str(obj.bit_resolution),' bits.']; + current_class = class(obj); + lbdesc = ['AWG ', current_class , '// k_over:',num2str(obj.kover),'. f_dac:',num2str(obj.fdac*1e-9),'GHz. Resolution:',num2str(obj.bit_resolution),' bits.']; signalclass_in = signalclass_in.logbookentry(lbdesc); % write to output diff --git a/Classes/01_transmit/M8196A.m b/Classes/01_transmit/M8196A.m new file mode 100644 index 0000000..8a0d7f8 --- /dev/null +++ b/Classes/01_transmit/M8196A.m @@ -0,0 +1,35 @@ +classdef M8196A < AWG + properties + + end + + methods + + function obj = M8196A() + %This is a ready to use AWG simulation that resembles the + %properties of the Keysighe M8196A % M8196A (92GBd) https://www.keysight.com/us/en/product/M8196A/92-gsa-s-arbitrary-waveform-generators.html + arguments + %options.bla = 1; + end + + obj = obj@AWG(); + obj.dac_max = 0.5; + obj.dac_min = -.5; + obj.bit_resolution = 5.5; + + obj.fdac = 92e9; + + obj.lpf_active = 1; + + obj.f_cutoff = 32e9; + obj.lpf_type = filtertypes.gaussian; + obj.H_lpf = Filter('filtdegree',5,"f_cutoff",obj.f_cutoff,"fsamp",obj.kover*obj.fdac,"filterType",obj.lpf_type); + + + end + + + end + + +end \ No newline at end of file diff --git a/Classes/01_transmit/M8199B.m b/Classes/01_transmit/M8199B.m new file mode 100644 index 0000000..448e692 --- /dev/null +++ b/Classes/01_transmit/M8199B.m @@ -0,0 +1,34 @@ +classdef M8199B < AWG + properties + + end + + methods + + function obj = M8199B() + %This is a ready to use AWG simulation that resembles the + %properties of the Keysighe M8199B https://www.keysight.com/us/en/assets/3120-1465/data-sheets/M8199A-128-256-GSa-s-Arbitrary-Waveform-Generator.pdf + arguments + %options.bla = 1; + end + + obj = obj@AWG(); + obj.dac_max = 0.5; + obj.dac_min = -.5; + obj.bit_resolution = 5.5; + + obj.fdac = 256e9; + + obj.lpf_active = 1; + + obj.f_cutoff = 80e9; + obj.lpf_type = filtertypes.gaussian; + obj.H_lpf = Filter('filtdegree',5,"f_cutoff",obj.f_cutoff,"fsamp",obj.kover*obj.fdac,"filterType",obj.lpf_type); + + end + + + end + + +end \ No newline at end of file diff --git a/Classes/01_transmit/PAMmapper.m b/Classes/01_transmit/PAMmapper.m index 50592c7..451c85b 100644 --- a/Classes/01_transmit/PAMmapper.m +++ b/Classes/01_transmit/PAMmapper.m @@ -19,10 +19,16 @@ classdef PAMmapper end - function signalclass_out = map(obj,signalclass_in) - signalclass_in.signal = obj.map_(signalclass_in.signal); - signalclass_in = signalclass_in.logbookentry(); - signalclass_out = signalclass_in; + function out = map(obj,signal_in) + + if isa(signal_in,'Signal') + signal_in.signal = obj.map_(signal_in.signal); + signal_in = signal_in.logbookentry(); + out = signal_in; + else + out = signal_in; + end + end function signalclass_out = demap(obj,signalclass_in) @@ -157,7 +163,7 @@ classdef PAMmapper case 2 % 4-ASK - data_out=[comp_real(:,:,2); ones(s1,s2)-comp_real(:,:,1)+comp_real(:,:,3)]; + data_out=[comp_real(:,:,2); ones(s1,s2) - comp_real(:,:,1) + comp_real(:,:,3)]; case 3 @@ -179,16 +185,43 @@ classdef PAMmapper end - function [data_out] = decide_pamlevel(obj,data_in) + function [data_out] = decide_pamlevel(obj,data_in,options) + arguments + obj + data_in Signal + options.symbol_levels = [] + end + + %A) normally return the preproduct of the decision a = squeeze(repmat(real(data_in.signal),[1 1 length(obj.thresholds)])); %Eingangssignal in 3 spalten b = squeeze(repmat(reshape(obj.thresholds(:).',[1 1 length(obj.thresholds)]),[1 length(data_in.signal) 1])); %Threshold in 3 Spalten comp_real = a > b; %check for each symbol/ sampling if it exeeds the obj.thresholdseshold 1, 2 or 3 + + data_out = data_in; + data_out.signal = sum(comp_real,2); + + %Option: return the actual level values/ just map onto given + %symbol levels + if ~isempty(options.symbol_levels) + data_out.signal = options.symbol_levels(data_out.signal+1); + end + + end - data_out = sum(comp_real,2); + function [out] = separate_pamlevels(obj,data_in) + + %A) normally return the preproduct of the decision + a = squeeze(repmat(real(data_in.signal),[1 1 length(obj.thresholds)])); %Eingangssignal in 3 spalten + b = squeeze(repmat(reshape(obj.thresholds(:).',[1 1 length(obj.thresholds)]),[1 length(data_in.signal) 1])); %Threshold in 3 Spalten + comp_real = a > b; %check for each symbol/ sampling if it exeeds the obj.thresholdseshold 1, 2 or 3 + comp_real_sum = sum(comp_real,2); - %data_out = (data_out*2)-3; + out = NaN(length(data_in),length(obj.thresholds)+1); + + for idx = 1:length(data_in) + out(idx,comp_real_sum(idx)+1) = data_in.signal(idx); + end - %data_out = pam_level_decision .* 1/sqrt(5); end end diff --git a/Classes/01_transmit/Pulseformer.m b/Classes/01_transmit/Pulseformer.m index a6f531b..e51c601 100644 --- a/Classes/01_transmit/Pulseformer.m +++ b/Classes/01_transmit/Pulseformer.m @@ -5,7 +5,7 @@ classdef Pulseformer properties(Access=public) fdac fsym - pulseform + pulse pulselength rrcalpha end @@ -18,7 +18,7 @@ classdef Pulseformer arguments options.fdac double options.fsym double - options.pulseform pulseform = pulseform.rrc + options.pulse pulseform = pulseform.rrc options.pulselength double {mustBeInteger} = 32 options.rrcalpha double = 0.05 end @@ -44,6 +44,9 @@ classdef Pulseformer lbdesc = 'Applied Pulseshaping'; signalclass_in = signalclass_in.logbookentry(lbdesc); + % write fs to signal + signalclass_in.fs = obj.fdac; + % write to output signalclass_out = signalclass_in; @@ -62,11 +65,11 @@ classdef Pulseformer % Detailed explanation goes here arguments(Input) obj - data_in double + data_in end arguments(Output) - data_out double + data_out end if ~rem(obj.fdac,obj.fsym) @@ -81,9 +84,9 @@ classdef Pulseformer sps= up; end - if obj.pulseform == pulseform.rrc + if obj.pulse == pulseform.rrc %Bau das Filter (hier rrc) - racos_len = obj.pulselength*2 ; + racos_len = obj.pulselength*2; alpha = obj.rrcalpha; h = rcosdesign(alpha,racos_len,sps); h = h./ max(h); @@ -115,12 +118,13 @@ classdef Pulseformer % %Apply Filter using Matlab build in fctn. -% data_out = upfirdn(data_in,h,up,dn); + h = rcosdesign(alpha,racos_len,sps); + h = h./ max(h); + %data_out_ = upfirdn(data_in,h,up,dn); % % %cut signal, which is longer due to fir filter st = round(up/dn*racos_len/2); %we need to cut y_out % en = round(st + (length(data_in)*up/dn) -1); -% % data_out = data_out(st:en); %scaling?! see pulsef module line 696 diff --git a/Classes/02_etc/Filter.m b/Classes/02_etc/Filter.m index a38cab9..35b017a 100644 --- a/Classes/02_etc/Filter.m +++ b/Classes/02_etc/Filter.m @@ -1,4 +1,4 @@ -classdef Filter +classdef Filter < handle %FILTER Summary of this class goes here % Detailed explanation goes here @@ -12,7 +12,7 @@ classdef Filter filtdegree passband_ripple stopband_ripple - fsamp + fs w end @@ -24,7 +24,8 @@ classdef Filter arguments options.filterType filtertypes = filtertypes.bessel_inp ; options.f_cutoff = 0; - options.fsamp = 0; + options.fs = 0; + options.signal_length = 0; options.filtdegree = 3; options.passband_ripple = 0.5; options.stopband_ripple = 0.5; @@ -36,7 +37,12 @@ classdef Filter for n = 1:numel(fn) obj.(fn{n}) = options.(fn{n}); end - + + % *Nur* wenn die beiden gegeben wurde, kann das Filter vorab + % erstellt werden - sonst zur Laufzeit + if obj.fs ~= 0 && obj.signal_length ~= 0 + [obj.H,obj.w] = obj.buildFilter(); + end end @@ -61,75 +67,77 @@ classdef Filter end function yout = process_(obj,xin) - + obj.signal_length = length(xin); - - [obj.H,obj.w] = obj.buildFilter(obj.filterType); - + + if isempty(obj.H) || obj.signal_length ~= length(obj.H) + [obj.H,obj.w] = obj.buildFilter(); + end + yout = obj.applyFilter(xin); - - end - + + end + function y_filtered = applyFilter(obj,xin) - + y_filtered = ifft(obj.H.*fft(xin)); - + end - - function [H,w] = buildFilter(obj,filterType) + + function [H,w] = buildFilter(obj) w = []; rp = obj.passband_ripple; %passband ripple rs = obj.stopband_ripple; %stopband ripple - switch filterType + switch obj.filterType case 1 % Bessel filter, impulse invariant transformed [B, A] = besself(obj.filtdegree, 2*pi*obj.f_cutoff); - [B ,A] = impinvar(B,A,obj.fsamp); + [B ,A] = impinvar(B,A,obj.fs); case 2 % Bessel filter, impulse bilinear transformed [Z, P, K] = besself(obj.filtdegree, 2*pi*obj.f_cutoff); - [Z ,P, K] = bilinear(Z,P,K,obj.fsamp); + [Z ,P, K] = bilinear(Z,P,K,obj.fs); [B ,A] = zp2tf(Z ,P ,K); case 3 % Butterworth filter if obj.lowpass == 1 %lowpass - [B, A] = butter(obj.filtdegree, obj.f_cutoff/(obj.fsamp/2),'low'); + [B, A] = butter(obj.filtdegree, obj.f_cutoff/(obj.fs/2),'low'); else % highpass - [B, A] = butter(obj.filtdegree, obj.f_cutoff/(obj.fsamp/2),'high'); + [B, A] = butter(obj.filtdegree, obj.f_cutoff/(obj.fs/2),'high'); end case 4 % Chebyshev 1 filter - [B, A] = cheby1(obj.filtdegree,rp, obj.f_cutoff/(obj.fsamp/2)); + [B, A] = cheby1(obj.filtdegree,rp, obj.f_cutoff/(obj.fs/2)); case 5 % Chebyshev 2 filter - [B, A] = cheby2(obj.filtdegree,rs, obj.f_cutoff/(obj.fsamp/2)); + [B, A] = cheby2(obj.filtdegree,rs, obj.f_cutoff/(obj.fs/2)); case 6 % Elliptic filter - [B, A] = ellip(obj.filtdegree,rp,rs,obj.f_cutoff/(obj.fsamp/2)); + [B, A] = ellip(obj.filtdegree,rp,rs,obj.f_cutoff/(obj.fs/2)); case 7 % Hamming filter g=(obj.filtdegree-1)/2; - wc=obj.f_cutoff/(obj.fsamp/2); + wc=obj.f_cutoff/(obj.fs/2); B = wc*sinc(wc*(-g:g)).*hamming(obj.filtdegree)'; A=1; @@ -137,7 +145,7 @@ classdef Filter % Raised Cosine filter - B = firrcos(obj.filtdegree,obj.f_cutoff,para.df,obj.fsamp); + B = firrcos(obj.filtdegree,obj.f_cutoff,df,obj.fs); A=1; case 9 @@ -145,7 +153,7 @@ classdef Filter % Sinc filter g=(obj.filtdegree-1)/2; - wc=obj.f_cutoff/(obj.fsamp/2); + wc=obj.f_cutoff/(obj.fs/2); B = wc*sinc(wc*(-g:g)); A=1; @@ -155,55 +163,122 @@ classdef Filter %check if order ist multiple of 1/2 blocklen=obj.signal_length; - faxis=linspace(-obj.fsamp/2,obj.fsamp/2,blocklen+1)';%generates arow vector faxis of blocklen+1 points linearly spaced between and including -para.fs/2 and para.fs/2 + faxis=linspace(-obj.fs/2,obj.fs/2,blocklen+1)';%generates arow vector faxis of blocklen+1 points linearly spaced between and including -fs/2 and fs/2 faxis=ifftshift(faxis(1:end-1)); H=exp(-((faxis)/(obj.f_cutoff*2)).^(2*obj.filtdegree)*log(2)*2^(2*obj.filtdegree-1)); -% figure() -% hold on -% xline(obj.f_cutoff*1e-9,'LineWidth',3,LineStyle='--'); -% xline(-obj.f_cutoff*1e-9,'LineWidth',3,LineStyle='--'); -% plot(faxis*1e-9,20*log10(abs(H)),'LineWidth',3); -% ax = gca; -% ylim([-6 0]) -% grid on -% xlabel('Freq in GHz') -% ylabel('Magnitude (dB)') - end % Build Filter from coefficients - if filterType ~= 10 + if obj.filterType ~= 10 [H,w] = freqz(B, A, obj.signal_length,'whole'); -% figure() -% hold on -% plot(w/(2*pi)*obj.fsamp*1e-9,20*log10(abs(H)),'LineWidth',3); -% ax = gca; -% %ylim([-6 0]) -% grid on -% xlabel('Freq in GHz') -% ylabel('Magnitude (dB)') - -% figure(30) -% freqz(B, A, 2048, obj.fsamp); - % hfvt = fvtool(B,A); end end function show(obj) - obj.signal_length = 1024; - [H,w] = obj.buildFilter(obj.filterType); - plot(w/pi,20*log10(abs(H))); - ax = gca; - ax.XTick = 0:.5:2; - grid on - xlabel('Normalized Frequency (\times\pi rad/sample)') - ylabel('Magnitude (dB)') - %freqz(H); + obj.signal_length = 512; + [H_,~] = obj.buildFilter(); + + figure() + sgtitle('- Filter Visualization -','FontSize',11) + hold on + fs = obj.fs; + frex = linspace(-fs/2,fs/2,numel(H_)).*1e-9; + + filt = fftshift(20*log10(abs(H_))); + filt_angle = fftshift(angle(H_)); + + filt_grpdly = diff(unwrap(filt_angle))./diff(frex.'); + + fc = obj.f_cutoff; + fcut_= [-fc*1e-9;fc*1e-9]; + + [~, index] = min(abs(filt - (-3))); + threedB = frex(index); + + cur_module_name = ''; + + filter_desc = char([num2str(obj.filtdegree),'th-order ',char(obj.filterType),'; f cut:',num2str(fc.*1e-9),' GHz']); + legend Interpreter none + + subplot(3,1,1); + legend Interpreter none + title('Magnitude') + hold on + p = plot(frex,filt,'LineWidth',1,'LineStyle','-','DisplayName',[cur_module_name,': ',filter_desc]); + xline([-threedB,threedB],'LineStyle','-','LineWidth',1,'HandleVisibility','off','Color',p.Color); + xline(fcut_,'LineStyle',':','LineWidth',1,'HandleVisibility','off','Color',p.Color); + yline(-3,'LineStyle',':','LineWidth',1,'HandleVisibility','off'); + ylim([-10 0]) + xlim([-200 200]) + ylabel('Magnitude (dB)') + legend + + subplot(3,1,2); + legend Interpreter none + title('Phase') + hold on + p = plot(frex,filt_angle,'LineWidth',1,'Color',p.Color,'LineStyle',':','DisplayName',[cur_module_name,': ',filter_desc]); + ylim([-pi,pi]); + yticks([-pi, 0, pi]); + yticklabels({'-$\pi$', '0', '$\pi$'}); + ylabel('Angle (rad)'); + + subplot(3,1,3); + legend Interpreter none + title('Group Delay') + hold on + p = plot(frex(2:end),filt_grpdly,'LineWidth',1,'Color',p.Color,'LineStyle',':','DisplayName',[cur_module_name,': ',filter_desc]); + ylabel('Angle (rad)'); + + + ax = gca; + + xlim([-200 200]) + grid on + xlabel('Freq in GHz') + + legend + + + + end + + function showHere(obj) + obj.signal_length = 512; + [H_,~] = obj.buildFilter(); + fs = obj.fs; + frex = linspace(-fs/2,fs/2,numel(H_)).*1e-9; + + filt = fftshift(20*log10(abs(H_))); + filt_angle = fftshift(angle(H_)); + + filt_grpdly = diff(unwrap(filt_angle))./diff(frex.'); + + fc = obj.f_cutoff; + fcut_= [-fc*1e-9;fc*1e-9]; + + [~, index] = min(abs(filt - (-3))); + threedB = frex(index); + + cur_module_name = ''; + + filter_desc = char([num2str(obj.filtdegree),'th-order ',char(obj.filterType),'; f cut:',num2str(fc.*1e-9),' GHz']); + legend Interpreter none + + legend Interpreter none + title('Magnitude') + hold on + p = plot(frex,filt,'LineWidth',1,'LineStyle','-','DisplayName',[cur_module_name,': ',filter_desc]); + xline([-threedB,threedB],'LineStyle','-','LineWidth',1,'HandleVisibility','off','Color',p.Color); + xline(fcut_,'LineStyle',':','LineWidth',1,'HandleVisibility','off','Color',p.Color); + yline(-3,'LineStyle',':','LineWidth',1,'HandleVisibility','off'); + legend + end end end diff --git a/Classes/02_optical/EML.m b/Classes/02_optical/EML.m index ebeb60b..96765da 100644 --- a/Classes/02_optical/EML.m +++ b/Classes/02_optical/EML.m @@ -14,6 +14,7 @@ classdef EML bias u_pi randomkey + randomstream %on instance creation field @@ -49,6 +50,8 @@ classdef EML obj.(fn{l}) = options.(fn{l}); end + obj.randomstream = RandStream('mlfg6331_64','Seed',obj.randomkey); + obj.field=sqrt(10^(obj.power/10-3)); %dbm to sqrt(mw) obj.noisefactor=sqrt(2*pi*obj.linewidth/obj.fsimu); @@ -105,8 +108,9 @@ classdef EML end function noi = createPhaseNoise(obj) + %create random vector - noi = randn(obj.signal_len/2,1); + noi = randn(obj.randomstream,obj.signal_len,1); %scale with noisefactor noi = noi * obj.noisefactor; @@ -114,9 +118,6 @@ classdef EML %cumsum to accumulate noise over time vector noi = cumsum(noi); - noi = [noi ; flip(noi)]; - - end function modulated_laserfield = externalmodulation(obj,laserfield,electrical_in) diff --git a/Classes/02_optical/Fiber.m b/Classes/02_optical/Fiber.m index b903d8c..eb0c67c 100644 --- a/Classes/02_optical/Fiber.m +++ b/Classes/02_optical/Fiber.m @@ -29,7 +29,7 @@ classdef Fiber options.D = 17 options.Dslope = 0.06 options.lambda0 = 1550 - options.gamma = 0.0013 + options.gamma = 0 options.dphimax = 5e-3 end diff --git a/Classes/03_receive/Photodiode.m b/Classes/03_receive/Photodiode.m index a1e5bdb..9b02436 100644 --- a/Classes/03_receive/Photodiode.m +++ b/Classes/03_receive/Photodiode.m @@ -8,6 +8,9 @@ classdef Photodiode dark_current temperature + randomkey + randomstream + end methods @@ -18,14 +21,16 @@ classdef Photodiode options.fsimu options.responsivity = 1; options.dark_current = 0; - options.temperature = 20; + options.temperature = 20; + options.randomkey = 1; end - obj.fsimu = options.fsimu; - obj.responsivity = options.responsivity; - obj.dark_current = options.dark_current; - obj.temperature = options.temperature; - + fn = fieldnames(options); + for l = 1:numel(fn) + obj.(fn{l}) = options.(fn{l}); + end + + obj.randomstream = RandStream('mlfg6331_64','Seed',obj.randomkey); end function signalclass_out = process(obj,signalclass_in) @@ -50,15 +55,16 @@ classdef Photodiode % Detailed explanation goes here k = Constant.Boltzmann; + q = Constant.ElementaryCharge; + T = obj.temperature + 273.15 ; %celsius + 273 = kelvin R = 50; %resistance of phdiode (50ohm is typical value) - % Magnitude squared detection - yout = sum( abs(xin) .^2*obj.responsivity, 2 ) ; - + % Magnitude squared detection (sum over pols) + yout = sum( abs(xin) .^2 * obj.responsivity, 2 ) ; % Shot Noise - shot_noise = sqrt(k * obj.fsimu .* yout) .* randn(size(yout,1),1); + shot_noise = sqrt(q * obj.fsimu .* yout) .* randn(obj.randomstream,size(yout,1),1); yout = yout + shot_noise; @@ -68,7 +74,7 @@ classdef Photodiode Bw = obj.fsimu; %is this correct? shouldnt it be the bandwidth of the actual component? e.g. 70GHz? therm_noise_pow = therm_current_psd * 2 * Bw; %squared - therm_noise = sqrt(therm_noise_pow) .* randn(size(yout,1),1); + therm_noise = sqrt(therm_noise_pow) .* randn(obj.randomstream,size(yout,1),1); yout = yout + therm_noise; diff --git a/Classes/03_receive/Scope.m b/Classes/03_receive/Scope.m index d7125f9..075f4c4 100644 --- a/Classes/03_receive/Scope.m +++ b/Classes/03_receive/Scope.m @@ -17,8 +17,10 @@ classdef Scope fixed_delay %fix the delay of the filter or use minimal delay for kausal system delay %specify a fixed delay of the filter + lpf_active filtertype lpf_bw + H_lpf block_dc @@ -47,8 +49,10 @@ classdef Scope options.fixed_delay = 0; options.delay = 0; + options.lpf_active = 0; options.filtertype = filtertypes.bessel_bilin; options.lpf_bw = 120e9; + options.H_lpf; options.block_dc = 1; end @@ -65,12 +69,27 @@ classdef Scope function signalclass_out = process(obj,signalclass_in) - % apply LPF - lpf = Filter('filtdegree',3,"f_cutoff",obj.lpf_bw,"fsamp",obj.fsimu,"filterType",obj.filtertype); - signalclass_in = lpf.process(signalclass_in); - % actual processing of the signal signalclass_in.signal = obj.process_(signalclass_in.signal); + signalclass_in.fs = obj.fadc; + + % apply LPF + % 4. Apply LPF on the signal + assert (signalclass_in.fs == obj.H_lpf.fs) + if obj.lpf_active + if isa(obj.H_lpf,'Filter') + %4.A) user already specified a complete filter class when + % initializing the AWG module + signalclass_in = obj.H_lpf.process(signalclass_in); + else + %4.B) just use a standard filter + obj.H_lpf = Filter('filtdegree',3,"f_cutoff",obj.lpf_bw,"fsamp",obj.fsimu,"filterType",obj.filtertype); + signalclass_in = obj.H_lpf.process(signalclass_in); + end + end + + % cast the electrical signal to information signal + signalclass_in = Informationsignal(signalclass_in,"fs",obj.fadc,"logbook",signalclass_in.logbook); % append to logbook lbdesc = ['Scope ']; diff --git a/Classes/04_DSP/EQ.m b/Classes/04_DSP/EQ.m index 3e2837d..c3facf0 100644 --- a/Classes/04_DSP/EQ.m +++ b/Classes/04_DSP/EQ.m @@ -52,6 +52,8 @@ classdef EQ coeff_number constellation_in + error_log + end methods @@ -258,6 +260,7 @@ classdef EQ e_dfe = b_.'*reference_vec; error = e_dc + e_ffe - e_dfe - ref_in(m-obj.k0); + %error = e_dc + e_.'*input_vec - b_.'*reference_vec - ref_in(m-obj.k0); if real(obj.FFEmu) @@ -279,7 +282,11 @@ classdef EQ end e_dc = e_dc - obj.DCmu*error; - error_log(cnt,trainloops) = e_dc; + + obj.error_log.e_ffe(cnt,trainloops) = e_ffe; + obj.error_log.e_dfe(cnt,trainloops) = e_dfe; + obj.error_log.e_(cnt,trainloops) = error; + cnt = cnt+1; if obj.Nb(1) > 0 b_ = b_ + obj.DFEmu*error*reference_vec; % Seems like normalized DFE has worse performance diff --git a/Classes/04_DSP/EQ_silas.m b/Classes/04_DSP/EQ_silas.m index 695684e..497dbe9 100644 --- a/Classes/04_DSP/EQ_silas.m +++ b/Classes/04_DSP/EQ_silas.m @@ -51,9 +51,9 @@ classdef EQ_silas < handle trainloops ddloops - eq_blocklength % block lengt of EQ (until now, only the dc subtraction is affected by this) + eq_parallelization_blocklength % block lengt of EQ (until now, only the dc subtraction is affected by this) eq_updatelatency % time in symbols until the calculated updates reach the signal again (until now, only the dc subtraction is affected by this) - + eq_avg_blocklength end @@ -79,8 +79,9 @@ classdef EQ_silas < handle options.mu_ffe_dd = [0.0004 0.0005 0.0006]; options.mu_dfe_dd = 0.0005; - options.eq_blocklength = 1; + options.eq_parallelization_blocklength = 1; options.eq_updatelatency = 1; + options.eq_avg_blocklength = 0; end fn = fieldnames(options); @@ -151,8 +152,8 @@ classdef EQ_silas < handle %% Adaptive Equalization Modes function trainingMode(obj) - - dc_block = ones(obj.eq_blocklength,1); + + dc_block = ones(obj.eq_parallelization_blocklength,1); for tloop = 1:obj.trainloops m = 1+obj.delay; @@ -176,7 +177,7 @@ classdef EQ_silas < handle % Calculate the Error obj.error = obj.e_dc + obj.e_ffe - obj.e_dfe - obj.d(obj.Nb(1)-1+m-obj.delay); - + if obj.mu_ffe_train ~= 0 %update FFE coefficients with LMS @@ -191,8 +192,8 @@ classdef EQ_silas < handle %update DC error dc_block(dc_cnt) = obj.error .* obj.mu_dc_train; - - if dc_cnt == obj.eq_blocklength + + if dc_cnt == obj.eq_parallelization_blocklength obj.e_dc = obj.e_dc - mean(dc_block(dc_cnt)); dc_cnt = 0; end @@ -208,8 +209,10 @@ classdef EQ_silas < handle %start the dd mode with coefficients from training coeff = [obj.e;obj.b]; + obj.e_dc = ones(obj.eq_updatelatency,1).*obj.e_dc; + dc_block = ones(obj.eq_parallelization_blocklength,1); + - dc_block = ones(obj.eq_blocklength,1); for ddloop = 1:obj.ddloops @@ -227,50 +230,96 @@ classdef EQ_silas < handle mu_dfe = ones(1,sum(obj.Cb))*obj.mu_dfe_dd; - y = zeros(1,floor(obj.x_length/obj.sps)); d_feedback = zeros(obj.Cb(1),1); d_vnle = obj.calcVNLENonlinVecs(d_feedback,obj.Ib2,obj.Ib3,obj.Nb,obj.d_norm); d_hat = zeros(obj.x_length,1); + m_reg = 0; + + if obj.eq_avg_blocklength > 0 + averaging_window = zeros(obj.eq_avg_blocklength,1); + end + for k = 1:obj.sps:obj.x_length dc_cnt = dc_cnt+1; m=m+1; %get Sigal input vectors with correct length for VNLE x = obj.x_in(obj.Ne(1)+k-1:-1:k).'; +% + if obj.eq_avg_blocklength > 0 %% Das läuft gut mit 400er Fenster!! + averaging_window = circshift(averaging_window,obj.sps); + averaging_window(1:obj.sps,1) = x(1:obj.sps); + avg_(k) = mean(averaging_window); + x = x-avg_(k); + end + %bring this signal to "special" VNLE format x_vnle = obj.calcVNLENonlinVecs(x,obj.Ie2,obj.Ie3,obj.Ne,obj.x_norm); %combine FFE with DFE to one vector (cursor between the two sequences) x_d = [x_vnle;-d_vnle]; + %Apply filter - y(m) = obj.e_dc + x_d.'* coeff; + %y(m) = (m_reg(end)*dc_cnt + obj.e_dc(end)) + x_d.'* coeff; + if obj.mu_dc_dd > 0 + y(m) = obj.e_dc(end) + x_d.'* coeff; + else + y(m) = x_d.'* coeff; + end + +% if obj.eq_avg_blocklength > 0 %% Das läuft nicht gut!! +% averaging_window = circshift(averaging_window,obj.sps); +% averaging_window(1:obj.sps,1) = y(m); +% avg_(m) = mean(averaging_window); +% y(m) = y(m)-avg_(m); +% end %Decision [~,symbol_idx] = min(abs(y(m) - obj.d_constellation)); % decision for closest constellation point d_hat(k) = obj.d_constellation(symbol_idx); %Error between FFE & DFE filtered signal and Decision - obj.error = y(m) - d_hat(k); + obj.error(k) = y(m) - d_hat(k); - %Update coefficients (both FFE and DFE) - % obj.e = obj.e - obj.error * mu_ffe * conj(x_vnle); - % obj.b = obj.b + obj.error * mu_dfe * conj(d_vnle); - % coeff = [obj.e;obj.b]; +% if obj.eq_avg_blocklength > 0 %% Das läuft nicht gut!! +% averaging_window = circshift(averaging_window,obj.sps); +% averaging_window(1:obj.sps,1) = y(m); +% avg_(m) = mean(averaging_window); +% y(m) = y(m)-avg_(m); +% end - coeff = coeff - mu_mat*obj.error*conj(x_d); - - - %update DC error - dc_block(dc_cnt) = obj.error .* obj.mu_dc_dd; + %Update FFE and DFE coefficients + coeff = coeff - mu_mat*obj.error(k) * conj(x_d); + + %Update DC error + dc_block(dc_cnt) = obj.error(k) ; + + if dc_cnt == obj.eq_parallelization_blocklength + if obj.eq_updatelatency > 1 + + obj.e_dc = circshift(obj.e_dc,1); + +% m_reg(end+1) = ((1:obj.eq_parallelization_blocklength)' \ (cumsum(dc_block))); +% +% obj.e_dc(1) = obj.e_dc(2) - sign(m_reg(end)) .* (sum(dc_block).* m_reg(end) .* obj.mu_dc_dd); + obj.e_dc(1) = obj.e_dc(2) - sum(dc_block) .* obj.mu_dc_dd; + + else + %m_reg(end+1) = ((1:obj.eq_parallelization_blocklength)' \ (cumsum(dc_block))); + +% obj.e_dc = obj.e_dc - sign(m_reg(end)) .* (sum(dc_block).* m_reg(end) .* obj.mu_dc_dd); + + obj.e_dc = obj.e_dc - sum(dc_block) .* obj.mu_dc_dd; + end - if dc_cnt == obj.eq_blocklength - obj.e_dc = obj.e_dc - mean(dc_block(dc_cnt)); dc_cnt = 0; + end + % Append new decision to decision feedback if obj.Nb(1) > 0 diff --git a/Classes/04_DSP/EQ_silas_plain.m b/Classes/04_DSP/EQ_silas_plain.m index 13ccb44..798a45c 100644 --- a/Classes/04_DSP/EQ_silas_plain.m +++ b/Classes/04_DSP/EQ_silas_plain.m @@ -1,369 +1,395 @@ -classdef EQ_silas_plain < handle - %EQ_SILAS FFE and DFE Equalizer Playground - - properties - % Important Signals - x_in %Input Sequence to be equalized - x_length - x_norm - - d %reference signal - d_norm - d_constellation %constellation points of the reference - - y_out %equalizer output signal - - % FFE coefficients always named with "e" - Ne - Ce %memory length FFE - Ie1 %Indice Combination of 1nd order FFE - Ie2 %Indice Combination of 2nd order FFE - Ie3 %Indice Combination of 3nd order FFE - e %coefficients for FFE - - % DFE coefficients always named with "b" - Nb - Cb %memory length DFE - Ib1 %Indice Combination of 1nd order DFE - Ib2 %Indice Combination of 2nd order DFE - Ib3 %Indice Combination of 3nd order DFE - b %coefficients for DFE - - error - e_ffe - e_dfe - e_dc - error_log - - mu_dc_train - mu_ffe_train - mu_dfe_train - - mu_dc_dd - mu_combined_dd - - delay - trainlength - sps - - trainloops - ddloops - - - - - end - - methods - function obj = EQ_silas_plain(options) - %EQ_SILAS Construct an instance of this class - arguments(Input) - options.Ne = [50 5 0] %Number of FFE coefficients (1st, 2nd and 3rd order) - options.Nb = [30 5 3] %Number of DFE coefficients (1st, 2nd and 3rd order) - options.trainloops = 2; - options.trainlength = 4096; - options.ddloops = 2; - - options.delay = 0; - options.sps = 2; - - options.mu_dc_train = 0.01; - options.mu_ffe_train = 0.005; - options.mu_dfe_train = 0.005; - - options.mu_dc_dd = 0.01; - options.mu_combined_dd = [0.0004 0.0005 0.0006 0.0007 ]; - end - - fn = fieldnames(options); - for n = 1:numel(fn) - obj.(fn{n}) = options.(fn{n}); - end - - % Generate helpful vectors and initialize the filters with - % correct length: - - obj.Ce = obj.calcVNLEMemoryLength(obj.Ne); - - [obj.Ie2,obj.Ie3] = obj.calcIndiceVectors(obj.Ne); - - obj.e = zeros(sum(obj.Ce),1); - - - obj.Cb = obj.calcVNLEMemoryLength(obj.Nb); - - [obj.Ib2,obj.Ib3] = obj.calcIndiceVectors(obj.Nb); - - obj.b = zeros(sum(obj.Cb),1); - - end - - function [signalclass_out,error_log] = process(obj,signalclass_in, reference_signalclass_in) - - % actual processing of the signal (steps 1. - 3.) - % 1 normalize RMS - signalclass_in = signalclass_in.normalize("mode","rms"); - - % Process the EQ optimization - obj.process_(signalclass_in.signal', reference_signalclass_in.signal'); - - signalclass_in.signal = obj.y_out'; - % append to logbook - lbdesc = ['EQ von Silas ist gelaufen ']; - signalclass_in = signalclass_in.logbookentry(lbdesc); - - % write to output - signalclass_out = signalclass_in; - - end - - function process_(obj,x_in,d_in) - - % 1) prepare signals - obj.e_dc = mean(x_in); - - % 1.1) Input Signal - obj.x_in = [zeros(1,floor(obj.Ne(1)/2)) x_in zeros(1,obj.Ne(1))]; - obj.x_length = length(x_in); - obj.x_norm = obj.calcPowerNormalization(x_in); - - % 1.2 Reference Signal // Constellation - obj.d = [zeros(1,obj.Nb(1)-1) d_in zeros(1,obj.Nb(1))]; - obj.d_constellation = unique(d_in); - obj.d_norm = obj.calcPowerNormalization(d_in); - - % 1.3 Training - obj.trainingMode(); - - % 1.4 Decision Directed Mode - obj.decisionDirectedMode(); - - end - - %% Adaptive Equalization Modes - - function trainingMode(obj) - - for tloop = 1:obj.trainloops - m = 1+obj.delay; - - for n = obj.sps*obj.delay+1:obj.sps:obj.sps*obj.trainlength - m = m+1; - - %get Sigal input vectors with correct length for VNLE - x_in_block = obj.x_in(obj.Ne(1)+n+(obj.sps-1):-1:n+obj.sps).'; - x_in_vnle_format = obj.calcVNLENonlinVecs(x_in_block,obj.Ie2,obj.Ie3,obj.Ne,[1,1,1]); - - %get Reference input vectors with correct length for VNLE - d_block = obj.d(obj.Nb(1)-obj.delay+m-2:-1:m-obj.delay-1).'; - d_vnle_format = obj.calcVNLENonlinVecs(d_block,obj.Ib2,obj.Ib3,obj.Nb,obj.d_norm); - - % Calculate the Error - obj.e_ffe = obj.e.' * x_in_vnle_format; - obj.e_dfe = obj.b.' * d_vnle_format; - - obj.error = obj.e_dc + obj.e_ffe - obj.e_dfe - obj.d(obj.Nb(1)-1+m-obj.delay); - - %update FFE coefficients with LMS - obj.e = obj.e - obj.error*conj(x_in_vnle_format)*obj.mu_ffe_train; - - %update DFE coefficients with LMS - obj.b = obj.b + obj.mu_dfe_train*obj.error*d_vnle_format; - - %update DC error - obj.e_dc = obj.e_dc - obj.error .* obj.mu_dc_train; - - end - - end - - end - - function decisionDirectedMode(obj) - - %start the dd mode with coefficients from training - coeff = [obj.e;obj.b]; - - for ddloop = 1:obj.ddloops - - m = 0; - - if all(obj.mu_combined_dd == obj.mu_combined_dd(1)) - mu_mat = obj.mu_combined_dd(1); - else - mu_mat = diag([ones(1,obj.Ce(1))*obj.mu_combined_dd(1)... %1st order ffe - ones(1,obj.Ce(2))*obj.mu_combined_dd(2)... %2nd order ffe - ones(1,obj.Ce(3))*obj.mu_combined_dd(3)... %3rd order ffe - ones(1,sum(obj.Cb))*obj.mu_combined_dd(4)]); %all order dfe - end - - y = zeros(1,floor(obj.x_length/obj.sps)); - d_feedback = zeros(obj.Cb(1),1); - d_vnle = obj.calcVNLENonlinVecs(d_feedback,obj.Ib2,obj.Ib3,obj.Nb,obj.d_norm); - d_hat = zeros(obj.x_length,1); - - for k = 1:obj.sps:obj.x_length - - m=m+1; - - %get Sigal input vectors with correct length for VNLE - x = obj.x_in(obj.Ne(1)+k-1:-1:k).'; - x_vnle = obj.calcVNLENonlinVecs(x,obj.Ie2,obj.Ie3,obj.Ne,[1,1,1]); - - %combine FFE with DFE to one vector (cursor between the two sequences) - x_d = [x_vnle;-d_vnle]; - - %Apply filter - y(m) = obj.e_dc + x_d.'* coeff; - - %Decision - [~,symbol_idx] = min(abs(y(m) - obj.d_constellation)); % decision for closest constellation point - d_hat(k) = obj.d_constellation(symbol_idx); - - %Error between FFE & DFE filtered signal and Decision - obj.error = y(m) - d_hat(k); - - %Update coefficients (both FFE and DFE) - coeff = coeff - mu_mat*obj.error*conj(x_d); - - if 1 %mu_mat ~= 0 - obj.e_dc = obj.e_dc - obj.mu_dc_dd * obj.error; - obj.error_log(ddloop,m) = obj.e_dc.^2; - end - - % Append new decision to decision feedback - if obj.Nb(1) > 0 - - %shift up one index - d_feedback(2:end) = d_feedback(1:end-1); - %replace 1st index with current estimation - d_feedback(1) = d_hat(k); - %build memorylike VNLE version - d_vnle = obj.calcVNLENonlinVecs(d_feedback,obj.Ib2,obj.Ib3,obj.Nb,obj.d_norm); - - end - - end - end - - obj.y_out = (circshift( y.' ,-(obj.delay))).'; - - end - - %% Functions needed During Adaption - function x_in_vnle_format = calcVNLENonlinVecs(~,x_in_block,I_2,I_3,N_,norm_) - % These are the second and third order input signal products of the VNLE EQ - % ∑ h1 x_in(k-n1) + ∑∑ h2 x_in(k-n1)*x_in(k-n2) + ∑∑∑ h3 x_in(k-n1)*x_in(k-n2)*x_in(k-n3) - - x1 = x_in_block; - x2 = []; - x3 = []; - - if N_(2) > 0 - delta_2 = round((N_(1)-N_(2))/2); - input_vec_se = x_in_block(delta_2:end)/norm_(2); %TODO normalization step - x2 = input_vec_se(I_2(:,1)).*input_vec_se(I_2(:,2)); - end - - if N_(3) > 0 - delta_3 = round((N_(1)-N_(3))/2); - input_vec_th = x_in_block(delta_3:end)/norm_(3); - x3 = input_vec_th(I_3(:,1)).*input_vec_th(I_3(:,2)).*input_vec_th(I_3(:,3)); - end - - x_in_vnle_format = [x1;x2;x3]; - - end - - - %% Functions needed for Preparation - function [C] = calcVNLEMemoryLength(~,N) - - %calculates the memory length of VNLE - C = zeros(size(N)); - - for o = 1:numel(N) - switch o - case 1 - C(o) = N(o); - case 2 - C(o) = N(o)*(N(o)+1) / 2; - case 3 - C(o) = N(o)*(N(o)+1)*(N(o)+2) / 6; - end - end - - end - - function [indvec2nd, indvec3rd] = calcIndiceVectors(~,N) - - % Init vectors of 2nd and 3rd order coefficient indices -> - % yield combination with - - for order = 2:numel(N) - n = N(order); - v = 1:n; % Ursprünglicher Vektor - row = 1; - - % Schleifen zur Generierung des Indize Vektors - switch order - - case 2 - - indvec2nd = zeros(n*(n+1)/2, order); - for i = 1:n - for j = i:n - indvec2nd(row, :) = [v(i) v(j)]; - row = row + 1; - end - end - - case 3 - - indvec3rd = zeros(n*(n+1)*(n+2)/6, 3); - for i = 1:n - for j = i:n - for k = j:n - indvec3rd(row, :) = [v(i) v(j) v(k)]; - row = row + 1; - end - end - end - end - end - - end - - function powerNorm = calcPowerNormalization(~,v) - - powerNorm(1) = sqrt(mean(abs(v ).^2)); - powerNorm(2) = sqrt(mean(abs(v.^2).^2)); - powerNorm(3) = sqrt(mean(abs(v.^3).^2)); - - end - - - end -end - - - - - - - - - - - - - - - - - - - - - - +classdef EQ_silas_plain < handle + %EQ_SILAS FFE and DFE Equalizer Playground + + properties + % Important Signals + x_in %Input Sequence to be equalized + x_length + x_norm + + d %reference signal + d_norm + d_constellation %constellation points of the reference + + y_out %equalizer output signal + + % FFE coefficients always named with "e" + Ne + Ce %memory length FFE + Ie1 %Indice Combination of 1nd order FFE + Ie2 %Indice Combination of 2nd order FFE + Ie3 %Indice Combination of 3nd order FFE + e %coefficients for FFE + + % DFE coefficients always named with "b" + Nb + Cb %memory length DFE + Ib1 %Indice Combination of 1nd order DFE + Ib2 %Indice Combination of 2nd order DFE + Ib3 %Indice Combination of 3nd order DFE + b %coefficients for DFE + + error + e_ffe + e_dfe + e_dc + error_log + + mu_dc_train + mu_ffe_train + mu_dfe_train + + mu_dc_dd + mu_combined_dd + + delay + trainlength + sps + + trainloops + ddloops + + + + + end + + methods + function obj = EQ_silas_plain(options) + %EQ_SILAS Construct an instance of this class + arguments(Input) + options.Ne = [50 5 0] %Number of FFE coefficients (1st, 2nd and 3rd order) + options.Nb = [30 5 3] %Number of DFE coefficients (1st, 2nd and 3rd order) + options.trainloops = 2; + options.trainlength = 4096; + options.ddloops = 2; + + options.delay = 0; + options.sps = 2; + + options.mu_dc_train = 0.01; + options.mu_ffe_train = 0.005; + options.mu_dfe_train = 0.005; + + options.mu_dc_dd = 0.01; + options.mu_combined_dd = [0.0004 0.0005 0.0006 0.0007 ]; + end + + fn = fieldnames(options); + for n = 1:numel(fn) + obj.(fn{n}) = options.(fn{n}); + end + + % Generate helpful vectors and initialize the filters with + % correct length: + + obj.Ce = obj.calcVNLEMemoryLength(obj.Ne); + + [obj.Ie2,obj.Ie3] = obj.calcIndiceVectors(obj.Ne); + + obj.e = zeros(sum(obj.Ce),1); + + obj.Cb = obj.calcVNLEMemoryLength(obj.Nb); + + [obj.Ib2,obj.Ib3] = obj.calcIndiceVectors(obj.Nb); + + obj.b = zeros(sum(obj.Cb),1); + + end + + function [signalclass_out,error_log] = process(obj,signalclass_in, reference_signalclass_in) + + assert(signalclass_in.fs / reference_signalclass_in.fs == obj.sps) + + % actual processing of the signal (steps 1. - 3.) + % 1 normalize RMS + signalclass_in = signalclass_in.normalize("mode","rms"); + + % Process the EQ optimization + obj.process_(signalclass_in.signal', reference_signalclass_in.signal'); + + signalclass_in.signal = obj.y_out'; + + %change sampling frequency of outgoing signal + signalclass_in.fs = reference_signalclass_in.fs; + + % append to logbook + lbdesc = ['EQ von Silas ist gelaufen ']; + signalclass_in = signalclass_in.logbookentry(lbdesc); + + % write to output + signalclass_out = signalclass_in; + + end + + function process_(obj,x_in,d_in) + + % 1) prepare signals + obj.e_dc = mean(x_in); + + % 1.1) Input Signal + obj.x_in = [zeros(1,floor(obj.Ne(1)/2)) x_in zeros(1,obj.Ne(1))]; + obj.x_length = length(x_in); + obj.x_norm = obj.calcPowerNormalization(x_in); + + % 1.2 Reference Signal // Constellation + obj.d = [zeros(1,obj.Nb(1)-1) d_in zeros(1,obj.Nb(1))]; + obj.d_constellation = unique(d_in); + obj.d_norm = obj.calcPowerNormalization(d_in); + + % 1.3 Training + obj.trainingMode(); + + % 1.4 Decision Directed Mode + obj.decisionDirectedMode(); + + end + + %% Adaptive Equalization Modes + + function trainingMode(obj) + + for tloop = 1:obj.trainloops + m = 1+obj.delay; + cnt = 1; + for n = obj.sps*obj.delay+1:obj.sps:obj.sps*obj.trainlength + m = m+1; + + %get Sigal input vectors with correct length for VNLE + x_in_block = obj.x_in(obj.Ne(1)+n+(obj.sps-1):-1:n+obj.sps).'; + x_in_vnle_format = obj.calcVNLENonlinVecs(x_in_block,obj.Ie2,obj.Ie3,obj.Ne,obj.x_norm); + + %get Reference input vectors with correct length for VNLE + d_block = obj.d(obj.Nb(1)-obj.delay+m-2:-1:m-obj.delay-1).'; + d_vnle_format = obj.calcVNLENonlinVecs(d_block,obj.Ib2,obj.Ib3,obj.Nb,obj.d_norm); + + % Calculate the Error + obj.e_ffe = obj.e.' * x_in_vnle_format; + obj.e_dfe = obj.b.' * d_vnle_format; + + obj.error = obj.e_dc + obj.e_ffe - obj.e_dfe - obj.d(obj.Nb(1)-1+m-obj.delay); + obj.error_log(tloop,cnt) = obj.error; + + %update FFE coefficients with LMS + obj.e = obj.e - obj.error*conj(x_in_vnle_format)*obj.mu_ffe_train; + + %update DFE coefficients with LMS + obj.b = obj.b + obj.mu_dfe_train*obj.error*d_vnle_format; + + %update DC error + obj.e_dc = obj.e_dc - obj.error .* obj.mu_dc_train; + cnt = cnt+1; + end + + end + + end + + function decisionDirectedMode(obj) + + %start the dd mode with coefficients from training + coeff = [obj.e;obj.b]; + + for ddloop = 1:obj.ddloops + + m = 0; + + if all(obj.mu_combined_dd == obj.mu_combined_dd(1)) + mu_mat = obj.mu_combined_dd(1); + else + mu_mat = diag([ones(1,obj.Ce(1))*obj.mu_combined_dd(1)... %1st order ffe + ones(1,obj.Ce(2))*obj.mu_combined_dd(2)... %2nd order ffe + ones(1,obj.Ce(3))*obj.mu_combined_dd(3)... %3rd order ffe + ones(1,sum(obj.Cb))*obj.mu_combined_dd(4)]); %all order dfe + end + + y = zeros(1,floor(obj.x_length/obj.sps)); + d_feedback = zeros(obj.Cb(1),1); + d_vnle = obj.calcVNLENonlinVecs(d_feedback,obj.Ib2,obj.Ib3,obj.Nb,obj.d_norm); + d_hat = zeros(obj.x_length,1); + + for k = 1:obj.sps:obj.x_length + + m=m+1; + + %get Sigal input vectors with correct length for VNLE + x = obj.x_in(obj.Ne(1)+k-1:-1:k).'; + x_vnle = obj.calcVNLENonlinVecs(x,obj.Ie2,obj.Ie3,obj.Ne,obj.x_norm); + + %combine FFE with DFE to one vector (cursor between the two sequences) + x_d = [x_vnle;-d_vnle]; + + %Apply filter + y(m) = obj.e_dc + x_d.'* coeff; + + %Decision + [~,symbol_idx] = min(abs(y(m) - obj.d_constellation)); % decision for closest constellation point + d_hat(k) = obj.d_constellation(symbol_idx); + + %Error between FFE & DFE filtered signal and Decision + obj.error = y(m) - d_hat(k); + + %Update coefficients (both FFE and DFE) + coeff = coeff - (mu_mat * (obj.error * conj(x_d))); + + if 1 %mu_mat ~= 0 + obj.e_dc = obj.e_dc - obj.mu_dc_dd * obj.error; + end + + % Append new decision to decision feedback + if obj.Nb(1) > 0 + + %shift up one index + d_feedback(2:end) = d_feedback(1:end-1); + %replace 1st index with current estimation + d_feedback(1) = d_hat(k); + %build memorylike VNLE version + d_vnle = obj.calcVNLENonlinVecs(d_feedback,obj.Ib2,obj.Ib3,obj.Nb,obj.d_norm); + + end + + end + end + + obj.y_out = (circshift( y.' ,-(obj.delay))).'; + + end + + %% Functions needed During Adaption + function x_in_vnle_format = calcVNLENonlinVecs(~,x_in_block,I_2,I_3,N_,norm_) + % These are the second and third order input signal products of the VNLE EQ + % ∑ h1 x_in(k-n1) + ∑∑ h2 x_in(k-n1)*x_in(k-n2) + ∑∑∑ h3 x_in(k-n1)*x_in(k-n2)*x_in(k-n3) + l1=length(x_in_block); + l2=length(I_2); + l3=length(I_3); + final_length = l1+l2+l3; + + x_in_vnle_format = zeros(final_length,1); + + idx = l1; + x_in_vnle_format(1:idx) = x_in_block; + + if N_(2) > 0 + delta_2 = round((N_(1)-N_(2)) / 2); + input_vec_se = x_in_block(delta_2:end) / norm_(2); %TODO normalization step + + % Extract columns from I_2 + col1 = input_vec_se(I_2(:,1)); + col2 = input_vec_se(I_2(:,2)); + + x2 = col1 .* col2; + x_in_vnle_format(idx+1:idx+l2) = x2; + end + + if N_(3) > 0 + delta_3 = round((N_(1)-N_(3))/2); + input_vec_th = x_in_block(delta_3:end) / norm_(3); + + % Extract columns from I_3 + col1 = input_vec_th(I_3(:,1)); + col2 = input_vec_th(I_3(:,2)); + col3 = input_vec_th(I_3(:,3)); + + % Perform matrix multiplication + x3 = col1 .* col2 .* col3; + + idx = idx+l2; + x_in_vnle_format(idx+1:idx+l3) = x3; + end + + + end + + + + %% Functions needed for Preparation + function [C] = calcVNLEMemoryLength(~,N) + + %calculates the memory length of VNLE + C = zeros(size(N)); + + for o = 1:numel(N) + switch o + case 1 + C(o) = N(o); + case 2 + C(o) = N(o)*(N(o)+1) / 2; + case 3 + C(o) = N(o)*(N(o)+1)*(N(o)+2) / 6; + end + end + + end + + function [indvec2nd, indvec3rd] = calcIndiceVectors(~,N) + + % Init vectors of 2nd and 3rd order coefficient indices -> + % yield combination with + + for order = 2:numel(N) + n = N(order); + v = 1:n; % Ursprünglicher Vektor + row = 1; + + % Schleifen zur Generierung des Indize Vektors + switch order + + case 2 + + indvec2nd = zeros(n*(n+1)/2, order); + for i = 1:n + for j = i:n + indvec2nd(row, :) = [v(i) v(j)]; + row = row + 1; + end + end + + case 3 + + indvec3rd = zeros(n*(n+1)*(n+2)/6, 3); + for i = 1:n + for j = i:n + for k = j:n + indvec3rd(row, :) = [v(i) v(j) v(k)]; + row = row + 1; + end + end + end + end + end + + end + + function powerNorm = calcPowerNormalization(~,v) + + powerNorm(1) = sqrt(mean(abs(v ).^2)); + powerNorm(2) = sqrt(mean(abs(v.^2).^2)); + powerNorm(3) = sqrt(mean(abs(v.^3).^2)); + + end + + + end +end + + + + + + + + + + + + + + + + + + + + + + diff --git a/Classes/04_DSP/EQ_silas_sliding_window_dc_removal.m b/Classes/04_DSP/EQ_silas_sliding_window_dc_removal.m new file mode 100644 index 0000000..0787434 --- /dev/null +++ b/Classes/04_DSP/EQ_silas_sliding_window_dc_removal.m @@ -0,0 +1,408 @@ +classdef EQ_silas_sliding_window_dc_removal < handle + %EQ_SILAS FFE and DFE Equalizer Playground + + properties + % Important Signals + x_in %Input Sequence to be equalized + x_length + x_norm + + d %reference signal + d_norm + d_constellation %constellation points of the reference + + y_out %equalizer output signal + + % FFE coefficients always named with "e" + Ne + Ce %memory length FFE + Ie1 %Indice Combination of 1nd order FFE + Ie2 %Indice Combination of 2nd order FFE + Ie3 %Indice Combination of 3nd order FFE + e %coefficients for FFE + + % DFE coefficients always named with "b" + Nb + Cb %memory length DFE + Ib1 %Indice Combination of 1nd order DFE + Ib2 %Indice Combination of 2nd order DFE + Ib3 %Indice Combination of 3nd order DFE + b %coefficients for DFE + + error + e_ffe + e_dfe + e_dc + + % coefficients + mu_dc_train + mu_ffe_train + mu_dfe_train + + mu_dc_dd + mu_ffe_dd + mu_dfe_dd + mu_combined_dd % [1st order FFE, 2nd order FFE, 3rd order FFE, all orders DFE] + + delay + trainlength + sps + + trainloops + ddloops + + eq_blocklength % block lengt of EQ (until now, only the dc subtraction is affected by this) + eq_updatelatency % time in symbols until the calculated updates reach the signal again (until now, only the dc subtraction is affected by this) + + + + end + + methods + function obj = EQ_silas_sliding_window_dc_removal(options) + %EQ_SILAS Construct an instance of this class + arguments(Input) + options.Ne = [50 5 0] %Number of FFE coefficients (1st, 2nd and 3rd order) + options.Nb = [30 5 3] %Number of DFE coefficients (1st, 2nd and 3rd order) + options.trainloops = 2; + options.trainlength = 4096; + options.ddloops = 2; + + options.delay = 0; + options.sps = 2; + + options.mu_dc_train = 0.01; + options.mu_ffe_train = 0.005; + options.mu_dfe_train = 0.005; + + options.mu_dc_dd = 0.01; + options.mu_ffe_dd = [0.0004 0.0005 0.0006]; + options.mu_dfe_dd = 0.0005; + + options.eq_blocklength = 1; + options.eq_updatelatency = 1; + end + + fn = fieldnames(options); + for n = 1:numel(fn) + obj.(fn{n}) = options.(fn{n}); + end + + % Generate helpful vectors and initialize the filters with + % correct length: + + obj.Ce = obj.calcVNLEMemoryLength(obj.Ne); + + [obj.Ie2,obj.Ie3] = obj.calcIndiceVectors(obj.Ne); + + obj.e = zeros(sum(obj.Ce),1); + + + obj.Cb = obj.calcVNLEMemoryLength(obj.Nb); + + [obj.Ib2,obj.Ib3] = obj.calcIndiceVectors(obj.Nb); + + obj.b = zeros(sum(obj.Cb),1); + + end + + function [signalclass_out] = process(obj,signalclass_in, reference_signalclass_in) + + % actual processing of the signal (steps 1. - 3.) + % 1 normalize RMS + %signalclass_in = signalclass_in.normalize("mode","rms"); + + % Process the EQ optimization + obj.process_(signalclass_in.signal', reference_signalclass_in.signal'); + + signalclass_in.signal = obj.y_out'; + % append to logbook + lbdesc = ['EQ von Silas ist gelaufen ']; + signalclass_in = signalclass_in.logbookentry(lbdesc); + + % write to output + signalclass_out = signalclass_in; + + end + + function process_(obj,x_in,d_in) + + % 1) prepare signals + obj.e_dc = mean(x_in); + + % 1.1) Input Signal + obj.x_in = [zeros(1,floor(obj.Ne(1)/2)) x_in zeros(1,obj.Ne(1))]; + obj.x_length = length(x_in); + obj.x_norm = obj.calcPowerNormalization(x_in); + + % 1.2 Reference Signal // Constellation + obj.d = [zeros(1,obj.Nb(1)-1) d_in zeros(1,obj.Nb(1))]; + obj.d_constellation = unique(d_in); + obj.d_norm = obj.calcPowerNormalization(d_in); + + % 1.3 Training + obj.trainingMode(); + + % 1.4 Decision Directed Mode + obj.decisionDirectedMode(); + + end + + %% Adaptive Equalization Modes + + function trainingMode(obj) + + dc_avg_block = zeros(obj.eq_blocklength,1); + + for tloop = 1:obj.trainloops + m = 1+obj.delay; + + for n = obj.sps*obj.delay+1:obj.sps:obj.sps*obj.trainlength + m = m+1; + + %get Sigal input vectors with correct length for VNLE + x_in_block = obj.x_in(obj.Ne(1)+n+(obj.sps-1):-1:n+obj.sps).'; + + x_in_vnle_format = obj.calcVNLENonlinVecs(x_in_block,obj.Ie2,obj.Ie3,obj.Ne,obj.x_norm); + + %get Reference input vectors with correct length for VNLE + d_block = obj.d(obj.Nb(1)-obj.delay+m-2:-1:m-obj.delay-1).'; + d_vnle_format = obj.calcVNLENonlinVecs(d_block,obj.Ib2,obj.Ib3,obj.Nb,obj.d_norm); + + obj.e_ffe = obj.e.' * x_in_vnle_format; + + dc_avg_block = circshift(dc_avg_block,1); + + dc_avg_block(1) = obj.e_ffe; + + if n > obj.eq_blocklength * obj.sps + obj.e_dc = obj.e_ffe - mean(dc_avg_block).* obj.mu_dc_train; + else + obj.e_dc = 0; + end + + obj.e_dfe = obj.b.' * d_vnle_format; + + % Calculate the Error + obj.error = obj.e_ffe - obj.e_dfe - obj.d(obj.Nb(1)-1+m-obj.delay); + + if obj.mu_ffe_train ~= 0 + %update FFE coefficients with LMS + obj.e = obj.e - obj.error*conj(x_in_vnle_format)*obj.mu_ffe_train; + else + %update FFE coefficients with NLMS + obj.e = obj.e - obj.error*x_in_vnle_format/(x_in_vnle_format.'*x_in_vnle_format); + end + + %update DFE coefficients with LMS + obj.b = obj.b + obj.mu_dfe_train*obj.error*d_vnle_format; + + end + + end + + end + + + function decisionDirectedMode(obj) + + %start the dd mode with coefficients from training + coeff = [obj.e;obj.b]; + + + for ddloop = 1:obj.ddloops + + dc_avg_block = zeros(obj.eq_blocklength,1); + m = 0; + dc_cnt = 0; + + mu_mat = diag([ones(1,obj.Ce(1))*obj.mu_ffe_dd(1)... %1st order ffe + ones(1,obj.Ce(2))*obj.mu_ffe_dd(2)... %2nd order ffe + ones(1,obj.Ce(3))*obj.mu_ffe_dd(3)... %3rd order ffe + ones(1,sum(obj.Cb))*obj.mu_dfe_dd]); %all order dfe + + mu_ffe = [ones(1,obj.Ce(1))*obj.mu_ffe_dd(1)... %1st order ffe + ones(1,obj.Ce(2))*obj.mu_ffe_dd(2)... %2nd order ffe + ones(1,obj.Ce(3))*obj.mu_ffe_dd(3)]; + + mu_dfe = ones(1,sum(obj.Cb))*obj.mu_dfe_dd; + + y = zeros(1,floor(obj.x_length/obj.sps)); + y_= zeros(1,floor(obj.x_length/obj.sps)); + d_feedback = zeros(obj.Cb(1),1); + d_vnle = obj.calcVNLENonlinVecs(d_feedback,obj.Ib2,obj.Ib3,obj.Nb,obj.d_norm); + d_hat = zeros(obj.x_length,1); + + for k = 1:obj.sps:obj.x_length + dc_cnt = dc_cnt+1; + m=m+1; + + %get Sigal input vectors with correct length for VNLE + x = obj.x_in(obj.Ne(1)+k-1:-1:k).'; + + x_vnle = obj.calcVNLENonlinVecs(x,obj.Ie2,obj.Ie3,obj.Ne,obj.x_norm); + + %combine FFE with DFE to one vector (cursor between the two sequences) + x_d = [x_vnle;-d_vnle]; + + %Apply filter + y_(m) = x_d.'* coeff ; + y(m) = y_(m) - obj.e_dc(end) ; + + dc_avg_block = circshift(dc_avg_block,1); + + dc_avg_block(1) = y(m); + + obj.e_dc(m) = (y(m) - mean(dc_avg_block)) .* obj.mu_dc_dd; + + %Decision + [~,symbol_idx] = min(abs(y(m) - obj.d_constellation)); % decision for closest constellation point + + d_hat(k) = obj.d_constellation(symbol_idx); + + %Error between FFE & DFE filtered signal and Decision + obj.error = y(m) - d_hat(k); + + %Update coefficients (both FFE and DFE) + % obj.e = obj.e - obj.error * mu_ffe * conj(x_vnle); + % obj.b = obj.b + obj.error * mu_dfe * conj(d_vnle); + % coeff = [obj.e;obj.b]; + + coeff = coeff - mu_mat*obj.error*conj(x_d); + + % Append new decision to decision feedback + if obj.Nb(1) > 0 + + %shift up one index + d_feedback(2:end) = d_feedback(1:end-1); + %replace 1st index with current estimation + d_feedback(1) = d_hat(k); + %build memorylike VNLE version + d_vnle = obj.calcVNLENonlinVecs(d_feedback,obj.Ib2,obj.Ib3,obj.Nb,obj.d_norm); + + end + + end + end + + obj.y_out = (circshift( y.' ,-(obj.delay))).'; + + end + + %% Functions needed During Adaption + function x_in_vnle_format = calcVNLENonlinVecs(~,x_in_block,I_2,I_3,N_,norm_) + % These are the second and third order input signal products of the VNLE EQ + % ∑ h1 x_in(k-n1) + ∑∑ h2 x_in(k-n1)*x_in(k-n2) + ∑∑∑ h3 x_in(k-n1)*x_in(k-n2)*x_in(k-n3) + + x1 = x_in_block; + x2 = []; + x3 = []; + + if N_(2) > 0 + delta_2 = round((N_(1)-N_(2))/2); + input_vec_se = x_in_block(delta_2:end)/norm_(2); %TODO normalization step + x2 = input_vec_se(I_2(:,1)).*input_vec_se(I_2(:,2)); + end + + if N_(3) > 0 + delta_3 = round((N_(1)-N_(3))/2); + input_vec_th = x_in_block(delta_3:end)/norm_(3); + x3 = input_vec_th(I_3(:,1)).*input_vec_th(I_3(:,2)).*input_vec_th(I_3(:,3)); + end + + x_in_vnle_format = [x1;x2;x3]; + + end + + + %% Functions needed for Preparation + function [C] = calcVNLEMemoryLength(~,N) + + %calculates the memory length of VNLE + C = zeros(size(N)); + + for o = 1:numel(N) + switch o + case 1 + C(o) = N(o); + case 2 + C(o) = N(o)*(N(o)+1) / 2; + case 3 + C(o) = N(o)*(N(o)+1)*(N(o)+2) / 6; + end + end + + end + + function [indvec2nd, indvec3rd] = calcIndiceVectors(~,N) + + % Init vectors of 2nd and 3rd order coefficient indices -> + % yield combination with + + for order = 2:numel(N) + n = N(order); + v = 1:n; % Ursprünglicher Vektor + row = 1; + + % Schleifen zur Generierung des Indize Vektors + switch order + + case 2 + + indvec2nd = zeros(n*(n+1)/2, order); + for i = 1:n + for j = i:n + indvec2nd(row, :) = [v(i) v(j)]; + row = row + 1; + end + end + + case 3 + + indvec3rd = zeros(n*(n+1)*(n+2)/6, 3); + for i = 1:n + for j = i:n + for k = j:n + indvec3rd(row, :) = [v(i) v(j) v(k)]; + row = row + 1; + end + end + end + end + end + + end + + function powerNorm = calcPowerNormalization(~,v) + + powerNorm(1) = sqrt(mean(abs(v ).^2)); + powerNorm(2) = sqrt(mean(abs(v.^2).^2)); + powerNorm(3) = sqrt(mean(abs(v.^3).^2)); + + end + + + end +end + + + + + + + + + + + + + + + + + + + + + + diff --git a/Functions/calc_ber.m b/Functions/calc_ber.m index 6ed7bc3..046cc0c 100644 --- a/Functions/calc_ber.m +++ b/Functions/calc_ber.m @@ -1,37 +1,58 @@ -function [bits,errors,ber,errorLocation] = calc_ber(data_in,data_ref,options) +function [bits,errors,ber,errorIndice] = calc_ber(data_in,data_ref,options) arguments(Input) data_in; data_ref; - options.skip = 0; + options.skip_front = 0; + options.skip_end = 0; options.returnErrorLocation = 0; end +options.skip_end = abs(options.skip_end); +options.skip_front = abs(options.skip_front); bits = 0; errors= 0; ber= 0; -errorLocation= []; +errorIndice= []; -data_ref = logical(data_ref)'; -data_in = logical(data_in)'; +% trim sequence to given start; stop. trim reference if signal is shorter +[data_in,data_ref]=trimseq(data_in,data_ref,options.skip_front,options.skip_end); -bits = 0; -bits = bits+numel(data_in(:,options.skip+1:end)); +if length(data_ref) == length(data_in) -try - if options.returnErrorLocation == 0 - errors = sum( data_in(:,options.skip+1:end) ~= data_ref(:,options.skip+1:end),"all" ); - else - errorLocation = sum(data_in(:,options.skip+1:end) ~= data_ref(:,options.skip+1:end),1); - errors = sum(errorLocation ,"all" ); + bits = 0; + bits = bits+numel(data_in(:,options.skip_front+1:end)); + + try + if options.returnErrorLocation == 0 + errors = sum( data_in ~= data_ref,"all" ); + else + errorIndice = sum(data_in ~= data_ref,1); + errors = sum(errorIndice ,"all" ); + [~,errorIndice] = find(errorIndice==1); + end + + catch + warning('BER calculation not optimal: Arrays have incompatible sizes for this operation.') + errors = NaN; end -catch - %warning('BER calculation not optimal: Arrays have incompatible sizes for this operation.') - errors = NaN; + % Determine BER + ber = sum(errors)/sum(bits); + +else + errormsg('Sequence length does not match'); end -% Determine BER -ber = sum(errors)/sum(bits); + function [data_,reference_]=trimseq(data,reference,skipstart,skip_end) + + data_ = logical(data(skipstart+1:end-skip_end,:))'; + + delta_bits = length(reference) - length(data_); + + skip_end = max(skip_end,delta_bits); + reference_ = logical(reference(skipstart+1:end-skip_end,:))'; + + end end \ No newline at end of file diff --git a/MPI_offline_analysis/mpi_analysis_pam4.m b/MPI_offline_analysis/mpi_analysis_pam4.m deleted file mode 100644 index dabfd70..0000000 --- a/MPI_offline_analysis/mpi_analysis_pam4.m +++ /dev/null @@ -1,120 +0,0 @@ - -clear; -col = linspecer(6); - -if ismac - foldername = '/Users/silasoettinghaus/Documents/MATLAB/Labor_Datensatz_PAM4_MPI/pam4_10km_1km'; -else - foldername = 'C:\Users\Silas\Documents\MATLAB\Labor_Datensatz_PAM4_MPI\pam4_10km_1km'; -end - -allfiles = dir(foldername); - -mudc_loop = [0.05 0.01 0.005]; - -parfsym = zeros(length(mudc_loop),length(allfiles)); -parsir = zeros(length(mudc_loop),length(allfiles)); -ber_wo_reduction = zeros(length(mudc_loop),length(allfiles)); -ber_w_dcremoval = zeros(length(mudc_loop),length(allfiles)); -mudc = zeros(length(mudc_loop),length(allfiles)); - - -for j = 1:length(mudc_loop) - mudc_current = mudc_loop(j); - parfor i = 1:length(allfiles) - - if allfiles(i).bytes ~= 0 - current_filename = allfiles(i).name; - else - continue - end - - % load dataset - simdata = load([foldername,filesep, current_filename]); - % get common variables of simulation loop - - %even if redundant, save stuff to results struct - parsir(j,i) = simdata.saveStructTemp.awg2scope_keysight_state.eigenlight_mpi - 3.3 + 7.2; - parfsym(j,i) = simdata.saveStructTemp.common.f_sym; - - fdac = simdata.saveStructTemp.common.f_DAC; - fadc = simdata.saveStructTemp.common.f_ADC; - fsym = simdata.saveStructTemp.common.f_sym; - - % 0) build RX Signal - - y_rx = Electricalsignal(simdata.saveStructTemp.dp_tsynch_out'); - - y_rx.fs = 2.*fsym; - - % 0) build tx reference Signal for eq training - y_digimod = Electricalsignal(simdata.saveStructTemp.digi_mod_out'); - y_digimod.fs = fsym; - - - % 1) normlaize - y_rx = y_rx.normalize("mode","rms"); - - eq = EQ_silas("Ne",[50,9,9],"Nb",[2,0,0],"trainlength",4096,... - "sps",2,... - "mu_dc_dd",mudc_current,... - "mu_dc_train",mudc_current,... - "mu_ffe_train",0,... - "mu_dfe_train",0.005,... - "mu_ffe_dd",[0.0004 0.0004 0.0004],... - "mu_dfe_dd",0.005,... - "ddloops",4,... - "trainloops",4,... - "eq_blocklength",10); - - [Eq_out] = eq.process(y_rx,y_digimod); - - - % 2.2) MPI reduction blocks - - % 3) digital demodulation - digimod = PAMmapper(2^simdata.saveStructTemp.common.M,0); - - d_hat = digimod.demap(Eq_out); - - d_ = digimod.demap(Electricalsignal(simdata.saveStructTemp.digi_mod_out')); - d = simdata.saveStructTemp.prms_out; - - % 4) BER calculation - [totalbits,errors_bm,ber_dhat,loc] = calc_ber(d_hat.signal(1:end,:) ,d(1:end,:)',"skip",0,"returnErrorLocation",1); - - ber_nodctap = mean(simdata.saveStructTemp.prms_compare_state.BER ); - - disp(['fsym: ',num2str(parfsym(j,i).*1e-9),'GBd, SIR: ',num2str(parsir(j,i)),' ->> BER_new: ',sprintf('%2E',ber_dhat),' BER_old: ',sprintf('%2E',ber_nodctap)]); - - % 5) save to struct - - ber_wo_reduction(j,i) = ber_nodctap; - ber_w_dcremoval(j,i) = ber_dhat; - mudc(j,i) = eq.mu_dc_dd; - - end -end - -rate = [92e9, 56e9]; -figure(99) - -for j = 1:length(mudc_loop) - for r = 1:length(rate) - a = [parsir(j,parfsym(j,:)==rate(r));ber_wo_reduction(j,parfsym(j,:)==rate(r))]; - a=sort(a,2); - plot(a(1,:),a(2,:),'DisplayName',[num2str(rate(r).*1e-9),' GBd; PAM4; w/o reduction'],'LineWidth',2,'Marker','o','MarkerSize',3); - %scatter(parsir(parfsym==rate(r)),ber_wo_reduction(parfsym==rate(r)),'DisplayName',[num2str(rate(r).*1e-9),' GBd; PAM4; w/o reduction'],'LineWidth',2,'Marker','o') - hold on - b = [parsir(j,parfsym(j,:)==rate(r));ber_w_dcremoval(j,parfsym(j,:)==rate(r))]; - b=sort(b,2); - plot(b(1,:),b(2,:),'DisplayName',[num2str(rate(r).*1e-9),' GBd; PAM4; adaptive DC removal'],'LineWidth',2,'Marker','+','MarkerSize',3); - %scatter(parsir(parfsym==rate(r)),ber_w_dcremoval(parfsym==rate(r)),'DisplayName',[num2str(rate(r).*1e-9),' GBd; PAM4; adaptive DC removal'],'LineWidth',2,'Marker','+'); - end -end -set(gca, 'YScale', 'log'); -yline(3.8e-3,'LineWidth',1.5); - - - - diff --git a/filter_debug.m b/projects/IMDD_base_system/filter_debug.m similarity index 94% rename from filter_debug.m rename to projects/IMDD_base_system/filter_debug.m index 9d095fc..ad0a4be 100644 --- a/filter_debug.m +++ b/projects/IMDD_base_system/filter_debug.m @@ -1,70 +1,70 @@ - -clear -%% Set Simulation Variables - -O = 18; %order of prbs -N = 2^(O-1); %length of prbs -[~,seed] = prbs(O,1); %initialize first seed of prbs - -% Modulation -M = 4; %PAM-M -bitpattern = zeros(N,log2(M)); - -% Symbol Rate -fsym = 112e9; - -% DAC Rate -fdac = 120e9; - -% Simulation oversampling rate "k"; -kover = 16; - -% ADC Rate -fadc = 256e9; - -% Simulation frequency in "analog domain" -fsimu = kover * fdac ; - - - -%% CONSTRUCT ALL CLASSES - -digimod = PAMmapper(M,0); - -pulsef = Pulseformer("pulseform","rrc","fdac",fdac,"fsym",fsym,"pulselength",32,"rrcalpha",0.05); - -awg = AWG('preset','M8199B','fdac',fdac,'kover',kover,'lpf_active',1,'f_cutoff',56e9,'lpf_type',filtertypes.gaussian,'bit_resolution',5.5); - -lp_laser = Filter('filtdegree',1,"f_cutoff",30e9,"fsamp",fdac,"filterType",filtertypes.bessel_inp); - - - -%% PROCESS - -% PRBS Generation -for i = 1:log2(M) - [bitpattern(:,i),seed] = prbs(O,N,seed); -end - -% Build Inf. signal class -bits = Informationsignal(bitpattern); - -% Digi Mod -mod_out = digimod.map(bits); -% merken für EQ training -reference = mod_out; - -% shape shape -mod_out = pulsef.process(mod_out); -test = applyPulseShaping(reference.signal,fsym,fdac); - -% AWG -> ELECTRICAL DOMAIN -X = lp_laser.process(mod_out); -%X.spectrum(fsimu,"displayname",'AWG out','figurename','after AWG'); - - - - - - - + +clear +%% Set Simulation Variables + +O = 18; %order of prbs +N = 2^(O-1); %length of prbs +[~,seed] = prbs(O,1); %initialize first seed of prbs + +% Modulation +M = 4; %PAM-M +bitpattern = zeros(N,log2(M)); + +% Symbol Rate +fsym = 112e9; + +% DAC Rate +fdac = 120e9; + +% Simulation oversampling rate "k"; +kover = 16; + +% ADC Rate +fadc = 256e9; + +% Simulation frequency in "analog domain" +fsimu = kover * fdac ; + + + +%% CONSTRUCT ALL CLASSES + +digimod = PAMmapper(M,0); + +pulsef = Pulseformer("pulseform","rrc","fdac",fdac,"fsym",fsym,"pulselength",32,"rrcalpha",0.05); + +awg = AWG('preset','M8199B','fdac',fdac,'kover',kover,'lpf_active',1,'f_cutoff',56e9,'lpf_type',filtertypes.gaussian,'bit_resolution',5.5); + +lp_laser = Filter('filtdegree',1,"f_cutoff",30e9,"fsamp",fdac,"filterType",filtertypes.bessel_inp); + + + +%% PROCESS + +% PRBS Generation +for i = 1:log2(M) + [bitpattern(:,i),seed] = prbs(O,N,seed); +end + +% Build Inf. signal class +bits = Informationsignal(bitpattern); + +% Digi Mod +mod_out = digimod.map(bits); +% merken für EQ training +reference = mod_out; + +% shape shape +mod_out = pulsef.process(mod_out); +test = applyPulseShaping(reference.signal,fsym,fdac); + +% AWG -> ELECTRICAL DOMAIN +X = lp_laser.process(mod_out); +%X.spectrum(fsimu,"displayname",'AWG out','figurename','after AWG'); + + + + + + + diff --git a/imddmodel.m b/projects/IMDD_base_system/imddmodel.m similarity index 96% rename from imddmodel.m rename to projects/IMDD_base_system/imddmodel.m index 7510896..b3879bf 100644 --- a/imddmodel.m +++ b/projects/IMDD_base_system/imddmodel.m @@ -1,226 +1,227 @@ -function output = imddmodel(sir,dcmode) - -for realiz = 1:3 - - rng(realiz); - - %% Set Simulation Variables - delay = 0;10+(10*realiz); %mpi delay in meter - - fiblen = 0; %main link in km - - laser_linewidth =0e6; - - O = 18; %order of prbs - N = 2^(O-1); %length of prbs - [~,seed] = prbs(O,1); %initialize first seed of prbs - - % Modulation - M = 4; %PAM-M - bitpattern = zeros(N,log2(M)); - - % Symbol Rate - fsym = 112e9; - - % DAC Rate - fdac = 120e9; - - % Simulation oversampling rate "k"; - kover = 16; - - % ADC Rate - fadc = 256e9; - - % Simulation frequency in "analog domain" - fsimu = kover * fdac ; - - - - %% B) CONSTRUCT ALL CLASSES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - digimod = PAMmapper(M,0); - - pulseform = Pulseformer("pulseform","rrc","fdac",fdac,"fsym",fsym,"pulselength",32,"rrcalpha",0.027); - - awg = AWG('fdac',fdac,'kover',kover,'lpf_active',1,'f_cutoff',56e9,'lpf_type',filtertypes.gaussian,'bit_resolution',5.5); - - lp_laser = Filter('filtdegree',1,"f_cutoff",50e9,"fsamp",fdac*kover,"filterType",filtertypes.bessel_inp); - - u_pi = 3.5; - vbias = (0.5*u_pi)-u_pi; - extmodlaser = EML("mode",eml_mode.im_cosinus,"power",5,"fsimu",fsimu,"lambda",1550,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth); - - fib = Fiber("fsimu",fdac*kover,"fiber_length",fiblen,"alpha",0.2,"D",16,"lambda0",thz2nm(193.1),"gamma",0); - - reflectionpoint = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",sir); - reflectionprop = Fiber("fsimu",fdac*kover,"fiber_length",delay/1000,"alpha",0.2,"D",16,"lambda0",1550,"gamma",0); - - opticatten = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",-8); - - phdiode = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20); - - lp_diode = Filter('filtdegree',1,"f_cutoff",50e9,"fsamp",fdac*kover,"filterType",filtertypes.bessel_inp); - - scp = Scope("fsimu",fdac*kover,"fadc",fadc,... - "delay",0,"fixed_delay",0,"lpf_bw",113e9,"filtertype",filtertypes.butterworth,... - "samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,... - "adcresolution",6,"quantbuffer",0.1,'block_dc',1); - - eq = EQ("K",2,"plottrain",0,"plotfinal",0,... - "training_length",4096,"training_loops",2,... - "Ne",[50,5,3],"Nb",[3,2,2],... - "DCmu",0.05,"DDmu",[0.0004 0.0005 0.0006 0.0007 ],"DFEmu",0.005,"FFEmu",0.005,... - "dd_loops",2,"epsilon",[10 100 1000 ],"M",2,... - "thres",[0.005 0.004 0.0005 ],"l1act",0,"delay",0,"rho",0.0005,"ideal_dfe",0,"DB_aim",0); - - eq2 = EQ_silas("Ne",[20,0,0],"Nb",[3,0,0],"trainlength",4096,"mu_dc_dd",0.005,"mu_dc_train",0.05,... - "mu_ffe_train",0.005,"mu_combined_dd",[0.0004 0.0006 0.0003 0.005],"ddloops",3,"dcmode",dcmode); - - - %% C) PROCESS TX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - % 1) PRBS Generation - for i = 1:log2(M) - [bitpattern(:,i),seed] = prbs(O,N,seed); - end - - % 2 ) Build Inf. signal class - bits = Informationsignal(bitpattern); - - % 3) Digi modulation -> PAM-M signal - digimod_out = digimod.map(bits); - - % 4) Pulse shaping -> racos - X = pulseform.process(digimod_out); - - % 5) AWG (lowpass, quantization, sample and hold) - X = awg.process(X); - - % 6) Lowpass behavior of laser and hf-cable? why twice? - X = lp_laser.process(X); - X = lp_laser.process(X); - - % 7) Normalize signal - X = X.normalize("mode","oneone"); - X.signal = X.signal .* 1.3800; - - - - %% D) PROCESS OPTICAL CHANNEL %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - % 1) Laser; Modulation -> OPTICAL DOMAIN - [X,extmodlaser] = extmodlaser.process(X); - - % 2) Main fiber propagation - X = fib.process(X); - - if delay ~= 0 - - % 3) Reflection - % Reflection is just an attenuation - R = reflectionpoint.process(X); - - % Propagate back and forth (actual fiber propagation) - R = reflectionprop.process(R); - - % Delay the reflected signal - [R,n] = R.delay("delay_meter",delay); - - % Add together - X.signal = X.signal(n:end); - R.signal = R.signal(n:end); - - %disp(['SIR ',num2str(10*log10(X.power/R.power))]); - - X = X+R; - - % 4) Equalize - % cut reference signal to correct length (nessecary due to MPI delay) - digimod_out.signal = digimod_out.signal(round(n * fsym/fsimu) : end,:); - bitpattern = bitpattern(round(n * fsym/fsimu):end,:); - end - - % 4) Attenuation - X = opticatten.process(X); - % X = edfaamp.process(X); - - % 5) Photo Diode -> ELECTRICAL DOMAIN - X = phdiode.process(X); - X = lp_diode.process(X); - - - %% E) PROCESS RX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - % 1) Oscilloscope (Sampling to f_adc; Quantization; Bandwidth Limitation) - X = scp.process(X); - - % 2) Resample to 2x symbol rate - X = X.resample("fs_out",2*fsym,"fs_in",fadc); - - % 3) Normalize - Eq_in = X.normalize("mode","rms"); - - - - % MPI reduction DC removal BEFORE EQ - wl = 3000; % symbols - Eq_in.signal = Eq_in.signal - 1/wl .* movsum( Eq_in.signal,[wl/2,wl/2]); - - % Equalize Signal - [Eq_out] = eq2.process(Eq_in,digimod_out); - - - - %% A1: MPI reduction DC removal - wl = 1000; % symbols - yk_dcsm = Eq_out; - yk_dcsm.signal = Eq_out.signal - 1/wl .* movsum( Eq_out.signal,[wl/2,wl/2]); - - %% A2: MPI reduction Level wise error removal - yk_lvsm = Eq_out; - yk_lvlp = Eq_out; - pre_decision_level_uni = digimod.decide_pamlevel(Eq_out); - pre_decision_level_bi = ( pre_decision_level_uni*2-3 ) .* 1/sqrt(5); - - e = Eq_out.signal - pre_decision_level_bi; - - lp_mpi = Filter('filtdegree',1,"f_cutoff",2e6,"fsamp",fsym,"filterType",filtertypes.bessel_inp); - - filtered = lp_mpi.process(e); - - wl = 30; % symbols - smoothed = ( 1/wl .* movsum(e,[wl/2,wl/2]) ); - - % remove interference - for level = 0:3 - yk_lvsm.signal(pre_decision_level_uni==level) = yk_lvsm.signal(pre_decision_level_uni==level) - smoothed(pre_decision_level_uni==level); - yk_lvlp.signal(pre_decision_level_uni==level) = yk_lvlp.signal(pre_decision_level_uni==level) - filtered(pre_decision_level_uni==level); - end - - %% PROCESS RX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - % 1) Digi Demod - d_bm = digimod.demap(Eq_out); - d_dcsm = digimod.demap(yk_dcsm); - d_lvsm = digimod.demap(yk_lvsm); - d_lvlp = digimod.demap(yk_lvlp); - - % 2) BER - - dbit = length(d_bm.signal)-length(bitpattern); - - [~,errors_bm,ber_bm(realiz),loc] = calc_ber(d_bm.signal(1:end-dbit,:) ,bitpattern(1:end,:),"skip",0,"returnErrorLocation",1); - [~,errors_dcsm,ber_dcsm(realiz)] = calc_ber(d_dcsm.signal(1:end-dbit,:) ,bitpattern(1:end,:),"skip",0,"returnErrorLocation",0); - [~,errors_lvsm,ber_lvsm(realiz)] = calc_ber(d_lvsm.signal(1:end-dbit,:) ,bitpattern(1:end,:),"skip",0,"returnErrorLocation",0); - [~,errors_lvlp,ber_lvlp(realiz)] = calc_ber(d_lvlp.signal(1:end-dbit,:) ,bitpattern(1:end,:),"skip",0,"returnErrorLocation",0); - - -end - -output.ber_bm = mean(ber_bm); -output.ber_dcsm = mean(ber_dcsm); -output.ber_lvsm = mean(ber_lvsm); -output.ber_lvlp = mean(ber_lvlp); - - +function output = imddmodel(sir,dcmode) + +for realiz = 1:3 + + rng(realiz); + + %% Set Simulation Variables + delay = 0;10+(10*realiz); %mpi delay in meter + + fiblen = 0; %main link in km + + laser_linewidth =0e6; + + O = 18; %order of prbs + N = 2^(O-1); %length of prbs + [~,seed] = prbs(O,1); %initialize first seed of prbs + + % Modulation + M = 4; %PAM-M + bitpattern = zeros(N,log2(M)); + + % Symbol Rate + fsym = 112e9; + + % DAC Rate + fdac = 120e9; + + % Simulation oversampling rate "k"; + kover = 16; + + % ADC Rate + fadc = 256e9; + + % Simulation frequency in "analog domain" + fsimu = kover * fdac ; + + + + %% B) CONSTRUCT ALL CLASSES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + digimod = PAMmapper(M,0); + + pulseform = Pulseformer("pulseform","rrc","fdac",fdac,"fsym",fsym,"pulselength",32,"rrcalpha",0.027); + + awg = AWG('fdac',fdac,'kover',kover,'lpf_active',1,'f_cutoff',56e9,'lpf_type',filtertypes.gaussian,'bit_resolution',5.5); + + lp_laser = Filter('filtdegree',1,"f_cutoff",50e9,"fsamp",fdac*kover,"filterType",filtertypes.bessel_inp); + + u_pi = 3.5; + vbias = (0.5*u_pi)-u_pi; + extmodlaser = EML("mode",eml_mode.im_cosinus,"power",5,"fsimu",fsimu,"lambda",1550,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth); + + fib = Fiber("fsimu",fdac*kover,"fiber_length",fiblen,"alpha",0.2,"D",16,"lambda0",thz2nm(193.1),"gamma",0); + + reflectionpoint = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",sir); + reflectionprop = Fiber("fsimu",fdac*kover,"fiber_length",delay/1000,"alpha",0.2,"D",16,"lambda0",1550,"gamma",0); + + opticatten = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",-8); + + phdiode = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20); + + lp_diode = Filter('filtdegree',1,"f_cutoff",50e9,"fsamp",fdac*kover,"filterType",filtertypes.bessel_inp); + + scp = Scope("fsimu",fdac*kover,"fadc",fadc,... + "delay",0,"fixed_delay",0,"lpf_bw",113e9,"filtertype",filtertypes.butterworth,... + "samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,... + "adcresolution",6,"quantbuffer",0.1,'block_dc',1); + + eq = EQ("K",2,"plottrain",0,"plotfinal",0,... + "training_length",4096,"training_loops",2,... + "Ne",[50,5,3],"Nb",[3,2,2],... + "DCmu",0.05,"DDmu",[0.0004 0.0005 0.0006 0.0007 ],"DFEmu",0.005,"FFEmu",0.005,... + "dd_loops",2,"epsilon",[10 100 1000 ],"M",2,... + "thres",[0.005 0.004 0.0005 ],"l1act",0,"delay",0,"rho",0.0005,"ideal_dfe",0,"DB_aim",0); + + eq2 = EQ_silas("Ne",[20,0,0],"Nb",[3,0,0],"trainlength",4096,"mu_dc_dd",0.005,"mu_dc_train",0.05,... + "mu_ffe_train",0.005,"mu_combined_dd",[0.0004 0.0006 0.0003 0.005],"ddloops",3,"dcmode",dcmode); + + + %% C) PROCESS TX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % 1) PRBS Generation + for i = 1:log2(M) + [bitpattern(:,i),seed] = prbs(O,N,seed); + end + + % 2 ) Build Inf. signal class + bits = Informationsignal(bitpattern); + + % 3) Digi modulation -> PAM-M signal + digimod_out = digimod.map(bits); + + % 4) Pulse shaping -> racos + X = pulseform.process(digimod_out); + + % 5) AWG (lowpass, quantization, sample and hold) + X = awg.process(X); + + % 6) Lowpass behavior of laser and hf-cable? why twice? + X = lp_laser.process(X); + X = lp_laser.process(X); + + % 7) Normalize signal + X = X.normalize("mode","oneone"); + X.signal = X.signal .* 1.3800; + + + + %% D) PROCESS OPTICAL CHANNEL %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % 1) Laser; Modulation -> OPTICAL DOMAIN + [X,extmodlaser] = extmodlaser.process(X); + + % 2) Main fiber propagation + X = fib.process(X); + + if delay ~= 0 + + % 3) Reflection + % Reflection is just an attenuation + R = reflectionpoint.process(X); + + % Propagate back and forth (actual fiber propagation) + R = reflectionprop.process(R); + + % Delay the reflected signal + [R,n] = R.delay("delay_meter",delay); + + % Add together + X.signal = X.signal(n:end); + R.signal = R.signal(n:end); + + %disp(['SIR ',num2str(10*log10(X.power/R.power))]); + + X = X+R; + + % 4) Equalize + % cut reference signal to correct length (nessecary due to MPI delay) + digimod_out.signal = digimod_out.signal(round(n * fsym/fsimu) : end,:); + bitpattern = bitpattern(round(n * fsym/fsimu):end,:); + end + + % 4) Attenuation + X = opticatten.process(X); + % X = edfaamp.process(X); + + % 5) Photo Diode -> ELECTRICAL DOMAIN + X = phdiode.process(X); + X = lp_diode.process(X); + + + %% E) PROCESS RX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % 1) Oscilloscope (Sampling to f_adc; Quantization; Bandwidth Limitation) + X = scp.process(X); + + % 2) Resample to 2x symbol rate + X = X.resample("fs_out",2*fsym,"fs_in",fadc); + + % 3) Normalize + Eq_in = X.normalize("mode","rms"); + + + + % MPI reduction DC removal BEFORE EQ + wl = 3000; % symbols + Eq_in.signal = Eq_in.signal - 1/wl .* movsum( Eq_in.signal,[wl/2,wl/2]); + + % Equalize Signal + [Eq_out] = eq2.process(Eq_in,digimod_out); + + + + + %% A1: MPI reduction DC removal + wl = 1000; % symbols + yk_dcsm = Eq_out; + yk_dcsm.signal = Eq_out.signal - 1/wl .* movsum( Eq_out.signal,[wl/2,wl/2]); + + %% A2: MPI reduction Level wise error removal + yk_lvsm = Eq_out; + yk_lvlp = Eq_out; + pre_decision_level_uni = digimod.decide_pamlevel(Eq_out); + pre_decision_level_bi = ( pre_decision_level_uni*2-3 ) .* 1/sqrt(5); + + e = Eq_out.signal - pre_decision_level_bi; + + lp_mpi = Filter('filtdegree',1,"f_cutoff",2e6,"fsamp",fsym,"filterType",filtertypes.bessel_inp); + + filtered = lp_mpi.process(e); + + wl = 30; % symbols + smoothed = ( 1/wl .* movsum(e,[wl/2,wl/2]) ); + + % remove interference + for level = 0:3 + yk_lvsm.signal(pre_decision_level_uni==level) = yk_lvsm.signal(pre_decision_level_uni==level) - smoothed(pre_decision_level_uni==level); + yk_lvlp.signal(pre_decision_level_uni==level) = yk_lvlp.signal(pre_decision_level_uni==level) - filtered(pre_decision_level_uni==level); + end + + %% PROCESS RX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % 1) Digi Demod + d_bm = digimod.demap(Eq_out); + d_dcsm = digimod.demap(yk_dcsm); + d_lvsm = digimod.demap(yk_lvsm); + d_lvlp = digimod.demap(yk_lvlp); + + % 2) BER + + dbit = length(d_bm.signal)-length(bitpattern); + + [~,errors_bm,ber_bm(realiz),loc] = calc_ber(d_bm.signal(1:end-dbit,:) ,bitpattern(1:end,:),"skip",0,"returnErrorLocation",1); + [~,errors_dcsm,ber_dcsm(realiz)] = calc_ber(d_dcsm.signal(1:end-dbit,:) ,bitpattern(1:end,:),"skip",0,"returnErrorLocation",0); + [~,errors_lvsm,ber_lvsm(realiz)] = calc_ber(d_lvsm.signal(1:end-dbit,:) ,bitpattern(1:end,:),"skip",0,"returnErrorLocation",0); + [~,errors_lvlp,ber_lvlp(realiz)] = calc_ber(d_lvlp.signal(1:end-dbit,:) ,bitpattern(1:end,:),"skip",0,"returnErrorLocation",0); + + +end + +output.ber_bm = mean(ber_bm); +output.ber_dcsm = mean(ber_dcsm); +output.ber_lvsm = mean(ber_lvsm); +output.ber_lvlp = mean(ber_lvlp); + + end \ No newline at end of file diff --git a/run_loop.m b/projects/IMDD_base_system/run_loop.m similarity index 96% rename from run_loop.m rename to projects/IMDD_base_system/run_loop.m index 6fe1b56..2041d50 100644 --- a/run_loop.m +++ b/projects/IMDD_base_system/run_loop.m @@ -1,208 +1,208 @@ - -clear - -sir_loop = [-25:2:-13]; -sir_loop = 0; -lw_loop = [1,2,3]; - -data = cell(length(sir_loop),length(lw_loop)); - -iterations=size(data); - -parfor ix = 1:numel(data) - - [u1,u2] = ind2sub(iterations,ix); - - output = imddmodel(sir_loop(u1),lw_loop(u2)); - - data{ix} = output; - -end - -for sir = 1:size(data,1) - for lw = 1:size(data,2) - ber_bm(sir,lw) = data{sir,lw}.ber_bm; - ber_dcsm(sir,lw) = data{sir,lw}.ber_dcsm; - ber_lvsm(sir,lw) = data{sir,lw}.ber_lvsm; - ber_lvlp(sir,lw) = data{sir,lw}.ber_lvlp; - end -end - -col = [ 0.6510 0.8078 0.8902 - 0.6980 0.8745 0.5412 - 0.9922 0.7490 0.4353]; -figure(22) -for l = 1:numel(lw_loop) - hold on - plot(sir_loop,ber_bm(:,l),'DisplayName',['DC mode:', num2str(lw_loop(l)),', Linewidth= 50 MHz'],'Marker','o','MarkerFaceColor',col(l,:),'Color',col(l,:),'LineWidth',2,'LineStyle','--'); -end -set(gca,'yscale','log'); -yline(3.8e-3,'LineWidth',2,'LineStyle','--','HandleVisibility','off'); -set(gca,'xdir','reverse'); - - - -%% Plot -col = linspecer(8); -figure(1) -hold on -m = ["x","o","pentagram","hexagram","*","+"]; -cnt = 1; -for d = 1%:size(data,2) - for dc = 1:size(data,1) - for lw = 1:size(data,3) - ber_bm(dc,d,lw) = (data{dc,d,lw}.ber_bm); - ber_dcsm(dc,d,lw) = data{dc,d,lw}.ber_dcsm; - ber_lvsm(dc,d,lw) = data{dc,d,lw}.ber_lvsm; - ber_lvlp(dc,d,lw) = (data{dc,d,lw}.ber_lvlp); - end - end - i = 1; - comm_dn = [];% ['EQ DC Tap: ',num2str(dc_tap_loop(d)),' m']; - title("Dependency on Laser Linewidth; B2B; Delay : 2*50m") - plot(lw_loop*1e-6,mean(squeeze(ber_bm(1:end,d,:)),1),"LineWidth",1.2,"Marker",m(1),"MarkerSize",5,'Color',col(cnt,:),'DisplayName',[' ',comm_dn]); - plot(lw_loop*1e-6,mean(squeeze(ber_dcsm(1:end,d,:)),1),"LineWidth",1,"Marker",m(1+4),"MarkerSize",5,'LineStyle','--','Color',col(cnt,:),'DisplayName',['DC smoothing ',comm_dn],'HandleVisibility','on'); - plot(lw_loop*1e-6,mean(squeeze(ber_lvsm(1:end,d,:)),1),"LineWidth",1.2,"Marker",m(1+1),"MarkerSize",5,'LineStyle',':','Color',col(cnt,:),'DisplayName',['Lvl Smoothing ',comm_dn],'HandleVisibility','on'); - plot(lw_loop*1e-6,mean(squeeze(ber_lvlp(1:end,d,:)),1),"LineWidth",1,"Marker",m(1+3),"MarkerSize",5,'LineStyle','-.','Color',col(cnt,:),'DisplayName',['Lvl Lowpass ',comm_dn],'HandleVisibility','on'); - set(gca,'yscale','log'); - set(gca,'xscale','log'); - %xticklabels([10 100 1000 10000]); - grid minor - yline(3.8e-3,'LineWidth',2,'LineStyle','--','HandleVisibility','off'); - ylim([1e-3 4e-2]); - ylabel("BER"); - xlabel("Linewidth in MHz") - legend - cnt = cnt+1; -end - - -figure(1) -contour(lw_loop,sir_loop,thres_a0,16:0.4:21,'LineWidth',2,'FaceAlpha',0.3,'ShowText','on',"LabelFormat","%0.1f dB"); -clim([16 21]); -ylabel("Bandwidth in Multiples of Linewidth"); -xlabel("Linewidth in MHz"); -% set(gca,'yscale','log'); -grid minor -title("MPI removal - Optimization of Lowpass Filter Bandwidth") - -figure(2) -contour(lw_loop,sir_loop([1:10,12:end]),thres_a0([1:10,12:end],:),16:0.3:21,'LineWidth',2,'FaceAlpha',0.3,'ShowText','on',"LabelFormat","%0.1f dB"); -clim([16 21]); -ylabel("Window Length"); -xlabel("Linewidth in MHz"); -set(gca,'yscale','log'); -grid minor -title("MPI removal - Optimization of Averaging Window Length") - - - -%% Plot Winlen Contour -hdfec = 3.8e-3.*ones(size(sir_loop)); - -for dc = 1:size(data,2) - for lw = 1:size(data,3) - for s = 1:size(data,1) - ber_lvlp(s,dc,lw) = data{s,dc,lw}.ber_dcsm; - % ber_dcsm(wl,lw,s) = data{wl,s,lw}.ber_dcsm; - end - - a_bm = InterX([sir_loop;squeeze(ber_lvlp(:,dc,lw))'],[sir_loop;hdfec]); - % a_dcsm = InterX([sirloop;squeeze(ber_dcsm(wl,lw,:))'],[sirloop;hdfec]); - try - thres_a0(dc,lw) = -a_bm(1); - %thres_a1(wl,lw) = -a_dcsm(1); - catch - thres_a0(dc,lw) = NaN; - %thres_a1(wl,lw) = NaN; - end - - end -end - -figure(1) -contour(lw_loop,bw_loop,thres_a0,14:0.2:21,'LineWidth',1.5,'FaceAlpha',0.3,'ShowText','on',"LabelFormat","%0.1f dB"); -a = flip(cbrewer2('seq','Spectral',32)); -a = [a(1:12,:); a(22:end,:)]; -colormap(a); -clim([16 21]); -ylabel("Window Length"); -xlabel("Linewidth in MHz"); -yticks(bw_loop); -yticklabels(bw_loop); -set(gca,'yscale','log'); -set(gca,'xscale','log'); -grid minor - - - - -%% Plot DC Tap Contour -hdfec = 3.8e-3.*ones(size(bw_loop)); -thres_a0 = zeros(size(ber,1),size(ber,3)); -thres_a1 = zeros(size(ber,1),size(ber,3)); -for dc = 1:size(ber,1) - for lw = 1:size(ber,3) - a_lvsm = InterX([bw_loop;ber(dc,:,lw)],[bw_loop;hdfec]); - thres_a0(dc,lw) = -a_lvsm(1); - - a1 = InterX([bw_loop;ber_a1(dc,:,lw)],[bw_loop;hdfec]); - thres_a1(dc,lw) = -a1(1); - end -end - -figure(1) -subplot(2,1,1) -contour(lw_loop,sir_loop,thres_a0,16:0.5:21,'LineWidth',3,'FaceAlpha',0.3,'ShowText','on',"LabelFormat","%0.1f dB"); -clim([16 21]); -ylabel("DC Tap"); -xlabel("Linewidth in MHz"); -set(gca,'yscale','log'); -grid minor -subplot(2,1,2) -contour(lw_loop,sir_loop,thres_a1,16:0.5:21,'LineWidth',3,'FaceAlpha',0.3,'ShowText','on',"LabelFormat","%0.1f dB"); -clim([16 21]); -ylabel("DC Tap"); -xlabel("Linewidth in MHz"); -set(gca,'yscale','log'); -grid minor - -%% Plot Curves of required SIR to see the minimum a bit better -col = flip(cbrewer2('seq','Spectral',16)); -col = col([1:4, 10:end],:); -figure(2) -hold on -for lw = [size(ber,3):-2:2 2 1] - plot(sir_loop,thres_a0(:,lw),'DisplayName',[' Linewidth: ',num2str(lw_loop(lw)*1e-6), ' MHz'],'Color',col(lw,:),'Marker','o','MarkerFaceColor',col(lw,:),'LineWidth',2); - set(gca,'xscale','log'); -end -xlabel("DC Tap Value"); -ylabel("Required SIR to rech FEC in dB"); - -%% -col = linspecer(7); -figure(3) -cnt=1; -for lw = [1,2,11] - - for i = 1:length(sir_loop)-1 - - subplot(1,3,cnt) - hold on - plot(-1.*bw_loop,ber(i,:,lw),"LineWidth",2,"Marker","o","MarkerSize",5,'Color',col(i,:),'DisplayName',['DC tap ',num2str(sir_loop(i))]); - plot(-1.*bw_loop,ber_a1(i,:,lw),"LineWidth",2,"Marker","x","MarkerSize",5,'LineStyle','--','Color',col(i,:),'DisplayName',['A1. DC tap ',num2str(sir_loop(i))]); - yline(3.8e-3,'LineWidth',2,'LineStyle','--','HandleVisibility','off'); - set(gca,'yscale','log'); - grid minor - - xlim([15,30]); - ylim([1e-4,1e-2]); - xlabel("SIR in dB"); - ylabel("BER"); - title(['BER for different SIR;',' Linewidth: ',num2str(lw_loop(lw)*1e-6), ' MHz']) - text(25,4.2e-3,"FEC $3.8 e^{-3}$"); - end - cnt = cnt+1; -end -legend - + +clear + +sir_loop = [-25:2:-13]; +sir_loop = 0; +lw_loop = [1,2,3]; + +data = cell(length(sir_loop),length(lw_loop)); + +iterations=size(data); + +parfor ix = 1:numel(data) + + [u1,u2] = ind2sub(iterations,ix); + + output = imddmodel(sir_loop(u1),lw_loop(u2)); + + data{ix} = output; + +end + +for sir = 1:size(data,1) + for lw = 1:size(data,2) + ber_bm(sir,lw) = data{sir,lw}.ber_bm; + ber_dcsm(sir,lw) = data{sir,lw}.ber_dcsm; + ber_lvsm(sir,lw) = data{sir,lw}.ber_lvsm; + ber_lvlp(sir,lw) = data{sir,lw}.ber_lvlp; + end +end + +col = [ 0.6510 0.8078 0.8902 + 0.6980 0.8745 0.5412 + 0.9922 0.7490 0.4353]; +figure(22) +for l = 1:numel(lw_loop) + hold on + plot(sir_loop,ber_bm(:,l),'DisplayName',['DC mode:', num2str(lw_loop(l)),', Linewidth= 50 MHz'],'Marker','o','MarkerFaceColor',col(l,:),'Color',col(l,:),'LineWidth',2,'LineStyle','--'); +end +set(gca,'yscale','log'); +yline(3.8e-3,'LineWidth',2,'LineStyle','--','HandleVisibility','off'); +set(gca,'xdir','reverse'); + + + +%% Plot +col = linspecer(8); +figure(1) +hold on +m = ["x","o","pentagram","hexagram","*","+"]; +cnt = 1; +for d = 1%:size(data,2) + for dc = 1:size(data,1) + for lw = 1:size(data,3) + ber_bm(dc,d,lw) = (data{dc,d,lw}.ber_bm); + ber_dcsm(dc,d,lw) = data{dc,d,lw}.ber_dcsm; + ber_lvsm(dc,d,lw) = data{dc,d,lw}.ber_lvsm; + ber_lvlp(dc,d,lw) = (data{dc,d,lw}.ber_lvlp); + end + end + i = 1; + comm_dn = [];% ['EQ DC Tap: ',num2str(dc_tap_loop(d)),' m']; + title("Dependency on Laser Linewidth; B2B; Delay : 2*50m") + plot(lw_loop*1e-6,mean(squeeze(ber_bm(1:end,d,:)),1),"LineWidth",1.2,"Marker",m(1),"MarkerSize",5,'Color',col(cnt,:),'DisplayName',[' ',comm_dn]); + plot(lw_loop*1e-6,mean(squeeze(ber_dcsm(1:end,d,:)),1),"LineWidth",1,"Marker",m(1+4),"MarkerSize",5,'LineStyle','--','Color',col(cnt,:),'DisplayName',['DC smoothing ',comm_dn],'HandleVisibility','on'); + plot(lw_loop*1e-6,mean(squeeze(ber_lvsm(1:end,d,:)),1),"LineWidth",1.2,"Marker",m(1+1),"MarkerSize",5,'LineStyle',':','Color',col(cnt,:),'DisplayName',['Lvl Smoothing ',comm_dn],'HandleVisibility','on'); + plot(lw_loop*1e-6,mean(squeeze(ber_lvlp(1:end,d,:)),1),"LineWidth",1,"Marker",m(1+3),"MarkerSize",5,'LineStyle','-.','Color',col(cnt,:),'DisplayName',['Lvl Lowpass ',comm_dn],'HandleVisibility','on'); + set(gca,'yscale','log'); + set(gca,'xscale','log'); + %xticklabels([10 100 1000 10000]); + grid minor + yline(3.8e-3,'LineWidth',2,'LineStyle','--','HandleVisibility','off'); + ylim([1e-3 4e-2]); + ylabel("BER"); + xlabel("Linewidth in MHz") + legend + cnt = cnt+1; +end + + +figure(1) +contour(lw_loop,sir_loop,thres_a0,16:0.4:21,'LineWidth',2,'FaceAlpha',0.3,'ShowText','on',"LabelFormat","%0.1f dB"); +clim([16 21]); +ylabel("Bandwidth in Multiples of Linewidth"); +xlabel("Linewidth in MHz"); +% set(gca,'yscale','log'); +grid minor +title("MPI removal - Optimization of Lowpass Filter Bandwidth") + +figure(2) +contour(lw_loop,sir_loop([1:10,12:end]),thres_a0([1:10,12:end],:),16:0.3:21,'LineWidth',2,'FaceAlpha',0.3,'ShowText','on',"LabelFormat","%0.1f dB"); +clim([16 21]); +ylabel("Window Length"); +xlabel("Linewidth in MHz"); +set(gca,'yscale','log'); +grid minor +title("MPI removal - Optimization of Averaging Window Length") + + + +%% Plot Winlen Contour +hdfec = 3.8e-3.*ones(size(sir_loop)); + +for dc = 1:size(data,2) + for lw = 1:size(data,3) + for s = 1:size(data,1) + ber_lvlp(s,dc,lw) = data{s,dc,lw}.ber_dcsm; + % ber_dcsm(wl,lw,s) = data{wl,s,lw}.ber_dcsm; + end + + a_bm = InterX([sir_loop;squeeze(ber_lvlp(:,dc,lw))'],[sir_loop;hdfec]); + % a_dcsm = InterX([sirloop;squeeze(ber_dcsm(wl,lw,:))'],[sirloop;hdfec]); + try + thres_a0(dc,lw) = -a_bm(1); + %thres_a1(wl,lw) = -a_dcsm(1); + catch + thres_a0(dc,lw) = NaN; + %thres_a1(wl,lw) = NaN; + end + + end +end + +figure(1) +contour(lw_loop,bw_loop,thres_a0,14:0.2:21,'LineWidth',1.5,'FaceAlpha',0.3,'ShowText','on',"LabelFormat","%0.1f dB"); +a = flip(cbrewer2('seq','Spectral',32)); +a = [a(1:12,:); a(22:end,:)]; +colormap(a); +clim([16 21]); +ylabel("Window Length"); +xlabel("Linewidth in MHz"); +yticks(bw_loop); +yticklabels(bw_loop); +set(gca,'yscale','log'); +set(gca,'xscale','log'); +grid minor + + + + +%% Plot DC Tap Contour +hdfec = 3.8e-3.*ones(size(bw_loop)); +thres_a0 = zeros(size(ber,1),size(ber,3)); +thres_a1 = zeros(size(ber,1),size(ber,3)); +for dc = 1:size(ber,1) + for lw = 1:size(ber,3) + a_lvsm = InterX([bw_loop;ber(dc,:,lw)],[bw_loop;hdfec]); + thres_a0(dc,lw) = -a_lvsm(1); + + a1 = InterX([bw_loop;ber_a1(dc,:,lw)],[bw_loop;hdfec]); + thres_a1(dc,lw) = -a1(1); + end +end + +figure(1) +subplot(2,1,1) +contour(lw_loop,sir_loop,thres_a0,16:0.5:21,'LineWidth',3,'FaceAlpha',0.3,'ShowText','on',"LabelFormat","%0.1f dB"); +clim([16 21]); +ylabel("DC Tap"); +xlabel("Linewidth in MHz"); +set(gca,'yscale','log'); +grid minor +subplot(2,1,2) +contour(lw_loop,sir_loop,thres_a1,16:0.5:21,'LineWidth',3,'FaceAlpha',0.3,'ShowText','on',"LabelFormat","%0.1f dB"); +clim([16 21]); +ylabel("DC Tap"); +xlabel("Linewidth in MHz"); +set(gca,'yscale','log'); +grid minor + +%% Plot Curves of required SIR to see the minimum a bit better +col = flip(cbrewer2('seq','Spectral',16)); +col = col([1:4, 10:end],:); +figure(2) +hold on +for lw = [size(ber,3):-2:2 2 1] + plot(sir_loop,thres_a0(:,lw),'DisplayName',[' Linewidth: ',num2str(lw_loop(lw)*1e-6), ' MHz'],'Color',col(lw,:),'Marker','o','MarkerFaceColor',col(lw,:),'LineWidth',2); + set(gca,'xscale','log'); +end +xlabel("DC Tap Value"); +ylabel("Required SIR to rech FEC in dB"); + +%% +col = linspecer(7); +figure(3) +cnt=1; +for lw = [1,2,11] + + for i = 1:length(sir_loop)-1 + + subplot(1,3,cnt) + hold on + plot(-1.*bw_loop,ber(i,:,lw),"LineWidth",2,"Marker","o","MarkerSize",5,'Color',col(i,:),'DisplayName',['DC tap ',num2str(sir_loop(i))]); + plot(-1.*bw_loop,ber_a1(i,:,lw),"LineWidth",2,"Marker","x","MarkerSize",5,'LineStyle','--','Color',col(i,:),'DisplayName',['A1. DC tap ',num2str(sir_loop(i))]); + yline(3.8e-3,'LineWidth',2,'LineStyle','--','HandleVisibility','off'); + set(gca,'yscale','log'); + grid minor + + xlim([15,30]); + ylim([1e-4,1e-2]); + xlabel("SIR in dB"); + ylabel("BER"); + title(['BER for different SIR;',' Linewidth: ',num2str(lw_loop(lw)*1e-6), ' MHz']) + text(25,4.2e-3,"FEC $3.8 e^{-3}$"); + end + cnt = cnt+1; +end +legend + diff --git a/setup_simulation.m b/projects/IMDD_base_system/setup_simulation.m similarity index 97% rename from setup_simulation.m rename to projects/IMDD_base_system/setup_simulation.m index b9242c6..324a409 100644 --- a/setup_simulation.m +++ b/projects/IMDD_base_system/setup_simulation.m @@ -1,399 +1,399 @@ - - -lw = [0.1e6 1e6 10e6]; - -for lp = 1 - - rng(9); - - %% A) Set Simulation Variables - - sir = -18; - delay = 50; %mpi delay in meter - - fiblen = 0; %main link in km - - laser_linewidth =5e6; - - O = 19; %order of prbs - N = 2^(O-1); %length of prbs - [~,seed] = prbs(O,1); %initialize first seed of prbs - - % Modulation - M = 4; %PAM-M - bitpattern = zeros(N,log2(M)); - - % Symbol Rate - fsym = 112e9; - - % DAC Rate - fdac = 120e9; - - % Simulation oversampling rate "k"; - kover = 16; - - % ADC Rate - fadc = 256e9; - - % Simulation frequency in "analog domain" - fsimu = kover * fdac ; - - - - %% B) CONSTRUCT ALL CLASSES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - digimod = PAMmapper(M,0); - - pulseform = Pulseformer("pulseform","rrc","fdac",fdac,"fsym",fsym,"pulselength",32,"rrcalpha",0.027); - - awg = AWG('fdac',fdac,'kover',kover,'lpf_active',1,'f_cutoff',56e9,'lpf_type',filtertypes.gaussian,'bit_resolution',5.5); - - lp_laser = Filter('filtdegree',1,"f_cutoff",50e9,"fsamp",fdac*kover,"filterType",filtertypes.bessel_inp); - - u_pi = 3.5; - vbias = (0.5*u_pi)-u_pi; - extmodlaser = EML("mode",eml_mode.im_cosinus,"power",5,"fsimu",fsimu,"lambda",1310,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth); - - fib = Fiber("fsimu",fdac*kover,"fiber_length",fiblen,"alpha",0.2,"D",16,"lambda0",thz2nm(193.1),"gamma",0); - - reflectionpoint = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",sir); - reflectionprop = Fiber("fsimu",fdac*kover,"fiber_length",delay/1000,"alpha",0.2,"D",16,"lambda0",1550,"gamma",0); - - opticatten = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",-8); - - phdiode = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20); - - lp_diode = Filter('filtdegree',1,"f_cutoff",50e9,"fsamp",fdac*kover,"filterType",filtertypes.bessel_inp); - - scp = Scope("fsimu",fdac*kover,"fadc",fadc,... - "delay",0,"fixed_delay",0,"lpf_bw",113e9,"filtertype",filtertypes.butterworth,... - "samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,... - "adcresolution",6,"quantbuffer",0.1,'block_dc',1); - - eq = EQ("K",2,"plottrain",0,"plotfinal",0,... - "training_length",4096,"training_loops",2,... - "Ne",[50,5,3],"Nb",[3,0,0],... - "DCmu",0.05,"DDmu",[0.0004 0.0005 0.0006 0.0007 ],"DFEmu",0.005,"FFEmu",0.00,... - "dd_loops",2,"epsilon",[10 100 1000 ],"M",2,... - "thres",[0.005 0.004 0.0005 ],"l1act",0,"delay",0,"rho",0.0005,"ideal_dfe",0,"DB_aim",0); - - eq = EQ_silas("Ne",[50,5,3],"Nb",[2,0,0],"trainlength",4096,... - "sps",2,... - "mu_dc_dd",0.05,... - "mu_dc_train",0.05,... - "mu_ffe_train",0,... - "mu_dfe_train",0.005,... - "mu_ffe_dd",[0.0004 0.0004 0.0004],... - "mu_dfe_dd",0.005,... - "ddloops",3,... - "trainloops",4,... - "dcmode",1); - - - %% C) PROCESS TX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - % 1) PRBS Generation - for i = 1:log2(M) - [bitpattern(:,i),seed] = prbs(O,N,seed); - end - - % 2 ) Build Inf. signal class - bits = Informationsignal(bitpattern); - - % 3) Digi modulation -> PAM-M signal - digimod_out = digimod.map(bits); - - % 4) Pulse shaping -> racos - X = pulseform.process(digimod_out); - - % 5) AWG (lowpass, quantization, sample and hold) - X = awg.process(X); - - % 6) Lowpass behavior of laser and hf-cable? why twice? - X = lp_laser.process(X); - % X = lp_laser.process(X); - - % 7) Normalize signal - X = X.normalize("mode","oneone"); - X.signal = X.signal .* 1.3800; - - - %% D) PROCESS OPTICAL CHANNEL %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - % 1) Laser; Modulation -> OPTICAL DOMAIN - [X,extmodlaser] = extmodlaser.process(X); - - % 2) Main fiber propagation - X = fib.process(X); - - if delay ~= 0 - % 3) Reflection - % Reflection is just an attenuation - R = reflectionpoint.process(X); - - % Propagate back and forth (actual fiber propagation) - R = reflectionprop.process(R); - - % Delay the reflected signal - [R,n] = R.delay("delay_meter",delay); - - % Add together - col = cbrewer2('qual','Paired',8); - - xax = (1:X.length);% / fsimu * (physconst("LightSpeed")/1.4677); - xax_sec = (1:X.length) / fsimu .* 1e6;%* (physconst("LightSpeed")/1.4677); - xax_sec = xax_sec(n:end).'; - thresh = pi/4; - phaseX = phase(X.signal(n:end)); - phaseR = phase(R.signal(n:end)); - phasediff = wrapToPi(phaseX-phaseR); - [pos_high] = find(abs(phasediff)>thresh); - [pos_low] = find(abs(phasediff)<=thresh); - - - figure(6);hold on;histogram(phaseX-phaseR,500,'EdgeAlpha',0) - - figure() - subplot(3,1,1) - scatter(xax_sec.',wrapToPi(phaseX),4,'.','MarkerEdgeColor',col(2,:),'DisplayName','Desired Signal') - hold on - scatter(xax_sec,wrapToPi(phaseR),4,'.','MarkerEdgeColor',col(1,:),'DisplayName','Reflected Signal') - hold on - lg = legend; - lg.Location = "southwest"; - xlim([xax_sec(1) xax_sec(end)]); - xlabel('time in $\mu$s') - ylim([-pi pi]); - yticks([-pi 0 pi]) - yticklabels({'$-\pi$',0, '$\pi$'}) - ylabel('$\phi$') - - subplot(3,1,2) - scatter(xax_sec(pos_low),phasediff(pos_low),4,'.','MarkerEdgeColor',col(5,:),'DisplayName','Phase Difference') - hold on - scatter(xax_sec(pos_high),phasediff(pos_high),4,'.','MarkerEdgeColor',col(2,:),'HandleVisibility','off') -% yline(thresh,'LineWidth',2,'LineStyle','-','HandleVisibility','off') -% yline(-thresh,'LineWidth',2,'LineStyle','-','HandleVisibility','off') - lg = legend; - lg.Location = "southwest"; - xlabel('time in $\mu$s') - xlim([xax_sec(1) xax_sec(end)]); - ylim([-pi pi]); - yticks([-pi 0 pi]) - yticklabels({'$-\pi$',0, '$\pi$'}) - ylabel('$\Delta \phi$') - - - X.signal = X.signal(n:end); - R.signal = R.signal(n:end); - X1 = X; - X = X+R; - - subplot(3,1,3) - scatter(xax_sec(pos_low),abs(X.signal(pos_low).^2)*1000,4,'.','MarkerEdgeColor',col(5,:),'DisplayName','Constructive Interference'); - hold on - scatter(xax_sec(pos_high),abs(X.signal(pos_high).^2)*1000,4,'.','MarkerEdgeColor',col(2,:),'DisplayName','Destructive Interference'); - %scatter(xax_sec(pos_high),abs(X1.signal(pos_high).^2)*1000,4,'.','MarkerEdgeColor',col(1,:),'MarkerFaceAlpha',0.3,'DisplayName','Destructive Interference'); - xlim([xax_sec(1) xax_sec(end)]); - xlabel('time in $\mu$s') - lg = legend; - lg.Location = "southwest"; - ylabel('Optical power in mW'); - - - - - disp(['SIR ',num2str(10*log10(X.power/R.power))]); - - % cut reference signal to correct length (nessecary due to MPI delay) - digimod_out.signal = digimod_out.signal(round(n * fsym/fsimu) : end,:); - - bitpattern = bitpattern(round(n * fsym/fsimu):end,:); - - end - - -% plot(angle(R.signal)) - - % 4) Attenuation - X = opticatten.process(X); - - - % 5) Photo Diode -> ELECTRICAL DOMAIN - X = phdiode.process(X); - X = lp_diode.process(X); - - - %% E) PROCESS RX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - % 1) Oscilloscope (Sampling to f_adc; Quantization; Bandwidth Limitation) - X = scp.process(X); - - % 2) Resample to 2x symbol rate - X = X.resample("fs_out",2*fsym,"fs_in",fadc); - - % 3) Normalize - Eq_in = X.normalize("mode","rms"); - - % 4) Equalize - - % MPI reduction DC removal BEFORE EQ - wl = 3000; % symbols - Eq_in.signal = Eq_in.signal - 1/wl .* movsum( Eq_in.signal,[wl/2,wl/2]); - - % Equalize Signal - [Eq_out] = eq.process(Eq_in,digimod_out); - - - - %% A1: MPI reduction DC removal - wl = 1000; % symbols - yk_dcsm = Eq_out; - yk_dcsm.signal = Eq_out.signal - 1/wl .* movsum( Eq_out.signal,[wl/2,wl/2]); - - %% A2: MPI reduction Level wise error removal - yk_lvsm = Eq_out; - yk_lvlp = Eq_out; - pre_decision_level_uni = digimod.decide_pamlevel(Eq_out); - pre_decision_level_bi = ( pre_decision_level_uni*2-3 ) .* 1/sqrt(5); - - e = Eq_out.signal - pre_decision_level_bi; - - lp_mpi = Filter('filtdegree',1,"f_cutoff",2e6,"fsamp",fsym,"filterType",filtertypes.bessel_inp); - - filtered = lp_mpi.process(e); - - wl = 30; % symbols - smoothed = ( 1/wl .* movsum(e,[wl/2,wl/2]) ); - - % remove interference - for level = 0:3 - yk_lvsm.signal(pre_decision_level_uni==level) = yk_lvsm.signal(pre_decision_level_uni==level) - smoothed(pre_decision_level_uni==level); - yk_lvlp.signal(pre_decision_level_uni==level) = yk_lvlp.signal(pre_decision_level_uni==level) - filtered(pre_decision_level_uni==level); - end - - % Calc EVM - - evm_bm = calc_evm(Eq_out.signal, pre_decision_level_bi); - evm_dcsm = calc_evm(yk_dcsm.signal, pre_decision_level_bi); - evm_lsm = calc_evm(yk_lvsm.signal, pre_decision_level_bi); - evm_llp = calc_evm(yk_lvlp.signal, pre_decision_level_bi); - - % figure(1);bar([evm_bm' evm_dcsm' evm_llp' evm_lsm']);ylim([0.01 0.1]);set(gca,'yscale','log'); - - %% PROCESS RX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - % 1) Digi Demod - d_bm = digimod.demap(Eq_out); - d_dcsm = digimod.demap(yk_dcsm); - d_lvsm = digimod.demap(yk_lvsm); - d_lvlp = digimod.demap(yk_lvlp); - - % 2) BER - dbit = length(d_bm.signal)-length(bitpattern); - - [~,errors_bm,ber_bm,loc] = calc_ber(d_bm.signal(1:end-dbit,:) ,bitpattern(1:end,:),"skip",0,"returnErrorLocation",1); - [~,errors_dcsm,ber_dcsm] = calc_ber(d_dcsm.signal(1:end-dbit,:) ,bitpattern(1:end,:),"skip",0,"returnErrorLocation",0); - [~,errors_lvsm,ber_lvsm] = calc_ber(d_lvsm.signal(1:end-dbit,:) ,bitpattern(1:end,:),"skip",0,"returnErrorLocation",0); - [~,errors_lvlp,ber_lvlp] = calc_ber(d_lvlp.signal(1:end-dbit,:) ,bitpattern(1:end,:),"skip",0,"returnErrorLocation",0); - - % Display BER - disp(['BER benchmark: ', sprintf('%2E',ber_bm), ' ERRORS: ' ,num2str(sum(errors_bm))]); - disp(['BER dc smooth (A1): ', sprintf('%2E',ber_dcsm), ' ERRORS: ' ,num2str(sum(errors_dcsm))]); - disp(['BER lv smooth (A2): ', sprintf('%2E',ber_lvsm), ' ERRORS: ' ,num2str(sum(errors_lvsm))]); - disp(['BER lv lowpas: ', sprintf('%2E',ber_lvlp), ' ERRORS: ' ,num2str(sum(errors_lvlp))]); - - - %% Generate some Plots - - if 1 - - % SCATTER - col = cbrewer2('Paired',8); - - xax = 1:Eq_out.length; - - figure(3) - sgtitle('') - - subplot(1,3,1) - hold on - eq_decision = digimod.decide_pamlevel(Eq_out); - true_symbols = digimod.decide_pamlevel(digimod_out); - - xindices = 1:Eq_in.length; - errorpos = find(loc~=0)*2; - errorpos(errorpos>length(Eq_in.signal)) = length(Eq_in.signal); - correct = find(loc==0)*2; - correct(correct>length(Eq_in.signal)) = length(Eq_in.signal); - scatter(xindices(correct),Eq_in.signal(correct),4,'.','MarkerEdgeColor',col(4,:),'DisplayName','After EQ'); - hold on - scatter(xindices(errorpos),Eq_in.signal(errorpos),6,'x','MarkerEdgeColor',col(6,:),'DisplayName','Wrong Decision'); - hold off - - xlim([1, xindices(end)]); - ylim([-3 3]); - xlabel('Sampling Index') - ylabel('Amplitude') - a = legend; - a.Location = "best"; - - subplot(1,3,2) - xindices = 1:Eq_out.length; - xax_sec = (1:Eq_out.length) / fsym .* 1e6; - scatter(xax_sec(loc==0),Eq_out.signal(loc==0),4,'.','MarkerEdgeColor',col(4,:),'DisplayName','After EQ'); - hold on - scatter(xindices(loc~=0),Eq_out.signal(loc~=0),8,'x','MarkerEdgeColor',col(6,:),'DisplayName','Wrong Decision'); - hold off - xlim([1, xax(end)]); - ylim([-2 2]); - xlabel('Sampling Index') - ylabel('Amplitude') - legend - a = legend; - a.Location = "best"; - hold off - - subplot(1,3,3) - xindices = 1:yk_lvsm.length; - scatter(xindices(loc==0),yk_lvsm.signal(loc==0),4,'.','MarkerEdgeColor',col(4,:),'DisplayName','After A2'); - hold on - scatter(xindices(loc~=0),yk_lvsm.signal(loc~=0),8,'x','MarkerEdgeColor',col(6,:),'DisplayName','Wrong Decision'); - hold off - xlim([1, xax(end)]); - ylim([-2 2]); - xlabel('Sampling Index') - ylabel('Amplitude') - legend - a = legend; - a.Location = "best"; - hold off - - - end - - if 0 - col = cbrewer2('Paired',8); - - figure(11) - clf - - subplot(2,1,1) - hold on - plot(digimod_out.signal(4150:4175),'DisplayName','Tx','Color',col(1,:),'LineWidth',3); - plot(Eq_out.signal(4150:4175),'DisplayName','Rx after EQ','Color',col(6,:),'LineWidth',1); - title('Modulated Sequence Zoom'); - legend - hold off - subplot(2,1,2) - hold on - stem(d_bm.signal(4150:4175,1),'DisplayName','Tx','Color',col(1,:),'LineStyle','-','LineWidth',5) - stem(bitpattern(4150:4175,1)','DisplayName','Rx','Color',col(6,:),'LineStyle','--','LineWidth',2) - title('Bitpattern Tx - Rx'); - legend - hold off - - end - -end - + + +lw = [0.1e6 1e6 10e6]; + +for lp = 1 + + rng(9); + + %% A) Set Simulation Variables + + sir = -18; + delay = 50; %mpi delay in meter + + fiblen = 0; %main link in km + + laser_linewidth =5e6; + + O = 19; %order of prbs + N = 2^(O-1); %length of prbs + [~,seed] = prbs(O,1); %initialize first seed of prbs + + % Modulation + M = 4; %PAM-M + bitpattern = zeros(N,log2(M)); + + % Symbol Rate + fsym = 112e9; + + % DAC Rate + fdac = 120e9; + + % Simulation oversampling rate "k"; + kover = 16; + + % ADC Rate + fadc = 256e9; + + % Simulation frequency in "analog domain" + fsimu = kover * fdac ; + + + + %% B) CONSTRUCT ALL CLASSES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + digimod = PAMmapper(M,0); + + pulseform = Pulseformer("pulseform","rrc","fdac",fdac,"fsym",fsym,"pulselength",32,"rrcalpha",0.027); + + awg = AWG('fdac',fdac,'kover',kover,'lpf_active',1,'f_cutoff',56e9,'lpf_type',filtertypes.gaussian,'bit_resolution',5.5); + + lp_laser = Filter('filtdegree',1,"f_cutoff",50e9,"fsamp",fdac*kover,"filterType",filtertypes.bessel_inp); + + u_pi = 3.5; + vbias = (0.5*u_pi)-u_pi; + extmodlaser = EML("mode",eml_mode.im_cosinus,"power",5,"fsimu",fsimu,"lambda",1310,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth); + + fib = Fiber("fsimu",fdac*kover,"fiber_length",fiblen,"alpha",0.2,"D",16,"lambda0",thz2nm(193.1),"gamma",0); + + reflectionpoint = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",sir); + reflectionprop = Fiber("fsimu",fdac*kover,"fiber_length",delay/1000,"alpha",0.2,"D",16,"lambda0",1550,"gamma",0); + + opticatten = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",-8); + + phdiode = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20); + + lp_diode = Filter('filtdegree',1,"f_cutoff",50e9,"fsamp",fdac*kover,"filterType",filtertypes.bessel_inp); + + scp = Scope("fsimu",fdac*kover,"fadc",fadc,... + "delay",0,"fixed_delay",0,"lpf_bw",113e9,"filtertype",filtertypes.butterworth,... + "samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,... + "adcresolution",6,"quantbuffer",0.1,'block_dc',1); + + eq = EQ("K",2,"plottrain",0,"plotfinal",0,... + "training_length",4096,"training_loops",2,... + "Ne",[50,5,3],"Nb",[3,0,0],... + "DCmu",0.05,"DDmu",[0.0004 0.0005 0.0006 0.0007 ],"DFEmu",0.005,"FFEmu",0.00,... + "dd_loops",2,"epsilon",[10 100 1000 ],"M",2,... + "thres",[0.005 0.004 0.0005 ],"l1act",0,"delay",0,"rho",0.0005,"ideal_dfe",0,"DB_aim",0); + + eq = EQ_silas("Ne",[50,5,3],"Nb",[2,0,0],"trainlength",4096,... + "sps",2,... + "mu_dc_dd",0.05,... + "mu_dc_train",0.05,... + "mu_ffe_train",0,... + "mu_dfe_train",0.005,... + "mu_ffe_dd",[0.0004 0.0004 0.0004],... + "mu_dfe_dd",0.005,... + "ddloops",3,... + "trainloops",4,... + "dcmode",1); + + + %% C) PROCESS TX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % 1) PRBS Generation + for i = 1:log2(M) + [bitpattern(:,i),seed] = prbs(O,N,seed); + end + + % 2 ) Build Inf. signal class + bits = Informationsignal(bitpattern); + + % 3) Digi modulation -> PAM-M signal + digimod_out = digimod.map(bits); + + % 4) Pulse shaping -> racos + X = pulseform.process(digimod_out); + + % 5) AWG (lowpass, quantization, sample and hold) + X = awg.process(X); + + % 6) Lowpass behavior of laser and hf-cable? why twice? + X = lp_laser.process(X); + % X = lp_laser.process(X); + + % 7) Normalize signal + X = X.normalize("mode","oneone"); + X.signal = X.signal .* 1.3800; + + + %% D) PROCESS OPTICAL CHANNEL %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % 1) Laser; Modulation -> OPTICAL DOMAIN + [X,extmodlaser] = extmodlaser.process(X); + + % 2) Main fiber propagation + X = fib.process(X); + + if delay ~= 0 + % 3) Reflection + % Reflection is just an attenuation + R = reflectionpoint.process(X); + + % Propagate back and forth (actual fiber propagation) + R = reflectionprop.process(R); + + % Delay the reflected signal + [R,n] = R.delay("delay_meter",delay); + + % Add together + col = cbrewer2('qual','Paired',8); + + xax = (1:X.length);% / fsimu * (physconst("LightSpeed")/1.4677); + xax_sec = (1:X.length) / fsimu .* 1e6;%* (physconst("LightSpeed")/1.4677); + xax_sec = xax_sec(n:end).'; + thresh = pi/4; + phaseX = phase(X.signal(n:end)); + phaseR = phase(R.signal(n:end)); + phasediff = wrapToPi(phaseX-phaseR); + [pos_high] = find(abs(phasediff)>thresh); + [pos_low] = find(abs(phasediff)<=thresh); + + + figure(6);hold on;histogram(phaseX-phaseR,500,'EdgeAlpha',0) + + figure() + subplot(3,1,1) + scatter(xax_sec.',wrapToPi(phaseX),4,'.','MarkerEdgeColor',col(2,:),'DisplayName','Desired Signal') + hold on + scatter(xax_sec,wrapToPi(phaseR),4,'.','MarkerEdgeColor',col(1,:),'DisplayName','Reflected Signal') + hold on + lg = legend; + lg.Location = "southwest"; + xlim([xax_sec(1) xax_sec(end)]); + xlabel('time in $\mu$s') + ylim([-pi pi]); + yticks([-pi 0 pi]) + yticklabels({'$-\pi$',0, '$\pi$'}) + ylabel('$\phi$') + + subplot(3,1,2) + scatter(xax_sec(pos_low),phasediff(pos_low),4,'.','MarkerEdgeColor',col(5,:),'DisplayName','Phase Difference') + hold on + scatter(xax_sec(pos_high),phasediff(pos_high),4,'.','MarkerEdgeColor',col(2,:),'HandleVisibility','off') +% yline(thresh,'LineWidth',2,'LineStyle','-','HandleVisibility','off') +% yline(-thresh,'LineWidth',2,'LineStyle','-','HandleVisibility','off') + lg = legend; + lg.Location = "southwest"; + xlabel('time in $\mu$s') + xlim([xax_sec(1) xax_sec(end)]); + ylim([-pi pi]); + yticks([-pi 0 pi]) + yticklabels({'$-\pi$',0, '$\pi$'}) + ylabel('$\Delta \phi$') + + + X.signal = X.signal(n:end); + R.signal = R.signal(n:end); + X1 = X; + X = X+R; + + subplot(3,1,3) + scatter(xax_sec(pos_low),abs(X.signal(pos_low).^2)*1000,4,'.','MarkerEdgeColor',col(5,:),'DisplayName','Constructive Interference'); + hold on + scatter(xax_sec(pos_high),abs(X.signal(pos_high).^2)*1000,4,'.','MarkerEdgeColor',col(2,:),'DisplayName','Destructive Interference'); + %scatter(xax_sec(pos_high),abs(X1.signal(pos_high).^2)*1000,4,'.','MarkerEdgeColor',col(1,:),'MarkerFaceAlpha',0.3,'DisplayName','Destructive Interference'); + xlim([xax_sec(1) xax_sec(end)]); + xlabel('time in $\mu$s') + lg = legend; + lg.Location = "southwest"; + ylabel('Optical power in mW'); + + + + + disp(['SIR ',num2str(10*log10(X.power/R.power))]); + + % cut reference signal to correct length (nessecary due to MPI delay) + digimod_out.signal = digimod_out.signal(round(n * fsym/fsimu) : end,:); + + bitpattern = bitpattern(round(n * fsym/fsimu):end,:); + + end + + +% plot(angle(R.signal)) + + % 4) Attenuation + X = opticatten.process(X); + + + % 5) Photo Diode -> ELECTRICAL DOMAIN + X = phdiode.process(X); + X = lp_diode.process(X); + + + %% E) PROCESS RX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % 1) Oscilloscope (Sampling to f_adc; Quantization; Bandwidth Limitation) + X = scp.process(X); + + % 2) Resample to 2x symbol rate + X = X.resample("fs_out",2*fsym,"fs_in",fadc); + + % 3) Normalize + Eq_in = X.normalize("mode","rms"); + + % 4) Equalize + + % MPI reduction DC removal BEFORE EQ + wl = 3000; % symbols + Eq_in.signal = Eq_in.signal - 1/wl .* movsum( Eq_in.signal,[wl/2,wl/2]); + + % Equalize Signal + [Eq_out] = eq.process(Eq_in,digimod_out); + + + + %% A1: MPI reduction DC removal + wl = 1000; % symbols + yk_dcsm = Eq_out; + yk_dcsm.signal = Eq_out.signal - 1/wl .* movsum( Eq_out.signal,[wl/2,wl/2]); + + %% A2: MPI reduction Level wise error removal + yk_lvsm = Eq_out; + yk_lvlp = Eq_out; + pre_decision_level_uni = digimod.decide_pamlevel(Eq_out); + pre_decision_level_bi = ( pre_decision_level_uni*2-3 ) .* 1/sqrt(5); + + e = Eq_out.signal - pre_decision_level_bi; + + lp_mpi = Filter('filtdegree',1,"f_cutoff",2e6,"fsamp",fsym,"filterType",filtertypes.bessel_inp); + + filtered = lp_mpi.process(e); + + wl = 30; % symbols + smoothed = ( 1/wl .* movsum(e,[wl/2,wl/2]) ); + + % remove interference + for level = 0:3 + yk_lvsm.signal(pre_decision_level_uni==level) = yk_lvsm.signal(pre_decision_level_uni==level) - smoothed(pre_decision_level_uni==level); + yk_lvlp.signal(pre_decision_level_uni==level) = yk_lvlp.signal(pre_decision_level_uni==level) - filtered(pre_decision_level_uni==level); + end + + % Calc EVM + + evm_bm = calc_evm(Eq_out.signal, pre_decision_level_bi); + evm_dcsm = calc_evm(yk_dcsm.signal, pre_decision_level_bi); + evm_lsm = calc_evm(yk_lvsm.signal, pre_decision_level_bi); + evm_llp = calc_evm(yk_lvlp.signal, pre_decision_level_bi); + + % figure(1);bar([evm_bm' evm_dcsm' evm_llp' evm_lsm']);ylim([0.01 0.1]);set(gca,'yscale','log'); + + %% PROCESS RX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % 1) Digi Demod + d_bm = digimod.demap(Eq_out); + d_dcsm = digimod.demap(yk_dcsm); + d_lvsm = digimod.demap(yk_lvsm); + d_lvlp = digimod.demap(yk_lvlp); + + % 2) BER + dbit = length(d_bm.signal)-length(bitpattern); + + [~,errors_bm,ber_bm,loc] = calc_ber(d_bm.signal(1:end-dbit,:) ,bitpattern(1:end,:),"skip",0,"returnErrorLocation",1); + [~,errors_dcsm,ber_dcsm] = calc_ber(d_dcsm.signal(1:end-dbit,:) ,bitpattern(1:end,:),"skip",0,"returnErrorLocation",0); + [~,errors_lvsm,ber_lvsm] = calc_ber(d_lvsm.signal(1:end-dbit,:) ,bitpattern(1:end,:),"skip",0,"returnErrorLocation",0); + [~,errors_lvlp,ber_lvlp] = calc_ber(d_lvlp.signal(1:end-dbit,:) ,bitpattern(1:end,:),"skip",0,"returnErrorLocation",0); + + % Display BER + disp(['BER benchmark: ', sprintf('%2E',ber_bm), ' ERRORS: ' ,num2str(sum(errors_bm))]); + disp(['BER dc smooth (A1): ', sprintf('%2E',ber_dcsm), ' ERRORS: ' ,num2str(sum(errors_dcsm))]); + disp(['BER lv smooth (A2): ', sprintf('%2E',ber_lvsm), ' ERRORS: ' ,num2str(sum(errors_lvsm))]); + disp(['BER lv lowpas: ', sprintf('%2E',ber_lvlp), ' ERRORS: ' ,num2str(sum(errors_lvlp))]); + + + %% Generate some Plots + + if 1 + + % SCATTER + col = cbrewer2('Paired',8); + + xax = 1:Eq_out.length; + + figure(3) + sgtitle('') + + subplot(1,3,1) + hold on + eq_decision = digimod.decide_pamlevel(Eq_out); + true_symbols = digimod.decide_pamlevel(digimod_out); + + xindices = 1:Eq_in.length; + errorpos = find(loc~=0)*2; + errorpos(errorpos>length(Eq_in.signal)) = length(Eq_in.signal); + correct = find(loc==0)*2; + correct(correct>length(Eq_in.signal)) = length(Eq_in.signal); + scatter(xindices(correct),Eq_in.signal(correct),4,'.','MarkerEdgeColor',col(4,:),'DisplayName','After EQ'); + hold on + scatter(xindices(errorpos),Eq_in.signal(errorpos),6,'x','MarkerEdgeColor',col(6,:),'DisplayName','Wrong Decision'); + hold off + + xlim([1, xindices(end)]); + ylim([-3 3]); + xlabel('Sampling Index') + ylabel('Amplitude') + a = legend; + a.Location = "best"; + + subplot(1,3,2) + xindices = 1:Eq_out.length; + xax_sec = (1:Eq_out.length) / fsym .* 1e6; + scatter(xax_sec(loc==0),Eq_out.signal(loc==0),4,'.','MarkerEdgeColor',col(4,:),'DisplayName','After EQ'); + hold on + scatter(xindices(loc~=0),Eq_out.signal(loc~=0),8,'x','MarkerEdgeColor',col(6,:),'DisplayName','Wrong Decision'); + hold off + xlim([1, xax(end)]); + ylim([-2 2]); + xlabel('Sampling Index') + ylabel('Amplitude') + legend + a = legend; + a.Location = "best"; + hold off + + subplot(1,3,3) + xindices = 1:yk_lvsm.length; + scatter(xindices(loc==0),yk_lvsm.signal(loc==0),4,'.','MarkerEdgeColor',col(4,:),'DisplayName','After A2'); + hold on + scatter(xindices(loc~=0),yk_lvsm.signal(loc~=0),8,'x','MarkerEdgeColor',col(6,:),'DisplayName','Wrong Decision'); + hold off + xlim([1, xax(end)]); + ylim([-2 2]); + xlabel('Sampling Index') + ylabel('Amplitude') + legend + a = legend; + a.Location = "best"; + hold off + + + end + + if 0 + col = cbrewer2('Paired',8); + + figure(11) + clf + + subplot(2,1,1) + hold on + plot(digimod_out.signal(4150:4175),'DisplayName','Tx','Color',col(1,:),'LineWidth',3); + plot(Eq_out.signal(4150:4175),'DisplayName','Rx after EQ','Color',col(6,:),'LineWidth',1); + title('Modulated Sequence Zoom'); + legend + hold off + subplot(2,1,2) + hold on + stem(d_bm.signal(4150:4175,1),'DisplayName','Tx','Color',col(1,:),'LineStyle','-','LineWidth',5) + stem(bitpattern(4150:4175,1)','DisplayName','Rx','Color',col(6,:),'LineStyle','--','LineWidth',2) + title('Bitpattern Tx - Rx'); + legend + hold off + + end + +end + diff --git a/projects/MPI_April/mpi_simulation_cspr.m b/projects/MPI_April/mpi_simulation_cspr.m new file mode 100644 index 0000000..a9962e1 --- /dev/null +++ b/projects/MPI_April/mpi_simulation_cspr.m @@ -0,0 +1,198 @@ +%% Settings +clear + +filename = '112G_2'; +load_sequence = 0; + +M = 4; +datarate = 112e9; + +kover = 8; +fsym = round(datarate*1e-9 / log2(M))*1e9; +%fsym = 50e9; +fdac = 256e9; +fadc = 256e9; + +laser_linewidth = 1e6; +mpi_ = 50; %meter +sir = 20; %decibel = attenuation of interference path + +lowpass_cutoff = fsym/2 * 1; +awg_bw = lowpass_cutoff; +mod_bw = lowpass_cutoff; +phd_bw = lowpass_cutoff; +scp_bw = lowpass_cutoff; + +mpi_path=50; + +vp = [0.25]; +vb = [0.5:0.05:0.8]; +vb = 0.7; +pn_key = [27]; +rop = -5; + +LP_awg = Filter('filtdegree',2,"f_cutoff",awg_bw,"fs",fdac,"filterType","butterworth"); +LP_laser = Filter('filtdegree',4,"f_cutoff",mod_bw,"fs",fdac*kover,"filterType",filtertypes.gaussian); +LP_opt = Filter('filtdegree',4,"f_cutoff",fsym/log2(M).*1.5,"fs",fdac*kover,"filterType",filtertypes.gaussian); +LP_phd = Filter('filtdegree',2,"f_cutoff",phd_bw,"fs",fdac*kover,"filterType",filtertypes.butterworth); +LP_scpe = Filter('filtdegree',2,"f_cutoff",scp_bw,"fs",fadc,"filterType","butterworth"); + +for pnk = 1:length(pn_key) + for m = 1:length(vb) + if load_sequence + + load(['C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\MPI_April\',filename,'.mat'],'X'); + + spectrum_plot(X.signal',X.fs,'spectrum'); + else + + % 1) PRBS Generation + O = 17; %order of prbs + N = 2^(O-1); %length of prbs + [~,seed] = prbs(O,1); %initialize first seed of prbs + bitpattern=[]; + for i = 1:log2(M) + [bitpattern(:,i),seed] = prbs(O,N,seed); + end + + % 2 ) Build Inf. signal class + bits = Informationsignal(bitpattern); + + % 3) Digi modulation -> PAM-M signal + digimod_out = PAMmapper(M,0).map(bits); + digimod_out.fs = fsym; + + X = Pulseformer("fsym",fsym,"fdac",fdac,"pulse","rrc","pulselength",16,"rrcalpha",0.01).process(digimod_out); + + % 5) AWG (lowpass, quantization, sample and hold) + %X = M8196A().process(X); + kover = 16; + X = AWG("fdac",fdac,"dac_min",-1,"dac_max",1,"H_lpf",LP_awg,"kover",kover).process(X); + + % 6) Lowpass behavior before laser + X = LP_laser.process(X); + + % 7) Normalize signal + X = X.normalize("mode","oneone"); + + %save(['C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\MPI_April\',char(filename)],'X'); + end + + % 1) Laser; Modulation -> OPTICAL DOMAIN + u_pi = 2; + vbias = -vb(m)*u_pi; + extmodlaser = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",X.fs,"lambda",1290,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",pn_key(pnk)); + E = X.*vp(1); + [Opt,extmodlaser] = extmodlaser.process(E); + +% figure(m) +% hold on +% scatter(E.signal(1:100000),(abs(Opt.signal(1:100000)).^2)*1e3,0.1,'.','DisplayName','Modulator TF') +% xlabel('Input in V') +% ylabel('abs(Output) in mW') + + Opt = LP_opt.process(Opt); + + cspr(m) = Opt.cspr; + + % 2) ping pong fiber propagation + Interference_sig = Fiber("fsimu",Opt.fs,"fiber_length",mpi_path*2/1000,"alpha",0,"D",0,"lambda0",1310,"gamma",0).process(Opt); + + Interference_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",-sir).process(Interference_sig); + + % In the meantime: delay the main signal + [Main_sig,n] = Opt.delay("delay_meter",mpi_path*2); + + % Add + Combined_sig = Main_sig + Interference_sig; + + % Cut (due to the delays there is a jump in the signals) + if n == 0;n = 1;end + Combined_sig.signal = Combined_sig.signal(ceil(n):end); + + % Fiber + Combined_sig = Fiber("fsimu",Combined_sig.fs,"fiber_length",0,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.08).process(Combined_sig); + + for i = 1:length(rop) + + % Set ROP + Rx_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",rop(i)).process(Combined_sig); + rop_save(m,i) = Rx_sig.power; + + % Square Law + Rx_sig = Photodiode("fsimu",Rx_sig.fs,"dark_current",2e-08,"responsivity",1,"temperature",20).process(Rx_sig); + + %Lowpass PhDiode + Rx_sig = LP_phd.process(Rx_sig); + + % Scope + Rx_sig = Scope("fsimu",Rx_sig.fs,"fadc",fadc,... + "delay",0,"fixed_delay",0,"lpf_bw",scp_bw,"filtertype",filtertypes.butterworth,... + "samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,... + "adcresolution",16,"quantbuffer",0.1,'block_dc',1,'lpf_active',1,'H_lpf',LP_scpe).process(Rx_sig); + + % Sample to 2x fsym + Rx_sig = Rx_sig.resample("fs_in",fadc,"fs_out",2*fsym); + + % Sync Rx signal with reference + [Rx_sig,D,cuts] = Rx_sig.tsynch("reference",digimod_out,"fs_ref",fsym); + + % % simple EQ (optimum mudc: 0.05 -> 0.005) + EQ_sig = EQ_silas_plain("Ne",[20,8,8],"Nb",[2,0,0],"trainlength",4096,"mu_dc_dd",0.005,"mu_dc_train",0.05,... + "mu_ffe_train",0.005,"mu_combined_dd",[0.0004 0.0006 0.0003 0.005],"ddloops",3,'trainloops',3,'sps',2).process(Rx_sig,digimod_out); + + % Demap + Rx_Bits = PAMmapper(M,0).demap(EQ_sig); + + Rx_symboldecision = PAMmapper(M,0).decide_pamlevel(EQ_sig,"symbol_levels",unique(digimod_out.signal)); + + levels = PAMmapper(M,0).separate_pamlevels(EQ_sig); + level_avg(:,m,i) = mean(levels,'omitnan'); + level_std(:,m,i) = std(levels,'omitnan'); + + % BER + [~,errors_bm,BER(m,i),errors] = calc_ber(Rx_Bits.signal,bitpattern,"skip_front",0,"skip_end",0,"returnErrorLocation",1); + + formatted_ber = sprintf('%.1e', BER(m,i)); + disp(formatted_ber); + + end + end +end +plot_analysis_window; +save(['C:\Users\Silas\Nextcloud4\Dokumente\02_Ablage_Office\MPI\Investigation_April_2024\','PAM_',num2str(M),'_mpi_',num2str(mpi_path),'_lw_',num2str(laser_linewidth)],"BER"); + +% BER plot +figure(340) +cols = linspecer(7); +for m = 1:size(BER,1) + hold on + plot(rop,BER(m,:),'DisplayName',['Bias: ',num2str(vb(m)), ' V'],'LineStyle','-','Color',cols(m,:),'LineWidth',1,'Marker','o','MarkerEdgeColor',[1 1 1],'MarkerFaceColor',cols(m,:)); + +end +set(gca,'YScale','log'); +legend +xlabel("ROP in dBm") +yline(3.8e-3,'DisplayName','FEC'); + + +figure(21) +hold on +plot(rop,mean(BER),'DisplayName',['Modulation: ',num2str(2*vp/extmodlaser.u_pi*100), ' $\%$'],'LineStyle','-','Color',cols(2,:),'LineWidth',1); +set(gca,'YScale','log'); +legend +ylabel("ROP in dBm") +yline(3.8e-3,'DisplayName','FEC'); + + + + + + + + +%check Rx and TX symbols +figure(101) +scatter(1:100,Rx_symboldecision.signal(1:100),10,'o'); +hold on +scatter(1:100,digimod_out.signal(1:100),5,'x'); diff --git a/projects/MPI_April/plot_analysis_window.m b/projects/MPI_April/plot_analysis_window.m new file mode 100644 index 0000000..6d3f443 --- /dev/null +++ b/projects/MPI_April/plot_analysis_window.m @@ -0,0 +1,70 @@ + +%straight outta workspace, show everything I want + +cols = linspecer(8); +colpairs = cbrewer2('paired',8); + +f=figure(20); +clf(f) +f.Name="Optspec"; +%Spectrum of Optical Signal after superposition +subplot(3,2,1:2) +spectrum_plot(Combined_sig.normalize("mode","rms").signal',Combined_sig.fs,'Optspec',['Papr: ',num2str(Combined_sig.papr), '; CSPR: ',num2str(Combined_sig.cspr)]); +hold on +ylim([-100,0]); + +%Phase Investigation +subplot(3,2,3) +hold on +phase_int = extmodlaser.phase(n:end); +phase_main = delayseq(extmodlaser.phase,n); +phase_main = phase_main(n:end); +phase_diff = phase_int - phase_main; +t = (1:length(phase_int))' ./ extmodlaser.fsimu ; +plot(t*1e6,phase_int,'Color',colpairs(1,:),'DisplayName','Interferer Phase'); +plot(t*1e6,phase_main,'Color',colpairs(2,:),'DisplayName','Signal Phase (delayed)'); +plot(t*1e6,phase_diff,'Color',colpairs(4,:),'DisplayName','Delta Phase'); +legend + +subplot(3,2,4) +%Received Signal after Phdiode +yyaxis left +t = (1:Rx_sig.length)' ./ Rx_sig.fs; +scatter(t*1e6,Rx_sig.normalize("mode","oneone").signal.*max(unique(digimod_out.signal)),1,'.','MarkerEdgeColor',cols(5,:)); +hold on +errors_t = errors./Rx_Bits.fs; +errors_sym = Rx_symboldecision.signal(errors); +scatter(errors_t*1e6,errors_sym,2,'x','LineWidth',1); +yline(PAMmapper(M,0).thresholds); + +%Again: Phase Diff +yyaxis right +t = (1:length(phase_int))' ./ extmodlaser.fsimu ; +plot(t*1e6,phase_diff,'Color',colpairs(4,:),'DisplayName','Delta Phase'); + +title(['Lwidth:',num2str(laser_linewidth) ,'; SIR: ', num2str(sir) ,'dB; ROP: ',num2str(rop(i)),' dBm']) +xlabel('t in $\mu$s') +legend + +%Equalized Signal +subplot(3,2,5) +t = (1:EQ_sig.length)' ./ EQ_sig.fs; +for h = 1:size(levels,2) + scatter(t*1e6,levels(:,h),1,'.'); + hold on +end +yline(PAMmapper(M,0).thresholds); +title(['BER: ',sprintf('%.1e', BER(m,i))]); + +%Histogram of EQzed Signal +subplot(3,2,6) +for h = 1:size(levels,2) + std_dev = std(levels(:,h),'omitnan'); + mean_val = mean(levels(:,h),'omitnan'); + hold on + histogram(levels(:,h),1000,"EdgeColor","none","Normalization","pdf"); + text(mean_val, 0.1, ['$\sigma^2:$ ',sprintf('%.2f', std_dev)], 'HorizontalAlignment', 'center'); +end +hold on +xline(PAMmapper(M,0).thresholds); +xticks(sort([unique(digimod_out.signal);PAMmapper(M,0).thresholds'])); diff --git a/projects/MPI_April/simulate_tx_sig.m b/projects/MPI_April/simulate_tx_sig.m new file mode 100644 index 0000000..0752c99 --- /dev/null +++ b/projects/MPI_April/simulate_tx_sig.m @@ -0,0 +1,55 @@ +function simulate_tx_sig(filename) + + %% Params + M=4; + digimod = PAMmapper(M,0); + + fdac = 120e9; + fsym = 112e9; + pulseform = Pulseformer("pulseform","rrc","fdac",fdac,"fsym",fsym,"pulselength",32,"rrcalpha",0.027); + + kover = 16; + awg = AWG('fdac',fdac,'kover',kover,'lpf_active',1,'f_cutoff',56e9,'lpf_type',filtertypes.gaussian,'bit_resolution',5.5); + %awg = M8199B(); + + lp_laser = Filter('filtdegree',2,"f_cutoff",50e9,"fsamp",fdac*kover,"filterType",filtertypes.butterworth); + + + + %% PROCESS TX + + % 1) PRBS Generation + O = 20; %order of prbs + N = 2^(O-1); %length of prbs + [~,seed] = prbs(O,1); %initialize first seed of prbs + for i = 1:log2(M) + [bitpattern(:,i),seed] = prbs(O,N,seed); + end + + % 2 ) Build Inf. signal class + bits = Informationsignal(bitpattern); + + % 3) Digi modulation -> PAM-M signal + digimod_out = digimod.map(bits); + + % 4) Pulse shaping -> racos + X = pulseform.process(digimod_out); + + % 5) AWG (lowpass, quantization, sample and hold) + X = awg.process(X); + + % 6) Lowpass behavior before laser + X = lp_laser.process(X); + + % 7) Normalize signal + X = X.normalize("mode","oneone"); + X.signal = X.signal; + + spectrum_plot(X.signal',X.fs,'spectrum'); + + if nargin == 1 + save(['C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\',char(filename)],'X'); + end + +end + diff --git a/projects/MPI_offline_analysis/mpi_analysis_pam4.m b/projects/MPI_offline_analysis/mpi_analysis_pam4.m new file mode 100644 index 0000000..a75b9f6 --- /dev/null +++ b/projects/MPI_offline_analysis/mpi_analysis_pam4.m @@ -0,0 +1,444 @@ + +clear; +col = linspecer(6); + +% GENERATE SIR CURVE +if ismac + foldername = '/Users/silasoettinghaus/Documents/MATLAB/Labor_Datensatz_PAM4_MPI/pam4_10km_1km'; +else + foldername = 'C:\Users\Silas\Documents\MATLAB\Datensätze\Labor_Datensatz_PAM4_MPI_OFC2023\pam4_10km_1km'; +end + + +% SETTINGS + +optimize_mudc = 0; +run_sir_sweep = 0; +run_feed_forward = 1; +run_baseline = 0; +run_ideal_dc_tap = 0; +plot_timesignal = 1; + +block_loop = [1]; + +for bl = 1:numel(block_loop) + + eq_parallelization_blocklength = block_loop(bl); + eq_updatelatency = 3; + eq_avg_blocklength =0; + + % FIND BEST MU DC + + if optimize_mudc + mudc_loop = [0, 0.0001,0.0005, 0.001,0.005, 0.01:0.01:0.1, 0.2:0.1:1]; + + current_filename = 'pam4__loop_14_92Gbd_19092023_1513.mat'; + current_filename = 'pam4__loop_29_56Gbd_19092023_1546.mat'; + recorded_data = load([foldername,filesep, current_filename]); + + [best_mudc,ber]= optimizeMuDc(recorded_data,eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength,mudc_loop); + +% figure(16) +% hold on +% plot(mudc_loop,ber); +% set(gca, 'YScale', 'log'); +% set(gca, 'XScale', 'log'); +% scatter(best_mudc,ber(mudc_loop==best_mudc),100,'Marker','x','LineWidth',2) +% yline(ber(1)); +% xlim([0,1]) + end + + if run_sir_sweep + + mudc = best_mudc; + [ber,sir,fsym] = runSIRsweep(eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength,mudc); + + plotSirSweep(ber,sir,fsym,eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength,mudc) + + crossing = calculateCrossing([92e9, 56e9],ber,sir,fsym); + + req_sir_56(bl) = crossing(2,2) ; + req_sir_92(bl) = crossing(2,1) ; + + end + + +end + +if run_feed_forward + + current_filename = 'pam4__loop_14_92Gbd_v19092023_1513.mat'; + recorded_data = load([foldername,filesep, current_filename]); + eq_avg_blocklength = [50,100,1000,3000]; + [best_block,ber]= optimizeAvgBlocklength(recorded_data,eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength,0); + + + best_block = 400; + + [ber_base,sir_base,fsym_base] = runSIRsweep(eq_parallelization_blocklength,eq_updatelatency,best_block,0); + + plotSirSweep(ber_base,sir_base,fsym_base,eq_parallelization_blocklength,eq_updatelatency,best_block,0) + + crossing = calculateCrossing([92e9, 56e9],ber_base,sir_base,fsym_base); + + [ber_base,sir_base,fsym_base] = runSIRsweep(eq_parallelization_blocklength,eq_updatelatency,0,0); + + plotSirSweep(ber_base,sir_base,fsym_base,eq_parallelization_blocklength,eq_updatelatency,0,0) + + crossing = calculateCrossing([92e9, 56e9],ber_base,sir_base,fsym_base); + +end + + +if run_baseline + + mudc = 0; + [ber_base,sir_base,fsym_base] = runSIRsweep(eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength,mudc); + + plotSirSweep(ber_base,sir_base,fsym_base,eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength,mudc) + + crossing = calculateCrossing([92e9, 56e9],ber_base,sir_base,fsym_base); + + req_sir_56_baseline = crossing(2,2) ; + req_sir_92_baseline = crossing(2,1) ; + +end + +if run_ideal_dc_tap + + current_filename = 'pam4__loop_29_56Gbd_19092023_1546.mat'; + recorded_data = load([foldername,filesep, current_filename]); + [mudc,~]= optimizeMuDc(recorded_data,eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength,mudc_loop); + + eq_parallelization_blocklength = 1; + eq_updatelatency = 1; + [ber_ideal,sir_ideal,fsym_ideal] = runSIRsweep(eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength,mudc); + + plotSirSweep(ber_ideal,sir_ideal,fsym_ideal,eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength,mudc) + + crossing = calculateCrossing([92e9, 56e9],ber_ideal,sir_ideal,fsym_ideal); + + req_sir_56_ideal = crossing(2,2) ; + req_sir_92_ideal = crossing(2,1) ; + +end + + + +if plot_timesignal + + current_filename = 'pam4__loop_26_92Gbd_19092023_1542.mat'; + recorded_data = load([foldername,filesep, current_filename]); + + sir = recorded_data.saveStructTemp.awg2scope_keysight_state.eigenlight_mpi - 3.3 + 7.2; + + fsym = 92e9; + eq_parallelization_blocklength = 1; + eq_updatelatency = 1; + eq_avg_blocklength = 0; + + % TRUE SYMBOLS + correct_symbols = recorded_data.saveStructTemp.digi_mod_out'; + + %3) PLAIN RECEIVED SIGNAL + disp('RX Signal') + y_rx = Electricalsignal(recorded_data.saveStructTemp.dp_tsynch_out'); + y_rx.fs = 2.*fsym; + rx_symbols = y_rx.resample("fs_in",2.*fsym,"fs_out",fsym); + rx_symbols = rx_symbols.normalize("mode","rms").signal; + disp(std(rx_symbols)); + scatterleveldependent(rx_symbols,correct_symbols,fsym); + + %1) PLOT WITH WITH ALGORITHM + if 0 + mudc_loop = [0, 0.0001,0.0005, 0.001,0.005, 0.01:0.01:0.1, 0.2:0.1:1]; + [mudc,~]= optimizeMuDc(recorded_data,eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength,mudc_loop); + end + + disp('ALGORITHM ') + mudc = 0.05; + results_alg = runPostProcessing(recorded_data,mudc,eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength); + disp(results_alg.ber); + rx_symbols = results_alg.EQ_out.signal; + disp(std(rx_symbols)); + scatterleveldependent(rx_symbols,correct_symbols,fsym); + + %2) JUST EQ + disp('Just EQ') + mudc = 0; + results = runPostProcessing(recorded_data,mudc,eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength); + disp(results.ber) + rx_symbols = results.EQ_out.signal; + disp(std(rx_symbols)); + scatterleveldependent(rx_symbols,correct_symbols,fsym); + + + + +end + + +function plotSirSweep(ber,sir,fsym,eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength,mudc) + + + rate = [92e9, 56e9]; + figure(95) + + for r = 1:length(rate) + + sorted = sortrows([sir(fsym == rate(r)); ber(fsym == rate(r))]',1)'; + hold on + plot(abs(sorted(1,:)),sorted(2,:),'DisplayName',[num2str(rate(r).*1e-9),' GBd; PAM4; mudc: ',num2str(mudc)],'LineWidth',1,'Marker','o','MarkerSize',5,'LineStyle','-','HandleVisibility','on'); + + end + + set(gca, 'YScale', 'log'); + yline(3.8e-3,'LineWidth',1, 'LineStyle','--','HandleVisibility','off'); + xlim([15,35]); + xlabel('SIR in dB'); + ylabel('BER'); + +end + +function [crossing] = calculateCrossing(rate,ber,sir,fsym) + + for r = 1:length(rate) + sorted = sortrows([sir(fsym == rate(r)); ber(fsym == rate(r))]',1)'; + berVals = sorted(2,:); + xAxis = sorted(1,:); + + hdfec = 3.8e-3 .* ones(1,length(sorted)); + crossing(:,r) = InterX([hdfec(:)';xAxis],[berVals;xAxis]) ; + end + +end + +function [ber,sir,fsym] = runSIRsweep(eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength,mudc) + + % GENERATE SIR CURVE + if ismac + foldername = '/Users/silasoettinghaus/Documents/MATLAB/Labor_Datensatz_PAM4_MPI/pam4_10km_1km'; + else + foldername = 'C:\Users\Silas\Documents\MATLAB\Labor_Datensatz_PAM4_MPI\pam4_10km_1km'; + end + + allfiles = dir(foldername); + + parfor i = 1:length(allfiles) + + if allfiles(i).bytes ~= 0 + current_filename = allfiles(i).name; + recorded_data = load([foldername,filesep, current_filename]); + else + continue + end + + results = runPostProcessing(recorded_data,mudc,eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength); + ber(i) = results.ber; + fsym(i) = results.fsym; + sir(i) = results.sir + + disp(['fsym: ',num2str(results.fsym*1e-9),'GBd, SIR: ',num2str(results.sir),' ->> BER: ',sprintf('%2E',results.ber)]); + + end +end + +function [best_blocklength,ber] = optimizeAvgBlocklength(recorded_data,eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength_loop,mudc) + + for j = 1:length(eq_avg_blocklength_loop) + eq_avg_blocklength = eq_avg_blocklength_loop(j); + + + results = runPostProcessing(recorded_data,mudc,eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength); + ber(j) = results.ber; + + disp(['fsym: ',num2str(results.fsym*1e-9),'GBd, SIR: ',num2str(results.sir),' ->> BER: ',sprintf('%2E',results.ber)]); + + end + + [~,pos]=min(ber); + best_blocklength = eq_avg_blocklength_loop(pos); + +end + +function [best_mudc,ber] = optimizeMuDc(recorded_data,eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength,mudc_loop) + + parfor j = 1:length(mudc_loop) + mudc = mudc_loop(j); + + + results = runPostProcessing(recorded_data,mudc,eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength); + ber(j) = results.ber; + + disp(['fsym: ',num2str(results.fsym*1e-9),'GBd, SIR: ',num2str(results.sir),' ->> BER: ',sprintf('%2E',results.ber)]); + + end + + [~,pos]=min(ber); + best_mudc = mudc_loop(pos); + +end + +function results = runPostProcessing(recorded_data,mudc,eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength) + + results.sir = recorded_data.saveStructTemp.awg2scope_keysight_state.eigenlight_mpi - 3.3 + 7.2; + results.fsym = recorded_data.saveStructTemp.common.f_sym; + + results.fdac = recorded_data.saveStructTemp.common.f_DAC; + results.fadc = recorded_data.saveStructTemp.common.f_ADC; + + + % 0) build RX Signal + y_rx = Electricalsignal(recorded_data.saveStructTemp.dp_tsynch_out'); + y_rx.fs = 2.*results.fsym; + + % 0) build tx reference Signal for eq training + y_digimod = Electricalsignal(recorded_data.saveStructTemp.digi_mod_out'); + y_digimod.fs = results.fsym; + + + % 1) normlaize + y_rx = y_rx.normalize("mode","rms"); + + eq = EQ_silas("Ne",[25,0,0],"Nb",[2,0,0],"trainlength",4096,... + "sps",2,... + "mu_dc_dd",mudc,... + "mu_dc_train",mudc,... + "mu_ffe_train",0,... + "mu_dfe_train",0.005,... + "mu_ffe_dd",[0.0004 0.0004 0.0004],... + "mu_dfe_dd",0.005,... + "ddloops",3,... + "trainloops",4,... + "eq_parallelization_blocklength",eq_parallelization_blocklength, ... + "eq_updatelatency",eq_updatelatency,... + "eq_avg_blocklength",eq_avg_blocklength); + + [Eq_out] = eq.process(y_rx,y_digimod); + + results.EQ_out = Eq_out; + + % 3) digital demodulation object + digimod = PAMmapper(2^recorded_data.saveStructTemp.common.M,0); + + %estimated/ equalized symbol sequence + d_estimated = digimod.demap(Eq_out); + + %correct data symbols + d_correct = recorded_data.saveStructTemp.prms_out; + + % 4) BER calculation + [totalbits,errors,results.ber,loc] = calc_ber(d_estimated.signal(1:end,:) ,d_correct(1:end,:)',"skip",0,"returnErrorLocation",1); +end + +function std_per_lvl = calcleveldependentstd(rx_symbols,correct_symbols) + error_of_rx_signal = rx_symbols - correct_symbols; + levels = unique(correct_symbols); + for l = 1:4 + level_amplitude = levels(l); + std_per_lvl(l) = var(( 1/64 .* movsum(error_of_rx_signal(correct_symbols==level_amplitude),[64/2,64/2]) )); + %std_per_lvl(l) = std(error_of_rx_signal(correct_symbols==level_amplitude)); + end +end + +function scatterleveldependent(rx_symbols,correct_symbols,f_sym) + col = cbrewer2('Paired',8); + ccnt = -1; + figure1 = figure(); + levels = unique(correct_symbols); + + start = 1; + ende = length(correct_symbols); + + start = 30000; + ende = 40000; + + for l = 1:4 + ccnt = ccnt+2; + + level_amplitude = levels(l); + + symbols_for_lvl = NaN(1,length(correct_symbols)); + symbols_for_lvl(correct_symbols==level_amplitude) = rx_symbols(correct_symbols==level_amplitude); + std_lvl(l) = std(symbols_for_lvl,'omitnan'); + xax_in_sec = ((1:length(correct_symbols)) / f_sym) * 1e6; + xax_in_sec = 1:length(correct_symbols); + + scatter(xax_in_sec(start:ende),symbols_for_lvl(start:ende),10,'.','MarkerFaceAlpha',0.5,'MarkerEdgeAlpha',0.5,'MarkerEdgeColor',col(ccnt,:)); + hold on; + end + + std_lvl = round(std_lvl,2); + disp(std_lvl); + + ccnt = 0; + + % Add the windowed/ smoothed curves + for l = 1:4 + ccnt = ccnt+2; + level_amplitude = levels(l); + + symbols_for_lvl = NaN(1,length(correct_symbols)); + + movmean = 1/250 .* movsum(rx_symbols(correct_symbols==level_amplitude),[250/2,250/2]); + + symbols_for_lvl(correct_symbols==level_amplitude) = movmean; + + nanx = isnan(symbols_for_lvl); + t = 1:numel(symbols_for_lvl); + symbols_for_lvl(nanx) = interp1(t(~nanx), symbols_for_lvl(~nanx), t(nanx)); + + xax_in_sec = ((1:length(correct_symbols)) / f_sym) * 1e6; + xax_in_sec = 1:length(correct_symbols); + + plot(xax_in_sec(start:ende),symbols_for_lvl(start:ende),'Color',col(ccnt,:)); + hold on + end + + %yline(max(rx_symbols(correct_symbols==levels(2)))) + + if 0 + annotation(figure1,'textbox',... + [0.660523809523809 0.844444444444448 0.133523809523809 0.0603174603174607],... + 'String',['\sigma = ',num2str(std_lvl(4))],... + 'LineWidth',1.8,... + 'LineStyle','none',... + 'FontSize',12,... + 'FitBoxToText','off'); + + % Create textbox + annotation(figure1,'textbox',... + [0.667666666666665 0.642857142857147 0.133523809523809 0.0603174603174607],... + 'String',['\sigma = ',num2str(std_lvl(3))],... + 'LineWidth',1.8,... + 'LineStyle','none',... + 'FontSize',12,... + 'FitBoxToText','off'); + + % Create textbox + annotation(figure1,'textbox',... + [0.671238095238093 0.442857142857148 0.133523809523809 0.0603174603174608],... + 'String',['\sigma = ',num2str(std_lvl(2))],... + 'LineWidth',1.8,... + 'LineStyle','none',... + 'FontSize',12,... + 'FitBoxToText','off'); + + % Create textbox + annotation(figure1,'textbox',... + [0.670047619047616 0.265079365079371 0.133523809523809 0.0603174603174608],... + 'String',['\sigma = ',num2str(std_lvl(1))],... + 'LineWidth',1.8,... + 'LineStyle','none',... + 'FontSize',12,... + 'FitBoxToText','off'); + end + + xlim([0, 2.6]) + ylim([-2 2]) + xlabel('Time in $\mu$s'); + ylabel('Normalized Amplitude'); + +end + + diff --git a/projects/MPI_offline_analysis/mpi_analysis_pam4_old_notworking.m b/projects/MPI_offline_analysis/mpi_analysis_pam4_old_notworking.m new file mode 100644 index 0000000..df67869 --- /dev/null +++ b/projects/MPI_offline_analysis/mpi_analysis_pam4_old_notworking.m @@ -0,0 +1,188 @@ + + + +% parfsym = zeros(length(mudc_loop),length(allfiles)); +% parsir = zeros(length(mudc_loop),length(allfiles)); +% ber_wo_reduction = zeros(length(mudc_loop),length(allfiles)); +% ber_w_dcremoval = zeros(length(mudc_loop),length(allfiles)); +% mudc = zeros(length(mudc_loop),length(allfiles)); +% blocklen = zeros(length(mudc_loop),length(allfiles)); + +% parameters.mudc = [0, 0.0001,0.0005, 0.001,0.005, 0.01:0.01:0.1, 0.2:0.1:1]; +% parameters.eq_parallelization_blocklength = [1]; +% parameters.eq_updatelatency = [1]; +% parameters.eq_avg_blocklength =[0]; +block_loop = [1]; + +for j = 1:length(mudc_loop) + + mudc_current = mudc_loop(j); + block_current = block_loop(1); + + parfor i = 1:length(allfiles) + + if 1 + + if allfiles(i).bytes ~= 0 + current_filename = allfiles(i).name; + else + continue + end + + else +% +% current_filename = 'pam4__loop_29_56Gbd_19092023_1546.mat'; +% current_filename = 'pam4__loop_30_92Gbd_19092023_1547.mat'; +% + current_filename = 'pam4__loop_13_56Gbd_19092023_1512.mat'; + current_filename = 'pam4__loop_14_92Gbd_19092023_1513.mat'; + + end + + % load dataset + recorded_data = load([foldername,filesep, current_filename]); + % get common variables of simulation loop + + %even if redundant, save stuff to results struct + parsir(j,i) = recorded_data.saveStructTemp.awg2scope_keysight_state.eigenlight_mpi - 3.3 + 7.2; + parfsym(j,i) = recorded_data.saveStructTemp.common.f_sym; + + fdac = recorded_data.saveStructTemp.common.f_DAC; + fadc = recorded_data.saveStructTemp.common.f_ADC; + fsym = recorded_data.saveStructTemp.common.f_sym; + + % 0) build RX Signal + + y_rx = Electricalsignal(recorded_data.saveStructTemp.dp_tsynch_out'); + + y_rx.fs = 2.*fsym; + + % 0) build tx reference Signal for eq training + y_digimod = Electricalsignal(recorded_data.saveStructTemp.digi_mod_out'); + y_digimod.fs = fsym; + + rx_symbols = y_rx.resample("fs_in",2.*fsym,"fs_out",fsym); + rx_symbols = rx_symbols.normalize("mode","rms").signal; + correct_symbols = recorded_data.saveStructTemp.digi_mod_out'; + +% scatterleveldependent(rx_symbols,correct_symbols,parfsym(j,i)); +% title([num2str(fsym.*1e-9),' GBd; SIR: ',num2str(parsir(j,i)),' dB']); + +% if parfsym(j,i) == 56e9 +% deviation56(i) = var(( 1/500 .* movsum(y_rx.signal,[500/2,500/2]) )); +% sir56(i) = simdata.saveStructTemp.awg2scope_keysight_state.eigenlight_mpi - 3.3 + 7.2; +% std56(i,:) = calcleveldependentstd(rx_symbols,correct_symbols); +% else +% deviation92(i) = var(( 1/500 .* movsum(y_rx.signal,[500/2,500/2]) )); +% sir92(i) = simdata.saveStructTemp.awg2scope_keysight_state.eigenlight_mpi - 3.3 + 7.2; +% std92(i,:) = calcleveldependentstd(rx_symbols,correct_symbols); +% end + + % 1) normlaize + y_rx = y_rx.normalize("mode","rms"); + + eq = EQ_silas("Ne",[25,0,0],"Nb",[2,0,0],"trainlength",4096,... + "sps",2,... + "mu_dc_dd",mudc_current,... + "mu_dc_train",mudc_current,... + "mu_ffe_train",0,... + "mu_dfe_train",0.005,... + "mu_ffe_dd",[0.0004 0.0004 0.0004],... + "mu_dfe_dd",0.005,... + "ddloops",3,... + "trainloops",4,... + "eq_parallelization_blocklength",block_current, ... + "eq_updatelatency",3,... + "eq_avg_blocklength",10000); + + [Eq_out] = eq.process(y_rx,y_digimod); + + % 2.2) MPI reduction blocks + + % 3) digital demodulation + digimod = PAMmapper(2^recorded_data.saveStructTemp.common.M,0); + + d_hat = digimod.demap(Eq_out); + + d_ = digimod.demap(Electricalsignal(recorded_data.saveStructTemp.digi_mod_out')); + d = recorded_data.saveStructTemp.prms_out; + +% scatterleveldependent(Eq_out.signal,correct_symbols,parfsym(j,i)); +% title([num2str(fsym.*1e-9),' GBd; SIR: ',num2str(parsir(j,i)),' dB']); + + % 4) BER calculation + [totalbits,errors_bm,ber_dhat,loc] = calc_ber(d_hat.signal(1:end,:) ,d(1:end,:)',"skip",0,"returnErrorLocation",1); + + ber_nodctap = mean(recorded_data.saveStructTemp.prms_compare_state.BER ); + + disp(['fsym: ',num2str(parfsym(j,i).*1e-9),'GBd, SIR: ',num2str(parsir(j,i)),' ->> BER_new: ',sprintf('%2E',ber_dhat),' BER_old: ',sprintf('%2E',ber_nodctap)]); + + % 5) save to struct + + ber_wo_reduction(j,i) = ber_nodctap; + ber_w_dcremoval(j,i) = ber_dhat; + mudc(j,i) = eq.mu_dc_dd; + blocklen(j,i) = eq.eq_parallelization_blocklength; + + end +end + + + +rate = [92e9, 56e9]; +figure(95) +col = cbrewer2('Paired',8); +for r = 1:length(rate) + for j = 1:length(mudc_loop) + + + if mudc_loop(j) == 0 + lnsty = '--'; + + mark = 'square'; + + else + lnsty = "-"; + + mark = 'v'; + end + + if r == 1 + handlvis = 'on'; + colr = col(6,:); + else + handlvis = 'off'; + colr = col(8,:); + end + + %plot curve of standard FFE+DFE + if j == 1 + load(['C:\Users\Silas\Nextcloud4\Dokumente\02_Ablage_Office\OFC2024\',num2str(rate(r).*1e-9),'GBd_no_dcs.mat'],'b'); + plot(abs(b(1,:)),b(2,:),'DisplayName',[num2str(rate(r).*1e-9),' GBd; PAM4; mudc: ',0],'LineWidth',1,'Marker','o','MarkerSize',5,'LineStyle',lnsty,'HandleVisibility',handlvis,'Color',colr); + end + + hold on + b = [parsir(j,parfsym(j,:)==rate(r));ber_w_dcremoval(j,parfsym(j,:)==rate(r))]; + b=sort(b,2); + plot(abs(b(1,:)),b(2,:),'DisplayName',[num2str(rate(r).*1e-9),' GBd; PAM4; mudc: ',num2str(mudc_loop(j))],'LineWidth',1,'Marker',mark,'MarkerSize',5,'LineStyle',lnsty,'HandleVisibility',handlvis,'Color',colr); + + end +end +set(gca, 'YScale', 'log'); +yline(3.8e-3,'LineWidth',1, 'LineStyle','--','HandleVisibility','off'); +xlim([15,35]); +xlabel('SIR in dB'); +ylabel('BER'); +annotation('textbox',... + [0.328380952380952 0.477777777777781 0.197809523809524 0.106349206349207],... + 'String',{'3.8 $\cdot 10^{-3}$'},... + 'LineStyle','none',... + 'Interpreter','latex',... + 'FitBoxToText','off'); + +% legend('Standard EQ','Adaptive DC Subtraction') + +figure(16) +hold on +plot(mudc_loop,ber_w_dcremoval(:,1));set(gca, 'YScale', 'log');set(gca, 'XScale', 'log'); +xlim([0,1])