GCC-PHAT cross correlation

Question

I am trying to implement GCC-PHAT in python.

The approach is similar to the following two links: link1 and link2

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):
    pad1 = np.zeros(len(s1))
    pad2 = np.zeros(len(s2))
    s1 = np.hstack([s1,pad1])
    s2 = np.hstack([pad2,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!

Answer 1

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.

Additional info

Although your code for GCC PhaT seems to be OK (haven't tested it but all necessary steps seem to be present) there is some "better" approach to calculate the "PhaT filter". According to the filters formula

$$ 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))).

As a final remark, please note that if there's any numerical errors resulting in imaginary values in the time-domain (after the inverse transform) you should get the real part of it and not the absolute value/magnitude.