classdef Godard_Timing_Recovery < handle

    properties(Access=public)
        mode
        num_blocks
        fft_length
        sps
        rolloff
        mu
        Ki
    end

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

                options.mode = 0;
                options.num_blocks = 1;
                options.fft_length = 1024;
                options.sps = 2;
                options.rolloff = 0;
                options.mu = 0;
                options.Ki = 1e-3;

            end

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

            %obj-Initialization here%

        end

        function [data_out, tau_hat] = process(obj, data_in)

            data_out = data_in;
            output_vector = zeros(length(data_in),1);
            block_length = length(data_in)/obj.num_blocks;
            for i = 1:obj.num_blocks
                i_start = ((i-1)*block_length)+1;
                i_end = i*block_length;
                x = data_in.signal(i_start:i_end);
                beta = obj.rolloff;
                eta  = obj.sps;
                R = fft(x,obj.fft_length);
                R_full = fft(x);
            
                if obj.mode == 0 % Classic Godard
                    % Calculate time shift
                    k0   = (1:obj.fft_length/2).';
                    idx1 = k0;
                    idx2 = k0+(obj.fft_length/2);
                    tau_hat = sum(imag(R(idx1) .* conj(R(idx2))));
                elseif obj.mode == 1 || obj.mode == 2 % Modified Godard 1
                    % Calculate Shift for the received signal
                    shiftBins = (1 - 1/eta) * obj.fft_length;
                    if abs(shiftBins - round(shiftBins)) > 1e-12
                        disp('Warning: shiftBins=(1-1/eta)*fft_length is non-integer. Choose compatible values for fft_length and eta.');
                    end
                    shiftBins = round(shiftBins);
                
                    % Calculate upper and lower bounds
                    kStart = ((1-beta)/(2*eta)) * obj.fft_length;
                    kEnd   = ((1+beta)/(2*eta)) * obj.fft_length - 1;
                    if abs(kStart - round(kStart)) > 1e-12 || abs(kEnd - round(kEnd)) > 1e-12
                        disp('Warning: kStart/kEnd are non-integer. Choose compatible values for fft_length, eta, and beta.');
                    end
                    kStart = round(kStart);
                    kEnd   = round(kEnd);
                
                    k0   = (kStart:kEnd).';
                    idx1 = k0 + 1;
                    idx2 = mod(k0 + shiftBins, obj.fft_length) + 1;
        
                    % Calculate time shift 
                    if obj.mode == 1
                        tau_hat = sum(imag(R(idx1) .* conj(R(idx2))));
                    elseif obj.mode == 2
                        tau_hat = sum(angle(R(idx1))-angle(R(idx2)));
                    end
                end

                
                % % Normalization
                % denom = floor(log10(abs(tau_hat)));
                % tau_hat = tau_hat / 10^denom;

                % Calculate mu using a first-order loop filter
                mu_block = obj.mu + obj.Ki * tau_hat;

                % % Shifting the signal in time domain using interpolation
                % x_original = linspace(0,length(x)-1,length(x)).';
                % x_new = x_original + mu_block;
                % output_vector(i_start:i_end) = interp1(x_original,x,x_new,'linear','extrap');
    
                % Shifting the signal in frequency domain
                k = (0:block_length-1).';
                phaseRamp = exp(-1j * 2*pi * (k/block_length) * mu_block);
                R_shifted = R_full .* phaseRamp;

                output_vector(i_start:i_end) = real(ifft(R_shifted));

            end
            data_out.signal = output_vector;
        end


    end
end