%% Parameter to simulate and save
params = struct;

params.M                        = [8];
params.datarate                 = [448];
params.rop                      = [-12:0];

precomp_mode    = 0; %0=do nothing ; 1= measure; 2=precomp active
postfilter      = 0; % noise whiten. approach -> Postfilter + MLSE

db_precode      = 1;
db_encode       = 1;
db_channelapproach = 0;


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\";
end

precomp_fn = "400G_simulative_setup";
usemrds = 0;

name = ['wh_',strrep(num2str(now),'.','')];

wh = DataStorage(params);

wh.addStorage("ber_ffe");

%% Init Params
link_length = 1000; %meter
pn_key = 2;
laser_linewidth = 0;

endcnt = prod(wh.dim);
cnt=0;

disp(['Start Simulation of ',num2str(endcnt),' loops...'])
tic

for M = wh.parameter.M.values
    for datarate = wh.parameter.datarate.values

        % SETUP HERE: %%
        kover = 16;
        M8199 = M8199B("kover",kover);
        fdac = M8199.fdac;
        fsym = round(datarate / log2(M)) * 1e9;
        rrcalpha = 0.05;
        Pform   = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"rrcalpha",rrcalpha);

        % MAIN SIGNAL

        %%%%% Symbol Generation %%%%%%
        [Digi_sig,Symbols,Bits] = PAMsource("fsym",fsym,"M",M,"order",18,"useprbs",0,...
            "fs_out",M8199.fdac,"applyclipping",1,"clipfactor",1.5,...
            "applypulseform",0,"pulseformer",Pform,"randkey",pn_key,...
            "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');

        if precomp_mode == 1
            freqresp = ChannelFreqResp("Nacq",1024,"Navg",64,"Ncp",63,'f_ref',Digi_sig.fs);
            Digi_sig = freqresp.buildOFDM();
        elseif precomp_mode == 2
            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.spectrum("fignum",11,"displayname",'after precomp');
        end
        
        %%%%% AWG %%%%%%
        El_sig = M8199.process(Digi_sig);

        % El_sig.spectrum("displayname",'el','fignum',123434);

        % El_sig.signal = awgn(El_sig.signal,-3,'measured',pn_key);

        %%%%% Lowpass el. components %%%%%%
        El_sig = Filter('filtdegree',2,"f_cutoff",60e9,"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);
        %  El_sig = El_sig.setPower(6,"dBm");

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


        % MAIN SIGNAL
        %%%%% MODULATE E/O CONVERSION %%%%%%
        vbias_rel = 0.5;
        u_pi = 2.9;
        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",pn_key).process(El_sig);

        % Opt_sig.eye(fsym,7);
        %
        % figure(10)
        % hold on
        % scatter(El_sig.signal(1:100000)+vbias,(abs(Opt_sig.signal(1:100000)).^2)*1e3,0.1,'.','DisplayName','Modulator TF')
        % ylim([0 4]);
        % xlim([-u_pi/2, u_pi/2]+vbias);
        % xlabel('Input in V')
        % ylabel('abs(Output) in mW')

        Optfilter = Filter('filtdegree',6,"f_cutoff",fsym.*0.7,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true);
        Opt_sig = Optfilter.process(Opt_sig);

        % Opt_sig.spectrum("fignum",122,"displayname",['Tx SPectrum; PAM ',num2str(M)]);

        Opt_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",0).process(Opt_sig);

        i_ = wh.parameter.rop.length;

        ber_ffe=zeros(i_);

        patten=zeros(i_);

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

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

            Scpe_sig.spectrum("displayname",'After Scope','fignum',123434);

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

            if db_channelapproach
                % ref symbols and transm. sequence are precoded
                [EQ_sig, Noi] = Eq.process(Scpe_sig,Duobinary().encode(Symbols));
            else
                [EQ_sig, Noi] = Eq.process(Scpe_sig,Symbols);
            end

            if db_encode || db_channelapproach
                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

            if postfilter
                % Noi.spectrum("displayname",'Noise Spectrum','fignum',1234);
                % EQ_sig.spectrum("displayname","Signal Spectrum","fignum",1234);
                
                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 0
                    [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

            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

        % wh.save('C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\MPI_August\auswertung\')

    end
end


disp('Simulation Done!')


cols = linspecer(8);

%cnt = cnt+1;
ber_ffe        = wh.getStoValue('ber_ffe',M,datarate,wh.parameter.rop.values);

% Create the initial plot

figure(44);
a = gca;
hold on; % Retain the plot so new points can be added without complete redraw

plot(wh.parameter.rop.values,ber_ffe,"LineWidth",0.5,"LineStyle","-","Marker",".","MarkerSize",15,"DisplayName","FFE only");
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