# GCC-PHAT cross correlation

I am trying to implement GCC-PHAT in python.

It seems the only difference between GCC-PHAT and normal cross-correlation using FFT is the division by the magnitude.

Here is my code:

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import rfft, irfft, fftfreq, fft, ifft

def xcorr_freq(s1,s2):
f_s1 = fft(s1)
f_s2 = fft(s2)
f_s2c = np.conj(f_s2)
f_s = f_s1 * f_s2c
denom = abs(f_s)
denom[denom < 1e-6] = 1e-6
f_s = f_s / denom  # This line is the only difference between GCC-PHAT and normal cross correlation
return np.abs(ifft(f_s))[1:]


I have checked by commenting out fs = fs / denom The function produces the same result as normal cross correlation for a wide band signal.

Here is a sample test code which shows the GCC-PHAT code above performs worse than normal cross-correlation:

LENG = 500
a = np.array(np.random.rand(LENG))
b = np.array(np.random.rand(LENG))
plt.plot(xcorr_freq(a,b))
plt.figure()
plt.plot(np.correlate(s1,s2,'full'))


Here is the result with GCC-PHAT: Here is the result with normal cross-correlation: Since GCC-PHAT should give better cross-correlation performance for wide band signal, I know there is something wrong with my code. Any help is very appreciated!

• Concerning your second code fragment: how did you compute s1 and s2? I think that with np.correlate you compute something totally different from what you get with xcorr_freq. Your code for xcorr_freq looks OK. Sep 4, 2014 at 9:07

Yep, as Matt L. state in their comment your code is correct... As can be seen in the figures you've provided your GCC PhaT works quite well!

If you have a look in the literature, you'll realise that the cross-correlation of two independent random vectors, like your vectors, doesn't peak somewhere, since they are "purely" random and their samples are not correlated in any way. Below you can see the results I get from "simple" cross correlation implemented in the frequency domain, like you did. The random vectors are drawn from a uniform distribution. According to its documentation np.correlate does perform correlation, albeit in the time domain. Now, you get "unexpectedly" different results because you have different expectations. This means that the results you get should be what they are depending on what the signals s1 and s2 are. In the FFT version of your correlation (not sure if you still have the "PhaT version" commented or not) you use a and b while you use s1 and s2 for your time domain version. Most probably the two pairs of variables do not hold the same vectors.

$$y_{PhaT} \left( l \right) = \mathcal{IFFT} \left\{ \frac{G_{x y}}{\left| G_{x y} \right|} \right\}$$
where $$y_{PhaT}$$ is the PhaT cross correlation function, $$l$$ is the lag index (the free variable of the function), $$\mathcal{IFFT}$$ denotes the Inverse Fourier Transform, $$G_{x y}$$ the cross spectrum and $$\left| ~ \cdot ~ \right|$$ denotes magnitude.
The idea is that you only retain the phase information of the function, since time is entirely encoded in the phase. Thus, instead of dividing by the magnitude to, correctly in the mathematical sense, replicate the formula, you can just use the phase of the cross spectrum, which you can get with np.angle(f_s). This way you avoid division with zero, or the need to clamp the values of the magnitude to avoid numerical overflows, or "explosions". Consequently, you can calculate the GCC PhaT functions as ifft(np.exp(1j * np.angle(f_s))).