Cross correlation of periodic signals
Michal Kotze
Good day to all. I have a problem when cross correlating periodic signals and finding the correct delay between the two signals. In the code below i generated a 100 hertz sine wave and 0.5 sec delayed version of it. When i cross correlate the two sequences i expect to find the delay between the two sequences to be 0.5 sec but MATLAB gives in correct value of 4.000000000000000e-002. (Finding the delay is done by finding the indices of the maximum of Rxy as illustrated in the code below. I know this is to do with using the max function and the floating point representation of values. I would like to know how to correct this so that i get the correct delay when cross correlating periodic signals.
Thanks
Regards
Michal Kotze
mko...@mtnloaded.co.za
%}
clear all
t = 0:0.0001:0.9999 ; %time vector
Fs = 10000 ; % Sample Rate
CH1 = sin(2.*pi.*100.*t) ; %Change sine to randn(1,10000) then get correct delay
CH2 = zeros(1,10000) ;
CH2(7001:10000) = CH1(1:3000) ; % Set up delayed vector of CH1
[Rxy,lags] = xcorr(CH1,CH2,'coef') ;
figure (1)
subplot(3,1,1)
plot((1:length(CH1))/Fs,CH1)
subplot(3,1,2)
plot((1:length(CH2))/Fs,CH2)
subplot(3,1,3)
plot(lags/Fs,Rxy,'g')
format ('long','e')
%__Determine delay__
[Y1,I1] = max(Rxy) ;
delay = abs(lags(I1)/Fs)
%___________________
Rune Allnor
> I know this is to do with using the max function and the floating point representation of values.
No, it isn't.
The problem is that the cross correlation produces a discrete
series, and a naive search for the maximum amplitude will be
constrained to the discrete delays. If the actual delay is
not an integer number of samples, there will be poblems.
You might want to refine your delay estimates by investigating
the phase of the cross spectrum.
Rune
Michal Kotze
No your wrong look at the plot of the cross correlation function you can clearly see at -0.7 the first maximum starts and repeats. However with floating point the values of these peaks are not the same (Up until a certain amount of digits they are the same). Use format long u will see the values of the peaks are 5.477225575051663e-001 ; 5.477225575051663e-001 ; 5.477225575051664e-001 as an example, thus using the max function it will find the max value but not at the correct delay due to the fact that the values at the peaks differ. Here is the update code and could you please explain to me how to determine the delay from the phase of the cross spectrum that i dont understand. Thanks for you reply i appreciate it alot.
Regards
Michal Kotze
clear all
t = 0:0.001:0.999 ; %time vector
Fs = 1000 ; % Sample Rate
blocksize = 2000 ;
CH1 = sin(2.*pi.*10.*t) ; %Change sine to randn(1,10000) then get correct delay
CH2 = zeros(1,1000) ;
CH2(701:1000) = CH1(1:300) ; % Set up delayed vector of CH1
%CH2=CH1 ;
[Rxy,lags] = xcorr(CH1,CH2,'coef') ;
figure (1)
subplot(3,1,1)
plot((1:length(CH1))/Fs,CH1)
subplot(3,1,2)
plot((1:length(CH2))/Fs,CH2)
subplot(3,1,3)
plot(lags/Fs,Rxy,'g')
format ('long','e')
XC_results = [Rxy >= 0.5469; Rxy ; lags./Fs]'
%__Determine delay__
[Y1,I1] = max(Rxy) ;
delay = abs(lags(I1)/Fs)
%___________________
%__Cross spectrum__
xfft = abs(fft(Rxy));
index = find(xfft == 0);
mag = (xfft);
mag = mag(1:floor(blocksize/2));
f = (0:length(mag)-1)*Fs/blocksize;
f = f(:);
figure (2)
subplot(2,1,1)
plot(f,mag)
grid on
ylabel('Magnitude ')
xlabel('Frequency (Hz)')
title('Frequency Components ')
theta = unwrap(angle(fft(Rxy)) );
theta = theta(1:floor(blocksize/2));
subplot(2,1,2)
plot(f,theta)
grid on
ylabel('Magnitude ')
xlabel('Frequency (Hz)')
title('Frequency Components ')
Mark Shore
What makes you think that your CH2 is just a delayed version of CH1? Where in CH1 do you have 0.7 seconds of continuous zeros? MATLAB is doing exactly what you've asked it to, and calculated the crosscorrelation of two arbitrary vectors.
Michal Kotze
> Your code is not at all doing what you think it is.
>
> What makes you think that your CH2 is just a delayed version of CH1? Where in CH1 do you have 0.7 seconds of continuous zeros? MATLAB is doing exactly what you've asked it to, and calculated the crosscorrelation of two arbitrary vectors.
t = 0:0.001:0.999 ; %time vector
Fs = 1000 ; % Sample Rate
blocksize = 2000 ;
CH1 = sin(2.*pi.*10.*t) ; %Change sine to randn(1,1000) then get correct delay
CH2 = zeros(1,1000) ;
CH2(701:1000) = CH1(1:300) ; % Set up delayed vector of CH1
Halo Mark thanks for your reply. If t = 0:0.001:0.999 then its length is 1000 points and since it represent a time vector of 1 second the sample rate must equal the number of points generated in one second (Fs = 1000). So it is correct to say one point represents 1/1000 = 1e-3 seconds thus the 700 point will represent in time domain 0.7 seconds. So setting CH2 to zero from 1:700 and from 701:1000 = CH1(1:300) ,Channel 2 is then the delayed version of channel 1 with the delay equal to 700 points and if you want to get this value in "time domain" you must divide by the sample rate (Fs). In this case 700/Fs = 700/1000 = 0.7. Copy my code and make CH1 = randn(1,1000) you will then calculate the correct value of 0.7. The problem is with periodic signals.
Greg Heath
> %{
> Good day to all. I have a problem when cross correlating periodic signals and finding the correct delay between the two signals. In the code below i generated a 100 hertz sine wave and 0.5 sec delayed version of it. When i cross correlate the two sequences i expect to find the delay between the two sequences to be 0.5 sec but MATLAB gives in correct value of 4.000000000000000e-002. (Finding the delay is done by finding the indices of the maximum of Rxy as illustrated in the code below. I know this is to do with using the max function and the floating point representation of values. I would like to know how to correct this so that i get the correct delay when cross correlating periodic signals.
>
> Thanks
>
> Regards
>
> Michal Kotze
> mko...@mtnloaded.co.za
>
> %}
> clear all
> t = 0:0.0001:0.9999 ; %time vector
> Fs = 10000 ; % Sample Rate
>
> CH1 = sin(2.*pi.*100.*t) ; %Change sine to randn(1,10000) then get correct delay
> CH2 = zeros(1,10000) ;
> CH2(7001:10000) = CH1(1:3000) ; % Set up delayed vector of CH1
>
> [Rxy,lags] = xcorr(CH1,CH2,'coef') ;
'coeff' % ?
> figure (1)
> subplot(3,1,1)
> plot((1:length(CH1))/Fs,CH1)
> subplot(3,1,2)
> plot((1:length(CH2))/Fs,CH2)
> subplot(3,1,3)
> plot(lags/Fs,Rxy,'g')
>
> format ('long','e')
>
> %__Determine delay__
> [Y1,I1] = max(Rxy) ;
> delay = abs(lags(I1)/Fs)
> %___________________
[rmax imax0] = max(rxy) % [5.477225575051667e-001 3900]
Imax0 = find(rxy == rmax);
Npeaks0 = length(Imax0) % 40 maximum peaks
delay0 = max(t)-t(imax0) % 0.61 sec
% Exact equality doesn't work. Need a tolerance
tol(1) = 1e-17
for j = 1:18
Imax = find(rxy >= rmax-tol(j));
Npeaks(j,1) = length(Imax);
delay(j,1) = max(t)-t(min(Imax));
tol(j+1,1) = 10*tol(j);
end
whos tol Imax Npeaks delay
format short g
[(1:j)' tol(1:end-1) Npeaks delay]
% 1 1e-017 40 0.61
% 2 1e-016 49 0.69
% 3 1e-015 71 0.7
% 4 1e-014 71 0.7
% 5 1e-013 71 0.7
% 6 1e-012 71 0.7
% 7 1e-011 71 0.7
% 8 1e-010 71 0.7
% 9 1e-009 71 0.7
% 10 1e-008 71 0.7
% 11 1e-007 71 0.7
% 12 1e-006 71 0.7
% 13 1e-005 71 0.7
% 14 0.0001 71 0.7
% 15 0.001 71 0.7
% 16 0.01 497 0.7003
% 17 0.1 1479 0.7503
% 18 1 18531 0.9999
Therefore only get the desired answer when
(approx)
~eps <= tol <= ~0.001
Hope this helps.
Greg
Michal Kotze
Regards and happy programming.
Michal Kotze
clear all
t = 0:0.001:0.999 ; %time vector
Fs = 1000 ; % Sample Rate
blocksize = 2000 ;
length(t)
CH1 = sin(2.*pi.*1.*t) ; %Change sine to randn(1,10000) then get correct delay
CH2 = zeros(1,1000) ;
CH2(701:1000) = CH1(1:300) ; % Set up delayed vector of CH1
%CH2=CH1 ;
[Rxy,lags] = xcorr(CH1,CH2,'coef') ;
figure (1)
subplot(3,1,1)
plot((1:length(CH1))/Fs,CH1)
subplot(3,1,2)
plot((1:length(CH2))/Fs,CH2)
subplot(3,1,3)
plot(lags/Fs,Rxy,'g')
format ('long','e')
%XC_results = [Rxy >= 0.5469; Rxy ; lags./Fs]'
%__Determine delay__
[Y1,I1] = max(Rxy) ;
delay_incorrect = abs(lags(I1)/Fs)
%___________________
%__Section rounding of maximum value of the Rxy__
digits(5)
Y1 = vpa(Y1)
Y1 = double(Y1);
delay_correct = max(abs((lags(find(Rxy >= (Y1)))/Fs)))
%_________________________________________________
Greg Heath
> Thanks Greg you seem to be the only one understanding the problem. I appreciate it alot.I found another method to determine the correct value by rounding of the maximum value of thecrosscorrelationfunction and finding the values were thecrosscorrelationis bigger or equal to it. See the section in code {%__Section rounding of maximum value of the Rxy__}.
That doesn't work when I run your code:
>> %__Determine delay__
[Y1,I1] = max(Rxy)
delay_incorrect = abs(lags(I1)/Fs)
%__Section rounding of maximum value of the Rxy__
digits(5)
Y1 = vpa(Y1)
Y1 = double(Y1)
delay_correct = max(abs((lags(find(Rxy >= (Y1)))/Fs)))
Y1 = 6.372047745990616e-001
I1 = 363
delay_incorrect = 6.370000000000000e-001
Y1 = .63720
Y1 = 6.372000000000000e-001
delay_correct = 6.370000000000000e-001
Is this what you got?
Greg
Greg Heath
OH! .. You changed the example
from Fs = 1e4, N = 1e4, f0 = 100
to Fs = 1e3, N = 1e3, f0 = 1
In the latter case I get
1 1e-017 1 0.637
2 1e-016 1 0.637
3 1e-015 1 0.637
4 1e-014 1 0.637
5 1e-013 1 0.637
6 1e-012 1 0.637
7 1e-011 1 0.637
8 1e-010 1 0.637
9 1e-009 1 0.637
10 1e-008 1 0.637
11 1e-007 1 0.637
12 1e-006 1 0.637
13 1e-005 2 0.638
14 0.0001 5 0.639
15 0.001 18 0.646
16 0.01 57 0.665
17 0.1 181 0.727
18 1 1691 0.999
and in the former case your method yields
Y1 = 5.477225575051666e-001
I1 = 8700
delay_incorrect = 1.300000000000000e-001
Y1 = .54772
Y1 = 5.477200000000000e-001
delay_correct = 7.000000000000000e-001
Interesting that I found the first max at I1 = 3900
yielding a delay of 0.61 compared with your
8700 and 0.13
Hope this helps.
Greg
Greg Heath
WHOOPS!
I forgot to ask you why do you think we obtained
0.637 for the 1Hz example instead of the 0.7 delay
seen in the plot?
HINT:
See Mark's post.
Hope this helps.
Greg
Hope this helps.
Greg
Greg Heath
Rune,
I have not been able to get anywhere with this. If Rxy = fft(rxy),
is it a simple matter of dividing angle(Rxy) at the peak of
abs(Rxy) by 2*pi*f?
Greg
Mark Shore
k=400; m=10; n=20; p=50;
% m sine wave cycles with n samples/cycle padded with k leading and trailing zeros
CH1=[zeros(1,k) repmat(sin(2*pi/n*(1:n)),1,m) zeros(1,k)];
% shift CH1 by p sample intervals
CH2=[CH1(p:end) CH1(1:p-1)]; % no 1D equivalent of circshift in MATLAB(?)
plot(CH1,'b'); hold on; plot(CH2,'r'); hold off % optional plot
He should also consider the limiting case of a vector with a single non-zero value, say [1 zeros(1,999)] and its cross-correlation with a shifted equivalent.
Michal Kotze
Greg Heath
> On Apr 9, 7:05 am, Rune Allnor <all...@tele.ntnu.no> wrote:
>
> > On 9 apr, 12:24, "Michal Kotze" <mkotz...@yahoo.de> wrote:
>
> > > I know this is to do with using the max function and the floating point representation of values.
>
> > No, it isn't.
>
> > series, and a naive search for the maximum amplitude will be
> > constrained to the discrete delays. If the actual delay is
> > not an integer number of samples, there will be poblems.
>
> > You might want to refine your delay estimates by investigating
> > the phase of thecrossspectrum.
>
> Rune,
>
> I have not been able to get anywhere with this. If Rxy = fft(rxy),
> is it a simple matter of dividing angle(Rxy) at the peak of
> abs(Rxy) by 2*pi*f?
OK. I see it now: Ideally
g = ifft(exp(i*angle(Rxy))
should yield a pulse whos position represents the time delay.
M = length(Rxy)
t = dt*(0:2*M-1);
[ gmax jmax0] = max(g) % [ 1 1500]
delay = max(t)-t(jmax0) % 0.5 sec
Greg
Greg Heath
For the test case
f0 = 1
x0 = sin(2*pi*f0*t0) ;
x1 = [x0 zeros(1,N)];
x2 = [zeros(1,N/2) x0 zeros(1,N/2)];
t = dt*(0:2*N-1);
Greg
Mark Shore
I thought a bit more carefully about how to use circshift.
k=400; m=10; n=20; p=50; % just examples
% m sine wave cycles with n samples/cycle padded with k leading and trailing zeros
CH1=[zeros(1,k) repmat(sin(2*pi/n*(1:n)),1,m) zeros(1,k)];
CH2=circshift(CH1',p)'; % shift CH1 by p sample intervals