Problem with Generalized Cross Correlation (PHAT filter)

Question

I have to calculate time delay between two signals.I am using the GCC-PHAT method. First, I cut the signals into smaller parts (~0.5s) and then try to calculate the delay. When I display those two signals on the graph, the delay is very obvious.

The problem is that for some parts of that long signal, my program calculates the correct delay, but for some parts the delay is totally incorrect. In both cases, when I display those two signals on the graph, the delay seems very obvious.

Relevant code:

duration=20000
a=hanning(duration+1)
[signal1,Fs] = audioread(mic1);
[signal2,Fs] = audioread(mic2);

start=1000;
stop=start+duration;
x1=[a.*signal1(start:stop)];
x2=[a.*signal2(start:stop)];

A=fft(x1);
B=fft(x2);

PHAT=A.*conj(B)./(abs((A).*conj(B)));   
result=abs(ifft(PHAT)); 

Examples: 1) Correct result:

2)Incorrect result:

I recorded those signals in a closed room, so there are probably also reverberations. Am I doing something wrong? Or is this method just not very effective in my circumstances?

I also tried to prefilter signals, but the results were the same. There is also no improvement, if I change the length of the signals.

Answer 1

Using the gcc-phat method, you can use the matlab function "gccphat(s2,s1,fs)" (R2016a) with s1 as the reference signal, s2 the other one and fs the sampling frequency.

Basically it will do a similar correlation calculation :

  Ncorr = 2*N-1;  % N size of the signal
  NFFT = 2^nextpow2(Ncorr);
  R12 = bsxfun(@times,fft(s2,NFFT),conj(fft(s1,NFFT)));
  r12_temp = fftshift(ifft(exp(1i*angle(R12))),1);
  r12 = r12_temp(NFFT/2+1-(N-1)/2:NFFT/2+1+(N-1)/2,:);

The only difference with your code is that the matlab function use the nearest power of two sequence length for the fft (as Peter.K suggest with "Nfft = length(x1) + length(x2) - 1") and doing an exponentiel stuff (and this part make the difference with a regular correlation, what you were doing)

But to return the time delay, the function is doing an other trick :

  lags = (-(Ncorr-1)/2:(Ncorr-1)/2).';
  lags = lags/fs;
  [~,idx] = max(abs(r12));
  tau = ts/(2*fs)+lags(idx); % ts as the total sample of the signal

I'm doing similar experimentation and this method is working great for long and short non-moving signal. the reverberation of the room may cause some issue, but if you can see an obvious delay when displaying the signal it should work fine.

Right now, i'm having difficulty to return a good time delay with a moving source (the function retourn 0 or +/- the length of the signal) And i'm still working on it... So if anybody got a lead, it will be much apreciate ;-)

Hope this will help you.

Peace

PS : Also, when doing an fft/ifft you really should use fftshift/ifftshift because the matlab fft/ifft is mixing some stuff during the process.

Answer 2

I suspect you have an issue with time domain aliasing.

Instead of

A=fft(x1);
B=fft(x2);

try

Nfft = length(x1) + length(x2) - 1;
A=fft(x1, Nfft);
B=fft(x2, Nfft);