clear
Ref = load("C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\MPI_Juni\EQ_optimization\ref_sig_with_mpi.mat",'Ref');
Rx_sig = load("C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\MPI_Juni\EQ_optimization\rx_sig_with_mpi.mat",'Rx_sig');
Bits = load("C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\MPI_Juni\EQ_optimization\ref_bits_prms.mat","Bits");

M = 4;
Ref = Ref.Ref;
Rx_sig = Rx_sig.Rx_sig;
Bits = Bits.Bits;

kover = 8;
fdac = 256e9;
fadc = 160e9;

Lp_phd      = Filter('filtdegree',2,"f_cutoff",70e9,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true);

Lp_scpe = Filter('filtdegree',4,"f_cutoff",63e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true);

Scp     = 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",16,"quantbuffer",0.1,'block_dc',1,'lpf_active',1,'H_lpf',Lp_scpe);

fadc = Scp.fadc;

% SETUP TX Model
Pmap    = PAMmapper(M,0);

% SETUP RX Model
Phd     = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20);

% Eq      = EQ_silas_ofc("Ne",[50,2,2],"Nb",[2,0,0],"trainlength",4096,...
%     "sps",2,...
%     "mu_dc_dd",[0.5, 0.5, 0.5, 0.5],...
%     "mu_dc_train",0.0,...
%     "mu_ffe_train",0.00,...
%     "mu_dfe_train",0.005,...
%     "mu_ffe_dd",[0.0004 0.0004 0.0004],...
%     "mu_dfe_dd",0.005,...
%     "ddloops",3,...
%     "trainloops",4,...
%     "eq_parallelization_blocklength",0, ...
%     "eq_updatelatency",1,...
%     "eq_avg_blocklength",0);


if 0
    Eq = EQ("Ne",[21,4,4],"Nb",[0,0,0],"training_length",4096,"training_loops",4,"dd_loops",4,"K",2,"DCmu",0,"DDmu",[0.0004 0.0005 0.0006 0.0003 ],"DFEmu",0.000,"FFEmu",0,"plotfinal",0,"ideal_dfe",0);
    
    
    % % RX Model
    Rx_bits = rx_model(Rx_sig,Ref,Phd,Lp_phd,Scp,Eq,Pmap);
    
    % % BER
    [~,errors_bm,ber,errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
    
    disp(['BER: ',sprintf('%.1E',ber),' - - ROP: ',num2str(Rx_sig.power),'dBm - - PAM-',num2str(M),' - - ']);
end

%%%%% Local Parameter %%%%%%%%
    fsym = Ref.fs; 
    
%%%%%% Square Law %%%%%% 
    Rx_sig = Phd.process(Rx_sig);
 
%%%%%% Lowpass PhDiode %%%%%%
    Rx_sig = Lp_phd.process(Rx_sig);
    % Rx_sig.spectrum("displayname","Received Signal after PhD","fignum",201); 

%%%%%% Scope %%%%%% 
    Scpe_sig = Scp.process(Rx_sig);
    % Scpe_sig.spectrum("displayname",'after scope','fignum',123);
    
%%%%%% Sample to 2x fsym %%%%%% 
    Scpe_sig =  Scpe_sig.resample("fs_in",Scp.fadc,"fs_out",2*fsym);
    % Scpe_sig.plot('fignum',12345,'displayname','bla')
    
%%%%%% Sync Rx signal with reference %%%%%% 
    [Scpe_sig,D,cuts] = Scpe_sig.tsynch("reference",Ref,"fs_ref",fsym);



Scpe_sig = Scpe_sig.normalize("mode","rms");


Eq     = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",20,"sps",2,"decide",1);
% Eq = FFE_DFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"ffe_mu_dd",1e-4,"dfe_mu_dd",5e-4,"ffe_mu_tr",0,"dfe_mu_tr",0,"ffe_order",21,"dfe_order",0,"sps",2,"decide",1);
% Eq = VNLE("epochs_tr",7,"epochs_dd",7,"len_tr",4096,"mu_dd",[0.0004 0.0005 0.0006],"mu_tr",0,"order",[21,5,5],"sps",2,"decide",1);

Eq_sig = Eq.process(Scpe_sig,Ref);

Rx_bits = Pmap.demap(Eq_sig);

% BER
[~,errors_bm,ber,errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);

disp(['BER: ',sprintf('%.1E',ber),' - - ROP: ',num2str(Rx_sig.power),'dBm - - PAM-',num2str(M),' - - ']);