classdef Electrical_Hybrid < handle

    properties(Access=public)
        file_path
        plot = 0

        % If true: perform digital residual-echo cancellation using known TX
        cancel_echo = 1

        % If true: in addition return v_hyb and v_echo_est
        return_intermediates = 1
    end

    methods(Access=public)

        function obj = Electrical_Hybrid(options)
            arguments
                options.file_path = ''
                options.plot = 0
                options.cancel_echo = 1
                options.return_intermediates = 1
            end

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

        function [v_fe_rec, v_hyb, v_echo_est] = process(obj, v_ne_tx, v_fe_tx)
            % v_ne_tx : Near-End TX signal object (Port 1)
            % v_fe_tx : Far-End TX signal object  (Port 2)
            %
            % Output:
            % v_hyb      : physical hybrid differential output (Port4 - Port3)
            % v_echo_est : estimated residual echo due to local TX only
            % v_fe_rec   : v_hyb - v_echo_est  (if cancel_echo enabled), else v_hyb

            % --- Load S-Parameters (.s4p) ---
            net = sparameters(obj.file_path);
            f_s = net.Frequencies(:);
            S4  = net.Parameters;

            % Extract needed S-parameters for differential output:
            % V3 = S31*V1 + S32*V2
            % V4 = S41*V1 + S42*V2
            % Vhyb = V4 - V3 = (S41-S31)*V1 + (S42-S32)*V2
            S31_s = squeeze(S4(3,1,:));
            S41_s = squeeze(S4(4,1,:));
            S32_s = squeeze(S4(3,2,:));
            S42_s = squeeze(S4(4,2,:));

            % --- Time-domain signals ---
            x1 = v_ne_tx.signal(:);
            x2 = v_fe_tx.signal(:);

            if length(x2) ~= length(x1)
                error('Near-end and far-end signals must have the same length.');
            end

            N  = length(x1);
            Fs = v_ne_tx.fs;

            % --- Use zero padding to avoid circular convolution artifacts ---
            Nfft = 2^nextpow2(2*N);  % robust choice

            % --- FFT ---
            V1 = fft(x1, Nfft);
            V2 = fft(x2, Nfft);

            % --- Frequency axis for interpolation (signed then abs) ---
            f_fft = (0:Nfft-1).' * (Fs/Nfft);
            f_signed = f_fft;
            idxNeg = f_signed > Fs/2;
            f_signed(idxNeg) = f_signed(idxNeg) - Fs;      % (-Fs/2, Fs/2]
            f_pos = abs(f_signed);

            % --- Interpolate S-parameters onto f_pos ---
            S31 = interp1(f_s, S31_s, f_pos, 'linear', 'extrap');
            S41 = interp1(f_s, S41_s, f_pos, 'linear', 'extrap');
            S32 = interp1(f_s, S32_s, f_pos, 'linear', 'extrap');
            S42 = interp1(f_s, S42_s, f_pos, 'linear', 'extrap');

            % Hermitian symmetry for real time-domain response:
            % For negative frequencies enforce conj symmetry.
            S31(idxNeg) = conj(S31(idxNeg));
            S41(idxNeg) = conj(S41(idxNeg));
            S32(idxNeg) = conj(S32(idxNeg));
            S42(idxNeg) = conj(S42(idxNeg));

            % --- Physical hybrid differential output ---
            % Vhyb = (S41-S31)*V1 + (S42-S32)*V2
            He = (S41 - S31);              % residual echo transfer from local TX
            Hr = (S42 - S32);              % transfer from far-end TX to output

            V_hyb = He .* V1 + Hr .* V2;

            % --- Residual echo estimate (digital canceller model) ---
            V_echo = He .* V1;

            % --- Back to time-domain (take first N samples after padding) ---
            v_hyb_full = ifft(V_hyb, 'symmetric');
            v_echo_full = ifft(V_echo, 'symmetric');

            v_hyb = v_hyb_full(1:N);
            v_echo_est = v_echo_full(1:N);

            % --- Optional cancellation ---
            if obj.cancel_echo
                v_fe_rec = v_hyb - v_echo_est;
            else
                v_fe_rec = v_hyb;
            end

            if ~obj.return_intermediates
                v_hyb = [];
                v_echo_est = [];
            end

            % --- Bring output in the correct form ---
            v_fe_rec = Informationsignal(v_fe_rec,"fs",v_fe_tx.fs);
            v_hyb = Informationsignal(v_hyb,"fs",v_fe_tx.fs);
            v_echo_est = Informationsignal(v_echo_est,"fs",v_fe_tx.fs);

            if obj.plot
                rfplot(net)
            end
        end
    end
end