imdd_silas/projects/MPI_April/mpi_simulation_cspr.m

%% Settings
clear

filename = '112G_2';
load_sequence = 0;

M = 4;
datarate = 112e9;

kover = 8;
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;
awg_bw = lowpass_cutoff;
mod_bw = lowpass_cutoff;
phd_bw = lowpass_cutoff;
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_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

            load(['C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\MPI_April\',filename,'.mat'],'X');

            spectrum_plot(X.signal',X.fs,'spectrum');
        else

            % 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);
            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);

        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');