R Own Implementation of Cross Correlation using Convolution

Question

I try to implement my own cross correlation function in R by translating it as a convolution problem.

Part I: So I have two arrays, e.g. two identical arrays, and I want to get the cross correlation in R, then I need the following code?!:

a1 = 1:9
a2 = 1:9

# Now translate into a conv. problem
a2 = rev(a1)
F2 = fft(a2)
F1 = fft(a1)
FR = F1 * F2
Re(fft(FR,inverse=TRUE))/length(FR)

The result is: 249 222 204 195 195 204 222 249 285

As I am working with two identical arrays, I would expect that I get the highest correlation value at position zero. If I calculate the same problem on wolfram alpha, I get the sequence 285 249 222 ... as expected.

Part II: Normalization In order to get normalized values, I need to subtract the mean and divide by the standard deviation:

a1 = ( a1 - mean(a1) ) / sd(a1)
a2 = rev(a1)

Then I get the following values: 3.2 -0.4 -2.8 -4.0 -4.0 -2.8 -0.4 3.2 8.0 although the correlation values should be between -1 and 1.

So what are my mistakes? :)

Answer 1

I think you're missing two things. First, you need to take a complex conjugate of fft(x) before taking the inverse FFT. More important, and subtle, is that the FFT assumes periodicity in your data. As a result, if you calculate the cross correlation directly you're calculating correlations with wrap around, which isn't what you want I suspect. For example, for 1:5 and a lag of 1, instead of calculating correlation between
1 2 3 4 5
0 1 2 3 4
you're calculating it between
1 2 3 4 5
5 1 2 3 4
You can account for this by padding your vectors with 0 up to a size of $2n-1$. There's a nice explanation of this in: Fast variogram computation with FFT (Marcotte 1996); pdf here. This paper focuses on 2D spatial data, but I think the idea is the same.

Here's some R code to calculate cross-correlation for lags $-(n-1)$ to $n-1$:

fftXcor <- function(x, y) {
    n <- length(x)
    # Normalize
    x <- as.numeric(scale(x)) 
    y <- as.numeric(scale(y))
    # Enlarge with 0's to size 2*n-1 to account for periodicity
    x <- c(x, rep(0, length(x) - 1))
    y <- c(y, rep(0, length(y) - 1))
    # FFT
    xfft <- fft(x)
    yfft <- fft(y)
    # Cross-correlation via convolution
    crosscor <- fft(Conj(xfft) * yfft, inverse=T) / length(x)
    crosscor <- Re(crosscor) / (n - 1)
    # Slice it up to make it for lags -n:n not 0:(2n-1)
    crosscor <- c(crosscor[(n+1):length(crosscor)], crosscor[1:n])
    # Store lag as names attribute of vector
    names(crosscor) <- (1-n):(n-1)
    return(crosscor)
}
x <- 1:9
y  <- 9:1
xc <- fftXcor(x, y)
# Just look at half the vector since it's symmetric
round(xc, 3)[9:17]
#      0      1      2      3      4      5      6      7      8 
# -1.000 -0.667 -0.350 -0.067  0.167  0.333  0.417  0.400  0.267 

# Compare to R's built in cross correlation function
xc_ccf <- ccf(x, y, lag.max=8, plot=F, type='correlation')
as.vector(xc / xc_ccf$acf)
# [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1