1
$\begingroup$

I am trying to use cross correlation instead of FFT to find lags between 2 signals. I am trying to understand cross correlation by using it on sin and cos and expecting phase lag to be pi/2 but it is not:

clc
clear

frequency =1000;
t= linspace(-1000,0.1,1000);


sine_wave = sin(2*pi*frequency*t);
cos_wave  = cos(2*pi*frequency*t);

[xc,lags] = xcorr(cos_wave,sine_wave);
[~,I] = max(abs(xc));

figure
stem(lags,xc,'filled')
legend(sprintf('Maximum at lag %d',lags(I)))
title('Sample Cross-Correlation Sequence')

I am getting a maximum at lag 22 (if you copy and paste the code provided you ll see it in the plot) and I understand that phase lag = 2*pif (time delay) so in my case:

phase lag = 2*pi*1000*22 which is by no means close to pi/2... What am I missing here? Thank you

$\endgroup$
1
$\begingroup$

The following modification to your code shows what you want. Note that you need sufficiently dense sampling of the sine and cosine, far above the nyquist sampling frequency, to get the correlation result you are going for. You also need to convert lags into radians. You also want to discard integer multiples of the wavelength from the lag, since the cross correlation should be periodic as such. I did this all in the code shown below. Check that the label applied to the graph shows the $\frac{\pi}{2}$ radians lag.

frequency = 1000;
% VERY IMPORTANT
% need at least 2 samples for every 1/frequency seconds, i.e. 1 sample per
%      1/(2*frequency) seconds
% let's actually use Oversample samples for every 1/(2*frequency) seconds
Oversample = 1024;
% that means one sample every 1/frequence/100 seconds per sample
% Let's just use NumPeriods periods = time from 0 to NumPeriods/frequency
WavePeriods = 10;
SamplePeriod = 1/frequency/2/Oversample;
% start from time=0
t= (0:SamplePeriod:WavePeriods/frequency);
t2= (0:SamplePeriod:2*WavePeriods/frequency);


sine_wave = sin(2*pi*frequency*t);
cos_wave  = cos(2*pi*frequency*t2);
[xc,lags] = xcorr(cos_wave,sine_wave);
[~,I] = max(abs(xc));

figure
% set (gca,'FontName','Symbol');
stem(lags,xc,'filled');
% discard integer multiples of the wavelength
NumWavelengthsLag = mod(lags(I) * SamplePeriod * frequency,1);
legend(sprintf('Maximum at lag %d wavelengths = pi * %d radians',NumWavelengthsLag,NumWavelengthsLag*2))
title('Sample Cross-Correlation Sequence')
$\endgroup$
  • $\begingroup$ Thank you bean, could you however explain some of your choices please: Why are you choosing t2 = 2*t? why not use t for both? I am very confused about this. Thank you $\endgroup$ – user11266 Oct 9 '14 at 16:02
  • $\begingroup$ Note that bean is not doing t2=2*t. t2 just has double the duration of t. $\endgroup$ – MBaz Oct 9 '14 at 16:33
  • $\begingroup$ Thank you MBaz for the response, you are correct, however why is that needed? $\endgroup$ – user11266 Oct 9 '14 at 16:37
  • $\begingroup$ You may notice, looking at the xcorr plot, that the envelope is a diamond shape. This is because the envelope here is indicative of the length of (nonzero) overlap between the two signals being correlated. This will tend to skew the phase estimate toward zero. By correlating a short segment of sine with a longer segment of cosine, there will be a considerable section of the xcorr where the amount of nonzero overlap is constant. This helps eliminate the phase estimate skew. Note that this explanation is a bit hand-wavy, but it gets the idea across. $\endgroup$ – bean Oct 9 '14 at 19:35
  • $\begingroup$ Thank you bean, I think I get your point. However,my main issue is that I have to find phase lag between a known sine wave(like above) and a periodic wave that has the same frequency as the sine wave and is generated through the same time vector (nothing else is known about it). With the case you helped with, i already know what my phase has to be so I can tell if the answer is correct or not, however it is not that easy when I need to find out what the phase is without knowing prior answer. Do you think you can help me with this? $\endgroup$ – user11266 Oct 13 '14 at 15:58
1
$\begingroup$

To start with there is a local maxima at 2 lags. If you zoom the x-axis (say axis[-200 200 -500 500]) you can see it more clearly. However, it is not the global maxima and the one at lag 22 is higher.

This is due to the digitization of your sample. If you plot one of your waves against t and zoom the x-axis to say a 200 width box you will see that the peaks clearly have different heights. You therefore get a ringing effect in the resulting cross-correlation.

You could reduce this effect by significantly increasing your sampling rate or decreasing the signal frequency. You currently have ~5 points per $\pi$ of phase which is very low.

You could also adjust your sampling rate/frequency such that an integer number of samples fit exactly in $2\pi$ of phase.

$\endgroup$
  • $\begingroup$ Thank you nivag, it would be wonderful if you help me out with modifying my code, I will post what I think you meant here in few as well and see if it works. $\endgroup$ – user11266 Oct 7 '14 at 16:35
  • $\begingroup$ Dear nivag I tried changing the number of points but I don't see any thing good, this is what i did: frequency =1000; t= linspace(-100000,0.1,100000); however i am not getting what I need. Honestly what I expected is max lag at 1/4000 (since that time lag gives me pi/2 phase lag) is this even what it should be? please point me in the right direction. Thank you $\endgroup$ – user11266 Oct 7 '14 at 16:53
  • $\begingroup$ That hasn't changed the number of points per cycle just the number of points per cycle, just the number of cycles. I think you may misunderstand how linspace works. Either look it up or use : indexing which is much easier IMHO. $\endgroup$ – nivag Oct 8 '14 at 8:45
  • $\begingroup$ So looking at this more closely you are actually massively under-sampling your signal, you only take a point every ~1000 waves. Your sampling frequency is 1 and your wave frequency is 1000. Try setting frequency = 10 and replace t with t=1:0.001:1e5 in your original code. $\endgroup$ – nivag Oct 8 '14 at 8:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.