function  [output] = imdd_model(simulation_mode,database,varargin)




%%% Change folder
curFolder = pwd;
funcFolder=fileparts(mfilename('fullpath'));
if ~isempty(funcFolder)
    cd(funcFolder);
end

%%% Run parameters
% TX
M               = 4;
fsym            = 180e9;
f_nyquist       = fsym/2;
apply_pulsef    = 1;
fdac            = 256e9;
fadc            = 256e9;
random_key      = 1;

db_precode      = 1;
emulate_precode = 0;
discard_precode = 1;
db_encode       = 0;
% duob_mode       = db_mode.db_emulate;
emulate_db = 1;

rcalpha         = 0.05;
kover           = 16;
vbias_rel       = 0.5;
u_pi            = 2.9;
vbias           = -vbias_rel*u_pi;
laser_wavelength    = 1310;
laser_linewidth     = 0;
tx_bw_nyquist   = 0.9;

% Channel
link_length     = 1;

% RX
rop = -8;
rx_bw_nyquist = 0.7;

% EQ
eq_mode = equalizer_structure.vnle_pf_mlse;
ffe_order=[50,0,0];
vnle_order=[50,7,7];
dfe_order = [2 0 0];

len_tr = 4096*2;
mu_ffe = [0.0004 0.0004 0.0004];
mu_dfe = 0.0004;
mu_dc = 0.00;

dfe_ = sum(dfe_order)>0;

%%% change specific parameter if given in varargin
% Parse optional input arguments
if ~isempty(varargin)
    var_s = varargin{1};
    if isstruct(var_s)
        fields = fieldnames(var_s);
        for i = 1:numel(fields)
            if isnumeric(fields{i})
                eval([fields{i}, ' = ',  num2str( var_s.(fields{i}) ), ';']);
                fprintf("%s <-- %.2f \n", fields{i}, var_s.(fields{i}));
            else
                eval([fields{i}, ' = ',  'var_s.(fields{',num2str(i),'})' , ';']);
            end

        end
    else
        error('Optional variables should be passed as a struct.');
    end
end

%%% run the simulation or measurement or ...


if simulation_mode

    fsym_ = floor( bitrate*1e-9./log2(M) ).*1e9;

    if fsym_ ~= fsym
        fsym = fsym_;
        fprintf('Adapted symbolrate to %d GBd, to match provided bitrate of %d GBit/s using PAM %d \n',fsym.*1e-9,bitrate.*1e-9, M);
    end

    tx_simulation;

    

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

    rx_simulation;

    Scpe_cell{1} = Scpe_sig;

else

    basePath = 'C:\Users\Silas\Documents\MATLAB\Datensätze\sioe_labor\';
    useGui = 0;
    % db = DBHandler("pathToDB",[basePath,'silas_labor.db']);
    filterParams = database.tables;
    % filterParams.Runs.run_id = 2958; % no db
    % filterParams.Runs.run_id = 2937; % no db
    filterParams.Configurations = struct( ...
        'bitrate',          bitrate, ...
        'db_mode',          db_precode+db_encode, ...
        'fiber_length',     link_length, ...
        'interference_attenuation', [], ...
        'interference_path_length', [], ...
        'is_mpi',           0, ...
        'pam_level',        M, ...
        'precomp_amp',      [], ...
        'rop_attenuation',  0, ...
        'symbolrate',       [], ...
        'v_awg',            [], ...
        'v_bias',           [], ...
        'wavelength',       1310 ...
        );

    selectedFields = {'Runs.run_id','Runs.tx_bits_path', 'Runs.tx_symbols_path', 'Runs.rx_sync_path','Runs.rx_raw_path',...
        'Configurations.db_mode','Configurations.pam_level','Configurations.bitrate','Configurations.symbolrate','Configurations.fiber_length','Configurations.wavelength','Configurations.precomp_amp','BERs.ber'};

    [dataTable,sql_query] = database.queryDB(filterParams, selectedFields);
    [~, uniqueIdx] = unique(dataTable.run_id);  % Get unique run_id indices
    dataTable = dataTable(uniqueIdx,:);  % Extract unique configurations for each run_id
    fprintf('Found %d entries for requested Configuration. IDs are: %s \n \n',size(dataTable,1),jsonencode(dataTable.run_id(1:min(size(dataTable,1),100))));

    Tx_bits     = load([basePath, char(dataTable.tx_bits_path(2))]);
    Tx_bits     = Tx_bits.Bits;

    Symbols  = load([basePath, char(dataTable.tx_symbols_path(2))]);
    Symbols  = Symbols.Symbols;

    Scpe_load     = load([basePath, char(dataTable.rx_sync_path(2))]);
    Scpe_cell     = Scpe_load.S;

    % Raw_signal = load([basePath, char(dataTable.rx_raw_path(1))]);
    % Raw_signal = Raw_signal.Scpe_sig_raw;

    fsym = Symbols.fs;

end

output = struct();

for occ = 1:1%length(Scpe_cell)

    Scpe_sig = Scpe_cell{occ};

    %%%%%% Sample to 2x fsym %%%%%%
    Scpe_sig =  Scpe_sig.resample("fs_out",2*fsym);

    %%%%%% Sync Rx signal with reference %%%%%%
    [Scpe_sig,S] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym);

    % Scpe_sig.spectrum("displayname",'Rx (Scpe+Sync+Resample)','fignum',100,'normalizeTo0dB',0);

    Scpe_sig = Filter('filtdegree',4,"f_cutoff",Symbols.fs.*0.6,"fs",Scpe_sig.fs,"filterType",filtertypes.gaussian,"active",true).process(Scpe_sig);

    Scpe_sig.spectrum("displayname",'Rx (Scpe+Sync+Resample+LPF)','fignum',110,'normalizeTo0dB',0);

    Scpe_sig.plot("displayname",'Filtered Scope Signal','fignum',111,'clear',1);

    Scpe_sig.eye(fsym,M);



    %%% EQUALIZING

    switch eq_mode

        case equalizer_structure.ffe

            %FFE
            if db_precode
                Bits_ = PAMmapper(M,0).demap(Symbols);
            else
                Bits_ = Bits;
            end

            eq_ffe = EQ("Ne",ffe_order,"Nb",dfe_order,"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1);
            % eq_ffe     = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",len_tr,"mu_dd",mu_ffe(1),"mu_tr",0,"order",ffe_order(1),"sps",2,"decide",0);
            % eq_ffe     = FFE_DCremoval("epochs_tr",5,"epochs_dd",5,"len_tr",len_tr,"mu_dd",mu_ffe(1),"mu_tr",0,"order",ffe_order(1),"sps",2,"decide",0,"dc_buffer_len",512,"mu_dc",0.05);

            [eq_sig,Eq_noise,ber_ffe(occ),totalErrors] = vnle( eq_ffe,M,Scpe_sig ,Symbols, Bits_);
            Eq_noise.spectrum("displayname",'Noise Spectrum after FFE','fignum',41,'normalizeTo0dB',0);

            fprintf('BER FFE: %.2e \n',ber_ffe(occ));

        case equalizer_structure.vnle





            %VNLE
            eq_vnle = EQ("Ne",vnle_order,"Nb",dfe_order,"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",0);

            [Eq_sig,Eq_noise] = eq_vnle.process(Scpe_sig,Symbols);

            Eq_sig = PAMmapper(M,0).quantize(Eq_sig);

            if emulate_precode && db_precode == 0
                % emulation
                Eq_sig = Duobinary().encode(Eq_sig);
                Eq_sig = Duobinary().decode(Eq_sig);

                Symbols= Duobinary().encode(Symbols);
                Symbols = Duobinary().decode(Symbols);
                Tx_bits = PAMmapper(M,0).demap(Symbols);
            elseif db_precode == 1 && db_encode == 0 && discard_precode == 1
                % normal dsp for precoded sequence
                Tx_bits = PAMmapper(M,0).demap(Symbols);

            elseif db_precode == 1 && db_encode == 1
                error('not implemented')

            elseif db_precode == 1 && db_encode == 0
                Eq_sig = Duobinary().encode(Eq_sig);
                Eq_sig = Duobinary().decode(Eq_sig);
            end

            % M = numel(unique(tx_symbols.signal));
            Rx_bits = PAMmapper(M,0).demap(Eq_sig);

            [~,numErrors,ber,~] = calc_ber(Rx_bits.signal,Tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);

            Eq_noise.spectrum("displayname",'Noise Spectrum after VNLE','fignum',41,'normalizeTo0dB',0);

            output.ber_vnle(occ) =  ber;

            if 1c
                figure(51);
                clf
                title(sprintf('DB coded PAM after EQ ; BER: %1.2e',M, ber ));
                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

            Eq_sig.eye(fsym,M,"displayname",'Eye after EQ','fignum',52);
            Eq_sig.spectrum("displayname",'Spectrum after EQ','fignum',53,'normalizeTo0dB',1);
            Scpe_sig.spectrum("displayname",'Spectrum before EQ','fignum',53,'normalizeTo0dB',1);
            fprintf('BER FFE: %.2e \n',ber);





        case equalizer_structure.vnle_pf_mlse


            %VNLE + PF + MLSE
            eq_mlse     = EQ("Ne",vnle_order,"Nb",dfe_order,"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1);
            % eq_mlse     = FFE_DCremoval("epochs_tr",5,"epochs_dd",5,"len_tr",len_tr,"mu_dd",mu_ffe(1),"mu_tr",0,"order",ffe_order(1),"sps",2,"decide",0,"dc_buffer_len",1,"mu_dc",0.05);
            % eq_mlse     = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",len_tr,"mu_dd",mu_ffe(1),"mu_tr",0,"order",ffe_order(1),"sps",2,"decide",0);
            % eq_mlse     = FFE_DCremoval("epochs_tr",5,"epochs_dd",5,"len_tr",len_tr,"mu_dd",mu_ffe(1),"mu_tr",0,"order",ffe_order(1),"sps",2,"decide",0,"dc_buffer_len",512,"mu_dc",0.05);

            pf_         = Postfilter("ncoeff",2);
            mlse_       = MLSE("DIR",[0,0],"duobinary_output",0,"M",[],"trellis_states",[]);

            %FFE or VNLE
            [Eq_sig,Eq_noise] = eq_mlse.process(Scpe_sig,Symbols);

            Eq_sig = pf_.process(Eq_sig,Eq_noise);

            %M = numel(unique(tx_symbols.signal));
            mlse_.DIR = pf_.burg_coeff;
            mlse_.trellis_states = PAMmapper(M,0).levels;
            mlse_.M = M;

            Eq_sig = mlse_.process(Eq_sig);
            pf_.showFilter(Eq_noise);
            % Eq_noise.spectrum("displayname",'Noise Spectrum after VNLE+PF','fignum',41,'normalizeTo0dB',0);

            if emulate_precode && db_precode == 0
                % emulation
                Eq_sig = Duobinary().encode(Eq_sig);
                Eq_sig = Duobinary().decode(Eq_sig);

                Symbols= Duobinary().encode(Symbols);
                Symbols = Duobinary().decode(Symbols);
                Tx_bits = PAMmapper(M,0).demap(Symbols);
            elseif db_precode == 1 && db_encode == 0 && discard_precode == 1
                % normal dsp for precoded sequence
                Tx_bits = PAMmapper(M,0).demap(Symbols);

            elseif db_precode == 1 && db_encode == 1
                error('not implemented')

            elseif db_precode == 1 && db_encode == 0
                Eq_sig = Duobinary().encode(Eq_sig);
                Eq_sig = Duobinary().decode(Eq_sig);
            end

            Rx_bits = PAMmapper(M,0).demap(Eq_sig);

            [~,numErrors,ber,~] = calc_ber(Rx_bits.signal,Tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);

            output.ber_mlse(occ) =  ber;
            output.pf_taps(occ,:) = pf_.burg_coeff;


            fprintf('BER MLSE: %.2e \n',ber);

        case equalizer_structure.db_precoded

            %EQ targets DB => less precompensation; pre-coded
            mlse_db_pre = MLSE("DIR",[1,1],"duobinary_output",1,"M",M,"trellis_states",PAMmapper(M,0).levels);
            eq_db_pre = EQ("Ne",vnle_order,"Nb",dfe_order,"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1);
            [eq_sig,Eq_noise,ber_db(occ),totalErrors] = duobinary_target(eq_db_pre, mlse_db_pre,M, Scpe_sig ,Symbols, Bits);
            Eq_noise.spectrum("displayname",'Noise Spectrum after DB','fignum',41,'normalizeTo0dB',0);
            %->append BER to DB

            fprintf('BER VNLE+DB: %.2e \n',ber_db(occ));

        case equalizer_structure.db_encoded

            %db signaling => db encoded
            mlse_db_enc = MLSE("DIR",[1,1],"duobinary_output",1,"M",M,"trellis_states",PAMmapper(M,0).levels);
            eq_db_enc = EQ("Ne",vnle_order,"Nb",dfe_order,"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1);
            [eq_sig,Eq_noise,ber_db_enc,totalErrors] = duobinary_signaling(eq_db_enc, mlse_db_enc,M, Scpe_sig ,Symbols, Bits);
            %->append BER to DB

            fprintf('BER DB: %.2e \n',ber_db_enc(occ));
    end
    autoArrangeFigures;


end



if ~isempty(curFolder)
    cd(curFolder);
end

end