# Computing phase difference between two sinusoidal signals

Based on this post here, I tried to calculate the phase difference between sine waves. Here is my code.

Fs = 64*10**6
N = 2048
f = 100000
t = np.linspace(0, (N-1)/Fs,N)
order = 3
amp = (2**(12)-1)/2
def butter_lowpass(cutoff, fs, order):
return butter(order, cutoff, fs=fs, btype='low', analog=False)

def butter_lowpass_filter(data, cutoff, fs, order):
b, a = butter_lowpass(cutoff, fs, order=order)
y = lfilter(b, a, data)
return y

noise = np.random.random(len(t))
y = amp * (np.sin( f * 2 * np.pi* t) +1 )
# y = y - np.mean(y)
plt.plot(t,y, label = 'y')

z = amp * (np.sin(f * 2 * np.pi * t + np.pi/4) + 1)
# z = z - np.mean(z)
plt.plot(t, z, label = 'z')

x  = y * z
# print(max(x))
x /= np.max(x)
# x = 2 * amp * x
x*= amp

# x = x - np.mean(x)
plt.plot(t,x, label = 'x')
filt = butter_lowpass_filter(x,f , Fs, order)

# filt = filt/np.max(filt)
# plt.plot(t, filt, label = 'filt')
phi = np.arccos(filt)
phi = np.degrees(phi)
# plt.plot(phi, 'ro') #to degrees
print(phi)
# # # plt.show()
# plt.legend()
# i = np.sum(x)


But the output phi is wrong. Where do I go wrong? I assume we need to take inverse cosine of the filtered signal while I am not sure regarding taking the mean of the array to calculate phi.

I don't want to use alternative methods like FFT, cross correlation because I want to implement this in STM32 microcontroller later. I'd appreciate if anyone could help me with this.

• The way You are generating sine tone in code does not look appropriate. Commented Apr 3, 2023 at 12:35
• np.sin((2 * np.pi * t * f)/Fs) will give You a tone of frequency "f". Your code is missing "Fs". Commented Apr 3, 2023 at 12:38
• I corrected my code. For now, I am just not taking in account of Fs in my code. But the phase calculation is wrong.
– Rima
Commented Apr 3, 2023 at 12:42
• Even the updated code doesn't look correct. Signal y in the code still looks to be an array of zeros Commented Apr 3, 2023 at 12:47
• t = np.linspace(0, (N-1)/Fs,N) I agree. I finally decided to use higher f so that I could see in the plot with the updated code. Even though my target is f = 1000-10k.
– Rima
Commented Apr 3, 2023 at 12:49

The phase between the sine waves can be computed as suggested, at any sampling rate as long as it is greater than twice the sum of the frequency of the sine waves. Further the result is sensitive to the peak amplitude so the waveforms need to be normalized or the phase result accordingly scaled.

This is all clear from the cosine product rule that the product of two sinusoids will have a low frequency component that is proportional to the cosine of the phase between them:

$$\cos(\alpha)\cos(\beta)= \frac{1}{2}(\alpha + \beta)+ \frac{1}{2}(\alpha - \beta)$$

So when the two frequencies are the same and there is a phase difference we have:

$$2\cos(\omega_1t+ \theta)\cos(\omega_1t) = \cos(2\omega_1t+ \theta)+ \cos(\theta)$$

After a low pass filter we get $$\cos(\theta)$$. Note that the time delay of the low pass filter has no impact on the result (the answer for a static phase offset is the DC mean itself). The loss of the filter will change the magnitude and this must be calibrated out accordingly.

If the system is sampled and the product is done in discrete time, the sampling rate must be high enough to support the double frequency component above as otherwise it can alias to be within the bandwidth of the subsequent low pass filter. Sampling at an even higher frequency will only provide more samples in the determining the mean (which is the purpose of the low pass filter), so will not be detrimental to the result other than needlessly increasing filter complexity to achieve the same cutoff frequency.

Here is a simulation using the OP's parameters demonstrating in simplest form the desired result and operation:

Python code:

Fs = 64*10**6
N = 2048
f = 100000
t = np.linspace(0, (N-1)/Fs,N)
amp = 1
phase = np.pi/4
y = amp * (np.sin( f * 2 * np.pi* t) )
z = amp * (np.sin(f * 2 * np.pi * t+ phase) )
x  = y * z
print(f"Actual phase : {phase:0.3f}")
print(f"Measured phase: {np.arccos(2*np.mean(x)):0.3f}")


NOTES on the results:

• y and x are the inputs as sinusoids with a 0.785 radian offset in phase between them. z is the result of the product $$xy$$.

• There should not be any DC offsets on the sinusoidal waveforms; given the result for a static phase offset is a DC term, this will only be an interfering offset (error).

• The average of the product is the desired result that will be proportional to the cosine of the phase angle. Here with the waveforms each normalized to 1, the solution is $$\phi = \overline{2xy}$$, meaning twice the average of the product. The purpose of any low pass filter is to remove the higher (double) frequency component and pass the DC value. If the phase is expected to change with time within a certain bandwidth, then the low pass filter used must be designed to pass that bandwidth.

• Note the measured result is 0.755, which is in error from the actual phase offset of 0.785. This is because the output has not gone through an integer number of cycles, so the final cycle is incompletely contributing to the mean. This however does not mean we should always attempt to average over an integer number of cycles of the double frequency output, but instead illustrates that the averaging duration in this case is not long enough for the frequency of the sine wave used. As the duration is increased (or similarly as the frequency of the sine wave is increased), the fractional error from any partial cycles is reduced. For applications where this can't be avoided, windowing techniques could be used to reduce the effects of this similar to that done for FFT processing.

Given the target application requires high efficiency, consider using a CIC filter instead of the IIR filters suggested in the OP's code which will combine performing a moving average together with a decimated (down-sampled output). The implementation to do a moving average over all 2048 samples is shown below:

The $$z^{-1}$$ blocks are one sample delays at the input rate on the left of the decimator, and output rate at the right of the decimator (running at 1/2048 of the input rate). The diagram to the left of the decimator just after the product represents and accumulator, and this accumulator must be designed such that it wraps on an overflow, and no more than one wrap is allowed during the period of a decimation cycle. The entire diagram shown must be implemented in fixed point arithmetic. The decimator indicates selecting every 2048th sample from the accumulator output and throwing away the rest. Every subsequent sample at the decimator output is subtracted from the previous. Alternatively this could be implemented with an "integrate and dump" where an accumulator grows over 2048 samples, dumps the result and clears. The result is equivalent to a moving average, exactly so if the output is divided by 2048 (convenient bit shift operation given that is a power of 2, but other scaling may be needed depending on input level). What is convenient with the CIC implementation shown, is that it is completely programmable in terms of the moving average number of sample duration just by changing the decimation: change 2048 above to 512 for example to then be a moving average over 512 samples (and scale accordingly).

Finally for ultimate efficiency a Cordic Rotator could be considered for efficient computation of the inverse cosine.

HOWEVER, if this is to be used in a Phase Lock Loop, then the loop filter itself will properly serve as the filter, and there is no reason for an inverse cosine computation: The loop will drive the phase difference to 90°, where $$\cos(\theta)=0$$, since that is the measure of loop error. The error in phase from locked condition will then be $$\sin(\theta)$$ and for small angles $$\sin(\theta) \approx \theta$$. Further any gain will be part of the loop gain itself, so the only concern is for the input sinusoids to be normalized to some constant value (as done with an AGC or automatic gain control). Everything gets simpler in terms of what else is needed to measure phase when the multiplier is used as the phase detector in a loop!

• Thanks for the detailed answer. I checked out changing some parameters and for N = Fs/f(just one period), for the range of phase, there is significant error between the expected and the calculated phase. This will heavily lead to error for phase correction. Regarding your last paragraph, I want to make a PLL and this phase comparator is a part of it(that's what I thought). You mean, i just multiply in phase detector part and need not calculate phase but rather use loop filter for this multiplier to drive the phase difference to be 0? can we still use this mean filter as a loop filter?
– Rima
Commented Apr 5, 2023 at 7:46
• @Rima glad this helped! When used in a PLL the scaling error does not matter nor do you need to compute the inverse cosine since for small angles the error around 0 (which occurs at 90 degree offset) will be linear. I won’t be able to detail a full PLL design in the comments here but you will want to use a “proportional integral loop filter” for a stable solution that drives to zero error and not any arbitrary low pass filter Commented Apr 5, 2023 at 12:01
• I checked out some online material and found the reference for PI filter. Also read about some other Phase detectors. This multiplier, as you said, will be linear around 90 degree offset but for small angles is nonlinear. Btw, I see that greater the ratio between Fs and f is, the lesser the deviation would be; however, my f could go as high as 10KHz so that causes a lot of deviation. If arcos is not required in PLL then should I just use a multiplier and then PI filter? I want to make sure to choose a more precise Phase detector, have a stable PLL for a wide range of phase difference.
– Rima
Commented Apr 5, 2023 at 12:22
• @Rima for small angles close to 90 degrees it is linear. The loop will lock to where the zero crossing is for the error output (your phase detector output). The multiplier based PLL will have a lock range for maximum frequency offset-- beyond that you need to incorporate a helping FLL to pull-in or Phase/Frequency detector (or initial frequency search for acquisition). We won't be able to get into a lot of the PLL design issues in the comments since back and forth discussion is discouraged and it will flag the moderators. Feel free to post a new question if your stuck. Commented Apr 5, 2023 at 12:29