More focus on Database analysis and direct DSP'ing of run_id's

Functions/Theory/modifiedGodardTimingRecovery.m

@@ -0,0 +1,62 @@
+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