classdef Electrical_Trace_BiDi < handle

    properties(Access=public)
        file_path
        fsym
        rolloff
        K_over
        plot
        test
        S_test
    end

    methods(Access=public)
        function obj = Electrical_Trace_BiDi(options)
            arguments(Input)

                options.file_path = 'C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\COM Test\Mellitzz\TA_6002_6003_FX_B6_C6_B7_C7_Terminated.s4p'
                options.fsym = 0;
                options.rolloff = 0;
                options.K_over = 1;
                options.plot = 0;
                options.test = 0;
                options.S_test = [0,1;1,0];
            end

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

 

        end

        function [b_1,b_2] = process(obj, a_1, a_2)

            % Rx Signal Calculation Using S-Paramters For a 2-Port Network
            % [B_1(f); B_2(f)] = [S_11(f), S_12(f); S_21(f), S_22(f)] * [A_1(f); A_2(f)]

            % Initialize Rx Time Domain Signals
            b_1 = a_2;
            b_2 = a_1;

            % Loading and Extract S-Paramters
            net = sparameters(obj.file_path);
            f_s = net.Frequencies(:);

            if ~obj.test
                S2  = net.Parameters(1:2, 1:2, :);
            else
                S2 = repmat(obj.S_test,1,1,size(f_s,1));
            end
            
            S11_s = squeeze(S2(1,1,:));
            S12_s = squeeze(S2(1,2,:));
            S21_s = squeeze(S2(2,1,:));
            S22_s = squeeze(S2(2,2,:));

            % Setup Time Domain Signals
            x1 = a_1.signal(:);
            x2 = a_2.signal(:);
            N  = length(x1);
            assert(length(x2)==N, 'a_1 and a_2 must have same length');

            % Extract Time Domain Signal Frequencies
            if obj.fsym == 0 && obj.rolloff == 0 && obj.K_over == 1
                Fs = a_1.fs;
            else
                Fs = (1+obj.rolloff)*obj.fsym;
            end
            
            % Calculate Frequency Domain Signals
            A1 = fft(x1);
            A2 = fft(x2);
            
            % Calculate Frequency Axis [-Fs/2,...,Fs/2]
            f_fft = (0:N-1).' * (Fs/N);
            f_signed = f_fft;
            idxNeg = f_signed > Fs/2;
            f_signed(idxNeg) = f_signed(idxNeg) - Fs;      % Now in (-Fs/2, Fs/2]
            
            % Absolute Value Used for Interpolation
            f_pos = abs(f_signed);
            
            % Interpolate S-Parameters onto f_pos
            S11_p = interp1(f_s, S11_s, f_pos, 'linear', 'extrap');
            S12_p = interp1(f_s, S12_s, f_pos, 'linear', 'extrap');
            S21_p = interp1(f_s, S21_s, f_pos, 'linear', 'extrap');
            S22_p = interp1(f_s, S22_s, f_pos, 'linear', 'extrap');
            
            % Apply Hermitian Symmetry: S(-f) = conj(S(f))
            S11 = S11_p; S12 = S12_p; S21 = S21_p; S22 = S22_p;
            S11(idxNeg) = conj(S11_p(idxNeg));
            S12(idxNeg) = conj(S12_p(idxNeg));
            S21(idxNeg) = conj(S21_p(idxNeg));
            S22(idxNeg) = conj(S22_p(idxNeg));
            
            % Compute Rx Frequency Domain Signals
            B1 = S11 .* A1 + S12 .* A2;
            B2 = S21 .* A1 + S22 .* A2;
            
            % Calculate Rx Time Domain Signals
            b_1.signal = ifft(B1, 'symmetric');
            b_2.signal = ifft(B2, 'symmetric');

            if obj.plot
                rfplot(net)
            end


        end
    end
end