folderpath = 'C:\Users\sioe\Documents\High_Speed_Measurement_2024\bias_5km\'; experiment_name = 'PAMX_5km_'; currentTime = datetime('now', 'Format', 'yyyyMMdd_HHmmss'); timeStr = char(currentTime); experiment_name = [experiment_name, timeStr]; ffe_only = 1; postfilter_approach = 0; db_channel_approach = 0; db_coding_approach = 0; db_precode = db_coding_approach || db_channel_approach; %%% SIR Sweep for MPI Experiment %%% params = struct; params.vbias = [1.7:0.02:3.2]; % PAM6=2.3V %PAM8=2.68V params.awg_vpp = [2.7]; params.precomp_amp_max = [5]; params.rop_atten = [0]; params.M = [4,6,8]; params.lambda = [1292,1310,1327]; %calcWavelengthPlan(16, 400e9 , 1310); % params.vbias = [2.23]; % PAM6=2.3V %PAM8=2.68V % params.awg_vpp = [2.7]; % params.precomp_amp_max = [5]; % params.rop_atten = [0]; % params.M = [6]; % params.lambda = [1310]; %calcWavelengthPlan(16, 400e9 , 1310); wh = DataStorage(params); wh.addStorage("ber_collect"); wh.addStorage("ber_ffe"); wh.addStorage("ber_mlse"); wh.addStorage("ber_db"); wh.addStorage("pd_in"); wh.addStorage("rop"); wh.addStorage("m"); wh.addStorage("rx_logbook"); precomp_path = "C:\Users\sioe\Documents\High_Speed_Measurement_2024\precomp\"; precomp_fn = "lab_high_speed"; precomp_mode = 2; %0=do nothing ; 1= measure; 2=precomp active precomp_amp_max = 4; random_key = 2; pd_in_set = 7; looptotal = prod(wh.dim); disp(['Start Measurement of ',num2str(looptotal),' loops...']) iterationTimes = zeros(looptotal, 1); % Preallocate for speed if ~exist('hWaitbar', 'var') || ~isvalid(hWaitbar) hWaitbar = waitbar(0, sprintf('Starting %d measurements',looptotal), 'Name', 'Processing Progress'); else waitbar(0, hWaitbar, sprintf('Starting %d measurements',looptotal)); end loopcnt = 0; estimatedTimeRemaining = 0; estimatedTotalTime = 0; for lambda = wh.parameter.lambda.values exfo = Exfo_laser("serialport_number",'COM8','mainframe_channel',1,'safety_mode',0); pdfa = Thor_PDFA("safety_mode",0); exfo.getLaserInfo; if ~(exfo.cur_wavelength == lambda) % 1) pdfa.disablePDFA; % 2) exfo.setWavelength(lambda); end % 3) pdfa.enablePDFA(); % 4) pdfa.setPumpLevel(100); % 5) SET to first vbias and wait 30 minutes v_bias_first = wh.parameter.vbias.values(1); dcs = DC_supply("active",[1,0],"voltage",[v_bias_first, 0]); dcs.set("voltage",[v_bias_first, 0]); % pause(30*60); %wait 30 minutes for stable bias for v_bias = wh.parameter.vbias.values for rop_atten = wh.parameter.rop_atten.values for precomp_amp_max = wh.parameter.precomp_amp_max.values for M = wh.parameter.M.values for awg_vpp = wh.parameter.awg_vpp.values iterationStartTime = tic; loopcnt = loopcnt+1; if M == 4 fsym = 168e9; pulsef = 1; elseif M == 6 fsym = 144e9; pulsef = 0; elseif M == 8 fsym = 120e9; pulsef = 0; end %%%%% Loop Preps %fsym = round(targetrate/log2(M)); loop_name = ['_fsym_',num2str(fsym)]; %%%%% SET Voltages %%%%%% dcs = DC_supply("active",[1,0],"voltage",[v_bias, 0]); dcs.set("voltage",[v_bias, 0]); %%%%% SET Attenuator %%%%%% voa = OptAtten("active",[1,2,1,1],"value",[rop_atten,pd_in_set,0,0],"wavelength",[1310,1310,1310,1310]); voa.set('active',[1,2,1,1],'value',[rop_atten,pd_in_set,0,0]); % voa.readvals(); %%%%% Construct AWG and Scope Modules %%%%%% fdac = 256e9; fadc = 256e9; SCP = ScopeKeysight("model","UXR1104B",'autoscale',1,"fadc","GSa_256","channel",[0,1,0,0],"recordLen",2000000,"removeDC",1); AWG = AwgKeysight("model","M8199B","fdac",fdac,"scaletodac",[1,1],"skews",[0,0],"voltages",[0,awg_vpp]); A2S = Awg2Scope(AWG,SCP,[0,2,0,0],"waitUntilClick",0); % %%%%% Symbol Generation %%%%%% Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"rrcalpha",0.05); [Digi_sig,Symbols,Bits] = PAMsource(... "fsym",fsym,"M",M,"order",19,"useprbs",1,... "fs_out",fdac,... "applyclipping",0,"clipfactor",1.5,... "applypulseform",pulsef,"pulseformer",Pform,... "randkey",random_key,... "db_precode",db_precode,"db_encode",db_coding_approach,... "mrds_code",0,"mrds_blocklength",512).process(); Digi_sig.spectrum("displayname","Normal Tx","fignum",20,"normalizeToNyquist",0,"normalizeTo0dB",0); %%%%% Precompensation Routine %%%%%% if precomp_mode == 1 % measure channel precomp_est = ChannelFreqResp("Nacq",1024,"Navg",64,"Ncp",63,'f_ref',fdac); Digi_sig = precomp_est.buildOFDM(); elseif precomp_mode == 2 % apply precomp precomp_est = ChannelFreqResp("Nacq",1024,"Navg",64,"Ncp",63,'f_ref',Digi_sig.fs); Digi_sig = precomp_est.precomp(Digi_sig,'maxampdb',precomp_amp_max,'loadPath',precomp_path,'fileName',precomp_fn); end %%%%% Resample to DAC rate %%%%%% Digi_sig = Digi_sig.resample("fs_out",AWG.fdac); %%%%% Plot and Save Routine 1 %%%%%%%%%%%%%%%%%%%%%%%%% Digi_sig.spectrum("displayname","Normal Tx","fignum",10,"normalizeToNyquist",0,"normalizeTo0dB",0); save([folderpath,[experiment_name,'_bits'],loop_name],"Bits"); save([folderpath,[experiment_name,'_symbols'],loop_name],"Symbols"); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%% AWG --> Scope %%%%%% [~,Scpe_sig_raw,~,D] = A2S.process("signal2",Digi_sig,"waitUntilClick",0); % Scpe_sig_raw.spectrum("displayname","Scope PSD","fignum",20); % Scpe_sig_raw.plot("displayname","Scope raw signal","fignum",25,"clear",1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%% Sample to 2x fsym %%%%%% Scpe_sig_resampled = Scpe_sig_raw.resample("fs_in",fadc,"fs_out",2*fsym); %%%%% Precompensation Routine %%%%%% if precomp_mode == 1 precomp_est.estimate(Scpe_sig_resampled,"save",true,"savePath",precomp_path,"fileName",precomp_fn); precomp_est.plot(); end voa.readvals(); rop = voa.power_state(1); pd_in = voa.power_state(2); disp(['ROP: ',num2str(rop),' dBm || PD in: ',num2str(pd_in), ' dBm']); %%%%%% Sync Rx signal with reference (S is a cell array with all occurences) %%%%%% [Scpe_sig_syncd,S,isFlipped] = Scpe_sig_resampled.tsynch("reference",Symbols,"fs_ref",fsym); %%%%%% SNR CHEAT - Avges the measured signal occurences found after correlation in "tsynch" %%%%%% average_signals = 0; if average_signals scope_mean = zeros(size(S{1}.signal)); for n=1:numel(S) scope_mean = scope_mean + S{n}.signal; end scope_mean = scope_mean ./ n; Scpe_sig_syncd.signal = scope_mean; end %%%%% Plot and Save Routines: SAVE RECEIVED SIGNALS %%%%%%%%%%%%%%%%%%%%%%%%% save([folderpath,experiment_name,'_rx_signal',loop_name],"S"); Scpe_sig_syncd.eye(fsym,M,"fignum",40,"displayname",' after Scope'); %%%%% EQUALIZE %%%%%% Eq = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",50,"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); % set to minus one not zero not avoid confusion if BER is acutally zero ber_ffe = -1; ber_mlse = -1; ber_db = -1; if ffe_only %%%%%%%%%%%%%%%%%%%%%%%%%%% ber = []; parfor i = 1:numel(S) [EQ_sig] = Eq.process(S{i},Symbols); % EQ_sig.plot("fignum",50,"displayname",'After EQ'); Rx_bits = PAMmapper(M,0).demap(EQ_sig); [~,errors_bm,ber(i),errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); disp(['FFE: ',sprintf('%.1E',ber(i)),'| ROP: ',num2str(rop),' dB | PD_in: ',num2str(pd_in),' dBm']); end try ber_ffe = mean(rmoutliers(ber)); catch ber_ffe = min(ber); end if 0 EQ_sig.plot("fignum",50,"displayname",'After EQ','clear',1); figure(56); clf title(sprintf('PAM %d ; BER: %1.2e',M, ber_ffe)); constellation = unique(Symbols.signal); received = NaN(numel(constellation),length(Symbols)); for lvl = 1:numel(constellation) %Separate the equalized signal into the %respective levels based on the actually %transmitted level! received(lvl,Symbols.signal==constellation(lvl)) = EQ_sig.signal(Symbols.signal==constellation(lvl)); intermediate = received(lvl,:); cnt(lvl) = numel(intermediate(~isnan(intermediate))); hold on histogram(received(lvl,:),1000,"EdgeAlpha",0,'DisplayName',['Lvl ',num2str(lvl),' | ',num2str(cnt(lvl)),' entries']); end legend end elseif postfilter_approach %%%%%%%%%%%%%%%%%%%%%%%%%%% [EQ_sig] = Eq.process(Scpe_sig_syncd,Symbols); EQ_sig.plot("fignum",50,"displayname",'After EQ','clear',1); Noi = EQ_sig-Symbols; Rx_bits = PAMmapper(M,0).demap(EQ_sig); [~,num_errors,ber_ffe,pos_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 0 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 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); [~,num_errors,ber_mlse,pos_errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); disp(['FFE: ',sprintf('%.1E',ber_ffe),' -> PF -> MLSE: ',sprintf('%.1E',ber_mlse),' dB | PD_in: ',num2str(pd_in),' dBm']); elseif db_channel_approach %%%%%%%%%%%%%%%%%%%%%%%%%%% [EQ_sig, Noi] = Eq.process(Scpe_sig_syncd,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); [~,num_errors,ber_db,pos_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),' | PD_in: ',num2str(pd_in),' dBm']); elseif db_coding_approach %%%%%%%%%%%%%%%%%%%%%%%%%%% [EQ_sig, Noi] = Eq.process(Scpe_sig_syncd,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); [~,errors_bm,ber_db,errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); EQ_sig.plot("fignum",50,"displayname",'After EQ'); disp([' DB Precode -> DB Code -> Channel -> FFE -> Decode/ Mod ',sprintf('%.1E',ber_db),' | PD_in: ',num2str(pd_in),' dBm']); end %%%%% Store measurement into measurement "warehouse" %%%%%% wh.addValueToStorage(ber,'ber_collect',v_bias,awg_vpp,precomp_amp_max,rop_atten,M,lambda); wh.addValueToStorage(ber_ffe,'ber_ffe',v_bias,awg_vpp,precomp_amp_max,rop_atten,M,lambda); wh.addValueToStorage(ber_mlse,'ber_mlse',v_bias,awg_vpp,precomp_amp_max,rop_atten,M,lambda); wh.addValueToStorage(ber_db,'ber_db',v_bias,awg_vpp,precomp_amp_max,rop_atten,M,lambda); wh.addValueToStorage(rop,'rop',v_bias,awg_vpp,precomp_amp_max,rop_atten,M,lambda); wh.addValueToStorage(pd_in,'pd_in',v_bias,awg_vpp,precomp_amp_max,rop_atten,M,lambda); Rx_bits.logbook.SignalCopy = []; wh.addValueToStorage(Rx_bits,'rx_logbook',v_bias,awg_vpp,precomp_amp_max,rop_atten,M,lambda); wh.addValueToStorage(M,'m',v_bias,awg_vpp,precomp_amp_max,rop_atten,M,lambda); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%% Plot stuff into Table (feel free to add own values in -> 'name',value <- notation. Must be closed when table is changed)%%%%%%%%%%%%%%%%%%%%%% showCurrentMeasurement('FFE', ber_ffe,'MLSE',ber_mlse, 'Fsym',fsym.*1e-9, 'ROP', rop, 'PD in', pd_in, 'PAM',M, 'Vbias', v_bias, 'AWG Vpp', awg_vpp, 'Precomp MaxAmp',precomp_amp_max); %%%%% Arrange Figures %%%%%%%%%%%%%%%%%%%%%% autoArrangeFigures(3,3,2); iterationTimes(loopcnt) = toc(iterationStartTime); averageTimePerIteration = mean(iterationTimes(1:loopcnt)); estimatedTotalTime = averageTimePerIteration * looptotal; estimatedTimeRemaining = estimatedTotalTime - sum(iterationTimes(1:loopcnt)); progressFraction = loopcnt / looptotal; waitbar(progressFraction, hWaitbar, ... sprintf('Loop: %d of %d \n Runtime: %.1f min | %.1f sec per Loop |Time to go: %.1f min ', ... loopcnt, looptotal, sum(iterationTimes(1:loopcnt))/60, averageTimePerIteration, estimatedTimeRemaining/60 )); wh.save([folderpath,experiment_name,'_wh']); end end end end end end close(hWaitbar); wh.save([folderpath,experiment_name,'_wh']); autoArrangeFigures(3,3,2) %%% LAMBDA PLOT if 0 lambda_vals = wh.parameter.lambda.values; ber_ffe = wh.getStoValue('ber_ffe',v_bias,awg_vpp,precomp_amp_max,rop_atten,M,lambda_vals); ber_mlse = wh.getStoValue('ber_mlse',v_bias,awg_vpp,precomp_amp_max,rop_atten,M,lambda_vals); rop_measured = wh.getStoValue('rop',v_bias,awg_vpp,precomp_amp_max,rop_atten,M,lambda_vals); pd_in_measured = wh.getStoValue('pd_in',v_bias,awg_vpp,precomp_amp_max,rop_atten,M,lambda_vals); figure(240); hold on legendname = ['5 km Fsym:',num2str(fsym.*1e-9),' GBd | PAM',num2str(M),' | PD: ',num2str(pd_in_set),'| Vbias: ', num2str(v_bias),'V | ']; ffeLine = plot(lambda_vals, ber_ffe, "LineWidth", 0.5, "LineStyle", "-", "Marker", ".", "MarkerSize", 15, "DisplayName", [legendname,' + VNLE']); % Continue with the rest of your plot settings yline(3.8e-3, 'DisplayName', 'HD-FEC', 'LineStyle', '--', 'HandleVisibility', 'off'); xlabel('Wavelength in nm'); ylabel('Bit Error Rate (BER)'); title('Bit Error Rate vs. ROP'); set(gca, 'yscale', 'log'); set(gca, 'Box', 'on'); grid on; grid minor; legend('Interpreter', 'none'); end %%% ROP PLOT if 0 rop_vals = wh.parameter.rop_atten.values; ber_ffe = wh.getStoValue('ber_ffe',v_bias,awg_vpp,precomp_amp_max,rop_vals,M); ber_mlse = wh.getStoValue('ber_mlse',v_bias,awg_vpp,precomp_amp_max,rop_vals,M); rop_measured = wh.getStoValue('rop',v_bias,awg_vpp,precomp_amp_max,rop_vals,M); pd_in_measured = wh.getStoValue('pd_in',v_bias,awg_vpp,precomp_amp_max,rop_vals,M); legendname = ['Thormax: ',num2str(fsym.*1e-9),' GBd | PAM',num2str(M),' | PD: ',num2str(pd_in_set),'| Vbias: ', num2str(v_bias),'V | ']; figure(230); hold on; % Retain the plot so new points can be added without complete redraw % Plot the data and get the line handle ffeLine = plot(rop_measured, ber_ffe, "LineWidth", 0.5, "LineStyle", "-", "Marker", ".", "MarkerSize", 15, "DisplayName", [legendname,' + FFE']); mlseLine = plot(rop_measured, ber_mlse, "LineWidth", 0.5, "LineStyle", "-", "Marker", ".", "MarkerSize", 15, "DisplayName", [legendname,' +MLSE']); % Store pd_in_measured in the ZData property ffeLine.ZData = pd_in_measured; % Customize the data tips % Set labels for existing data tip rows ffeLine.DataTipTemplate.DataTipRows(1).Label = 'ROP'; ffeLine.DataTipTemplate.DataTipRows(2).Label = 'FFE'; ffeLine.DataTipTemplate.DataTipRows(2).Format = '%.2e'; % Format BER as "3e-4" % Add a new data tip row for PDin pdinRow = dataTipTextRow('PDin', 'ZData'); ffeLine.DataTipTemplate.DataTipRows(3) = pdinRow; % Store pd_in_measured in the ZData property mlseLine.ZData = pd_in_measured; % Customize the data tips % Set labels for existing data tip rows mlseLine.DataTipTemplate.DataTipRows(1).Label = 'ROP'; mlseLine.DataTipTemplate.DataTipRows(2).Label = 'MLSE'; mlseLine.DataTipTemplate.DataTipRows(2).Format = '%.2e'; % Format BER as "3e-4" % Add a new data tip row for PDin pdinRow = dataTipTextRow('PDin', 'ZData'); mlseLine.DataTipTemplate.DataTipRows(3) = pdinRow; % Continue with the rest of your plot settings 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('Interpreter', 'none'); end