Forums

Generalized Cross Correlation

Started by ryujin_ssdt March 8, 2007
Following this thread http://www.dsprelated.com/showmessage/1737dd5/1.php
and the code from Davide Renzi at
http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=8581&objectType=file
I tried to implement a simple demo to learn the differences between all GCC
methods for time delay estimation.

I have some questions and request:

1 - Can someone run this script and tell me if I am on the rigth track? I
can see some nice correlation peaks that may indicate the delay but  I
cannot understand the Roth graph and the HT method does not produce any
results.

2 - The SCOT and cps-m methods seem to have a mirror image (I get two
peaks). Any explanations? maybe cyclic correlation?

3 - How can I make the graphs display the x axis in time (i.e. 3ms, 4ms).
Is any other way to extract the time delay value from the resulting cross
correlations?

4 - Since I use 8192 point FFT for cross correlating the result has 8192
sample values... can I simple graph from 0 to 1024 without lossing any
information?

5 - What is the best way to compare all these methods?? obtaining the time
delay estimate for different SNR values maybe? I tried it with the SCOT
method and got the same exact graph for all SNR (1~50dB) I tried... 

Here is the code:

FFTLength = 8192;

% Generate random data
x = randn(1, 1024);

% What is this filter for??
y = filter([1 1],1,x);
 
% Delay y a bit.
y = [zeros(1,150) y];
y = y(1:length(x));
 
% Get power spectra of x and y and cross power spectra
Rxx = xcorr(x);
Ryy = xcorr(y);
Rxy = xcorr(x,y);
Sxx = fft(Rxx,FFTLength);
Syy = fft(Ryy,FFTLength);
Sxy = fft(Rxy,FFTLength);
 
% Unfiltered Correlation (plain correlation) 
W = ones(1,FFTLength);
% Apply the filter
R = Sxy.*W;
% Obtain the GCC
G = fftshift(real(ifft(R)),1);
figure(1);
subplot(321);
plot(G);
title('Plain Time Cross Correlation');
 
% ROTH Filter
W = 1./Sxx;
% Apply the filter
R = Sxy.*W;
% Obtain the GCC
G = fftshift(real(ifft(R)),1);
figure(1);
subplot(322);
plot(G);
title('Roth Kernel GCC');

% SCOT Filter
W = 1./(Sxx .* Syy).^0.5;
% Apply the filter
R = Sxy.*W;
% Obtain the GCC
G = fftshift(real(ifft(R)),1);
figure(1);
subplot(323);
plot(G);
title('Scot Kernel GCC');

% PHAT Filter
W = 1./abs(Sxy);
% Apply the filter
R = Sxy.*W;
% Obtain the GCC
G = fftshift(real(ifft(R)),1);
figure(1);
subplot(324);
plot(G);
title('Phat Kernel GCC');

% cps-m Filter  (SCOT filter modified)
factor = .75;       % common value between .5 and 1
W = 1./(Sxx .* Syy).^factor;
% Apply the filter
R = Sxy.*W;
% Obtain the GCC
G = fftshift(real(ifft(R)),1);
figure(1);
subplot(325);
plot(G);
title('cpc-m Kernel GCC');

% Hannah and Thomson filter
gamma = Sxy./sqrt(Sxx .* Syy);
W = abs(gamma).^2 ./ (abs(Sxy).*(1-abs(gamma).^2));
% Apply the filter
R = Sxy.*W;
% Obtain the GCC
G = fftshift(real(ifft(R)),1);
figure(1);
subplot(326);
plot(G);
title('HT Kernel GCC');

>Following this thread
http://www.dsprelated.com/showmessage/1737dd5/1.php
>and the code from Davide Renzi at >http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId…81&objectType=file >I tried to implement a simple demo to learn the differences between all
GCC
>methods for time delay estimation. > >I have some questions and request: > >1 - Can someone run this script and tell me if I am on the rigth track?
I
>can see some nice correlation peaks that may indicate the delay but I >cannot understand the Roth graph and the HT method does not produce any >results. > >2 - The SCOT and cps-m methods seem to have a mirror image (I get two >peaks). Any explanations? maybe cyclic correlation? > >3 - How can I make the graphs display the x axis in time (i.e. 3ms,
4ms).
>Is any other way to extract the time delay value from the resulting
cross
>correlations? > >4 - Since I use 8192 point FFT for cross correlating the result has 8192 >sample values... can I simple graph from 0 to 1024 without lossing any >information? > >5 - What is the best way to compare all these methods?? obtaining the
time
>delay estimate for different SNR values maybe? I tried it with the SCOT >method and got the same exact graph for all SNR (1~50dB) I tried... > >Here is the code: > >FFTLength = 8192; > >% Generate random data >x = randn(1, 1024); > >% What is this filter for?? >y = filter([1 1],1,x); > >% Delay y a bit. >y = [zeros(1,150) y]; >y = y(1:length(x)); > >% Get power spectra of x and y and cross power spectra >Rxx = xcorr(x); >Ryy = xcorr(y); >Rxy = xcorr(x,y); >Sxx = fft(Rxx,FFTLength); >Syy = fft(Ryy,FFTLength); >Sxy = fft(Rxy,FFTLength); > >% Unfiltered Correlation (plain correlation) >W = ones(1,FFTLength); >% Apply the filter >R = Sxy.*W; >% Obtain the GCC >G = fftshift(real(ifft(R)),1); >figure(1); >subplot(321); >plot(G); >title('Plain Time Cross Correlation'); > >% ROTH Filter >W = 1./Sxx; >% Apply the filter >R = Sxy.*W; >% Obtain the GCC >G = fftshift(real(ifft(R)),1); >figure(1); >subplot(322); >plot(G); >title('Roth Kernel GCC'); > >% SCOT Filter >W = 1./(Sxx .* Syy).^0.5; >% Apply the filter >R = Sxy.*W; >% Obtain the GCC >G = fftshift(real(ifft(R)),1); >figure(1); >subplot(323); >plot(G); >title('Scot Kernel GCC'); > >% PHAT Filter >W = 1./abs(Sxy); >% Apply the filter >R = Sxy.*W; >% Obtain the GCC >G = fftshift(real(ifft(R)),1); >figure(1); >subplot(324); >plot(G); >title('Phat Kernel GCC'); > >% cps-m Filter (SCOT filter modified) >factor = .75; % common value between .5 and 1 >W = 1./(Sxx .* Syy).^factor; >% Apply the filter >R = Sxy.*W; >% Obtain the GCC >G = fftshift(real(ifft(R)),1); >figure(1); >subplot(325); >plot(G); >title('cpc-m Kernel GCC'); > >% Hannah and Thomson filter >gamma = Sxy./sqrt(Sxx .* Syy); >W = abs(gamma).^2 ./ (abs(Sxy).*(1-abs(gamma).^2)); >% Apply the filter >R = Sxy.*W; >% Obtain the GCC >G = fftshift(real(ifft(R)),1); >figure(1); >subplot(326); >plot(G); >title('HT Kernel GCC');
Hello, I run your program a little bit. Actually ,time delay is wrong. What may be the reason ? %%Building our signals , the same type signal tho=2; A=1; t=-15:0.1:15; SNR=20; x= rect( t+1,tho,A ); % First rectangle signal plot(t,x) hold on y=rect(t-5,tho,A); %Second rectangle signal ,delayed 5 unit 250 sample plot(t,y,'r') %%Correlation Part % Get power spectra of x and y and cross power spectra Rxx = xcorr(x); Ryy = xcorr(y); Rxy = xcorr(x,y); Sxx = fft(Rxx,FFTLength); Syy = fft(Ryy,FFTLength); Sxy = fft(Rxy,FFTLength); % Unfiltered Correlation (plain correlation) W = ones(1,FFTLength); % Apply the filter R = Sxy.*W; % Obtain the GCC G = fftshift(real(ifft(R)),1); figure; plot(G); %peak value of GCC give us time delay f=max(G); find(G==f) disp('time delay-cross correlation') title('Plain Time Cross Correlation'); time delay is 60 but the program gives values around 240. thanks --------------------------------------- Posted through http://www.DSPRelated.com
On Wednesday, February 24, 2016 at 1:13:00 AM UTC+13, mekazanc wrote:
> >Following this thread > http://www.dsprelated.com/showmessage/1737dd5/1.php > >and the code from Davide Renzi at > >http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId...81&objectType=file > >I tried to implement a simple demo to learn the differences between all > GCC > >methods for time delay estimation. > > > >I have some questions and request: > > > >1 - Can someone run this script and tell me if I am on the rigth track? > I > >can see some nice correlation peaks that may indicate the delay but I > >cannot understand the Roth graph and the HT method does not produce any > >results. > > > >2 - The SCOT and cps-m methods seem to have a mirror image (I get two > >peaks). Any explanations? maybe cyclic correlation? > > > >3 - How can I make the graphs display the x axis in time (i.e. 3ms, > 4ms). > >Is any other way to extract the time delay value from the resulting > cross > >correlations? > > > >4 - Since I use 8192 point FFT for cross correlating the result has 8192 > >sample values... can I simple graph from 0 to 1024 without lossing any > >information? > > > >5 - What is the best way to compare all these methods?? obtaining the > time > >delay estimate for different SNR values maybe? I tried it with the SCOT > >method and got the same exact graph for all SNR (1~50dB) I tried... > > > >Here is the code: > > > >FFTLength = 8192; > > > >% Generate random data > >x = randn(1, 1024); > > > >% What is this filter for?? > >y = filter([1 1],1,x); > > > >% Delay y a bit. > >y = [zeros(1,150) y]; > >y = y(1:length(x)); > > > >% Get power spectra of x and y and cross power spectra > >Rxx = xcorr(x); > >Ryy = xcorr(y); > >Rxy = xcorr(x,y); > >Sxx = fft(Rxx,FFTLength); > >Syy = fft(Ryy,FFTLength); > >Sxy = fft(Rxy,FFTLength); > > > >% Unfiltered Correlation (plain correlation) > >W = ones(1,FFTLength); > >% Apply the filter > >R = Sxy.*W; > >% Obtain the GCC > >G = fftshift(real(ifft(R)),1); > >figure(1); > >subplot(321); > >plot(G); > >title('Plain Time Cross Correlation'); > > > >% ROTH Filter > >W = 1./Sxx; > >% Apply the filter > >R = Sxy.*W; > >% Obtain the GCC > >G = fftshift(real(ifft(R)),1); > >figure(1); > >subplot(322); > >plot(G); > >title('Roth Kernel GCC'); > > > >% SCOT Filter > >W = 1./(Sxx .* Syy).^0.5; > >% Apply the filter > >R = Sxy.*W; > >% Obtain the GCC > >G = fftshift(real(ifft(R)),1); > >figure(1); > >subplot(323); > >plot(G); > >title('Scot Kernel GCC'); > > > >% PHAT Filter > >W = 1./abs(Sxy); > >% Apply the filter > >R = Sxy.*W; > >% Obtain the GCC > >G = fftshift(real(ifft(R)),1); > >figure(1); > >subplot(324); > >plot(G); > >title('Phat Kernel GCC'); > > > >% cps-m Filter (SCOT filter modified) > >factor = .75; % common value between .5 and 1 > >W = 1./(Sxx .* Syy).^factor; > >% Apply the filter > >R = Sxy.*W; > >% Obtain the GCC > >G = fftshift(real(ifft(R)),1); > >figure(1); > >subplot(325); > >plot(G); > >title('cpc-m Kernel GCC'); > > > >% Hannah and Thomson filter > >gamma = Sxy./sqrt(Sxx .* Syy); > >W = abs(gamma).^2 ./ (abs(Sxy).*(1-abs(gamma).^2)); > >% Apply the filter > >R = Sxy.*W; > >% Obtain the GCC > >G = fftshift(real(ifft(R)),1); > >figure(1); > >subplot(326); > >plot(G); > >title('HT Kernel GCC'); > > > > Hello, > > I run your program a little bit. > Actually ,time delay is wrong. > What may be the reason ? > > %%Building our signals , the same type signal > tho=2; > A=1; > t=-15:0.1:15; > SNR=20; > > x= rect( t+1,tho,A ); % First rectangle signal > > plot(t,x) > > hold on > > y=rect(t-5,tho,A); %Second rectangle signal ,delayed 5 unit 250 sample > > plot(t,y,'r') > > > > > > > %%Correlation Part > % Get power spectra of x and y and cross power spectra > Rxx = xcorr(x); > Ryy = xcorr(y); > Rxy = xcorr(x,y); > Sxx = fft(Rxx,FFTLength); > Syy = fft(Ryy,FFTLength); > Sxy = fft(Rxy,FFTLength); > > % Unfiltered Correlation (plain correlation) > W = ones(1,FFTLength); > % Apply the filter > R = Sxy.*W; > % Obtain the GCC > G = fftshift(real(ifft(R)),1); > figure; > plot(G); > %peak value of GCC give us time delay > f=max(G); > find(G==f) > disp('time delay-cross correlation') > title('Plain Time Cross Correlation'); > > time delay is 60 but the program gives values around 240. > thanks > > > > --------------------------------------- > Posted through http://www.DSPRelated.com
Try swapping the inputs around