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