From d8b4d6fe7c45a1bd9979a42d0cb9c1ec02044740 Mon Sep 17 00:00:00 2001 From: sioe Date: Tue, 15 Oct 2024 10:16:35 +0200 Subject: [PATCH] small updates from work pc try to implement kalman filter for MPI mitigation --- Classes/00_signals/Signal.m | 2 +- Classes/04_DSP/Equalizer/EQ.m | 7 +- Classes/04_DSP/Equalizer/FFE.m | 33 --- Classes/04_DSP/Equalizer/FFE_DCremoval.m | 14 +- Classes/04_DSP/Equalizer/FFE_FFDCAVG.m | 35 --- Classes/04_DSP/Equalizer/FFE_Kalman.m | 245 +++++++++++++++++ .../04_DSP/Equalizer/FFE_Kalman_Feedback.m | 197 ++++++++++++++ .../04_DSP/Equalizer/FFE_adaptive_decision.m | 13 +- .../400G_FTN_setups/imdd_dsp_approaches.m | 14 +- .../400G_FTN_setups/imdd_dsp_lab_copy_200G.m | 247 ++++++++++++++++++ projects/400G_FTN_setups/imdd_mpi_dsp.m | 73 ++++-- projects/MPI_August/mpi_uncorrelated.m | 40 ++- 12 files changed, 799 insertions(+), 121 deletions(-) create mode 100644 Classes/04_DSP/Equalizer/FFE_Kalman.m create mode 100644 Classes/04_DSP/Equalizer/FFE_Kalman_Feedback.m create mode 100644 projects/400G_FTN_setups/imdd_dsp_lab_copy_200G.m diff --git a/Classes/00_signals/Signal.m b/Classes/00_signals/Signal.m index e955679..77cd287 100644 --- a/Classes/00_signals/Signal.m +++ b/Classes/00_signals/Signal.m @@ -284,7 +284,7 @@ classdef Signal N = 2^(nextpow2(length(obj.signal))-8); [p_lin,w] = pwelch(obj.signal,hanning(N),N/2,N,obj.fs,"centered","power","mean"); - normalize = 1; + normalize = 0; if normalize p_lin = p_lin./ max(p_lin); p_dbm = 10*log10(p_lin); %dB to dBm in case of "power" diff --git a/Classes/04_DSP/Equalizer/EQ.m b/Classes/04_DSP/Equalizer/EQ.m index 888f85a..a92f827 100644 --- a/Classes/04_DSP/Equalizer/EQ.m +++ b/Classes/04_DSP/Equalizer/EQ.m @@ -104,10 +104,10 @@ classdef EQ end - function [signalclass_out,error_log] = process(obj,signalclass_in, reference_signalclass_in) + function [signalclass_out,noi] = process(obj,signalclass_in, reference_signalclass_in) % actual processing of the signal (steps 1. - 3.) - [signalclass_in.signal,error_log] = obj.process_(signalclass_in.signal', reference_signalclass_in.signal'); + [signalclass_in.signal] = obj.process_(signalclass_in.signal', reference_signalclass_in.signal'); signalclass_in.signal = signalclass_in.signal'; @@ -120,6 +120,9 @@ classdef EQ % write to output signalclass_out = signalclass_in; + noi = signalclass_out; + noi = signalclass_out - reference_signalclass_in; + end function [yout,error_log] = process_(obj,data_in,ref_in) diff --git a/Classes/04_DSP/Equalizer/FFE.m b/Classes/04_DSP/Equalizer/FFE.m index fdd82e5..2f414ec 100644 --- a/Classes/04_DSP/Equalizer/FFE.m +++ b/Classes/04_DSP/Equalizer/FFE.m @@ -101,29 +101,6 @@ classdef FFE < handle end x = [zeros(floor(obj.order/2),1); x; zeros(obj.order,1)]; - - if showviz - f = figure(111); - subplot(2,2,1:2); - hold on - a = scatter(1:numel(x),x,1,'.'); - a2 = scatter(1,1,1,'.'); - a3 = scatter(1,1,2,'.'); - a4 = xline(1); - ylim([-3 3]) - xlim([0 length(x)]); - % subplot(2,2,3) - % dplot = x(1:1+500); - % b = scatter(1:numel(dplot),dplot,5,'x'); - % xline(1) - % xline(obj.order) - % ylim([-3 3]) - % xlim([0 500]); - subplot(2,2,3:4) - c = stem(obj.e); - ylim([-1 1]) - drawnow - end for epoch = 1 : epochs symbol = 0; @@ -153,16 +130,6 @@ classdef FFE < handle obj.e = obj.e - err(symbol) * U / normalizationfactor; % Weight update rule of NLMS end - if mod(sample,100) == 1 && showviz - a2.XData = 1:2*numel(y); - a2.YData = repelem(y, 2); - a3.XData = 1:2*numel(d_hat); - a3.YData = repelem(d_hat, 2); - a4.Value = sample; - % b.YData = x(symbol:symbol+500); - c.YData = obj.e; - drawnow; - end obj.error(epoch,symbol) = err(symbol) * err(symbol)'; % Instantaneous square error diff --git a/Classes/04_DSP/Equalizer/FFE_DCremoval.m b/Classes/04_DSP/Equalizer/FFE_DCremoval.m index 30f689f..dafe65c 100644 --- a/Classes/04_DSP/Equalizer/FFE_DCremoval.m +++ b/Classes/04_DSP/Equalizer/FFE_DCremoval.m @@ -59,7 +59,7 @@ classdef FFE_DCremoval < handle end - function [X] = process(obj, X, D) + function [X,Noi] = process(obj, X, D) % actual processing of the signal (steps 1. - 3.) % 1 normalize RMS @@ -86,6 +86,8 @@ classdef FFE_DCremoval < handle lbdesc = [num2str(obj.order),' tap FFE']; X = X.logbookentry(lbdesc); % append to logbook + Noi = X; + Noi = X - D; end @@ -145,11 +147,21 @@ classdef FFE_DCremoval < handle e_dc_buffer(1) = e_dc_est - obj.mu_dc * err(symbol); e_dc_buffer = circshift(e_dc_buffer,1); end + + e_dc_save(symbol) = e_dc_est; obj.error(epoch,symbol) = err(symbol) * err(symbol)'; % Instantaneous square error end end + + if ~training + figure(1122) + hold on + scatter(1:numel(e_dc_save),e_dc_save,1,'.'); + scatter(1:numel(y),y,1,'.'); + + end end diff --git a/Classes/04_DSP/Equalizer/FFE_FFDCAVG.m b/Classes/04_DSP/Equalizer/FFE_FFDCAVG.m index c0da9a6..93ac1fd 100644 --- a/Classes/04_DSP/Equalizer/FFE_FFDCAVG.m +++ b/Classes/04_DSP/Equalizer/FFE_FFDCAVG.m @@ -101,29 +101,6 @@ classdef FFE_FFDCAVG < handle end x = [zeros(floor(obj.order/2),1); x; zeros(obj.order,1)]; - - if showviz - f = figure(111); - subplot(2,2,1:2); - hold on - a = scatter(1:numel(x),x,1,'.'); - a2 = scatter(1,1,1,'.'); - a3 = scatter(1,1,2,'.'); - a4 = xline(1); - ylim([-3 3]) - xlim([0 length(x)]); - % subplot(2,2,3) - % dplot = x(1:1+500); - % b = scatter(1:numel(dplot),dplot,5,'x'); - % xline(1) - % xline(obj.order) - % ylim([-3 3]) - % xlim([0 500]); - subplot(2,2,3:4) - c = stem(obj.e); - ylim([-1 1]) - drawnow - end for epoch = 1 : epochs @@ -185,19 +162,7 @@ classdef FFE_FFDCAVG < handle normalizationfactor = (U.' * U); obj.e = obj.e - err(symbol) * U / normalizationfactor; % Weight update rule of NLMS end - - if mod(sample,100) == 1 && showviz - a2.XData = 1:2*numel(y); - a2.YData = repelem(y, 2); - a3.XData = 1:2*numel(d_hat); - a3.YData = repelem(d_hat, 2); - - a4.Value = sample; - % b.YData = x(symbol:symbol+500); - c.YData = obj.e; - drawnow; - end obj.error(epoch,symbol) = err(symbol) * err(symbol)'; % Instantaneous square error diff --git a/Classes/04_DSP/Equalizer/FFE_Kalman.m b/Classes/04_DSP/Equalizer/FFE_Kalman.m new file mode 100644 index 0000000..62e7b33 --- /dev/null +++ b/Classes/04_DSP/Equalizer/FFE_Kalman.m @@ -0,0 +1,245 @@ +classdef FFE_Kalman < handle + % Implementation of plain and simple FFE. + % 1) Training mode (stable performance when you use NLMS) + % 2) Decision directed mode + + % Eq = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",0); + + properties + sps % usually 2 + order + e + error + + len_tr + mu_tr + epochs_tr + + mu_dd + epochs_dd + + constellation + + decide + end + + methods + function obj = FFE_Kalman(options) + arguments(Input) + + options.sps = 2; + options.order = 15; + + options.len_tr = 4096; + options.mu_tr = 0; + options.epochs_tr = 5; + + options.mu_dd = 1e-5; + options.epochs_dd = 5; + + options.decide = false; + + end + + fn = fieldnames(options); + for n = 1:numel(fn) + obj.(fn{n}) = options.(fn{n}); + end + + obj.e = zeros(obj.order,1); + obj.error = 0; + + end + + function [X,Noi] = process(obj, X, D) + + % actual processing of the signal (steps 1. - 3.) + % 1 normalize RMS + X = X.normalize("mode","rms"); + + obj.constellation = unique(D.signal); + + % Training Mode + training = 1; + showviz = 0; + obj.equalize(X.signal, D.signal,obj.mu_tr,obj.epochs_tr,obj.len_tr,training,showviz); + + % Decision Directed Mode + n = X.length; + training = 0; + showviz = 0; + [signal,decision]=obj.equalize(X.signal, D.signal,obj.mu_dd,obj.epochs_dd,n,training,showviz); + + % Output Signal + if obj.decide + X.signal = decision; + else + X.signal = signal; + end + + X.fs = D.fs; %change sampling frequency of outgoing signal from fdac e.g. 2 sps to symbol spaced = fsym + lbdesc = [num2str(obj.order),' tap FFE']; + X = X.logbookentry(lbdesc); % append to logbook + + Noi = X; + Noi = X - D; + + + end + + function [y,d_hat] = equalize(obj,x,d,mio,epochs,N,training,showviz) + + arguments + obj + x + d + mio + epochs + N + training + showviz + end + + x = [zeros(floor(obj.order/2),1); x; zeros(obj.order,1)]; + + for epoch = 1 : epochs + + + + symbol = 0; + % Initialization of Kalman filter variables + A = 1; % State transition matrix + H = 1; % Observation matrix + Q = 1e-4; % Process noise covariance + R = 1e-1; % Measurement noise covariance + P = 1; % Initial error covariance + mpi_est = 0; % Initial estimate for MPI noise + K = 0; % Kalman gain + subtract_mpi_est = 1; + + for sample = 1 : obj.sps : N + symbol = symbol + 1; + + % Get the current input sample and the equalizer output + U = x(obj.order + sample - 1 : -1 : sample); + + if subtract_mpi_est || ~training + y(symbol,1) = (obj.e.' * U) - mpi_est .* 1 ; % Subtract MPI estimate + else + y(symbol,1) = obj.e.' * U; + end + + % Decision and error calculation + if training + d_hat(symbol,1) = d(symbol); + else + [~, symbol_idx] = min(abs(y(symbol) - obj.constellation)); % Closest constellation point + d_hat(symbol,1) = obj.constellation(symbol_idx); + end + + err(symbol) = y(symbol) - d_hat(symbol); % Instantaneous residual error + + true_err(symbol) = y(symbol) - d(symbol); + + % Kalman filter update to track the MPI noise + % Prediction step + P = A * P * A' + Q; % Update error covariance + K = P * H' / (H * P * H' + R); % Kalman gain + + % Update step + mpi_est_new = mpi_est + K * (err(symbol) - H * mpi_est); % MPI noise estimation + alpha=0; + mpi_est = alpha * mpi_est + (1 - alpha) * mpi_est_new; + + P = (1 - K * H) * P; % Update error covariance + + % Subtract MPI noise from the signal + if subtract_mpi_est || ~training + y(symbol) = y(symbol);% - mpi_est; + else + + end + + % Equalizer weight update (LMS or NLMS) + if mio ~= 0 + obj.e = obj.e - (mio * err(symbol) * U); % LMS weight update + else + normalizationfactor = (U.' * U); + obj.e = obj.e - err(symbol) * U / normalizationfactor; % NLMS weight update + end + + % Store MPI estimate for visualization + mpi_estimates(symbol) = mpi_est; + + Kgain(symbol) = K; + P_(symbol) = P; + H_(symbol) = H; + + + end %symbols + end %epoch + + if ~training + + if 1 + % figure; + % subplot(2,2,1) + % hold on + % scatter(1:numel(y),y,1,'.'); + % plot(1:numel(mpi_estimates), mpi_estimates); + % subplot(2,2,3) + % plot(1:numel(true_err), true_err); + % subplot(2,2,2) + % scatter(1:numel(y),y'-mpi_estimates,1,'.'); + % xlabel('Sample Index'); + % ylabel('MPI Noise Estimate'); + % title('MPI Noise Estimation Over Time (Kalman Filter)'); + % grid on; + + figure(111); + + subplot(3,1,1) + hold on + cla + scatter(1:numel(y),y,1,'.'); + plot(1:numel(mpi_estimates), mpi_estimates,'DisplayName','mpi est'); + ylim([-2,2]); + + subplot(3,1,2) + cla + plot(1:numel(true_err), true_err,'DisplayName','true err'); + ylim([-1,1]); + + subplot(3,1,3) + plot(1:numel(true_err), true_err-mpi_estimates,'DisplayName',['diff rms: ',num2str(rms(true_err-mpi_estimates))]); + ylim([-1,1]); + legend + + figure(333) + hold on + von = 5000; + bis = 15000; + plot(von:bis,true_err(von:bis),'DisplayName','Error') + plot(von:bis,movmean(true_err(von:bis), [30,30]),'DisplayName','Error') + plot(von:bis,mpi_estimates(von:bis),'DisplayName','Est','LineWidth',2); + legend + + figure(222) + hold on + crr = xcorr(mpi_estimates,true_err,'normalized'); + plot(crr); + + end + + % y = y-mpi_estimates'; + + + end + + + + end + + end +end + diff --git a/Classes/04_DSP/Equalizer/FFE_Kalman_Feedback.m b/Classes/04_DSP/Equalizer/FFE_Kalman_Feedback.m new file mode 100644 index 0000000..68875aa --- /dev/null +++ b/Classes/04_DSP/Equalizer/FFE_Kalman_Feedback.m @@ -0,0 +1,197 @@ +classdef FFE_Kalman_Feedback < handle + % Implementation of plain and simple FFE. + % 1) Training mode (stable performance when you use NLMS) + % 2) Decision directed mode + + % Eq = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",0); + + properties + sps % usually 2 + order + e + error + + len_tr + mu_tr + epochs_tr + + mu_dd + epochs_dd + + constellation + + decide + end + + methods + function obj = FFE_Kalman_Feedback(options) + arguments(Input) + + options.sps = 2; + options.order = 15; + + options.len_tr = 4096; + options.mu_tr = 0; + options.epochs_tr = 5; + + options.mu_dd = 1e-5; + options.epochs_dd = 5; + + options.decide = false; + + end + + fn = fieldnames(options); + for n = 1:numel(fn) + obj.(fn{n}) = options.(fn{n}); + end + + obj.e = zeros(obj.order,1); + obj.error = 0; + + end + + function [X,Noi] = process(obj, X, D) + + % actual processing of the signal (steps 1. - 3.) + % 1 normalize RMS + X = X.normalize("mode","rms"); + + obj.constellation = unique(D.signal); + + % Training Mode + training = 1; + showviz = 0; + obj.equalize(X.signal, D.signal,obj.mu_tr,obj.epochs_tr,obj.len_tr,training); + + % Decision Directed Mode + n = X.length; + training = 0; + showviz = 0; + [signal,decision]=obj.equalize(X.signal, D.signal,obj.mu_dd,obj.epochs_dd,n,training); + + % Output Signal + if obj.decide + X.signal = decision; + else + X.signal = signal; + end + + X.fs = D.fs; %change sampling frequency of outgoing signal from fdac e.g. 2 sps to symbol spaced = fsym + lbdesc = [num2str(obj.order),' tap FFE']; + X = X.logbookentry(lbdesc); % append to logbook + + Noi = X; + Noi = X - D; + + + end + + + + + function [y, d_hat] = equalize(obj, x, d, mio, epochs, N, training) + + arguments + obj + x + d + mio + epochs + N + training + end + + % Initialize Kalman filter variables + A = 0.7; % State transition matrix + H = 1; % Observation matrix + Q = 1e-2; % Process noise covariance + R = 1e-3; % Measurement noise covariance + P = 1; % Initial error covariance + mpi_est = 0; % Initial estimate for MPI noise + K = 0; % Kalman gain + + % Error buffers for parallelized structure + parallel_depth = 100; % Depending on your parallelization depth + err_buffer = zeros(parallel_depth, 1); % Store collected errors + mpi_est_buffer = zeros(parallel_depth, 1); % Store MPI estimates + + x = [zeros(floor(obj.order/2), 1); x; zeros(obj.order, 1)]; + + for epoch = 1 : epochs + symbol = 0; + + for sample = 1 : obj.sps : N + symbol = symbol + 1; + + U = x(obj.order + sample - 1 : -1 : sample); % Input window for FFE + + % Subtract the estimated MPI noise from the input signal before FFE + x_mpi_reduced = U; + y(symbol, 1) = obj.e.' * x_mpi_reduced; % Output of FFE + + y(symbol, 1) = y(symbol, 1) - mpi_est; + + % Decision and error calculation + if training + [~, symbol_idx] = min(abs(d(symbol) - obj.constellation)); % Closest constellation point + d_hat(symbol, 1) = d(symbol); + else + [~, symbol_idx] = min(abs(y(symbol) - obj.constellation)); % Closest constellation point + d_hat(symbol, 1) = obj.constellation(symbol_idx); + end + + % Calculate the residual error + err(symbol) = y(symbol) - d_hat(symbol); % Instantaneous residual error + + % Collect error for feedback after 'parallel_depth' samples + err_buffer(mod(symbol, parallel_depth) + 1) = err(symbol); + + % Update Kalman filter based on the accumulated errors every parallel_depth samples + if mod(symbol, parallel_depth) == 0 + % Compute the mean error over the last 'parallel_depth' samples + avg_error = mean(err_buffer); + + % Prediction step (Kalman filter) + P = A * P * A' + Q; % Update the error covariance + K = P * H' / (H * P * H' + R); % Compute Kalman gain + + % Update MPI noise estimate + mpi_est = mpi_est + K * (avg_error - H * mpi_est); % Update based on averaged error + P = (1 - K * H) * P; % Update error covariance + + end + + % Store MPI estimate for visualization or debugging + k_buffer(symbol) = K; + mpi_est_buffer(symbol) = mpi_est; + + % FFE weight update using NLMS + if mio ~= 0 + obj.e = obj.e - (mio * err(symbol) * U); % LMS weight update + else + normalizationfactor = (U.' * U); + obj.e = obj.e - err(symbol) * U / normalizationfactor; % NLMS weight update + end + + end + end + + if ~training + + figure(111); + + subplot(3,1,1) + hold on + cla + scatter(1:numel(y),y,1,'.'); + plot(1:numel(mpi_est_buffer), mpi_est_buffer,'DisplayName','mpi est'); + + end + end + + end + +end + + diff --git a/Classes/04_DSP/Equalizer/FFE_adaptive_decision.m b/Classes/04_DSP/Equalizer/FFE_adaptive_decision.m index f1f9c9d..14de94c 100644 --- a/Classes/04_DSP/Equalizer/FFE_adaptive_decision.m +++ b/Classes/04_DSP/Equalizer/FFE_adaptive_decision.m @@ -55,7 +55,7 @@ classdef FFE_adaptive_decision < handle end - function [X] = process(obj, X, D) + function [X,Noi] = process(obj, X, D) % actual processing of the signal (steps 1. - 3.) % 1 normalize RMS @@ -80,11 +80,14 @@ classdef FFE_adaptive_decision < handle else X.signal = signal; end + X.fs = D.fs; %change sampling frequency of outgoing signal from fdac e.g. 2 sps to symbol spaced = fsym lbdesc = [num2str(obj.order),' tap FFE']; X = X.logbookentry(lbdesc); % append to logbook - + Noi = X; + Noi = X - D; + end function [y,d_hat] = equalize(obj,x,d,mio,epochs,N,training,showviz) @@ -124,19 +127,21 @@ classdef FFE_adaptive_decision < handle d_hat(symbol,1) = obj.constellation(symbol_idx); end - y_buffer(symbol_idx,1) = y(symbol); y_buffer(symbol_idx,:) = circshift(y_buffer(symbol_idx,:),1); + adap_constellation = mean(y_buffer,2,"omitnan"); + delta_y_y(symbol) = y(symbol) - adap_constellation(symbol_idx); + err(symbol) = y(symbol) - d_hat(symbol); % Instantaneous error + if mio ~= 0 obj.e = obj.e - (mio * err(symbol) * U) ; % Weight update rule of LMS else normalizationfactor = (U.' * U); obj.e = obj.e - err(symbol) * U / normalizationfactor; % Weight update rule of NLMS end - obj.error(epoch,symbol) = err(symbol) * err(symbol)'; % Instantaneous square error diff --git a/projects/400G_FTN_setups/imdd_dsp_approaches.m b/projects/400G_FTN_setups/imdd_dsp_approaches.m index 8192c39..53430bb 100644 --- a/projects/400G_FTN_setups/imdd_dsp_approaches.m +++ b/projects/400G_FTN_setups/imdd_dsp_approaches.m @@ -2,16 +2,16 @@ %% Parameter to simulate and save params = struct; -params.M = [6]; -params.datarate = [448]; -params.rop = [0]; +params.M = [4]; +params.datarate = [224]; +params.rop = [-12:-5]; precomp_mode = 0; %0=do nothing ; 1= measure; 2=precomp active -postfilter = 1; % noise whiten. approach -> Postfilter + MLSE +postfilter = 0; % noise whiten. approach -> Postfilter + MLSE -db_precode = 0; +db_precode = 1; db_encode = 0; -db_channelapproach = 0; +db_channelapproach = 1; if ismac @@ -125,7 +125,7 @@ for M = wh.parameter.M.values Opt_sig = Fiber("fsimu",Opt_sig.fs,"fiber_length",link_length/1000,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.07).process(Opt_sig); % Receiver ROP curve - for i = 1:i_ + parfor i = 1:i_ rop=wh.parameter.rop.values(i); % Set ROP diff --git a/projects/400G_FTN_setups/imdd_dsp_lab_copy_200G.m b/projects/400G_FTN_setups/imdd_dsp_lab_copy_200G.m new file mode 100644 index 0000000..998d9dd --- /dev/null +++ b/projects/400G_FTN_setups/imdd_dsp_lab_copy_200G.m @@ -0,0 +1,247 @@ + +%% Parameter to simulate and save +params = struct; + +params.M = [4]; +params.datarate = [184]; +params.rop = [-12:-5]; + +precomp_mode = 0; %0=do nothing ; 1= measure; 2=precomp active +postfilter = 0; % noise whiten. approach -> Postfilter + MLSE + +db_precode = 0; +db_encode = 0; +db_channelapproach = 0; + +if ismac + precomp_path = "/Users/silasoettinghaus/Documents/MATLAB/imdd_simulation/projects/standard_system"; +else + precomp_path = "C:\Users\sioe\Documents\MATLAB\imdd_simulation\projects\standard_system\"; +end + +precomp_fn = "400G_simulative_setup"; +usemrds = 0; + +name = ['wh_',strrep(num2str(now),'.','')]; + +wh = DataStorage(params); + +wh.addStorage("ber_ffe"); + +%% Init Params +link_length = 10000; %meter +pn_key = 2; +laser_linewidth = 0; + +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 + + % SETUP HERE: %% + kover = 16; + Awg = M8196A("kover",kover); + fdac = Awg.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",Awg.fdac,"applyclipping",1,"clipfactor",1.5,... + "applypulseform",0,"pulseformer",Pform,"randkey",pn_key,... + "db_precode",db_precode,"db_encode",db_encode,... + "mrds_code",usemrds,"mrds_blocklength",512).process(); + + % Digi_sig.eye(fsym,M); + Digi_sig.spectrum("fignum",123434,"displayname",'Digital Tx Signal'); + + if precomp_mode == 1 + freqresp = ChannelFreqResp("Nacq",1024,"Navg",64,"Ncp",63,'f_ref',Digi_sig.fs); + Digi_sig = freqresp.buildOFDM(); + elseif precomp_mode == 2 + Digi_sig = ChannelFreqResp("Nacq",1024,"Navg",64,"Ncp",63,'f_ref',Digi_sig.fs).precomp(Digi_sig,'maxampdb',1,'loadPath',precomp_path,'fileName',precomp_fn); + Digi_sig.spectrum("fignum",11,"displayname",'after precomp'); + end + + %%%%% AWG %%%%%% + El_sig = Awg.process(Digi_sig); + + % El_sig.spectrum("displayname",'el','fignum',123434); + + % El_sig.signal = awgn(El_sig.signal,-3,'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",3).process(El_sig); + % El_sig = El_sig.setPower(6,"dBm"); + + 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.5; + 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.eye(fsym,7); + % + % figure(10) + % hold on + % scatter(El_sig.signal(1:100000)+vbias,(abs(Opt_sig.signal(1:100000)).^2)*1e3,0.1,'.','DisplayName','Modulator TF') + % ylim([0 4]); + % xlim([-u_pi/2, u_pi/2]+vbias); + % xlabel('Input in V') + % ylabel('abs(Output) in mW') + + Optfilter = Filter('filtdegree',6,"f_cutoff",fsym.*0.7,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true); + Opt_sig = Optfilter.process(Opt_sig); + + % Opt_sig.spectrum("fignum",122,"displayname",['Tx SPectrum; PAM ',num2str(M)]); + + Opt_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",0).process(Opt_sig); + + i_ = wh.parameter.rop.length; + + ber_ffe=zeros(i_); + + patten=zeros(i_); + + %%%%% Interference Signal Fiber Prop %%%%%% + + Opt_sig = Fiber("fsimu",Opt_sig.fs,"fiber_length",link_length/1000,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.07).process(Opt_sig); + + % Receiver ROP curve + parfor 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); + patten(i) = Rx_sig.power; + + %%%%%% 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); + + %%%%%% Scope %%%%%% + fadc = 160e9; + Lp_scpe = Filter('filtdegree',4,"f_cutoff",63e9,"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",5.5,"quantbuffer",0.1,'block_dc',1,'lpf_active',1,'H_lpf',Lp_scpe).process(Rx_sig); + + if precomp_mode == 1 + freqresp.estimate(Scpe_sig,"save",true,"savePath",precomp_path,"fileName",precomp_fn); + freqresp.plot(); + end + + Scpe_sig.spectrum("displayname",'After Scope','fignum',123434); + + %%%%%% Sample to 2x fsym %%%%%% + Scpe_sig = Scpe_sig.resample("fs_in",fadc,"fs_out",2*fsym); + + %%%%%% Sync Rx signal with reference %%%%%% + [Scpe_sig,S] = 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",0); + %Eq = VNLE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",[0.0004 0.0005 0.0006],"mu_tr",0,"order",[50,7,7],"sps",2,"decide",1); + + if db_channelapproach + % ref symbols and transm. sequence are precoded + [EQ_sig, Noi] = Eq.process(Scpe_sig,Duobinary().encode(Symbols)); + else + [EQ_sig, Noi] = Eq.process(Scpe_sig,Symbols); + end + + if db_encode || db_channelapproach + EQ_sig = MLSE("DIR",[1,1],"duobinary_output",1,"M",M,"trellis_states",PAMmapper(M,0).levels).process(EQ_sig); + EQ_sig = Duobinary().decode(EQ_sig); + end + + if postfilter + % Noi.spectrum("displayname",'Noise Spectrum','fignum',1234); + % EQ_sig.spectrum("displayname","Signal Spectrum","fignum",1234); + + nc = 2; + burg_coeff = arburg(Noi.signal,nc); + + EQ_sig = EQ_sig.filter(burg_coeff,1); + + % EQ_sig.spectrum("displayname","Signal Spectrum after Postfilter","fignum",1234); + tic + EQ_sig = MLSE("DIR",burg_coeff,"duobinary_output",0,"M",M,"trellis_states",PAMmapper(M,0).levels).process(EQ_sig); + toc + % EQ_sig.spectrum("displayname","Signal Spectrum after MLSE","fignum",1234); + + if 1 + Noi.spectrum('displayname','Noise PSD','fignum',123) + [h,w] = freqz(1,burg_coeff,length(Noi),"whole",Noi.fs); + h = h/max(abs(h)); + hold on + w_ = (w - Noi.fs/2); + plot(w_.*1e-9,20*log10(fftshift(h)),'DisplayName',['', num2str(nc), ' coefficients for burg alg.']); + end + end + + Rx_bits = PAMmapper(M,0).demap(EQ_sig); + [~,errors_bm,ber_ffe(i),errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + disp(['BER: ',sprintf('%.1E',ber_ffe(i)),' - - ROP: ',num2str(patten(i)),'dBm - - PAM-',num2str(M),' - - ',num2str(fsym*1e-9),' GBd']); + + end + + for i = 1:i_ + rop=wh.parameter.rop.values(i); + + wh.addValueToStorage(ber_ffe(i),'ber_ffe',M,datarate,rop); + + end + + toc + + % wh.save('C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\MPI_August\auswertung\') + + end +end + + +disp('Simulation Done!') + + +cols = linspecer(8); + +%cnt = cnt+1; +ber_ffe = wh.getStoValue('ber_ffe',M,datarate,wh.parameter.rop.values); + +% Create the initial plot + +figure(44); +a = gca; +hold on; % Retain the plot so new points can be added without complete redraw + +plot(wh.parameter.rop.values,ber_ffe,"LineWidth",0.5,"LineStyle","-","Marker",".","MarkerSize",15,"DisplayName","FFE only"); +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/400G_FTN_setups/imdd_mpi_dsp.m b/projects/400G_FTN_setups/imdd_mpi_dsp.m index 6825996..a812568 100644 --- a/projects/400G_FTN_setups/imdd_mpi_dsp.m +++ b/projects/400G_FTN_setups/imdd_mpi_dsp.m @@ -2,21 +2,21 @@ %% Parameter to simulate and save params = struct; -params.M = [6]; -params.datarate = [448]; +params.M = [4]; +params.datarate = [300]; params.rop = [0]; precomp_mode = 0; %0=do nothing ; 1= measure; 2=precomp active -postfilter = 1; % noise whiten. approach -> Postfilter + MLSE +postfilter = 0; % noise whiten. approach -> Postfilter + MLSE db_precode = 0; db_encode = 0; db_channelapproach = 0; -laser_linewidth = 1e6; +laser_linewidth = 5e6; random_key_sequence = 2; -random_key_laser_phase = 4; -sir = 25; +random_key_laser_phase = 11; +sir = 20; if ismac precomp_path = "/Users/silasoettinghaus/Documents/MATLAB/imdd_simulation/projects/standard_system"; @@ -25,6 +25,7 @@ else end precomp_fn = "400G_simulative_setup"; + usemrds = 0; name = ['wh_',strrep(num2str(now),'.','')]; @@ -56,7 +57,7 @@ for M = wh.parameter.M.values % MAIN SIGNAL %%%%% Symbol Generation MAIN %%%%%% - [Digi_sig,Symbols,Bits] = PAMsource("fsym",fsym,"M",M,"order",18,"useprbs",0,... + [Digi_sig,Symbols,Bits] = PAMsource("fsym",fsym,"M",M,"order",18,"useprbs",1,... "fs_out",M8199.fdac,"applyclipping",1,"clipfactor",1.5,... "applypulseform",0,"pulseformer",Pform,"randkey",random_key_sequence,... "db_precode",db_precode,"db_encode",db_encode,... @@ -80,7 +81,7 @@ for M = wh.parameter.M.values Digi_sig_I = ChannelFreqResp("Nacq",1024,"Navg",64,"Ncp",63,'f_ref',Digi_sig_I.fs).precomp(Digi_sig_I,'maxampdb',1,'loadPath',precomp_path,'fileName',precomp_fn); Digi_sig.spectrum("fignum",11,"displayname",'after precomp'); end - + %%%%% AWG MAIN %%%%%% @@ -122,7 +123,7 @@ for M = wh.parameter.M.values Optfilter = Filter('filtdegree',6,"f_cutoff",fsym.*0.7,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true); 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); - + %%%%% Interference Signal Fiber Prop 2x fiber length %%%%%% Opt_sig_I_prop = Fiber("fsimu",Opt_sig_I.fs,"fiber_length",2*link_length/1000,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.07).process(Opt_sig_I); @@ -171,7 +172,7 @@ for M = wh.parameter.M.values freqresp.plot(); end - Scpe_sig.spectrum("displayname",'After Scope','fignum',123434); + %Scpe_sig.spectrum("displayname",'After Scope','fignum',123434); %%%%%% Sample to 2x fsym %%%%%% Scpe_sig = Scpe_sig.resample("fs_in",fadc,"fs_out",2*fsym); @@ -183,33 +184,43 @@ for M = wh.parameter.M.values Eq = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",0); %Eq = VNLE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",[0.0004 0.0005 0.0006],"mu_tr",0,"order",[50,7,7],"sps",2,"decide",1); + Eq = EQ("Ne",[50,7,7],"Nb",[0,0,0],"training_length",4096*2,"training_loops",5,"dd_loops",5,"K",2,"DCmu",0.05,"DDmu",[0.0004 0.0004 0.0004 0.0004 ],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1); + + Eq = FFE_Kalman("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",0); + + % Eq = FFE_Kalman_Feedback("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",0); + + % 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",80); + + 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",0,"mu_dc",0.05,"dc_buffer_len",100); + if db_channelapproach % ref symbols and transm. sequence are precoded [EQ_sig, Noi] = Eq.process(Scpe_sig,Duobinary().encode(Symbols)); - else - [EQ_sig, Noi] = Eq.process(Scpe_sig,Symbols); - end - - if db_encode || db_channelapproach EQ_sig = MLSE("DIR",[1,1],"duobinary_output",1,"M",M,"trellis_states",PAMmapper(M,0).levels).process(EQ_sig); EQ_sig = Duobinary().decode(EQ_sig); - end - if postfilter - % Noi.spectrum("displayname",'Noise Spectrum','fignum',1234); - % EQ_sig.spectrum("displayname","Signal Spectrum","fignum",1234); - + elseif db_encode + + [EQ_sig, Noi] = Eq.process(Scpe_sig,Symbols); + EQ_sig = MLSE("DIR",[1,1],"duobinary_output",1,"M",M,"trellis_states",PAMmapper(M,0).levels).process(EQ_sig); + EQ_sig = Duobinary().decode(EQ_sig); + + elseif postfilter + + [EQ_sig, Noi] = Eq.process(Scpe_sig,Symbols); + nc = 2; burg_coeff = arburg(Noi.signal,nc); - + EQ_sig = EQ_sig.filter(burg_coeff,1); - + % EQ_sig.spectrum("displayname","Signal Spectrum after Postfilter","fignum",1234); tic EQ_sig = MLSE("DIR",burg_coeff,"duobinary_output",0,"M",M,"trellis_states",PAMmapper(M,0).levels).process(EQ_sig); toc % EQ_sig.spectrum("displayname","Signal Spectrum after MLSE","fignum",1234); - + if 1 Noi.spectrum('displayname','Noise PSD','fignum',123) [h,w] = freqz(1,burg_coeff,length(Noi),"whole",Noi.fs); @@ -218,12 +229,24 @@ for M = wh.parameter.M.values w_ = (w - Noi.fs/2); plot(w_.*1e-9,20*log10(fftshift(h)),'DisplayName',['', num2str(nc), ' coefficients for burg alg.']); end + + + else + + [EQ_sig, Noi] = Eq.process(Scpe_sig,Symbols); + + if 0 + Noi.spectrum('displayname','Noise PSD','fignum',123) + EQ_sig.plot("displayname",'After EQ','fignum',1112); + end end Rx_bits = PAMmapper(M,0).demap(EQ_sig); [~,errors_bm,ber_ffe(i),errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); disp(['BER: ',sprintf('%.1E',ber_ffe(i)),' - - ROP: ',num2str(patten(i)),'dBm - - PAM-',num2str(M),' - - ',num2str(fsym*1e-9),' GBd']); + + end for i = 1:i_ @@ -255,7 +278,9 @@ figure(44); a = gca; hold on; % Retain the plot so new points can be added without complete redraw -plot(wh.parameter.rop.values,ber_ffe,"LineWidth",0.5,"LineStyle","-","Marker",".","MarkerSize",15,"DisplayName","FFE only"); +dispname = ['Linewidth: ',num2str(laser_linewidth.*1e-6),' MHz']; + +plot(wh.parameter.rop.values,ber_ffe,"LineWidth",0.5,"LineStyle","-","Marker",".","MarkerSize",15,"DisplayName",dispname); yline(3.8e-3,'DisplayName','HD-FEC','LineStyle','--','HandleVisibility','off'); xlabel('Received Optical Power (dBm)'); ylabel('Bit Error Rate (BER)'); diff --git a/projects/MPI_August/mpi_uncorrelated.m b/projects/MPI_August/mpi_uncorrelated.m index ed50896..64e8311 100644 --- a/projects/MPI_August/mpi_uncorrelated.m +++ b/projects/MPI_August/mpi_uncorrelated.m @@ -5,14 +5,15 @@ params = struct; params.M = [4]; params.datarate = [224]; params.sir = [35]; %decibel = attenuation of interference path -params.laser_linewidth = [10e6]; +params.laser_linewidth = [1e6]; params.pn_key = [1]; params.rop = [-12:1:-1]; - +params.rop = 0; usemrds = 0; +wl = 512; name = ['wh_',strrep(num2str(now),'.','')]; @@ -144,7 +145,8 @@ for M = wh.parameter.M.values % Receiver ROP curve - parfor i = 1:i_ + for i = 1:i_ + rop=wh.parameter.rop.values(i); % Set ROP @@ -234,33 +236,43 @@ for M = wh.parameter.M.values % 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); + [Scpe_sig] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym); %%%%% EQUALIZE %%%%%% 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 = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",0); % 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); + [EQ_sig,Noi] = 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_adaptive_decision("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",0,"buffer_length",112); + [EQ_sig,Noi] = 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']); + + Noi.spectrum('displayname','Noise PSD','fignum',123) + [h,w] = freqz(1,burg_coeff,length(Noi),"whole",Noi.fs); + h = h/max(abs(h)); + hold on + w_ = (w - Noi.fs/2); + plot(w_.*1e-9,20*log10(fftshift(h)),'DisplayName',['', num2str(nc), ' coefficients for burg alg.']); + % % % 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);