# GCC-PHAT (Generalized cross correlation MATLAB)

I've been trying to implement generalized cross correlation with a PHAT weighting function for a while now, and cannot get it to work. I have tried correlating using MATLAB's xcorr.m file, and it works for getting a correct delay tau (on simulated sinusoidal signals).

The code for the signals are thus:

Fs = 8000;
dt = 1/Fs;%0.125e-3
f1 = 100;
tdelay = 0.625e-03;%try different values
t3 = (0:dt:(1)-dt)';
x3 = cos(2*pi*f1*t3);
x4 = cos(2*pi*f1*(t3-tdelay));


As can be seen in part of the source code for xcorr.m, the cross correlation is implemented thus:

%Transform both vectors

X = fft(x,2^nextpow2(2*M-1));
Y = fft(y,2^nextpow2(2*M-1));

% Compute cross-correlation

c = ifft(X.*conj(Y));


According to definitions of GCC-PHAT, the only addition I needed to make was to divide the product by its own magnitude, before taking the ifft. Here is my version with this change.

%Transform both vectors
X = fft(x,2^nextpow2(2*M-1));
Y = fft(y,2^nextpow2(2*M-1));

% Compute cross-correlation

R = X.*conj(Y);
c = ifft(R./abs(R));


However I always end up with a tau of zero with the PHAT weighting! On closer examination of the array produced as a result of this division, it seems as if the first value of R is a real value (with no imaginary component), and so when divided by its magnitude it becomes 1. All the other values in array R are complex, so do not end up as 1 when divided by their own magnitude and so end up with a value of <1.

This can be seen below, for the 1st 10 values of R.

    K>> R(1:10,1)

ans =

0.000000000000000 + 0.000000000000000i
-0.494299608718696 - 0.003002230689022i
-0.002678647083223 - 0.000032538742345i
-0.488954228290329 - 0.008909374553649i
-0.010656518992354 - 0.000258902698589i
-0.478379290671260 - 0.014528074329782i
-0.023760667475633 - 0.000865926459320i
-0.462803929640386 - 0.019677623220519i
-0.041707017319469 - 0.002026674993917i
-0.442565618721743 - 0.024194329448597i

K>> abs(R(1:10,1))

ans =

0.000000000000000
0.494308725968464
0.002678844707371
0.489035391682370
0.010659663580139
0.478599844010494
0.023776441018801
0.463222070011989
0.041756229537848
0.443226457301486

K>> R(1:10,1)./abs(R(1:10,1))

ans =

1.000000000000000 + 0.000000000000000i
-0.999981555555690 - 0.006073594357736i
-0.999926227844417 - 0.012146557900713i
-0.999834033705084 - 0.018218261306199i
-0.999705001216859 - 0.024288074069393i
-0.999539169638280 - 0.030355367874845i
-0.999336589392997 - 0.036419515378054i
-0.999097322000240 - 0.042479891383435i
-0.998821440083941 - 0.048535871565711i
-0.998509027227826 - 0.054586834901284i


As can be seen above, the largest value ends up being at the 1st index, when we divide R by its own magnitude. So once the ifft is taken, the highest value is ALWAYS at the beginning of the array, which gives a lag and timediff of zero...even when I set the delay between the 2 identical signals at the beginning to a nonzero value (e.g. delay = 0.75e-03).

What am I doing wrong?? Any help appreciated.

Rory

P.S. if anyone is wondering why I'm bothering with the PHAT weighting, its because it should give far better results (in theory) in a real life scenario, for TDOA.

• This sentence doesn't make sense to me: All the other values in array R are complex, so do not end up as 1 when divided by their own magnitude and so end up with a value of <1. Surely any number, complex or otherwise, when divided by its own magnitude will end up as a unit-long value? – Peter K. Jul 5 '16 at 13:57
• @PeterK. have tried to clarify what I was saying by editing my question; I have added in the 1st 10 values of R, its absolute value and also of R divided by its absolute value. – CN railfan Jul 5 '16 at 14:38
• What is M in your code? Perhaps the lengths of the FFTs are too short? No need to make it a power of two, just make it at least as long as length(x3)+length(x4)-1. – Peter K. Jul 5 '16 at 14:49
• @PeterK. M is the length of x3 or x4, the MATLAB function I am adapting declares M as the longest of either of these. The original function also uses a power of 2 in the ffts, but yes you have a point that length(3) + length(x4) -1 would be long enough for the ffts. – CN railfan Jul 5 '16 at 14:54
• Would you support this: area51.stackexchange.com/proposals/94189/matlab Just up vote question with less than 10 votes. – Royi Jan 2 '17 at 12:10

I remember I was facing a similar problem some time back. I am pretty sure it has to do with the DC and nyquist values when performing the equalisation and then the inverse. The first element in the array which is becoming 1.0 is the DC value.

I have this code which I’m pretty sure works for GCC-PHAT:

function [r,tau] = gccphat( x, y, fs )
M = max(numel(x),numel(y));

%%Transform both vectors
% X = fft(x,2^nextpow2(2*M-1));
% Y = fft(y,2^nextpow2(2*M-1));
%
% % Compute cross-correlation
%
% R = X.*conj(Y);
% c = ifft(R./abs(R));

%%
N = 2*M-1;
Nfft = 2^nextpow2(N);

R = bsxfun(@times, ...
fft(y,Nfft), ...
conj(fft(x,Nfft)));
rtmp = fftshift( ...
ifft(exp(1i*angle(R))) ,1);
r = rtmp(Nfft/2+1-(M-1)/2:Nfft/2+1+(M-1)/2,:);

lags = (-(N-1)/2:(N-1)/2).';
lags = lags/fs;
[~,idx] = max(abs(r));
tau = N/(2*fs)+lags(idx);

end


The algorithm you use works fine for me in R.

The image below shows the standard cross correlation function (CCF) and the generalized CCF. Both give the right answer, but the generalized has a much more distinct peak.

R Code Below

#31956

Fs <- 8000
dt <- 1/Fs
f1 <- 100
tdelay <- 0.625e-03
t3 <- seq(0,1-dt,dt)
x3 <- cos(2*pi*f1*t3)
x4 <- cos(2*pi*f1*(t3-tdelay))

xcorr_31956 <- function(x,y,normalize = FALSE)
{
xfft <- fft(x, 4*length(x))
yfft <- fft(y, 4*length(x))

R <- xfft*Conj(yfft);
if (normalize)
{
R <- R/abs(R)
}
c <- fft(R, inverse=TRUE);

return(c)
}

xc <- xcorr_31956(x3,x4, FALSE)
xc_phat <- xcorr_31956(x3,x4, TRUE)

par(mfrow=c(2,1))
plot(seq(0,length(xc)-1),abs(xc), type="l", xlim=c(0,20), col="blue", lwd=2)
ix <- which.max(abs(xc))
points(ix-1,abs(xc[ix]), col="red", lwd=5);
title('Standard CCF')

plot(seq(0,length(xc)-1),abs(xc_phat), type="l", xlim=c(0,20), col="blue", lwd=2)
ix_phat <- which.max(abs(xc_phat))
points(ix_phat-1,abs(xc_phat[ix_phat]), col="red", lwd=5);
title('Generalized')

print(paste("Delay is",tdelay*Fs), quote = FALSE)
print(paste("Estimate is",which.max(abs(xc)) - 1), quote = FALSE)

• Hi Peter thanks for your answer. Can I clarify that the ' <- ' symbols in the above code should actually be ' = ' symbols? (So as to assign variables a value, e.g. Fs = 8000). – CN railfan Jul 5 '16 at 14:59
• R uses <- as the "assignment" operator, so yes, Fs <- 8000 is equivalent to matlab's Fs = 8000;. – Peter K. Jul 5 '16 at 15:11
• Thanks for clarifying. I will adapt this for MATLAB and let you know how I got on. cheers! – CN railfan Jul 5 '16 at 15:12
• @RoryHand : I basically just cut-and-pasted your code and then edited it to work in R. – Peter K. Jul 5 '16 at 15:16