3
$\begingroup$

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? :)

$\endgroup$
2
$\begingroup$

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  
$\endgroup$
  • $\begingroup$ Regarding the padding, I agree. However, conjugation in frequency domain or a reverse operation in time domain before should be the same, right? So if I remove the Conj function in your code above and use xc<-fftXcor(x, rev(y)), it should create the same result, but I get the highest value at the end. The sequence is shifted and I do not understand why regardless the padding. $\endgroup$ – PeteChro Apr 20 '15 at 20:50
  • $\begingroup$ @PeteChro Unfortunately I don't have a good answer to this either. It's unclear to me why padding and conjugation in the time and frequency domain, respectively, would be equivalent. Hopefully someone with more experience can chime in! $\endgroup$ – Matt SM Apr 22 '15 at 0:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.