Cross correlation of periodic signals
in [Matlab]
|
Prev: READING TEXT BOX VALUE....
Next: Bilevel optimization
From: Michal Kotze on 9 Apr 2010 06:24 %{ 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(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') ; 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) %___________________
From: Rune Allnor on 9 Apr 2010 07:05 On 9 apr, 12:24, "Michal Kotze" <mkotz...(a)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
From: Michal Kotze on 9 Apr 2010 10:03 Rune Allnor <allnor(a)tele.ntnu.no> wrote in message <f11579f8-18c1-4640-b9ef-81976f32a77d(a)z4g2000yqa.googlegroups.com>... > On 9 apr, 12:24, "Michal Kotze" <mkotz...(a)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 ')
From: Mark Shore on 9 Apr 2010 10:31 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.
From: Michal Kotze on 9 Apr 2010 10:59
"Mark Shore" <mshore(a)magmageosciences.ca> wrote in message <hpndn9$lcd$1(a)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. |