function tau_error = modifiedGodardTimingRecovery(rx, N, eta, beta)
% modifiedGodardTimingRecovery
%
% This function estimates the symbol timing error using the modified Godard
% approach in the frequency domain as described in:
%
% "Modified Godard Timing Recovery for Non-Integer Oversampling Receivers"
% Appl. Sci. 2017, 7, 655. :contentReference[oaicite:0]{index=0}​:contentReference[oaicite:1]{index=1}
%
% Inputs:
%   rx   - Received time-domain signal (vector)
%   N    - FFT size (should be an even integer)
%   eta  - Effective oversampling factor used for timing recovery (eta > 1)
%   beta - Roll-off related parameter (0 < beta <= 1)
%
% Output:
%   tau_error - Estimated timing error (in sample units)
%
% Implementation Notes:
%   1. The function computes an N-point FFT of the first N samples of rx.
%   2. It then determines an offset (Delta) defined as:
%          offset = round((1 - 1/eta) * N)
%   3. To avoid index overflow, the summation is taken over indices k from 1 to
%      floor(N/2) - offset.
%   4. The timing error is estimated as:
%          tau_error = ( (1+beta)/(2*eta*N - 1) * sum(phase difference) ) / (2*pi)
%      where the phase difference is (angle(R(k)) - angle(R(k+offset)))
%
% Make sure that the input signal rx contains at least N samples.

% Check input length
if length(rx) < N
    error('Input signal length must be at least N.');
end

% Compute the N-point FFT of the first N samples of rx
R = fft(rx(1:N), N);

% Determine the offset based on the oversampling factor (eta)
offset = round((1 - 1/eta) * N);

% Define the summation range to avoid index overflow
k_min = 1;
k_max = floor(N/2) - offset;
if k_max < k_min
    error('Chosen parameters result in an empty summation range. Adjust N, eta, or beta.');
end

% Compute the sum of phase differences over the selected frequency bins
phase_diff_sum = 0;
for k = k_min:k_max
    phase_k = angle(R(k));
    phase_k_offset = angle(R(k + offset));
    phase_diff_sum = phase_diff_sum + (phase_k - phase_k_offset);
end

% Normalization factor as per the modified Godard algorithm
norm_factor = (1 + beta) / (2 * eta * N - 1);

% Estimate the timing error in sample units
tau_error = (norm_factor * phase_diff_sum) / (2 * pi);
end