7
$\begingroup$

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.

$\endgroup$
6
  • $\begingroup$ 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? $\endgroup$
    – Peter K.
    Jul 5, 2016 at 13:57
  • $\begingroup$ @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. $\endgroup$
    – CN railfan
    Jul 5, 2016 at 14:38
  • $\begingroup$ 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. $\endgroup$
    – Peter K.
    Jul 5, 2016 at 14:49
  • $\begingroup$ @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. $\endgroup$
    – CN railfan
    Jul 5, 2016 at 14:54
  • $\begingroup$ Would you support this: area51.stackexchange.com/proposals/94189/matlab Just up vote question with less than 10 votes. $\endgroup$
    – Royi
    Jan 2, 2017 at 12:10

3 Answers 3

2
$\begingroup$

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
$\endgroup$
0
$\begingroup$

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.

CCF comparison


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)
$\endgroup$
4
  • $\begingroup$ 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). $\endgroup$
    – CN railfan
    Jul 5, 2016 at 14:59
  • $\begingroup$ R uses <- as the "assignment" operator, so yes, Fs <- 8000 is equivalent to matlab's Fs = 8000;. $\endgroup$
    – Peter K.
    Jul 5, 2016 at 15:11
  • $\begingroup$ Thanks for clarifying. I will adapt this for MATLAB and let you know how I got on. cheers! $\endgroup$
    – CN railfan
    Jul 5, 2016 at 15:12
  • $\begingroup$ @RoryHand : I basically just cut-and-pasted your code and then edited it to work in R. $\endgroup$
    – Peter K.
    Jul 5, 2016 at 15:16
0
$\begingroup$

You can add a constant factor to the weighting function to avoid amplifying the error when the denominator is small,like that c=R./(abs(R)+a); where a can be a constant

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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