Minor changes from Star PC.

Better plot in channel freq response

imdd_mpi simulation is a good 400G model

DSP offline analysis script

Classes/00_signals/Signal.m

@@ -337,15 +337,24 @@ classdef Signal
                 obj
                 options.fignum
                 options.displayname = "";
+                options.color = [];
+                options.normalizeToNyquist = 0;
+                options.normalizeTo0dB = 0;
             end
 
             % spectrum_plot(obj.signal,options.fsamp,options.figurename,options.displayname);
 
             N = 2^(nextpow2(length(obj.signal))-8);
 
-            [p_lin,w] = pwelch(obj.signal,hanning(N),N/2,N,obj.fs,"centered","power","mean");
+            if options.normalizeToNyquist==0
-            normalize = 0;
+                [p_lin,w] = pwelch(obj.signal,hanning(N),N/2,N,obj.fs,"centered","power","mean");
-            if normalize 
+                w=w.*1e-9;
+            else
+                [p_lin,w] = pwelch(obj.signal,hanning(N),N/2,N,"centered","power","mean");
+            end
+
+            
+            if options.normalizeTo0dB
                 p_lin = p_lin./ max(p_lin);
                 p_dbm = 10*log10(p_lin); %dB to dBm in case of "power"
                 ylab = "normalized to 0 dB";
@@ -357,14 +366,27 @@ classdef Signal
             figure(options.fignum);  % If figure does not exist, create new figure
             ax = gca;
             hold on
-            plot(w.*1e-9,p_dbm,'DisplayName',options.displayname,'LineWidth',1);
+
-            xlabel("Frequency in GHz");
+            if isempty(options.color)
-            %ylabel("Power/frequency (dB/Hz)");
+                plot(w,p_dbm,'DisplayName',options.displayname,'LineWidth',1);
+            else
+                plot(w,p_dbm,'DisplayName',options.displayname,'LineWidth',1,'Color',options.color);
+            end
+
+            if options.normalizeToNyquist==0
+                xlabel("Frequency in GHz");
+                %xlim([-obj.fs/2 obj.fs/2].*1e-9)
+                edgetick = 2^(nextpow2(obj.fs*1e-9));
+                xticks([-edgetick:16:edgetick]);
+                xlim([100*round( min(w.*1e-9)/100,1)-10,100*round( max(w.*1e-9)/100,1)+10])
+            else
+                xlabel("Normalized Frequency");
+                xlim([-pi, pi]);
+            end
+            
+            ylabel("Power/frequency (dB/Hz)");
             ylabel(ylab);
-            xlim([-obj.fs/2 obj.fs/2].*1e-9)
+            
-            edgetick = 2^(nextpow2(obj.fs*1e-9));
-            % xticks([-edgetick:16:edgetick]);
-            xlim([100*round( min(w.*1e-9)/100,1)-10,100*round( max(w.*1e-9)/100,1)+10])
             ylim([min(floor( min(p_dbm))-3 , ax.YLim(1)),  max(ceil( max(p_dbm) )+(3), ax.YLim(2))]);
             yticks([-200:10:10]);
             grid on

Classes/01_transmit/ChannelFreqResp.m

@@ -240,6 +240,13 @@ classdef ChannelFreqResp < handle
             plot(obj.faxis/1e9, unwrap(angle(1./Havg)),"LineWidth",2,"Color",[0.3467    0.5360    0.6907]) ;
             xlim([0.2 .5*max(obj.faxis)*1e-9]); grid on;
 
+            %%% plot for publication
+            figure(98989);hold all;box on;title('Magnitude Freq. Response');
+            xlim([0.2 .5*max(obj.faxis)*1e-9]);
+            ylim([-20, 2]);
+            plot(obj.faxis/1e9, 20*log10(abs(Havg)),'LineWidth',2);
+            grid on;
+
 
         end
 

Classes/01_transmit/PAMmapper.m

@@ -293,7 +293,7 @@ classdef PAMmapper
         end
 
         function [out] = separate_pamlevels(obj,data_in)
-            
+            %data_in is Signal class
             %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
@@ -308,6 +308,18 @@ classdef PAMmapper
 
         end
 
+        function [Signal_out] = quantize(obj,Signal_in)
+            constellation = obj.get_levels();
+            constellation = constellation ./ rms(constellation);
+            Signal_out = Signal_in;
+            dist = abs(Signal_in.signal - constellation);
+            [~,symbol_idx]  = min(dist,[],2); % decision for closest constellation point
+            Signal_out.signal = constellation(symbol_idx);
+
+            Signal_out.signal = reshape(Signal_out.signal,size(Signal_in.signal));
+
+        end
+
     end
 
 end

projects/400G_FTN_setups/imdd_mpi_dsp.m

@@ -3,25 +3,27 @@
 params = struct;
 
 params.M                        = [4];
-params.datarate                 = [300];
+params.datarate                 = [448];
 params.rop                      = [0];
+params.sir                      = 40;%15:1:40;
+params.random_key_laser_phase   = 10:20;
 
 precomp_mode    = 0; %0=do nothing ; 1= measure; 2=precomp active
 postfilter      = 0; % noise whiten. approach -> Postfilter + MLSE
 
-db_precode      = 0;
+db_precode      = 1;
 db_encode       = 0;
-db_channelapproach = 0;
+db_channelapproach = 1;
 
-laser_linewidth = 5e6;
+laser_linewidth = 50e5;
-random_key_sequence = 2;
+random_key_sequence = 15;
-random_key_laser_phase = 11;
+random_key_laser_phase = 66;
 sir = 20;
 
 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\";
+    precomp_path = "C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\standard_system\";
 end
 
 precomp_fn = "400G_simulative_setup";
@@ -32,7 +34,8 @@ name = ['wh_',strrep(num2str(now),'.','')];
 
 wh = DataStorage(params);
 
-wh.addStorage("ber_ffe");
+wh.addStorage("ber_vnle");
+wh.addStorage("ber_mlse");
 
 %% Init Params
 link_length = 1000; %meter
@@ -43,6 +46,7 @@ cnt=0;
 disp(['Start Simulation of ',num2str(endcnt),' loops...'])
 tic
 
+for random_key_laser_phase = wh.parameter.random_key_laser_phase.values
 for M = wh.parameter.M.values
     for datarate = wh.parameter.datarate.values
 
@@ -58,27 +62,28 @@ for M = wh.parameter.M.values
 
         %%%%% Symbol Generation MAIN %%%%%%
         [Digi_sig,Symbols,Bits] = PAMsource("fsym",fsym,"M",M,"order",18,"useprbs",1,...
-            "fs_out",M8199.fdac,"applyclipping",1,"clipfactor",1.5,...
+            "fs_out",M8199.fdac,"applyclipping",0,"clipfactor",1.5,...
             "applypulseform",0,"pulseformer",Pform,"randkey",random_key_sequence,...
             "db_precode",db_precode,"db_encode",db_encode,...
             "mrds_code",usemrds,"mrds_blocklength",512).process();
 
         %%%%% Symbol Generation INTERFERENCE %%%%%%
         [Digi_sig_I,Symbols_I,Bits_I] = PAMsource("fsym",fsym,"M",M,"order",18,"useprbs",0,...
-            "fs_out",M8199.fdac,"applyclipping",1,"clipfactor",1.5,...
+            "fs_out",M8199.fdac,"applyclipping",0,"clipfactor",1.5,...
             "applypulseform",0,"pulseformer",Pform,"randkey",random_key_sequence+1,...
             "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');
+        % Digi_sig.normalize("mode","rms").spectrum("displayname",'Tx Signal','fignum',10);
 
         if precomp_mode == 1 %measure
             freqresp = ChannelFreqResp("Nacq",1024,"Navg",64,"Ncp",63,'f_ref',Digi_sig.fs);
             Digi_sig = freqresp.buildOFDM();
+            Digi_sig_I = freqresp.buildOFDM();
         elseif precomp_mode == 2 %apply
-            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 = ChannelFreqResp("Nacq",1024,"Navg",64,"Ncp",63,'f_ref',Digi_sig.fs).precomp(Digi_sig,'maxampdb',3,'loadPath',precomp_path,'fileName',precomp_fn);
-            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_I = ChannelFreqResp("Nacq",1024,"Navg",64,"Ncp",63,'f_ref',Digi_sig_I.fs).precomp(Digi_sig_I,'maxampdb',3,'loadPath',precomp_path,'fileName',precomp_fn);
             Digi_sig.spectrum("fignum",11,"displayname",'after precomp');
         end
 
@@ -88,7 +93,7 @@ for M = wh.parameter.M.values
         El_sig = M8199.process(Digi_sig);
 
         %%%%% Lowpass el. components %%%%%%
-        El_sig = Filter('filtdegree',2,"f_cutoff",60e9,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(El_sig);
+        El_sig = Filter('filtdegree',2,"f_cutoff",100e9,"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);
@@ -96,11 +101,13 @@ for M = wh.parameter.M.values
         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)));
 
+        % El_sig.spectrum("displayname",'Transmit PDS','fignum',10);
+
         %%%%% AWG INTERFERENCE %%%%%%
         El_sig_I = M8199.process(Digi_sig_I);
 
         %%%%% Lowpass el. components %%%%%%
-        El_sig_I = Filter('filtdegree',2,"f_cutoff",60e9,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(El_sig_I);
+        El_sig_I = Filter('filtdegree',3,"f_cutoff",100e9,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(El_sig_I);
 
         %%%%% Electrical Driver Amplifier %%%%%%
         El_sig_I = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",3).process(El_sig_I);
@@ -115,146 +122,237 @@ for M = wh.parameter.M.values
         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",random_key_laser_phase).process(El_sig);
-        Optfilter = Filter('filtdegree',6,"f_cutoff",fsym.*0.7,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true);
+        Optfilter = Filter('filtdegree',3,"f_cutoff",110e9,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true);
         Opt_sig = Optfilter.process(Opt_sig);
         Opt_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",0).process(Opt_sig);
 
         [Opt_sig_I] = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig_I.fs,"lambda",1290,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",random_key_laser_phase+1).process(El_sig_I);
-        Optfilter = Filter('filtdegree',6,"f_cutoff",fsym.*0.7,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true);
+        Optfilter = Filter('filtdegree',3,"f_cutoff",110e9,"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);
+        Opt_sig_I = 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);
-
-        %%%%% Set SIR %%%%%%
-        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 = Opt_sig_I_prop + Opt_sig;
-
-        %%%%% 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);
-
-
 
+        % ber=zeros(i_);
+        % patten=zeros(i_);
         i_ = wh.parameter.rop.length;
+        j_ = wh.parameter.sir.length;
 
-        ber_ffe=zeros(i_);
+        ber_vnle=zeros(i_,j_);
+        ber_mlse=zeros(i_,j_,3);
 
-        patten=zeros(i_);
+        for j = 1:j_
     
-        % Receiver ROP curve
+            sir = wh.parameter.sir.values(j);
-        for i = 1:i_
-            rop=wh.parameter.rop.values(i);
     
-            % Set ROP
+            %%%%% Set SIR %%%%%%
-            Rx_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",rop).process(Opt_sig);
+            Opt_sig_I_atten = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",Opt_sig.power-sir).process(Opt_sig_I);
-            patten(i) = Rx_sig.power;
     
-            %%%%%% Square Law %%%%%%
+            %%%%% ADD Interference and Main Signal  %%%%%%
-            Rx_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20,"nep",1.8e-11).process(Rx_sig);
+            Opt_sig_MPI = Opt_sig_I_atten + Opt_sig;
     
-            %%%%%% Lowpass PhDiode %%%%%%
+            %%%%% Interference Signal Fiber Prop %%%%%%
-            Rx_sig = Filter('filtdegree',2,"f_cutoff",70e9,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true).process(Rx_sig);
+            Opt_sig_MPI    = 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_MPI);
     
-            %%%%%% Scope %%%%%%
-            fadc = 256e9;
-            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",10,"quantbuffer",0.1,'block_dc',1,'lpf_active',1,'H_lpf',Lp_scpe).process(Rx_sig);
     
-            if precomp_mode == 1
+                % Receiver ROP curve
-                freqresp.estimate(Scpe_sig,"save",true,"savePath",precomp_path,"fileName",precomp_fn);
+                for i = 1:i_
-                freqresp.plot();
+                    rop=wh.parameter.rop.values(i);
-            end
         
-            %Scpe_sig.spectrum("displayname",'After Scope','fignum',123434);
+                    % Set ROP
+                    Rx_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",rop).process(Opt_sig_MPI);
+                    % patten(i) = Rx_sig.power;
         
-            %%%%%% Sample to 2x fsym %%%%%%
+                    %%%%%% Square Law %%%%%%
-            Scpe_sig =  Scpe_sig.resample("fs_in",fadc,"fs_out",2*fsym);
+                    Rx_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20,"nep",1.8e-11).process(Rx_sig);
         
-            %%%%%% Sync Rx signal with reference %%%%%%
+                    %%%%%% Lowpass PhDiode %%%%%%
-            [Scpe_sig,S] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym);
+                    Rx_sig = Filter('filtdegree',2,"f_cutoff",70e9,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true).process(Rx_sig);
         
-            %%%%% EQUALIZE %%%%%%
+                    %%%%%% Scope %%%%%%
-            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);
+                    fadc = 256e9;
-            %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);
+                    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",10,"quantbuffer",0.1,'block_dc',1,'lpf_active',1,'H_lpf',Lp_scpe).process(Rx_sig);
         
-             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);
+                    if precomp_mode == 1
+                        freqresp.estimate(Scpe_sig,"save",true,"savePath",precomp_path,"fileName",precomp_fn);
+                        freqresp.plot();
+                    end
 
-            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);
+                    % Scpe_sig_normalized = Scpe_sig.normalize("mode","rms");
         
-           %  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);
+                    % Scpe_sig.normalize("mode","rms").spectrum("displayname",'After Scope','fignum',10);
         
-            % 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);
+                    %%%%%% Sample to 2x fsym %%%%%%
+                    Scpe_sig =  Scpe_sig.resample("fs_in",fadc,"fs_out",2*fsym);
         
-           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);
+                    %%%%%% Sync Rx signal with reference %%%%%%
+                    [Scpe_sig,S] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym);
         
-            if db_channelapproach
+                    %%%%% EQUALIZE %%%%%%
-                % ref symbols and transm. sequence are precoded
+                     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_sig, Noi] = Eq.process(Scpe_sig,Duobinary().encode(Symbols));
+                    % 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_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 db_encode
+                  % Eq     = EQ("Ne",[50,7,7],"Nb",[0,0,0],"training_length",4096*2,"training_loops",5,"dd_loops",5,"K",2,"DCmu",0.0,"DDmu",[0.0004 0.0004 0.0004 0.0004 ],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1);
         
-                [EQ_sig, Noi] = Eq.process(Scpe_sig,Symbols);
+                    % 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_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     = 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);
+                   % 
+                   % Eq = EQ("Ne",[50,7,7],"Nb",[0,0,0],"training_length",4096*2,"training_loops",5,"dd_loops",5,"K",2,"DCmu",0.0,"DDmu",[0.0004 0.0004 0.0004 0.0004 ],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1);
+        
+                    if db_channelapproach
+                        % ref symbols and transm. sequence are precoded
+                        [EQ_sig, Noi] = Eq.process(Scpe_sig,Duobinary().encode(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);
+                        Rx_bits = PAMmapper(M,0).demap(EQ_sig);
+                        [~,~,ber_vnle(i,j),~] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
+        
+                    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);
+        
+                        % EQ_sig.plot("displayname",'After VNLE','fignum',90,'clear',1);
+
+                        % Quantization is too far from orig. symbols ->
+                        % error psd is quite different
+                        % Sym_ = PAMmapper(M,0).quantize(EQ_sig);
+                        % Noi_ = Sym_-EQ_sig;
+                        % Noi_.normalize('mode','rms').spectrum('displayname','Noise PSD','fignum',1234,'normalizeTo0dB',1,'normalizeToNyquist',1,'color',cols(nc+1,:));
+
+                        Rx_bits = PAMmapper(M,0).demap(EQ_sig);
+                        [~,~,ber_vnle(i,j),~] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
+
+                        EQ_sig.normalize('mode','rms').spectrum('displayname','EQ Out','fignum',1234,'normalizeTo0dB',1,'normalizeToNyquist',1,'color',cols(nc,:));
+                        
+                        Noi.normalize('mode','rms').spectrum('displayname','Noise PSD optimal','fignum',1234,'normalizeTo0dB',1,'normalizeToNyquist',1,'color',cols(nc+1,:));
+
+                        for nc = 1:3
+                        
+                            burg_coeff = arburg(Noi.signal,nc);
+            
+                            EQ_sig_filt = EQ_sig.filter(burg_coeff,1);
+            
+                            % EQ_sig.spectrum("displayname","Signal Spectrum after Postfilter","fignum",1234);
+                            
+                            EQ_sig_mlse = MLSE("DIR",burg_coeff,"duobinary_output",0,"M",M,"trellis_states",PAMmapper(M,0).levels).process(EQ_sig_filt);
+                            
+                            % EQ_sig.spectrum("displayname","Signal Spectrum after MLSE","fignum",1234);
+            
+                            if 1
+                                cols = linspecer(12);
+
+
+
+                                EQ_sig_filt.normalize('mode','rms').spectrum('displayname','Noise PSD','fignum',1234,'normalizeTo0dB',1,'normalizeToNyquist',1,'color',cols(nc+2,:));
+
+                                % [h,w] = freqz(1,burg_coeff,length(Noi),"whole",Noi.fs);
+                                % h = h/max(abs(h));
+                                % hold on
+                                % w_ = (w - Noi.fs/2);
+                                % figure(123)
+                                % plot(w_.*1e-9,20*log10(fftshift(h)),'DisplayName',['', num2str(nc), ' coefficients for burg alg.']);
+
+                                [h,w] = freqz(1,burg_coeff,length(Noi),"whole");
+                                h = h/max(abs(h));
+                                hold on
+                                w_ = (w - pi);
+                                plot(w_,20*log10(fftshift(h)),'DisplayName',['', num2str(nc), ' coefficients for burg alg.']);
+
+                                
+                            end
+    
+                            Rx_bits = PAMmapper(M,0).demap(EQ_sig_mlse);
+                            [~,errors_bm,ber_mlse(i,j,nc),errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
+                            % disp(['BER: ',sprintf('%.1E',ber_mlse(i,j)),' - - ROP: ',num2str(patten(i)),'dBm - - PAM-',num2str(M),' - - ',num2str(fsym*1e-9),' GBd']);
+
+                        end
+
+                    else
+
+                                                % S = Scpe_sig.signal;
+                        % N1 = 101;
+                        % 
+                        % % Initialize the running sum with the first window's sum
+                        % running_sum = mean( S(1:N1) );
+                        % 
+                        % % Calculate the first output value
+                        % S_(1) = S(1) - running_sum;
+                        % 
+                        % % Recursive running sum filter
+                        % for n = 2 : length(S) - N1
+                        %     % Update running sum by removing the oldest sample and adding the newest
+                        %     avg_win(n) = mean( S(n:n+N1) );
+                        %     S_(n) = S(n) - avg_win(n);
+                        % end
+                        % 
+                        % % movmean
+                        % S__ = S - movmean(S,[floor(N1/2),ceil(N1/2)]);
+                        % 
+                        % % recursion
+                        % % Initialize the moving sum for the first window
+                        % half_window = (N1 - 1) / 2;
+                        % moving_sum = sum(S(1:N1));
+                        % 
+                        % % Calculate the first element of R1
+                        % S___(half_window+1) = S(half_window+1) - (moving_sum / N1);
+                        % 
+                        % % Loop over the signal and apply the recursive moving average subtraction
+                        % for n = (half_window+2):(length(S)-half_window)
+                        %     % Update the moving sum by subtracting the oldest value and adding the new one
+                        %     moving_sum = moving_sum - S(n-half_window-1) + S(n+half_window);
+                        % 
+                        %     % Calculate the new value of R1
+                        %     S___(n) = S(n) - (moving_sum / N1);
+                        % end
+
+
+                        [EQ_sig, Noi] = Eq.process(Scpe_sig,Symbols);
+
+                        if 0
+                            Noi.spectrum('displayname','Noise PSD','fignum',123,'normalizeTo0dB',1,'normalizeToNyquist',1);
+                            EQ_sig.plot("displayname",'After EQ','fignum',1113);
+                        end
+                        % 
+                        Rx_bits = PAMmapper(M,0).demap(EQ_sig);
+                        [~,errors_bm,ber_vnle(i,j),errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
+
+                    end
 
-                [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);
-                    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
 
-            else
+       for j = 1:j_
+           sir = wh.parameter.sir.values(j);
+            for i = 1:i_
+                rop=wh.parameter.rop.values(i);
     
-                 [EQ_sig, Noi] = Eq.process(Scpe_sig,Symbols);
+                wh.addValueToStorage(ber_vnle(i,j),'ber_vnle',M,datarate,rop,sir,random_key_laser_phase);
+                wh.addValueToStorage(ber_mlse(i,j,:),'ber_mlse',M,datarate,rop,sir,random_key_laser_phase);
     
-                 if 0
-                     Noi.spectrum('displayname','Noise PSD','fignum',123)
-                     EQ_sig.plot("displayname",'After EQ','fignum',1112);
-                 end
             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
 
@@ -262,25 +360,36 @@ for M = wh.parameter.M.values
 
     end
 end
-
+end
 
 disp('Simulation Done!')
     
-
+ber_mlse=[];
+ber_vnle=[];
 cols = linspecer(8);
+random_key_laser_phase_ = wh.parameter.random_key_laser_phase.values;
+cnt = 0;
+for r = random_key_laser_phase_
+    cnt = cnt+1;
+    ber_mlse(cnt,:,1:3)        = wh.getStoValue('ber_mlse',M,datarate,wh.parameter.rop.values(1),wh.parameter.sir.values,r);
+    ber_vnle(cnt,:,1)        = wh.getStoValue('ber_vnle',M,datarate,wh.parameter.rop.values(1),wh.parameter.sir.values,r);
+end
 
-%cnt = cnt+1;
+ber_mlse=squeeze(mean(ber_mlse,1));
-ber_ffe        = wh.getStoValue('ber_ffe',M,datarate,wh.parameter.rop.values);
+ber_vnle = mean(ber_vnle,1);
-
 % Create the initial plot
 
 figure(44);
 a = gca;
 hold on; % Retain the plot so new points can be added without complete redraw
 
-dispname = ['Linewidth: ',num2str(laser_linewidth.*1e-6),' MHz'];
+dispname = ['Lw: ',num2str(laser_linewidth.*1e-6),' MHz'];
+
+plot(wh.parameter.sir.values,ber_vnle,"LineWidth",0.5,"LineStyle","-","Marker",".","MarkerSize",15,"DisplayName",['VNLE ',dispname]);
+plot(wh.parameter.sir.values,ber_mlse(:,1),"LineWidth",0.5,"LineStyle","-","Marker",".","MarkerSize",15,"DisplayName",['MLSE 1 ',dispname]);
+plot(wh.parameter.sir.values,ber_mlse(:,2),"LineWidth",0.5,"LineStyle","-","Marker",".","MarkerSize",15,"DisplayName",['MLSE 2',dispname]);
+plot(wh.parameter.sir.values,ber_mlse(:,3),"LineWidth",0.5,"LineStyle","-","Marker",".","MarkerSize",15,"DisplayName",['MLSE 3',dispname]);
 
-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)');

projects/Lab_2024/offline_dsp_analysis/load_n_dsp.m

@@ -2,76 +2,123 @@
 % load data points
 
 foldername = 'C:\Users\Silas\Nextcloud4\Dokumente\02_Ablage_Office\Lab_Data_24\sir_sweep_pam4';
-filename = 'PAM4_10km';
+filename = 'PAM4_10km_';
 
-stuff = load([foldername,filesep,filename,'_wh']);
+
+stuff = load([foldername,filesep,filename,'wh']);
+%stuff = load([foldername,filesep,'PAM4_v2_10km_wh']);
 wh = stuff.obj;
 i_atten_vals = wh.parameter.i_atten.values;
 v_bias = wh.parameter.vbias.values(1);
 awg_vpp = wh.parameter.awg_vpp.values(1);
-eq_mode = wh.parameter.eq_mode.values(1);
+eqmode = wh.parameter.eq_mode.values(1);
-
 
+M = 4;
 
 % Tx Bits
 % Tx Symbols
 % Rx Signal
 
-Bits = load([foldername, filesep, filename, '_bits'],'Bits');
+Bits = load([foldername, filesep, filename, 'bits'],'Bits');
 Bits = Bits.Bits;
-Symbols = load([foldername, filesep, filename, '_symbols_92gbd'],'Symbols');
+Symbols = load([foldername, filesep, filename, 'symbols'],'Symbols');
 Symbols = Symbols.Symbols;
 
+% all_signals = load([foldername, filesep, filename, 'best_rx_signals']);
+% all_signals = all_signals.all_signals;
+
 Eq = EQ("Ne",[50,7,7],"Nb",[0,0,0],"training_length",4096*2,"training_loops",5,"dd_loops",5,"K",2,"DCmu",0.0,"DDmu",[0.0004 0.0004 0.0004 0.0004 ],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1);
 
+cnt = 0;
+for atten = 40:-4:0
+    cnt = cnt+1;
 
-for atten = 0:4:40
+    Signal_cell = load([foldername, filesep, filename, 'rx_signal_iatten_',num2str(atten)]);
-    Signal_cell = load([foldername, filesep, filename, '_rx_signal_iatten_',num2str(atten)]);
-    Scpe_sig = Signal_cell.S{1};
     
-    bers = wh.getStoValue('ber',v_bias,awg_vpp,eq_mode,atten);
+    % Scpe_sig = all_signals(cnt);
+
+    ber_ffe_only_ = zeros(numel(Signal_cell.S),1);
+    ber_mlse_ = zeros(numel(Signal_cell.S),1);
+    ber_db_ = zeros(numel(Signal_cell.S),1);
+    
+    parfor i = 1:numel(Signal_cell.S)
+        
+        Scpe_sig = Signal_cell.S{i};
+
+        sir_measured = wh.getStoValue('sir',v_bias,awg_vpp,eqmode,atten);
+        
+        if eqmode == 2
+            [EQ_sig] = Eq.process(Scpe_sig,Symbols);
+        
+            EQ_sig.plot("fignum",50,"displayname",'After EQ','clear',1);
+        
+            Noi = EQ_sig-Symbols;
+        
+            Rx_bits = PAMmapper(M,0).demap(EQ_sig);
+            [~,errors_bm,ber_ffe_only_(i),errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
+        
+            nc = 2;
+            burg_coeff = arburg(Noi.signal,nc);
+        
+            EQ_sig = EQ_sig.filter(burg_coeff,1);
+        
+            if 1
+                Noi.spectrum('displayname',['SIR: ', num2str(sir_measured), ' dB '],'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
+        
+            if 0
+                figure(53);
+                clf
+                constellation = unique(Symbols.signal);
+                received = NaN(numel(constellation),length(Symbols));
+                for lvl = 1:numel(constellation)
+                    received(lvl,Symbols.signal==constellation(lvl)) = EQ_sig.signal(Symbols.signal==constellation(lvl));
+                    hold on
+                    histogram(received(lvl,:),1000,"EdgeAlpha",0);
+                end
+            end
+        
+            EQ_sig = MLSE("DIR",burg_coeff,"duobinary_output",0,"M",M,"trellis_states",PAMmapper(M,0).levels).process(EQ_sig);
+        
+            Rx_bits = PAMmapper(M,0).demap(EQ_sig);
+            [~,errors_bm,ber_mlse_(i),errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
+        
+            disp(['FFE: ',sprintf('%.1E',ber_ffe_only_(i)),' -> PF -> MLSE: ',sprintf('%.1E',ber_mlse_(i)),'  dB ']);
+        elseif eqmode == 3
 
     
-    [EQ_sig] = Eq.process(Scpe_sig,Symbols);
+            [EQ_sig, Noi] = Eq.process(Scpe_sig,Duobinary().encode(Symbols));
 
-    EQ_sig.plot("fignum",50,"displayname",'After EQ','clear',1);
+            % EQ_sig.plot("fignum",50,"displayname",'After EQ','clear',1);
 
-    Noi = EQ_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);
 
-    Rx_bits = PAMmapper(M,0).demap(EQ_sig);
+            Rx_bits = PAMmapper(M,0).demap(EQ_sig);
-    [~,errors_bm,ber_ffe_only,errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
+            [~,errors_bm,ber_db_(i),errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
 
-    nc = 2;
+            disp([' DB Precode -> Channel -> FFE -> Decode/ Mod  ',sprintf('%.1E',ber_db_(i)),' ']);
-    burg_coeff = arburg(Noi.signal,nc);
 
-    EQ_sig = EQ_sig.filter(burg_coeff,1);
 
-    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
 
-    if 0
-        figure(53);
-        constellation = unique(Symbols.signal);
-        received = NaN(numel(constellation),length(Symbols));
-        for lvl = 1:numel(constellation)
-            received(lvl,Symbols.signal==constellation(lvl)) = EQ_sig.signal(Symbols.signal==constellation(lvl));
-            hold on
-            histogram(received(lvl,:),1000,"EdgeAlpha",0);
         end
+
     end
 
-    EQ_sig = MLSE("DIR",burg_coeff,"duobinary_output",0,"M",M,"trellis_states",PAMmapper(M,0).levels).process(EQ_sig);
+    if eqmode == 2
+        [ber_ffe_2(cnt),idx] = min(ber_ffe_only_);
+        [ber_mlse_2(cnt),idx] = min(ber_mlse_);
 
-    Rx_bits = PAMmapper(M,0).demap(EQ_sig);
+    elseif eqmode == 3
-    [~,errors_bm,ber,errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
+        [ber_db(cnt),idx] = min(ber_db_);
+    end
     
-    disp(['FFE: ',sprintf('%.1E',ber_ffe_only),' -> PF -> MLSE: ',sprintf('%.1E',ber),'  dB | PD_in: ',num2str(pd_in),' dBm']);
+    all_signals(:,cnt) = Signal_cell.S{idx};
 
 end
 
@@ -83,21 +130,23 @@ cols = linspecer(8);
 i_atten_vals = wh.parameter.i_atten.values;
 v_bias = wh.parameter.vbias.values(1);
 awg_vpp = wh.parameter.awg_vpp.values(1);
-eq_mode = wh.parameter.eq_mode.values(1);
+eqmode = wh.parameter.eq_mode.values(1);
 
-bers = wh.getStoValue('ber',v_bias,awg_vpp,eq_mode,i_atten_vals);
+bers = wh.getStoValue('ber',v_bias,awg_vpp,eqmode,i_atten_vals);
+sirs = wh.getStoValue('sir',v_bias,awg_vpp,eqmode,i_atten_vals);
 
 figure(90);
 hold on; % Retain the plot so new points can be added without complete redraw
 
 % Plot the data and get the line handle
-hLine = plot(i_atten_vals, bers, "LineWidth", 0.5, "LineStyle", "-", "Marker", ".", "MarkerSize", 15, "DisplayName", 'nbla');
+if eqmode == 2
+    hLine1 = plot(sirs, ber_ffe_2, "LineWidth", 0.5, "LineStyle", "-", "Marker", ".", "MarkerSize", 15, "DisplayName", 'VNLE 2');
+    hLine2 = plot(sirs, ber_mlse_2, "LineWidth", 0.5, "LineStyle", "-", "Marker", ".", "MarkerSize", 15, "DisplayName", 'VNLE+MLSE 2');
+ elseif eqmode == 3
+    hLine2 = plot(sirs, ber_db, "LineWidth", 0.5, "LineStyle", "-", "Marker", ".", "MarkerSize", 15, "DisplayName", 'DB ');
+end
 
-% Customize the data tips
+hLine3 = plot(sirs, bers, "LineWidth", 0.5, "LineStyle", "-", "Marker", ".", "MarkerSize", 15, "DisplayName", 'from measurement');
-% Set labels for existing data tip rows
-hLine.DataTipTemplate.DataTipRows(1).Label = 'Fsym';
-hLine.DataTipTemplate.DataTipRows(2).Label = 'BER';
-hLine.DataTipTemplate.DataTipRows(2).Format = '%.2e'; % Format BER as "3e-4"
 
 
 % Continue with the rest of your plot settings