diff --git a/Classes/00_signals/Signal.m b/Classes/00_signals/Signal.m index 8224f4b..cb5c411 100644 --- a/Classes/00_signals/Signal.m +++ b/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"); - normalize = 0; - if normalize + if options.normalizeToNyquist==0 + [p_lin,w] = pwelch(obj.signal,hanning(N),N/2,N,obj.fs,"centered","power","mean"); + 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"); - %ylabel("Power/frequency (dB/Hz)"); + + if isempty(options.color) + 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 diff --git a/Classes/01_transmit/ChannelFreqResp.m b/Classes/01_transmit/ChannelFreqResp.m index fe644a6..e77cf43 100644 --- a/Classes/01_transmit/ChannelFreqResp.m +++ b/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 diff --git a/Classes/01_transmit/PAMmapper.m b/Classes/01_transmit/PAMmapper.m index 64ef02f..c233a13 100644 --- a/Classes/01_transmit/PAMmapper.m +++ b/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 diff --git a/projects/400G_FTN_setups/imdd_mpi_dsp.m b/projects/400G_FTN_setups/imdd_mpi_dsp.m index a812568..d9464c2 100644 --- a/projects/400G_FTN_setups/imdd_mpi_dsp.m +++ b/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; -random_key_sequence = 2; -random_key_laser_phase = 11; +laser_linewidth = 50e5; +random_key_sequence = 15; +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_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 = 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',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); - - %%%%% 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); - - + 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); + % 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_ + + sir = wh.parameter.sir.values(j); + + %%%%% Set SIR %%%%%% + Opt_sig_I_atten = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",Opt_sig.power-sir).process(Opt_sig_I); + + %%%%% ADD Interference and Main Signal %%%%%% + Opt_sig_MPI = Opt_sig_I_atten + Opt_sig; + + %%%%% Interference Signal Fiber Prop %%%%%% + 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); + + + % 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_MPI); + % 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 = 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 + freqresp.estimate(Scpe_sig,"save",true,"savePath",precomp_path,"fileName",precomp_fn); + freqresp.plot(); + end - % Receiver ROP curve - for i = 1:i_ - rop=wh.parameter.rop.values(i); + % Scpe_sig_normalized = Scpe_sig.normalize("mode","rms"); + + % Scpe_sig.normalize("mode","rms").spectrum("displayname",'After Scope','fignum',10); + + %%%%%% 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); + + % 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 = 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); + % + % 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); - % 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; + % 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,:)); - %%%%%% Square Law %%%%%% - Rx_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20,"nep",1.8e-11).process(Rx_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); - %%%%%% Lowpass PhDiode %%%%%% - Rx_sig = Filter('filtdegree',2,"f_cutoff",70e9,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true).process(Rx_sig); + 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,:)); - %%%%%% 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); + 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); - 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); + EQ_sig_filt.normalize('mode','rms').spectrum('displayname','Noise PSD','fignum',1234,'normalizeTo0dB',1,'normalizeToNyquist',1,'color',cols(nc+2,:)); - %%%%%% Sync Rx signal with reference %%%%%% - [Scpe_sig,S] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym); + % [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.']); - %%%%% 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); + [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.']); - 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); + + 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']); - 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); + end - % 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); + else - % 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); + % 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 = 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)); - 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); + [EQ_sig, Noi] = Eq.process(Scpe_sig,Symbols); - elseif db_encode + 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); - [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); + end - 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); - 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 - - [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 + for j = 1:j_ + sir = wh.parameter.sir.values(j); + for i = 1:i_ + rop=wh.parameter.rop.values(i); + + 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); + 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 + 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_ffe = wh.getStoValue('ber_ffe',M,datarate,wh.parameter.rop.values); - +ber_mlse=squeeze(mean(ber_mlse,1)); +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)'); diff --git a/projects/Lab_2024/offline_dsp_analysis/load_n_dsp.m b/projects/Lab_2024/offline_dsp_analysis/load_n_dsp.m index 6e8c897..2c8baab 100644 --- a/projects/Lab_2024/offline_dsp_analysis/load_n_dsp.m +++ b/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)]); - Scpe_sig = Signal_cell.S{1}; + Signal_cell = load([foldername, filesep, filename, 'rx_signal_iatten_',num2str(atten)]); + + % Scpe_sig = all_signals(cnt); - bers = wh.getStoValue('ber',v_bias,awg_vpp,eq_mode,atten); + 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, Noi] = Eq.process(Scpe_sig,Duobinary().encode(Symbols)); + + % EQ_sig.plot("fignum",50,"displayname",'After EQ','clear',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); + + Rx_bits = PAMmapper(M,0).demap(EQ_sig); + [~,errors_bm,ber_db_(i),errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + + disp([' DB Precode -> Channel -> FFE -> Decode/ Mod ',sprintf('%.1E',ber_db_(i)),' ']); - [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,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','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); - [~,errors_bm,ber,errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - - disp(['FFE: ',sprintf('%.1E',ber_ffe_only),' -> PF -> MLSE: ',sprintf('%.1E',ber),' dB | PD_in: ',num2str(pd_in),' dBm']); + elseif eqmode == 3 + [ber_db(cnt),idx] = min(ber_db_); + end + + 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 -% 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" +hLine3 = plot(sirs, bers, "LineWidth", 0.5, "LineStyle", "-", "Marker", ".", "MarkerSize", 15, "DisplayName", 'from measurement'); % Continue with the rest of your plot settings diff --git a/projects/standard_system/400G_simulative_setup.mat b/projects/standard_system/400G_simulative_setup.mat index 0202264..4716e8f 100644 Binary files a/projects/standard_system/400G_simulative_setup.mat and b/projects/standard_system/400G_simulative_setup.mat differ