From cb54f4480e7ac83bde161545c0f2d7b58a7a46fb Mon Sep 17 00:00:00 2001 From: Silas Oettinghaus Date: Mon, 22 Apr 2024 11:14:06 +0200 Subject: [PATCH] Update --- Classes/00_signals/Signal.m | 28 ++ Classes/04_DSP/EQ.m | 4 +- Classes/04_DSP/EQ_silas.m | 71 ++-- .../EQ_silas_sliding_window_dc_removal.m | 39 +- Functions/calc_ber.m | 2 +- projects/MPI_April/auswertung_1.m | 51 +++ projects/MPI_April/mpi_simulation_cspr.m | 366 ++++++++++-------- projects/MPI_April/mzm_er.m | 9 + projects/MPI_April/plot_analysis_window.m | 8 +- projects/MPI_April/save_to_warehouse.m | 61 +++ 10 files changed, 428 insertions(+), 211 deletions(-) create mode 100644 projects/MPI_April/auswertung_1.m create mode 100644 projects/MPI_April/mzm_er.m create mode 100644 projects/MPI_April/save_to_warehouse.m diff --git a/Classes/00_signals/Signal.m b/Classes/00_signals/Signal.m index 05a309d..a50fd5e 100644 --- a/Classes/00_signals/Signal.m +++ b/Classes/00_signals/Signal.m @@ -364,6 +364,34 @@ classdef Signal end + function eye(obj,fsym) + + disp('h'); + + fsig = obj.fs; + q = fsig/fsym; + + if q > 10 + sig = abs(obj.signal).^2; + else + sig = abs(obj.resample("fs_in",fsig,"fs_out",fsym*10).signal); + q = 10; + end + + figure(12) + clf + cursor = 200*q; + + for s = 1:700 + plot(sig(cursor-q:cursor+q),'Color','black','LineWidth',0.1,'LineStyle','-'); + hold on + cursor=cursor+q; + s=s+1; + end + + + end + end end diff --git a/Classes/04_DSP/EQ.m b/Classes/04_DSP/EQ.m index c3facf0..b999c69 100644 --- a/Classes/04_DSP/EQ.m +++ b/Classes/04_DSP/EQ.m @@ -384,7 +384,7 @@ classdef EQ end for dd_loop = 1:obj.dd_loops - cnt = obj.training_length+1; + cnt = 1; m = 0; output_vec = zeros(1,floor(length(data_in)/obj.K)); % initilaization of the output vector dd_DFE = zeros(obj.Nb(1),1); @@ -452,7 +452,7 @@ classdef EQ % e_save(:,save_ind) = coeff; % save_ind = save_ind+1; - if mu_mat ~= 0 + if 1%mu_mat ~= 0 e_dc = e_dc - obj.DCmu*error; error_log(cnt,dd_loop) = e_dc; cnt = cnt+1; diff --git a/Classes/04_DSP/EQ_silas.m b/Classes/04_DSP/EQ_silas.m index 497dbe9..39cbfa3 100644 --- a/Classes/04_DSP/EQ_silas.m +++ b/Classes/04_DSP/EQ_silas.m @@ -111,12 +111,16 @@ classdef EQ_silas < handle % actual processing of the signal (steps 1. - 3.) % 1 normalize RMS - %signalclass_in = signalclass_in.normalize("mode","rms"); + signalclass_in = signalclass_in.normalize("mode","rms"); % Process the EQ optimization obj.process_(signalclass_in.signal', reference_signalclass_in.signal'); signalclass_in.signal = obj.y_out'; + + %change sampling frequency of outgoing signal + signalclass_in.fs = reference_signalclass_in.fs; + % append to logbook lbdesc = ['EQ von Silas ist gelaufen ']; signalclass_in = signalclass_in.logbookentry(lbdesc); @@ -198,6 +202,7 @@ classdef EQ_silas < handle dc_cnt = 0; end + end end @@ -209,11 +214,9 @@ classdef EQ_silas < handle %start the dd mode with coefficients from training coeff = [obj.e;obj.b]; - obj.e_dc = ones(obj.eq_updatelatency,1).*obj.e_dc; + obj.e_dc = ones(obj.eq_updatelatency,1);%.*obj.e_dc; dc_block = ones(obj.eq_parallelization_blocklength,1); - - for ddloop = 1:obj.ddloops m = 0; @@ -261,7 +264,6 @@ classdef EQ_silas < handle %combine FFE with DFE to one vector (cursor between the two sequences) x_d = [x_vnle;-d_vnle]; - %Apply filter %y(m) = (m_reg(end)*dc_cnt + obj.e_dc(end)) + x_d.'* coeff; if obj.mu_dc_dd > 0 @@ -270,13 +272,6 @@ classdef EQ_silas < handle y(m) = x_d.'* coeff; end -% if obj.eq_avg_blocklength > 0 %% Das läuft nicht gut!! -% averaging_window = circshift(averaging_window,obj.sps); -% averaging_window(1:obj.sps,1) = y(m); -% avg_(m) = mean(averaging_window); -% y(m) = y(m)-avg_(m); -% end - %Decision [~,symbol_idx] = min(abs(y(m) - obj.d_constellation)); % decision for closest constellation point d_hat(k) = obj.d_constellation(symbol_idx); @@ -284,15 +279,8 @@ classdef EQ_silas < handle %Error between FFE & DFE filtered signal and Decision obj.error(k) = y(m) - d_hat(k); -% if obj.eq_avg_blocklength > 0 %% Das läuft nicht gut!! -% averaging_window = circshift(averaging_window,obj.sps); -% averaging_window(1:obj.sps,1) = y(m); -% avg_(m) = mean(averaging_window); -% y(m) = y(m)-avg_(m); -% end - %Update FFE and DFE coefficients - coeff = coeff - mu_mat*obj.error(k) * conj(x_d); + coeff = coeff - (mu_mat * (obj.error(k) * conj(x_d))); %Update DC error dc_block(dc_cnt) = obj.error(k) ; @@ -312,7 +300,9 @@ classdef EQ_silas < handle % obj.e_dc = obj.e_dc - sign(m_reg(end)) .* (sum(dc_block).* m_reg(end) .* obj.mu_dc_dd); - obj.e_dc = obj.e_dc - sum(dc_block) .* obj.mu_dc_dd; + obj.e_dc = obj.e_dc - sum(dc_block) .* obj.mu_dc_dd; + +% obj.e_dc = obj.e_dc - obj.mu_dc_dd * obj.error(k); %newapril end dc_cnt = 0; @@ -343,24 +333,43 @@ classdef EQ_silas < handle function x_in_vnle_format = calcVNLENonlinVecs(~,x_in_block,I_2,I_3,N_,norm_) % These are the second and third order input signal products of the VNLE EQ % ∑ h1 x_in(k-n1) + ∑∑ h2 x_in(k-n1)*x_in(k-n2) + ∑∑∑ h3 x_in(k-n1)*x_in(k-n2)*x_in(k-n3) + l1=length(x_in_block); + l2=length(I_2); + l3=length(I_3); + final_length = l1+l2+l3; - x1 = x_in_block; - x2 = []; - x3 = []; + x_in_vnle_format = zeros(final_length,1); + + idx = l1; + x_in_vnle_format(1:idx) = x_in_block; if N_(2) > 0 - delta_2 = round((N_(1)-N_(2))/2); - input_vec_se = x_in_block(delta_2:end)/norm_(2); %TODO normalization step - x2 = input_vec_se(I_2(:,1)).*input_vec_se(I_2(:,2)); + delta_2 = round((N_(1)-N_(2)) / 2); + input_vec_se = x_in_block(delta_2:end) / norm_(2); %TODO normalization step + + % Extract columns from I_2 + col1 = input_vec_se(I_2(:,1)); + col2 = input_vec_se(I_2(:,2)); + + x2 = col1 .* col2; + x_in_vnle_format(idx+1:idx+l2) = x2; end if N_(3) > 0 delta_3 = round((N_(1)-N_(3))/2); - input_vec_th = x_in_block(delta_3:end)/norm_(3); - x3 = input_vec_th(I_3(:,1)).*input_vec_th(I_3(:,2)).*input_vec_th(I_3(:,3)); - end + input_vec_th = x_in_block(delta_3:end) / norm_(3); - x_in_vnle_format = [x1;x2;x3]; + % Extract columns from I_3 + col1 = input_vec_th(I_3(:,1)); + col2 = input_vec_th(I_3(:,2)); + col3 = input_vec_th(I_3(:,3)); + + % Perform matrix multiplication + x3 = col1 .* col2 .* col3; + + idx = idx+l2; + x_in_vnle_format(idx+1:idx+l3) = x3; + end end diff --git a/Classes/04_DSP/EQ_silas_sliding_window_dc_removal.m b/Classes/04_DSP/EQ_silas_sliding_window_dc_removal.m index 0787434..cf1bde1 100644 --- a/Classes/04_DSP/EQ_silas_sliding_window_dc_removal.m +++ b/Classes/04_DSP/EQ_silas_sliding_window_dc_removal.m @@ -293,24 +293,43 @@ classdef EQ_silas_sliding_window_dc_removal < handle function x_in_vnle_format = calcVNLENonlinVecs(~,x_in_block,I_2,I_3,N_,norm_) % These are the second and third order input signal products of the VNLE EQ % ∑ h1 x_in(k-n1) + ∑∑ h2 x_in(k-n1)*x_in(k-n2) + ∑∑∑ h3 x_in(k-n1)*x_in(k-n2)*x_in(k-n3) + l1=length(x_in_block); + l2=length(I_2); + l3=length(I_3); + final_length = l1+l2+l3; - x1 = x_in_block; - x2 = []; - x3 = []; + x_in_vnle_format = zeros(final_length,1); + + idx = l1; + x_in_vnle_format(1:idx) = x_in_block; if N_(2) > 0 - delta_2 = round((N_(1)-N_(2))/2); - input_vec_se = x_in_block(delta_2:end)/norm_(2); %TODO normalization step - x2 = input_vec_se(I_2(:,1)).*input_vec_se(I_2(:,2)); + delta_2 = round((N_(1)-N_(2)) / 2); + input_vec_se = x_in_block(delta_2:end) / norm_(2); %TODO normalization step + + % Extract columns from I_2 + col1 = input_vec_se(I_2(:,1)); + col2 = input_vec_se(I_2(:,2)); + + x2 = col1 .* col2; + x_in_vnle_format(idx+1:idx+l2) = x2; end if N_(3) > 0 delta_3 = round((N_(1)-N_(3))/2); - input_vec_th = x_in_block(delta_3:end)/norm_(3); - x3 = input_vec_th(I_3(:,1)).*input_vec_th(I_3(:,2)).*input_vec_th(I_3(:,3)); - end + input_vec_th = x_in_block(delta_3:end) / norm_(3); - x_in_vnle_format = [x1;x2;x3]; + % Extract columns from I_3 + col1 = input_vec_th(I_3(:,1)); + col2 = input_vec_th(I_3(:,2)); + col3 = input_vec_th(I_3(:,3)); + + % Perform matrix multiplication + x3 = col1 .* col2 .* col3; + + idx = idx+l2; + x_in_vnle_format(idx+1:idx+l3) = x3; + end end diff --git a/Functions/calc_ber.m b/Functions/calc_ber.m index 046cc0c..fb5a9af 100644 --- a/Functions/calc_ber.m +++ b/Functions/calc_ber.m @@ -48,7 +48,7 @@ end data_ = logical(data(skipstart+1:end-skip_end,:))'; - delta_bits = length(reference) - length(data_); + delta_bits = length(reference) - length(data); skip_end = max(skip_end,delta_bits); reference_ = logical(reference(skipstart+1:end-skip_end,:))'; diff --git a/projects/MPI_April/auswertung_1.m b/projects/MPI_April/auswertung_1.m new file mode 100644 index 0000000..e3ca3a7 --- /dev/null +++ b/projects/MPI_April/auswertung_1.m @@ -0,0 +1,51 @@ + + + + +vp = wh.parameter.vp.values(1); +vb = wh.parameter.vb.values(1); +rop = wh.parameter.rop.values; + + +f=figure(1113); +tiledlayout(2,4) +for sir = [20,36] + for lw = wh.parameter.laser_linewidth.values + nexttile + cols = linspecer(9); + cnt = 1; + for bias = wh.parameter.vb.values + curber = []; + curstd = []; + rop_meas = []; + txpow_meas = wh.getStoValue('mod_out_pow',sir,lw,1,vp,bias,rop(1)); + for pn_key = wh.parameter.pn_key.values + curber(end+1,:) = wh.getStoValue('ber',sir,lw,pn_key,vp,bias,rop); + rop_meas(end+1,:) = wh.getStoValue('rop_save',sir,lw,pn_key,vp,bias,rop); + + curstd(end+1,:,:) = wh.getStoValue('level_std',sir,lw,pn_key,vp,bias,rop); + plot(rop_meas(end,:)-txpow_meas, curber(end,:) ,'LineStyle',':','Color',cols(cnt,:),'LineWidth',0.1,'Marker','o','MarkerEdgeColor',[1 1 1],'MarkerFaceColor',cols(cnt,:),'HandleVisibility','off','MarkerSize',1); + hold on + end + [wrst,idx]=max(curber); + for i = 1:numel(idx) + rop_wrst(i)=rop_meas(idx(i),i); + end + plot( rop_wrst-txpow_meas, wrst ,'LineStyle',':','Color',cols(cnt,:),'LineWidth',1,'Marker','^','MarkerSize',5,'MarkerEdgeColor',[1 1 1],'MarkerFaceColor',cols(cnt,:),'HandleVisibility','off'); + plot(mean(rop_meas,1)-txpow_meas, mean(curber,1),'DisplayName',['Vbias: ',num2str(bias),' V'],'LineStyle','-','Color',cols(cnt,:),'LineWidth',1,'Marker','o','MarkerEdgeColor',[1 1 1],'MarkerFaceColor',cols(cnt,:)); + hold on + cnt = cnt+1; + end + + title(['SIR: ',num2str(sir),'; Lw: ',num2str(lw*1e-6),' MHz; Vpeakpeak: ',num2str(2*vp)]) + set(gca,'YScale','log'); + legend('Location','southwest') + xlabel("measured ROP in dBm") + ylabel("BER") + yline(3.8e-3,'DisplayName','FEC'); + xlim([-5,3]); + ylim([1e-4,3e-1]); + drawnow + end +end + diff --git a/projects/MPI_April/mpi_simulation_cspr.m b/projects/MPI_April/mpi_simulation_cspr.m index a9962e1..ca1f15c 100644 --- a/projects/MPI_April/mpi_simulation_cspr.m +++ b/projects/MPI_April/mpi_simulation_cspr.m @@ -4,20 +4,16 @@ clear filename = '112G_2'; load_sequence = 0; -M = 4; -datarate = 112e9; +M = 8; +datarate = 448e9; -kover = 8; +kover = 4; fsym = round(datarate*1e-9 / log2(M))*1e9; %fsym = 50e9; fdac = 256e9; fadc = 256e9; -laser_linewidth = 1e6; -mpi_ = 50; %meter -sir = 20; %decibel = attenuation of interference path - -lowpass_cutoff = fsym/2 * 1; +lowpass_cutoff = fsym/2 * 1.5; awg_bw = lowpass_cutoff; mod_bw = lowpass_cutoff; phd_bw = lowpass_cutoff; @@ -25,174 +21,218 @@ scp_bw = lowpass_cutoff; mpi_path=50; -vp = [0.25]; -vb = [0.5:0.05:0.8]; -vb = 0.7; -pn_key = [27]; -rop = -5; - -LP_awg = Filter('filtdegree',2,"f_cutoff",awg_bw,"fs",fdac,"filterType","butterworth"); +LP_awg = Filter('filtdegree',2,"f_cutoff",awg_bw,"fs",fdac*kover,"filterType","butterworth"); LP_laser = Filter('filtdegree',4,"f_cutoff",mod_bw,"fs",fdac*kover,"filterType",filtertypes.gaussian); LP_opt = Filter('filtdegree',4,"f_cutoff",fsym/log2(M).*1.5,"fs",fdac*kover,"filterType",filtertypes.gaussian); LP_phd = Filter('filtdegree',2,"f_cutoff",phd_bw,"fs",fdac*kover,"filterType",filtertypes.butterworth); LP_scpe = Filter('filtdegree',2,"f_cutoff",scp_bw,"fs",fadc,"filterType","butterworth"); -for pnk = 1:length(pn_key) - for m = 1:length(vb) - if load_sequence +% 1) PRBS Generation +O = 18; %order of prbs +N = 2^(O-1); %length of prbs +[~,seed] = prbs(O,1); %initialize first seed of prbs +bitpattern=[]; - load(['C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\MPI_April\',filename,'.mat'],'X'); +for i = 1:log2(M) + [bitpattern(:,i),seed] = prbs(O,N,seed); +end - spectrum_plot(X.signal',X.fs,'spectrum'); - else +% 2 ) Build Inf. signal class +bits = Informationsignal(bitpattern); - % 1) PRBS Generation - O = 17; %order of prbs - N = 2^(O-1); %length of prbs - [~,seed] = prbs(O,1); %initialize first seed of prbs - bitpattern=[]; - for i = 1:log2(M) - [bitpattern(:,i),seed] = prbs(O,N,seed); +% 3) Digi modulation -> PAM-M signal +digimod_out = PAMmapper(M,0).map(bits); +digimod_out.fs = fsym; + +X = Pulseformer("fsym",fsym,"fdac",fdac,"pulse","rrc","pulselength",16,"rrcalpha",0.1).process(digimod_out); + +% 5) AWG (lowpass, quantization, sample and hold) +X = AWG("fdac",fdac,"dac_min",-1,"dac_max",1,"H_lpf",LP_awg,"kover",kover,"bit_resolution",5,"lpf_active",1).process(X); + +% 6) Lowpass behavior before laser +X = LP_laser.process(X); + +% 7) Normalize signal +X = X.normalize("mode","oneone"); + +sir = [20:2:36]; %decibel = attenuation of interference path + laser_linewidth = [1e5 1e6 3e6 10e6]; + pn_key = [1:10]; + vp = [0.25,0.5,0.75,1]; + vb = [1:0.1:1.8]; + + rop = -9:3; +% sir = 25; +% laser_linewidth = 1e6; +% pn_key = 1; +% vp = 0.8; +% vb = 1; + +cnt = 1; +for s = 1:length(sir) + for l = 1:length(laser_linewidth) + for pnk = 1:length(pn_key) + for n = 1:length(vp) + for m = 1:length(vb) + + cnt = cnt+1; + + % 1) Laser; Modulation -> OPTICAL DOMAIN + u_pi = 2; + vbias = -vb(m); + extmodlaser = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",X.fs,"lambda",1290,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth(l),"randomkey",pn_key(pnk)); + E = X.*vp(1); + + % E.signal = min(max(E.signal,-0.2),0.2); + + [Opt,extmodlaser] = extmodlaser.process(E); + + % figure(m) + % hold on + % scatter(E.signal(1:100000),(abs(Opt.signal(1:100000)).^2)*1e3,0.1,'.','DisplayName','Modulator TF') + % xlabel('Input in V') + % ylabel('abs(Output) in mW') + + + % ER = 10*log10(max(abs(Opt.signal).^2)/min(abs(Opt.signal).^2)); + + Opt = LP_opt.process(Opt); + + cspr(s,l,pnk,n,m) = Opt.cspr; + mod_out_pow(s,l,pnk,n,m) = Opt.power; + + % 2) ping pong fiber propagation + Interference_sig = Fiber("fsimu",Opt.fs,"fiber_length",mpi_path*2/1000,"alpha",0,"D",0,"lambda0",1310,"gamma",0).process(Opt); + + Interference_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",-sir(s)).process(Interference_sig); + + % In the meantime: delay the main signal + [Main_sig,dly] = Opt.delay("delay_meter",mpi_path*2); + + % Add + Combined_sig = Main_sig + Interference_sig; + + % Cut (due to the delays there is a jump in the signals) + if dly == 0;dly = 1;end + Combined_sig.signal = Combined_sig.signal(ceil(dly):end); + + % Fiber + Combined_sig = Fiber("fsimu",Combined_sig.fs,"fiber_length",2,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.08).process(Combined_sig); + + parfor i = 1:length(rop) + + % Set ROP + Rx_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",rop(i)).process(Combined_sig); + rop_save(s,l,pnk,n,m,i) = Rx_sig.power; + + % Square Law + Rx_sig = Photodiode("fsimu",Rx_sig.fs,"dark_current",2e-08,"responsivity",1,"temperature",20).process(Rx_sig); + + %Lowpass PhDiode + Rx_sig = LP_phd.process(Rx_sig); + + % Scope + Rx_sig = Scope("fsimu",Rx_sig.fs,"fadc",fadc,... + "delay",0,"fixed_delay",0,"lpf_bw",scp_bw,"filtertype",filtertypes.butterworth,... + "samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,... + "adcresolution",16,"quantbuffer",0.1,'block_dc',1,'lpf_active',1,'H_lpf',LP_scpe).process(Rx_sig); + + % Sample to 2x fsym + Rx_sig = Rx_sig.resample("fs_in",fadc,"fs_out",2*fsym); + + % Sync Rx signal with reference + [Rx_sig,D,cuts] = Rx_sig.tsynch("reference",digimod_out,"fs_ref",fsym); + + % % % simple EQ (optimum mudc: 0.05 -> 0.005) + % EQ_sig = EQ_silas_plain("Ne",[20,8,8],"Nb",[2,0,0],"trainlength",4096,"mu_dc_dd",0.005,"mu_dc_train",0.05,... + % "mu_ffe_train",0.005,"mu_combined_dd",[0.0004 0.0006 0.0003 0.005],"ddloops",3,'trainloops',3,'sps',2).process(Rx_sig,digimod_out); + % +% EQ_sig = EQ("K",2,"plottrain",0,"plotfinal",0,... +% "training_length",4096,"training_loops",3,... +% "Ne",[20,8,8],"Nb",[2,0,0],... +% "DCmu",0.005,"DDmu",[0.0004 0.0006 0.0003 0.005],"DFEmu",0.005,"FFEmu",0.00,... +% "dd_loops",3,"epsilon",[10 100 1000 ],"M",2,... +% "thres",[0.005 0.004 0.0005 ],"l1act",0,"delay",0,"rho",0.0005,"ideal_dfe",0,"DB_aim",0).process(Rx_sig,digimod_out); + + EQ_sig = EQ_silas("Ne",[20,2,0],"Nb",[2,0,0],"trainlength",4096,... + "sps",2,... + "mu_dc_dd",0.00,... + "mu_dc_train",0.00,... + "mu_ffe_train",0.00,... + "mu_dfe_train",0.005,... + "mu_ffe_dd",[0.0004 0.0006 0.0003],... + "mu_dfe_dd",0.005,... + "ddloops",3,... + "trainloops",3,... + "eq_parallelization_blocklength",1, ... + "eq_updatelatency",1,... + "eq_avg_blocklength",0).process(Rx_sig,digimod_out); + + % Demap + Rx_Bits = PAMmapper(M,0).demap(EQ_sig); + + Rx_symboldecision = PAMmapper(M,0).decide_pamlevel(EQ_sig,"symbol_levels",unique(digimod_out.signal)); + + levels = PAMmapper(M,0).separate_pamlevels(EQ_sig); + %levels = PAMmapper(M,0).separate_pamlevels(Rx_sig.resample("fs_in",Rx_sig.fs,"fs_out",fsym)); + + level_avg(s,l,pnk,n,m,i,:) = mean(levels,'omitnan'); + level_std(s,l,pnk,n,m,i,:) = std(levels,'omitnan'); + + % BER + [~,errors_bm,BER(s,l,pnk,n,m,i),errors] = calc_ber(Rx_Bits.signal,bitpattern,"skip_front",0,"skip_end",0,"returnErrorLocation",1); + + formatted_ber = sprintf('%.1e', BER(s,l,pnk,n,m,i)); + disp(['SIR: ',num2str(sir(s)),'; Lw:',num2str(laser_linewidth(l)),'; Key:',num2str(pn_key(pnk)),'; Vpeak: ',num2str(vp(n)),'; Vbias',num2str(vbias),'; BER: ',formatted_ber,'; run: ',num2str(cnt),' / 12961']); + +% plot_analysis_window; +% drawnow; + + end + end end - - % 2 ) Build Inf. signal class - bits = Informationsignal(bitpattern); - - % 3) Digi modulation -> PAM-M signal - digimod_out = PAMmapper(M,0).map(bits); - digimod_out.fs = fsym; - - X = Pulseformer("fsym",fsym,"fdac",fdac,"pulse","rrc","pulselength",16,"rrcalpha",0.01).process(digimod_out); - - % 5) AWG (lowpass, quantization, sample and hold) - %X = M8196A().process(X); - kover = 16; - X = AWG("fdac",fdac,"dac_min",-1,"dac_max",1,"H_lpf",LP_awg,"kover",kover).process(X); - - % 6) Lowpass behavior before laser - X = LP_laser.process(X); - - % 7) Normalize signal - X = X.normalize("mode","oneone"); - - %save(['C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\MPI_April\',char(filename)],'X'); - end - - % 1) Laser; Modulation -> OPTICAL DOMAIN - u_pi = 2; - vbias = -vb(m)*u_pi; - extmodlaser = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",X.fs,"lambda",1290,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",pn_key(pnk)); - E = X.*vp(1); - [Opt,extmodlaser] = extmodlaser.process(E); - -% figure(m) -% hold on -% scatter(E.signal(1:100000),(abs(Opt.signal(1:100000)).^2)*1e3,0.1,'.','DisplayName','Modulator TF') -% xlabel('Input in V') -% ylabel('abs(Output) in mW') - - Opt = LP_opt.process(Opt); - - cspr(m) = Opt.cspr; - - % 2) ping pong fiber propagation - Interference_sig = Fiber("fsimu",Opt.fs,"fiber_length",mpi_path*2/1000,"alpha",0,"D",0,"lambda0",1310,"gamma",0).process(Opt); - - Interference_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",-sir).process(Interference_sig); - - % In the meantime: delay the main signal - [Main_sig,n] = Opt.delay("delay_meter",mpi_path*2); - - % Add - Combined_sig = Main_sig + Interference_sig; - - % Cut (due to the delays there is a jump in the signals) - if n == 0;n = 1;end - Combined_sig.signal = Combined_sig.signal(ceil(n):end); - - % Fiber - Combined_sig = Fiber("fsimu",Combined_sig.fs,"fiber_length",0,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.08).process(Combined_sig); - - for i = 1:length(rop) - - % Set ROP - Rx_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",rop(i)).process(Combined_sig); - rop_save(m,i) = Rx_sig.power; - - % Square Law - Rx_sig = Photodiode("fsimu",Rx_sig.fs,"dark_current",2e-08,"responsivity",1,"temperature",20).process(Rx_sig); - - %Lowpass PhDiode - Rx_sig = LP_phd.process(Rx_sig); - - % Scope - Rx_sig = Scope("fsimu",Rx_sig.fs,"fadc",fadc,... - "delay",0,"fixed_delay",0,"lpf_bw",scp_bw,"filtertype",filtertypes.butterworth,... - "samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,... - "adcresolution",16,"quantbuffer",0.1,'block_dc',1,'lpf_active',1,'H_lpf',LP_scpe).process(Rx_sig); - - % Sample to 2x fsym - Rx_sig = Rx_sig.resample("fs_in",fadc,"fs_out",2*fsym); - - % Sync Rx signal with reference - [Rx_sig,D,cuts] = Rx_sig.tsynch("reference",digimod_out,"fs_ref",fsym); - - % % simple EQ (optimum mudc: 0.05 -> 0.005) - EQ_sig = EQ_silas_plain("Ne",[20,8,8],"Nb",[2,0,0],"trainlength",4096,"mu_dc_dd",0.005,"mu_dc_train",0.05,... - "mu_ffe_train",0.005,"mu_combined_dd",[0.0004 0.0006 0.0003 0.005],"ddloops",3,'trainloops',3,'sps',2).process(Rx_sig,digimod_out); - - % Demap - Rx_Bits = PAMmapper(M,0).demap(EQ_sig); - - Rx_symboldecision = PAMmapper(M,0).decide_pamlevel(EQ_sig,"symbol_levels",unique(digimod_out.signal)); - - levels = PAMmapper(M,0).separate_pamlevels(EQ_sig); - level_avg(:,m,i) = mean(levels,'omitnan'); - level_std(:,m,i) = std(levels,'omitnan'); - - % BER - [~,errors_bm,BER(m,i),errors] = calc_ber(Rx_Bits.signal,bitpattern,"skip_front",0,"skip_end",0,"returnErrorLocation",1); - - formatted_ber = sprintf('%.1e', BER(m,i)); - disp(formatted_ber); - + save('big_run_easy_EQ_one_nonlin'); + disp('saved_run2'); end end end -plot_analysis_window; -save(['C:\Users\Silas\Nextcloud4\Dokumente\02_Ablage_Office\MPI\Investigation_April_2024\','PAM_',num2str(M),'_mpi_',num2str(mpi_path),'_lw_',num2str(laser_linewidth)],"BER"); - -% BER plot -figure(340) -cols = linspecer(7); -for m = 1:size(BER,1) - hold on - plot(rop,BER(m,:),'DisplayName',['Bias: ',num2str(vb(m)), ' V'],'LineStyle','-','Color',cols(m,:),'LineWidth',1,'Marker','o','MarkerEdgeColor',[1 1 1],'MarkerFaceColor',cols(m,:)); - -end -set(gca,'YScale','log'); -legend -xlabel("ROP in dBm") -yline(3.8e-3,'DisplayName','FEC'); - - -figure(21) -hold on -plot(rop,mean(BER),'DisplayName',['Modulation: ',num2str(2*vp/extmodlaser.u_pi*100), ' $\%$'],'LineStyle','-','Color',cols(2,:),'LineWidth',1); -set(gca,'YScale','log'); -legend -ylabel("ROP in dBm") -yline(3.8e-3,'DisplayName','FEC'); - - - - -%check Rx and TX symbols -figure(101) -scatter(1:100,Rx_symboldecision.signal(1:100),10,'o'); -hold on -scatter(1:100,digimod_out.signal(1:100),5,'x'); +% save(['C:\Users\Silas\Nextcloud4\Dokumente\02_Ablage_Office\MPI\Investigation_April_2024\','PAM_',num2str(M),'_mpi_',num2str(mpi_path),'_lw_',num2str(laser_linewidth)],"BER"); +% +% % BER plot +% figure(340) +% cols = linspecer(7); +% for m = 1:size(BER,1) +% hold on +% plot(rop,BER(m,:),'DisplayName',['Bias: ',num2str(vb(m)), ' V'],'LineStyle','-','Color',cols(m,:),'LineWidth',1,'Marker','o','MarkerEdgeColor',[1 1 1],'MarkerFaceColor',cols(m,:)); +% +% end +% set(gca,'YScale','log'); +% legend +% xlabel("ROP in dBm") +% yline(3.8e-3,'DisplayName','FEC'); +% +% +% figure(21) +% hold on +% plot(rop,mean(BER),'DisplayName',['Modulation: ',num2str(2*vp/extmodlaser.u_pi*100), ' $\%$'],'LineStyle','-','Color',cols(2,:),'LineWidth',1); +% set(gca,'YScale','log'); +% legend +% ylabel("ROP in dBm") +% yline(3.8e-3,'DisplayName','FEC'); +% +% +% +% +% +% +% +% +% %check Rx and TX symbols +% figure(101) +% scatter(1:100,Rx_symboldecision.signal(1:100),10,'o'); +% hold on +% scatter(1:100,digimod_out.signal(1:100),5,'x'); diff --git a/projects/MPI_April/mzm_er.m b/projects/MPI_April/mzm_er.m new file mode 100644 index 0000000..b18bf13 --- /dev/null +++ b/projects/MPI_April/mzm_er.m @@ -0,0 +1,9 @@ + + +field=sqrt(10^(3/10-3)); %dbm to sqrt(mw) + +vpi = 2; +vbias = 1; +vin = 0.1; + +(field * cos(pi/2*(vin+vbias)/vpi)).^2 diff --git a/projects/MPI_April/plot_analysis_window.m b/projects/MPI_April/plot_analysis_window.m index 6d3f443..626938b 100644 --- a/projects/MPI_April/plot_analysis_window.m +++ b/projects/MPI_April/plot_analysis_window.m @@ -16,9 +16,9 @@ ylim([-100,0]); %Phase Investigation subplot(3,2,3) hold on -phase_int = extmodlaser.phase(n:end); -phase_main = delayseq(extmodlaser.phase,n); -phase_main = phase_main(n:end); +phase_int = extmodlaser.phase(ceil(n):end); +phase_main = delayseq(extmodlaser.phase,ceil(n)); +phase_main = phase_main(ceil(n):end); phase_diff = phase_int - phase_main; t = (1:length(phase_int))' ./ extmodlaser.fsimu ; plot(t*1e6,phase_int,'Color',colpairs(1,:),'DisplayName','Interferer Phase'); @@ -53,7 +53,7 @@ for h = 1:size(levels,2) scatter(t*1e6,levels(:,h),1,'.'); hold on end -yline(PAMmapper(M,0).thresholds); +%yline(PAMmapper(M,0).thresholds); title(['BER: ',sprintf('%.1e', BER(m,i))]); %Histogram of EQzed Signal diff --git a/projects/MPI_April/save_to_warehouse.m b/projects/MPI_April/save_to_warehouse.m new file mode 100644 index 0000000..f79e282 --- /dev/null +++ b/projects/MPI_April/save_to_warehouse.m @@ -0,0 +1,61 @@ + +sir = [20:2:36]; %decibel = attenuation of interference path +laser_linewidth = [1e5 1e6 3e6 10e6]; +pn_key = [1:10]; +vp = [0.25]; +vb = [1:0.1:1.8]; +rop = -9:3; + +params = struct; +params.sir = [20:2:36]; %decibel = attenuation of interference path +params.laser_linewidth = [1e5 1e6 3e6 10e6]; +params.pn_key = [1:10]; +params.vp = [0.25,0.5,0.75,1]; +params.vb = [1:0.1:1.8]; +params.rop = -9:3; +wh = DataStorage(params); +wh.addStorage("ber"); +wh.addStorage("level_avg"); +wh.addStorage("level_std"); +wh.addStorage("rop_save"); +wh.addStorage("cspr"); +wh.addStorage("mod_out_pow"); + +cnt = 1; +for s = 1:length(sir) + for l = 1:length(laser_linewidth) + for pnk = 1:length(pn_key) + for n = 1 + for m = 1:length(vb) + for i = 1:length(rop) + cnt = cnt +1; + + wh.addValueToStorage(BER(s,l,pnk,n,m,i),'ber',sir(s),laser_linewidth(l),pn_key(pnk),vp(n),vb(m),rop(i)); + wh.addValueToStorage(level_avg(s,l,pnk,n,m,i,:),'level_avg',sir(s),laser_linewidth(l),pn_key(pnk),vp(n),vb(m),rop(i)); + wh.addValueToStorage(level_std(s,l,pnk,n,m,i,:),'level_std',sir(s),laser_linewidth(l),pn_key(pnk),vp(n),vb(m),rop(i)); + wh.addValueToStorage(rop_save(s,l,pnk,n,m,i),'rop_save',sir(s),laser_linewidth(l),pn_key(pnk),vp(n),vb(m),rop(i)); + wh.addValueToStorage(cspr(s,l,pnk,n,m),'cspr',sir(s),laser_linewidth(l),pn_key(pnk),vp(n),vb(m),rop(i)); + wh.addValueToStorage(mod_out_pow(s,l,pnk,n,m),'mod_out_pow',sir(s),laser_linewidth(l),pn_key(pnk),vp(n),vb(m),rop(i)); + + end + end + end + end + end +end + +% 4) Hey! the warehouse is here and (hopefully) filled with data :-) + +% Create a save dialog +defaultDir = 'C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\MPI_April'; +defaultExt = '*.mat'; +[filename, pathname] = uiputfile(fullfile(defaultDir, defaultExt),'', 'wh.mat'); + +% Check if the user pressed Cancel +if isequal(filename, 0) || isequal(pathname, 0) + disp('Save operation canceled.'); +else + % Save the variable to the selected file + save(fullfile(pathname, filename), 'wh'); + disp(['Variable "wh" saved to: ', fullfile(pathname, filename)]); +end \ No newline at end of file