Since the OP mentions an index of frequency I assume a missing step not shown is an FFT is computed on the M samples, and the algorithm is to determine the maximum frequency.
I would not recommend such an approach to determine Doppler offset and fine frequency tuning except for initial acquisition (course tuning), given the computation involved in the FFT compared to what is alternatively explained below. Also I would expect that the maximum peak of an FFT on any given capture would be quite variable for modulated waveforms.
A favored algorithm is the cross product frequency discriminator which is the imaginary term of the complex conjugate product of adjacent time domain samples of the complex baseband waveform (or can also be applied to correlator outputs etc), where baseband waveform is the modulated waveform with only a residual carrier offset to be determined. This works as the imaginary term of the complex conjugate cross product is proportional to the phase between two complex phasors. If the two phasors are values separated in time, then the result is a change in phase over a change in time, which is frequency (as $d\phi/dt$). What is beautiful about this approach is that the imaginary term of the complex conjugate of two phasors given as $I_1+jQ_1$ and $I_2 + jQ_2$ is simply $I_1 Q_2 - I_2Q_1$, or $I[n] Q[n-1] - I[n-1]Q[n]$ for each time domain sample $n$ of the complex waveform. Such a result can be averaged for very fine frequency offset correction.
The imaginary term is the scaled sine of the angle (scaled by the square of the input amplitude), and for small angles $\sin(\phi) \approx \phi$, so can be effectively used in a recovery loop which would drive that angle to zero. It is typical to AGC (automatic gain control) the waveform first so that the signal amplitude is not controlling the loop time constant as it otherwise would.
This is often used in a decision directed carrier recovery loop such as shown in the block diagram below specific to a QPSK demodulator:
Note that the autocorrelation function is composed of conjugate multiplications of a signal with delayed versions of itself (!), so the imaginary term of the complex autocorrelation result is itself a frequency discriminator (cool!), as demonstrated in the plot below. There are many more computations than necessary in an autocorrelation (which can also be efficiently computed using FFT's as $\text{fft}(x)\text{fft}^*(x)$ where $\text{fft}^*(x)$ is the complex conjugate of the fft result) but wanted to include this in case similar processing was being done for other reasons such as spread-spectrum receivers or pilot detection, or if this can otherwise facilitate acquisition. For ongoing carrier tracking, I would recommend the simple $I[n] Q[n-1] - I[n-1]Q[n]$ frequency error discriminator.
