Over-the-Air Synchronization of Time, Phase and Frequency in a Distributed Receiving System

This paper addresses over-the-air synchronization of two distributed receivers connected to a fusion center via digital links. We distinguish between received signals’ synchronization and receiving channels’ synchronization. The applications are digital receive beamforming (BF) and radio source localization, respectively. We introduce the system and signal models and propose a synchronization procedure. The scheme for signal synchronization includes estimating and compensating for the instantaneous phase and time offsets between the received user signals. This is achieved by using wideband and narrowband pilots sent by a beacon transmitter, along with a wideband preamble sent by the user transmitter. Optionally, the scheme allows for equalization of the user signal, employing a narrowband pilot signal sent by the user transmitter. For channel synchronization, the procedure involves estimating and compensating for the instantaneous phase and time offsets between the receiving channels using pilot signals sent by the beacon transmitter. We assume variable frequency and constant time offsets during the observation interval. In addition, we propose a novel adaptive algorithm for instantaneous phase offset estimation. We demonstrate the performance of the synchronization procedure and the proposed algorithms through Monte Carlo simulations and experiments using USRPs (Universal Software Radio Peripheral). With a 20 dB signal-to-noise ratio and realistic frequency offset amplitude and dynamics, the 0.9-quantiles for time estimation error and instantaneous phase estimation error are one 20th of the sampling interval and 10 degrees, respectively. The total error is reduced almost to the noise-only error.


I. INTRODUCTION
Spatially distributed radio transceivers offer great benefits in distributed beamforming (DBF), such as improved directivity and energy/spectral efficiency, as documented in previous research [1], [2], [3], [4].They also enhance radio source localization, leading to improved accuracy, as demonstrated in earlier studies [5], [6].These distributed transceivers enable the use of single-antenna elements as components of a virtual array, thereby relaxing constraints on antenna The associate editor coordinating the review of this manuscript and approving it for publication was Hans-Peter Bernhard .
apertures and the number of antenna elements, which are characteristic of classical centralized antenna arrays.
Each transceiver uses its own local oscillator (LO) and A/D (analog-to-digital) or D/A (digital-to-analog) converter, resulting in phase, frequency, and time offsets between them.To ensure proper array processing, synchronization in phase, frequency, and time among the array elements is essential.This synchronization challenge is a significant hurdle in practical implementations of systems with distributed transceivers.
One logical approach to tackle this challenge is to distribute a common reference signal to every transceiver.The simplest method for achieving this is through the use of measured (matched) coaxial cables.However, synchronization quality tends to degrade with increasing cable length [7].Another approach, which permits synchronization over greater distances, involves the use of optical fibers [8].This method can be effective with large co-located antenna arrays but may be expensive or cumbersome when deploying many distributed single-antenna nodes, especially in scenarios where nodes are attached to different buildings/posts or different walls/ceilings in indoor spaces.
It is important to highlight that the use of any type of cables is impractical for mobile receivers, although this paper does not address this scenario.

A. THE STATE-OF-THE ART
A decade ago, a cable-free solution for distributing references was proposed [1].This solution is based on RF signals and is flexible and more scalable for synchronizing distributed transceivers over-the-air.Initially, most of the authors who dealt with over-the-air synchronization addressed the problem in the context of distributed transmit beamforming (BF) [1], [2], [3], [4], [9], [10], [11], [12], [13].Their synchronization methods relied on a master-slave protocol, and BF coefficients were calculated based on Channel State Information (CSI) obtained through cooperation with the intended receiver [1].
One solution involved a protocol for time and frequency synchronization of nodes connected to a digital backbone [2].In this paper, the authors also introduced a joint maximum likelihood (ML) time and frequency estimator and achieved phase synchronization during BF coefficient calculations.The same group of authors implemented a system that achieved timing synchronization within an Orthogonal Frequency Division Multiplexing (OFDM) cyclic prefix and produced phase synchronization within a few degrees [9].
Another approach addressed phase synchronization of distributed transmitters connected to an Ethernet backbone [13].Additionally, an approach for frequency synchronization of distributed transmitters and receivers was proposed in [14].An all-wireless distributed BF system was discussed in [4], where the authors employed a Costas loop for frequency synchronization and one-bit feedback [10] to adjust Tx phases, allowing signals to add together constructively at the Rx.
Similar implementations were presented in [11], where the Kalman filter was applied instead of the Costas loop, and in [12], where the receiver played the role of the master in the frequency estimation process instead of one of the transmitters.
A new receive BF method that does not use CSI was proposed in [15].The method is based on the assumption of perfect time synchronization and on the use of a reference signal transmitted along with data signals.In [16], the authors present an implementation of carrier and timing synchronization for distributed arrays.In [17], a new scheme for remote carrier phase synchronization was introduced, with emphasis on reducing the used bandwidth for synchronization.The central unit transmits two tones with similar frequencies, and the remote unit responds with a tone at the frequency, which is the average of the frequencies of the first two tones.
In [18], the authors proposed a method to synchronize the carriers of remote clock sources in the context of distributed cooperative remote sensing applications.To cope with the interference between the transmitted and received synchronization signals at the same transceiver, which is significant in a single carrier frequency scenario, the method uses two frequency carriers per direction of propagation with equal frequency midpoints to achieve phase synchronization between a master and a slave.
In [19], the authors describe some use cases that motivate high-precision over-the-air time synchronization, such as industrial automation and smart grids, and discuss requirements and some possible enabler technologies in wireless systems.
Different sources of time synchronization error for an integration of a wired and a 5G wireless network were modeled in [20], and their impact on the end-to-end error was analyzed through simulations.
A method for over-the-air frequency synchronization between transceiver panels contained within a geographically distributed massive MIMO system was presented in [21].This setup plays a critical role in BF towards a user terminal, enabling coherent reception.
For a distributed system that has already been frequencysynchronized, [22] introduced and experimentally verified a method for estimating phase and sampling time offsets for each channel when the antennas are in dominant LoS conditions and at known positions.
A technique for spontaneous synchronization of the pulse repetition frequency between two independent radars was discussed in [23].The simulations show that their frequencies can converge to each other even if one of the radars is noncooperative.
Interested readers are also referred to recently published review papers on the subject of synchronization [3], [24], and [25].
Our motivation is to propose a procedure that jointly estimates the time and frequency offsets while continuously updating the phase offset.This approach offers advantages over existing methods because it compensates for the coupling effect between the time and frequency offset and successfully copes with abrupt changes in the phase offset.Furthermore, it is suitable for hardware deployment, which remains cost-effective even if the antennas are highly dispersed throughout the served area.

B. THE CONTRIBUTIONS OF THE PAPER
In this paper, we propose a procedure for over-the-air synchronization of distributed receiving systems.We address time, frequency, and phase synchronization of such systems for the purpose of BF toward an intended transmitter (the user) or for radio source localization.
The system consists of distributed network nodes, each with only one antenna (or a small number) and a digital channel (RF chain).For each user Tx, a subset of these receiving nodes is chosen to ensure they are in LoS conditions with the user for most of the time.To cover the entire area served by the network in this LoS manner, multiple beacons (referent transmitters that are part of our system) might be required.
The receivers are connected to a fusion center via digital links, and they are assumed to be incapable of accurately conveying reference signals [2].Although we deal with a two-channel system, the analysis and results can readily be extended to more than two receiving channels by considering pairs of channels.
In the context of localization, we focus on the joint estimation and compensation of the instantaneous values of phase and time offsets between the receiving channels.For BF we have two stages, where in stage 1 we perform joint estimation and compensation of the instantaneous phase and time offsets between the received user signals, and in stage 2, optionally, equalization of the user signal.In the localization case, we use a beacon transmitter that sends both a wideband and a narrowband pilot signal.In the BF case, besides the beacon transmitter sending a wideband and a narrowband pilot signal, the user sends a wideband preamble, and if stage 2 takes place, it sends a narrowband pilot.
In our work, we chose the length of the observation interval to be short enough that the time offsets are approximately constant within it, yet long enough to ensure a stable estimate.A more detailed discussion is provided in Section VI.However, it is important to note that the frequency offsets vary over time.This approximation works well for the Universal Software Radio Peripheral (USRP) platforms [26] and the observation intervals used in our experiments.Therefore, we do not explicitly estimate the frequency offsets; instead, we propose an adaptive algorithm for instantaneous phase offset estimation.
In our procedure, the phase offset is continuously estimated, in contrast to the approaches proposed in the references above, except for [4] and [15].However, these references do not consider time shift estimation.We also observe that most of the cited papers deal with transmit BF, where it is more demanding to continuously track the phase offset (requiring an additional antenna and frequency band for reception of a reference signal), and it is impossible to use the estimated phase offset in real-time, requiring some form of prediction.Considering the round-trip delay, including the algorithm execution time, if the variability of the frequency offset allows for accurate prediction this far into the future, our procedure can also be applied to transmit distributed BF.
We consider a spatially coherent signal scenario [5] with a dominant LOS component, such as one within the mmWave range [6], although our procedure is not limited to such scenarios.In this setting, we calculate the coefficients of the user signal at the receiving channels based on the estimation of the phase offset and Time Difference Of Arrival (TDOA).These tasks are significantly less demanding than estimating classical CSI.For time shift estimation, we propose a noncoherent ML-type algorithm [6] due to the presence of relative phase offsets in the signals [27].To estimate the residual constant phase offset, we propose a numerically efficient algorithm based on signal correlation.We evaluate the proposed algorithms through simulations and experiments using USRP platforms [26].
The rest of the paper is organized as follows: The system model and problem formulation are presented in Section II.In Section III, we introduce our signal model.We describe our synchronization procedure in Section IV and the used algorithms in Section V.In Section VI, we provide a description of our simulations and experiments and present various numerical results for the performance of the proposed algorithms and the procedure.We conclude the paper in Section VII with some final remarks.In the Appendix, we provide a detailed mathematical presentation of the proposed synchronization procedure.

II. SYSTEM MODEL AND PROBLEM FORMULATION
Our receiving system consists of M distributed receivers, denoted as Rx 1 , Rx 2 , . . ., Rx M , connected via digital links to a fusion center (see Fig. 1).In our model, the digital links are wired (e.g., Ethernet), but wireless links can be used instead to allow receiver movement.Each receiver (Rx) operates independently with its own LO and A/D converter.Therefore, there is neither time, frequency, nor phase synchronization among the receiving channels.However, we assume that a rough time synchronization exists in the sense that the observation intervals of the channels overlap.Signals from two different transmitters (Tx) at separate locations also overlap within the same observation interval.128200 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.The user Tx emits an RF signal carrying useful data, but these data are only present in the signal segment processed during the data transfer phase (see Section IV).In this section and in Section III, we only deal with the first part of the signal, which is processed in the initial phase (see Section IV).That being said, the user signal includes a wideband preamble and, optionally, a narrowband pilot.We denote the corresponding baseband frequency of the signal as f U , as shown in Fig. 2. The synchronization Tx (the beacon) emits wideband and narrowband pilot signals at corresponding baseband frequencies f W and f N , respectively.If the user narrowband signal is present (explained in Section IV), it can occupy a separate portion of the spectrum beside the wideband preamble/useful data, similar to the beacon signals.
Alternatively, it can be one subcarrier if the user data is in OFDM format (this case is shown in Fig. 2 and considered in Section III).
The location of the user Tx is unknown to the receiving system.The requirements regarding knowledge of the beacon's location are discussed in Section IV.There is a LOS connection between every Tx-Rx pair.The signals received at every receiver are IQ (in-phase quadraturephase) demodulated and A/D converted.The resulting digital samples are then delivered to the fusion center for joint processing.We also assume that all A/D and D/A converters in the system have the same sampling frequency, denoted as fs .
For simplicity, we consider a two-channel receiving system (M = 2), as shown in Fig. 3.The results and conclusions can be readily extended to settings with more than two receiving channels by considering additional pairs of receiving channels.
The variables in Fig. 3 have the following meanings: the phase of the local oscillator of any transceiver is given by ω c t + ϕ(t), relative to that transceiver's local time.The term ϕ(t) also includes the (variable) frequency misalignment of the LO relative to a nominal carrier frequency, ω c .The frequency and phase misalignments of the beacon Tx, user Tx, Rx1, and Rx2 are denoted as ϕ B (t), ϕ U (t), ϕ 1 (t), and ϕ 2 (t), respectively.The lack of time synchronization of the beacon, user Tx, Rx1, and Rx2 is modeled by τ B , τ U , τ T1 , and τ T2 , respectively.These values represent the delays of their local times relative to a nominal time axis.The signal propagation delays are denoted as τ PB1 , τ PB2 , τ PU1 , and τ PU2 .The length of the observation interval is chosen to allow the time shifts to be modeled as constant within it and to account for the phase shifts they cause, which can be compensated by the phase estimation algorithms.The beacon and user Tx waveforms are denoted as s B (t) and s U (t), respectively.In general, where s W (t) and s N (t) are the wideband pilot and narrowband pilot waveforms, respectively, and where s D (t) and s UN (t) are the user data signal and user narrowband pilot waveforms, respectively.The goal is to achieve time, frequency, and phase synchronization of the distributed receiving channels by utilizing the beacon signals, enabling user localization, or achieving signal synchronization of the user data received in the channels by utilizing both the beacon and user signals, which enables digital receive BF.

III. SIGNAL MODEL
Let s B = [s B0 , s B1 , . . ., s BN −1 ] T be a complex sequence at the input of the D/A converter of the beacon.The continuous waveform at the output of the D/A converter is denoted as s B (t), with respect to the local time of the beacon.The variable t is normalized with 1/ fs .This signal is IQ upconverted, and the resulting RF signal is then transmitted via the antenna.In the following equations, we express the baseband equivalent of the emitted signal with respect to the nominal time as where ω c is the angular carrier frequency normalized by fs .Unless otherwise stated, all time and frequency variables throughout the paper are normalized in the same way as t and ω c , respectively.Non-normalized, natural variables are denoted with a tilde (˜) above the variable.
Similarly, if we denote s U (t) as the continuous waveform at the output of the D/A converter of the user Tx according to the local time of the user, the baseband equivalent of the signal it emits according to the nominal time is given by These two signals propagate to the receiving antennas, and the baseband equivalents of the signals received by channels 1 and 2, each represented in their respective local times, are defined by where a m is a constant real scalar factor representing the signal attenuation due to propagation, m ∈ {1, 2}, and η m (t) is independent circularly symmetric complex Gaussian noise in the frequency range [−1/2, 1/2).We assume that a m = 1 (∀m), so that the noise variance for a given channel models the propagation attenuation.We recall that a time shift always implies that there is a corresponding phase shift as well, since we deal with RF signals [27].
Starting from ( 5)-( 6), we derive a discrete-time matrix baseband model of the signals at the receivers.We have where s m = [s m (0), s m (1), . . ., s m (N − 1)] T , N is the sequence length, and the sequences s EkvB , s EkvU , s EkvW , s EkvN , s EkvD , s EkvUN , and η m are defined similarly.
For any complex sequence s = [s 0 , s 1 , . . ., s N −1 ] T , the DFT (Discrete Fourier Transform) of s is: and the ideal periodic continuous-time signal (waveform) whose period corresponds to the sequence s, s ideal (t), is given by The ideal signal is bandlimited to [−1/2, 1/2) and is equivalent to the signal at the output of an ideal D/A converter.This means that the samples at its input are identical to the corresponding values of the analog signal at its output at the sampling instants when s is at its input.Using this method, we can uniquely represent any sequence with a function defined for any real argument.This is especially useful since we are dealing with non-integer time shifts.In practice, the signal segment used for time shift estimation contains only a few periods of the sequence.
The symbol m represents a matrix operator that models phase shifts in the channels.It contains phase terms defined by the vector −ϕ m on its main diagonal.This vector comprises phase shifts for each sample.The symbol D τ represents a time-delay-by-τ operator, which models both the complex envelope and carrier phase shifts, essentially representing the RF signal delay, while D BB τ is a timedelay-by-τ operator that models only the complex envelope shift, representing the baseband signal delay.Finally, F is a modified DFT matrix.It is structured such that F −1 = F H , and its rows are arranged based on natural RF frequencies.In a more formal representation, where exp(•) is the element-by-element exponential function, and Diag{•} is a diagonal matrix with the given elements on its main diagonal.As a short digression, Fig. 4 illustrates time shifts in a signal received by two receivers Rx 1 and Rx 2 , with propagation delays τ P1 and τ P2 , and delays of their time axes τ T1 and τ T2 , respectively.The Tx time axis is assumed nominal.Note that, as opposed to a Tx, a delay of the observation interval in an Rx for some τ , means a delay of −τ for the received signal in that channel.
The signal model ( 7)-( 8) is suitable for scenarios where the original transmitted sequence is unknown to the receiving system or not used in the estimation.In such cases, the time and phase offsets of the Tx are irrelevant, and we focus solely on the signal at the output of the Tx, denoting it with the 'Ekv' index.However, if the receiving system has knowledge of the original sequence and utilizes it in the algorithms, we must explicitly model the Tx time and phase offsets.In this scenario, the discrete-time matrix baseband model of the signals at the receivers is given by where B = Diag exp (jϕ B ) , and U is defined similarly.
128202 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

IV. THE PROPOSED SYNCHRONIZATION PROCEDURE
The procedure we propose takes two signals as input, compares them, and estimates the relative time and phase deviations between them, producing the compensated version of the second signal as the output.Recall that the time shift is approximately constant within the observation interval, so its estimate is a scalar, whereas the phase shift is variable, so its estimate is a function of time (a vector of time samples of the same length as the signal segment).The input signals contain the same components, which may have experienced different distortions and noise.The idea is to use various algorithms, some for estimating time shifts and others for frequency/phase shifts, in four steps to progressively estimate and eliminate all these distortions (except for the noise).
To successfully estimate a time shift, we need to compare wideband signals; otherwise, we would not be able to distinguish between time and phase shifts.For this purpose, we use a wideband pilot, which is a component of the signal.Since the observation interval is short enough that the time shift is a scalar, we can obtain an estimate of the time shift from just the beginning of the signal (in the time domain) and use the estimate to compensate for it during the entire signal segment.This is why we refer to the wideband pilot in the user signal as a ''preamble'' throughout the rest of the paper.Moreover, the procedure is designed to execute in real-time during exploitation, so we refer to the processing of the preamble as ''the initial phase.''Once it is complete, there is no need for the wideband pilot until the next signal segment is acquired.In the meantime, only compensation takes place, which we refer to as ''the data transfer phase.''Therefore, for one observation interval, we distinguish between the initial phase and the data transfer phase.
To estimate a phase shift, on the other hand, we can use a signal component of any bandwidth.In fact, using a narrowband pilot is preferable because the narrower the band, the less noise energy we have to deal with.Since we are considering a phase shift that varies with time (it implicitly contains the frequency shift), we need the corresponding pilot to last for the entire duration of the signal segment.
Our system contains a network of antenna nodes spread across the served area.Each of them has its own distortions due to its local oscillator/clock.We propose applying the procedure to pairs of receiver nodes (the Rx channels) to synchronize them, i.e., to remove their relative distortions.We do this in stage 1 of the procedure, as shown in Fig. 5.In step 1, we roughly estimate (and compensate for) the time shift.In step 2, the algorithm works sample by sample to estimate the time-variable component of the phase shift.In step 3, we make an accurate estimate of the remaining time shift.Due to the nature of the problem, the step 3 algorithm introduces an additional phase distortion in the signal, but fortunately, it is time-constant.In step 4, we accurately estimate the remaining (time-constant) phase shift.All the steps except step 2 use block-estimation algorithms (as opposed to sample-by-sample), and we execute them in the initial phase only, but use their results throughout the signal segment.
At a first glance, it may appear that not all of the steps are necessary.However, because the effects of a time shift and a phase shift are coupled, each step serves a specific purpose.Namely, if we did not execute step 1 before step 2, the step 2 algorithm would not be able to compare corresponding samples -those with phase distortions that correspond to approximately the same sample of the useful signal component.Consequently, it would not produce sensible results.The step 3 algorithm can only provide an accurate time estimate if we remove the frequency distortion first.For the same reason, step 1 cannot make an accurate time estimate.Consequently, step 2 compares samples of signal 1 to the values of signal 2 at time instants that are slightly in the past or future.Given that the narrowband pilot has its own frequency, this introduces a bias in the form of a constant phase shift, which is later removed in step 4.
Note that we require these pilot signals to be present with a high enough SNR and favorable channel conditions throughout the entire duration when synchronization is necessary for the system (network).It might be overly optimistic to expect this level of consistency from user-terminal signals alone.This is why we incorporate a beacon Tx as a critical component of our system.We have control over its placement and the structure of its signals, providing several advantages.One of the key benefits of this approach is that it allows us to ensure network-wide coverage.Multiple beacons may be necessary to cover the entire service area of the network, guaranteeing that we have sufficiently large subsets of Rx nodes in LoS conditions with at least one beacon Tx.
If we apply our procedure to nodes Rx 1 and Rx 2 (without loss of generality), in the form of stage 1 depicted in Fig. 5, we align the received beacon signals in both time and phase at point A 4 .However, because the propagation delays from the beacon to Rx 1 and Rx 2 differ, we must compensate for the TDOA to synchronize these two channels.Consequently, we achieve a synchronized pair of receivers at A 5 .
Recall that we receive signals from both the beacon and all user terminals (served by a subset of Rx nodes) simultaneously.This implies that all distortions encountered in stage 1 are also compensated for in the user signals.Consequently, the user signals at A 5 are ready for use as inputs to a localization algorithm.However, if we need to perform BF for an individual user, we must first align the user's signals.For this purpose, we can employ a state-ofthe-art BF algorithm or apply steps 3 and 4 of our procedure to the user's signals (instead of the beacon signals).We point out that we require knowledge of the beacon's TDOA, and most likely its exact position, to obtain the A 5 signals from A 4 .Nonetheless, we can relax this requirement and use the A 4 signals (see Fig. 5) as inputs to the BF algorithm.We can even skip steps 3 and 4 in stage 1 and directly align the user signals, using the A 2 signals as inputs to BF, to save processing time.However, estimates based on beacon signals are likely to be of higher quality, and the number of search grid points in BF could potentially be reduced.Therefore, it involves a trade-off, as the total numerical complexity might be worse if we use the A 2 signals.In Section VI and the Appendix of this paper, we use the A 2 signals as inputs to BF.Thus, the implemented procedure path for BF is A 0 -A 1 -A 2 -A * 3 -A * 4 , marked by a solid arrow in Fig. 5.The other possibilities are indicated by dashed arrows.
In the BF application, once the user signals are aligned, we sum them to obtain a beamforming gain, resulting in increased SNR in B 0 , as shown in Fig. 5.At this point, we still cannot decode the user data because the user terminal has its own distortions relative to the Rx system (the user Tx is not synchronized with the Rx).To address this, we run the procedure again, denoted as stage 2 in Fig. 5.This time, we apply it to the original transmitted training signal (preamble + narrowband pilot) and the combined signal from B 0 as the inputs.Finally, at B 4 , the Rx signal is aligned with the Tx signal, and the decoding process can commence.
Next, we provide a comprehensive description of the procedure option marked by a solid arrow in Fig. 5.The system and signal models outlined above pertain to the initial signal phase, and a significant portion of this section is dedicated to it.Consequently, we will begin by providing a detailed explanation of the estimation and compensation procedure for the initial phase of the signal synchronization scenario, i.e., the BF scenario.Subsequently, we will highlight the differences that arise during the data transfer phase compared to the initial phase.Finally, we will identify the variations in the procedure for channel synchronization (localization scenario) relative to signal synchronization (BF scenario).

A. INITIAL PHASE FOR BF
The initial phase involves processing the first signal segment within the observation interval.In the context of BF, this segment includes the wideband preamble and, optionally, the user's narrowband pilot, transmitted by the user Tx, as well as the narrowband and wideband pilots transmitted by the beacon, as shown in Fig. 6.In this figure, we see that stage 2 is also performed, and the user's narrowband pilot occupies a portion of the data signal spectrum.Both the preamble and wideband pilot consist of periodic signals, where one period corresponds to the sequences s U and s B , respectively.During the time offset estimation, a correlation window is constructed, consisting of N samples of one of the signals, which is then shifted along N samples and compared with the corresponding samples of the other signal for a given shift.It is crucial that the sent signals span the entire signal segment used for correlation; otherwise, there will be a portion of the signal with noise only.Hence, a minimum of three signal periods is sent, and the correlation window includes samples that roughly correspond to the second period.In contrast, for localization, the user's Tx does not send a preamble, while the beacon signal maintains the same structure.
We distinguish between two stages in the synchronization procedure.These two stages can be performed in two separate subsystems.
In stage 1, named the alignment stage, the goal is to align the Rx 2 to the Rx 1 signals, which involves synchronizing their time, frequency, and phase, utilizing both the beacon and user signals.The knowledge of the original sequences is not required or used.Therefore, the wideband preamble does not have to be known, allowing a part of the data signal to serve as the preamble.The synchronization procedure consists of four steps.In each step, we first estimate the offsets in one of the signals received in channel 2 relative to the corresponding signal received in channel 1.Then, we compensate for the signals in channel 2 using the estimated offset.At the end of this stage, the system does not enable correct decoding of the data yet, as it does not provide any type of synchronization with the user Tx.If data decoding needs to be enabled, we proceed to stage 2. The corresponding signals from the channels are added together, and the resulting combined signal serves as the input to stage 2, thereby providing a higher effective SNR.
In stage 2, named the decoding stage, the goal is to achieve equalization of the chain [the user Tx] -[the propagation channel] -[receiving channel 1] (at this point, the user signal of channel 2 is aligned, up to the estimation error, with the signal of channel 1).We emphasize that only the user signals are exploited.The knowledge of the original sequence transmitted by the user Tx is required and used.For this reason, a known preamble occupies the first part of the data signal.Therefore, since stage 1 and stage 2 use the same portion of the user Tx signal, if stage 2 takes place, then the preamble would also be known in stage 1, allowing us to use algorithms that exploit this knowledge for even better performance.The synchronization procedure consists again of four steps.In each step, the combined signals from stage 1 are compared with the corresponding original sequence.Then, the combined signals are compensated for the estimated offset.At the end of this stage, the system does enable correct decoding of the data.
A detailed mathematical presentation of the steps in stages 1 and 2 is given in the Appendix.

B. DATA TRANSFER PHASE FOR BF
The data transfer phase involves processing the next signal segment.In the BF setting, this segment of the signal contains the user data signal, optionally the user narrowband pilot, and the beacon narrowband pilot.In the localization setting, the user signal transmits only the data signal, while the beacon signal remains the same as in the BF setting.The length of the observation interval is limited by the assumption that time shifts remain constant within it.The following observation intervals maintain the same structure.
The procedure is the same as in stages 1 and 2 of the initial phase, with the only exception that in steps 1, 3, and 4, no estimation is performed, but compensation is.The estimates from the corresponding steps of the initial phase are used for compensating the signals in these steps.This use of these estimates is possible because the time shifts remain constant in the observation interval.Furthermore, the wideband pilot and user preamble are not needed in this phase.However, phase offset estimation in step 2 still takes place due to a time-variable frequency offset, which is why narrowband pilots from both the beacon and the user Tx are necessary.It is important to note that time compensation is still performed in two separate steps (step 1 and step 3).Merging steps 1 and 3 would invalidate the phase estimate calculated in step 4 of stage 1 (or stage 2) of the initial phase, as it is coupled with the constant phase error at the output of step 2 (ϕ II err or ϕ UII err , see the Appendix).This, in turn, is coupled with the time estimation error at the output of step 1 (τ I err or τ UI err , see the Appendix).

C. LOCALIZATION
Unlike in the case of BF, for localization the goal is to achieve time, frequency, and phase synchronization of the distributed receiving channels by using the pilot signals from the beacon (channel synchronization instead of signal synchronization).The user Tx signal is not utilized for either time or phase synchronization; only the beacon signals are employed.Therefore, the user narrowband pilot is not needed at all.The signals are displayed in Fig. 7.
In the initial phase, the processing in steps 1 and 2 is identical to that used in BF.In steps 3 and 4, the processing is similar to that in BF, but there are two significant differences.The first difference is that we use the wideband pilot for estimation instead of the user Tx preamble.Compensation for the estimated shifts is applied to the user Tx signal, just as in BF.In step 3, compensation is also applied to the wideband pilot since it is used for estimation in step 4. The second difference is that we additionally compensate the user Tx signal from channel 2 for −TDOA B , effectively delaying it by TDOA B .In BF, time compensation of the signals includes accounting for the time shifts of the receiving channels and the TDOA from the user Tx to the receiving channels.However, for localization, our goal is to compensate only for the time shifts in the receiving channels.By delaying the channel 2 signal by TDOA B , we remove the compensation for TDOA from the beacon to the receiving channels.We deliberately avoid compensating for TDOA B as it contains information about the source location.This is why we rely on estimates from channel synchronization rather than signal synchronization.A detailed mathematical presentation of steps 3 and 4 for the localization case is given in the Appendix.
Similar to BF, the procedure in the data transfer phase (named as such due to the analogy with BF) is the same as in the initial phase, except that in steps 1, 3, and 4, no estimation is performed -only compensation.The estimates from the corresponding steps in the initial phase are also used for compensating the signals in these steps.
Hence, at this point, we can affirm that the distributed receiving system has been transformed into a classical (collocated) receiving system concerning synchronization, as discussed in [6].Furthermore, in the case of coherent localization, it requires spatial coherence within the scenario and could also benefit from enhanced time offset estimation accuracy.Coherent time estimation, under certain conditions, provides improvements over noncoherent estimation by a factor proportional to f c [27], [28].As previously mentioned, our procedure uses noncoherent algorithms for time shifts due to the presence of relative phase offsets in the signals.Some potential avenues for improving accuracy include increasing the SNR and/or signal bandwidth (though there are limitations to both), or expanding the statistical signal sample.

V. ALGORITHMS A. STEP 1 ALGORITHM
In step 1, we use the algorithm for constant time shift (CTSA), which is of ML-type and based on correlation of signals [6].
128206 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
Let us assume that we want to apply the algorithm to two signal vectors, x 1 , x 2 ∈ C N ×1 , containing a useful signal with two different time shifts.The algorithm estimates the time shift of x 2 relative to x 1 .For the sake of numerical optimization, the estimation of the time shift is divided into two substeps.In substep 1, the integer part of the time shift is estimated by Obviously, the integer time shift is implemented as simple shifting of the samples of one of the signals.
Then, in substep 2, the fractional part of the time shift is estimated with a resolution depending on the expected accuracy and according to where Finally, the estimated time shift of step 1 is given by Note that the algorithm does not utilize information from the carrier phase; in other words, the algorithm is noncoherent.Also, the algorithm can be further numerically optimized by calculating Fx 1 and Fx τ I int 1 once during preprocessing and then computing the criterion functions over the search grid in the frequency domain.This approach allows us to avoid the calculation of the Inverse Discrete Fourier Transform (IDFT) and Discrete Fourier Transform (DFT) on the search grid.

B. STEP 2 ALGORITHM
In step 2, we employ an adaptive algorithm for instantaneous phase offset estimation, referred to as APSA.The vector of filter coefficients is denoted as w(n) ∈ C 2×1 , with w(n) = [w 1 (n), w 2 (n)] T and the vector of signals at the two inputs of the filter at time instant n as x(n) = [x 1 (n), x 2 (n)] T .The algorithm's purpose is to estimate the instantaneous phase shift between the signals at its input.We define the optimization criterion for the filter coefficients as the minimization of the mean power of the signal at the output of the filter.This criterion is chosen under the constraint that the absolute values of the coefficients are not allowed to converge to zero, e.g.|w 1 (n)| = 1.This constraint ensures that the output signals are in counter phase.The signal at the filter's output is then calculated as follows: and its mean power is where E is the expectation operator, (•) * denotes complex conjugation, and R = E x(n)x H (n) is the covariance matrix of the input signal.The filter coefficients are adapted in the way that the successive corrections are in the direction of the negative of the mean power gradient (the method of steepest descent [29]), or where µ ϕ is the step size parameter, and ∇ξ (w(n)) is the gradient of ξ (w(n)) [29].If the matrix R is approximated by its instantaneous estimate R(n) = x(n)x H (n), we obtain the LMS (Least Mean Square) algorithm [29] whose coefficients are updated as follows: where The same result is obtained starting from the standard formula for an LMS algorithm, i.e., where d(n) is the desired signal at the filter output.When we set d = 0, which means (with constrained values of the coefficients) that signals in the output branches of the filter (see Fig. 8) are in counter phase, we get: The compensation of x 2 for the instantaneous phase offset is performed according to the following formula: where the negative sign is added because, while the algorithm is designed to minimize the mean power at the filter output, our objective is actually to maximize it.This maximization indicates that the output signals, when added together, are in phase, meaning there is no phase offset between them.
In order for the algorithm to converge, 0 < µ ϕ < 2/ x H (n)x(n) has to hold, where x H (n)x(n) is the power of the signal at the filter input at n. Since the input signal power is variable, it is convenient to write (26) as where μϕ = µ ϕ x H (n)x(n), and where the condition for convergence of the algorithm is 0 < μϕ < 2.

C. STEP 3 ALGORITHM
In step 3, we employ CTSA with improved resolution in substep 2 compared to step 1.

D. STEP 4 ALGORITHM
In step 4, we use an algorithm for mean phase offset estimation, known as CPSA.If x 1 , x 2 ∈ C N ×1 are two signal vectors at the input to CPSA, the algorithm estimates the mean phase offset between them as follows: The algorithm is suitable for scenarios where the phase offset between the signals is constant in time.
Note that our procedure is modular, allowing the replacement of an algorithm in a step with another algorithm of the same type.For instance, in steps 1 and 3 of stage 1, we could use an algorithm that relies on knowledge of the original sequence.This approach offers higher accuracy at lower SNRs but comes at the cost of increased numerical complexity [27].

VI. SIMULATIONS AND EXPERIMENTS A. GENERAL DESCRIPTION
In the remainder of the paper, for the sake of simplicity, we present results on signal synchronization (referred to as the ''BF case'' here).We used a 4QAM signal as the user data signal.This signal was interpolated by a factor of four, and its spectrum was centered around the frequency f U = 0.175, as shown in Fig. 9.
The user's narrowband pilot signal was represented as a cisoid with a frequency of f N = 400/1024, which corresponds to an integer multiple of the Discrete Fourier Transform (DFT) resolution.The wideband pilot, on the other hand, was a Gaussian sequence.This sequence was also interpolated by a factor of four and had its spectrum centered around f W = −0.25.The narrowband pilot had a frequency of f N = −50/1024.
Each of these signals had a length of N = 1024 samples.To create the preamble, we extended each of these signals with one copy of itself on each side.As a result, the preamble length was increased to N prmb = 3072.
During the time offset estimation in step 3, we shifted the correlation window by N samples of channel 1's signal  and compared it with the corresponding samples of channel 2's signal.This process aimed to calculate the value of the criterion function for that specific shift.
128208 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
For the correlation window, we chose the segment of the signal that approximately corresponded to the second period of the signal, specifically the samples 1025 to 2048.This choice helped avoid edge effects if the time shifts did not exceed N .
In step 4, which also involves block-type processing for estimation, we used the samples from 1025 to 2048.
In step 2, where the algorithm operates at the sample level, we used the entire received signal.
In each step, compensation was applied to the entire signal, allowing any part of it to be shown in the results.When compensation was performed on the segment of the signal used for estimation in steps 1, 3, and 4, as described in Section IV, we referred to it as the initial phase.When compensation was applied to the other part of the signal, it was termed as the data transfer phase.
In scenarios where the frequency offset values were high and/or changed rapidly, it was advantageous in step 1 to use a shorter segment for calculating the maximum of the criterion function.This adjustment was due to the fact that frequency offset can distort the waveform, and longer correlation windows tend to accentuate this effect.
Shorter windows experience the distortion as a relatively constant phase shift, which the algorithm can effectively handle.To determine the optimal segment size, we conducted simulations, and a segment size of 50 samples yielded the best overall performance.
Similar to step 3, the correlation window was shifted over one signal period.However, instead of spanning 1024 samples, the correlation window was 50 samples in length.We focused on samples 1525 to 1574, which represented a signal segment approximately in the middle of the second period.
In cases where the correlation window does not include an entire signal period, and the signal energy is not evenly distributed over that period, it becomes essential to normalize the correlation in the criterion function.
The (normalized) carrier frequency was f c = 3960.We configured the step size parameter for the APSA as μϕ = 1.
Due to space limitations, most of the results for stage 2 are omitted.They closely resemble those for stage 1 but exhibit a 6 dB improvement.This improvement is attributed to two key factors: 1) a 3 dB enhancement achieved by adding the signals from the channels before stage 2 (the signal levels in the channels were equal), and 2) another 3 dB improvement as the combined signal was compared with the Tx signal that did not contain noise, unlike in stage 1, where both signals were noisy.
At the conclusion of stage 1, the user signals in channels 1 and 2 exhibit disparities due to independent noise, residual time, and variable phase offsets.Despite these factors, the normalized autocorrelation function of the user signals for these shifts remains close to 1. Consequently, due to the precision of stage 1, it can be inferred that the user signals are nearly identical, particularly for reasonable SNR values.This characteristic enabled us to determine estimation errors by considering the means of the time and phase shifts of the signals in channels 1 and 2 as the initial shifts of the combined signals (relative to the original sequence) in stage 2. In cases where the signals possessed different power levels, we applied appropriate weighting.
To facilitate the comparison between simulation and experimental results, whenever feasible, we present the corresponding outcomes in the same figure.
In the following two subsections, we provide the remaining details regarding the simulations and the experiments.

B. SIMULATIONS
We conducted the simulations using Matlab.We modeled the frequency offset of the transceiver, denoted as f (n), as follows: where α is the autocorrelation of f (n), where the normalized autocorrelation (the autocorrelation coefficient) was ) and ζ (n 2 ) were independent; σ is the standard deviation of f (0) in [1/sample]; β is a scalar factor chosen to ensure E f 2 (n) = σ 2 .The frequency offsets of different transceivers were independent.The parameter ρ = R f (100) specified the normalized autocorrelation of the frequency offset for a shift of 100 samples, see Fig. 10.
In the simulations, we used the same values for σ and ρ for all four transceivers, although the parameters for different transceivers can be different in general.The phase offsets were modeled by where ϕ(0) had a uniform distribution on [−π, π), i.e., ϕ(0) ∼ U [−π, π), and was independent across the transceivers.Fig. 13d shows the phase offset between the channels for different values of ρ.The time offsets were modeled according to where τ B1 , τ , and τ U1 were i.i.d. with distribution U (−100, 100).The difference between the user Tx TDOA and the beacon TDOA was τ rel ∼ U [−0.05, 0.05).

C. EXPERIMENTS
The experiments were conducted in an urban environment, in front of the Innovation center of the School of Electrical Engineering in Belgrade, Fig. 11.We employed two USRP platforms NI-2920 as transmitters, and two USRP platforms NI-2932 as receivers.We connected the USRPs via Ethernet cables to a PC (the blue cables in Fig. 11).Further, we connected the USRPs via matched cables to omnidirectional antennas (the black cables in Fig. 11).The distances between the transceivers were of the order of several meters, with LOS between every Tx-Rx antenna couple.
We prepared the signals of the transmitters in the same way as in the simulations (see Fig. 9 b), except that the signal of length 3072 (used in the simulations) was repeated 10 times.We refer to these 3072-sample-long parts of the signal as signal segments 1-10.The natural carrier frequency was fc = 990 MHz, and the receiver bandwidth was B = 0.25 MHz.Thus, as in the simulations, we had f c = fc /B = 3960.

D. RESULTS OF SIMULATIONS AND EXPERIMENTS
In this subsection, we present the signals in the channels at different stages of the synchronization process for both a simulation and an experimental run.The parameters used for the simulation were as follows: SNR = 20 dB, σ = 1.5 × 10 −3 , ρ = 0.1.
First, we discuss some results related to step 2 of the APSA algorithm.Fig. 12 displays the coefficients of the APSA applied to the signals in channel 2. Since the phase of the signals of channel 2 relative to the signals of channel 1 was of importance, the coefficient for channel 1 was kept at 1.The algorithm converges in just a few samples.These coefficients are time-varying due to the variable phase and frequency offset in the signals.
If we look at the period of the sinusoids in Fig. 12a, which is around 1000 samples, we can estimate that, in the simulation, the frequency offset between the channels was approximately 1 × 10 −3 1/sample.As time progresses, the coefficients deviate from a cisoid due to the time-varying frequency offset.
When there is a time shift in addition to the frequency and phase offset between the signals, the algorithm interprets it as an additional constant phase shift since it deals with complex sinusoids (narrowband pilots).This results in a constant phase error in the estimate of the instantaneous phase offset, as shown in Fig. 13a.
The actual frequency offset between the channels in the simulation is shown in Fig. 13b.Based on the observations in Figs.13c and 13d, we can conclude that for the experiment, ρ is greater than 0.99.Also, the average frequency offset between the channels is calculated as 106 rad/(2π rad • 20000 sample) = 0.843 × 10 −3 1/sample.This value corresponds to the slope of the curve in Fig. 13c.This finding is further supported by observing the period of the coefficients in Fig. 12b.The value of the frequency offset approximately corresponds to the parameter σ = 0.5 × 10 −3 used in the simulations.Fig. 14 displays the time domain of the real parts of the narrowband pilots in the channels at both the input and the output of the APSA.We observe that the user narrowband pilots were not used for the estimation in stage 1.Now we show some results that illustrate the performance of the overall estimation and compensation process.Fig. 15 presents the real parts of the user data signals in the channels, before and after compensation of the offsets, respectively.In the simulation, the obtained normalized BF gain (NBFG) at the end of stage 1 was 0.94.We note that the NBFG is a scalar measure for quality of compensation of offsets between two signals.It is defined by where s 1,in , s 2,in are signals before and s 1,out , s 2,out are signals after offset compensation.In our case, the respective signals are the user data signals in the channels at the beginning and the end of stage 1, respectively.The ideal compensation of offsets corresponds to NBFG = 1, and the obtained NBFG in the experiment was 0.995.In order to confirm that the time offsets were constant in time, we used the signal segment 1 for estimation in steps 1, 3, and 4, and these estimates were employed for offset compensation and NBFG calculation for samples 28673-29696 (the third period of the signal segment 10).Fig. 16 displays the constellation of the user data signals in the channels at the end of stage 1, and Fig. 17 presents the constellation of the combined user data signal for steps 0 to 4 of stage 2. It is presented alongside the same signal at the user Tx (the original sequence).
The analysis of these figures helps in understanding the synchronization and compensation process: • In stage 1, the goal was compensate for offsets in channel 2 relative to channel 1.As a result, no synchronization with the user Tx was achieved, making data decoding impossible.If channel 1 had no offsets, data decoding would have been possible.
• Stage 2 comprised several steps to address different types of offsets: -In step 1, the time offset (intersymbol interference) was partially removed, but the variable frequency offset remained.-Step 2 removed the frequency offset up to an estimation error, but a constant phase offset remained.-Step 3 removed the time offset up to an estimation error, and only a phase offset remained, distinct from the one after step 2. -Finally, in step 4, the constant phase offset was removed up to an estimation error.Decoding of the data was possible at the end of stage 2 since alignment with the original sequence was achieved.
• Fig. 18 allows for a comparison of the constellation at the end of stage 2 in the simulated scenario and a setting without time, frequency, and phase offsets.This comparison helps estimate the Signal-to-Noise Ratio (SNR) in the experiment.Based on the comparison with simulated scenarios at different SNRs, it can be inferred that the experiment's SNR was around 40 dB.
These figures collectively demonstrate the synchronization and compensation process and its effectiveness in addressing different types of offsets.
The experimental results presented in this section provide confirmation that the procedure is effective in a real-world multipath environment.The subsequent part of the section focuses exclusively on simulation results.

1) PHASE DISTORTIONS
Figs. 19 and 20 show the phase distortion of the user data signal for the steps of stages 1 and 2, respectively.The results were obtained from one simulation run and only one signal period (samples 1025-2048).The parameters were: SNR = 20 dB, σ = 1.5 × 10 −3 , ρ = 0.1.With the increase of σ , the phase distortions (the mean slope of the curves) of the signals at the beginning of the procedure increased.Further, the increase of ρ, the tendency to abrupt changes of frequency offset (the roughness of the curves) decreased.Note that in each stage the phase offset after step 2 was constant in time up to the estimation error.

2) PERFORMANCE OF THE ALGORITHMS IN STAGE 1
In the remaining part of this section of the paper, we demonstrate the performance of the proposed algorithms and the procedure through cumulative distribution functions (CDFs) of absolute error values.In this subsection, we present CDFs that depict the quality of the algorithms in each step.The results obtained using beacon signals are shown with solid We initially selected a set of parameters for which the errors followed Gaussian cumulative distribution functions (CDFs).These starting parameters included SNR = 20 dB, σ = 0.5×10 −3 , and ρ = 0.99.We then individually modified the values of each parameter to demonstrate their impact on the estimation accuracy.The results are presented in Figs.21  and 22.
As anticipated, steps 3 and 4 exhibited significantly better accuracy than steps 1 and 2, respectively.The findings indicate that in favorable scenarios, SNR primarily influences estimation accuracy.For SNR = 20 dB, the median time error in step 3 was approximately one 50th of a sample, the median instantaneous phase error for the step 2 algorithm was around one 14th of a radian (equivalent to roughly 4 degrees), and the median mean phase error for the step 4 algorithm was roughly one 110th of a radian (about 0.5 degrees).
Next, we present the results for a group of parameters in which the time estimation errors deviated from Gaussian CDFs and exhibited a threshold effect.The initial parameters for this group were SNR = 10 dB, σ = 5 × 10 −3 , and ρ = 0.1.Subsequently, we adjusted each parameter individually to assess its impact on the estimation accuracy.These results are displayed in Figs.23 and 24.
Compared to the previous group of parameters, we observe several differences.In step 1, improving σ had a more significant impact in terms of reducing deviation in the CDF curves than enhancing the SNR.In step 3, improvements in both σ and ρ were more beneficial in terms of reducing deviation than increasing the SNR.
Regarding phase estimation, the estimation in step 2 is relatively insensitive to the threshold effect observed in step 1 because the algorithm in step 2 treats time offsets as phase offsets.In step 4 (Fig. 24 b), there is a correlation with the error observed in step 3 (Fig. 23 b).

3) PERFORMANCE OF THE PROCEDURE
In this subsection, we present CDF curves illustrating the estimation error for the steps of stage 1 in the procedure.The parameters for these simulations were SNR = 10 dB, σ = 1.5 × 10 −3 , and ρ = 0.1.For the beacon signals, the accuracy of time estimation after steps 1 and 3 matched the accuracy of the algorithms in those steps.The same applied to the user signals after step 3.However, we point out that the error in step 1 of stage 1 for the user signals was higher than The results presented in this subsection focus on phase estimation for steps 2 and 4. Each figure compares the phase error after a given step with the error of the algorithm used in that step, as described in subsection VI-D2.
In the case of phase error after step 2, it includes a constant phase shift, ϕ II err , due to the time error in step 1, τ I err .This shift arises because the algorithm in step 2 (the APSA) interprets the time shift of a cisoid as the corresponding phase shift: ϕ II err = (ω c + ω N ) τ I err , (mod 2π).This phase error is not taken into account when assessing the quality of the algorithm but must be included when evaluating the quality of the overall procedure, as demonstrated in Fig. 25 a).Due to the higher time error in step 1 for the user signal, the step 2 phase error is also higher for the user signal.It is important to note that the maximum absolute value of a phase error cannot exceed π.

FIGURE 18.
Constellation diagrams for the user data signal at the end of stage 2 for a) the scenario in Fig. 17 -simulation and b) a scenario with the same SNR and no offsets present in the signal -simulation.The red circles denote the combined user signal at the Rx, and the green circles denote the original sequence at the Tx.TABLE 1.The 0.9-quantile of time and instantaneous phase offsets at the end of stage 1, and the corresponding NBFG (the 0.1-quantile) as a function of SNR.
The algorithm in step 4 (the CPSA) estimates the mean phase error across the samples of a 1024-sample-long signal.However, what is important for evaluating the quality of the procedure after step 4 is the instantaneous phase shift, which represents the phase distortion of every single sample of the signal.This instantaneous phase error is significantly higher than the error of the algorithm, as shown in Fig. 25 b), and it is approximately equal to the error of the algorithm in step 2. This observation confirms our previous statement that step 4 compensates for the constant phase offset induced by step 2.

4) PERFORMANCE AS A FUNCTION OF SNR
The parameters for this subsection were chosen so that σ and ρ were equal to the corresponding ones in the USRP experiments: SNR ∈ {10, 15, 20, 25} dB, σ = 0.5 × 10 −3 , and ρ = 0.99.The results illustrate the quality of the procedure.The SNR curves for the 0.9-quantile of the absolute value of the time and instantaneous phase estimation error of the steps of stage 1 are shown in Figs. 26 a) and b), respectively.The time offsets in step 0 and after step 1 were limited to (−100, 100).The phase offsets change also in steps 1 and 3 as a part of the time offset compensation.These changes of the phase offset are of random nature so they are not shown in Fig. 26 b), and they are corrected in steps 2 and 4, respectively.Note that in Fig. 26 b) the phase offset values of the user signals after step 2 of stage 1 are equal as in step 0, since the user narrowband pilot is not used for the estimation in stage 1.For 20 dB, the 0.9-quantiles of the time estimation error and the instantaneous phase estimation error at the end of the stage 1 were one 20th part of a sample and 10 degrees, respectively.The corresponding NBFGs, i.e., curves for the 0.1-quantile at the end of stage 1 and stage 2 are shown in Fig. 27.We summarize the results from Figs. 26 and 27 in Table 1.
128214 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

VII. CONCLUSION AND FUTURE WORK
In this paper, we propose a procedure for over-the-air synchronization of a distributed dual-channel receiving system for the purpose of digital beamforming (BF) and radio source localization.The receivers are connected to a fusion center via digital links.We assume dominant LOS conditions, constant time offsets, and variable frequency offsets.For localization, the procedure utilizes wideband and narrowband pilots sent by a beacon to synchronize the receiving channels.For BF, in addition to the beacon signals, the procedure employs a user's wideband preamble and a narrowband pilot to synchronize the received user data signals.
For BF, the procedure is split into two stages with nearly identical steps, which are applied to different signals.aligns the signals in time using a non-coherent ML-type algorithm.
Step 2 removes the time-variable component of the phase offset, for which we propose a new adaptive algorithm.
Step 3 then performs fine time alignment, using an algorithm similar to the one in step 1.Finally, step 4 removes the constant component of the phase offset using a numerically efficient algorithm based on signal correlation.
The localization procedure contains only one stage, corresponding to Stage 1 of BF, where, instead of aligning the signals, it aligns the receiving channels, preserving TDOA in the received user signals.
The proposed procedures are modular, allowing every algorithm to be replaced by another algorithm of the same type.
We evaluated the performance of the procedure and the algorithms using both Monte Carlo simulations and experiments conducted on software-defined radio platforms, specifically USRPs.For simplicity, our analysis focused solely on signal synchronization, i.e., the BF case.
The experiments took place in an urban setting with LOS conditions.The USRPs featured variable frequency offsets and approximately constant time offsets throughout the observation interval.
Our adaptive algorithm for instantaneous phase offset estimation demonstrated the capability to effectively track abrupt changes in phase and frequency offsets.The accuracies obtained suggest that our procedure and algorithms are particularly well-suited for receive beamforming and noncoherent/semi-coherent localization [6].In cases where the frequency offset is predictable enough, they can also be applied to transmit beamforming.In spatially coherent environments, the procedure is well-suited for coherent localization.The results from our experiments closely matched the outcomes of the simulations, confirming the accuracy of the signal model and the underlying assumptions.
Future work may involve the optimal extension of our procedure to accommodate multichannel receiving systems and account for non-constant time offsets.Our procedure, which can effectively manage abrupt phase changes, could also benefit from testing in high-mobility settings involving Doppler shifts, including scenarios where network nodes are mounted on vehicles or drones.Another important aspect to explore is the behavior of the procedure in extreme multipath environments.It would be interesting to assess how the conditions in which this synchronization procedure is executed impact the performance of coherent localization.Additionally, a further extension could encompass a set of receivers, each equipped with an analog beamforming antenna array.

APPENDIX DETAILED DESCRIPTION OF THE SYNCHRONIZATION PROCEDURE
This section includes a detailed mathematical presentation of the steps of the proposed procedure for the BF and localization cases.Throughout this section we use the following equalities: which express the commutativity of the operator D. The same applies to .A time shift of τ applied to a phase shifted analog signal e jϕ(t) s(t), in the discrete-time matrix model is expressed as where −τ denotes compensated for −τ , i.e., delayed by τ .We note that is a unitary matrix, i.e., and if the Hermitian operator and time shifting are applied one after another, the following equality holds: To put it more generally, some of the signals in this section (vectors/matrices) have superscripts denoting time and phase shifts, where the superscripts mean that the signal is compensated for these shifts strictly in the order of the superscripts.For example, s ϕ(t),τ means that s is first compensated for a phase ϕ(t), i.e., a phase delay of −ϕ(t) is 128216 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.induced to the signal, and then compensated for a time shift τ , that is, a time delay of −τ is induced to the signal.

A. BF 1) STAGE 1
This stage corresponds to the path A 0 -A 1 -A 2 -A * 3 -A * 4 in Fig. 5. Let us assume that signals s 1 and s 2 arrive at the fusion center, see Fig. 28.In this stage, the original sequence/waveform is considered unknown, so the offsets in the transmitters are irrelevant because they affect both s 1 and s 2 in the same way.Therefore, we write: The goal is to estimate and compensate for the time and instantaneous phase offsets in s 2 relative to s 1 .(Since the frequency offset is a time variable, we deal with the instantaneous phase offset instead.)Upon arrival at the fusion center, each s 1 and s 2 are passed through three filters to obtain the wideband pilots, s W1 and s W2 , the narrowband pilots, s N1 and s N2 , and the user signals, s U1 and s U2 .The vectors of Gaussian noise at the output of the filters are denoted by η Wm , η Nm , and η Um , respectively.In the course of this procedure, four estimation-compensation pairs, i.e., steps, take place.
In step 1, CTSA (details provided in Section V) uses s W1 and s W2 to estimate the time shift of the observation intervals in the channels, 12 = τ T2 − τ T1 (s W2 is delayed by − 12 relative to s W1 ).However, these signals also contain a time shift due to different propagation times from the beacon to the respective receiving antennas, TDOA B = τ PB2 − τ PB1 , so we, in effect, estimate − 12 + TDOA B .The estimate at the output of the algorithm is where τ I err is the error due to phase distortions and the Gaussian noise.This is a rough estimate, with an accuracy of the order of one fifth of a sample (for reasonable SNR and N ).All the signals in channel 2 are compensated for the estimated time shift, τ I , according to ) Inducing a time or phase offset does not change the statistical properties of the Gaussian noise.
In step 2, the adaptive algorithm for phase shift (APSA) (see Section V) uses s N1 and s τ 1 N2 to estimate the instantaneous phase offset between the channels.Although the user Tx sends its own narrowband pilot, s UN , it is not used in stage 1.Because the beacon is a part of our system, it is likely that its signals have higher SNR than the user signals, so its narrowband pilot is used for the compensation of the user signals as well.Its use is possible because the knowledge of the original sequence is not taken advantage of, and the receiving system (the only source of phase offset that counts in this stage) is common for both transmitters.Since the phase shift is estimated for the time offset τ I relative to the signals at the Rx system input, the phase compensation in the user signal has to be performed for exactly the same time offset.This is why the user signal in step 1 is compensated for the time shift estimated using the beacon signal.After step 2, the beacon and user signals are processed separately.
Since the narrowband pilot is a cisoid, we choose to represent the time shifts in s N1 and s τ 1 N2 by phase terms, s N1 = 1 e jϕ N1 + η N1 , (56) where The APSA at its output yields where II is a diagonal matrix containing the estimated instantaneous phase offsets as a part of e j ϕ II (t) , II err is a matrix with errors due to noise, τ I 2 H 1 is the desired term, and e −jϕ II err is a constant phase error term, where ϕ II err = (ω c + ω N ) τ I err (mod 2π), with ω N being the angular baseband frequency of the narrowband pilot.The phase error appears because at the end of step 1, there is a residual time offset of τ I err between the signals in the channels.Therefore, instead of the corresponding samples, the phase offset is calculated between the samples separated by τ I err , which is, for a cisoid, equivalent to the phase shift represented by e −jϕ II err .However, the main task of step 2 is to resolve the time-variable component of the phase offset.The residual constant phase offset is resolved later in 128218 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply. .
In order for this compensation to work properly, all the signals in a channel have to be aligned in time.This is already achieved in step 1 because all the signals in channel 2 are shifted by the same value τ I .
In step 3, the CTSA (Section V) estimates the time shift between s τ I , ϕ II (t) U2 and s U1 , see Fig. 29.Since an unknown phase offset is induced in step 2 (no phase synchronization, [27]), the information from the carrier phase cannot be used (there would be no improvement if we used a coherent algorithm).The algorithm at its output yields where τ UIII is the estimated time shift, and τ UIII err is the error due to noise.The error of the estimation, τ UIII err , is of the order of the 50th part of a sample (for reasonable SNR and N ).The index ''U'' denotes that the estimate is obtained using user signals.
Here we make a small digression.A time shift of an RF signal is inherently accompanied by a corresponding phase shift [27], regardless of whether this phase shift is used in the estimation.Let us assume that an RF signal, with a complex envelope s(t), and a carrier frequency ω c , experiences a time shift of τ .Then the complex envelope takes the form s(t − τ )e −jω c τ .Therefore, although the envelope time shift of τ UIII err can be neglected, the phase shift caused by τ UIII is significant.If the signal s U2 has a flat spectrum in its subband, this phase shift is given by ϕ τ UIII err = (ω c + ω U ) τ UIII err (mod 2π), where ω U is the central angular baseband frequency of the user signal.Now, s τ I , ϕ II (t) U2 is compensated for the estimated time offset, τ UIII , according to Note that constant phase terms, such as e jϕ II err , do not change when being time shifted.In step 4, the algorithm for constant phase shift (CPSA) (Section V) estimates the mean phase shift between s τ I , ϕ II (t), τ UIII U2 and s U1 .The signal at its output is given by where ϕ UIV err is the phase estimation error due to noise.We denote by C {•} the contribution of its argument to the total estimated constant component of the phase.Now, s is compensated for the estimated phase offset, ϕ UIV , that is, where, for any vector/diagonal matrix V , C {V } = e jC{V } .We note that the terms C − τ UIII II err and H

II err
τ UIII approximately cancel each other out.Also, if the frequency offset does not change noticeably over τ UIII (τ UIII is usually less than a sample thanks to the existence of step 1), and The remaining error of these two approximations, along with ϕ UIV err , is contained in Uerr .Further, neglecting the small envelope shift by τ UIII err , represented by the term D BB τ UIII err , we obtain The obtained user signal is now added to the corresponding channel 1 signal, producing the combined user signal s US , in order to achieve the BF gain (higher effective SNR), Due to lack of synchronization with the user Tx, decoding of the user data is not possible yet.If this decoding is needed, we proceed to stage 2, where this synchronization is achieved.

2) STAGE 2
In this stage, which corresponds to the path B 0 -B 1 -B 2 -B 3 -B 4 in Fig. 5, we estimate and compensate for the offsets in the combined user signal relative to the corresponding known original sequence/waveform, see Fig. 30, i.e., we perform equalization of the chain [the user Tx] -[the propagation channel] -[Rx channel 1].In this stage, it is necessary that the user Tx sends a known preamble, or in other words, that the portion of the user signal we deal with in stage 1 is known to the Rx system.
First, s US is passed through two filters to obtain the combined user narrowband pilot, s UNS , and the combined user data signal, s DS .The vectors of Gaussian noise at the output of the filters for the combined signals are denoted η UNS and η DS .Unlike in stage 1, the offsets in the transmitters cannot be neglected, as they have an impact on the original transmitted sequences on their way to the receivers.Therefore, we write As in stage 1, the estimation and compensation processes are carried out in four steps.
In step 1, the CTSA estimates the time shift between s DS and s D .The estimate at the output of the algorithm is given by where τ UI err is the error due to phase distortions and noise.This is a rough estimate with an accuracy of the order of one fifth of a sample (for reasonable SNR and N ).
In step 2, the APSA estimates the instantaneous phase offset between s UNS and s UN .Since the narrowband pilot is a cisoid, we choose to represent the time shift in s UNS as a phase term, that is, where UII is the diagonal matrix containing the estimated instantaneous phase offsets as a part of e j ϕ UII (t) , U IIerr is the error due to noise, and similarly to stage 1, e −jϕ UII err is a constant phase error term, where ϕ UII err = (ω c + ω UN ) τ I Uerr (mod 2π), with ω UN being the angular baseband frequency of the user narrowband pilot.Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.In order to simplify our notation, we used the same symbols for some variables within different stages, e.g., τ UIII , ϕ UIV , since it is clear from the context what they refer to.
Finally, after stage 2, decoding of the data is enabled.If decoding of the beacon data is needed, the procedure (stage 1, stage 2) is the same as for any user signal.

B. LOCALIZATION
The localization case corresponds to the path A 0 -A 1 -A 2 -A 3 -A 4 -A 5 in Fig. 5.As stated in Section IV, the processing up to A 2 is the same as for the BF case.
where τ III err is the error due to the noise.The estimation error, τ III err , is of the order of the 50th part of a sample (for reasonable SNR and N ).Unlike for BF, the procedure here needs to know the location of the beacon.Since the beacon is a part of our system, this condition is almost certainly satisfied, so the propagation times from the beacon can be calculated, and therefore the shift of the time axes of the channels can be calculated with the accuracy of the order of the 50th part of a sample.Additionally, as previously stated, SNR for the beacon signals can be made relatively high, and larger statistical sample can be used, so that the synchronization accuracy could be further improved.The wideband pilot and the user signal in channel 2 after the compensation in step 3 are obtained according to . (86) The estimate of step 4 is:

FIGURE 2 .
FIGURE 2. Signals at the Rx channels with their baseband frequencies.

FIGURE 4 .
FIGURE 4.An illustration of time shifts in a signal: a) a signal; b) the signal at Rx 1 and its time axis; c) the signal at Rx 2 and its time axis.

FIGURE 5 .
FIGURE 5. A flowchart of the proposed procedure for signal synchronization and channel synchronization.APSA, CPSA, and CTSA are the acronyms of the algorithms we propose in a later section.

FIGURE 6 .
FIGURE 6. Signals at the Rx channels for BF.

FIGURE 7 .
FIGURE 7. Signals at the Rx channels for localization.

FIGURE 8 .
FIGURE 8.A schematic diagram of the adaptive algorithm for instantaneous phase offset estimation in a two-channel Rx system.

FIGURE 9 .
FIGURE 9. Rx 1 signal amplitude spectrum from a) a simulation and b) an experiment (similar for Rx 2 ).

FIGURE 10 .
FIGURE 10.Normalized frequency offset autocorrelation for different values of ρ.

FIGURE 11 .
FIGURE 11.Our experimental setup from two different angles.

FIGURE 12 .
FIGURE 12. APSA coefficients for channel 2 in step 2 of stage 1 for a) a simulation and b) an experiment.

FIGURE 13 .
FIGURE 13. a) Actual and estimated phase offset between the Rx channels (stage 1, step 2) -simulation, b) the actual frequency offset between the Rx channels (stage 1, step 2) -simulation, c) the estimated phase offset between the Rx channels (stage 1, step 2) -experiment, d) realizations of the actual phase offsets for different ρ -simulation.

FIGURE 14 .
FIGURE 14.Real part of the narrowband pilot in the Rx channels a) at the input of the APSA (stage 1, step 2) -simulation, b) at the output of the APSA -simulation, c) at the input of the APSA -experiment, and d) at the output of the APSA -experiment.

FIGURE 15 .
FIGURE 15.Real part of the user data signal in the Rx channels a) before the compensation for the offsets (stage 1, step 0) -simulation, b) after the compensation for the offsets (stage 1, step 4) -simulation, c) before the compensation for the offsets -experiment, and d) after the compensation for the offsets -experiment.

FIGURE 16 .
FIGURE 16.Constellation diagrams for the user data signal at the end of stage 1 a) a simulation and b) experiment.

FIGURE 17 .
FIGURE 17.Constellation diagrams for the user data signal across the steps of stage 2 a) simulation -left: a), c), e), g), i), and experimentright: b), d), f), h), j).The red circles denote the combined user signal at the Rx, and the green circles denote the original sequence at the Tx.

FIGURE 19 .
FIGURE 19.Phase distortion of the user data signal throughout steps 0-4 (a)-e)) of stage 1.The blue line denotes the signal in channel 1, and the red line denotes the signal in channel 2.

FIGURE 20 .
FIGURE 20.Phase distortion of the user signal in steps 0-4 (a)-e)) of stage 2. The blue line denotes the original signal at the Tx, and the red line represents the combined signal at the Rx.
Stage 1 aligns the signal in Rx 2 to the signal in Rx 1 , whereas Stage 2 aligns the sum of these signals to the original sequence of the selected Tx.Each stage consists of four steps.Step 1 coarsely

FIGURE 21 .
FIGURE 21.CDF curves of the time estimation error of the algorithm in a) step 1 and b) step 3 of stage 1.

FIGURE 22 .
FIGURE 22. CDF curves of the phase estimation error of the algorithm in a) step 2 and b) step 4 of stage 1.

FIGURE 23 .
FIGURE 23.CDF curves of the time estimation error of the algorithm in a) step 1 and b) step 3 of stage 1 with a threshold effect.

FIGURE 24 .
FIGURE 24.CDF curves of the phase estimation error of the algorithm in a) step 2 and b) step 4 of stage 1 with a threshold effect.

FIGURE 26 .
FIGURE 26.0.9-quantile of the a) time and b) instantaneous phase offset of the steps of stage 1 as a function of SNR, for σ = 0.5 × 10 −3 , ρ = 0.99.

FIGURE 28 .
FIGURE 28.A schematic diagram for steps 1 and 2 of stage 1 for both BF and localization cases.The solid lines denote vectors, and the dashed lines denote scalars.

step 4 .,=
We assume that the frequency offset does not change significantly over τ I err , as it is of the order of one fifth of a sample.The rest of the signals in channel 2 are compensated for the estimated phase offset according tos τ I , ϕ II (t) W2 = H II s τ I W2 = e jϕ II errH II err 1 F H D −τ T1 +τ PB1 +τ I err Fs EkvW + η τ I , ϕ II (t) W2 e jϕ II err H II err 1 F H D −τ T1 +τ PU2 −TDOA B +τ I err Fs EkvU + η τ I , ϕ II (t) U2

FIGURE 29 .
FIGURE 29.A schematic diagram for steps 3 and 4 of stage 1 for BF.The solid lines denote vectors, and the dashed lines represent scalars.

2 τ UI 1 − 1 −
The combined user signals are compensated for the estimated time shift, τ UI , bys τ UI DS = F H D − τ UI Fs DS = τ UI err U F H D τ UI err Fs D + η = F H D − τ UI Fs UNS = 2 τ I τ UI err U F H D τ UI err Fs UN + η τ UI UNS .

FIGURE 30 .
FIGURE 30.A schematic diagram for stage 2 for BF.The solid lines denote vectors, and the dashed lines represent scalars.

FIGURE 31 . 1 ≈
FIGURE 31.A schematic diagram for step 3 and step 4 for localization.The solid lines denote vectors, and the dashed lines denote scalars.