# Phase difference between signals sampled at different frequencies

I want to know that if it is possible to measure the relative phase difference between a signal that has been sampled at two different locations with different sampling frequencies? Also can that method be extended to undersampled cases as well?

Edit: Adding Matlab script to test possible solution (Eq.3) provided by Dan Boschen


clear all
close all
clc

Len = 768/121e6;
Fs1  = 157e6;
t1 = 0:1/(13*Fs1) :Len-1/Fs1; %Time vector for Channel 1
Fs2 = 121e6;
t2 = 0:1/(13*Fs2) :Len-1/Fs2; %Time vector for Channel 1

f=25e6; % Incoming signal frequency

phase_diff_in=0; % Modelling the actual phase difference taking In-Phase for now

% Creating signals
sign1 = cos(2*pi*f*t1);
sign2 = cos(2*pi*f*t2 + deg2rad(phase_diff_in) );
sign1 = sign1(1:13:end);
sign2 = sign2(1:13:end);

sig_ref=cos(2*pi*Fs1*t2);% Fs1 sampled by Fs2
sig_ref =sig_ref(1:13:end);

% Test of phase difference formula in time domain
phi1=acos(sign1(1:256));% In first window of 256 points
phi2=acos(sign2(1:256));
phi3=acos(sig_ref(1:256));

T1=1/Fs1;
n=0:255;
phase_diff=2*pi*n*f*( ((T1*phi3(n+1))/(2*pi*n)) -T1)...
- (phi2(n+1) - phi1(n+1));
phase_diff=wrapToPi(phase_diff);


As far as I understood the phase difference in this case should have been 0 but that is not the case. The phase difference (in deg) is as shown below:

Update: Simulating the code added by Dan

Fs1  = 157e6;
Fs2 = 121e6;
f=500e6;%25e6
samples = 400;
Len = samples;
Phi = 45;
phase_out=phase_scale(Fs1,Fs2,f,Phi,Len);
figure;
plot(phase_out)
mean(phase_out)



for the case when f=25e6 and phi=45 the following was obtained:

And for the case when f=500e6 and phi=45 the following was obtained:

The error increases significantly as the frequency is increased further.

Update #2: Simulation results after the code modifications by Dan

for the case when f=25MHz and phi=45 the following is obtained:

Which shows that the phase difference was measured very accurately.

Also for the subnyquist case as well @f=600MHz and phi=75, the following is obtained:

which shows that this works in the subnyquist cases as well. Hence the given solution works under the assumptions stated by Dan in 'Practical Limitations' section of the answer.

• You title says "Phase difference between signals" (plural), your question says phase difference between a signal (singular)--which really wouldn't make sense but wanted to ask what you are really trying to do? (Purpose?) It can help simplify the answer. Jan 9 '20 at 2:11
• Sorry for the confusion. By 'signals' I meant two different sampled versions of the same signal. The phase difference is introduced due to the sampling taking place at different locations (sensors being physically separated) Jan 9 '20 at 4:48
• ok they are physically in different locations with sampling clocks that cannot be synchronized (meaning even if you want to use two different frequencies, which is fine, you don't have the means to phase lock them to each other due to limitations of the set-up, correct?) Jan 9 '20 at 4:53
• Offset error corrected and tested over broad range of input frequencies so hopefully it is done, please let me know if it works for you! Is this a beam forming application? Jan 10 '20 at 17:43

## SOLUTION

Bottom Line

$$(\theta_2-\theta_1) = 2\pi f(T_2-T_1)n -(\phi_2[n]-\phi_1[n]) \tag{1}$$

$$f$$: frequency in Hz of two tones of the same frequency and fixed phase offset

$$(\theta_2-\theta_1)$$: phase difference in radians of tones being sampled

$$T_1$$: period of sampling clock 1 with sampling rate $$f_{s1}$$ in seconds

$$T_2$$: period of sampling clock 2 with sampling rate $$f_{s1}$$ in seconds

$$\phi_1[n]$$: phase result from sampling tone with $$f_{s1}$$ in radians/sample

$$\phi_2[n]$$: phase result from sampling tone with $$f_{s2}$$ in radians/sample

This shows how any standard approach of finding the phase between two tones of the same frequency that are sampled with the same sampling rate (common phase detectors approaches including multiplication, correlation etc) can be extended to handle the case when the two sampling rates are different.

Simpler explanation first:

Consider the exponential frequency form of equation (1):

$$e^{j(\theta_2-\theta_1)} = e^{j2\pi f(T_2-T_1)n}e^{-j(\phi_2[n]-\phi_1[n])} \tag{2}$$

The term $$e^{j2\pi f(T_2-T_1)n}$$ is the predicted difference in frequency between the two tones that would result from sampling a single tone with two different sampling rates (when observing both on the same normalized frequency scale).

The observed difference in frequency between the two tones would be $$e^{j(\phi_2[n]-\phi_1[n])}$$.

Both terms have the same frequency with a fixed phase offset. This phase offset is to the actual difference in phase between the two continuous-time tones. By conjugate multiplication we subtract the two, removing the phase slope and the fixed phase difference results.

Derivation

The approach is to carefully work with all units with a time axis of samples. The frequency domain is thus in units of normalized frequency: cycles/sample or radians/sample corresponding to cycles/sec or radians/sec when the time axis is seconds. Therefore our sampling rate, regardless of what it is in time given in seconds, will be always equal to $$1$$ cycle/sample (or $$2\pi$$ radians/sample if working in normalized radian frequency).

The two signals with the same analog frequency once sampled each with a different rate in the time domain, will be two signals each with a different normalized frequency.

This simplifies the problem to gives us the following result:

Given our original signals as normalized sinusoidal tones at the same frequency with different phase offsets:

$$x_1(t) = \cos(2\pi f t + \theta_1) \tag{3}$$ $$x_1(t) = \cos(2\pi f t + \theta_2) \tag{4}$$

Once sampled but still with time in seconds: $$x_1(nT_1) = \cos(2\pi f n T_1 + \theta_1) \tag{5}$$ $$x_2(nT_2) = \cos(2\pi f n T_2 + \theta_2) \tag{6}$$

Equation (5) and Equation (6) expressed time in units of samples is:

$$x_1[n] = \cos(2\pi f T_1 n+ \theta_1) \tag{7}$$ $$x_2[n] = \cos(2\pi f T_2 n+ \theta_2) \tag{8}$$

Convert to complex exponential form so that we can easily extract the phase terms using complex conjugate multiplication, (for a single tone we just need to split the input signal into quadrature components; $$\cos(\phi) \rightarrow [\cos(\phi),\sin(\phi)]\rightarrow \cos(\phi)+j\sin(\phi) = e^{j\phi}$$, this is described using the Hilbert Transform as $$h\{\}$$)

$$h\{x_1[n]\} =e^{-j(\phi_1[n])} = e^{2\pi f T_1 n+ \theta_1} = e^{2\pi f T_1 n}e^{\theta_1} \tag{9}$$ $$h\{x_2[n]\} = e^{-j(\phi_2[n])} =e^{2\pi f T_2 n+ \theta_2} =e^{2\pi f T_2 n}e^{\theta_2} \tag{10}$$

The complex conjugate multiplication gives us the difference phase term we seek and its relation to our measured results:

$$e^{-j(\phi_2[n]-\phi_1[n])} = e^{2\pi f T_2 n}e^{\theta_2}e^{-2\pi f T_1 n}e^{-\theta_1} \tag{11}$$

Resulting in

$$e^{j(\theta_2-\theta_1)} = e^{j2\pi f(T_2-T_1)n}e^{-j(\phi_2[n]-\phi_1[n])} \tag{12}$$

Note that $$e^{-j(\phi_2[n]-\phi_1[n])}$$ represents the measurement which for single tones will result in a frequency and this frequency is predicted to be $$\omega = 2\pi f(T_2-T_1)n$$, given by the $$e^{j2\pi f(T_2-T_1)n}$$ term. If we remove the frequency offset (by the multiplication above), the result is the phase difference of the original signal.

Taking the natural log of both sides reveals the result in units of phase (radians):

$$(\theta_2-\theta_1) = 2\pi f(T_2-T_1)n-(\phi_2[n]-\phi_1[n]) \tag{13}$$

So in summary, $$\phi_1[n]$$, $$\phi_2[n]$$ come from our measurements given as $$cos(\phi_1[n])$$, $$cos(\phi_2[n])$$ and we establish the difference that we need to get our answer through the complex conjugate multiplication of the Hilbert Transform of those measurements.

## Demonstration

I demonstrate this with the script below similar to the OP's configuration with the results plotted below, which now includes the decimation and was tested for both f = 25 MHz and f = 400 MHz (undersampled) with similar results This shows each step to demonstrate the process above, and the operations can be further combined. The Hilbert Transform in implementation would be any approach of choice to delay the sampled tones 90° (A fractional delay all-pass filter is a reasonable choice).

Len = 10000;
phase_diff_in = 45;
f=400e6; % Incoming signal frequency
D = 13
Fs1 = 157e6*D;
Fs2 = 121e6*D;
t1 = [0:Len-1]/Fs1;  % Time vector channel 1
t2 = [0:Len-1]/Fs2;  % Time vector channel 2
phi1 = 2*pi*f*t1;
sign1 = cos(phi1);
sign2 = cos(phi2);

% emulation of perfect Hilbert Transform for each tone:
c1_in = 2*(sign1 - 0.5*exp(j*phi1));
c2_in = 2*(sign2 - 0.5*exp(j*phi2));

% create expected phase slope to remove
n = [0:Len-1];
comp_in = exp(-j*2*pi*f*(1/Fs2-1/Fs1)*n);

% decimation
c1 = c1_in(1:D:end);
c2 = c2_in(1:D:end);
comp = comp_in(1:D:end);
pdout = c1.*conj(c2);
result = pdout.*comp;

%determine phase_diff

Below shows the result for the copies of the input signal at frequency $$f$$ sampled by $$f_{s1}$$ as sig1 and $$f_{s2}$$ as sig2 for the case of zero degree phase between them. The real of the complex conjugate product pdout is the bold red sinusoid, and we note that it has zero phase offset.