Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Cross correlation of periodic signals

Subject: Cross correlation of periodic signals

From: Michal Kotze

Date: 9 Apr, 2010 10:24:05

Message: 1 of 16

%{
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
mkotze@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)
%___________________

Subject: Cross correlation of periodic signals

From: Rune Allnor

Date: 9 Apr, 2010 11:05:42

Message: 2 of 16

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.

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

Subject: Cross correlation of periodic signals

From: Michal Kotze

Date: 9 Apr, 2010 14:03:21

Message: 3 of 16

Rune Allnor <allnor@tele.ntnu.no> wrote in message <f11579f8-18c1-4640-b9ef-81976f32a77d@z4g2000yqa.googlegroups.com>...
> 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.
>
> 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


 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 ')

Subject: Cross correlation of periodic signals

From: Mark Shore

Date: 9 Apr, 2010 14:31:06

Message: 4 of 16

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.

Subject: Cross correlation of periodic signals

From: Michal Kotze

Date: 9 Apr, 2010 14:59:23

Message: 5 of 16

"Mark Shore" <mshore@magmageosciences.ca> wrote in message <hpndn9$lcd$1@fred.mathworks.com>...
> 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.

Subject: Cross correlation of periodic signals

From: Greg Heath

Date: 10 Apr, 2010 04:34:14

Message: 6 of 16

On Apr 9, 6:24 am, "Michal Kotze" <mkotz...@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...@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

Subject: Cross correlation of periodic signals

From: Michal Kotze

Date: 10 Apr, 2010 09:28:05

Message: 7 of 16

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

Subject: Cross correlation of periodic signals

From: Greg Heath

Date: 10 Apr, 2010 14:42:33

Message: 8 of 16

On Apr 10, 5:28 am, "Michal Kotze" <mkotz...@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

Subject: Cross correlation of periodic signals

From: Greg Heath

Date: 10 Apr, 2010 15:40:23

Message: 9 of 16

On Apr 10, 10:42 am, Greg Heath <he...@alumni.brown.edu> wrote:
> On Apr 10, 5:28 am, "Michal Kotze" <mkotz...@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

Subject: Cross correlation of periodic signals

From: Greg Heath

Date: 10 Apr, 2010 16:05:18

Message: 10 of 16

On Apr 10, 11:40 am, Greg Heath <he...@alumni.brown.edu> wrote:
> On Apr 10, 10:42 am, Greg Heath <he...@alumni.brown.edu> wrote:
>
>
>
>
>
> > On Apr 10, 5:28 am, "Michal Kotze" <mkotz...@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

Subject: Cross correlation of periodic signals

From: Greg Heath

Date: 11 Apr, 2010 01:19:13

Message: 11 of 16

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.
>
> 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,

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

Subject: Cross correlation of periodic signals

From: Mark Shore

Date: 11 Apr, 2010 12:40:22

Message: 12 of 16

A more appropriate signal to examine cross-correlation the way I'm guessing Michal wants would be something along the lines of :

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.

Subject: Cross correlation of periodic signals

From: Michal Kotze

Date: 12 Apr, 2010 08:04:03

Message: 13 of 16

Thanks Greg and Mark for your replies and for looking into my problem. Greg and Mark I will answer you questions soon but i am a bit busy now at university. Check the post round about Wednesday, 14 April 2010. Regards Michal Kotze

Subject: Cross correlation of periodic signals

From: Greg Heath

Date: 12 Apr, 2010 08:45:38

Message: 14 of 16

On Apr 10, 9:19 pm, Greg Heath <he...@alumni.brown.edu> wrote:
> 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.
>
> > The problem is that thecrosscorrelationproduces 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 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

Subject: Cross correlation of periodic signals

From: Greg Heath

Date: 12 Apr, 2010 08:49:52

Message: 15 of 16

On Apr 12, 4:45 am, Greg Heath <he...@alumni.brown.edu> wrote:
> On Apr 10, 9:19 pm, Greg Heath <he...@alumni.brown.edu> wrote:
>
>
>
>
>
> > 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.
>
> > > The problem is that thecrosscorrelationproduces 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 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

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

Subject: Cross correlation of periodic signals

From: Mark Shore

Date: 13 Apr, 2010 17:14:07

Message: 16 of 16

"Mark Shore" <mshore@magmageosciences.ca> wrote in message <hpsfvm$rtp$1@fred.mathworks.com>...
> A more appropriate signal to examine cross-correlation the way I'm guessing Michal wants would be something along the lines of :
>
> 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.

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
 plot(CH1,'b'); hold on; plot(CH2,'r'); hold off % optional plot

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us