Auto-correlation/correlation in C++

Question

The following function computes the correlation between two vectors.

It doesn't give the same result as matlab function for small values:

I am really don't know if the bug becomes from this function or not ? the maximum lags by default is N-1 ? is this reasonable ?

inline int pow2i(int x) { return ((x < 0) ? 0 : (1 << x)); }`

     vec xcorr(vec x, vec y,bool autoflag)
    {
      int maxlag=0;
      int N = std::max(x.size(), y.size());
      //Compute the FFT size as the "next power of 2" of the input vector's length (max)
      int b = ceil(log2(2.0 * N - 1));
      int fftsize = pow2i(b);
    
      int e = fftsize - 1;
      cx_vec temp2;
    
      if (autoflag == true) {
        //Take FFT of input vector
        cx_vec X = cx_vec(x,zeros(x.size()));
        X= fft(X,fftsize);
        //Compute the abs(X).^2 and take the inverse FFT.
        temp2 = ifft(X%conj(X));
      }
      else{
       //Take FFT of input vectors
      cx_vec X=cx_vec(x,zeros(x.size()));
      cx_vec Y=cx_vec(y,zeros(y.size()));
      X = fft(X,fftsize);
      Y = fft(Y,fftsize);
      //cout<< "Y " << Y << endl;
      //cout<< "X " << X<< endl;
      temp2 =ifft(X%conj(Y));
      //cout<< "temp 2 " << temp2 << endl;
     }
      maxlag=N-1;
      vec out=real(join_cols(temp2(span(e - maxlag + 1, e)),temp2(span(0,maxlag))));
      return out;
    }

Possibly related https://stackoverflow.com/questions/588004/is-floating-point-math-broken — JohnFilleau, Jul 05 '20 at 15:48
Floating-point nitpicks aside, are you certain that rounding the FFT window size up is appropriate for this specific task? For small inputs, why don't you just directly compute the correlation in the time domain and skip the overhead of the Fourier transform? — nanofarad, Jul 05 '20 at 15:52
vec is the vector definition using armadillo framework http://arma.sourceforge.net/docs.html#Col — Mouna Karoui, Jul 05 '20 at 15:55
@nanofarad: I need something general that work for all cases, "For small inputs, why don't you just directly compute the correlation in the time domain and skip the overhead of the Fourier transform?" --> how can I do this ? — Mouna Karoui, Jul 05 '20 at 15:58
@MounaKaroui Equation 3 [here](https://en.wikipedia.org/wiki/Cross-correlation#Cross-correlation_of_deterministic_signals). This is the original formula, i.e. how cross-correlation is **defined**. Any weird games with the FFT are primarily attempts to improve performance for large vectors, or to show off one's signal processing knowledge. As long as you have the time and memory, and avoid overflow, the original time-domain formula will always work for any size input. — nanofarad, Jul 05 '20 at 16:03

score 2 · Answer 1 · answered Jul 08 '20 at 19:46

2

Just implement autocorrelation in time-domain, as one of the comments mentioned.

Armadillo does not have cross-correlation (and autocorrelation) implemented, but one easy way to implement them is using convolution, which armadillo does have. You just need to invert the other of the elements in the second vector and arma::conv will be essentially be computing the cross-correlation.

That is, you can easily compute the autocorrelation with of an arma::vec a with

arma::vec result = arma::conv(a, arma::reverse(a));

This gives the same result that xcorr in MATLAB/Octave returns (when you pass just a single vector to xcorr it computes the autocorrelation). Note that you might want to divide result by N or by N-1.

answered Jul 08 '20 at 19:46

darcamo

3,294
1
16
27

Hi, I tried your solution. It 's true gives the same values as MATLAB for not so small values However, It doesn't give the same values when I tried with the following input {1.056e-11, 1.0843e-11,1.113e-11,1.1431e-11,1737e-11}. – Mouna Karoui Jul 09 '20 at 06:52
And I can't understand the reason of this difference, from the previous comments, I tried to learn about floating point. But I didn't get any solution for that problem for my case. I need to convert a MATLAB code to C++ and being sure that it works well without any bugs. – Mouna Karoui Jul 09 '20 at 07:02
I get the same response no matter if I use C++ with armadillo or octave for any input, including this vector that you mentioned. The math is right and any discrepancies must be due to numerical precision. – darcamo Jul 10 '20 at 17:56

Auto-correlation/correlation in C++

1 Answers1