From bb228ae2bdf201924fe17b7ab0f6ec66fefd3dc1 Mon Sep 17 00:00:00 2001 From: Silas Oettinghaus Date: Mon, 2 Sep 2024 09:00:41 +0200 Subject: [PATCH] Code for Deliverable 05 --- Classes/00_signals/Opticalsignal.m | 11 +- Classes/00_signals/Signal.m | 15 +- Classes/01_transmit/AWG.m | 1 + Classes/01_transmit/PAMsource.m | 16 +- Classes/04_DSP/FFE_DCremoval.m | 2 +- Classes/04_DSP/FFE_FFDCAVG.m | 2 +- Classes/04_DSP/MRDS_coding.m | 269 ++++++++++++++++++ Classes/04_DSP/VNLE.m | 3 + Datatypes/coding_mode.m | 8 + .../MPI_August/auswertung/clipfactor_curve.m | 65 +++++ .../auswertung/dc_removal_sir_curve.m | 59 ++++ .../auswertung/dc_removal_window_length.m | 71 +++++ projects/MPI_August/auswertung/hpf_curve.m | 76 +++++ .../hpf_curve_different_filter_bw.m | 59 ++++ .../MPI_August/auswertung/mrds_sir_curve.m | 59 ++++ projects/MPI_August/auswertung/rop_curve.m | 35 ++- .../MPI_August/auswertung/sensitivity_curve.m | 81 ++++++ projects/MPI_August/auswertung/sir_curve.m | 2 +- projects/MPI_August/mpi_uncorrelated.m | 216 ++++++++++---- projects/MPI_August/simualte_MPI_data.m | 117 ++++++++ 20 files changed, 1094 insertions(+), 73 deletions(-) create mode 100644 Classes/04_DSP/MRDS_coding.m create mode 100644 Datatypes/coding_mode.m create mode 100644 projects/MPI_August/auswertung/clipfactor_curve.m create mode 100644 projects/MPI_August/auswertung/dc_removal_sir_curve.m create mode 100644 projects/MPI_August/auswertung/dc_removal_window_length.m create mode 100644 projects/MPI_August/auswertung/hpf_curve.m create mode 100644 projects/MPI_August/auswertung/hpf_curve_different_filter_bw.m create mode 100644 projects/MPI_August/auswertung/mrds_sir_curve.m create mode 100644 projects/MPI_August/auswertung/sensitivity_curve.m create mode 100644 projects/MPI_August/simualte_MPI_data.m diff --git a/Classes/00_signals/Opticalsignal.m b/Classes/00_signals/Opticalsignal.m index ee6137f..25900a4 100644 --- a/Classes/00_signals/Opticalsignal.m +++ b/Classes/00_signals/Opticalsignal.m @@ -53,19 +53,22 @@ classdef Opticalsignal < Signal function cspr = cspr(obj) - carrier_power_dbm = pow2db( abs(mean(obj.signal)).^2 )+30; % dB -> +30 -> dBm + carrier_power_dbm = pow2db( abs(mean(obj.signal)).^2 ) +30; % dB -> +30 -> dBm signal_power_dbm = obj.power; + cspr = carrier_power_dbm-signal_power_dbm; carr_lin = abs(mean(obj.signal)).^2; sign_lin = mean(abs(obj.signal).^2); cspr_lin = pow2db(carr_lin/sign_lin); + %carrier power is the mean value of the overall signal -> dc part c = abs(mean(obj.signal)).^2; - s = mean(abs(obj.signal).^2); - cspr_ = 10*log10(c / s); + %signal power is now only the "fluctuation"/ i.e. the ac part + s = mean( abs(obj.signal-mean(obj.signal)).^2 ); + cspr = 10*log10(c / s); - cspr =carrier_power_dbm-signal_power_dbm; + end diff --git a/Classes/00_signals/Signal.m b/Classes/00_signals/Signal.m index bafbef7..48ddbf0 100644 --- a/Classes/00_signals/Signal.m +++ b/Classes/00_signals/Signal.m @@ -264,8 +264,7 @@ classdef Signal % spectrum_plot(obj.signal,options.fsamp,options.figurename,options.displayname); - N = 2^(nextpow2(length(obj.signal))-6); - + N = 2^(nextpow2(length(obj.signal))-8); [p_lin,w] = pwelch(obj.signal,hanning(N),N/2,N,obj.fs,"centered","power","mean"); p_dbm = 10*log10(p_lin)+30; %dB to dBm in case of "power" @@ -278,7 +277,7 @@ classdef Signal ylabel("Power (dBm)"); xlim([-obj.fs/2 obj.fs/2].*1e-9) edgetick = 2^(nextpow2(obj.fs*1e-9)); - xticks([-edgetick:16:edgetick]); + % xticks([-edgetick:16:edgetick]); xlim([-244, 244]) ylim([-120,-0]); yticks([-200:10:10]); @@ -329,7 +328,7 @@ classdef Signal case power_notation.mW pow_pk = pow_pk .* 1e3; %mW case power_notation.W - %pow = pow % Watt + pow_pk = pow_pk; % Watt end end @@ -562,6 +561,9 @@ classdef Signal maxA = max(sig(100:end-100)); minA = min(sig(100:end-100)); + % maxA = 0.0015; + % minA = 0; + difference= maxA-minA; @@ -603,6 +605,7 @@ classdef Signal min_ = min(obj.signal(100:end-100)); max_ = abs(max(obj.signal(100:end-100))); end + xlabel('Time in ps') % add information @@ -704,8 +707,12 @@ classdef Signal yticks(linspace(0,histpoints,16)); + y_tickstring = sprintfc('%.2f', y_tickstring); yticklabels(y_tickstring); + xticks(linspace(0,histpoints_horizontal,8)) + x_tickstring = sprintfc('%.2f', linspace(0, 2/fsym, 8) .* 1e12); + xticklabels(x_tickstring); end diff --git a/Classes/01_transmit/AWG.m b/Classes/01_transmit/AWG.m index 79c5469..4a60a1d 100644 --- a/Classes/01_transmit/AWG.m +++ b/Classes/01_transmit/AWG.m @@ -135,6 +135,7 @@ classdef AWG < handle %%%%%%%%% PRECOMP SINC ROLLOFF %%%%%%%%% if 1 % X: design FIR filter for sinc precomp + % https://www.dsprelated.com/showarticle/1191.php ntaps = 13; npts = 32; % least-squares FIR design diff --git a/Classes/01_transmit/PAMsource.m b/Classes/01_transmit/PAMsource.m index 279a913..8df05af 100644 --- a/Classes/01_transmit/PAMsource.m +++ b/Classes/01_transmit/PAMsource.m @@ -9,6 +9,9 @@ classdef PAMsource fsym randkey + mrds_code + mrds_blocklength + applypulseform pulseformer @@ -30,6 +33,9 @@ classdef PAMsource options.fsym = 112e9; options.randkey = 0; + options.mrds_code = 0; + options.mrds_blocklength = 512; + options.applypulseform = 1; options.pulseformer Pulseformer; @@ -84,22 +90,28 @@ classdef PAMsource symbols = PAMmapper(obj.M,0).map(bits); symbols.fs = obj.fsym; + + if obj.mrds_code + symbols = MRDS_coding("blocklength",obj.mrds_blocklength).encode(symbols); + end if obj.applyclipping sym_min = min(symbols.signal); sym_max = max(symbols.signal); end - %%%%% Pulseforming %%%%%% + %%%%% Pulse-forming %%%%%% if obj.applypulseform digi_sig = obj.pulseformer.process(symbols); else digi_sig = symbols; end - %%%%% Resample to f DAC %%%%%% + %%%%% Re-sample to f DAC %%%%%% digi_sig = digi_sig.resample("fs_in",digi_sig.fs,"fs_out",obj.fs_out,"n",10,"beta",5); + % digi_sig.spectrum("fignum",111,"displayname","after pulseforming"); + %%%%% Hard clip digital signal to PAM range before DAC %%%%%% if obj.applyclipping try diff --git a/Classes/04_DSP/FFE_DCremoval.m b/Classes/04_DSP/FFE_DCremoval.m index 8ea8f6a..30f689f 100644 --- a/Classes/04_DSP/FFE_DCremoval.m +++ b/Classes/04_DSP/FFE_DCremoval.m @@ -136,7 +136,7 @@ classdef FFE_DCremoval < handle %Update the dc estimation every n-th symbol. This is a %trivial implementation of parallel EQ´s where the %errors are not apparent in every step. See Silas OFC - %2023 "MPI mitigation adaptive DC removal" + %2023 "MPI mitigation adaptive DC removal if mod(symbol,length(e_dc_buffer)) == 0 e_dc_buffer(1) = e_dc_est - obj.mu_dc * err(symbol); e_dc_buffer = circshift(e_dc_buffer,1); diff --git a/Classes/04_DSP/FFE_FFDCAVG.m b/Classes/04_DSP/FFE_FFDCAVG.m index f72e027..c0da9a6 100644 --- a/Classes/04_DSP/FFE_FFDCAVG.m +++ b/Classes/04_DSP/FFE_FFDCAVG.m @@ -129,7 +129,7 @@ classdef FFE_FFDCAVG < handle for epoch = 1 : epochs symbol = 0; - err_buffer = zeros(numel(obj.constellation),90); + err_buffer = zeros(numel(obj.constellation),112); dc_err = zeros(numel(obj.constellation),1); dc_sto = NaN(numel(obj.constellation),N); diff --git a/Classes/04_DSP/MRDS_coding.m b/Classes/04_DSP/MRDS_coding.m new file mode 100644 index 0000000..b7be55b --- /dev/null +++ b/Classes/04_DSP/MRDS_coding.m @@ -0,0 +1,269 @@ +classdef MRDS_coding + %MRDS implementation according to: + + % Optical Multi-Path Interference Mitigation for PAM4-IMDD Systems Using Balanced Coding + % Journal of Lightwave Technology; 2024 + + properties(Access=public) + + blocklength + delta + end + + methods (Access=public) + function obj = MRDS_coding(options) + %NAME Construct an instance of this class + % Detailed explanation goes here + + arguments + options.blocklength = 8; + end + + % + fn = fieldnames(options); + for n = 1:numel(fn) + try + obj.(fn{n}) = options.(fn{n}); + end + end + + + + end + + function process(~) + error("MRDS_coding has no process function. Use .encode(signal) and .dc_remove(signal) and .decode(signal)"); + end + + function signalclass_out = encode(obj,signalclass_in) + + data_in = signalclass_in.signal'; + + if mean(unique(data_in)) < 0.01 % --> check for bipolar + data_in = int32(data_in.*sqrt(5)); + data_in = double(data_in); + else % unipolar + data_in = int32(data_in.*sqrt(5)*2-3); % make bipolar [-3, -1, 1, 3] + data_in = double(data_in); + end + + data_out = obj.mrds_encoding(data_in, obj.blocklength); + data_out = data_out./sqrt(5); % normalized to Power=1 + + signalclass_in.signal = data_out'; + + % append to logbook + lbdesc = ['MRDS Coded']; + signalclass_in = signalclass_in.logbookentry(lbdesc); + + % write to output + signalclass_out = signalclass_in; + + end + + function signalclass_out = decode(obj,signalclass_in) + + data_in = signalclass_in.signal'; + + if mean(unique(data_in)) < 0.01 % --> check for bipolar + data_in = int32(data_in.*sqrt(5)); + data_in = double(data_in); + else % unipolar + data_in = int32(data_in.*sqrt(5)*2-3); % make bipolar [-3, -1, 1, 3] + data_in = double(data_in); + end + + data_oh = obj.oh_decider(data_in, obj.blocklength); % decider for overhead + + data_out = obj.mrds_decoding(data_oh, obj.blocklength); + data_out = data_out./sqrt(5); % normalized to Power=1 + + signalclass_in.signal = data_out'; + + % append to logbook + lbdesc = ['MRDS Coded']; + signalclass_in = signalclass_in.logbookentry(lbdesc); + + % write to output + signalclass_out = signalclass_in; + + end + + + function signalclass_out = dc_remove(obj,signalclass_in,options) + + arguments + obj + signalclass_in + options.oversampling_factor = 1; + end + + + signalclass_in.signal = obj.dcr(signalclass_in.signal, obj.blocklength, options.oversampling_factor); + + % append to logbook + lbdesc = ['MRDS DC Removed']; + signalclass_in = signalclass_in.logbookentry(lbdesc); + + % write to output + signalclass_out = signalclass_in; + + end + + end + + methods (Access=private) + % Cant be seen from outside! So put all your functions here that can/ + % shall not be called from outside + + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % Function 1 - encoding + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + function [data_out] = mrds_encoding(~,data_in, blocklength) + % data_in: bipolar PAM4 sequence with levels [-3, -1, 1, 3] + % with length power of two + % blocklength: power of two <= length of data_in + + oh_length = log2(blocklength); + data_out = zeros(1,length(data_in)+length(data_in)/blocklength*oh_length); + l = 0; + + for j = 1:blocklength:length(data_in) + + data = data_in(j:j+blocklength-1); + + z_N = zeros(1,blocklength); % RDS + z_rds = 0; + for i = 1:blocklength + z_rds = z_rds + data(i); + z_N(i) = z_rds; + end + + z = z_rds/2; % find inversion point k + k = find(z_N == z); + [~, index] = min(abs(blocklength/2 - k)); + k = k(index); + if isempty(k) == 1 + k = blocklength; + end + + if k == blocklength + overhead = ones(1,oh_length)*3; + else + overhead = (decimalToBinaryVector(k-1,log2(blocklength))-0.5)*6; % calculate OH + data(k+1:end) = data(k+1:end)*(-1); % invert + end + + data_oh = [data overhead]; + data_out(l*(blocklength+oh_length)+1:(l+1)*(blocklength+oh_length)) = data_oh; + l = l+1; + + end + + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % Function 2 - decoding + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + function [data_out] = mrds_decoding(~,data_in, blocklength) + + oh_length = log2(blocklength); + data_out = zeros(1,length(data_in)-length(data_in)/(blocklength+oh_length)*oh_length); + l = 1; + + for j = 1:(blocklength+oh_length):length(data_in) + + data = data_in(j:j+blocklength+oh_length-1); + + overhead = data(blocklength+1:end)/6+0.5; + k = binaryVectorToDecimal(overhead)+1; + + if k == blocklength + data = data(1:blocklength); + else + data = data(1:blocklength); + data(k+1:end) = data(k+1:end)*(-1); % invert + end + + data_out(l:l+blocklength-1) = data; + l = l+blocklength; + + end + + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % Function 3 - decider for overhead values + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + function [data_out] = oh_decider(~,data_in, blocklength) + + oh_length = log2(blocklength); + threshold = 0; + data_out = data_in; + + for j = 1:length(data_in)/(blocklength+oh_length) + for k = 1:oh_length + if data_in(j*blocklength+(j-1)*oh_length+k) >= threshold + data_out(j*blocklength+(j-1)*oh_length+k) = 3; + else + data_out(j*blocklength+(j-1)*oh_length+k) = -3; + end + end + end + + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % Function 4 - matched DC removal (DCR) + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + function [data_out] = dcr(~,data_in, blocklength, oversampling_factor) + + oh_length = log2(blocklength); + + winlength = blocklength+oh_length; + + if oversampling_factor > 1 + winlength = winlength*oversampling_factor; + end + + for j = 1:winlength:length(data_in) + + try + + data = data_in(j:j+winlength-1); + + rmean = mean(data); + + data_out(j:j+winlength-1) = data - rmean; + + catch + + if j+winlength > length(data_in) + data = data_in(j:length(data_in)); + else + error('indice problem.') + end + + rmean = mean(data); + + data_out(j:length(data_in)) = data - rmean; + + end + + + + end + + data_out = data_out'; + + end + + + + end +end diff --git a/Classes/04_DSP/VNLE.m b/Classes/04_DSP/VNLE.m index 60dc09d..0901fa5 100644 --- a/Classes/04_DSP/VNLE.m +++ b/Classes/04_DSP/VNLE.m @@ -3,6 +3,9 @@ classdef VNLE < handle % 1) Training mode (stable performance when you use NLMS) % 2) Decision directed mode + % Eq = VNLE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",[0.0004 0.0005 0.0006],"mu_tr",0,"order",[25,2,2],"sps",2,"decide",1); + % Somehow it is not possible to use only 1 nonlinear order + properties sps % usually 2 order diff --git a/Datatypes/coding_mode.m b/Datatypes/coding_mode.m new file mode 100644 index 0000000..f00000e --- /dev/null +++ b/Datatypes/coding_mode.m @@ -0,0 +1,8 @@ +classdef coding_mode < int32 + + enumeration + encode (1) + decode (2) + end + +end \ No newline at end of file diff --git a/projects/MPI_August/auswertung/clipfactor_curve.m b/projects/MPI_August/auswertung/clipfactor_curve.m new file mode 100644 index 0000000..ed22dd8 --- /dev/null +++ b/projects/MPI_August/auswertung/clipfactor_curve.m @@ -0,0 +1,65 @@ +M = wh.parameter.M.values(1); +datarate = wh.parameter.datarate.values(1); +sir = wh.parameter.sir.values(1); +laser_linewidth = wh.parameter.laser_linewidth.values(1); +pn_key = wh.parameter.pn_key.values(1); +rop = wh.parameter.rop.values; +cfs = wh.parameter.vbias_rel.values; + +cols = linspecer(6); +c_cnt = 0; +for c = cfs + c_cnt = c_cnt+1; + for pnk = pn_key + + %get ROP curve data for cur. sir and realization + + ber_ffe = wh.getStoValue('ber_ffe',M,datarate,sir,laser_linewidth,pnk,rop,c); + + try + rrop_ffe(c_cnt,1) = getIntersection(ber_ffe,rop); + end + + if 1 + % Create rop plot + figure(43); + hold on; % Retain the plot so new points can be added without complete redraw + plot(rop,ber_ffe',"LineWidth",1,"LineStyle","-","Marker",".","MarkerSize",10,"DisplayName","ber ffe",'Color',cols(5,:)); + yline(3.8e-3,'DisplayName','HD-FEC'); + xlabel('Signal to Interference Ratio (dB)'); + ylabel('Bit Error Rate (BER)'); + + set(gca,'yscale','log'); + grid on; + legend + end + + end +end + +figure(2243); +hold on; % Retain the plot so new points can be added without complete redraw +plot(cfs,mean(rrop_ffe,2,"omitnan")',"LineWidth",1,"LineStyle","--","Marker",".","MarkerSize",10,"DisplayName","FFE",'Color',cols(5,:)); + +xlabel('Clipfactor'); +ylabel('Receiver Sensitivity'); +%title(['Bit Error Rate vs. SIR; SIR: ',num2str(s),' dB']); +grid on; +legend + + + +function i = getIntersection(ber,rop) + +%get intersection between rop curve and hd-fec limit +hdfec = 3.8e-3 .* ones(size(ber)); +i = InterX([rop;hdfec'],[rop;ber']); + +if isempty(i) + i = NaN; +else + i = i(1); +end +end + + diff --git a/projects/MPI_August/auswertung/dc_removal_sir_curve.m b/projects/MPI_August/auswertung/dc_removal_sir_curve.m new file mode 100644 index 0000000..99ab0ab --- /dev/null +++ b/projects/MPI_August/auswertung/dc_removal_sir_curve.m @@ -0,0 +1,59 @@ +M = wh.parameter.M.values(1); +datarate = wh.parameter.datarate.values(1); +sirs = wh.parameter.sir.values; +laser_linewidths = wh.parameter.laser_linewidth.values; +pnoi = wh.parameter.pn_key.values; +rop = wh.parameter.rop.values; +winlengths = wh.parameter.winlength.values(1); + +cols = linspecer(numel(winlengths)); +figure(44); +hold on + + +cntlw = 0; + +for laser_linewidth = laser_linewidths + cntlw = cntlw+1; + cntsir = 0; + for sir = sirs + cntsir = cntsir+1; + for pnk = pnoi + + + + %cnt = cnt+1; + ber_ffe = wh.getStoValue('ber_ffe',M,datarate,sir,laser_linewidth,pnk,rop,winlength); + rrop_ffe(cntsir) = getIntersection(ber_ffe,rop); + + + end + end + + plot(sirs,rrop_ffe,'DisplayName',['Linewidth: ',num2str(laser_linewidth.*1e-6),' MHz'],"LineWidth",1,"LineStyle","-","Marker",".","MarkerSize",10) + xlabel('Signal to Interference Ratio (dB)'); + ylabel('Sensitivity in dBm'); + title('Receiver Sensitivity vs. SIR'); + set(gca,'Box','on'); + grid on; + grid minor + legend + +end + + + + + +function i = getIntersection(ber,rop) + +%get intersection between rop curve and hd-fec limit +hdfec = 3.8e-3 .* ones(size(ber)); +i = InterX([rop;hdfec'],[rop;ber']); + +if isempty(i) + i = NaN; +else + i = i(1); +end +end \ No newline at end of file diff --git a/projects/MPI_August/auswertung/dc_removal_window_length.m b/projects/MPI_August/auswertung/dc_removal_window_length.m new file mode 100644 index 0000000..e619384 --- /dev/null +++ b/projects/MPI_August/auswertung/dc_removal_window_length.m @@ -0,0 +1,71 @@ +M = wh.parameter.M.values(1); +datarate = wh.parameter.datarate.values(1); +sirs = wh.parameter.sir.values; +laser_linewidth = wh.parameter.laser_linewidth.values(1); +pnoi = wh.parameter.pn_key.values; +rop = wh.parameter.rop.values; +winlengths = wh.parameter.winlength.values; + +cols = linspecer(numel(winlengths)); +figure(44); +hold on + +cnt = 0; +for winlength = winlengths + cnt = cnt+1; + for sir = sirs + + for pnk = pnoi + + + + %cnt = cnt+1; + ber_ffe = wh.getStoValue('ber_ffe',M,datarate,sir,laser_linewidth,pnk,rop,winlength); + rrop_ffe(cnt,pnk) = getIntersection(ber_ffe,rop); + + + if sum(ber_ffe) + + plot(rop,ber_ffe',"LineWidth",0.5,"LineStyle","-","Marker",".","MarkerSize",15,"DisplayName",['win len: ',num2str(winlength), ' '],"Color",cols(cnt,:)); + yline(3.8e-3,'DisplayName','HD-FEC','LineStyle','--','HandleVisibility','off'); + xlabel('Received Optical Power (dBm)'); + ylabel('Bit Error Rate (BER)'); + title('Bit Error Rate vs. ROP'); + set(gca,'yscale','log'); + set(gca,'Box','on'); + grid on; + grid minor + legend + + end + + end + end + +end + +figure(123) +hold on +plot(winlengths,mean(rrop_ffe,2),'DisplayName',['Linewidth: 10 MHz'],"LineWidth",1,"LineStyle","-","Marker",".","MarkerSize",10) +xlabel('Window Length'); +ylabel('Sensitivity in dBm'); +title('Receiver Sensitivity vs. Window Length'); +set(gca,'Box','on'); +grid on; +grid minor +legend + + + +function i = getIntersection(ber,rop) + +%get intersection between rop curve and hd-fec limit +hdfec = 3.8e-3 .* ones(size(ber)); +i = InterX([rop;hdfec'],[rop;ber']); + +if isempty(i) + i = NaN; +else + i = i(1); +end +end \ No newline at end of file diff --git a/projects/MPI_August/auswertung/hpf_curve.m b/projects/MPI_August/auswertung/hpf_curve.m new file mode 100644 index 0000000..e0141e8 --- /dev/null +++ b/projects/MPI_August/auswertung/hpf_curve.m @@ -0,0 +1,76 @@ +M = wh.parameter.M.values(1); +datarate = wh.parameter.datarate.values(1); +sirs = wh.parameter.sir.values; +laser_linewidth = wh.parameter.laser_linewidth.values(1); +pnoi = wh.parameter.pn_key.values; +rop = wh.parameter.rop.values; +hpf = wh.parameter.hpf.values(2); + +cols = linspecer(16); + +for lw = laser_linewidth + cnt = 0; + for sir = sirs + cnt = cnt+1; + for pnk = pnoi + + + + %cnt = cnt+1; + ber_ffe = wh.getStoValue('ber_ffe',M,datarate,sir,lw,pnk,rop,fc); + rrop_ffe(cnt,pnk) = getIntersection(ber_ffe,rop); + + % if 0 + % % Create the initial plot + % figure(44); + % a = gca; + % cnt = numel(a.Children)+1; + % hold on; % Retain the plot so new points can be added without complete redraw + % end + + + + if sum(ber_ffe) + figure(44); + plot(rop,ber_ffe',"LineWidth",0.5,"LineStyle","-","Marker",".","MarkerSize",15,"DisplayName","FFE only",'Color',cols(cnt,:)); + yline(3.8e-3,'DisplayName','HD-FEC','LineStyle','--','HandleVisibility','off'); + xlabel('Received Optical Power (dBm)'); + ylabel('Bit Error Rate (BER)'); + title('Bit Error Rate vs. ROP'); + set(gca,'yscale','log'); + set(gca,'Box','on'); + grid on; + grid minor + legend + + end + + end + end + + + figure(123) + hold on + plot(sirs,mean(rrop_ffe,2),'DisplayName',['lw: ',num2str(lw), ' MHz'],"LineWidth",1,"LineStyle","-","Marker",".","MarkerSize",10) + xlabel('Signal to Interference (dBm)'); + ylabel('Sensitivity in dBm'); + title('Receiver Sensitivity vs. SIR'); + set(gca,'Box','on'); + grid on; + grid minor + legend + +end + +function i = getIntersection(ber,rop) + +%get intersection between rop curve and hd-fec limit +hdfec = 3.8e-3 .* ones(size(ber)); +i = InterX([rop;hdfec'],[rop;ber']); + +if isempty(i) + i = NaN; +else + i = i(1); +end +end \ No newline at end of file diff --git a/projects/MPI_August/auswertung/hpf_curve_different_filter_bw.m b/projects/MPI_August/auswertung/hpf_curve_different_filter_bw.m new file mode 100644 index 0000000..1b9838a --- /dev/null +++ b/projects/MPI_August/auswertung/hpf_curve_different_filter_bw.m @@ -0,0 +1,59 @@ +M = wh.parameter.M.values(1); +datarate = wh.parameter.datarate.values(1); +sirs = wh.parameter.sir.values; +laser_linewidth = wh.parameter.laser_linewidth.values(1); +pnoi = wh.parameter.pn_key.values; +rop = wh.parameter.rop.values; +hpfs = wh.parameter.hpf.values; + +cols = linspecer(16); +figure(44); +hold on +for hpf = hpfs + cnt = 0; + for sir = sirs + cnt = cnt+1; + for pnk = pnoi + + + + %cnt = cnt+1; + ber_ffe = wh.getStoValue('ber_ffe',M,datarate,sir,lw,pnk,rop,hpf); + % rrop_ffe(cnt,pnk) = getIntersection(ber_ffe,rop); + + + if sum(ber_ffe) + + plot(rop,ber_ffe',"LineWidth",0.5,"LineStyle","-","Marker",".","MarkerSize",15,"DisplayName",['fc: ',num2str(hpf.*1e-6), ' MHz']); + yline(3.8e-3,'DisplayName','HD-FEC','LineStyle','--','HandleVisibility','off'); + xlabel('Received Optical Power (dBm)'); + ylabel('Bit Error Rate (BER)'); + title('Bit Error Rate vs. ROP'); + set(gca,'yscale','log'); + set(gca,'Box','on'); + grid on; + grid minor + legend + + end + + end + end + + + + +end + +function i = getIntersection(ber,rop) + +%get intersection between rop curve and hd-fec limit +hdfec = 3.8e-3 .* ones(size(ber)); +i = InterX([rop;hdfec'],[rop;ber']); + +if isempty(i) + i = NaN; +else + i = i(1); +end +end \ No newline at end of file diff --git a/projects/MPI_August/auswertung/mrds_sir_curve.m b/projects/MPI_August/auswertung/mrds_sir_curve.m new file mode 100644 index 0000000..5911ede --- /dev/null +++ b/projects/MPI_August/auswertung/mrds_sir_curve.m @@ -0,0 +1,59 @@ +M = wh.parameter.M.values(1); +datarate = wh.parameter.datarate.values(1); +sirs = wh.parameter.sir.values; +laser_linewidths = wh.parameter.laser_linewidth.values; +pnoi = wh.parameter.pn_key.values; +rop = wh.parameter.rop.values; +winlengths = wh.parameter.winlength.values(1); + +cols = linspecer(numel(winlengths)); +figure(44); +hold on + + +cntlw = 0; + +for laser_linewidth = laser_linewidths + cntlw = cntlw+1; + cntsir = 0; + for sir = sirs + cntsir = cntsir+1; + for pnk = pnoi + + + + %cnt = cnt+1; + ber_ffe = wh.getStoValue('ber_ffe',M,datarate,sir,laser_linewidth,pnk,rop,winlengths); + rrop_ffe(cntsir) = getIntersection(ber_ffe,rop); + + + end + end + + plot(sirs,rrop_ffe,'DisplayName',['Linewidth: ',num2str(laser_linewidth.*1e-6),' MHz'],"LineWidth",1,"LineStyle","-","Marker",".","MarkerSize",10) + xlabel('Signal to Interference Ratio (dB)'); + ylabel('Sensitivity in dBm'); + title('Receiver Sensitivity vs. SIR'); + set(gca,'Box','on'); + grid on; + grid minor + legend + +end + + + + + +function i = getIntersection(ber,rop) + +%get intersection between rop curve and hd-fec limit +hdfec = 3.8e-3 .* ones(size(ber)); +i = InterX([rop;hdfec'],[rop;ber']); + +if isempty(i) + i = NaN; +else + i = i(1); +end +end \ No newline at end of file diff --git a/projects/MPI_August/auswertung/rop_curve.m b/projects/MPI_August/auswertung/rop_curve.m index c768bd7..02e898e 100644 --- a/projects/MPI_August/auswertung/rop_curve.m +++ b/projects/MPI_August/auswertung/rop_curve.m @@ -5,31 +5,48 @@ laser_linewidth = wh.parameter.laser_linewidth.values(1); pn_key = wh.parameter.pn_key.values; rop = wh.parameter.rop.values; -cols = linspecer(4); +cols = linspecer(8); cnt = 0; for pnk = pn_key - cnt = cnt+1; + %cnt = cnt+1; ber_ffe = wh.getStoValue('ber_ffe',M,datarate,sir,laser_linewidth,pnk,rop); ber_dcavg = wh.getStoValue('ber_dcavg',M,datarate,sir,laser_linewidth,pnk,rop); ber_adapt = wh.getStoValue('ber_adapt',M,datarate,sir,laser_linewidth,pnk,rop); ber_derem = wh.getStoValue('ber_dcrem',M,datarate,sir,laser_linewidth,pnk,rop); % Create the initial plot - figure(43); + figure(44); + a = gca; + cnt = numel(a.Children)+1; hold on; % Retain the plot so new points can be added without complete redraw - plot(rop,ber_ffe',"LineWidth",1,"LineStyle","-","Marker",".","MarkerSize",10,"DisplayName","ber ffe",'Color',cols(1,:)); - plot(rop,ber_dcavg',"LineWidth",1,"LineStyle","-","Marker",".","MarkerSize",10,"DisplayName","ber dcavg",'Color',cols(2,:)); - plot(rop,ber_adapt',"LineWidth",1,"LineStyle","-","Marker",".","MarkerSize",10,"DisplayName","ber adapt",'Color',cols(3,:)); - plot(rop,ber_derem',"LineWidth",1,"LineStyle","-","Marker",".","MarkerSize",10,"DisplayName","ber derem",'Color',cols(4,:)); + + if sum(ber_ffe) + plot(rop,ber_ffe',"LineWidth",0.5,"LineStyle","-","Marker",".","MarkerSize",15,"DisplayName","FFE only",'Color',cols(cnt,:)); + end + + if sum(ber_dcavg) + plot(rop,ber_dcavg',"LineWidth",0.5,"LineStyle","-","Marker",".","MarkerSize",15,"DisplayName","FFE + FF DC Avg.",'Color',cols(cnt+1,:)); + end + + if sum(ber_adapt) + plot(rop,ber_adapt',"LineWidth",0.5,"LineStyle","-","Marker",".","MarkerSize",15,"DisplayName","FFE + Adaptive Levels",'Color',cols(cnt+2,:)); + end + + if sum(ber_derem) + plot(rop,ber_derem',"LineWidth",0.5,"LineStyle","-","Marker",".","MarkerSize",15,"DisplayName","FFE + Adaptive DC-Subtraction",'Color',cols(cnt+3,:)); + end end -yline(3.8e-3,'DisplayName','HD-FEC'); -xlabel('Signal to Interference Ratio (dB)'); +yline(3.8e-3,'DisplayName','HD-FEC','LineStyle','--','HandleVisibility','off'); +xlabel('Received Optical Power (dBm)'); ylabel('Bit Error Rate (BER)'); title('Bit Error Rate vs. ROP'); set(gca,'yscale','log'); +set(gca,'Box','on'); grid on; +grid minor legend + diff --git a/projects/MPI_August/auswertung/sensitivity_curve.m b/projects/MPI_August/auswertung/sensitivity_curve.m new file mode 100644 index 0000000..6c99780 --- /dev/null +++ b/projects/MPI_August/auswertung/sensitivity_curve.m @@ -0,0 +1,81 @@ +M = wh.parameter.M.values(1); +datarate = wh.parameter.datarate.values(1); +sir = wh.parameter.sir.values(4); +laser_linewidths = wh.parameter.laser_linewidth.values; +pn_key = wh.parameter.pn_key.values; +rop = wh.parameter.rop.values; + +cols = linspecer(4); + +for laser_linewidth = laser_linewidths + s_cnt = 0; + + for s = sir + s_cnt = s_cnt+1; + for pnk = pn_key + + %get ROP curve data for cur. sir and realization + + ber_ffe = wh.getStoValue('ber_ffe',M,datarate,s,laser_linewidth,pnk,rop); + ber_dcavg = wh.getStoValue('ber_dcavg',M,datarate,s,laser_linewidth,pnk,rop); + ber_adapt = wh.getStoValue('ber_adapt',M,datarate,s,laser_linewidth,pnk,rop); + ber_derem = wh.getStoValue('ber_dcrem',M,datarate,s,laser_linewidth,pnk,rop); + + rrop_ffe(s_cnt,pnk) = getIntersection(ber_ffe,rop); + rrop_dcavg(s_cnt,pnk) = getIntersection(ber_dcavg,rop); + rrop_adapt(s_cnt,pnk) = getIntersection(ber_adapt,rop); + rrop_derem(s_cnt,pnk) = getIntersection(ber_derem,rop); + + if 1 + % Create rop plot + figure(43); + hold on; % Retain the plot so new points can be added without complete redraw + plot(rop,ber_ffe',"LineWidth",1,"LineStyle","-","Marker",".","MarkerSize",10,"DisplayName","ber ffe",'Color',cols(1,:)); + % scatter(rrop_ffe(s_cnt,pnk),3.8e-3,'LineWidth',2); + plot(rop,ber_dcavg',"LineWidth",1,"LineStyle","-","Marker",".","MarkerSize",10,"DisplayName","ber dcavg",'Color',cols(2,:)); + plot(rop,ber_adapt',"LineWidth",1,"LineStyle","-","Marker",".","MarkerSize",10,"DisplayName","ber adapt",'Color',cols(3,:)); + plot(rop,ber_derem',"LineWidth",1,"LineStyle","-","Marker",".","MarkerSize",10,"DisplayName","ber dcrem",'Color',cols(4,:)); + yline(3.8e-3,'DisplayName','HD-FEC'); + xlabel('Signal to Interference Ratio (dB)'); + ylabel('Bit Error Rate (BER)'); + title(['Bit Error Rate vs. SIR; SIR: ',num2str(s),' dB']); + set(gca,'yscale','log'); + grid on; + legend + end + + + end + end + + figure(223); + hold on; % Retain the plot so new points can be added without complete redraw + % plot(sir,mean(rrop_ffe,2,"omitnan")',"LineWidth",1,"LineStyle",":","Marker",".","MarkerSize",10,"DisplayName","FFE",'Color',cols(1,:)); + % plot(sir,mean(rrop_dcavg,2,"omitnan")',"LineWidth",1,"LineStyle","-","Marker",".","MarkerSize",10,"DisplayName",['Linewidth: ',num2str(laser_linewidth.*1e-6),' MHz'],'Color',cols(2,:)); + % plot(sir,mean(rrop_adapt,2,"omitnan")',"LineWidth",1,"LineStyle","-","Marker",".","MarkerSize",10,"DisplayName",['Linewidth: ',num2str(laser_linewidth.*1e-6),' MHz']); + plot(sir,mean(rrop_derem,2,"omitnan")',"LineWidth",1,"LineStyle","-","Marker",".","MarkerSize",10,"DisplayName",['Linewidth: ',num2str(laser_linewidth.*1e-6),' MHz']); + xlabel('Signal to Interference Ratio (dB)'); + ylabel('Receiver Sensitivity'); + title(['Bit Error Rate vs. SIR']); + grid on; + legend + + + +end + + +function i = getIntersection(ber,rop) + +%get intersection between rop curve and hd-fec limit +hdfec = 3.8e-3 .* ones(size(ber)); +i = InterX([rop;hdfec'],[rop;ber']); + +if isempty(i) + i = NaN; +else + i = i(1); +end +end + + diff --git a/projects/MPI_August/auswertung/sir_curve.m b/projects/MPI_August/auswertung/sir_curve.m index 60eaa03..ad6dc9a 100644 --- a/projects/MPI_August/auswertung/sir_curve.m +++ b/projects/MPI_August/auswertung/sir_curve.m @@ -17,7 +17,7 @@ for pnk = pn_key ber_derem = wh.getStoValue('ber_dcrem',M,datarate,sir,laser_linewidth,pnk,rop); % Create the initial plot - figure(43); + figure(44); hold on; % Retain the plot so new points can be added without complete redraw plot(sir,ber_ffe',"LineWidth",1,"LineStyle","-","Marker",".","MarkerSize",10,"DisplayName","ber ffe",'Color',cols(1,:)); plot(sir,ber_dcavg',"LineWidth",1,"LineStyle","-","Marker",".","MarkerSize",10,"DisplayName","ber dcavg",'Color',cols(2,:)); diff --git a/projects/MPI_August/mpi_uncorrelated.m b/projects/MPI_August/mpi_uncorrelated.m index 685b870..ed50896 100644 --- a/projects/MPI_August/mpi_uncorrelated.m +++ b/projects/MPI_August/mpi_uncorrelated.m @@ -4,16 +4,20 @@ params = struct; params.M = [4]; params.datarate = [224]; -params.sir = [15:45]; %decibel = attenuation of interference path -params.laser_linewidth = [1e6]; +params.sir = [35]; %decibel = attenuation of interference path +params.laser_linewidth = [10e6]; + params.pn_key = [1]; -params.rop = [-7]; +params.rop = [-12:1:-1]; + + + +usemrds = 0; name = ['wh_',strrep(num2str(now),'.','')]; wh = DataStorage(params); - wh.addStorage("ber_ffe"); wh.addStorage("ber_dcavg"); wh.addStorage("ber_adapt"); @@ -22,74 +26,96 @@ wh.addStorage("ber_dcrem"); %% Init Params link_length = 10000; %meter + + endcnt = prod(wh.dim); cnt=0; disp(['Start Simulation of ',num2str(endcnt),' loops...']) tic + + + + + for M = wh.parameter.M.values for datarate = wh.parameter.datarate.values for pn_key = wh.parameter.pn_key.values - - %% SETUP HERE: %% + % SETUP HERE: %% kover = 8; M8199 = M8199A("kover",kover); fdac = M8199.fdac; - fsym = round(datarate / log2(M))*1e9; + fsym = round(datarate / log2(M)) * 1e9; rrcalpha = 0.05; Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"rrcalpha",rrcalpha); % MAIN SIGNAL + %%%%% Symbol Generation %%%%%% - [Digi_sig,Symbols,Bits] = PAMsource("fsym",fsym,"M",M,"order",18,"useprbs",0,"fs_out",M8199.fdac,"applyclipping",1,"clipfactor",1.3,"applypulseform",1,"pulseformer",Pform,"randkey",pn_key).process(); + [Digi_sig,Symbols,Bits] = PAMsource("fsym",fsym,"M",M,"order",18,"useprbs",0,... + "fs_out",M8199.fdac,"applyclipping",1,"clipfactor",1.4,"applypulseform",1,"pulseformer",Pform,"randkey",pn_key,"mrds_code",usemrds,"mrds_blocklength",wl).process(); + %%%%% AWG %%%%%% El_sig = M8199.process(Digi_sig); - El_sig = El_sig.*0.7222; - El_sig.signal = awgn(El_sig.signal,20,'measured',pn_key); - %%%%% Lowpass before Modulator %%%%%% - El_sig = Filter('filtdegree',2,"f_cutoff",60e9,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(El_sig); - El_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",10).process(El_sig); + %El_sig.signal = awgn(El_sig.signal,20,'measured',pn_key); + %%%%% Lowpass el. components %%%%%% + El_sig = Filter('filtdegree',2,"f_cutoff",60e9,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(El_sig); + + %%%%% Electrical Driver Amplifier %%%%%% + El_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",9).process(El_sig); + + fprintf('Driver output power: %s dBm\n', num2str(El_sig.power)); + fprintf('Driver output peak voltage: %s Vpp \n', num2str(max(El_sig.signal)-min(El_sig.signal))); % INTERFERENCE SIGNAL + %%%%% Symbol Generation %%%%%% - [Digi_sig_i,Symbols_i,Bits_i] = PAMsource("fsym",fsym,"M",M,"order",18,"useprbs",0,"fs_out",M8199.fdac,"applyclipping",1,"clipfactor",1.3,"applypulseform",1,"pulseformer",Pform,"randkey",pn_key*2).process(); + [Digi_sig_i,Symbols_i,Bits_i] = PAMsource("fsym",fsym,"M",M,"order",18,"useprbs",0,... + "fs_out",M8199.fdac,"applyclipping",1,"clipfactor",1.4,"applypulseform",1,"pulseformer",Pform,"randkey",pn_key*2,"mrds_code",usemrds,"mrds_blocklength",wl).process(); + %%%%% AWG %%%%%% El_sig_i = M8199.process(Digi_sig_i); - El_sig_i = El_sig_i.*0.7222; - El_sig_i.signal = awgn(El_sig_i.signal,20,'measured',pn_key*2); + %El_sig_i.signal = awgn(El_sig_i.signal,20,'measured',pn_key*2); + %%%%% Lowpass before Modulator %%%%%% El_sig_i = Filter('filtdegree',2,"f_cutoff",60e9,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(El_sig_i); - El_sig_i = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",10).process(El_sig_i); + El_sig_i = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",9).process(El_sig_i); + fprintf('Driver output power: %s dBm\n', num2str(El_sig.power)); + fprintf('Driver output peak voltage: %s Vpp \n', num2str(max(El_sig.signal)-min(El_sig.signal))); + for laser_linewidth = wh.parameter.laser_linewidth.values % MAIN SIGNAL %%%%% MODULATE E/O CONVERSION %%%%%% - vbias_rel = 0.5; + vbias_rel = 0.6; u_pi = 2.9; vbias = -vbias_rel*u_pi; + [Opt_sig] = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs,"lambda",1290,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",pn_key).process(El_sig); + Optfilter = Filter('filtdegree',6,"f_cutoff",fsym.*0.7,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true); Opt_sig = Optfilter.process(Opt_sig); - er = Opt_sig.extinctionratio(fsym,M); + Opt_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",0).process(Opt_sig); % INTERFERENCE SIGNAL %%%%% MODULATE E/O CONVERSION %%%%%% - u_pi = 2.9; - vbias = -vbias_rel*u_pi; + [Opt_sig_i] = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs,"lambda",1290,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",pn_key+1).process(El_sig_i); + Opt_sig_i = Optfilter.process(Opt_sig_i); + Opt_sig_i = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",0).process(Opt_sig_i); j_ = wh.parameter.sir.length; i_ = wh.parameter.rop.length; - + ber_ffe=zeros(j_,i_); ber_dcavg=zeros(j_,i_); ber_adapt=zeros(j_,i_); @@ -97,8 +123,8 @@ for M = wh.parameter.M.values patten=zeros(j_,i_); - parfor j = 1:j_ - + for j = 1:j_ + sir = wh.parameter.sir.values(j); @@ -117,7 +143,8 @@ for M = wh.parameter.M.values % Opt = channel_model_mpi(Opt_sig,link_length,mpi_path,sir); % Receiver ROP curve - for i = 1:i_ + + parfor i = 1:i_ rop=wh.parameter.rop.values(i); % Set ROP @@ -129,47 +156,132 @@ for M = wh.parameter.M.values %%%%%% Lowpass PhDiode %%%%%% Rx_sig = Filter('filtdegree',2,"f_cutoff",70e9,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true).process(Rx_sig); - % Rx_sig.spectrum("displayname","Received Signal after PhD","fignum",201); + + % Rx_sig.signal = Rx_sig.signal - mean(Rx_sig.signal); + % + % Rx_sig = Rx_sig.setPower(0,"dBm"); + + if 0 + [b, a] = butter(1, 0.1e9/(Rx_sig.fs/2), 'high'); + % Rx_sig.spectrum("fignum",1111,"displayname","before HPF"); + Rx_sig.signal = filter(b, a, Rx_sig.signal); + % Rx_sig.spectrum("fignum",1111,"displayname","after HPF"); + end + + % figure(111) + % hold on + % [h, f] = freqz(b, a, length(Rx_sig), Rx_sig.fs); + % plot(f, 20*log10(abs(h))); + % title('Magnitude Response of Digital High-Pass Filter'); + % xlabel('Frequency (Hz)'); + % ylabel('Magnitude (dB)'); + % grid on; %%%%%% Scope %%%%%% fadc = 256e9; - Lp_scpe = Filter('filtdegree',4,"f_cutoff",63e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true); + Lp_scpe = Filter('filtdegree',4,"f_cutoff",100e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true); Scpe_sig = Scope("fsimu",fdac*kover,"fadc",fadc,... "delay",0,"fixed_delay",0,"filtertype",filtertypes.butterworth,... "samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,... - "adcresolution",16,"quantbuffer",0.1,'block_dc',1,'lpf_active',1,'H_lpf',Lp_scpe).process(Rx_sig); - + "adcresolution",10,"quantbuffer",0.1,'block_dc',1,'lpf_active',1,'H_lpf',Lp_scpe).process(Rx_sig); + + + % Scpe_sig.spectrum("displayname","Received Signal after Scope","fignum",201); + %%%%%% Sample to 2x fsym %%%%%% Scpe_sig = Scpe_sig.resample("fs_in",fadc,"fs_out",2*fsym); + if usemrds + Scpe_sig = MRDS_coding("blocklength",512).dc_remove(Scpe_sig,"oversampling_factor",2); + end + % + + rmean = zeros(1,2); + if 0 + + data_in = Scpe_sig.signal; + winlength = wl; + + for k = 1:winlength:length(data_in) + + try + + data = data_in(k:k+winlength-1); + + rmean(1) = mean(data); + rmean = circshift(rmean,1); + Scpe_sig.signal(k:k+winlength-1) = data - rmean(1); + + catch + + if k+winlength > length(data_in) + data = data_in(k:length(data_in)); + else + error('indice problem.') + end + + rmean(1) = mean(data); + + Scpe_sig.signal(k:length(data_in)) = data - rmean(1); + + end + end + % Scpe_sig.signal = data_out; + end + + + % Scpe_sig.plot("fignum",313,"displayname",'after dc removal'); + %%%%%% Sync Rx signal with reference %%%%%% [Scpe_sig,D,cuts] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym); + %%%%% EQUALIZE %%%%%% - Eq = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",1); - [EQ_sig] = Eq.process(Scpe_sig,Symbols); - Rx_bits = PAMmapper(M,0).demap(EQ_sig); - [~,errors_bm,ber_ffe(j,i),errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - disp(['BER: ',sprintf('%.1E',ber_ffe(j,i)),' - - ROP: ',num2str(patten(j,i)),'dBm - - PAM-',num2str(M),' - - ',num2str(fsym*1e-9),' GBd']); - Eq = FFE_FFDCAVG("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",1,"mu_buff",0.7); - [EQ_sig] = Eq.process(Scpe_sig,Symbols); - Rx_bits = PAMmapper(M,0).demap(EQ_sig); - [~,errors_bm,ber_dcavg(j,i),errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - disp(['BER: ',sprintf('%.1E',ber_dcavg(j,i)),' - - ROP: ',num2str(patten(j,i)),'dBm - - PAM-',num2str(M),' - - ',num2str(fsym*1e-9),' GBd']); + if ~usemrds + % + Eq = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",1); + % Eq = VNLE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",[0.0004 0.0005 0.0006],"mu_tr",0,"order",[25,2,2],"sps",2,"decide",1); + [EQ_sig] = Eq.process(Scpe_sig,Symbols); + Rx_bits = PAMmapper(M,0).demap(EQ_sig); + [~,errors_bm,ber_ffe(j,i),errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + disp(['BER: ',sprintf('%.1E',ber_ffe(j,i)),' - - ROP: ',num2str(patten(j,i)),'dBm - - PAM-',num2str(M),' - - ',num2str(fsym*1e-9),' GBd']); + % + % % + % Eq = FFE_FFDCAVG("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",1,"mu_buff",0.7); + % [EQ_sig] = Eq.process(Scpe_sig,Symbols); + % Rx_bits = PAMmapper(M,0).demap(EQ_sig); + % [~,errors_bm,ber_dcavg(j,i),errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + % disp(['BER: ',sprintf('%.1E',ber_dcavg(j,i)),' - - ROP: ',num2str(patten(j,i)),'dBm - - PAM-',num2str(M),' - - ',num2str(fsym*1e-9),' GBd']); + % + % + % Eq = FFE_adaptive_decision("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",1,"buffer_length",112); + % [EQ_sig] = Eq.process(Scpe_sig,Symbols); + % Rx_bits = PAMmapper(M,0).demap(EQ_sig); + % [~,errors_bm,ber_adapt(j,i),errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + % disp(['BER: ',sprintf('%.1E',ber_adapt(j,i)),' - - ROP: ',num2str(patten(j,i)),'dBm - - PAM-',num2str(M),' - - ',num2str(fsym*1e-9),' GBd']); + % + % + % Eq = FFE_DCremoval("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",1,"mu_dc",0.07,"dc_buffer_len",112); + % [EQ_sig] = Eq.process(Scpe_sig,Symbols); + % Rx_bits = PAMmapper(M,0).demap(EQ_sig); + % [~,errors_bm,ber_dcrem(j,i),errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + % disp(['BER: ',sprintf('%.1E',ber_dcrem(j,i)),' - - ROP: ',num2str(patten(j,i)),'dBm - - PAM-',num2str(M),' - - ',num2str(fsym*1e-9),' GBd']); + + + else + Eq = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",1); + + [EQ_sig] = Eq.process(Scpe_sig,Symbols); + Rx_bits = MRDS_coding("blocklength",wl).decode(EQ_sig); + Rx_bits = PAMmapper(M,0).demap(Rx_bits); + [~,errors_bm,ber_ffe(j,i),errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + disp(['BER: ',sprintf('%.1E',ber_ffe(j,i)),' - - ROP: ',num2str(patten(j,i)),'dBm - - PAM-',num2str(M),' - - ',num2str(fsym*1e-9),' GBd']); + + end - Eq = FFE_adaptive_decision("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",1,"buffer_length",85); - [EQ_sig] = Eq.process(Scpe_sig,Symbols); - Rx_bits = PAMmapper(M,0).demap(EQ_sig); - [~,errors_bm,ber_adapt(j,i),errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - disp(['BER: ',sprintf('%.1E',ber_adapt(j,i)),' - - ROP: ',num2str(patten(j,i)),'dBm - - PAM-',num2str(M),' - - ',num2str(fsym*1e-9),' GBd']); - Eq = FFE_DCremoval("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",1,"mu_dc",0.07); - [EQ_sig] = Eq.process(Scpe_sig,Symbols); - Rx_bits = PAMmapper(M,0).demap(EQ_sig); - [~,errors_bm,ber_dcrem(j,i),errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - disp(['BER: ',sprintf('%.1E',ber_dcrem(j,i)),' - - ROP: ',num2str(patten(j,i)),'dBm - - PAM-',num2str(M),' - - ',num2str(fsym*1e-9),' GBd']); end end @@ -179,7 +291,6 @@ for M = wh.parameter.M.values for i = 1:i_ rop=wh.parameter.rop.values(i); - wh.addValueToStorage(ber_ffe(j,i) ,'ber_ffe',M,datarate,sir,laser_linewidth,pn_key,rop); wh.addValueToStorage(ber_dcavg(j,i),'ber_dcavg',M,datarate,sir,laser_linewidth,pn_key,rop); wh.addValueToStorage(ber_adapt(j,i),'ber_adapt',M,datarate,sir,laser_linewidth,pn_key,rop); @@ -188,12 +299,13 @@ for M = wh.parameter.M.values end end - + toc disp(['Simulated: ',num2str(cnt/endcnt*100),' %']); wh.save('C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\MPI_August\auswertung\') end + end end end @@ -201,6 +313,8 @@ end +rop_curve; + diff --git a/projects/MPI_August/simualte_MPI_data.m b/projects/MPI_August/simualte_MPI_data.m new file mode 100644 index 0000000..0818bd2 --- /dev/null +++ b/projects/MPI_August/simualte_MPI_data.m @@ -0,0 +1,117 @@ + +usemrds = 0; + +%% Init Params +link_length = 10000; %meter + +Optfilter = Filter('filtdegree',6,"f_cutoff",fsym.*0.7,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true); + +tic + +% SETUP HERE: %% +kover = 8; +M8199 = M8199A("kover",kover); +fdac = M8199.fdac; +fsym = round(datarate / log2(M)) * 1e9; +rrcalpha = 0.05; +Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"rrcalpha",rrcalpha); + +% MAIN SIGNAL + +%%%%% Symbol Generation %%%%%% +[Digi_sig,Symbols,Bits] = PAMsource("fsym",fsym,"M",M,"order",18,"useprbs",0,... + "fs_out",M8199.fdac,"applyclipping",1,"clipfactor",1.4,"applypulseform",1,"pulseformer",Pform,"randkey",pn_key,"mrds_code",usemrds).process(); + +%%%%% AWG %%%%%% +El_sig = M8199.process(Digi_sig); +%El_sig.signal = awgn(El_sig.signal,20,'measured',pn_key); + +%%%%% Lowpass el. components %%%%%% +El_sig = Filter('filtdegree',2,"f_cutoff",60e9,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(El_sig); + +%%%%% Electrical Driver Amplifier %%%%%% +El_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",9).process(El_sig); + +fprintf('Driver output power: %s dBm\n', num2str(El_sig.power)); +fprintf('Driver output peak voltage: %s Vpp \n', num2str(max(El_sig.signal)-min(El_sig.signal))); + +% INTERFERENCE SIGNAL + +%%%%% Symbol Generation %%%%%% +[Digi_sig_i,Symbols_i,Bits_i] = PAMsource("fsym",fsym,"M",M,"order",18,"useprbs",0,... + "fs_out",M8199.fdac,"applyclipping",1,"clipfactor",1.4,"applypulseform",1,"pulseformer",Pform,"randkey",pn_key*2,"mrds_code",usemrds).process(); + +%%%%% AWG %%%%%% +El_sig_i = M8199.process(Digi_sig_i); +%El_sig_i.signal = awgn(El_sig_i.signal,20,'measured',pn_key*2); + +%%%%% Lowpass before Modulator %%%%%% +El_sig_i = Filter('filtdegree',2,"f_cutoff",60e9,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(El_sig_i); +El_sig_i = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",9).process(El_sig_i); + + +fprintf('Driver output power: %s dBm\n', num2str(El_sig.power)); +fprintf('Driver output peak voltage: %s Vpp \n', num2str(max(El_sig.signal)-min(El_sig.signal))); + + +% MAIN SIGNAL +%%%%% MODULATE E/O CONVERSION %%%%%% +vbias_rel = 0.6; +u_pi = 2.9; +vbias = -vbias_rel*u_pi; + +[Opt_sig] = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs,"lambda",1290,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",pn_key).process(El_sig); + +Opt_sig = Optfilter.process(Opt_sig); + +Opt_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",0).process(Opt_sig); + +% INTERFERENCE SIGNAL +%%%%% MODULATE E/O CONVERSION %%%%%% + +[Opt_sig_i] = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs,"lambda",1290,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",pn_key+1).process(El_sig_i); + +Opt_sig_i = Optfilter.process(Opt_sig_i); + +Opt_sig_i = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",0).process(Opt_sig_i); + + + + +for j = 1:j_ + + sir = wh.parameter.sir.values(j); + + %%%%% Interference Signal Fiber Prop %%%%%% + Opt_sig_i_prop = Fiber("fsimu",Opt_sig_i.fs,"fiber_length",link_length/1000,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.07).process(Opt_sig_i); + + Opt_sig_i_prop = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",Opt_sig.power-sir).process(Opt_sig_i_prop); + + %%%%% ADD Interference and Main Signal %%%%%% + Opt_sig_comb = Opt_sig_i_prop + Opt_sig; + + %%%%% Interference Signal Fiber Prop %%%%%% + Opt_sig_comb = Fiber("fsimu",Opt_sig_comb.fs,"fiber_length",link_length/1000,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.07).process(Opt_sig_comb); + + % Receiver ROP curve + for i = 1:i_ + rop=wh.parameter.rop.values(i); + + % Set ROP + Rx_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",rop).process(Opt_sig_comb); + + %%%%%% Square Law %%%%%% + Rx_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20,"nep",1.8e-11).process(Rx_sig); + + %%%%%% Lowpass PhDiode %%%%%% + Rx_sig = Filter('filtdegree',2,"f_cutoff",70e9,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true).process(Rx_sig); + + end +end + + +rop_curve; + + + +