classdef EQ
    %EQ Summary of this class goes here
    %   Detailed explanation goes here

    properties
        Ne %Number of feed forward coefficients (1st, 2nd and 3rd order)
        Nb %Number of decision feedback coefficients (1st, 2nd and 3rd order)
        K  %Number of samples per symbol
        delay %Delay of incoming signal
        training_length %Number of training symbols
        training_loops %Number of loops through sequence for training mode
        ideal_dfe %Error free DFE decisions

        DB_aim %Aim at duobinary output sequence

        M %Order of the PAM constellation (only relevant in case of DB aim)

        FFEmu %mu parameter for FFE part in training mode (0 means normalized LMS)
        DFEmu % mu parameter for DFE part in training mode
        dd_loops % Number of loops through sequence for DD mode
        DDmu % mu parameters for DD mode (individual value for each order)
        DCmu % mu parameter for the dc tap

        l1act %Activate/deactive l1 regularization
        rho %Parameter for speed of coeff shrinking
        epsilon %Reciprocal value of the magnitude of the coeff to converge to zero (1st,2nd,3rd order)
        thres %Theshold for neglecting coefficienties (1st,2nd,3rd order)
        static_act %Activate/deactive static coefficient reduction

        mode2nd %0: no reduction | 1: polynomial | 2: restricted to interval
        len_2nd %length of the interval (only for 2nd order mode = 2)

        mode3rd %0: no reduction | 1: polynomial | 2: restricted to interval
        len_3rd %length of the interval (only for 3rd order mode = 3/4)

        plottrain
        plotfinal
        load_decisions

        save_taps


        %during simulation
        k0
        b
        b2
        b3
        e
        e2
        e3

        coeff_number
        constellation_in

    end

    methods
        function obj = EQ(options)
            %EQ Construct an instance of this class
            %   Detailed explanation goes here
            arguments(Input)
                options.Ne = [10 0 0] %Number of feed forward coefficients (1st, 2nd and 3rd order)
                options.Nb = [10 0 0]%Number of decision feedback coefficients (1st, 2nd and 3rd order)
                options.K  = 1 %Number of samples per symbol
                options.delay = 0 %Delay of incoming signal
                options.training_length = 1024 %Number of training symbols
                options.training_loops = 1 %Number of loops through sequence for training mode
                options.ideal_dfe = 0 %Error free DFE decisions

                options.DB_aim %Aim at duobinary output sequence

                options.M = 1 %Order of the PAM constellation (only relevant in case of DB aim)

                options.FFEmu = 0 %mu parameter for FFE part in training mode (0 means normalized LMS)
                options.DFEmu = 0.005 % mu parameter for DFE part in training mode
                options.dd_loops = 1% Number of loops through sequence for DD mode
                options.DDmu = [0.0004 0.0004 0.0004 0.0004 ] % mu parameters for DD mode (individual value for each order)
                options.DCmu = 0.005 % mu parameter for the dc tap

                options.l1act = 0 %Activate/deactive l1 regularization
                options.rho = 5e-4%Parameter for speed of coeff shrinking
                options.epsilon = [10 100 1000] %Reciprocal value of the magnitude of the coeff to converge to zero (1st,2nd,3rd order)
                options.thres = [5e-3 4e-3 5e-4]%Theshold for neglecting coefficienties (1st,2nd,3rd order)
                options.static_act = 0 %Activate/deactive static coefficient reduction

                options.mode2nd = 1%0: no reduction | 1: polynomial | 2: restricted to interval
                options.len_2nd = 1 %length of the interval (only for 2nd order mode = 2)

                options.mode3rd = 1%0: no reduction | 1: polynomial | 2: restricted to interval
                options.len_3rd = 1%length of the interval (only for 3rd order mode = 3/4)

                options.plottrain = 0
                options.plotfinal = 0
                options.load_decisions = 0
                options.save_taps = 0
            end

            fn = fieldnames(options);
            for n = 1:numel(fn)
                obj.(fn{n}) = options.(fn{n});
            end

        end

        function signalclass_out = process(obj,signalclass_in, reference_signalclass_in)

             % actual processing of the signal (steps 1. - 3.)
            signalclass_in.signal = obj.process_(signalclass_in.signal', reference_signalclass_in.signal'); 

            signalclass_in.signal = signalclass_in.signal';
             % append to logbook
             lbdesc = ['EQ '];
            signalclass_in  = signalclass_in.logbookentry(lbdesc);

             % write to output
            signalclass_out = signalclass_in;

        end

        function yout = process_(obj,data_in,ref_in)
            %METHOD1 Summary of this method goes here
            %   Detailed explanation goes here

            if obj.DB_aim
                ref_DB = zeros(size(ref_in));
                for k = 1:length(ref_in)
                    if k == 1
                        ref_DB(k) = ref_in(k);
                    else
                        ref_DB(k) = ref_in(k) + ref_in(k-1);
                    end
                end
                ref_in = ref_DB;
            end

            ref = [zeros(1,obj.Nb(1)-1) ref_in zeros(1,obj.Nb(1))];

            if isreal(ref)
                cplx = 0;
            else
                cplx = 1;
            end

            if obj.static_act
                obj.mode2nd = obj.mode2nd + 1;
                obj.mode3rd = obj.mode3rd + 1;
            else
                obj.mode2nd = 1;
                obj.mode3rd = 1;
            end

            if obj.mode2nd == 1
                N2 = (obj.Ne(2)*(obj.Ne(2)+1))/2; % Number of coefficients for second order
            elseif obj.mode2nd == 2
                N2 = obj.Ne(2);
            elseif obj.mode2nd == 3
                N2 = (obj.len_2nd+1)*(2*obj.Ne(2)-obj.len_2nd)/2;
            elseif obj.mode2nd == 4
                N2 = (ceil(obj.Ne(2)/2)+1)*(2*obj.Ne(2)-ceil(obj.Ne(2)/2))/2;
            end

            if obj.mode3rd == 1
                if cplx
                    N3 = obj.Ne(3)^2*(obj.Ne(3)+1)/2;
                else
                    N3 = obj.Ne(3)*(obj.Ne(3)+1)*(obj.Ne(3)+2)/6; % Number of coefficients for third order
                end
            elseif obj.mode3rd == 2
                N3 = obj.Ne(3);
            elseif obj.mode3rd == 3
                N3 = obj.Ne(3)^2;
            elseif obj.mode3rd == 4
                N3 = round(1/6*(obj.len_3rd+1)*(obj.len_3rd+2)*(3*obj.Ne(3)-2*obj.len_3rd));
            elseif obj.mode3rd == 5
                N3 = 2*obj.Ne(3)*obj.len_3rd-obj.len_3rd*(obj.len_3rd+1)+obj.Ne(3);
            end

            Nb2 = (obj.Nb(2)*(obj.Nb(2)+1))/2;
            Nb3 = obj.Nb(3)*(obj.Nb(3)+1)*(obj.Nb(3)+2)/6;

            data_in = data_in/sqrt(mean(abs(data_in).^2));  % power normalization of input sequence

            if obj.FFEmu == 0
                norm_fac2 = sqrt(mean(abs(data_in.^2).^2));     % power normalization for second and third order terms
                norm_fac3 = sqrt(mean(abs(data_in.^3).^2));     % (not necessary, but seems to be more stable if applied --> same as different mu values for linear and nl terms)
            else
                norm_fac2 = 1;
                norm_fac3 = 1;
            end

            norm_fac_DFE2 = sqrt(mean(abs(ref_in.^2).^2)); % same for DFE input (reference)
            norm_fac_DFE3 = sqrt(mean(abs(ref_in.^3).^2));

            data = [zeros(1,floor(obj.Ne(1)/2)) data_in zeros(1,obj.Ne(1))];

            delta_2 = round((obj.Ne(1)-obj.Ne(2))/2);
            delta_3 = round((obj.Ne(1)-obj.Ne(3))/2);

            delta_DFE2 = 1;%round((obj.Nb-obj.Nb(2))/2);
            delta_DFE3 = 1;%round((obj.Nb-obj.Nb(3))/2);

            % calculate the indices for the combination of second and third order symbols
            % - done in advance because it's the same for each iteration, so time
            % can be saved

            [ind_mat_2nd,ind_mat_3rd] = obj.calc_ind(obj.Ne(2),N2,obj.Ne(3),N3,obj.mode2nd,obj.mode3rd,obj.len_2nd,obj.len_3rd,cplx);

            [ind_mat_DFE_2nd,ind_mat_DFE_3rd] = obj.calc_DFE_ind(obj.Nb(2),Nb2,obj.Nb(3),Nb3);

            if obj.l1act
                epsilon_ = diag([ones(1,obj.Ne(1))*obj.epsilon(1) ones(1,N2)*obj.epsilon(2) ones(1,N3)*obj.epsilon(3)]);
            end

            obj.k0 = obj.delay; % input delay compared to training sequence<

            if 1 %  obj.active
                %% Calculation of the filter coefficients in training based LMS mode
                e_ = zeros(obj.Ne(1)+N2+N3,1); % initialization of filter coefficients
                %             e(ceil(obj.Ne(1)/2)) = 1;     % set central tap to 1 (better starting point since it's closer to the expected solution)
                b_ = zeros(obj.Nb(1)+Nb2+Nb3,1);
                e_dc = mean(data_in);       % initilaization of the dc tap with the mean value of the data
                %             e_save = NaN(361,8.6e5);
                %             save_ind = 1;
                for trainloops = 1:obj.training_loops
                    m = obj.k0+1; % starting symbol index at the delay compared to the training sequence
                    for n = obj.K*obj.k0+1:obj.K:obj.K*obj.training_length
                        m = m+1;
                        X_1 = data(obj.Ne(1)+n+(obj.K-1):-1:n+obj.K).';
                        [X_2,X_3] = obj.calc_nl_vecs(X_1,ind_mat_2nd,ind_mat_3rd,norm_fac2,norm_fac3,delta_2,delta_3,cplx);

                        D_1 = ref(obj.Nb(1)-obj.k0+m-2:-1:m-obj.k0-1).';
                        [D_2,D_3] = obj.calc_nl_vecs(D_1,ind_mat_DFE_2nd,ind_mat_DFE_3rd,norm_fac_DFE2,norm_fac_DFE3,delta_DFE2,delta_DFE3,cplx);

                        input_vec = [X_1;X_2;X_3];
                        reference_vec = [D_1;D_2;D_3];
                        error = e_dc + e_.'*input_vec - b_.'*reference_vec  - ref_in(m-obj.k0); % error = e_dc + e.'*input_vec - b.'*reference_vec  - ref_in(m-obj.k0);

                        if real(obj.FFEmu)
                            if obj.l1act
                                sgn_e = e_;
                                sgn_e(e_~=0) = e_(e_~=0)./abs(e_(e_~=0));
                                e_ = e_ - obj.rho*sgn_e./(1+epsilon_*abs(e_)) - error*input_vec*obj.FFEmu;
                            else
                                e_ = e_ - error*conj(input_vec)*obj.FFEmu; %e = e - error*conj(input_vec)*obj.FFEmu; %e = e - error*input_vec*obj.FFEmu;
                            end
                        else
                            if obj.l1act
                                sgn_e = e_;
                                sgn_e(e_~=0) = e_(e_~=0)./abs(e_(e_~=0));
                                e_ = e_ - obj.rho*sgn_e./(1+epsilon_*abs(e_)) - error*input_vec/(input_vec.'*input_vec);
                            else
                                e_ = e_ - error*input_vec/(input_vec.'*input_vec);
                            end
                        end
                        %                       e_save(:,save_ind) = e;
                        %                       save_ind = save_ind+1;
                        e_dc = e_dc - obj.DCmu*error;

                        if obj.Nb(1) > 0
                            b_ = b_ + obj.DFEmu*error*reference_vec; % Seems like normalized DFE has worse performance
                        end
                    end
                end
                %%

                % Plot the intermediate coefficients after training mode
                obj.b = b_(1:obj.Nb(1));
                obj.b2 = b_(obj.Nb(1)+1:obj.Nb(1)+Nb2);
                obj.b3 = b_(obj.Nb(1)+Nb2+1:end);
                obj.e = e_(1:obj.Ne(1));
                obj.e2 = e_(obj.Ne(1)+1:obj.Ne(1)+N2);
                obj.e3 = e_(obj.Ne(1)+N2+1:end);

                if obj.plottrain
                    figure(8052)
                    subplot(2,3,1); stem(abs(obj.e),'Markersize',2);
                    title('FFE coeff linear')
                    xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12)
                    subplot(2,3,2); stem(obj.e2,'Markersize',2);
                    title('FFE coeff nl 2nd')
                    xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12)
                    subplot(2,3,3); stem(obj.e3,'Markersize',2);
                    title('FFE coeff nl 3rd')
                    xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12)
                    subplot(2,3,4);stem(obj.b,'Markersize',2);
                    title('DFE coeff linear')
                    xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12)
                    subplot(2,3,5);stem(obj.b2,'Markersize',2);
                    title('DFE coeff nl 2nd')
                    xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12)
                    subplot(2,3,6);stem(obj.b3,'Markersize',2);
                    title('DFE coeff nl 3rd')
                    xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12)
                    set(gcf,'Position',[200 500 700 400])
                end

                if obj.l1act
                    neg_lin = find(abs(obj.e) < obj.thres(1));
                    neg_2nd = find(abs(obj.e2) < obj.thres(2));
                    neg_3rd = find(abs(obj.e3) < obj.thres(3));

                    neg = [neg_lin;neg_2nd+obj.Ne(1);neg_3rd+obj.Ne(1)+N2]; % indices of the neglected coefficients

                    rel_lin = find(abs(obj.e) >= obj.thres(1));
                    rel_2nd = find(abs(obj.e2) >= obj.thres(2));
                    rel_3rd = find(abs(obj.e3) >= obj.thres(3));

                    obj.coeff_number = length(rel_lin)+2*length(rel_2nd)+3*length(rel_3rd);

                    rel = [rel_lin;rel_2nd+obj.Ne(1);rel_3rd+obj.Ne(1)+N2]; % indices of the relevant coefficients
                    e_(neg) = 0;

                    ind_mat_2nd(neg_2nd,:) = [];

                    ind_mat_3rd(neg_3rd,:) = [];
                end

                %% decision directed mode
                if ~obj.DB_aim
                    constellation_in_ = unique(ref_in); % getting the symbol constellation from reference data
                else
                    if obj.M == 2
                        constellation_in_ = [-3 -2 -1 0 1 2 3]/sqrt(5)*2;
                    elseif obj.M == 2.5
                        constellation_in_ = [-5 -4 -3 -2 -1 0 1 2 3 4 5]/sqrt(10)*2;
                    elseif obj.M == 3
                        constellation_in_ = [-7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7]/sqrt(21)*2;
                    else
                        constellation_in_ = unique(ref_in);
                    end
                end

                obj.constellation_in = constellation_in_;

                if obj.l1act
                    coeff = [e_(rel);b_];  % combine FFE and DFE coefficient vectors for DD mode
                else
                    coeff = [e_;b_];
                end

                for dd_loop = 1:obj.dd_loops
                    m = 0;
                    output_vec = zeros(1,floor(length(data_in)/obj.K)); % initilaization of the output vector
                    dd_DFE = zeros(obj.Nb(1),1);
                    D_2 = zeros(Nb2,1);
                    D_3 = zeros(Nb3,1);

                    if all(obj.DDmu == obj.DDmu(1))
                        mu_mat = obj.DDmu(1);
                    else
                        if obj.l1act
                            mu_mat = diag([ones(1,length(rel_lin))*obj.DDmu(1) ones(1,length(rel_2nd))*obj.DDmu(2) ones(1,length(rel_3rd))*obj.DDmu(3) ones(1,obj.Nb)*obj.DDmu(4)]);
                        else
                            mu_mat = diag([ones(1,obj.Ne(1))*obj.DDmu(1) ones(1,N2)*obj.DDmu(2) ones(1,N3)*obj.DDmu(3) ones(1,obj.Nb(1)+Nb2+Nb3)*obj.DDmu(4)]);
                        end
                    end

                    if obj.load_decisions
                        pathn = evalin('base','modeldir');
                        temp = load([pathn, 'MLSE_out', '.mat']) ;
                        %eval(['dd_out_vals = temp.', 'a', ';']) ;
                        dd_out_vals=temp.a;
                        dd_out = zeros(size(data_in));
                        dd_out(1:2:length(data_in)) = dd_out_vals;
                    else
                        dd_out = zeros(size(data_in));
                    end

                    for k = 1:obj.K:length(data_in)
                        m=m+1; % Symbol index
                        X_1 = data(obj.Ne(1)+k-1:-1:k).';
                        [X_2,X_3] = obj.calc_nl_vecs(X_1,ind_mat_2nd,ind_mat_3rd,norm_fac2,norm_fac3,delta_2,delta_3,cplx);

                        if obj.l1act
                            input_vec = [X_1(rel_lin);X_2;X_3;-dd_DFE;-D_2;-D_3];
                        else
                            input_vec = [X_1;X_2;X_3;-dd_DFE;-D_2;-D_3];
                        end

                        output_vec(m) = e_dc + input_vec.'*coeff;

                        if ~obj.load_decisions
                            [~,dd_idx] = min(abs(output_vec(m) - constellation_in_)); % decision for closest constellation point
                            dd_out(k) = constellation_in_(dd_idx);
                        end

                        if obj.Nb(1) > 0
                            dd_DFE(2:end) = dd_DFE(1:end-1);
                            dd_DFE(1) = dd_out(k);

                            if obj.ideal_dfe && m > obj.k0
                                dd_DFE(1) = ref_in(m-obj.k0);
                            end
                            [D_2,D_3] = obj.calc_nl_vecs(dd_DFE,ind_mat_DFE_2nd,ind_mat_DFE_3rd,norm_fac_DFE2,norm_fac_DFE3,delta_DFE2,delta_DFE3,cplx);
                        end
                        %                     if dd_loop ~= 21
                        error = output_vec(m) - dd_out(k);
                        %                     else
                        %                         error = 0;
                        %                     end

                        coeff = coeff - mu_mat*error*conj(input_vec);
                        %                     e_save(:,save_ind) = coeff;
                        %                     save_ind = save_ind+1;

                        if mu_mat ~= 0
                            e_dc = e_dc - obj.DCmu*error;
                        end

                    end
                end


                % shifting the output sequence by k0 symbols
                yout = (circshift(output_vec.',-(obj.k0))).'; %(circshift(dd_out.',-(obj.k0))).';

                e_ = coeff(1:end-obj.Nb(1)-Nb2-Nb3);
                b_ = coeff(end-obj.Nb(1)-Nb2-Nb3+1:end);

                if obj.l1act
                    obj.e = e_(1:length(rel_lin));
                    obj.e2 = e_(length(rel_lin)+1:length(rel_lin)+length(rel_2nd));
                    obj.e3 = e_(length(rel_lin)+length(rel_2nd)+1:end);
                else
                    obj.e = e_(1:obj.Ne(1));
                    obj.e2 = e_(obj.Ne(1)+1:obj.Ne(1)+N2);
                    obj.e3 = e_(obj.Ne(1)+N2+1:end);
                end
                obj.b = b_(1:obj.Nb(1));
                obj.b2 = b_(obj.Nb(1)+1:obj.Nb(1)+Nb2);
                obj.b3 = b_(obj.Nb(1)+Nb2+1:end);

                % plot the final coefficients after DD mode
                if obj.plotfinal
                    figure(8054)
                    if obj.l1act
                        subplot(2,3,1); stem(rel_lin,obj.e,'Markersize',2);
                        title('FFE coeff linear')
                        xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12)
                        subplot(2,3,2); stem(rel_2nd,obj.e2,'Markersize',2);
                        title('FFE coeff nl 2nd')
                        xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12)
                        subplot(2,3,3); stem(rel_3rd,obj.e3,'Markersize',2);
                        title('FFE coeff nl 3rd')
                        xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12)
                        subplot(2,3,4);stem(obj.b,'Markersize',2);
                        title('DFE coeff linear')
                        xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12)
                        subplot(2,3,5);stem(obj.b2,'Markersize',2);
                        title('DFE coeff nl 2nd')
                        xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12)
                        subplot(2,3,6);stem(obj.b3,'Markersize',2);
                        title('DFE coeff nl 3rd')
                        xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12)
                    else
                        subplot(2,3,1); stem(obj.e/max(e_),'Markersize',2);
                        title('FFE coeff linear')
                        xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12)
                        subplot(2,3,2); stem(obj.e2,'Markersize',2);
                        title('FFE coeff nl 2nd')
                        xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12)
                        subplot(2,3,3); stem(obj.e3,'Markersize',2);
                        title('FFE coeff nl 3rd')
                        xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12)
                        subplot(2,3,4);stem(obj.b,'Markersize',2);
                        title('DFE coeff linear')
                        xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12)
                        subplot(2,3,5);stem(obj.b2,'Markersize',2);
                        title('DFE coeff nl 2nd')
                        xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12)
                        subplot(2,3,6);stem(obj.b3,'Markersize',2);
                        title('DFE coeff nl 3rd')
                        xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12)
                    end
                    set(gcf,'Position',[1000 500 700 400])
                end

                % save frequency response to the work space
                if obj.save_taps
                    % save the FFE coefficients to the work space
                    % pathn = evalin('base','modeldir');
                    % eval([obj.field_ffe, ' = obj.e ;']) ;
                    % eval([obj.field_dfe, ' = b ;']) ;
                    % eval(['save(''', pathn, '\',obj.filen,''', ''', obj.field_ffe,''', ''',obj.field_dfe,''') ;']) ;

                    save("coefficients",obj.e, obj.b);
                end

            else
                yout = data_in;
            end

        end


        function [X_2,X_3] = calc_nl_vecs(obj,X_1,ind_mat_2,ind_mat_3,norm_fac2,norm_fac3,delta_2,delta_3,cplx)
            % calculation of the vectors containing all combinations of input symbols
            % of second and third order based on the linear symbols

            if ind_mat_2(1) > 0
                input_vec_se = X_1(delta_2:end)/norm_fac2;%(K*(k0-1):end)
                X_2 = input_vec_se(ind_mat_2(:,1)).*input_vec_se(ind_mat_2(:,2));
            else
                X_2 = [];
            end

            if ind_mat_3(1) > 0
                if cplx
                    input_vec_th = X_1(delta_3:end)/norm_fac3;
                    X_3 = input_vec_th(ind_mat_3(:,1)).*input_vec_th(ind_mat_3(:,2)).*conj(input_vec_th(ind_mat_3(:,3)));
                else
                    input_vec_th = X_1(delta_3:end)/norm_fac3;
                    X_3 = input_vec_th(ind_mat_3(:,1)).*input_vec_th(ind_mat_3(:,2)).*input_vec_th(ind_mat_3(:,3));
                end
            else
                X_3 = [];
            end
        end

        function [ind_mat_2nd,ind_mat_3rd] = calc_ind(obj,Ne2,N2,Ne3,N3,mode2nd,mode3rd,len_2nd,len_3rd,cplx)

            if Ne2 > 0
                ind_mat_2nd = NaN(N2,2);
                count=1;
                if mode2nd == 1
                    for t = 1:Ne2
                        for u = t:Ne2
                            ind_mat_2nd(count,:) = [t u];
                            count = count + 1 ;
                        end
                    end
                elseif mode2nd == 2
                    for t = 1:Ne2
                        ind_mat_2nd(t,:) = [t t];
                    end
                elseif mode2nd == 3
                    for t = 1:Ne2
                        for u = t:Ne2
                            if u-t<=len_2nd
                                ind_mat_2nd(count,:) = [t u];
                                count = count + 1 ;
                            end
                        end
                    end
                elseif mode2nd == 4
                    for t = 1:Ne2
                        for u = t:Ne2
                            if u-t<=ceil(Ne2/2)
                                ind_mat_2nd(count,:) = [t u];
                                count = count + 1 ;
                            end
                        end
                    end
                end

            else
                ind_mat_2nd = 0;
            end

            if Ne3 > 0
                ind_mat_3rd = NaN(N3,3);
                count=1;
                if mode3rd == 1
                    if cplx
                        for t = 1:Ne3
                            for u = t:Ne3
                                for v = 1:Ne3
                                    ind_mat_3rd(count,:) = [t u v];
                                    count = count + 1 ;
                                end
                            end
                        end
                    else
                        for t = 1:Ne3
                            for u = t:Ne3
                                for v = u:Ne3
                                    ind_mat_3rd(count,:) = [t u v];
                                    count = count + 1 ;
                                end
                            end
                        end
                    end
                elseif mode3rd == 2
                    for t = 1:Ne3
                        ind_mat_3rd(t,:) = [t t t];
                    end
                elseif mode3rd == 3
                    for t = 1:Ne3
                        for u = t:Ne3
                            ind_mat_3rd(count,:) = [t t u];
                            if t ~= u
                                count = count + 1;
                                ind_mat_3rd(count,:) = [t u u];
                            end
                            count = count + 1;
                        end
                    end
                elseif mode3rd == 4
                    for t = 1:Ne3
                        for u = t:Ne3
                            for v = u:Ne3
                                if u-t<=len_3rd && v-t<=len_3rd
                                    ind_mat_3rd(count,:) = [t u v];
                                    count = count + 1 ;
                                end
                            end
                        end
                    end
                elseif mode3rd == 5
                    for t = 1:Ne3
                        for u = t:Ne3
                            if u-t<=len_3rd
                                ind_mat_3rd(count,:) = [t t u];
                                if t ~= u
                                    count = count + 1;
                                    ind_mat_3rd(count,:) = [t u u];
                                end
                                count = count + 1;
                            end
                        end
                    end

                    ind_mat_3rd2 = NaN(N3,3);
                    count = 1;
                    %         for t = 1:Ne3
                    %             ind_mat_3rd2(count,:) = [t t t];
                    %             count = count + 1;
                    %         end
                    for t = 1:Ne3
                        %             ind_mat_3rd2(count,:) = [t t t];
                        %             count = count + 1;
                        for u = t:min(Ne3,t+len_3rd)
                            for v = unique([t u])
                                ind_mat_3rd2(count,:) = [t v u];
                                count = count + 1;
                                %                     ind_mat_3rd2(count,:) = [t u u];
                                %                     count = count + 1;
                            end
                        end
                    end
                end
            else
                ind_mat_3rd = 0;
            end

        end

        function [ind_mat_2nd,ind_mat_3rd] = calc_DFE_ind(obj,Ne2,N2,Ne3,N3)

            if Ne2 > 0
                ind_mat_2nd = NaN(N2,2);
                count=1;
                for t = 1:Ne2
                    for u = t:Ne2
                        ind_mat_2nd(count,:) = [t u];
                        count = count + 1 ;
                    end
                end
            else
                ind_mat_2nd = 0;
            end

            if Ne3 > 0
                ind_mat_3rd = NaN(N3,3);
                count=1;
                for t = 1:Ne3
                    for u = t:Ne3
                        for v = u:Ne3
                            ind_mat_3rd(count,:) = [t u v];
                            count = count + 1 ;
                        end
                    end
                end
            else
                ind_mat_3rd = 0;
            end

        end

    end
end