2
$\begingroup$

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

% Adding a reference cosine
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);
figure;plot(rad2deg(phase_diff),'-*r')

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:

enter image description here

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:

enter image description here

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

enter image description here

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:

enter image description here

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:

enter image description here

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Dan Boschen Jan 9 at 2:11
  • $\begingroup$ 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) $\endgroup$ – malik12 Jan 9 at 4:48
  • $\begingroup$ 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?) $\endgroup$ – Dan Boschen Jan 9 at 4:53
  • $\begingroup$ 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? $\endgroup$ – Dan Boschen Jan 10 at 17:43
4
$\begingroup$

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;   
phi2 = 2*pi*f*t2 + deg2rad(phase_diff_in);
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
phase_out = rad2deg(unwrap(angle(result)));
mean_phase = mean(phase_out); 

Below is the result for two test cases, 0° as the OP was trying in his example and then a 45° phase shift.

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.

Phase Detector Output

To confirm the calculations, the plot below compares it directly to the real of the compensation term $cos(2\pi f(T_2-T_1)) to see that they are the same frequency consistent with the equation.

With angle diff = 0

And repeating with $\theta_2-\theta_1 = 45°$

With angle diff = 45

The result of the raw phase data for every sample shows that each sample individually has extremely low noise (limited by numerical precision, so the result can be determined with very few samples!). Such performance will depend on the actual quality of the Hilbert transform to accurately delay the input tone by 90° to create a qaudrature copy. Under conditions of noise the result can be averaged to the extent of waveform stability for a very robust solution.

Phase vs Sample

Extended testing of performance with the undersampling case shows excellent results (f = 400e6):

Testing every difference angle in 1 degree steps: phase sweep

RMS error of 10,000 samples (Note vertical axis is in increments of 0.5e-11)

Std deviation vs Input phase

Result of a greatly extended frequency sweep of the input frequency from 1e6 to 4000e6 in steps of 1e6 with a 45 degree phase shift with 10,000 points measured at each frequency showed a consistent result for phase determination at all frequencies (oversampling and undersampling). This is with the OP's configuration with the two frequencies including the decimate by 13. (Thus the sampling rate of each of the input tones after decimation for this test was at fs = 157e6 and 121e6, thus the far right of this graph with the frequency of the tone being sampled being 4e9 is significantly under-sampled. The RMS error is proportional to the frequency of the tone as shown, but even under this extreme condition, the error is still less than 5e-10 degrees. (8.7e-12 radians or -221 dB).

extended frequency sweep


Practical Limitations

The accuracy of the above result is limited by knowledge of the exact frequencies and phase relationship given by $f_{s1}$ and $f_{s2}$, and knowledge of the frequency $f$ of the tone being sampled.

(As written the solution also assumes that the two sampling clocks both start at time $t=0$, but the sampling offset can be added starting with equation (8) with a similar result; bottom line is the starting phase relationship between the two sampling clocks must be known or measured as it will introduce an additional offset).

The reality is that no two free-running clocks will stay in perfect synchronization; there will be an inevitable drift in the actual frequency and phase difference between the sampling clocks that are not locked to a common reference (see Segal's Law https://en.wikipedia.org/wiki/Segal%27s_law). One of the clocks must be declared our reference of time (and our measurement will be to the accuracy of that one clock). If the clocks are not physically co-located, two-way time transfer techniques (see https://tf.nist.gov/time/twoway.htm) can be used to measure one clock versus the other. If they are physically co-located, then the simple thing to do would be to sample one clock with the other.

Below I show how this approach can completely eliminate one of the sampling clocks from the equation for our solution: (I haven't yet tested this so may contain math errors)

Consider sampling $f_{s1}=\frac{1}{T_1}$ with $f_{s2}=\frac{1}{T_2}$. This will ultimately remove $f_{s2}$ from the equation entirely by using $f_{s1}$ as the common reference (we essentially measured $f_{s2}$ with $f_{s1}$ by sampling $f_{s1}$ with $f_{s2}$ allowing us to put the samples of $f_{s2}$ in units of $f_{s1}$ counts.):

$f_{s1}$ as a cosine:

$$x_{s1}(t) = cos(2\pi f_{s1}t) \tag{14}$$

When sampled with $f_{s2}$ given the constraint they both start at t=0 becomes:

$$x_{s_1}(nT_2) = cos(2\pi f_{s1}nT_2) = cos(2\pi nT_2/T_1) \tag{15}$$

Which in units of samples is:

$$x_{s_1}[n] =cos(2\pi T_2/T_1 n) \tag{16}$$

Resulting in a third phase measurement in units of samples that we can get by sampling $f_{s1}$ with $f_{s2}$ (importantly to be done at the same time $x_1(t)$ and $x_2(t)$ are sampled!):

$$\phi_3[n] = 2\pi T_2/T_1 n \tag{17}$$

Thus if we don't know $T_2$ but have $\phi_3$ we can substitute the above equation to get:

$$T_2 = \frac{T_1 \phi_3[n]}{2\pi n} \tag{18}$$

substituting into (4):

$$ \phi_2[n]- \phi_1[n] = 2\pi nf\bigg(\frac{T_1 \phi_3[n]}{2\pi n}-T_1\bigg) + (\theta_2-\theta_1) \tag{19} $$

Resulting in the following solution for the desires phase difference of the original input signals:

$$ \theta_2-\theta_1= 2\pi f\bigg(\frac{T_1 \phi_3[n]}{2\pi n}-T_1\bigg)n - (\phi_2[n]- \phi_1[n]) \tag{20}

$$

Where

$f$: frequency of tone being sampled

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

$\phi_1[n]$: result from sampling tone with $f_{s1}$, values will be $cos(\phi_1[n]) $

$\phi_2[n]$: result from sampling tone with $f_{s2}$, values will be $cos(\phi_2[n]) $

$\phi_3[n]$: result of sampling $f_{s1}$ with $f_{s2}$, values will be $cos(\phi_3[n]) $

Thus by only knowing $T_1$ which is $1/f_{s1}$, we can measure $f$ from the samples of $x_1(t)$ directly, measure $\phi_1[n]$ by sampling $x_1(t)$ with $f_{s1}$, measure $\phi_2[n]$ by sampling $x_2(t)$ with $f_{s_2}$ and measure $\phi_3[n]$ by sampling $f_{s1}$ with $f_{s2}$ and from those measurements resolve $\theta_2-\theta_1$.

Similarly if your application is for a phase offset that would not be changing, then you can measure $f_{s2}$ error using the slope of the result without having to sample $f_{s1}$ with $f_{s2}$.

The true results will depend on the actual clock accuracy of $f_{s1}$ but we have completely removed $f_{s2}$ from the equation. If you can consider $f_{s1}$ to be your true reference of time, meaning it is accurate enough for the precision and accuracy of your measurement, then the result will be the phase difference of your two waveforms being sampled. This means that ultimately you need something to be your common reference of time.

|improve this answer|||||
$\endgroup$
  • $\begingroup$ The sampling clocks are not synchronized. Moreover, the clock periods have non-integer relation hence it will be difficult to "know the time difference between the sampling clocks at a specific absolute time". e.g. how can it be done if the sampling frequencies are 121 MHz and 157 MHz. $\endgroup$ – malik12 Jan 8 at 5:08
  • $\begingroup$ At t= 0 the phase between the two sampling clocks is 0. The phase between them for all times after that is a linear ramp versus time---so as long as you can anchor the start of that ramp you have the information you need. $\endgroup$ – Dan Boschen Jan 8 at 5:12
  • $\begingroup$ In practical reality- you will need to have the two sampling clock sources locked to a common reference---- (such as a 10 MHz clock), and that clock becomes your absolute reference for all your phase answers (errror in that clock will be proportional errors in your phase as your frequencies as given in my equations will actually be in error). Just consider $f_2$ and $f_1$ in my equation and if your frequencies were off by 1ppm (for example) what would happen to the answer. $\endgroup$ – Dan Boschen Jan 8 at 5:17
  • $\begingroup$ Kindly elaborate what you mean by " locked to a common reference---- (such as a 10 MHz clock)" $\endgroup$ – malik12 Jan 8 at 10:21
  • $\begingroup$ What I mean is the two sampling clocks have the same reference for their time. For example signal generators that create arbitrary tones have the option of providing an external 10MHz input so that they can be run off of a common reference. This way they will both have the same percentage error in their actual frequencies. If they weren’t tied to a common reference they would have independent and unknown errors. This is the practical part as clocks drift over time, so if the two were not tied together you and you have no other measurement you wouldn’t know what that drift is. $\endgroup$ – Dan Boschen Jan 8 at 13:00
1
$\begingroup$

For the first part of your question, perhaps this will shed some light:

Phase difference measurement of a signal sampled with two different sampling frequencies

The answer to your second part of your question is yes for a single pure tone. It will appear as a lower frequency alias in the DFT, but if you know the actual frequency range, you can calculate the correct frequency.

One caveat to that. If it is an alias of the DC or Nyquist frequencies, it may or may not show up. Those are the DFT's potential "blind spots".


I believe the third part of my linked answer is the most efficient and the most accurate, especially considering the possibility of being undersampled.

It goes like this:

Find values of $M$ and $N$ such that:

$$ \frac{M}{N} = \frac{T_1}{T_2} - \epsilon $$

That gets you this:

$$ (MT_2 \approx NT_1 )= T_{DFT \; frame} $$

Since you know $f$, you know $k_1$ in a $M$ sample DFT on the first signal, and $k_2$ in a $N$ sample DFT on the second signal for the same time duration. You only have to calculate two bin values in each DFT, the $k$ and $k+1$ which bookend where $f$ falls in each respectively. Use the phase calculation of my two bin solution to solve for relative parameters (not the prior article version which doesn't unfurl the two bins into a real vector). You already know the frequency, so you don't have to estimate it.

To keep the variable names distinct, let's call them:

$$ \begin{aligned} S_1[n] &= A_1 \cos( \omega_1 n + \tau_1 ) \\ S_2[m] &= A_2 \cos( \omega_2 m + \tau_2 ) \\ \end{aligned} $$

The latter part of the two bin solution will solve for the $A$ and $\tau$ parameters. The $\omega$s are known ahead of time.

Using the sampling rate (in samples per second), these can be converted to real world values and the phase values directly compared for that interval. If $\epsilon$ is large, it can be incorporated into this calculation.

This solution renders magnitude differences irrelevant. It also allows you to compensate for possible aliasing in the case of an undersampled signal.

There is an inherent latency of the DFT frame duration.

[The bolding is for the benefit of the OP and others, not Dan]


What I glossed over since it was much discussed is that any time differential of the arrivals of the signal will translate directly into an error in the phase difference. If the phase difference is across many samples, this is just a source of a little inaccuracy. If the phase difference is sample sized or even subsample (detectable by the DFT method) it presents a real problem. One solution for calibration, which may or may not be appropriate for the OP, is addressed in the first link.


It is also possible to select the duration of a whole number of cycles, and select $M$ and $N$ from that. Only one DFT bin each needs to be calculated then, and the basis vectors can be predefined.

Selecting a whole number of cycles plus a half, where $MT_2$ is very close to $NT_1$ and employing the two bin phase is more noise resistant, but it requires two DFT bin calculations per signal.


Reply to Dan's comment:

I'm not really good on the Hilbert. My understanding is that the discrete version is just an approximation of the continuous. Discrete differentiation compared to derivatives is a whole 'nuther topic.

Shall we set up some sample data and have a contest?

The frequency of the signal tone is known a priori, no need to estimate it.


Okay, this took a lot longer than it should have. Debugging prints left as comments.

The results:

 1     6     5   0.833333   0.770701   0.062633   0.955414   1.033058
 2    13    10   0.769231   0.770701   0.001470   2.070064   2.066116
 3    19    15   0.789474   0.770701   0.018773   3.025478   3.099174
 4    25    19   0.760000   0.770701   0.010701   3.980892   3.925620
 5    31    24   0.774194   0.770701   0.003493   4.936306   4.958678
 6    38    29   0.763158   0.770701   0.007543   6.050955   5.991736
 7    44    34   0.772727   0.770701   0.002027   7.006369   7.024793
 8    50    39   0.780000   0.770701   0.009299   7.961783   8.057851
 9    57    44   0.771930   0.770701   0.001229   9.076433   9.090909

Too much time spent on this already, the code will have to speak for itself.

[redacted - look below]


Comment on Dan's solution that is too long for a comment:

Instead of doing a Hilbert, you can get a sine from a cosine by shifting the signal by a quarter cycle. This retains the proper amplitude as well so you get a circular spiral. Then you can frequency shift it:

$$ A_1 e^{i (\omega_1 n + \phi_1) } \cdot e^{i \Delta \omega n } = A_1 e^{i [(\omega_1 + \Delta \omega ) n + \phi_1 ]} $$

I did this last summer on a FMCW project.

Shift the other signal in the other direction to generate the "conjugate of what it would have been" and multiply.

$$ A_1 e^{i [(\omega_1 + \Delta \omega ) n + \phi_1 ]} A_2 e^{-i (\omega_2 n + \phi_2) } = A_1 A_2 e^{i [(\omega_1 - \omega_2 + \Delta \omega ) n + ( \phi_1 - \phi_2)] } $$

Let $ \Delta \omega = \omega_2 - \omega_1 $ and you get $ A_1 A_2 e^{i ( \phi_1 - \phi_2 )} $.

You can now read the phase difference directly from the arg. You are "stretching time" on the signals to get them to match sample by sample, so your samples aren't matched in actual time. So, to get the best reading for a particular time interval, I would select samples so the shorter interval is centered in the longer interval. Then you will want to average the point by point $\Delta \phi$ readings to get a single value for $\phi_1 - \phi_2$. It might be advantageous to just use a shorter interior centered interval within the result interval for the averaging.


It turned out to be necessary to apply the off-bin phase adjustment in order to get the results reasonably accurate for a small number of cycles per frame. Increasing the cpf will still improve accuracy, but at the cost of latency. Overlapping sliding windows are no problem for steps sizes shorter than the latency.

The phase adjustment formula and derivation can be found here:

(One of my best answers ever, BTW, yet no upvotes.)

Here are the adjusted results:

Selected 9 57 44

28.5 28.5 -1.88182802674e-14
22.0 22.0 -1.60982338571e-15

Omegas 1.00050721452 1.29817878248

Peak and Fs 9.0 9.07643312102 9.09090909091

 0  0.2006  1.1933  0.9927
 1  0.4977  1.4948  0.9971
 2  0.7956  1.7982  1.0026
 3  1.0950  2.1022  1.0072
 4  1.3962  2.4054  1.0093
 5  1.6987  2.7067  1.0081
 6  2.0017 -3.2775 -5.2792
 7  2.3041 -2.9805 -5.2846
 8  2.6050 -2.6845 -5.2895
 9  2.9042 -2.3880 -5.2922

The $2\pi$ adjustment was purposefully not applied.

Here is the new code. It should be easy for any one to slip in a test algorithm of their own. The units comments should be instructional, even for non-programmer.

import numpy as np

#====================================================================
def main():

#---- Set Parameters

        Fs1 = 157e6
        Fs2 = 121e6
        f   =  25e6 # Incoming signal frequency

#---- Calculate Derived Values

                          # = samples per second / cycles per second 
        theSamplesPerCycle1 = Fs1 / f
        theSamplesPerCycle2 = Fs2 / f

#---- Display M and N combinations

        Q_21 = Fs2 / Fs1

        for cpf in range( 1, 10 ):
          N = int( theSamplesPerCycle1 * cpf + 0.5 )
          M = int( theSamplesPerCycle2 * cpf + 0.5 )

          Q_MN = float( M ) / float( N )

          E = abs( Q_MN - Q_21 )

          k1 = float( N ) / theSamplesPerCycle1
          k2 = float( M ) / theSamplesPerCycle2

          print "%2d %5d %5d %10.6f %10.6f %10.6f %10.6f %10.6f" %\
                 ( cpf,  N,  M,  Q_MN,  Q_21,     E,    k1,    k2 )


        print 

#---- Determine the DFT sizes

        theCyclesPerFrame = 9

        N = int( theSamplesPerCycle1 * theCyclesPerFrame + 0.5 )
        M = int( theSamplesPerCycle2 * theCyclesPerFrame + 0.5 )

        print "Selected", theCyclesPerFrame, N, M
        print

#---- Build the Basis Vectors of the DFT bin

        C_N, S_N = BuildDftVectors( theCyclesPerFrame, N )
        C_M, S_M = BuildDftVectors( theCyclesPerFrame, M )

        print C_N.dot( C_N ), S_N.dot( S_N ), C_N.dot( S_N )
        print C_M.dot( C_M ), S_M.dot( S_M ), C_M.dot( S_M )
        print

#---- Calculate Normalized Frequencies

        # radians per sample = radians per cycle 
        #                    / samples per cycle

        omega1 = 2.0 * np.pi / theSamplesPerCycle1
        omega2 = 2.0 * np.pi / theSamplesPerCycle2

        print "Omegas", omega1, omega2
        print

#---- Set the Adjustment Parameters

        # cycles per frame = samples per frame
        #                  / samples per cycle

        f1 = N / theSamplesPerCycle1
        p1 = float( theCyclesPerFrame )

        f2 = M / theSamplesPerCycle2
        p2 = float( theCyclesPerFrame )

        print "Peak and Fs", p1, f1, f2
        print

#---- Do Some Runs

        for theTestRun in range( 10 ):
          theSignal1 = BuildSignal( 1000, 1.1, omega1, 0.2 + 0.3 * theTestRun )
          theSignal2 = BuildSignal( 1000, 1.2, omega2, 1.2 + 0.3 * theTestRun )

          RunTest_Cedron( theTestRun, theSignal1, theSignal2,  \
                          omega1, omega2, Fs1, Fs2,            \
                          f1, p1, f2, p2,                      \
                          C_N, S_N, C_M, S_M )

#====================================================================
def BuildSignal( argSampleCount, argAmplitude, argOmega, argPhi ):

        x = np.zeros( argSampleCount )

        for n in range( argSampleCount ):
          x[n] = argAmplitude * np.cos( argOmega * n + argPhi )

        return x

#====================================================================
def RunTest_Cedron( argTestRun, argSignal1, argSignal2, \
                    omega1, omega2, Fs1, Fs2,           \
                    f1, p1, f2, p2,                     \
                    C_N, S_N, C_M, S_M ):

        theInterval1 = argSignal1[0:len( C_N )]
        theInterval2 = argSignal2[0:len( C_M )]

        thePhase1 = FindPhaseOf( theInterval1, C_N, S_N, f1, p1 )
        thePhase2 = FindPhaseOf( theInterval2, C_M, S_M, f2, p2 )

        theDeltaPhase = thePhase2 - thePhase1

        # samples = radians   / radians per sample
#        theShift1 = thePhase1 / omega1        
#        theShift2 = thePhase2 / omega2

        # seconds = samples / samples per second
        # theDelta1 = ( theShift1 / Fs1 ) * 1000000.0
#        theDelta2 = ( theShift2 / Fs2 ) * 1000000.0
#        theDiff   = theDelta1 - theDelta2

        print "%2d %7.4f %7.4f %7.4f" %\
              ( argTestRun, thePhase1, thePhase2, theDeltaPhase )

        return theDeltaPhase

#====================================================================
def FindPhaseOf( argInterval, C, S, f, p ):

#---- Calculate the DFT Bin Value

        real = argInterval.dot( C )
        imag = argInterval.dot( S )

        theBinPhase = np.arctan2( imag, real )

#---- Apply Off-bin Phase Approximation

        MN = float( len( C ) )

        theDeltaPhase = -( f - p ) * ( MN - 1.0 ) / MN * np.pi 

#---- Return the Bin's angle

        return theBinPhase + theDeltaPhase

#====================================================================
def BuildDftVectors( argCyclesPerFrame, argSamplesPerFrame ):

        C = np.zeros( argSamplesPerFrame )
        S = np.zeros( argSamplesPerFrame )

        theSlice = 2.0 * np.pi / float( argSamplesPerFrame )

        theStep  = argCyclesPerFrame * theSlice
        theAngle = 0.0

        for n in range( argSamplesPerFrame ):
          C[n] =  np.cos( theAngle )
          S[n] = -np.sin( theAngle )
          theAngle += theStep

        return C, S

#====================================================================
main()
|improve this answer|||||
$\endgroup$
  • $\begingroup$ Thank you for your detailed input on this matter as well. I will also simulate your solution and update accordingly. $\endgroup$ – malik12 Jan 10 at 10:08
  • $\begingroup$ @malik12 You're welcome. Dan B. also needed benchmark for speed and accuracy. Obviously longer frames (a larger cycle count) would be more accurate here. A choice of 9 would have been better as they are off integer by about the same amount so if you are taking their difference that error will largely wash out. Otherwise you know $f$ and $p$ and can apply the phase adjustment provided in your other question here: dsp.stackexchange.com/questions/63076/… $\endgroup$ – Cedron Dawg Jan 10 at 14:11
  • $\begingroup$ It would be nice to see the actual math that does the Hilbert transform. No matter what, it is going to be more computationally expensive than merely shifting your signal by a fixed number of samples. However, the shift trick will not work well in signals with a small number of samples per cycle. Note, you can skip the $\Delta\omega$ multiply step and read the intercept off the line formed by the arg of the products $(\omega_1 - \omega_2) n + ( \phi_1 - \phi_2)$. The product is actually just a bunch of dot products of the signals and their shifts. $\endgroup$ – Cedron Dawg Jan 10 at 15:54
  • $\begingroup$ @DanBoschen In this case, the frequency is known and that fact should be exploited if possible. Two metrics at play: Accuracy and Efficiency. $\endgroup$ – Cedron Dawg Jan 10 at 15:55
  • $\begingroup$ Yes of course- what I mean is for the case of a known single tone, the Hilbert Transform is just a delay by a quarter cycle- not instead of. The math is well documented if you are interested. Given it is a single tone I would likely implement with a delay line rather than the broadband filter solution (that is what you would do to create an analytic signal for a modulated waveform that occupies multiple frequencies- genetically for this single tone case “Hilbert” means create a quadrature tone from the real tone, with many ways to do that well established. $\endgroup$ – Dan Boschen Jan 10 at 16:09

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.