Write a function to implement the GCC-PHAT algorithm (generalized cross-correlation with phase transform). The input signal and reference signal lists must be the same length N. The returned delay list is of length (2 * N) - 1, with the element representing zero delay at index N. Indices less than N represent negative delays, and indices greater than N represent positive delays.

P. McGill
10Apr25

gccPhat[sig_List,ref_List]:=Module[{sigt,reft,prod,denom,delay}, (* version 1.0.0 of function *)
If[Length@sig!=Length@ref,Print["Error: Signal and reference must be same length"];Return[]];
sigt=Fourier[Join[sig,Table[0,Length[sig]]]];(* transform to frequency domain *)
reft=Fourier[Join[ref,Table[0,Length[ref]]]];
prod=sigt Conjugate[reft]; (* multiply transforms *)
denom=Clip[Abs[prod],{(10^(-6)),Infinity}]; (* find normalization and clip small values *)
delay=Abs[InverseFourier[prod/denom]];
Return[Join[delay[[Length[delay]/2+2;;]],delay[[;;Length[delay]/2]]]]
] (* end Module *)

Create a simulated earthquake signal.

swave=Table[(E^(-(x/25))) Sin[x],{x,0,50 Pi}];
pwave=Table[0.1 Sin[x],{x,0,10 Pi}];
simQuake=Join[Table[0,50],pwave,swave];

ListLinePlot[simQuake,PlotRange->All]

Length[simQuake]

Create a delayed version of the simulated earthquake.

dlySimQuake=Drop[Join[Table[0,25],simQuake],-25];

ListLinePlot[{simQuake,dlySimQuake},PlotRange->All]

Test the gccPhat function.

xcor=gccPhat[dlySimQuake,simQuake];

Length@xcor

Positive delay.

ListLinePlot[xcor[[190;;290]],PlotRange->All,DataRange->{-50,50}]

No delay.

xcor2=gccPhat[simQuake,simQuake];
ListLinePlot[xcor2[[190;;290]],PlotRange->All,DataRange->{-50,50}]

Negative delay.

xcor3=gccPhat[simQuake,dlySimQuake];
ListLinePlot[xcor3[[190;;290]],PlotRange->All,DataRange->{-50,50}]

Input lists not the same length.

xcor3=gccPhat[simQuake,dlySimQuake[[2;;]]];

Plot a standard correlation for comparison.

ListLinePlot[ListCorrelate[simQuake,dlySimQuake,1],PlotRange->All]

Add noise to the simulated earthquakes to see how that affects GCC-PHAT performance.

simQuakeN=simQuake+Table[RandomReal[{-0.1,0.1}],Length[simQuake]];
dlySimQuakeN=dlySimQuake+Table[RandomReal[{-0.1,0.1}],Length[dlySimQuake]];

ListLinePlot[simQuakeN,PlotRange->All]

xcor4=gccPhat[dlySimQuakeN,simQuakeN];
ListLinePlot[xcor4[[190;;290]],PlotRange->All,DataRange->{-50,50}]