# Understanding and converting Matlab's coarse frequency compensation function

I am trying to understand and convert to C Matlab's coarse frequency compensator the documentation for which can be found here: Documentation

Matlab's design itself is based on this IEEE paper: Carrier Frequency Recovery in All-Digital Modems for Burst-Mode Transmissions

I want to translate the correlation based algorithm for QPSK modulation; for example the invocation of such an object in Matlab would look like this:

cfc = comm.CoarseFrequencyCompensator('Modulation', 'QPSK', ...
'Algorithm', 'Correlation-based', ...
'MaximumFrequencyOffset', 12.5e3, ...
'SampleRate',400000);


The documentation states that the estimate is given by:

$$estimate=\frac{f_{samp}}{\pi(L+1)}arg\{\sum_{k=1}^{L}R(k)\}$$

where

$$L = round(\frac{f_{samp}}{f_{max}}) - 1$$

$$R(k) = \frac{1}{N-k}\sum_{i=k+1}^{N}r_{i}r^{*}_{i-k}$$, autocorrelation function,

$$r_k = e^{j(2\pi\Delta f k T_{s} + \theta)}$$, model of received signal

$$r^{*}_{k} = e^{-j(2\pi\Delta f k T_{s} + \theta)}$$, is the complex conjugate of $$r_{k}$$

$$N = \frac{L+1}{M}$$, number of samples

$$M = 4$$, because for QPSK there are four possible symbols

My C code:

double complex autocorrelation_func(double complex* input,
double complex* conjinput,
int k, int input_size)
{
double complex auto_corr_sum = 0+0*I;

for(int i = 0 + k; i < input_size; i++ )
{
auto_corr_sum += input[i]*conjinput[i];
}
auto_corr_sum *= (1/(input_size-k));

return auto_corr_sum;
}

double coarse_freq_comp(double complex* input,
double complex* output,
int input_size, int mod_order,
double max_freq_offset,
double sample_rate )
{
int auto_corr_seq = round(sample_rate/max_freq_offset)-1;
int num_samples = (auto_corr_seq+1)/mod_order;
double mult_const = sample_rate/((auto_corr_seq+1)*M_PI);
double complex input_conjugated[input_size];
double complex auto_corr_sum = 0 + 0*I;
double offset_est = 0;

for(int ndx = 0; ndx < input_size; ndx++)
{
input_conjugated[ndx] = conj(input[ndx]);
}

for(int i = 0; i < num_samples; i++)
{
auto_corr_sum += autocorrelation_func(input, input_conjugated, i, input_size);
}

offset_est = mult_const*catan(auto_corr_sum);

return offset_est;
}


Questions:

How do I apply the estimate to the input signal to get a compensated output signal?

What is wrong with my C code? I am not getting the same estimate value that I get for a given complex array in Matlab. I must be doing something wrong but this isn't obvious to me.