Cross correlation of periodic signals [Matlab]

Prev: READING TEXT BOX VALUE....
Next: Bilevel optimization

From: Greg Heath on 10 Apr 2010 00:34

On Apr 9, 6:24 am, "Michal Kotze" <mkotz...(a)yahoo.de> wrote:
> %{
> 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...(a)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

From: Michal Kotze on 10 Apr 2010 05:28

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 the cross correlation function and finding the values were the cross correlation is bigger or equal to it. See the section in code {%__Section rounding of maximum value of the Rxy__}.

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)))
%_________________________________________________

From: Greg Heath on 10 Apr 2010 10:42

On Apr 10, 5:28 am, "Michal Kotze" <mkotz...(a)yahoo.de> wrote:
> 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__}.
>
> 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)))
> %_________________________________________________

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

From: Greg Heath on 10 Apr 2010 11:40

On Apr 10, 10:42 am, Greg Heath <he...(a)alumni.brown.edu> wrote:
> On Apr 10, 5:28 am, "Michal Kotze" <mkotz...(a)yahoo.de> wrote:
>
>
>
>
>
> > 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__}.
>
> > 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)))
> > %_________________________________________________
>
> 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

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

From: Greg Heath on 10 Apr 2010 12:05

On Apr 10, 11:40 am, Greg Heath <he...(a)alumni.brown.edu> wrote:
> On Apr 10, 10:42 am, Greg Heath <he...(a)alumni.brown.edu> wrote:
>
>
>
>
>
> > On Apr 10, 5:28 am, "Michal Kotze" <mkotz...(a)yahoo.de> wrote:
>
> > > 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__}.
>
> > > 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)))
> > > %_________________________________________________
>
> > 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
>
> 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

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

First | Prev | Next | Last
Pages: 1 2 3 4
Prev: READING TEXT BOX VALUE....
Next: Bilevel optimization