GCC-PHAT auf MATLAB (immer noch schwierig)

Question

Ich versuche, Zeitverzögerung Schätzung zwischen zwei akustischen Signalen mit GCC-PHAT auf MATLAB zu berechnen. Es ist jedoch immer noch eine Herausforderung, GCC-PHAT mit fft/ifft zu codieren, wobei die Zeitverzögerungsschätzung immer Null ist. Ich habe stundenlang gegoogelt, war aber noch nicht erfolgreich. Mein Code wird an einem einfachen Beispiel mit zwei sinusförmigen Wellen wie folgt demonstriert.

Fs = 8000; 
dt = 1/Fs; 
f1 = 100; 
tdelay = 0.002; 
t3 = (0:dt:(1)-dt)'; 
x = cos(2*pi*f1*t3); 
y = cos(2*pi*f1*(t3-tdelay)); 

n=2*length(x)-1; 
N=2^nextpow2(n); 

xfft = fft(x,N); 
yfft = fft(y,N); 

df=Fs/N; 
sampleIndex = -N/2:N/2-1; 
f=sampleIndex*df; 

R=xfft.*conj(yfft); 
c=R./abs(R); 
ic=fftshift(ifft(c)); 
figure; plot(f,ic) 
[~,ind]=max(abs(ic)); 
lagDiff=f(ind);      
timeDiff=lagDiff/Fs; 

Ich verstehe, dass die GCC-PHAT durch R codiert werden, und es funktioniert wie in diesem Beitrag erwähnt https://dsp.stackexchange.com/questions/31956/gcc-phat-generalized-cross-correlation-matlab

jedoch die R-Code nicht robust ist, wo es in einigen Fällen versagt. Zum Beispiel sind für den Code in der oben genannten Verknüpfung die Ergebnisse bei f=100Hz und tdelay=0.01 oder bei f=50Hz und tdelay=0.01 falsch.

Kann jemand, der erfolgreich GCC-PHAT auf Matlab gelöst hat oder bekannt ist, die gccphat(), eine lizenzierte Funktion, näher zu diesem Thema erklären? Ich habe viele Beiträge über dieses scheinbar einfache Problem gesehen, aber es ist immer noch ungelöst.

Answer 1

Die neue Phased Array Toolbox von Matlab ermöglicht die Verwendung der Funktion GCCPHAT für Ihre Signale.

timeDelay = gccphat([Y1(span),Y2(span),Y3(span),Y4(span),Y5(span),Y6(span),Y7(span)],Y0(span),FS); 

Die timeDelay-Variable hat 7 Werte in Bezug auf die Abtastfrequenz FS.

Answer 2

Ihr Code und Matlabs gccphat() unterscheidet sich bei der Extraktion der Phase. Während Sie den allgemeinen Ansatz implementiert haben, berechnen Sie die Exponentialfunktion des Winkels. Daher vermeidet man numerische Probleme im Zusammenhang mit der Division durch 0 in:

c = R./abs (R);