Theoretical background

In [Reference Kuramoto and Araki11], Kuramoto described the coupling between a set of harmonic oscillators with

(1)\begin{equation}\varphi'_m = \omega_{m,0} + \varepsilon_m\sum_{n=1}^M g(\psi_{nm}), \end{equation}

where

(2)\begin{equation} \psi_{nm} = \varphi_n - \varphi_m. \end{equation}

φ_m is the phase of harmonic oscillator m and $\varphi'_m$ is its instantaneous frequency. Moreover, ɛ_m describes the coupling strength between oscillators, $\omega_{m,0}$ is the initial angular frequency of oscillator m, and ψ_nm agrees with the phase difference between oscillators n and m. Equation (1) describes the generalized version of the Kuramoto model, where g(x) is an arbitrary 2π-periodic function. The standard model as proposed by Kuramoto approximates g(x) with

(3)\begin{equation} g_\mathrm{S}(x) = \sin x. \end{equation}

For two oscillators, it simplifies to

(4)\begin{align} \varphi'_1 &= \omega_{1,0} + \varepsilon_1\ g(\psi_{21}), \end{align}

(5)\begin{align} \varphi'_2 &= \omega_{2,0} + \varepsilon_2\ g(\psi_{12}). \end{align}

For $g(x)=g_\text{S}(x)$ and by setting $\psi'_{nm} = 0$, it follows:

(6)\begin{equation} \psi'_{nm} = \arcsin\left(\frac{\omega_{n,0}-\omega_{m,0}}{\varepsilon_{1}+\varepsilon_{2}}\right). \end{equation}

Since the argument of the arcsine must be $\in[-1, 1]$, the following synchronization criteria can be derived [Reference Dallmann14]:

(7)\begin{align}\varepsilon_{1} &= \gamma_{1}\left|\omega_{2,0}-\omega_{1,0}\right| =2\pi\gamma_{1}\left|\,f_{2,0}-f_{1,0}\right|,\; \gamma_{1} \geq 1, \end{align}

(8)\begin{align}\varepsilon_{2} &= \gamma_{2}\left|\omega_{1,0}-\omega_{2,0}\right|=2\pi\gamma_{2}\left|\,f_{1,0}-f_{2,0}\right|,\; \gamma_{2} \geq 1. \end{align}

This means that synchronization can be successful if the coupling factor ɛ_l is $\gamma_{l} \geq 1$ times larger than the initial angular frequency difference between the oscillators (for $l\in\{1, 2\}$). In principle, these equations could be directly used for frequency alignment of CW signals, since they are meant for synchronization of harmonic oscillator frequencies. However, in practice it is unlikely that this approach can be used: The Kuramoto model describes a continuous process, which means that the oscillator frequencies must be adjusted within a fraction of the signal period $T_m=2\pi/\varphi'_m$ to come close to continuous behavior. For a CW signal of several GHz, this update rate is much too high to be realized with a digital circuit. In addition, continuous adjustment of the oscillator frequency is inconsistent with the way radar and communication systems are operated: A CW radar must take multiple samples while transmitting a CW signal of constant frequency to be able to determine the velocity of targets. Similarly, pilot symbols in communication protocol packets remain at a constant frequency and vary from packet to packet only. For this reason, continuous adjustment of the CW signal frequency is neither desirable nor technically feasible.

Synchronization of pilot tones

Instead, we propose a self-synchronization method based on CW signals whose frequency is kept constant during time duration T. In the following, a CW signal of constant frequency and duration T will be denominated as a pilot tone. If pilot tones of different frequencies follow one after another, the phase of the signal generated by system m can be described after a pilot tones by

(9)\begin{equation} \varphi_{m,a} = \int\limits_0^{aT} \omega_m(t) \,\mbox{d}t = T\sum_{p=0}^a \omega_{m,p} = 2\pi T\sum_{p=0}^a f_{m,p}, \end{equation}

where $f_{m,p}$ is the frequency of the pth pilot tone of system m. Rewritten as an update equation, it becomes

(10)\begin{equation} \varphi_{m,a} = \varphi_{m,a-1} + 2\pi T f_{m,a},\quad \varphi_{m,0} = 0. \end{equation}

In the following, it will be assumed without loss of generality that system 1 is the system under control and system 2 is the interferer. If applying the Kuramoto model described by (4), the frequency of the oscillator of system 1 is adjusted according to

(11)\begin{equation} f_{1,a+1} = f_{1,0} + \epsilon_1\ g(\psi_{21,a}), \end{equation}

where $\varepsilon_1 = 2\pi\epsilon_1$. By substitution of (10), $\psi_{21,a}$ can be determined to

(12)\begin{align} \psi_{21,a} &= \varphi_{2,a}-\varphi_{1,a} \end{align}

(13)\begin{align} &= \varphi_{2,a-1}-\varphi_{1,a-1} + 2\pi T f_{2,a}- 2\pi T f_{1,a} \end{align}

(14)\begin{align} &= \psi_{21,a-1} + 2\pi T f_{\Delta,a}. \end{align}

Here, $f_{\Delta,a}$ is the CFO between the pilot tones of system 2 and system 1. The calculation of $\psi_{21,a}$ corresponds to the sampling of the phase of a time-varying signal. To fulfill the Nyquist sampling criterion in this case, $T f_{\Delta,a} \leq 0.5$ is required. This is achieved for (14) with an arbitrary scaling factor κ:

(15)\begin{equation} \psi_{21,a} = \psi_{21,a-1} + 2\pi\frac{T}{\kappa} f_{\Delta,a}. \end{equation}

Now, the Nyquist criterion $T f_{\Delta,a}/\kappa \leq 0.5$ holds. Since $f_{\Delta,a}\leq 0.5 f_\mathrm{ADC}$ must be equally satisfied, where $f_\mathrm{ADC}$ is the sampling rate of the radars analog-to-digital converter (ADC), it can be derived:

(16)\begin{equation} \kappa \geq T f_\mathrm{ADC}\quad\Rightarrow\quad \kappa = N_\mathrm{s} T f_\mathrm{ADC}, \; N_\mathrm{s} \geq 1. \end{equation}

By substitution into (15), it follows

(17)\begin{equation} \psi_{21,a} = \psi_{21,a-1} + 2\pi\frac{f_{\Delta,a}}{N_\mathrm{s}\, f_\mathrm{ADC}}. \end{equation}

$N_\mathrm{s}$ will be denoted as sampling factor in the following.

The synchronization procedure therefore works out as follows: During time duration T, system 1 records a pilot tone of the interferer and estimates the CFO $\hat{f}_{\Delta,1}=\hat{f}_2 - f_1$, e.g. by Fourier transformation of the recorded signal. For CW radars, this Fourier transformation is already part of the signal processing chain, since Doppler frequencies are determined in the same way. Afterward, the frequency $f_{1,a+1}$ of the next pilot tone is calculated using (11) after determining the phase difference $\psi_{21,a}$ with (17). After setting the new pilot tone frequency, the process starts again from the beginning.

It shall be noted that (11) and (14) are comparable to the equations for PRF synchronization in [Reference Dallmann14], with the biggest difference being the quantities used to determine $\psi_{21,a}$. However, if compared to [Reference Dallmann15], it becomes apparent that one term of $\psi_{21,a}$ is missing, which is required to compensate the occurrence of more than one or less than one pulse of the interfering radar during one pulse repetition period. Since it is assumed that the signal duration T is the same for both systems, the acquisition of the pilot tone takes the same amount of time, rendering the additional term unnecessary here. However, slightly different pilot tone durations are acceptable as long as two different pilot tones transmitted by the interfering radar do not mix up at the system under control, e.g. due to Fourier transform. This can be prevented with various strategies like synchronizing pilot tone durations, adding pauses in between the tones or advanced processing techniques to separate two subsequent tones received by the radar.

Standard vs. generalized Kuramoto model

g(x) in (1) directly affects the synchronization process, since it determines how much the frequency of a pilot tone has to change depending on the phase difference between the tones. It therefore makes sense to identify other functions in addition to $g_\mathrm{S}(x)$ that could positively influence the synchronization process. Therefore, we investigate how synchronization can be accelerated.

Let $\psi_\mathrm{21,sync}=-\psi_\mathrm{12,sync}$ be the phase difference at which two pilot tone frequencies become synchronized, meaning that $\psi'_\mathrm{21,sync}=\psi'_\mathrm{12,sync}=0$. If synchronization shall be achieved fast, two requirements must be fulfilled:

First, to quickly reach $\psi_\mathrm{21,sync}$, the phase difference ψ ₂₁ needs to change fast over time, meaning that $\psi'_{21}$ must be maximized. From (4) and (5) follows:

(18)\begin{equation} \psi'_{21}=\varphi'_2-\varphi'_1 = \omega_{2,0}-\omega_{1,0} - (\varepsilon_1+\varepsilon_2)\ g(\psi_{21}). \end{equation}

Here, a strong phase change $\psi'_{21}$ is only achieved if $|g(\psi_{21})|$ is close to its maximum value for most values of ψ ₂₁.

Second, for values close to $\psi_\mathrm{21,sync}$, the frequency difference $\psi'_{21}$ needs to change fast to quickly reach $\psi'_\mathrm{21,sync}$. From (18) can be derived that

(19)\begin{equation} \psi''_\mathrm{21,sync}=-(\varepsilon_1+\varepsilon_2)\ g'(\psi_\mathrm{21,sync}). \end{equation}

Here, a strong frequency change $\psi''_\mathrm{21,sync}$ is only achieved if $|g'(\psi_\mathrm{21,sync})|$ is maximized. Since this requirement is independent from where $\psi_\mathrm{21,sync}$ is located, the strongest increase in the function even does not have to be at x = 0, as long as the function is still odd as required in “Theoretical background” section.

According to the last two requirements, $|g(\psi_{21})|$ should take the maximum value for all values of ψ ₂₁ except for $\psi_\mathrm{21,sync}$, where $|g'(\psi_\mathrm{21,sync})|$ must be maximum. Therefore, all requirements can be ideally achieved with the sign function

(20)\begin{equation} g_\mathrm{sgn}(x) = \left\{\begin{array}{ll} +1 & x \geq 0 \\ -1 & x \lt 0 \end{array}\right.,\;x\in[-\pi,\pi). \end{equation}

However, substitution of (20) into (18) shows that this function leads to an oscillation between $\omega_{2,0}-\omega_{1,0} \pm (\varepsilon_1+\varepsilon_2)$, causing a degradation in synchronization accuracy. Therefore, the choice of a suitable $g(\psi_{nm})$ points to a compromise between fast and accurate synchronization.

To investigate this, two different functions for the generalized Kuramoto model will be used, which are

(21)\begin{equation} g_P(x) = \sin(x)^P \end{equation}

and

(22)\begin{equation} g_{\Sigma,P}(x) = \sum_{q=-Q}^{Q}\sin(x-mw)^P, \end{equation}

where P is odd and w, Q are given in the Appendix. Examples for both functions are shown in Fig. 1 and Fig. 2, where the function amplitude is given as dimensionless quantity. It becomes apparent that for $g_P(x)$, the peak around $\pm\pi/2$ gets narrower and steeper with increasing P, thus fulfilling the second requirement. For $g_{\Sigma,P}(x)$, a plateau is formed with steep slopes with increasing P, thus fulfilling both requirements. Moreover, the similarity to $g_\mathrm{sgn}$ increases with higher values of P. It should be noted that also other functions that fulfill one or both requirements could have been selected, such as the Fourier series of a square wave function. The two functions in this paper were chosen mainly for convenience, as they can be set with a single parameter and allow a better comparison with each other due to their similarity.

Figure 1. $g_P(x)$ for $x\in [-180^{\circ}, 180^{\circ} ]$ and different values of P.

Figure 2. $g_{\Sigma,P}(x)$ for $x\in [-180^{\circ}, 180^{\circ} ]$ and different values of P.

In the following, both functions will be used for synchronization processes between two CW radars.

Measurements and results

To test the proposed synchronization method, two monostatic FMCW radars manufactured by Fraunhofer Institute for High Frequency Physics and Radar Techniques FHR covering a chirp bandwidth from 68 GHz to 93.5 GHz and featuring an ADC sampling rate of $f_\mathrm{ADC}=1$ MHz were used. For this experiment, they were operated as CW radars around 79 GHz. As shown on both images in Fig. 3, the radars, recognizable by their white lens antennas, were mounted to supports and pointed in the same direction toward a trihedral located at a distance of approx. 1.5 m. Radar 1 operates by recording measurement data during time T, applying a Hamming window onto the recorded data set, performing a Fourier transform and using thresholding for target detection. Moreover, static targets were filtered out by subtracting the mean from the time domain signal before Fourier transformation to suppress reflections caused by antenna mismatch and static surrounding. The frequency of the highest peak found by a maximum search is used as an estimate for the CFO $\hat{f}_{\Delta,a}$. This value is subsequently used to calculate $f_{1,a+1}$ and to adjust the CW frequency of radar 1. Radar 2 operates the same way, meaning that both radars transmit simultaneously during measurements. After updating their CW frequencies, both radars repeat the process from the beginning. Due to hardware limitations, only the absolute frequency value $|f_{\Delta,a}|$ can be measured, making it necessary to hand over the sign to the signal processing methods of both radars as a-priori information after each measurement. This sign is directly calculated from the values $f_{1,a+1}$ used to set the CW frequencies of both radars.

In the following, three different quantities will be discussed, the definition of which we give here. The CFO between both radars after synchronization $f_{\Delta\text{,sync}}$ is defined to the largest difference between frequencies observed after synchronization took place. The synchronization accuracy given in$\ \mathrm{ppb}$ is the aforementioned CFO in relation to the carrier frequency at 79 GHz. The number of pilot tones or iterations until synchronization is achieved $I_\mathrm{sync}$ is measured to be the number of tones until the CFO falls below $f_{\Delta,sync}$ for the first time.

Standard Kuramoto model

Initially, the standard model with $g(x) = g_\mathrm{S}(x)$ will be used. Measurements were carried out using the setup shown on the left side of Fig. 3. $f_{1,0}$ of radar 1 was set to 79 GHz and $f_{2,0}$ of radar 2 was set 0.4 MHz above $f_{1,0}$. The sampling factor was chosen to $N_\mathrm{s}=64$. The pilot tone duration was set to T = 4 ms, resulting in $T f_\mathrm{ADC}=4000$ samples per pilot tone. While keeping $\gamma_1=1.2$, the value of γ ₂ was varied between 0 and 9.6. In Fig. 4, the variation of $f_{1,a}$ (dashed lines) and $f_{2,a}$ (solid lines) is shown for 64 pilot tones. For $\gamma_2=1.2$, both radars adjust their CW frequencies using the same gamma factor. Due to this, they synchronize to a frequency lying approximately in the middle of $f_{1,0}$ and $f_{2,0}$. For $\gamma_2=2.4$ to $\gamma_2=9.6$, radar 2 forces the synchronization process to frequencies which are closer to $f_{1,0}$. The downside of choosing high values for γ ₂ is a high frequency fluctuation, resulting in low synchronization accuracy. However, it comes with the advantage that synchronization is achieved after a lower number of pilot tones. For $\gamma_2=0$, radar 2 does not adjust its CW frequency. Nevertheless, synchronization still takes place, since radar 1 adjusts its frequency to $f_{2,0}$. This demonstrates that self-synchronization can still be successful if one system is non-cooperative.

Figure 3. Measurement setups. Left: Setup realized for the standard Kuramoto model. Right: Setup realized for the generalized Kuramoto model.

Figure 4. Synchronization process for different γ ₂.

Next, $f_{1,0}$, $f_{2,0}$, and T were kept the same and $\gamma_1=\gamma_2=1.2$ was chosen. $N_\mathrm{s}$ was varied between the values 16, 64, and 256. In Fig. 5, the absolute frequency difference $|f_{\Delta,a}|$ is shown on logarithmic scale for 256 pilot tones on the y-axis, whereas the number of iterations is shown as dimensionless quantity on the x-axis. It can be seen that higher synchronization accuracies are achieved if a higher sampling factor is chosen: For $N_\mathrm{s}=16$ the CFO $f_{\Delta,sync}$ is below 19.73 kHz, agreeing with a synchronization accuracy of 256.23 ppb, for $N_\mathrm{s}=64$ it is below 3.77 kHz, agreeing with a synchronization accuracy of 48.96 ppb and for $N_\mathrm{s}=256$ it is below 147.6 Hz, agreeing with a synchronization accuracy of 1.92 ppb. However, better synchronization accuracies come at the price of a slower synchronization process: For $N_\mathrm{s}=16$, synchronization is achieved after $I_\mathrm{sync}=7$ pilot tones, for $N_\mathrm{s}=64$ after $I_\mathrm{sync}=47$ pilot tones and for $N_\mathrm{s}=256$ after $I_\mathrm{sync}=109$ pilot tones.

Figure 5. Synchronization process for different $N_\mathrm{s}$.

Afterward, $N_\mathrm{s}=256$ was chosen while keeping $f_{1,0}$, $f_{2,0}$, γ ₁, and γ ₂ the same as before. Time duration T was varied between $0.5\ \mathrm{ms}$, $1.0\ \mathrm{ms}$, $2.0\ \mathrm{ms}$ and $4.0\ \mathrm{ms}$. The use of longer pilot tones affects the estimation of the CFO $|f_{\Delta,a}|$ because it reduces the frequency point spacing in the frequency domain. Longer pilot tones therefore allow $|f_{\Delta,a}|$ to be estimated more accurately. This can be seen in Fig. 6, where using T = 0.5 ms resulted in a lower synchronization accuracy compared to $T\geq 1.0$ ms. Nevertheless, the synchronization accuracy could not be improved for values above T = 1.0 ms, possibly because the chosen sampling factor $N_\mathrm{s}=256$ is the dominant influence on synchronization accuracy from there on. Thus, sampling factors smaller than $N_\mathrm{s}=256$ would possibly allow for time durations of T = 0.5 ms and smaller. Since the radars are limited to a minimum time duration of T = 0.4 ms, this could not be investigated.

Figure 6. Synchronization process for different T.

Generalized Kuramoto model

After studying the standard model, its generalized version was investigated by application of $g_{P}(x)$ and $g_{\Sigma,P}(x)$. As parameters, $f_{1,0}=79$ GHz, $f_{2,0}=f_{1,0}+0.4$ MHz, $N_\mathrm{s}=64$, $\gamma_1=\gamma_2=1.2$, and $T=4\ \mathrm{ms}$ were chosen. Since the experiments were conducted at a later date, the last setup was reproduced as shown of the right side of Fig. 3. First, $g_{P}(x)$ was tested for different values of P, resulting in the curves shown in Fig. 7, showing the absolute frequency difference $|f_{\Delta,a}|$ on logarithmic scale for 256 iterations on the y-axis and the number of iterations as dimensionless quantity on the x-axis again. In comparison to the standard model, the number of required pilot tones can be reduced by increasing P, achieving a reduction from $I_\mathrm{sync}=47$ to $I_\mathrm{sync}=35$ pilot tones from P = 1 to P = 13. In contrast, $f_{\Delta,sync}=6.49$ kHz increases to $f_{\Delta,sync}=16.43$ kHz, agreeing with synchronization accuracies of 82.15 ppb and 207.97 ppb, respectively. The results show that although faster synchronization can be achieved, the added value of using $g_{P}(x)$ is limited. The reason can be seen from the progression of the curves: For higher values of P, the curves start the synchronization process later, but achieve synchronization faster afterward, leading to a small improvement in total. Since $g_{P}(x)$ only fulfills the second requirement described in “Standard vs. generalized Kuramoto model” section, a significant part of the signals shown in Fig. 1 is close to zero and therefore prevents a fast progression of ψ ₂₁, delaying the synchronization process.

Figure 7. Synchronization process for $g_{P}(x)$.

In comparison, $g_{\Sigma,P}(x)$ outperforms $g_{P}(x)$ as is shown in Fig. 8. Here, a reduction from $I_\mathrm{sync}=51$ to $I_\mathrm{sync}=27$ pilot tones from P = 1 to P = 9 is achieved. $f_{\Delta,sync}=8.26$ kHz increases to $f_{\Delta,sync}=16.49$ kHz with synchronization accuracies of 104.56 ppb and 208.73 ppb, respectively. This shows that a comparable accuracy could be reached with a lower number of pilot tones if compared with $g_{P}(x)$ for P = 13. Although the synchronization speed is similar for P = 9 and P = 13, this behavior seems to be an exception, as a further reduction in the number of required pilot tones can be seen for P = 17. Due to the better performance of $g_{\Sigma,P}(x)$ in comparison to $g_{P}(x)$ since the aforementioned first criterion is not fulfilled for $g_{P}(x)$, the focus is placed on $g_{\Sigma,P}(x)$ in the following.

Figure 8. Synchronization process for $g_{\Sigma,P}(x)$.

Figure 9 allows to study the performance of $g_{\Sigma,P}(x)$. It shows both $I_\mathrm{sync}$ as well as the synchronization accuracy for different P values on the y-axes and the number of iterations as dimensionless quantity on the x-axis. The nonlinear dependency of the number or required pilot tones on P hints at the added value of using the generalized model: For low values of P, a considerable synchronization speed-up can be achieved without a significant deterioration in synchronization accuracy. For higher values of P, this advantage is lost: Comparing Fig. 5 with $N_\mathrm{S}=16$ from last section with the results shown for P = 61 in Fig. 9, it becomes apparent that in both cases synchronization is achieved after $I_\mathrm{sync}=9$ tones, but the change of the sampling factor from $N_\mathrm{S}=64$ to $N_\mathrm{S}=16$ affects synchronization accuracy significantly less. This also becomes apparent from evaluating the product of $I_\mathrm{sync}$ and synchronization accuracy, which is highest for P = 11. It is therefore recommended to apply to the generalized Kuramoto model only small P values around P = 11 to gain an advantage in terms of synchronization speed and accuracy.

Figure 9. Change in synchronization speed and accuracy for $g_{\Sigma,P}(x)$.

A comparison of synchronization performance for $g_{P}(x)$ and $g_{\Sigma,P}(x)$ is shown in Table 1. The table also contains results from the last section for comparison. It should be noted that by comparing the three cases where $N_\mathrm{S}=64$ and P = 1, a significant variation in synchronization accuracy can be determined. This shows that the synchronization is sensitive to environmental influences, leading to performance fluctuations. However, during the experiments, consideration was given to changing $N_\mathrm{S}$ and P over a range of values that would result in significant changes in synchronization accuracy, which is why the discussed results retain their validity.

Table 1. Performance comparison of different synchronization parameters